Connecting Hodge and Sakaguchi-Kuramoto through a mathematical framework for coupled oscillators on simplicial complexes

Phase synchronizations in models of coupled oscillators such as the Kuramoto model have been widely studied with pairwise couplings on arbitrary topologies, showing many unexpected dynamical behaviors. Here, based on a recent formulation the Kuramoto model on weighted simplicial complexes with phases supported on simplices of any order k, we introduce linear and non-linear frustration terms independent of the orientation of the k + 1 simplices, as a natural generalization of the Sakaguchi-Kuramoto model to simplicial complexes. With increasingly complex simplicial complexes, we study the the dynamics of the edge simplicial Sakaguchi-Kuramoto model with nonlinear frustration to highlight the complexity of emerging dynamical behaviors. We discover various dynamical phenomena, such as the partial loss of synchronization in subspaces aligned with the Hodge subspaces and the emergence of simplicial phase re-locking in regimes of high frustration. Synchronization dynamics in the presence of higher order interactions is well represented through variations of the Kuramoto model and subject of current interest. Here, the authors study and characterize the behavior of the simplicial Kuramoto model with weights on any simplices and in the presence of linear and nonlinear frustration, defined as the simplicial Sakaguchi-Kuramoto model.

S ynchronization is an ubiquitous phenomenon observed in many complex systems across spatial and temporal scales 1 , from the firing patterns of neurons and the communication of fireflies, to the flow of traffic [2][3][4] . One of the most popular dynamical systems, capable of reproducing a wide range of observed synchronization behaviors, is the Kuramoto model of coupled oscillators [5][6][7] . Whilst the model was originally formulated in terms of all-to-all interacting oscillators, the interactions between oscillators are commonly considered in-homogeneous and represented with a graph, whose structure affects the resulting dynamics. For example, while the full synchronization of the oscillator population is usually a strong attractor for the dynamics irrespective of the underlying graph 1,7 , the transient dynamics on the path toward synchronization can reveal the modular structure of the oscillators' interactions 8 .
Beyond the structure of oscillator interactions, other variations of the Kuramoto model have been studied extensively, including: time-delayed interactions 9,10 , oriented or signed interactions 11,12 , time-varying parameters, stochasticity, and more, see 1 for a comprehensive review. Of particular interest for this study, the introduction of a frustration parameter 13 in the nonlinear term of the Kuramoto model, then known as the Sakaguchi-Kuramoto model, can produce rich dynamics [14][15][16][17] and appears in many applications 18,19 .
Recently, the study of higher-order interactions between elements of a system, that is, models with interactions involving more than two nodes, has garnered momentum and interest 20 . Higher-order interactions are typically represented with hypergraphs or simplicial complexes, both of which generalize the graph representation of pairwise interactions to instead encode three-, four-and higher-way interactions. Naturally, extensions of well known dynamical systems have been proposed to investigate the effect of higher order interactions on their behavior [21][22][23][24][25][26] .
The Kuramoto model-being a paradigmatic model for synchronization phenomena-is no exception. In this case, however, there are two main avenues to extend classical oscillator models to higher-order. The first approach maintains the usual setup of phases defined on nodes of a systems and upgrades the interactions to the polyadic case, e.g., using simplicial complexes as the underlying connectivity structure. Recent works investigated variations of this node Kuramoto model with higher-order interactions introducing various types of coupling terms 27,28 . These models display a rich variety of synchronization and desynchronization phenomena, as well as multi-stable behavior. The second approach instead promotes phase variables from nodes to higher-order simplices, thus defining phases for edges, triangles, and all higher-order interactions, coupled by boundary operators as generalized incidence matrices. Pioneering work in this direction, 24 showed that the edge dynamics projected onto the nodes and faces possesses explosive synchronization properties when specific nonlinear and non-local couplings are introduced between the two projections. More recently, a version of the same model with local coupling between orders was introduced 29 .
In this paper, we extend the simplicial Kuramoto model in 24 to include: (i) weights on any simplices with a precise mathematical formulation based on discrete differential geometry; and (ii) linear and nonlinear frustrations. We will refer to the linear frustrations as natural frequencies yielding non-fully synchronized stationary states, and to the non-linear frustrations as the higherorder generalization of the Sakaguchi-Kuramoto model 13 . The difficulty of introducing proper nonlinear frustration comes from the orientation of the simplices which make even a naive frustration orientation dependent. Here, inspired by previous work on higher-order random walks 23 , we lift the simplices to double their numbers with opposite signs, obtaining an equivalent formulation without frustration and an orientation independent frustration. We then study the resulting frustrated simplicial Kuramoto model on edges with numerical simulations of oscillators which internal frequency vector, the linear frustration term, lies in the kernel of the Hodge Laplacian, using several measures to quantify the resulting dynamics, such as Hodge decomposition, the order parameter and the largest Lyapunov exponent.

Results
Simplicial complexes and Hodge Laplacian. The central elements of the mathematical formulation of the Kuramoto model on simplicial complexes are the boundary operators and the related Hodge Laplacians, which are, respectively, generalizations to higher order structures of the graph incidence matrices and of the Laplacian operator. We briefly review the main concepts we will use in our work following 30 , see also 25 , with additional details in Methods. A k-simplex is defined by a set of k + 1 nodes (a 1-simplex is an edge, a 2-simplex is a triangle, etc.). A simplicial complex is defined as a set of simplices in which every face of a simplex is also a simplex. For our purposes, the relevant connectivity between k−simplices will be that induced by sharing a (k −1)-simplex as a face, e.g., triangles sharing an edge, or by being faces of a (k + 1)-simplex, e.g., edges belonging to the same triangle. A k-chain within a simplicial complex is a linear combination of k-simplices. We denote by n k the number of k-simplices of a complex, which is also the dimension of the k-chains and kcochains vector spaces, dual to k-chains. The coboundary operator N k and its dual N Ã k on a simplicial complex are defined using the generalized incidence matrices B T k 2 M n k n kþ1 which encode the topology of a simplicial complex, and the weight matrices W k , which are diagonal matrices of the k-simplices weights [Eq. 1] The weight matrices W k can be chosen in an ad-hoc fashion and no formal relations need to exist between the different order k. The only relative constraint is for the weights to be positive in order to remain in the realm of unsigned graphs. Note that our notation follows 30 which differs from the convention commonly used for these operators. Both act on k-cochains, defined as linear functional on the space of k-chains. The Hodge Laplacian of order k can then be written as [Eq. (2)] For k = 0, W 0 = I and W 1 = I, we obtain the graph Laplacian L 0 = D − A with A the simplicial complex 1-skeleton, namely the graph node adjacency matrix, and D the diagonal matrix of the nodes degree. The choice W 0 = D −1 defines the normalized graph Laplacian L norm The graph Laplacian L 0 can produce two types of dynamics. When acting on the left of a distribution f, it yields the consensus dynamics _ f ¼ L 0 f for any choice of W 1 while by acting on the right, it corresponds to the diffusion dynamics _ p ¼ pL 0 . Equally, both types of dynamics also exist for the edge Laplacian L 1 23,31 , defined as [Eq. 4] We refer the reader to the Methods section for the diffusion formulation of the weighted simplicial Kuramoto model. For the remainder of this paper we will use the standard consensus formulation, but we emphasize that our formulation is not restricted to consensus dynamics.
Simplicial Kuramoto model. The Kuramoto model 5 is typically formulated for a node phase dynamical variable θ 2 R n 0 , with natural frequencies ω ¼ ðω 1 ; ; ω n 0 Þ 2 R n 0 that are sitting on the nodes of a graph G = (V, E) (|V| = n 0 , |E| = n 1 ) and interact through the graph adjacency matrix A ij 2 R n 0 n 0 For simplicity, we will consider a unit coupling σ = 1 throughout the remainder of this paper and will thus omit it from here on. The unweighted node Kuramoto model can be equivalently formulated in vector form using the n 0 × n 1 incidence matrix B T 0 32 and a vector of internal frequencies ω as which is approximated by the Laplacian dynamics When ω i = ω for all i, it is customary to study the node Kuramoto model in a frame rotating at ωt and thus ignore the internal frequencies, yielding _ θ ¼ L 0 θ. The Kuramoto model is therefore a nonlinear extension of the consensus dynamics introduced above.
The weighted simplicial Kuramoto model is then given for a time-dependent k-cochain θ (k) , see 24,25 for the original equations, as or equally with the weight and incidence matrices We emphasize that the positions of the weight matrices are not arbitrary but constrained by the geometrical nature of the coboundary operators. The definition and interpretation of the weights themselves are defined by the system under study, which may themselves include geometrical constraints. The weighted model can be seen as an extension of the Kuramoto model on weighted graphs, where the weights represent heterogeneous couplings. The weights on simplices of different order can be coupled, but this is not a necessary requirement and each order can capture independent characteristics of the system studied. We briefly explore the effect of weights on the dynamics of the Sakaguchi-Kuramoto model that forms the basis of our numerical experiments and is introduced in the next section. In the limit where θ is close to the subspace kerðL k Þ, we recover the linear consensus dynamics _ θ ðkÞ ¼ L k θ ðkÞ . For k = 0 and a connected graph, the kernel consists of a constant vector, or full synchronization.
Similarly to the node Kuramoto, the internal frequencies of the oscillators can be introduced via a change of rotating frame θ (k) → θ (k) − h (k) t for any vector h ðkÞ 2 kerðL k Þ. Indeed, such a vector will leave invariant the nonlinear terms, due to the presence of the boundary operator, and thus only adds a constant drift to the phases. Again, if we consider k = 0 and a connected graph, the kernel of L 0 is the constant vector, corresponding to the stationary state of consensus dynamics, and thus, by extension, the node Kuramoto model in full synchronization. For higher-order Kuramoto models k > 0, the dimension of the kerðL k Þ corresponds to the number of k-dimensional holes, i.e., holes bounded by k-simplices, or -equivalentlyto the Betti number β k of the simplicial complex. For k = 0, the Betti number β 0 corresponds to the number of connected components, and the stationary states are given by the piece-wise constant vectors to which each component will synchronize. We did not introduce by hand any internal frequencies at this stage, as we will see in remainder of this section that they naturally emerge as a form of linear frustration.
Simplicial Sakaguchi-Kuramoto model. The frustration in the Kuramoto model was first introduced in the Kuramoto-Sakaguchi model 13 , and has been studied in the context of graph theory, where the graph topology can give rise to rich repertoires of stationary states such as chimera states 14,15 and remote synchronization 17 . The frustrated node Kuramoto model is usually written as [Eq. 9] where α 2 R n 1 is the edge frustration vector, often taken to be constant α ij = α 1 . This equation cannot be directly formulated using the incidence matrices because the relative sign between the difference of phases θ i − θ j and α ij must be independent from the orientation of edges. In the adjacency matrix formulation, the orientation of edges is 'hidden', because L 0 ¼ B T 0 B 0 and A = D − L 0 are independent of edge orientation, and the choice of ordering θ i − θ j , instead of θ j − θ i , is possible irrespective of the edge orientation. If one writes B T 0 sinðB 0 θ þ α 1 Þ, the resulting order in the difference of phases depends on the choice of edge orientation and will not be 'node-centered', i.e., the θ i term will not always appear in front.
Nevertheless, it is possible to introduce a frustration in the general formulation of the Kuramoto model in Eq. (7) with coboundary operators such that it reduces to the frustrated Kuramoto in Eq. (9) for k = 0 and remains orientation invariant for k + 1 simplices. Our construction uses two ingredients: (i) lift matrices 23 , defined as [Eq. 10] for any order k, and (ii) the projection onto the positive or negative entries of any matrix X, defined element-wise as [Eq. 11] where | ⋅ | denotes the absolute value function. The lift matrices create duplicates of simplices of order k with an orientation opposite to the original one, whilst the projection sets half of the doubled simplices to zero, i.e., removes them, based on their signs. One can define the lift of the coboundary operator as [Eq. 12] The projection to positive or negative entries is often used to define directed node graph Laplacians 30,33 by transforming the edge orientation to an edge direction with either [Eq. 13] With D out/in the diagonal matrices of out-or in-degrees and A dir the corresponding directed adjacency matrix, L 0,out models the directed diffusion dynamics written explicitly as L 0,out = D out − A dir and L 0,in corresponds to the directed consensus dynamics with L 0,in = D in − A dir .
As we have seen in the construction of the weighted simplicial Kuramoto model in Eq. (7), we are using the formulation that yields consensus dynamics. We will thus consider the associated projection onto the negative entries of the lifted simplicial Laplacian as [Eq. 14] First, we note that b L k ¼ L k , thus the application of the lift and the projection has a trivial effect on the Hodge Laplacian, but crucially it allows us to introduce the frustration via the linear frustration operator [Eq. 15] acting on any cochain x and arbitrary frustration cochain α k . We can now formulate the frustrated simplicial Kuramoto model as By construction, our formulation is independent of the orientation of the k + 1-simplices but not of the k-simplices because only the action of the k + 1 lift is non trivial as it acts inside the nonlinear part of the equation. From this point of view, α k is a linear frustration, whilst α k+1 is a nonlinear frustration. When the internal frequencies are all equal, α k can be any arbitrary vector not necessarily in kerðL k Þ, which corresponds to equal internal frequencies in the node Kuramoto. This can lead to a variety of dynamics, including partially synchronized dynamics or even non-stationary dynamics if its amplitude is large enough 24 .
For k = 0, we recover the frustrated node Kuramoto model in Eq. (9) as where the natural frequencies vector α 0 naturally appears from the frustration operator, whilst the rest vanishes with N −1 = 0. The case where k = 1 constitutes the main equation we consider in the rest of this paper, namely the frustrated edge simplicial Kuramoto model which is invariant under change of face orientations, but not under change of edge orientation. To lighten the notation, we will write θ for θ (1) in the remainder of the paper as we only consider systems with edge oscillators.
Hodge decomposition of the dynamics. The Hodge decomposition is an important tool to study the properties of simplicial complexes. Here, we use it to decompose the dynamics of the oscillators on the simplicial Kuramoto model to understand their properties in relation to the amount of frustration applied. The Hodge decomposition theorem states that the space of k-cochains can be decomposed into three orthogonal spaces 34,35 which can be seen as analogues to the gradient, harmonic and curl space respectively. When k = 1 the three orthogonal spaces are exactly the gradient, harmonic and curl space respectively. Any k-cochain θ (k) can thus be projected onto each subspace where θ (k−1) and θ (k+1) are the corresponding potentials. Here, instead of computing these potentials, as done for example in 24 , we project the k-cochain θ (k) onto each subspace using the projection operators where the matrices p grad and p curl are the orthonormal bases of the ranges of N k and N Ã kþ1 and p harm the orthonormal basis of the kernel of L k .
Simplicial order parameter. Probably the most popular and fundamental tool to measure the level of synchronization in a coupled dynamical system is the order parameter. It is usually defined as where j ¼ ffiffiffiffiffiffi ffi À1 p , and was introduced for the original Kuramoto model on a complete graph, i.e., where all oscillators are coupled. The generalization of the order parameter to any graph structure 32 can be expressed as where 1 n 1 is the unit vector of dimension n 1 , see Methods for the details. This formulation allows one to write the node Kuramoto model with uniform natural frequencies as a gradient flow of the form Notice that the usual minus sign is not needed as the order parameter is a concave function. Notice that only the cosine term is needed to express the gradient flow, while the other constant terms are needed for the normalization. For a simpler derivation, we will thus modify the normalization to define the simplicial order parameter (SOP) as where the normalization is C k ¼ 1 n kÀ1 Á W À1 kÀ1 1 n kÀ1 þ 1 n kþ1 Á W À1 kþ1 1 n kþ1 which corresponds to the weighted sum of nodes and faces of the simplex, or the combined number of nodes and faces for unweighted simplicial complexes, see Method for details.
As expected, R k = 1 if θ (k) is in the harmonic space, which corresponds to full synchronization. Notice that for k > 1, the harmonic space is in general not spanned by the constant vector, and full synchronization does not correspond to equal θ k j values on the k-simplices. The SOP generalizes the notion of full synchronization to the instantaneous phase vector to be in the harmonic space, where the phases are in general not equal, except in the node Kuramoto case. This type of harmonic synchronization is therefore akin to a simplicial phase locking, in which each higher-order phase evolves with a different proper frequency but overall the whole dynamics lives within the harmonic space, i.e., kerðL k Þ for the corresponding k. In addition, if the dimension of the harmonic space is larger than one, the fully synchronized state is in fact a linear combination of the basis vectors of the harmonic space. For α 2 = 0 and if α 1 is harmonic, i.e., α 1 2 KerðL 1 Þ, the particular linear combinations will be dictated by the choice of α 1 , or, if absent, by the choice of initial conditions. Thus our formulation extends the notion of full synchronization beyond constant phases to include a generalized harmonic phase locking.
As in the node Kuramoto case, this order parameter acts as a potential for the gradient flow formulation of the full k-order Kuramoto dynamics as Note that this formulation does not contain the harmonic natural frequencies which can be recovered, as before, via a change of rotating frame. Finally, we notice that in the case of the standard node Kuramoto, this measure corresponds to the weighted generalization of Eq. (24) with a different normalization factor  The single face triangle complex has no hole and thus the harmonic space of the Hodge Laplacian is of dimension zero. Therefore, to be in the full synchronization regime, defined as θ i = θ j &for all; i, j in the absence of harmonic space, one would expect that that state is only accessible for the non-frustrated model with α 1 = α 2 = 0 in Eq. (19). We will show it is not the case and the dynamics can still reach fullsynchronization. For simplicity, and without loss of generality, we will use α 1 as a constant vector in time with the same value on all three edges. In addition, whilst the model is invariant to face orientation, it is not invariant to edge orientation. We thus have two non-equivalent choices for edge orientation: (i) a fully oriented complex, or (ii) one edge oriented in the opposite direction, as shown in Fig. 1(d).
In (i) the fully oriented complex, all edges are equivalent and the frustrated simplicial Kuramoto model reduces to the scalar equation If |α 1 | < 1, any initial condition will converge, as time → ∞ to full synchronization with phase θ 1 ¼ 1 3 sin À1 ðÀα 1 Þ À α 2 À Á in the stationary state. Otherwise, the stationary solution will be periodic around a linearly increasing trend. In (ii), the case of a flipped edge orientation, only two edges are equivalent, yielding the following coupled differential equations _ We solve these equations numerically for values of α 1 ∈ [0, 2.5] and α 2 2 0; π 2 Â Ã and show in Fig. 1(a-c) three measures that help us characterize the ensuing dynamics.
In Fig. 1(a), we plot the SOP of Eq. (26) where we observe a full synchronization regime, R 2 1 ðθ ð1Þ Þ ¼ 1, for the non frustrated case with α 1 = 0, α 2 = 0, but also for a large region of the α 1 and α 2 parameter space. To understand this regime further in term of the Hodge decomposition of the stationary state, we show in Fig. 1(b-c) the slope of a linear fit, representing the drift, of the temporal evolution of the projection of the solution onto the gradient and the curl subspaces at stationarity. More precisely, we estimate the parameter a h , or slope, of the linear regression of the projections of θ(t) onto the grad, curl and harmonic spaces, as defined in Eq. (22). In addition, the white regions correspond to projections that are constant in time, while regions with vanishing slopes but non-constant projections are in dark blue-corresponding to a value of zero. The latter correspond to oscillating solutions around a constant value. As we will see below, these measures provide a finer, while still tractable, analysis of the solutions in the frustration parameter space than the order parameter. We highlight some important observations from these plots. First, the region where the gradient component of the solution is non-constant matches with the region where we observe a large drop in synchronization. This suggests that when the gradient, and curl, component of the phases becomes too large, the synchronization is abruptly reduced. Notice that this region is bounded below by α 1 = 1, as for the fully oriented case. Second, the region where the projection of the curl is not constant is strictly contained within the region of non-constant gradient. This is a general result that we show below.
In Fig. 1(e-f) we show two typical trajectories of non-constant gradient and non-constant gradient and non-constant curl (corresponding to the magenta and green markers on Fig. 1(a-c) respectively). We observe that the trajectory is a Lissajous curve when the curl component is constant and a more complex trajectory otherwise. The Lissajous behavior is simply explained by imposing a constant curl, θ 2 = 2θ 1 + δ for a constant δ, which reduces the coupled differential equations to a one dimensional dynamical system _ θ 1 ¼ Àα 1 À sinðα 2 Þ À sinðθ 1 Þ; parameterizing the speed of motion on this curve, represented in Fig. 1(e) as dots equally spaced in time.
Finally, in the regime with non-vanishing curl, upper right of Fig. 1(a-c), there exists a sharp transition along α 1 between almost vanishing and positive gradient slope while the curl projection grows continuously. In Fig. 1(g-h), we show two trajectories on each side of this transition, corresponding to the blue (zero gradient slope) and orange (non-zero gradient slope) dots in Fig. 1(a-c). The two trajectories are partially overlapping where the segments of the trajectories parallel to the sinðθ 1 Þ axis are switching sign of sinðθ 2 Þ. A more precise understanding of this transition in the context of dynamical system theory could be of interest but is beyond the scope of this work.
Although simple, this simplicial complex already displays interesting and non-trivial dynamical behavior of the simplicial Sakaguchi-Kuramoto model when the frustrations are turned on. However, this example does not contain a hole, i.e., there is no harmonic component to the dynamics. We explore the role of the harmonic component of the dynamic in the next section.
Synchronization and edge orientation. For our second example, we use a slightly larger simplicial complex which we display in Fig. 2(a) to study the properties of the dynamics in the presence of a hole. We previously mentioned, if α 1 2 KerðL 1 Þ and α 2 = 0, the dynamics will fully synchronize with the stationary state θ (1) = α 1 . Setting α 2 > 0 will perturb the stationary state by increasing the gradient and curl components, but may remain in a simplicial phase-lock for a wide range of parameters, including at high frustration.
In Fig. 2(b, e) without loss of generality, we fix α 1 = 0 and scan α 2 2 0; π 2 Â Ã and observe that for a given choice of edge orientation, the level of synchronization, as measured by the SOP, decreases with α 2 , while its standard deviation remain null, which is an indicator of simplicial phase-locking. The projections onto the gradient and the curl spaces are constant, while the harmonic projection is not. We also notice that the dynamics are very sensitive to changes in orientation. Reversing the orientation of the blue edge, Fig. 2(a), has a dramatic impact on the solution, with the gradient component becoming non-constant for some choices of α 2 , see Fig. 2(c, f). Reversing the orientation of both the blue and red edges make both the gradient and curl components non-constant (see Fig. 2(d, g). Similar to the first example of the single face triangle complex, we observe that the projection onto the curl is non-constant only if the projection onto the gradient is also non-constant.
We now show the existence of two critical values for α 2 corresponding to changes of regime: α 2,g when the gradient becomes non-constant and α 2,c when the curl becomes nonconstant, and that α 2,c ≥ α 2,g . For simplicity, we set α 1 = 0, but any α 1 2 KerðL 1 Þ can be considered. First, for small α 2 < α 2,g , i.e., when the gradient and curl component are constant in the simplicial phase-lock regime, the solution is of the form for a scalar Ω, h 2 KerðL 1 Þ and ϵ 2 KerðL 1 Þ ? small, and thus The term ðN Ã k V kþ1 Þ À 1 n kþ1 counts the number of k + 1-simplices adjacent to each k-simplex and is therefore a generalized degree. ðN Ã 0 V 1 Þ À 1 n 1 is simply the weighted node degree and ðN Ã 1 V 2 Þ À 1 n 2 the weighted edge degree.
In this linear approximation, there always exists a solution of Eq. (31) for ϵ, but, in the nonlinear regime, the presence of the sine function may prevent any solution to exist and the system will leave the phase-locked regime for α 2 > α 2,g . The exact value of α 2,g is difficult to find analytically, as we see in Fig. 2, it depends not only on the structure of the simplicial complex but on the edge orientation as well. In addition, if the dimension of the kernel of L 1 is larger than 1, the direction of the vector h in the harmonic space may also depend on α 2 . For small α 2 , the value of Ω is represented in Fig. 2(b-d) by the value of the harmonic slope (in green), and increases quasi-linearly as a function of α 2 as expected from Eq. (30). For larger values of α 2 , the previous linearization is not valid and cannot be used to correctly approximate the dynamics. However, the Hodge decomposition is still valid and the corresponding projections operators defined in Eq. (22) allow us to decompose the simplicial Kuramoto equation in its gradient and rotational parts as [Eqs. 32, 33] where θ (1) = θ g + θ c + θ h from the Hodge decomposition. Notice that the Hodge decomposition directly implies that N 1 θ (1) = N 1 θ c as θ g is in the range of N 0 , and N 1 N 0 = 0, thus only the curl component of θ, θ c , can survive in the sine terms containing N 1 , a similar argument applies for θ g . These two equations are coupled by the term P grad ðN Ã 1 V 2 Þ À sinðV 2 N 1 θÞ which vanishes if α 2 = 0 because P grad ðN Ã 1 V 2 Þ À sinðV 2 N 1 θÞ ¼ P grad N Ã 1 sinðN 1 θÞ ¼ 0, since the range of N 1 is orthogonal to the gradient space. The fact that, in the absence of nonlinear frustration, the curl, grad and harmonic projection of the dynamics are decoupled was already noted in 24 and is a direct result of the orthogonality of these three spaces. The nonlinear frustration makes the dynamics of the gradient projection depend on the solution of the curl projection. This coupling relies on the presence of the lift and projections which are necessary to preserve the independence on the face orientation of the dynamics.
The presence of the coupling in the gradient equation explains why the dynamics of the curl projection can be non-constant only if the dynamics of the grad projection is also non-constant. Indeed, for α 2,g < α 2 < α 2,c , the curl projection equation is stationary with a time-independent θ c,∞ solution, i.e., _ θ c;1 ¼ 0. The coupling term in the grad projection is then a constant, which we denote as δ, and the Kuramoto dynamics reduces to _ θ g ¼ Àδ À N 0 sinðN Ã 0 θ g Þ: These dynamics correspond to the edge Kuramoto model on the complex without faces and a non-harmonic natural frequency δ. It has a transition from synchronization to non-synchronization regime at α 2,g . For α 2 > α 2,c , the dynamics are non-stationary in the curl projection and in the gradient due to the presence of δ.
While the order of the transitions to non-stationarity for the different components hold whenever they exist, their existence and exact behavior is dependent on the orientation of the edges and the localization of holes. Remarkably, even the simple example we used here displays an abundance of varying behaviors, of which we have only described representative examples: we observe no transitions in Fig. 2(b), only a gradient transition as in Fig. 2(c) or two transitions as in Fig. 2(d) with a near singular re-phase-locked synchronization gap. In Fig. 2(c, d), we also observe a re-synchronization to a phase-locked regime for α 2 > α 2,c until π 2 , possibly a result of the small size of the complex and high degree of symmetry. Indeed, as we observe in the next section, this regime does not exist for larger, more irregular complexes, see Fig. 3, and is replaced by a more chaotic regime.
Larger simplicial complex. Until now, we have studied the frustrated simplicial Kuramoto dynamics on small simplicial complexes in order to study and understand in detail the effects of the frustration on the dynamics. As a final example for this paper, we consider a larger simplicial complex constructed from a Delaunay triangulation of random points on a plane around two circular holes as illustrated in Fig. 3(a). In Fig. 3(b,c), we show the same analysis as in Fig. 2 with the slope of the projections and the SOP. As expected, we observe more complex dynamics from the shape of these curves, obtained after averaging over 10 simulations with random initial conditions. In particular, we do not observe any re-synchronization for large α 2 but rather an even more complex set of dynamics as shown by the standard deviation of theSOP in Fig. 3(c).
To better quantify these complex dynamics, we compute the largest Lyapunov exponent 36,37 of the trajectories of each edge phase and show in Fig. 3(d) the mean and quartile of them for each value of α 2 . As soon as the dynamics are no longer constant, i.e., α 2 > α 2,g , the largest Lyapunov exponent is on average positive, but increases significantly for larger α 2 , clearly indicative of chaotic dynamics. We visually noticed two different regime of chaotic dynamics, which we highlighted with α 2,a and α 2,b which corresponds to the start of the decrease of the slope of the gradient and the curl projection, respectively. For α 2 > α 2,b , we also observe some sensitivity to initial condition on the value of the slope of the projection, the line thickness is the standard deviation across 10 simulations with random initial conditions. For example in (b), only the harmonic projection is non-constant in time. If the slope is present with a value of 0, the solution is oscillating around a fixed point, for example in (d) for gradient slope at large α 2 . In (e-g), we show the average and standard deviation across time of the simplicial order parameter: R 2 1 and σðR 2 1 Þ respectively for the simulations of (b-d). The order parameter is 1 for α 2 = 0 and decreases as the nonlinear frustration increases. The standard deviation of the order parameter allows us to detect in which regime the solution is non-constant in the gradient or curl space. With this simplicial complex, the solution is non-constant only when the grad and/or curl are non-constant, even for high α 2 .
These two regimes would be interesting to study in more detail, since a decrease of the slope of the projection can either suggests more synchronization, as in the examples of Fig. 2, or more random or chaotic dynamics, as in Fig. 3. Understanding the transition between these two regimes in term of the complexity of the simplicial complex, where complexity is for example measured by the number of holes, their relative localization and the symmetries of the simplicial complex, is an open problem. The Lyapunov exponent seems however a promising measure to identify the switching between the two regimes: for the example of Fig. 2(d) it remains at low values (see brown line in Fig. 3(d)) and does not increase for large α 2 .
Effect of weights. In this last example, we briefly explore the effect of weights on the Sakaguchi-Kuramoto dynamics using the general formulation of Eq. (19). Our formulation allows mathematically for arbitrary weights on simplices of any order, however, the choice of weights may be constrained by the nature of the modeled system, e.g., geometric constraints of lengths, areas, and volumes. Whilst the exact interpretation of weights depends mainly on the context, we nevertheless provide general guidelines regarding their effect. In the node Kuramoto model, weights on nodes can be interpreted as a modification of the underlying graph Laplacian and thus the dynamics it represents. For example, using the inverse degrees yields the normalized Laplacian. Edge weights are a natural mechanism to introduce heterogeneous interactions between oscillators and is known to give rise to interesting dynamics, e.g., metastable Chimera states 15 . In the frustrated edge Kuramoto model of Eq. (17), the focus of this paper, edge weights have a similar interpretation to edge weights in the node Kuramoto model: they quantify the strengths of the interactions for the node Kuramoto. Node and face weights both modulate the interaction strength between oscillators. Face weights parameterize the strength of the triple coupling between phases where the nonlinear frustration acts, so vanishing faces weights are another mechanism to control the effect of the nonlinear frustration. Node weights play a similar role, but are decoupled from the effect of the frustration. We finally point out that while the weights at the different simplicial orders can related to each others, it is not necessarily the case, except in limit cases such as an edge with zero weights cannot serve as a support for a face. A systematic exploration, both analytical and computational, of the role and effect of each type weights and their combinations is well beyond the scope of this paper. We present here a phenomenological description of the dynamics in a simple example as a preliminary to future work: we considered a simplicial complex comprised of two triangles, one full and one empty, sharing one face and varied the the weight w of the full face, see Fig. 4(a) for the complexes and Fig. 4(b) for the results of the simulation for various frustrations and w. In the limit where w = 0, the face vanishes and we have two holes, and no frustration from α 2 and the solution lies entirely in KerðL 1 Þ. For low weights, only a small region is synchronized with a large increase in the projection in the gradient and harmonic spaces, but no curl component. The projection on the curl space is non-vanishing only for larger weights on the face. Finally for w = 1, we recover a similar behavior to that of the simple triangle of Fig. 1. The structure of the projection suggest non trivial relations between the different components, as well as three regimes: non-vanishing and vanishing curl, and two gradients regime nested in the vanishing curl one.

Conclusion
In this work, we extend a previously introduced Kuramoto model on simplicial complexes 24,25 to include weights on any simplices We consider the simplicial complex of (a) obtained from a Delaunay triangulation of random points on a plane around two circular holes. Faces are represented in gray, and edge orientations with arrows. We set the linear frustration parameter α 1 = 0 and scan across the nonlinear frustration parameter α 2 . We plot the slope of the projection of the harmonic, gradient and curl component of the solution in (b). We highlight the various transitions (see main text) with vertical lines. If the slope value is absent, it corresponds to a constant projection. In (c), we show the average, R 2 1 , and standard deviation, σðR 2 1 Þ, across time of the simplicial order parameter R 2 1 . In addition to the two critical α 2 , we manually highlighted two possible points corresponding to the onset of chaotic dynamics regimes with α 2,a and α 2,b . The standard deviation of the curves in the top of (b), which are calculated over 10 simulations with random initial conditions, is zero, except for large α 2 where it is small and is represented by the thickness of the curves. In (d), the dark line represents the average largest Lyapunov exponent λ over edges, and the shaded gray region between the lower an upper quartile of the corresponding distribution. For the sake of comparison, we have included a light brown curve that is the mean largest Lyapunov exponent for the simplicial complex in Fig. 2(d).
as well as a linear and a non-linear frustration term to define the simplicial Sakaguchi-Kuramoto. This formulation naturally allows us to generalize the notion of synchronization, internal frequencies, and edge frustration, of the standard Kuramoto model.
Without frustration, the Kuramoto dynamics can be decomposed into three independent sub-systems aligned with the orthogonal spaces given by the Hodge decomposition. However, we have demonstrated that by adding frustration to the dynamics, the harmonic, gradient and curl subspaces become hierarchically coupled, see Eq. (32). The dynamics in the harmonic space is coupled to both the gradient and curl subspaces even in the absence of harmonic linear frustration. In the linear regime of small nonlinear frustration, the amplitude of the dynamics in the harmonic space is proportional to the amount of frustration. In the nonlinear regime of simplicial Sakaguchi-Kuramoto, we showed that the dynamics is highly varied, from constant to chaotic solutions. We have discovered that the edge orientation is of fundamental importance in the resulting Kuramoto dynamics and the change of orientation of one edge can be enough to significantly alter the dynamics. Understanding the precise relationship between the choice of orientation for a given simplicial complex and the resulting type of dynamics has remained elusive so far but would be an interesting topic to gain further understanding of these systems particularly in the context of control.
We foresee various interesting directions for further interrogation of our frustrated simplicial Kuramoto and also additional adaptions. Firstly, while we only provide a simple example of the effect of weights on the dynamics, we believe that a full exploration of the weights definition and their effect will open interesting avenue of research not only theoretically but also for applications of the simplicial Sakaguchi-Kuramoto model. Secondly, whilst we used consensus dynamics in our formulation, we also mentioned earlier the dual formulation of the diffusion We scan the frustration parameters for four different values of w > 1 in each column, discarding the trivial case w = 1 is trivial. We display the order parameter as well as gradient, curl and harmonic slopes on separate rows. The structure of the projection suggests non trivial relations between the different components, as well as three regimes: non-vanishing and vanishing curl, and two gradients regime nested in the vanishing curl one.
Kuramoto (see Method). Indeed, examining how the dynamics of the consensus and diffusion formulations deviate in the weighted setting could be of interest. Thirdly, we did not explore the possible interplay between the linear and nonlinear frustration. The linear frustration is known to have a transition between stationary and non-stationary regimes 24 , but may also affect the types of dynamics with nonlinear frustration, which can even be made non-constant. Finally, as we have shown with our three examples, the topology of the simplicial complex is crucial to determine the type of dynamics and in particular its complexity. A more complete characterization in term of graph theoretical or topological measures would be of interest to identify the criteria necessary for the transition between non-stationary and chaotic dynamics. In fact, particular geometries of simplicial complexes may support more specific types of dynamics, with maybe partial, cluster or metastable synchronizations.

Methods
Review of discrete geometry. We provide here a short review of discrete geometry, following the exposition in 30 to which we redirect the reader for an in depth and detailed exposition of all the notions introduced here. A simplicial complex is a collection of k-cells (node, edges, face, etc...), on which are defined k-chains as vectors with coefficients on each cell. The space of k-chains is denoted as C k , and the incidence matrix B T k maps k + 1-chains to k-chains, i.e., B T k : C kþ1 ! C k : Dual to the space C k of k-chains is the space C k of k-cochains, defined using the scalar product as duality pairing, i.e., for τ k ∈ C k and c k ∈ C k , the pairing is From this pairing, the dual of the incidence matrix is defined as but reduces to coboundary operator We have thus defined the incidence matrix and its dual, but acting on chains and cochains. The map between them is a metric represented by a diagonal matrix, or weight matrix W k : C k ! C k as c k = W k τ k , and it's inverse W À1 k : C k ! C k : Then, to obtain the dual of the incidence matrix to form diffusion or Kuramoto equations, we need to define the dual of simplicial complex and the Hodge operator. If n is the largest dimension of the k-cells in the complex, the dual of a complex is a complex where the k-cells of the primal complex are the (n − k)-cells of the dual complex. With this definition, one can see that the dual incidence matrices M k are defined as N T k ¼ M nÀkþ1 , that is the incidence matrix on k-cells is the transpose of the incidence matrix on n − k + 1-cells in the dual complex.
The Hodge operator maps a k-cochain x of the primal complex to the (n − k)cochain x * on the dual as k x; and the (n − k)-cochains y * on the dual complex to the k-chain y on the primal complex as We are now in the position to define N Ã k , the dual of the coboundary operator N k , which maps k + 1-cochains into k-cochains as : Finally, the Hodge Laplacian is Lift, projection and frustrations. We provide here more details on the derivation of the frustrated simplicial Kuramoto model and the resulting orientation invariance. The most general projected lifted L k Laplacian with the lift from Eq. (12) and the projections onto negative values is But we have the following propositions.
Proposition 1. The following holds for any k Proof. Using the relations V T k V k ¼ 2 and V T k V À k ¼ 1 For the down term of the Hodge Laplacian, we have which results in a non-lifted down term in the complete Laplacian. For the up Laplacian, the same computation applies, thus we have the result. □ Then, from the form of the frustration operator acting in the nonlinear term corresponding to the up Laplacian, we can only use the lift on the k + 1 simplices to get the lifted Laplacian in Eq. (14) of the main text.

Proposition 2.
The frustrated simplicial Kuramoto model is independent on the orientation of the k + 1 simplices.
Proof. The term of interest is L up ðN k ; θ ðkÞ Þ ¼ ðN Ã k V kþ1 Þ À sin V kþ1 N k θ ðkÞ þ α kþ1 À Á : If we change the orientation of a k + 1 simplex indexed by i the corresponding coboundary operator e N j has ð e N k Þ i ¼ ÀðN k Þ i . Then where P i permutes the rows i and 2i of the lifted matrix. Hence, we have the orientation invariance L up ð e N k ; θ ðkÞ Þ ¼ ðN Ã k V kþ1 P i Þ À sin P i V kþ1 N k θ ðkÞ þ α kþ1 À Á ¼ L up ðN k ; θ ðkÞ Þ; as the permutation of rows commute with the point-wise sine function. □ Simplicial order parameter. Following 32 , the generalization of node order parameter to any graph in Eq. (24) is obtained by rewriting the node order parameter for a complete graph with the graph incidence matrix ¼ n 0 þ 21 n 1 cosðB 0 θÞ ¼ n 2 0 À 2n 1 þ 21 n 1 cosðB 0 θÞ ¼ n 2 0 R 2 0;g where we used n 1 ¼ 1 2 ðn 2 0 À n 0 Þ, which holds for complete graphs. We then chose a simpler normalization to write the order parameter as an average over edges as with 1 n 1 the constant vector of ones of dimension n 1 , so that for full-synchronization, we still have R 2 0 ¼ 1. In order to obtain a proper generalization of the order parameter on simplicial complexes, one first has to notice that the order parameter generates the Kuramoto model as a gradient flow. A gradient flow is constructed from two elements: a convex potential function H(x) for x ∈ V with V a vector space, and a gradient structure K: V → V * where V * is the dual of V. The gradient is a symmetric operator while an anti-symmetric operator would result in a Hamiltonian equation, with its dynamics restricted to the level sets of the potential function. The gradient flow is then given as [Eq. 34] where δ δx : F ðVÞ ! V Ã is a variational derivative acting on the space F ðVÞ of functions of V to its dual V * .
In our case, the potential function is R 2 k defined in the main text in Eq. (26), the variational derivative is simply the gradient with respect to θ. The gradient results in a chain, as W À1 kÀ1 N Ã kÀ1 : C k ! C kÀ1 and W À1 kþ1 N k : C k ! C kþ1 ; and the gradient structure K is simply the weight matrix W k that converts the resulting chains to cochains.
Diffusion simplicial Kuramoto. All the derivations in this paper were done following the standard formulation of the Kuramoto model that corresponds to the consensus dynamics from the graph Laplacian. The effects of this choice are actually only noticeable once one considers weighted graphs or simplicial complexes in a geometrical setting as presented here. Indeed, without weights, both the diffusion and consensus dynamics are equivalent, as the graph Laplacian is symmetric. Hence, there is an obvious formulation of Kuramoto model with the diffusion interpretation, that is obtained simply by acting with the Hodge Laplacian on the θ (k) cochains from the right. The diffusion simplicial Kuramoto model is then or, explicitly with the weight matrices, For k = 0, the fully synchronized state is not the constant state as with the consensus dynamics but is proportional to diagðW À1 0 Þ, often given by the node degree for normalized graph Laplacian.
In order to include the frustration operator, one has to define it's 'dual' version acting on N Ã k and we obtain the frustrated diffusion simplicial Kuramoto model [Eq. 37] Such systems require non-trivial weight matrices to be different from the simplicial Kuramoto model presented in the main text and could be of interest for future studies.

Data availability
Data sharing not applicable to this paper as no datasets were generated or analyzed during the current study.