Higher-order simplicial synchronization of coupled topological signals

Simplicial complexes capture the underlying network topology and geometry of complex systems ranging from the brain to social networks. Here we show that algebraic topology is a fundamental tool to capture the higher-order dynamics of simplicial complexes. In particular we consider topological signals, i.e., dynamical signals defined on simplices of different dimension, here taken to be nodes and links for simplicity. We show that coupling between signals defined on nodes and links leads to explosive topological synchronization in which phases defined on nodes synchronize simultaneously to phases defined on links at a discontinuous phase transition. We study the model on real connectomes and on simplicial complexes and network models. Finally, we provide a comprehensive theoretical approach that captures this transition on fully connected networks and on random networks treated within the annealed approximation, establishing the conditions for observing a closed hysteresis loop in the large network limit.

In networks, dynamical processes are typically defined over signals associated to the nodes of the network. In particular, the Kuramoto model [39][40][41][42][43] investigates the synchronization of phases associated to the nodes of the network. This scenario can change significantly in the case of simplicial complexes [16,17,19]. In fact, simplicial complexes can sustain dynamical signals defined on simplices of different dimension, including nodes, links, triangles and so on, called topological signals. For instance, topological signals defined on links can represent fluxes of interest in neuroscience and in biological transportation networks. The interest on topological signals is rapidly growing with new results related to signal processing [17,19] and higher-order topological synchronization [16,28,44]. In particular, higher-order topological * Correspoding author. Email:juanga@colorado.edu † Correspoding author. Email:jtorres@onsager.ugr.es ‡ Correspoding author. Email:ginestra.bianconi@gmail.com synchronization [16] demonstrates that topological signals (phases) associated to higher dimensional simplices can undergo a synchronization phase transition. These results open a new uncharted territory for the investigation of higher-order synchronization.
Higher-order topological signals defined on simplices of different dimension can interact with one another in non-trivial ways. For instance in neuroscience the activity of the cell body of a neuron can interact with synaptic activity which can be directly affected by gliomes in the presence of brain tumors [45]. In order to shed light on the possible phase transitions that can occur when topological signals defined on nodes and links interact, here we build on the mathematical framework of higherorder topological synchronization proposed in Ref. [16] and consider a synchronization model in which topological signals of different dimension are coupled. We focus in particular on the coupled synchronization of topological signals defined on nodes and links, but we note that the model can be easily extended to topological signals of higher dimension. The reason why we focus on topological signals defined on nodes and links is three-fold. First of all we can have a better physical intuition of topological signals defined on nodes (traditionally studied by the Kuramoto model) and links (like fluxes) that is relevant in brain dynamics [45,46] and biological transport networks [10,38]. Secondly, although the coupled synchronization dynamics of nodes and links can be considered as a special case of coupled synchronization dynamics of higher-order topological signals on a generic simplicial complex, this dynamics can be observed also on networks including only pairwise interactions. Indeed nodes and links are the simplices that remain unchanged if we reduce a simplicial complex to its network skeleton. Since currently there is more availability of network data than simplicial complex data, this fact implies that the coupled dynamics studied in this work has wide applicability as it can be tested on any network data and shows a simplicial complex in which not only nodes but also links sustain dynamical variables whose coupled synchronization dynamics is captured by the higher-order topological Kuramoto model. network model. Thirdly, defining the coupled dynamics of topological signals defined on nodes and links can open new perspectives in exploiting the properties of the line graph of a given network which is the network whose nodes corresponds to the links or the original network [47].
In this work, we show that by adopting a global adaptive coupling of dynamics inspired by Refs. [48][49][50] the coupled synchronization dynamics of topological signals defined on nodes and links is explosive [51], i.e., it occurs at a discontinuous phase transition in which the two topological signals of different dimension synchronize at the same time. We also illustrate numerical evidence of this discontinuity on real connectomes and on simplicial complex models including the configuration model of simplicial complexes [52] and the non-equilibrium simplicial complex model called Network Geometry with Flavor [12,13]. We provide a comprehensive theory of this phenomenon on fully connected networks offering a complete analytical understanding of the observed transition. This approach can be extended to random networks treated within the annealed network approximation. The analytical results reveal that the investigated transition is discontinuous.

II. RESULTS
A. Higher-order topological Kuramoto model of topological signals of a given dimension Let us consider a simplicial complex K formed by N [m] simplices of dimension m, i.e., N [0] nodes, N [1] links, N [2] triangles, and so on. In order to define the higher-order synchronization of topological signals we will make use of algebraic topology (see the Appendix for a brief introduction) and specifically we indicate with B [m] the m-th incidence matrix representing the m-th boundary operator. The higher-order Kuramoto model generalizes the classic Kuramoto model to treat synchronization of topological signals of higher-dimension. The classic Kuramoto model describes the synchonization transition for phases associated to nodes, i.e., simplices of dimension n = 0 (see Figure 1). The Kuramoto model is typically defined on a network but it can treat also synchronization of the phases associated to the nodes of a simplicial complex. Each node i has associated an internal frequency ω i drawn from a given distribution, for instance a normal distribution ω i ∼ N (Ω 0 , 1/τ 0 ). In absence of any coupling, i.e., in absence of pairwise interactions, every node oscillates at its own frequency. However in a network or in a simplicial complex skeleton the phases associated to the nodes follow the dynamical evolution dictated by the equationθ = ω − σB [1] sin B [1] θ , where here and in the following we use the notation sin(x) to indicate the column vector where the sine function is taken element wise. Note that here we have chosen to write this system of equations in terms of the incidence matrix B [1] . However if we indicate with a the adjacency matrix of the network and with a ij its matrix elements, this system of equations is equivalent tȯ valid for every node i of the network. For coupling constant σ = σ c the Kuramoto model [39][40][41] displays a continuous phase transition above which the order parameter is non-zero also in the limit The higher-order topological Kuramoto model [16] describes synchronization of phases associated to the n dimensional simplices of a simplicial complex. Although the definition of the model applies directly to any value of n, here we consider specifically the case in which the higher-order Kuramoto model is defined on topological signals (phases) associated to the links where φ r indicates the phase associated to the r-th link of the simplicial complex (see Figure 1). The higher order Kuramoto dynamics defined on simplices of dimension n > 0 is the natural extension of the standard Kuramoto model defined by Eq. (2). Let us indicate with ω the internal frequencies associated to the links of the simplicial complex, sampled for example from a normal distribution,ω ∼ N (Ω 1 , 1/τ 1 ). The higher-order topological Kuramoto model is defined aṡ φ =ω − σB [1] sin(B [1] φ) − σB [2] sin(B [2] φ). (6) Once the synchronization dynamics is defined on higherorder topological signals of dimension n (here taken to be n = 1) an important question is whether this dynamics can be projected on (n + 1) and (n − 1) simplices. Interestingly, algebraic topology provides a clear solution to this question. Indeed for n = 1, when the dynamics describes the evolution of phases associated to the links, one can consider the projection φ [−] and φ [+] respectively on nodes and on triangles defined as Note that in this case B [1] acts as a discrete divergence and B [2] acts as a discrete curl. Interestingly, since the incidence matrices satisfy B [1] B [down2] = 0 and B [2] B [1] = 0 (see Methods V) these two projected phases follow the uncoupled dynamicṡ where [1] and L down [2] = B [2] B [2] . These two projected dynamics undergo a continuous synchronization transition at σ c = 0 [16] with order parameters In Ref. [16] an adaptive coupling between these two dynamics is considered formulating the explosive higherorder topological Kuramoto model in which the topological signal follows the set of coupled equationṡ The projected dynamics on nodes and triangles are now coupled by the modulation of the coupling constant σ with the order parameters R down 1 and R up 1 , i.e. the two projected phases follow the coupled dynamicṡ This explosive higher-order topological Kuramoto model has been shown in Ref. [16] to lead to a discontinuous synchronization transition on different models of simplicial complexes and on clique complexes of real connectomes. However, signals of different dimension can be coupled to each other in non-trivial ways. In this work we will show how topological signals of different dimensions can be coupled together leading to an explosive synchronization transition. Specifically we focus on the coupling of the traditional Kuramoto model [Eq. (2)] to a higher-order topological Kuramoto model defined for phases associated to the links [Eq. (6)]. The coupling between these two dynamics is here performed considering the modulation of the coupling constant σ with the global order parameters of the node dynamics [defined in Eq. (4)] and the link dynamics [defined in Eq. (9)]. Specifically, we consider two models denoted as Model NL (nodes and links) and model NLT (nodes, links, and triangles). Model NL couples the dynamics of the phases of the nodes θ and of the links φ according to the following dynamical equationṡ φ =ω − σR 0 B [1] sin(B [1] φ) − σB [2] sin(B [2] φ). (13) The projected dynamics for φ [−] and φ [+] then obeẏ Therefore the projection on the nodes φ [−] of the phases φ associated to the links [Eq. (14)] is coupled to the dynamics of the phases θ [Eq. (12)] associated directly to nodes. However the projection on the triangles φ [+] of the phases φ associated to the links is independent of φ [−] and of θ as well. Model NLT also describes the coupled dynamics of topological signals defined on nodes and links but the adaptive coupling captured by the model is different. In this case the dynamical equations are taken to beθ For Model NLT the projected dynamics for φ [−] and for Therefore, as in Model NL, the dynamics of the projection φ [−] of the phases φ associated to the links [Eq. (18)] is coupled to the dynamics of the phases θ associated directly to nodes [Eq. (16)] and vice versa. Moreover, the dynamics of the projection of the phases φ on the triangles φ [+] [Eq. (19)] is now also coupled with the dynamics of φ [−] [Eq. (18)] and vice versa. Here and in the following we will use the convenient notation (using the parameter m) to indicate both models NL and NLT with the same set of dynamical equations given bẏ which reduce to Eqs. (13) for m = 1 and to Eqs. (17) for m = 2. We make two relevant observations: • First, the proposed coupling between topological signals of different dimension can be easily extended to signals defined on higher-order simplices providing a very general scenario for coupled dynamical processes on simplicial complexes.
• Second, the considered coupled dynamics of topological signals defined on nodes and links can be also studied on networks with exclusively pairwise interactions where we assume that the number of simplices of dimension n > 1 is zero. Therefore in this specific case this topological dynamics can have important effects also on simple networks.
We have simulated Model NL and Model NLT on two main examples of simplicial complex models: the configuration model of simplicial complexes [52] and the Network Geometry with Flavor (NGF) [12,13] (see Figure  2). In the configuration model we have considered powerlaw distribution of the generalized degree with exponent γ < 3, and for the NGF model with have considered simplicial complexes of dimensions d = 3 whose skeleton is a power-law network with exponent γ = 3. In both cases we observe an explosive synchronization of the topological signals associated to the nodes and to the links. On finite networks, the discontinuous transition emerge together with the hysteresis loop formed by the forward and backward synchronization transition. However the two models display a notable difference. In Model NL we observe a discontinuity for R 0 and R down 1 at a nonzero coupling constant σ = σ c , however R up 1 follows an independent transition at zero coupling (see Figure 2, panels in the second and fourth column). In Model NLT, on the contrary, all order parameters R 0 , R down 1 , and R up 1 display a discontinuous transition occurring for the same non zero value of the coupling constant σ = σ c (see Figure 2 panels in the first and third column). This is a direct consequence of the fact that in Model NL the adaptive coupling leading to discontinuous phase transition only couples the phases φ [−] and θ, while for Model NLT the coupling involves also the phases φ [+] .
Additionally we studied both Model NL and Model NTL on two real connectomes: the human connectome of Ref. [53] and the c. elegans connectome from Ref. [54] (see Figure 3). Interestingly also for these real datasets we observe that in Model NL the explosive synchronization involves only the phases θ and φ [−] while in Model NLT we observe that also φ [+] undergoes an explosive synchronization transition at the same value of the coupling constant σ = σ c .

A. Theoretical solution of the NL model
As mentioned earlier the higher-order topological Kuramoto model coupling the topological signals of nodes and links can be defined on simplicial complexes and on networks as well. In this section we exploit this property of the dynamics to provide an analytical understanding of the synchronization transition on uncorrelated random networks.
It is well known that the Kuramoto model is challenging to study analytically. Indeed the full analytical understanding of the model is restricted to the fully connected case, while on a generic sparse network topology the analytical approximation needs to rely on some approximations. A powerful approximation is the annealed network approximation [41] which consists in approximating the adjacency matrix of the network with its expectation in a random uncorrelated network ensemble. In order to unveil the fundamental theory that determines the coupled dynamics of topological signals described by the higher-order Kuramoto model here we combine the annealed approximation with the Ott-Antonsen method [43]. This approach is able to capture the coupled dynamics of topological signals defined on nodes and links. In particular the solution found to describe the dynamics of topological signals defined on the links is highly non trivial and it is not reducible to the equations valid for the standard Kuramoto model. Conveniently, the calculations performed in the annealead approximation can be easily recasted in the exact calculation valid in the fully connected case previous a rescaling of some of the parameters. The analysis of the fully connected network reveals that the discontinuous sychronization transition of the considered model is characterized by a non-trivial backward transition with a well defined large network limit. On the contrary the forward transition is highly dependent on the network size and vanishes in the large network limit, indicating that the incoherent state remains stable for every value of the coupling constant σ in the large network limit. This implies that on a fully connected network the NL model does not display a closed hysteresis loop as it occurs also for the model proposed in Ref. [21]. This scenario is here shown to extend also to sparse networks with finite second moment of the degree distribution while scale-free networks display a well defined hysteresis loop in the large network limit.

B. Annealed dynamics
For the dynamics of the phases θ associated to the nodes -Eq. (20) -it is possible to proceed as in the tra-ditional Kuramoto model [42,55,56]. However the annealed approximation for the dynamics of the phases φ defined in Eq. (21) needs to be discussed in detail as it is not directly reducible to previous results. To address this problem our aim is to directly define the annealed approximation for the dynamics of the projected variables φ [−] which, here and in the following are indicated as in order to simplify the notation. Moreover we will indicate with N = N [0] the number of nodes in the network or in the simplicial complex skeleton. Here we focus on the NL Model defined on networks, i.e., we assume that there are no simplices of dimension two. We provide an analytical understanding of the coupled dynamics of nodes and links in the NL Model by determining the equations that capture the dynamics in the annealed approximation and predict the value of the complex order parameters (with R 0 , R down 1 , Θ and Ψ real) as a function of the coupling constant σ.
We notice that Eq. (14), valid for Model NL, can be written asψ This equation can be also written elementwise aṡ where the vectorω is given bŷ Let us now consider in detail these frequencies in the case in which the generic internal frequencyω of a link follows a Gaussian distribution, specifically in the case in whichω ∼ N (Ω 1 , 1/τ 1 ) for every link . Using the definition of the boundary operator on a link it is easy to show that the expectation ofω i is given by Given that each node has degree k i , the covariance matrix C is given by the graph Laplacian L [0] of the network, i.e. where we have indicated with . . . c the connected correlation. Therefore the variance ofω in the annealed approximation is Moreover, the projected frequencies are actually correlated and for i = j we have It follows that the frequenciesω are correlated Gaussian variables with average given by Eq. (27) and correlation matrix given by the graph Laplacian. The fact that the frequenciesω i are correlated is an important feature of the dynamics of ψ and, with few exceptions (e.g., [57]), this feature has remained relatively unexplored in the case of the standard Kuramoto model. Additionally let us note that the average ofω over all the nodes of the network is zero. In fact where with 1 we indicate the N -dimensional column vector of elements 1 i = 1. By using the symmetry of the adjacency matrix, i.e. the fact that a ij = a ji , Eq. (31) implies that the sum ofψ i over all the nodes of the network is zero, i.e.
We now consider the annealed approximation consisting in substituting the adjacency matrix element a ij with its expectation in an uncorrelated network ensemble where k i indicates the degree of node i and k is the average degree of the network. Note that the considered random networks can be both sparse [58] or dense [59] as long as they display the structural cutoff, i.e. k i k N for every node i of the network. In the annealed approximation we can put Also, in the annealed approximation the dynamical Eq. (20) and Eq. (24) reduce tȯ where indicates the Hadamard product (element by element multiplication) and where two auxiliary complex order parameters are defined aŝ withR 0 ,Θ,R down 1 andΨ real.

C. The dynamics on a fully connected network
On a fully connected network in which each node has degree k i = N − 1 the dynamics of the NL Model is well defined provided its parameter are properly rescaled. In particular we require a standard rescaling of the coupling constant with the network size, given by which guarantees that the interaction term in the dynamical equations has a finite contribution to the velocity of the phases. The Model NL on fully connected networks requires also some specific model dependent rescalings associated to the dynamics on networks. Indeed in order to have a finite expectation ω i of the projected frequenciesω i and a finite of the covariance matrix C, [given by Eqs. (27) and (28), respectively] we require that on a fully connected network both Ω 1 and τ 1 are rescaled according to Considering these opportune rescalings and noticing that the order parameters obeyR 0 = R 0 ,R down 1 = R down 1 , Θ =Θ, and Ψ =Ψ, we obtain that Model NL dictated by Eqs. (34)- (35) can be rewritten here aṡ with R 0 , R down 1 , Θ and Ψ given by Eq. (23) and D. Solution of the dynamical equations in the annealed approximation

General framework for obtaining the solution of the annealed dynamical equations
In this section we will provide the analytic solutions for the order parameter of the higher-order topological synchronization studied within the annealed approximation, i.e., captured by Eqs. (34) and (35). In particular first we will find an expression of the order parameters R 0 of the dynamics associated to the nodes (Eq. (34)) and subsequently in the next paragraph we will derive the expression for the order parameter R down 1 associated to the projection on the nodes of the topological signal defined on the links (Eq. 35)). By combining the two results it is finally possible to uncover the discontinuous nature of the transition.

Dynamics of the phases of the nodes
When we investigate Eq. (34) we notice that this equation can be easily reduced to the equation for the standard Kuramoto model treated within the annealed approximation [42] if one performs a rescaling of the coupling constant σ R 0 → σ. Therefore we can treat this model similarly to the known treatment of the standard Kuramoto model [40][41][42]. Specifically, starting from Eq. (34) and using a rescaling of the phases θ according to it is possible to show that we can set Θ = 0 and therefore Eq. (34) reduces to the well-known annealed expression for the standard order Kuramoto model given bẏ Assuming that the system of equations reaches a steady state in which both R down 1 andR 0 become time independent, the order parameters of this system of equations in the coherent stateR 0 > 0 and R down 1 > 0 can be found to obey [40,42,51,55] and g(ω) is the Gaussian distribution with expectation Ω 0 and standard deviation 1.

Dynamics of the phases of the links projected on the nodes
In this paragraph we will derive the expression of the order parameters R down 1 andR down 1 which, together with Eqs. (44), will provide the annealed solution of our model. To start with we assume that the frequenciesω are known. In this case we can express the order parameters R down 1 andR down 1 as a function of the probability density function ρ (i) (ψ, t|ω) that node i is associated to a projected phase of the link equal to ψ. Since in the annealed approximation ψ i has a dynamical evolution dictated by Eq. (35) the probability density function obeys the continuity equation with associated velocity v i given by where we have defined κ i as In this case the complex order parameters are given bŷ In order to solve the continuity equation we follow Ott-Antonsen [43] and we express ρ (i) (ψ, t|ω) in the Fourier basis as Making the ansatẑ we can derive the equation for the evolution of b i = b i (ω i , t) given by Since we showed before that the average value ofψ i over nodes is zero, we look for non-rotating stationary solutions of Eq. (52), ∂ t b i = 0. As long as R 0 > 0 these stationary solutions are given by where d i is given by By inserting this expression into Eq. (49) we get the expression of the order parameters given the projected fre-quenciesω, in the coherent phase in which R 0 > 0 where, indicating by θ(x) the Heaviside function, we have defined Finally, if the projected frequenciesω are not known we can average the result over the marginal frequency distribution of the projected frequencyω i given by G i (ω) gettinĝ and an analogous equations for R down 1 sin(Ψ) (not shown). We note that in the case of distributions g(ω) and G i (ω) that are symmetric around their means the above equations always admit the solution Ψ =Ψ = 0. Such values of the phases are also confirmed by direct numerical integration of the NL model. These equations together with Eqs. (44) capture the steady-state behavior of the higher-order Kuramoto model coupling topological signals defined on nodes and links within the annealed approximation in the coherent synchronized phase. Note that by derivation, these equations cannot capture the asynchronous phase which is instead always a trivial solution of the dynamical equations corresponding to R 0 = R down 1 = 0. Finally we observe that for the NL Model as well as for the standard Kuramoto model on random networks, it is expected that the annealed approximation is more accurate for networks that are connected and are sufficiently dense.
To illustrate the applicability of the theoretical analysis, we consider two examples of connected networks with N = 1600 nodes: a Poisson network with average degree c = 12 and an uncorrelated scale-free network with minimum degree m = 6 and power-law exponent γ = 2.5 In Fig. 4  from direct numerical integration of Eqs. (20) and (25) and the steady state solutions obtained from the numerical solution of Eqs. (55). The backward transition is fully captured by our theory, while the next paragraphs will clarify the theoretical expectations for the forward transition.

E. Solution on the fully connected network
The integration of Eq. (57) requires the knowledge of the marginal distributions G i (ω) which does not have in general a simple analytical expression. However, in the fully connected networks with Gaussian distribution of the internal frequency of nodes and links this calculation simplifies significantly. Indeed, when the link frequencies are sampled from a Gaussian distribution with mean Ω 1 /N and standard deviation 1/(τ 1 √ N − 1), the marginal frequency distribution G i (ω) of the internal frequencyω i of a node i in a fully connected network is given by (see Methods for details) wherec = N N −1 . By considering Ω 0 = Ω 1 = ω i = 0, and performing a direct integration of Eqs. (57) we obtain (see Methods section for details) the closed system of equations for R 0 and R down where the scaling function h(x) is given by corresponding to the upper solution approach one (full phase synchronization), while those for the lower solution approach zero asymptotically, thus indicating that the incoherent state never loses stability. Indeed, it can be easily checked (see Methods for details) that for large σ the unstable solution of Eqs. (59) has asymptotic behavior with J 0 and J 1 constants given by Therefore the unstable branch approaches the trivial solution R 0 = R down 1 = 0 only asymptotically for σ → ∞. This implies that the trivial solution remains stable for every possible value of σ although as σ increases it describes the stationary state of an increasingly smaller set of initial conditions. This scenario is confirmed by numerical simulations (see Figure 5) showing that the backward transition is captured very well by our theory and does not display notable finite size effects. The forward transition, instead, displays remarkable finite size effects. Indeed, as σ increases, the system remains in the incoherent state until it explosively synchronizes at a positive value of σ and reaches the stable synchronized branch. However the incoherent state is stable in the limit N → ∞, and this forward transition is the result of finite size fluctuations that push the system above the unstable branch, causing the observed explosive transition. This is consistent with the fact that for larger values of N , which have smaller finite size fluctuations, the system remains in the incoherent state for larger values of σ.
Therefore, while a closed hysteresis loop is not present in the NL model defined on fully connected networks, we observe fluctuation-driven hysteresis, in which finite-size fluctuations of the zero solution drive the system towards the synchronized solution, creating an effective hysteresis loop.

F. Hysteresis on homogeneous and scale-free networks
In this section we discuss how the scenario found for the fully connected network can be extended to random networks with given degree distribution. We will start from the self-consistent Eqs. (57) obtained within the annealed approximation model. These equations display a saddle point bifurcation with the emergence of two non-trivial solutions describing a stable and an unstable branch of these self-consistent equations. These solutions always exist in combination with the trivial solution R 0 = R down 1 = 0 describing the asynchronous state. Two scenarios are possible: either the unstable branch converges to the trivial solution only in the limit σ → ∞ or it converges to the trivial solution at a finite value of σ. In the first case, the scenario is the same as the one observed for the fully connected network, and the trivial solution remains stable for any finite value of σ. In this case the forward transition is not obtained in the limit N → ∞ and the transition observed on finite networks is only caused by finite size effects. In the second case the trivial solution loses its stability at a finite value of σ. Therefore the forward transition is not subjected to strong finite size effects and we expect to see a forward transition also in the N → ∞ limit. in order to determine which network topologies can sustain a non-trivial hysteresis loop we expand Eqs. (57) for 0 < R 0 1, 0 <R 0 1, and 0 < R down 1 1 under the hypothesis that the distributions g(ω) and G i (ω) are symmetric and unimodal. Under these hypothesis it is easy to show that Eqs. (57) predict an unstable solution in which R 0 and R down 1 scale with σ according to with J 0 and J 1 constants given by As long as the network does not have vanishing J 0 and J 1 the unstable branch converges to the trivial solution R 0 = R down 1 only in the limit σ → ∞. This happens for instance for Gaussian distribution of the internal frequency of the links and converging second moment k 2 of the degree distribution. However, when the second moment diverges, i.e., the network is scale free with k 2 → ∞ as N → ∞, then R 0 and R 1 can converge to the trivial solution R 0 = R down 1 = 0 also for finite σ. This analysis suggests that the scenario described for the fully connected network remains valid for sparse (connected) networks as long as the degree distribution does not have a diverging second moment, while a stable hysteresis loop can be observed for scale-free networks.

IV. CONCLUSIONS
Until recently the synchronization phenomenon has been explored only in the context of topological signals associated to the nodes of a network. However, the growing interest in simplicial complexes opens the perspective of investigating synchronization of higher order topological signals, for instance associated to the links of the discrete networked structure. Here we uncover how topological signals associated to nodes and links can be coupled to one another giving rise to an explosive synchronization phenomenon involving both signals at the same time. The model has been tested on real connectomes and on major examples of simplicial complexes (the configuration model [52] of simplicial complex and the Network Geometry with Flavor [13]). Moreover, we provide an analytical solution of this model that provides a theoretical understanding of the mechanism driving the emergence of this discontinuous phase transition and the mechanism responsible for the emergence of a closed hysteresis loop. This work can be extended in different directions including the study of the de-synchronization dynamics of this coupled higher-order synchronization and the duality of this model with the same model defined on the line graph of the same network.

A. Definition of simplicial complexes
Simplicial complexes represent higher-order networks whose interactions include two or more nodes. These many-body interactions are captured by simplices. A ndimensional simplex α is a set of n + 1 nodes For instance a node is a 0-dimensional simplex, a link is a 1-dimensional simplex, a triangle is a 2-dimensional simplex, a tetrahedron is a 3-dimensional simplex, and so on. A face of a simplex is the simplex formed by a proper subset of the nodes of the original simplex. For instance the faces of a tetrahedron are 4 nodes, 6 links and 4 triangles. A simplicial complex is a set of simplices closed under the inclusion of the faces of each simplex. Any simplicial complex can be reduced to its simplicial complex skeleton, which is the network formed by the simplicial complex nodes and links. Simplices have a relevant topological and geometrical interpretation and constitute the topological structures studied by discrete algebraic topology. Therefore representing the many-body interactions of a complex system with a simplicial complex opens the very fertile opportunity to use the tools of algebraic topology [5,60] to study the topology of the system under investigation. In this work we show that algebraic topology can also shed significant light on the role that topology has on higher-order synchronization.

B. Oriented simplices and boundary map
In algebraic topology simplices are oriented. For instance a link α = [i, j] has the opposite sign of the link [j, i], i.e., Similarly to higher order simplices we associate an orientation such that where σ(π) indicates the parity of the permutation π.
It is good practice to use as orientation of the simplices the orientation induced by the labelling of the nodes, i.e., giving, for example, a positive orientation to any simplex where This will ensure that the spectral properties of the higherorder Laplacians that will be defined later are independent of the labelling of the nodes. Given a simplicial complex, a n-chain consists of the elements of a free abelian group C n with basis formed by the set of all oriented nsimplices. Therefore every element of C n can be uniquely expressed as a linear combination of the basis elements (n-simplices) with coefficients in Z 2 . The boundary operator ∂ n is a linear operator ∂ n : C n → C n−1 whose action is determined by the action on each n-simplex of the simplicial complex given by As a concrete example, in Figure 6 we demonstrate the action of the boundary operator on links and triangles.

FIG. 6:
The boundary operators and their representation in terms of the incidence matrices. Panel (a) and (b) describe the action of the boundary operator on an oriented link and on an oriented triangle respectively. Panel (c) shows a toy example of a simplicial complex and panel (d) indicates its incidence matrices B [1] and B [2] representing the boundary operators ∂1 and ∂2 respectively.
A celebrated property of the boundary operator is that the boundary of a boundary is null, i.e.
Let us indicate with N [n] the number of simplices of the simplicial complex with generic dimension n. Given a basis for the linear space of n-chains C n and for the linear space of (n − 1)-chains C n−1 formed by an ordered list of the n simplices and (n − 1) simplices of the simplicial complex, the boundary operator ∂ n can be represented as N [n−1] × N [n] incidence matrix B [n] . In Figure 6 we show a 2-dimensional simplicial complex and its corresponding incidence matrices B [1] and B [2] . Given that the boundary matrices obey Eq. (72) it follows that the incidence matrices obey for any n > 0.

C. Higher order Laplacians
Using the incidence matrices it is natural to generalize the definition of the graph Laplacian for n > 0. The higher-order Laplacian can be proven to be independent of the orientation of the simplices as long as the simplicial complex has an orientation induced by a labelling of the nodes. The most celebrated property of higher-order Laplacian is that the degeneracy of the zero eigenvalue of the n Laplacian L [n] is equal to the Betti number β n and that their corresponding eigenvectors localize around the corresponding n-dimensional cavities of the simplicial complex. The higher-order Laplacians can be used to define higher-order diffusion [17] and can display a higher-order spectral dimension on network geometries. Here we are particularly interested in the use of incidence matrices and higher-order Laplacians to define higher-order topological synchronization. Here we study Eqs. (44), (57) assuming that the distributions g(ω) and G i (ω) are unimodal functions symmetric about their means. Setting Ψ =Ψ = 0 and considering the change of variables z = ω/(σR 0 R down 1 ), y =ω/(σR 0 ), Eqs. (44) can be written as while Eqs. (57) reduce to We notice that the equations for R 0 ,R 0 and R down 1 do not depend on the order parameterR down 1 so we can obtain a fully analytical solution of the model without solving the last equation. The above equations depend on the distribution g(ω) and the set of marginal distributions G i (ω i ). However we can show that, provided k 2 / k is finite, the solution of these equations does not converge to the trivial solution R 0 =R 0 = R down 1 = 0 for any finite value of σ. Indeed we are now going to show that the unstable branch of the solution these equations converges to the trivial solution only in the limit σ → ∞. Assuming 0 < R 0 1, 0 <R 0 1 and 0 < R down 1 1 we can expand the functions g(zσk iR0 R down 1 ) and G i (yσR 0 k i ) as Stopping at the first order of this expansion we get This equations lead to the following scaling of R 0 and R down 1 with σ with This confirms the theoretical framework revealing that in this dynamics there is always a trivial solution R 0 = R 0 = R down 1 = 0 while Eqs. (44), (57) are characterized by a saddle-point instability so that for σ > σ c two additional solutions emerge, a stable solution and an unstable solution. The stable solution describes the synchronized phase and captures the backward transition. As long as the second moment of the degree distribution does not diverge, the unstable solution converges to the trivial solution R 0 =R 0 = R down 1 = 0 only for σ → ∞. The asymptotic scaling for R 0 and R down 1 given by Eq. (80) can be adapted to capture the asymptotic scaling of the fully connected case with a suitable rescaling of the model parameters of the model, obtaining Eqs. (61), (63).

E. Marginal distributions in the fully connected case
The distribution G 1 (ω) ofω is a Gaussian distribution with averages given by Eq. (27) and covariance matrix C given by Eq. (28). The covariance matrix has N − 1 eigenvalues given by λ = 1/τ 2 1 and one zero eigenvalue λ = 0 corresponding to the eigenvector This means that we should always have N n=1 [ω n − ω n ] √ N = 0, a constraint that we can introduce as a delta function in the expression for the joint distributionĜ(ω) of the vectorω. Here we note that under these hypotheses and assuming that the distribution of the frequencies of the links is a Gaussian with average Ω 1 /N and standard deviation 1/(τ 1 √ N − 1) the marginal probability G i (ω) of ω i can be expressed as Eq. (58).
Given that the covariance matrix has a zero eigenvalue we can express the joint Gaussian distributionĜ(ω) aŝ where δ(x) indicates the delta function and where F(ω) and C are given by The marginal probability G i (ω) is given by