Dynamical systems on hypergraphs

Networks are a widely used and efficient paradigm to model real-world systems where basic units interact pairwise. Many body interactions are often at play, and cannot be modelled by resorting to binary exchanges. In this work, we consider a general class of dynamical systems anchored on hypergraphs. Hyperedges of arbitrary size ideally encircle individual units so as to account for multiple, simultaneous interactions. These latter are mediated by a combinatorial Laplacian, that is here introduced and characterised. The formalism of the master stability function is adapted to the present setting. Turing patterns and the synchronisation of non linear (regular and chaotic) oscillators are studied, for a general class of systems evolving on hypergraphs. The response to externally imposed perturbations bears the imprint of the higher order nature of the interactions.

Starting on these premises, it is clear that many body interactions constitute a relevant and transversal research field that is still in its embryonic stage, in particular as concerns studies that relate to hypergraphs. Indeed, novel light could be shed on a large plethora of systems, usually defined on standard networks, by accounting for generalised hypergraph architectures. This paper aims at taking one first step in this direction, by expanding along different axis. We will begin by adapting to the hypergraph setting the Master Stability Function [31] formalism. We will then consider the condition for the emergence of Turing patterns [32] for reaction-diffusion systems on hypergraphs, the synchronisation of nonlinear oscillators [7] and of chaotic orbits. It is here anticipated that for theoretical progress to be made one needs to characterise the spectral properties of a properly defined operator, which implements diffusion on hypergraphs.
The Master Stability Function (MSF), is a powerful technique developed in [31] to analyse synchronisation and it basically amounts to performing a linear stability analysis around a given equilibrium orbit, for a system of coupled interacting units. A straightforward application of linear stability analysis is for instance found in the context of the celebrated Turing instability, once the reference orbit is indeed a homogeneous fixed point.
In his seminal paper [32], Alan Turing set the mathematical basis of pattern formation. Initially proposed to explain the richness and diversity of forms displayed in Nature, the theory elaborated by Turing is nowadays an universally accepted paradigm of self-organisation [33][34][35]. The onset of pattern originates from the loss of stability of an homogeneous equilibrium, as triggered by diffusion. Turing instabilities have been initially studied for systems defined on continuous spatial domains and regular lattices [36]. More recently, the realm of application of Turing ideas has been extended to account for reaction-diffusion dynamics hosted on a complex network [37] and other related structures, such as multilayer networks [38,39] or multigraphs [40] just to mention a few. It is hence a natural question to generalise these studies to the broad framework of hypergraphs.
Turing patterns emerge from the destabilisation of a homogeneous equilibrium, that is a stationary solution of the examined model. In many real cases, however the system is not bound to evolve close to a stationary solution, but instead displays periodic oscillations. Examples ranges from biology to ecology, passing through physics [7,41]: individual nonlinear oscillators can synchronise and thus exhibit a coherent collective behaviour. Synchronisation, the spontaneous ability of coupled oscillators to operate in unison, has been studied for systems interacting via a complex and heterogeneous network of interlaced connections. To the best of our knowledge, however, this analysis has never been attempted for systems defined on hypergraphs of the type here considered. Let us observe that, although similar in their conception, the works [42,43] deal with hypernetworks, namely a network where several different links can connect two nodes, also called multigraph in the literature. The interactions are hence pairwise.
The formalism of the MSF can be also applied to chaotic oscillators. The synchronisation of chaotic systems defined on hypergraphs has been studied in [29] by using the formalism of the MSF under two main assumptions: (i) the work has been limited to p-hypergraphs, namely assuming that all the hyperedges have the same size p; (ii) the coupling function was assumed to be invariant with respect to permutations of the nodes, within each hyperedge. In this paper, we will relax both assumptions to deal with general hypergraphs with heterogenous hyperedge size distribution and without putting forward any hypothesis on the form of the coupling function.
In a recent work [30], the synchronisation phenomenon has been studied resorting again to the MSF, but employing however a Laplace operator [44] which cannot account in full for the high order interaction at play. The employed operator is defined from the hyper-adjacency matrix, which is solely capable to encode for the number of incident hyperedges without gauging their sizes. Moreover authors assumed the coupling function to depend on the average (arithmetic or geometric) value of the involved variables. Again, both assumptions are relaxed in the present work, because our Laplace operator takes into account both the number of incident hyperedges but also their size. We will moreover make use of a generic coupling function.
The paper is organised as follows. We first review the formalism of hypergraphs and introduce a new combinatorial Laplace matrix for hypergraphs. We then turn to discussing the spectra of the newly introduced Laplacian by emphasising its localisation properties. Then we present three applications, following the logic path outlined above, and elaborate on the impact of high-order structures. We finally conclude and sum up of our results.

II. HYPERGRAPHS.
Let us consider an hypergraph H(V, E), where V = {v 1 , . . . , v n } denotes the set of n nodes and E = {E 1 , . . . , E m } the set of m hyperedges, that is for all α = 1, . . . , m: E i ⊂ V , i.e. an unordered collections of vertices. Note that if E α = {u, v}, i.e. |E α | = 2, then the hyperedge is actually a "standard" edge denoting a binary interaction among u and v. If all hyperedges have size 2 then the hypergraph is actually a network. If an hyperedge contains all its subsets, then we recover a simplicial complex.
We can define the incidence matrix of the hypergraph [64], e iα , which carries information on how nodes are shared among edges (see middle panel Fig. 1). More precisely With such a matrix one can construct the n × n adjacency matrix of the hypergraph, A = ee T , whose entry A ij represents the number of hyperedges containing both nodes i and j. Note that often the adjacency matrix is defined by setting to 0 the main diagonal. Let us also define the m × m hyperedges matrix C = e T e, whose entry C αβ counts the number of nodes in E α ∩ E β .
The adjacency matrix of the hypergraph allows one to define a Laplace matrix [30,44], whose entries are given by k i δ ij − A ij , where k i = j A ij denotes the number of hyperedges incident with node i. This matrix generalises the (combinatorial) Laplace matrix for networks. However it does not account in full for the higher-order structures encoded in the hypergraph. Notably, the sizes of the incident hyperedges are neglected.
To overcome this limitation, authors of [16] studied a random walk process defined on a generic hypergraph using a new (random walk) Laplace matrix. It is worth mentioning that the transition rates of the associated process, linearly correlates with the size of the involved hyperedges. Stated differently, exchanges are favoured among nodes belonging to the same hyperedge (weighted according to its associated size). This allows in turn to describe the tightness of high-order interactions among "close nodes". More precisely: where the entries of K H are given by andĈ is a matrix whose diagonal coincides with that of C and it is zero otherwise. From this random walk Laplace operator, one can straightforwardly derive the (combinatorial) Laplace matrix, that will be employed in this paper to investigate the effect of diffusion on higher-order structures. In the above equation, matrix D contains on the diagonal the values k H i = =i k H i and zeros otherwise. It is clear from its very definition that K H takes into account both the number and the size of the hyperedges incident with the nodes. It can also be noted that K H can be considered as a weighted adjacency matrix whose weights have been self-consistently defined to account for the higher-order structures encoded in the hypergraph (see right panel of Fig. 1).
It is worth emphasising that the dynamics defined on this weighted network is equivalent [45] to the dynamics on the corresponding hypergraph. This observation allows us to transport existing tools targeted to networks' analysis to the realm where nodes are made to interact via hypergraphs. In particular, studying linear dynamical systems evolving on a hypergraph amounts to operating with standard n × n matrices, where n stands for the number of nodes. In this respect, the analysis is straightforward, and avoid the complications that are to be faced when dealing with simplicial complexes, where tensors are instead involved (see Section IV).
Given a hypergraph one can construct the projected network, that is the network obtained by mapping the nodes belonging to a hyperedge into a clique of suitable size (see left panel Fig. 1). If the hypergraph contains only simple hyperedges, then this projection is invertible and, given a network, one can construct a unique hypergraph whose projection coincides with the network itself [16]. Let us observe that the projected network keeps track of the many body interactions only though the cliques, i.e. relying on binary exchanges.
Let us conclude this section by remarking that the operator L H , given by Eq. (3), admits (1, . . . , 1) T as eigenvector associated to the zero eigenvalue. This latter homogeneous solution can be stable, so resilient, to external perturbation for a system evolving on a hypergraph and subject to nonlinear reaction terms. Instabilities can alternatively develop, depending on the specific explored setting. These issues will be addressed in the following by assuming higher-order interactions encoded by the hypergraph, to link co-evolving populations. Inspecting the stability of this generalised class of reaction-diffusion systems, amounts to studying the spectra of the coupling operator. For this reason we shall begin hereafter by analysing the spectra of a hypergraph Laplacian. In the middle panel, a hypergraph is displayed. Hyperedges are coloured according to their size (blue for size 2, red for size 3 and green for size 4). The hypergraph's characteristics are encoded in the incidence matrix eiα. Here, the information on how nodes are shared among hyperedges is stored. For ease of visualisation, we coloured the entries of eiα by using the same colour-code that was used to highlight the size of the hyperedges. From the hypergraph, we can construct the projected network, specified by the adjacency matrix A  (25) belongs to a hyperedge of size 3 and to another one of size 4. It is therefore the most important of the collection and because of this it receives the largest weight, k H 25 = 5. Observe also the link (28): it belongs to two hyperedges of size 3 and it is assigned a weight k H 28 = 4, larger than the one associated to the links that insist on the hyperedge of size 4.

III. LOCALISATION OF EIGENVECTORS
One can prove [16] that L H is symmetric, non-negatively defined and its largest eigenvalue equals 0. Moreover, let (Λ α H ) 1≤α≤n be the set of its eigenvalues of L H , then Λ n H ≥ . . . Λ 2 H > Λ 1 H = 0, and its eigenvectors, φ α 1≤α≤n form an orthonormal basis, φ α · φ β = δ α β . As already observed, φ 1 ∝ (1, . . . , 1). Finally L H reduces to the Laplace matrix defined on networks once all the hyperedges have size 2. In the following we will denote by (Λ α ) 1≤α≤n the eigenvalues of the Laplace operator of the projected network, L, and ψ α 1≤α≤n the associated eigenvectors. Based on the well known properties of L and assuming the network to be connected, we have Λ n ≥ . . . Λ 2 > Λ 1 = 0 and the eigenvectors do form an orthonormal basis. Localisation of eigenmodes is a phenomenon relevant to many fields of science, e.g. the Anderson localisation in disordered systems [46,47], with a particular relevance to dynamics. For this reason we decided to start our analysis by studying the localisation properties of the Laplacian eigenvectors for the hypergraph (3) and compare them with the corresponding quantities obtained for the projected network. Results reported in Fig. 2 show that the localisation is more evident for a hypergraph, than for the associated projected network. In the left panel of Fig. 2, we present the eigenvectors for the Laplace matrix stemming from the hypergraph (ordered for increasing eigenvalue Λ α H ) as a function of the nodes indexes (ordered for increasing k H i ). In the right panel, the same quantity is displayed for the Laplace matrix computed from the projected network. In this latter case, the nodes are ordered for increasing degree. By visual inspection (entries larger than 0.015 are coloured in black while the remaining ones are drawn in white), one can clearly appreciate the dark squarish zones, associated to small or medium rank eigenvectors, which appear in the left panel of Fig. 2: eigenvectors are found with relatively large entries across many nodes, i.e. a strong localisation. On the right panel, similar structures are present but much weaker. A substantially analogous behaviour is observed for high ranked eigenvectors, e.g. α 400 in the left panel and α ∼ 500 in the right one, for which only few entries display very large values, pointing hence to an even stronger localisation (see the thin dark "line" in the top right corners in both panels).
To illustrate our results, we employed as projected network a Scale Free network made by n = 500 nodes, built by using the configuration model with γ = −2 and k min = 2 [5]. The associated hypergraph is obtained by transforming all the maximal m-cliques into hyperedges of size m. The distribution of hyperedges sizes is reported in Fig. 3.

FIG. 2: Laplacian eigenvectors.
We report the absolute values of the components of the eigenvectors, φ α , ordered for increasing eigenvalues and nodes degree (right panel) and nodes hyper degree (left panel). Entries larger than 0.015 are pictured in black, while the remaining ones in white. The projected network is a scale free network made of 500 nodes and generated by using the configuration method with γ = −2 and kmin = 2. The corresponding hypergraph is obtained from the latter by transforming all the m-cliques into hyperedges of size m.
A more quantitative measure of the localisation, can be obtained by using the Inverse Participation Ratio (IPR) [48]. For a n-dimensional vector, v, this is defined as The above quantity ranges in [1/n, 1], where the lower bound is attained for a vector with uniform entries. The upper limit is hit when all entries are 0 but one, which equals 1. In Fig. 4 we report the IPR computed for the eigenvectors of the hypergraph (blue dots) and the projected network (black dots) used in Fig. 2. We can observe that in the case of the hypergraph, the IPR is always larger than the homologous quantity computed for the projected network, except for very high ranked eigenvectors (say, the last 5 ones).
In the next section we will show that the localisation which manifests on hypergraph, leaves macroscopic imprints on the dynamics of systems subject to many-body, higher-order interactions. This issue will be discussed in the following Section.

IV. DYNAMICAL SYSTEMS ON HYPERGRAPHS
In the remaining part of this paper we will consider the behaviour of dynamical systems defined on hypergraphs. In particular, we will analyse the consequences of dealing with higher-order couplings, exploiting to this end the spectral We report the distribution of hyperedges sizes for a hypergraph whose Scale Free projected network is made by 500 nodes and built using the configuration method with γ = −2 and kmin = 2. One can observe the presence of relatively large hyperedges responsible for high-order interactions. characteristics highlighted above. More specifically, assume n copies of the same low dimensional dynamical system to be hosted on each node of a given collection. This defines the local dynamics of the inspected system. Units belonging to different nodes are assumed to interact through higher-order structures identified as hyperedges. Many body interactions promote a preferential interaction among nodes belonging to the same large hyperedge. The nodes can be imagined to identify different spatial locations. For this reason we will denote by aspatial the system composed by one isolated dynamical unit, and use spatial to refer to its multi-dimensional version made of mutually entangled components.
As already mentioned the newly introduced (combinatorial) Laplace matrix (3) admits a homogeneous eigenvector associated to the zero eigenvalue. This will allow us to probe (in)stability of interconnected systems evolving close to reference orbits. For the sake of completeness, we will consider three distinct applications that cover several relevant research domains. We will begin by imposing a generalised diffusive coupling among nodes as exemplified by the aforementioned Laplace matrix (3). Working in this framework, we will study the emergence of Turing patterns, that is the conditions that promote the emergence of a stable heterogeneous solution. We will then turn to considering the synchronisation between nonlinear oscillators, diffusively coupled via higher-order combinatorial Laplacians. Finally, we will analyse the synchronisation of chaotic oscillators, in the setting of interest where higher-order interactions are at play. The formalism of the Master Stability Function, will be used to tackle the problem analytically. Projected networks will be employed as reference benchmarks to bring into evidence the role of hypergraphs and related higher order interactions.
Consider a d-dimensional system described by local, i.e. aspatial, equations: and fix a reference orbit, s(t). Let us observe that the latter can also be a fixed point. Assume further n identical copies of the above system coupled through a hypergraph, namely each copy is attached to a node of a hypergraph. Moreover, each unit belongs to one (or more) hyperedge. Units sharing the same hyperedge are tightly coupled, due to existing many body interactions. In formulas: We report the IPR of the eigenvectors of the hypergraph (blue dots) and of the associated projected network (black dots) used in Fig. 2. We can observe that in the case of the hypergraph the IPR is always larger than that obtained for the corresponding projected network, while the eigenvectors associated to the largest eigenvalues are more localised for the case of the network.
where x i denotes the state of the i-th unit, i.e. anchored to the i-th node, ε the strength of the coupling and G a generic nonlinear coupling function. The elements C α α of matrix C denote the size of the hyperedge E α . The factor −1 account for the fact that j should be different from i. Recalling the definition of e iα one can rewrite the previous formula as where we have used the definition of k H i = j k H ij and L H ij given by (3). Let us stress once again that all the high-order structure is encoded in a n × n matrix and there is no need for tensors as in the case of simplicial complexes: this simplifies the resulting analysis.
By exploiting the fact that j L H ij = 0 for all i = 1, . . . , n, it is immediate to conclude that the aspatial reference solution s(t) is also a solution of Eq. (6). A natural question that arises is hence to study the stability of the homogeneous solution for the system in its coupled variant.
To answer to this question one introduce the deviations from the reference orbit, i.e. δx i = x i − s. Assuming this latter to be small, one can derive a self-consistent set of linear equations for tracking the evolution of the perturbation in time. To this end, we make use of the above expression in Eq. (6) and perform a Taylor expansion by neglecting terms of order larger than two, to eventually get: where DF(s(t)) (resp. DG(s(t))) denotes the Jacobian matrix of the function F (resp. G) evaluated on the trajectory s(t).
Remember that L H is symmetric. Hence, there exists a basis formed by orthonormal eigenvectors, φ α , associated to eigenvalues Λ α H (see Section III). We can then project δx i on this basis and obtains for all α: where δy α is the projection of δx i on the α-th eigendirection. Let us finally conclude this section by observing that from Eq. (8) one can derive the Master Stability Function, i.e. the most general framework to address questions that pertain to the stability of the reference orbit. Despite its generality, the latter can only be handled numerically, except very few exceptions. In the following we begin by studying the setting where s(t) is a constant solution. In this case the Eq. (8) simplifies because the right hand side is no longer time dependent and the problem reduces to a classical study of Turing instability. Indeed, the rightmost term in Eq. (6) can be seen as a sort of generalised Fickean diffusion (see Section IV A). If the reference orbit is instead periodic in time, one can investigate the conditions which drive the synchronisation of regular oscillators. In this case the Master Stability Function can be analysed by resorting to the Floquet machinery. In the following, we have however chosen to study the synchronisation of Stuart-Landau oscillators via higher-order couplings (see Section IV B). Working in this setting, the Master Stability Function becomes again time independent and the analysis closely resembles the one carried out for addressing the onset of Turing instabilities. As a final step, we will turn to studying the case where s(t) is a chaotic trajectory (see Section IV C).

A. Turing patterns on hypergraphs
The Turing instability takes place for spatially extended systems: a stable homogeneous equilibrium becomes unstable upon injection of a heterogeneous, i.e. spatially dependent, perturbation once diffusion and reaction terms are simultaneously at play. Let us first consider two generic nonlinear functions f (u, v) and g(u, v) describing the local dynamics Then assume to replicate such system on all the nodes of a given hypergraphs, and label u i and v i the corresponding concentration. Here the index i refers to the specific node to which the dynamical variables are bound. Finally, assume that two nodes, i and j, communicate if they belong to the same hyperedge and moreover the strength of the interaction (which results in an effective transport across the involved nodes) is mediated by both the number of shared hyperedges and their sizes. Indeed, nodes belonging to the same hyperedge exhibit a higher-order interaction and we consequently assume that spreading among them is more probable than with nodes associated with other hyperedges or smaller ones. From a microscopic point of view, imagine to deal with a walker belonging to a given node. The walker assigns to all its neighbours a weight that gauges the size of the hyperedges and the number of incident hyperedges, and then she performs a jump with a probability proportional to this weight. This represents a higher-order extension Ficks' law: the rate of change of u i is proportional tȯ where use has been made of matrix C, as introduced above. Recalling the definition of e iα one can rewrite the previous formula asu where we have used the definition of k H i = j k H ij and L H ij . So in conclusion a reaction-diffusion processes on hypergraphs, where the diffusion takes into account the higherorder interactions among nodes in the same hyperedge, can be described by the following system where D u and D v are effective diffusion coefficients of species u and v. At first sight, the above model seems to solely account for binary interactions. However, higher-order interactions are also present, as encoded in the matrix L H . This is thus a compact formalism allowing to overcome the computational issues intrinsic to simplicial complexes. Finally, let us observe that if the hypergraph is a network, i.e. the hyperedges have size 2, |E α | = 2 ∀α, then L H reduces to the standard Laplace matrix. Thus Eqs. (10) converges to the standard reaction-diffusion system defined on a network. The condition for the emergence of a Turing instability can be detected by performing a linear stability analysis about the homogeneous equilibrium. More precisely, the latter is assumed to be stable with respect to homogeneous perturbations, while it loses its stability for heterogeneous perturbations once diffusion is at play, D u > 0 and D v > 0. The linear stability analysis can be performed by following the standard procedure: (i) by linearising the model (10) around the equilibrium, (u i , v i ) = (ū,v) for all i; (ii) by expanding the perturbations on the eigenbasis of L H , and (iii) by calculating the dispersion relation, i.e. the linear growth rate λ α = λ(Λ α H ) of the eigenmode α, as a function of the Laplacian eigenvalue Λ α H . The linear growth rate is the real part of the largest root of the second order equation where J 0 = ∂uf ∂vf ∂ug ∂vg is the Jacobian matrix of the reaction part evaluated at the equilibrium (u i , v i ) = (ū,v), tr (resp. det) is its trace (resp. determinant). The concept of dispersion relation is close to that of Lyapunov exponent: the existence of eigenvalues Λ α H for which the dispersion relation takes positive values, implies that the system goes unstable via a typical path first identified by Alan Turing in his seminal work. At variance, if the dispersion relation is negative the system cannot undergo a Turing instability: any tiny perturbation fades away and the system settles back to the homogeneous equilibrium.
To provide a concrete example, we assume the reaction kinetic to be modelled by the Brusselator scheme [49,50]. This is a nonlinear model defined by f (u, v) = 1 − (b + 1)u + cu 2 v and g(u, v) = bu − cu 2 v, where b and c act as tunable parameters. We first show an example of Turing pattern emerging in both the hypergraph and its related projected network (the same used in the previous section). In the main panels of Fig. 5 the dispersion relations are reported: a subset of eigenvalues exist which is associated to positive values of the dispersion relation, for both the hypergraph -panel (a)-and the projected network -panel (b). In the insets of Fig. 5 we display the ensuing patterns. Nodes are ordered for increasing hyper degree (resp. degree) for the hypergraph (resp. the projected network). One can clearly observe that, in the case of the hypergraph, patterns are strongly localised in nodes associated to larger hyper degree.  (panel (b)). We report the time evolution of the concentration of the species ui(t) in each node as a function of time, by using an appropriate colour code (yellow associated to large values, blue to small ones). In the former case, nodes are ordered for increasing hyper degree while in the second panel for increasing degree. One can hence conclude that nodes associated to large hyper degrees display a large concentration amount for species ui. This yields a very localised pattern. The hypergraph and the projected network are the same used in Fig. 2.
From Fig. 5 one can also observe that the domain of definition of the eigenvalues for the hypergraph cover a much wider range, as compared to that associated to the projected network. This observation can open the way to settings where patterns emerge only for systems defined on top of hypergraphs and not on the corresponding projected networks. In this case, patterns are the result of the higher-order interaction among nodes. To challenge this scenario, let us consider a small network built by using the Barabási-Albert algorithm [1] with 20 nodes. For each iteration of the generative algorithm, 3 new nodes are attached to the already existing ones, according to a preferential attachment scheme; because of the small size of the network, our goal here is not to resolve the scale free nature of the network but to obtain a hierarchical structure where 3-cliques, and larger ones, are mutually connected. We identify the complete cliques and build the associated hypergraph by assuming each m-clique to form a hyperedge with size m. We then turn to considering the resulting hypergraph and the associated projected network as the underlying support for the dynamics (10). The dispersion relation can be computed (see main panels of Fig. 6): observe that the homogenous equilibrium is stable even in presence of diffusion on the network while it looses stability in the case of the hypergraph. In this latter setting Turing patterns are hence expected to develop. This can be checked by computing the time evolution of the species density u i (t) both on the hypergraph and the projected network. By inspection of Fig. 6 one can appreciate that heterogenous patterns develop in the former case (see inset in the panel (a) of Fig. 6). Patterns are instead lacking in the latter scenario, i.e. when the Brussellator model hosted on the projected network (see inset in the panel (b) of Fig. 6).  (panel (b)). In the former case eigenvalues are found for which the dispersion relation is positive (red dots) while for the system defined on the projected network this conclusion does not hold. The blue line represents again the dispersion relation for the Brusselator model on a continuous support. Insets: Because the condition for Turing instability is satisfied for the hypergraph Laplace matrix, Turing Patterns emerge on the hypergraph (panel (a)). At variance, patterns do not manifest on the projected network (panel (b)). We report the time evolution of the concentration of species ui(t), on each node, as a function of time by using a proper colour code (yellow associated to large values, blue to small ones). In the former case nodes are ordered for increasing hyper degree while in the latter for increasing degree. The hypergraph and the projected network are obtained by means of the Barabási-Albert algorithm [1] with 20 nodes. Every iteration, 3 newly added nodes are attached to the existing ones.

B. Synchronisation of Stuart-Landau oscillators on hypergraphs
In the previous section we studied the emergence of Turing patterns for reaction-diffusion systems defined on a hypergraph so as to account for many body interactions. These patterns originate from a symmetry breaking instability induced by an externally imposed perturbation acting on systems initially close to a stationary homogeneous equilibrium. In many relevant problems, systems display periodic solutions. It is therefore important to investigate the stability of isolated periodic orbits and, even more essential, to study the dynamics of extended systems which combine several replica of the same nonlinear oscillators. Imagine that individual oscillators are evolving in phase and introduce a non homogeneous perturbation. If the system is globally stable the perturbation gets eventually re-absorbed and the oscillators display a synchronous dynamics [7]. Otherwise the perturbation develops in time and the system evolves towards a distinct, heterogeneous, attractor.
To study the synchronisation via a hypergraph, we consider individual units obeying to a Stuart-Landau (SL) equation [51,52]. This is a paradigmatic model of nonlinear oscillators, often invoked for modelling a wide range of phenomena, from nonlinear waves to second-order phase transitions, from superconductivity and superfluidity to Bose-Einstein condensation [53] Besides, the SL equation can be considered as a normal form for systems close to a supercritical Hopf-bifurcation. In this respect, the results here presented are more general than the specific setting explored.
Consider an ensemble made of n nonlinear oscillators and label with W i their associated complex amplitude. Each oscillator obeys a complex Stuart-Landau equation where c 2 is a real parameter and i = √ −1. Let us observe that the former admits the limit cycle solution W LC (t) = e −ic2t .
We then assume the oscillators to be coupled via a many body diffusive-like interaction which can be described by the discrete Laplacian matrix (3), returning thus the system where c 1 is a second real parameters and K is a suitable parameter setting the coupling strength. Based on the properties of the Laplace matrix, one can prove that the limit cycle solution, W LC (t), is also a solution of Eq. (12).
To characterise the stability of the latter to heterogeneous perturbation, we rewrite W j using polar coordinates as: Assuming |ρ i (t)| and |θ i (t)| to be small, one can insert the (13) into Eq. (12) and then linearise the resulting equation, to get: Remark that, even if we are perturbing around a limit cycle, namely a time dependent solution, the coefficients of the linearised equations do not depend on time, owing to the specific structure of the GL equation. This observation will simplify the successive analysis, which will follow closely that discussed in the preceding section for the case of a Turing instability. In the next section we will instead deal with a problem for which the linearised dynamics yields a time dependent Jacobian.
To proceed further we expand the perturbations ρ j and θ j on the Laplacian eigenvectors basis inserting the latter into (14), and by using the orthonormality of the eigenvectors, we obtain: We put forward the ansatz of exponential growth for each mode, that is ρ α ∼ e λαt and θ α ∼ e λαt and we eventually obtain a condition formally equivalent to the dispersion relation Let us observe that λ(Λ 1 ) = 0, signifying that the reference orbit is a limit cycle and thus neutral stable. On the other hand if λ(Λ α ) is positive for some α > 1, the perturbation grows exponentially in time, and the initial homogeneous state proves unstable. Conversely, if λ(Λ α ) < 0, for every α, the perturbation fades away and the system converges back to the fully synchronised state. Expanding (17) for small KΛ α we get λ(Λ α ) ∼ −KΛ α (1 + c 1 c 2 ) + . . . . By recalling that λ(0) = 0, K > 0 and Λ α > 0 for α > 1, one can conclude [54] that λ(Λ α ) > 0 for some α if and only if 1 + c 1 c 2 < 0, that is a necessary and sufficient condition for the loss of stability of the fully synchronised solution.
The numerical results reported in Fig. 7 complement the analytical theory discussed above. In panel (a) of Fig. 7 we present the dispersion relation and the heterogeneous patterns emerging for both the hypergraph and the associated projected network, for K = 1, c 1 = 0.5 and c 2 = −10. The dispersion relation is positive over a finite domain and the patterns (represented by W j (t)) that develop as follow the instability are pretty localised. In panel (b) of Fig. 7, the parameters are set to the values K = 1, c 1 = 1 and c 2 = −0.9. The dispersion relation is non positive and the system displays synchronised oscillations: the imposed perturbation dies out and the oscillators evolve at unison. . In both insets nodes are ordered for increasing hyper degree, resp. degree, for hypergraph, resp. projected network. The localisation is stronger when the dynamics is hosted on the hypergraph. In panel (b) the chosen parameters (K = 1, c1 = 1 and c2 = −0.9) result in the emergence of a globally synchronised state, the dispersion relation being always negative. This can be appreciated by looking at the insets ((b1) for the projected network and (b2) for the hypergraph) where we report Wj as function of time.

C. Master Stability Function on hypergraphs
In the previous section we have analysed the synchronisation of an ensemble made of Stuart-Landau (SL) oscillators defined on a hypergraphs. To this end we employed a straightforward generalisation of the techniques presented in section IV A, when investigating the emergence of Turing patterns. The use of the dispersion relation has been made possible because, for coupled SL equations, the variational problem yields a time independent Jacobian, once evaluated on the periodic homogeneous solution (14). This is not true for generic nonlinear oscillators. To overcome this problem one can however resort to the formalism of the Master Stability Function (MSF) [31], as introduced above. The aim of this section is thus to study the MSF in its full generality for systems defined on hypergraph. In particular, we will set to analyse the synchronisation of nonlinear chaotic oscillators coupled through a hypergraph and compare the outcome of the analysis to that obtained when operating the system on the corresponding projected network.
Let us consider again Eq. (8) and replace now in the latter equation εΛ α by a generic parameter κ > 0 and thus also the projection δy α by a generic "perturbation" vector δy The largest Lyapunov exponent of Eq. (18) is called the Master Stability Function [31]. Let us denote it by λ(κ) to emphasise its dependence on the parameter κ > 0. If for all κ, λ(κ) < 0, then δy decays to 0. At variance, if there exists κ > 0 such that λ(κ) > 0, then δy will grow. Back to Eq. (8) one can conclude that if for a given ε there exists α such that λ(εΛ α ) > 0, then the associated δy α grows in time. Thus individual units deviate from the reference solution s(t). On the other hand if for all α one has λ(εΛ α ) < 0 then the system reaches a globally synchronised state: all units will follow at the unison the same chaotic orbit. Let us observe that λ(0) > 0 being the reference orbit, s(t), a chaotic one.
To proceed in the analysis we assume linear coupling functions [55]: in this way the MSF simplifies, since DG is a constant matrix. Moreover, we will assume the matrix DG to have only one non zero element, say DG ba which denotes the existence of a coupling between the a-th and the b-th component of x.
Let us observe that the variational equation still contains explicitly the time variable via the Jacobian of the reaction part, DF(s(t)), which is indeed evaluated on the chaotic orbit. Hence to compute the MSF we have to solve a non autonomous system of ODEs, to study the evolution of the norm of δy(t) and then use the definition of the maximum Lyapunov exponent λ(κ) = lim t→∞ log ||δy(t)|| t . This can result in a tricky exercise. Indeed if λ(κ) > 0 then the norm can quickly increase to produce an overflow. On the other hand, if λ(κ) < 0, then ||δy(t)|| shrinks below round-off error. For this reason we employed in our analysis the Mean Exponential Growth of Nearby Orbits (MEGNO) algorithm [56,57]. This is an improved chaos indicator that allows to rapidly discriminate between chaotic and regular orbits. The method makes it possible for the Lyapunov exponent to be consequently recovered. For these reasons, MEGNO has been largely used in the framework of planetary systems [58,59], satellites and spatial debris [60][61][62] and also generic nonlinear dynamical systems [57]. The method overcomes the above mentioned limitation by performing a sort of time average of the norm of the deviation vector (see Appendix A).
Without loss of generality we will use the Lorenz model [63] for a demonstrative application: In the following we will fix the model parameters to the "standard values", β = 2, σ = 10 and ρ = 28 for which the system exhibits the chaotic orbit with a "butterfly shape". Once we couple the above ODE using high-order interactions, i.e. the hypergraph, we get d dt where the constant 3 × 3 matrix E encodes for the coupling among the three variables and its entries take values 0 or 1. For instance if E 21 = 1 and otherwise E ij = 0, (noted for short 1 → 2) then the growth rate of the second variable, y, depends on the first one, x, that isẏ i ∼ −ε j L H ij x j (discarding the reaction part). We are now in a position to adapt the above described theory, i.e. linearise about the reference orbit and project the perturbation on the eigenbase of the Laplace matrix, to Eq. (8) for the case of the Lorenz system. We will in particular compute the MSF to check the stability of the homogeneous states obtained by replicating chaotic Lorenz trajectory on each node of the collection. In the main panel of Fig. 8 we report the MSF for the coupling scheme, 1 → 1, that in the classification proposed in [55], corresponds to class Γ 1 , namely the MSF is monotone decreasing and it has a single root. We consider in particular two values of the coupling strength ε = 3 (panel (a)) and ε = 10 (panel (b)). For (sufficiently) small coupling strength (panel (a)), the MSF evaluated on the discrete spectrum of the hypergraph Laplace matrix (green dots) is always negative and thus the system synchronises to the chaotic reference orbit, as shown in the inset (a 2 ). On the other hand the MSF for the projected network (red dots) takes positive values: the chaotic oscillators cannot synchronise, as we can appreciate from inspection of inset (a 1 ). For large enough coupling strength (panel (b)), both spectra yeld a negative MSF (green and red dots in panel (b)) and hence, in both cases, the systems do synchronise (see insets (b 1 ) and (b 2 )).
From these results one can draw a first conclusion. Once we fix the coupling strength ε, the sign of the MSF depends on the spectrum of the Laplace matrix for the hypergraph, L H . Similarly for the projected network. However, as we observed in Section III the eigenvalues of the hypergraph Laplacian extend over a large portion of the real axis, as compared to what it happens when considering the projected network. Hence the coupling scheme 1 → 1 favours the synchronisation on the hypergraph, provided the coupling strength is sufficiently small. Said figuratively, one can act on the "knob" ε and have the spectra to slide on the MSF: by progressively reducing the value of ε one can force the spectrum of the projected network to enter the zone where the MSF is positive, whereas for the same value of ε the spectrum of the hypergraph is still associated to a negative MSF.
In Fig. 9 we report a similar analysis for the coupling schemes 1 → 2 (a panel) and 3 → 3 (panel (b)). In the classification proposed in [55] the former corresponds the class Γ 2 , two zeros, while the latter to Γ 3 , three zeros. From the results shown in Fig. 9, one can conclude that the system behaves similarly for couplings 3 → 3 and 1 → 1: if the . For a small coupling strength (panel (a)), the MSF is negative in correspondence of the eigenvalues of the Laplace matrix defined on the hypergraph (green dots), while the MSF can assume positive values once evaluated on the spectrum of the Laplace matrix for the projected network (red dots). In the former case ,the system synchronises (see inset (a2)) while in the latter it does not (see inset (a1)). For larger coupling strengths (panel (b)) the MSF is negative for both the hypergraph and the projected network and thus, in both cases, the systems do synchronise (see insets (b1) and (b2)).
coupling is sufficiently large (here ε = 20), synchronisation is found on the hypergraph but not on the corresponding projected network. This generalises our previous observation to all couplings belonging to an odd class Γ 2m+1 .
The reported behaviour is reversed once we consider couplings that belong to an even class. As we can appreciate from inspection of Fig. 9 panel (a) one can choose a sufficiently small coupling to have the MSF negative on the projected network (red dots), while it takes positive values, when the problem formulated on the hypergraph (green dots).

V. CONCLUSION
In this work we took a step forward in modelling dynamical systems on networks. The aim of the work is to account for high-order interactions among coupled units. In particular we focused on the hypergraphs, a very versatile setting where to model systems endowed with many-body interactions. Indeed one can easily represent such high-order interactions via the hyperedge, so as to overcome the limitations intrinsic to dealing with binary exchanges.
Starting from a microscopic process which takes place on the hypergraph, i.e. a random walk biases toward the size and the number of hyperedges a node belongs to, we defined a new combinatorial Laplace operator which generalises the concept of diffusive interaction to a multidimensional setting. This operator reduces to the standard combinatorial Laplacian once the hypergraph converges back to an ordinary network. In this respect, the newly introduced Laplacian can be rationalised as a natural extension of the usual operator.
In this framework we considered dynamical systems defined on top of hypergraphs and analysed the stability of the associated homogeneous equilibria. In particular we extended the Master Stability Function to this formalism and investigated the specificity of Turing patterns for the generalised proxy of reaction-diffusion systems on hypergraphs. We also analysed the synchronisation of periodic and chaotic orbits, shedding light on the role exerted by high-order couplings.
In all the inspected cases, the spectral properties of the novel Laplace operator are central in shaping the ensuing patterns, which appear remarkably localised, as illustrated with reference to the Turing setting. Further, hypergraphs can enhance or impede the synchronisation, as compared to what it happens on the corresponding projected network, depending on the specificity of the imposed couplings.
FIG. 9: Master Stability Function and patterns for the Lorenz system II. We report the MSF for the Lorenz model using the linear couplings, 1 → 2 (a panel) and 3 → 3 (b panel). In the former case we can observe that, for the chosen value of the coupling strength, ε = 2.4, the projected network yields a negative MSF (red dots) and thus the systems synchronises (inset a1), Conversely, the hypergraph possesses unstable eigenmodes (green dots) and the system goes consequently unstable (inset a2). The opposite behaviour is displayed in the case of the 3 → 3 coupling for ε = 20 (panel (b)): here the projected network exhibits unstable eigenmodes (red dots) while the hypergraph shows a negative MSF (green dots). The inset (b1) testifies on the absence of synchronisation for the projected network, while in the inset (b2) synchronisation is shown to occur on the hypergraph. Calculating the MSF proves less stable for the setting analysed in panel (b). The blue solid line refers to the average computed over 200 independent runs of the MEGNO algorithm. The shaded region (light blue) refer to the associated standard deviation.

Appendix A: Compute the MSF using MEGNO
To compute the MSF one has to solve Eq. (8). By discarding the partitioning into reaction and coupling parts, one can rewrite the previous equation as that is a time dependent ODE, often named variational equation. The latter should thus be solved together with the evolution of the reference trajectory where again we combine in F the reaction and the coupling parts. Then calling δx(t) = δx(t; δx 0 ) the solution of the variational equation with initial datum δx(0) = δx 0 , the Mean Exponential Growth factor by Nearby Orbits (MEGNO) [56,57], can be defined as: where [δ(τ )] 2 = ||δx(τ )|| 2 = (δx(τ ), δx(τ )), i.e. the norm of the vector δx, being (·, ·) the scalar product. We also emphasised that the MEGNO is being computed with respect to the reference orbit s(t). A trivial computation gives: where H = (J T + J )/2 is the Hermitian part of J . Together with the MEGNO one usually defines also the (time)-averaged MEGNO: Y (t) could in principle display large oscillations for large t, so limiting its effective predictive power. At variance, it can be shown that the average-MEGNO is well behaved and allows to study the dynamics for long times. Indeed the main feature of the average-MEGNO (and/or the MEGNO) is to allow to distinguish between regular orbits, for whichȲ (t) → 0, and irregular orbits, for whichȲ (t) grows unbounded. More precisely,Ȳ (t) ∼ λt/2 where λ is the largest Lyapunov characteristic number (or maximal Lyapunov exponent) of the orbit s(t). Let us observe that for regular orbits, MEGNO is able to differentiate between periodic ones, Y (t) → 0, and quasi-periodic ones, Y (t) → 2.