# Blog

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

11:31 22 agosto in Artículos por Website

### 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 following30, see also25, 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 nk the number of k-simplices of a complex, which is also the dimension of the k-chains and k-cochains vector spaces, dual to k-chains. The coboundary operator Nk and its dual $${N}_{k}^{* }$$ on a simplicial complex are defined using the generalized incidence matrices $${B}_{k}^{T}\in {M}^{{n}_{k}\times {n}_{k+1}}$$ which encode the topology of a simplicial complex, and the weight matrices Wk, which are diagonal matrices of the k-simplices weights [Eq. 1]

$${N}_{k}={B}_{k}\,,\,\,{N}_{k}^{* }={W}_{k}{B}_{k}^{T}{W}_{k+1}^{-1}.$$

(1)

The weight matrices Wk 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 follows30 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)]

$${L}_{k}={L}_{k}^{down}+{L}_{k}^{up}$$

(2)

$$:={N}_{k-1}{N}_{k-1}^{* }+{N}_{k}^{* }{N}_{k}.$$

(3)

For k = 0, W0 = I and W1 = I, we obtain the graph Laplacian L0 = 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 W0 = D−1 defines the normalized graph Laplacian $${L}_{0}^{norm}=I-{D}^{-1}A$$.

The graph Laplacian L0 can produce two types of dynamics. When acting on the left of a distribution f, it yields the consensus dynamics $$\dot{f}={L}_{0}f$$ for any choice of W1 while by acting on the right, it corresponds to the diffusion dynamics $$\dot{p}=p{L}_{0}$$. Equally, both types of dynamics also exist for the edge Laplacian L123,31, defined as [Eq. 4]

$${L}_{1}={B}_{0}{W}_{0}{B}_{0}^{T}{W}_{1}^{-1}+{W}_{1}{B}_{1}^{T}{W}_{2}^{-1}{B}_{1}.$$

(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 model5 is typically formulated for a node phase dynamical variable $$\theta \in {{\mathbb{R}}}^{{n}_{0}}$$, with natural frequencies $$\omega =({\omega }_{1},\ldots ,{\omega }_{{n}_{0}})\in {{\mathbb{R}}}^{{n}_{0}}$$ that are sitting on the nodes of a graph G = (V, E) (V = n0, E = n1) and interact through the graph adjacency matrix $${A}_{ij}\in {{\mathbb{R}}}^{{n}_{0}\times {n}_{0}}$$

$${\dot{\theta }}_{i}={\omega }_{i}-\sigma \mathop{\sum}\limits_{j}{A}_{ij}\sin ({\theta }_{i}-{\theta }_{j}).$$

(5)

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 n0 × n1 incidence matrix $${B}_{0}^{T}$$32 and a vector of internal frequencies ω as [Eq. 6]

$$\dot{\theta }=\omega -{B}_{0}^{T}\sin ({B}_{0}\theta ),$$

(6)

which is approximated by the Laplacian dynamics $$\dot{\theta }=\omega -{B}_{0}^{T}{B}_{0}\theta =\omega -{L}_{0}\theta$$ in the limit B0θ 1. 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 $$\dot{\theta }={L}_{0}\theta$$. 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), see24,25 for the original equations, as

$${\dot{\theta }}^{(k)}=-{N}_{k-1}\sin \left({N}_{k-1}^{* }{\theta }^{(k)}\right)-{N}_{k}^{* }\sin \left({N}_{k}{\theta }^{(k)}\right),$$

(7)

or equally with the weight and incidence matrices

$${\dot{\theta }}^{(k)}= -{B}_{k-1}\sin \left({W}_{k-1}{B}_{k-1}^{T}{W}_{k}^{-1}{\theta }^{(k)}\right)\\ -{W}_{k}{B}_{k}^{T}{W}_{k+1}^{-1}\sin \left({B}_{k}{\theta }^{(k)}\right).$$

(8)

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 $${\dot{\theta }}^{(k)}={L}_{k}{\theta }^{(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)}\in \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 L0 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 –equivalently– to 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 model13, 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 states14,15 and remote synchronization17. The frustrated node Kuramoto model is usually written as [Eq. 9]

$${\dot{\theta }}_{i}={\omega }_{i}-\mathop{\sum }\limits_{j}{A}_{ij}\sin ({\theta }_{i}-{\theta }_{j}+{\alpha }_{ij}),$$

(9)

where $$\alpha \in {{\mathbb{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}_{0}^{T}{B}_{0}$$ and A = D − L0 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}_{0}^{T}\sin ({B}_{0}\theta +{\alpha }_{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 matrices23, defined as [Eq. 10]

$${V}_{k}=\left(\begin{array}{l}{I}_{{n}_{k}}\\ -{I}_{{n}_{k}}\end{array}\right),$$

(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]

$${X}_{ij}^{\pm }=\frac{1}{2}\left({X}_{ij}\pm \left|{X}_{ij}\right|\right)\,\,\forall ij,$$

(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]

$${N}_{k}\to {V}_{k+1}{N}_{k}{V}_{k}^{T}\,.$$

(12)

The projection to positive or negative entries is often used to define directed node graph Laplacians30,33 by transforming the edge orientation to an edge direction with either [Eq. 13]

$${L}_{0,{{{{{{{\rm{out}}}}}}}}}={N}_{0}^{* }{N}_{0}^{+}\,,\,{{{{{{{\rm{or}}}}}}}}\,{L}_{0,{{{{{{{\rm{in}}}}}}}}}={({N}_{0}^{-})}^{* }{N}_{0}\,.$$

(13)

With Dout/in the diagonal matrices of out- or in- degrees and Adir the corresponding directed adjacency matrix, L0,out models the directed diffusion dynamics written explicitly as L0,out = Dout − Adir and L0,in corresponds to the directed consensus dynamics with L0,in = Din − Adir.

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]

$${\widehat{L}}_{k}={N}_{k-1}^{-}{N}_{k-1}^{* }+{({N}_{k}^{* }{V}_{k+1}^{T})}^{-}{V}_{k+1}{N}_{k}.$$

(14)

First, we note that $${\widehat{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]

$${{{{{{{{\mathcal{F}}}}}}}}}_{k}^{{\alpha }_{k}}({N}_{k}):x\;\mapsto {N}_{k}x+{\alpha }_{k},$$

(15)

acting on any cochain x and arbitrary frustration cochain αk. We can now formulate the frustrated simplicial Kuramoto model as

$${\dot{\theta }}^{(k)}= -{{{{{{{{\mathcal{F}}}}}}}}}_{k}^{{\alpha }_{k}}({N}_{k-1})\left[\sin \left({N}_{k-1}^{* }{\theta }^{(k)}\right)\right]\\ -{({N}_{k}^{* }{V}_{k+1}^{T})}^{-}\sin \left({{{{{{{{\mathcal{F}}}}}}}}}_{k+1}^{{\alpha }_{k+1}}({V}_{k+1}{N}_{k})\left[{\theta }^{(k)}\right]\right)$$

(16)

$$= -{\alpha }_{k}-{N}_{k-1}\sin \left({N}_{k-1}^{* }{\theta }^{(k)}\right)\\ -{({N}_{k}^{* }{V}_{k+1})}^{-}\sin \left({V}_{k+1}{N}_{k}{\theta }^{(k)}+{\alpha }_{k+1}\right).$$

(17)

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 enough24.

For k = 0, we recover the frustrated node Kuramoto model in Eq. (9) as

$${\dot{\theta }}^{(0)}=-{\alpha }_{0}-{({N}_{0}^{* }{V}_{1}^{T})}^{-}\sin \left({V}_{1}{N}_{0}{\theta }^{(0)}+{\alpha }_{1}\right),$$

(18)

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

$${\dot{\theta }}^{(1)}= -{\alpha }_{1}-{N}_{0}\sin \left({N}_{0}^{* }{\theta }^{(1)}\right)\\ -{({N}_{1}^{* }{V}_{2})}^{-}\sin \left({V}_{2}{N}_{1}{\theta }^{(1)}+{\alpha }_{2}\right),$$

(19)

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 spaces34,35

$${C}^{(k)}={{{{{{{\rm{Im}}}}}}}}({N}_{k-1})\oplus \ker ({L}_{k})\oplus {{{{{{{\rm{Im}}}}}}}}({N}_{k}^{* }),$$

(20)

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 $${\theta }^{(k)}={\theta }_{{{{{{{{\rm{g}}}}}}}}}^{(k)}+{\theta }_{{{{{{{{\rm{h}}}}}}}}}^{(k)}+{\theta }_{{{{{{{{\rm{c}}}}}}}}}^{(k)}$$ as follow

$${\theta }_{{{{{{{{\rm{g}}}}}}}}}^{(k)} ={N}_{k}{\theta }^{(k-1)}\\ {L}_{k}{\theta }_{{{{{{{{\rm{h}}}}}}}}}^{(k)} =0\\ {\theta }_{{{{{{{{\rm{c}}}}}}}}}^{(k)} ={N}_{k+1}^{* }{\theta }^{(k+1)}.$$

(21)

where θ(k−1) and θ(k+1) are the corresponding potentials. Here, instead of computing these potentials, as done for example in24, we project the k-cochain θ(k) onto each subspace using the projection operators

$${P}_{{{{{{{{\rm{grad}}}}}}}}} ={p}_{{{{{{{{\rm{grad}}}}}}}}}^{T}{p}_{{{{{{{{\rm{grad}}}}}}}}}\\ {P}_{{{{{{{{\rm{curl}}}}}}}}} ={p}_{{{{{{{{\rm{curl}}}}}}}}}^{T}{p}_{{{{{{{{\rm{curl}}}}}}}}}\\ {P}_{{{{{{{{\rm{harm}}}}}}}}} ={p}_{{{{{{{{\rm{harm}}}}}}}}}^{T}{p}_{{{{{{{{\rm{harm}}}}}}}}},$$

(22)

where the matrices $${p}_{{{{{{{{\rm{grad}}}}}}}}}$$ and $${p}_{{{{{{{{\rm{curl}}}}}}}}}$$ are the orthonormal bases of the ranges of Nk and $${N}_{k+1}^{* }$$ and pharm the orthonormal basis of the kernel of Lk.

### 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

$${R}_{0,c}^{2}(\theta ):=\frac{1}{{n}_{0}}{\left|\mathop{\sum }\limits_{i = 1}^{{n}_{0}}\exp \left(j{\theta }_{i}\right)\right|}^{2},$$

(23)

where $$j=\sqrt{-1}$$, 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 structure32 can be expressed as

$${R}_{0,g}^{2}(\theta ):=1+\frac{2}{{n}_{0}^{2}}{1}_{{n}_{1}}\cdot (\cos ({N}_{1}^{* }\theta )-1),$$

(24)

where $${1}_{{n}_{1}}$$ is the unit vector of dimension n1, 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

$$\dot{\theta }=\frac{1}{2}{n}_{0}^{2}{\nabla }_{\theta }{R}_{0}^{2}(\theta ).$$

(25)

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

$$\begin{array}{ll}{R}_{k}^{2}({\theta }^{(k)})=&\frac{1}{{C}_{k}}\left({1}_{{n}_{k-1}}\cdot {W}_{k-1}^{-1}\cos ({N}_{k-1}^{* }{\theta }^{(k)})\right.\\ &\left.+{1}_{{n}_{k+1}}\cdot {W}_{k+1}^{-1}\cos ({N}_{k}{\theta }^{(k)})\right),\end{array}$$

(26)

where the normalization is $${C}_{k}={1}_{{n}_{k-1}}\cdot {W}_{k-1}^{-1}{1}_{{n}_{k-1}}+{1}_{{n}_{k+1}}\cdot {W}_{k+1}^{-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, Rk = 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 $${\theta }_{j}^{k}$$ 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., $${\alpha }_{1}\in {{{{{{{\rm{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

$${\dot{\theta }}^{(k)}={C}_{k}{W}_{k}{\nabla }_{{\theta }^{(k)}}{R}_{k}^{2}({\theta }^{(k)}).$$

(27)

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

$${{{{{{{{\mathcal{R}}}}}}}}}_{0}^{2}({\theta }^{(0)})=\frac{1}{{C}_{0}}{1}_{{n}_{1}}\cdot {W}_{1}^{-1}\cos ({N}_{1}^{* }{\theta }^{(0)}),$$

(28)

where $${C}_{0}=\mathop{\sum }\nolimits_{i = 0}^{{n}_{1}}{({W}_{1}^{-1})}_{ii}$$ is the weighted sum of edges, or for an unweighted graph, the number of edges.

### Frustrated simplicial Kuramoto model on a face

To showcase the properties of the frustrated simplicial Kuramoto model, we begin with one of the simplest examples: the single face complex of a triangle graph. 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 full-synchronization. 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

$$\dot{\theta }=-{\alpha }_{1}-\sin (3\theta +{\alpha }_{2}).$$

If α1 < 1, any initial condition will converge, as time →  to full synchronization with phase $${\theta }_{\infty }=\frac{1}{3}\left({\sin }^{-1}(-{\alpha }_{1})-{\alpha }_{2}\right)$$ 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

$$\begin{array}{l}{\dot{\theta }}_{1}=-{\alpha }_{1}-\sin (-{\theta }_{1}+{\theta }_{2})-\sin (2{\theta }_{1}-{\theta }_{2}+{\alpha }_{2})\\ {\dot{\theta }}_{2}=-{\alpha }_{1}+2\sin (-{\theta }_{1}+{\theta }_{2})-\sin (2{\theta }_{1}-{\theta }_{2}+{\alpha }_{2}).\end{array}$$

We solve these equations numerically for values of α1 [0, 2.5] and $${\alpha }_{2}\in \left[0,\frac{\pi }{2}\right]$$ 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}_{1}^{2}({\theta }^{(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 ah, 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

$${\dot{\theta }}_{1}=-{\alpha }_{1}-\sin ({\alpha }_{2})-\sin ({\theta }_{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 ({\theta }_{1})$$ axis are switching sign of $$\sin ({\theta }_{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 $${\alpha }_{1}\in {{{{{{{\rm{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 $${\alpha }_{2}\in \left[0,\frac{\pi }{2}\right]$$ 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 non-constant, and that α2,cα2,g. For simplicity, we set α1 = 0, but any $${\alpha }_{1}\in {{{{{{{\rm{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

$${\theta }_{\infty }(t)={{\Omega }}th+\epsilon ,$$

for a scalar Ω, $$h\in {{{{{{{\rm{Ker}}}}}}}}({L}_{1})$$ and $$\epsilon \in {{{{{{{\rm{Ker}}}}}}}}{({L}_{1})}^{\perp }$$ small, and thus

$${{\Omega }}h=-{L}_{1}\epsilon -{({N}_{1}^{* }{V}_{2})}^{-}{{{{{{{\bf{1}}}}}}}}{\alpha }_{2}.$$

(29)

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.

From Eq. (29), Ω, h and ϵ are defined as [Eqs. 30, 31]

$${{\Omega }}h=-{P}_{{{{{{{{\rm{harm}}}}}}}}}{({N}_{1}^{* }{V}_{2})}^{-}{{{{{{{\bf{1}}}}}}}}{\alpha }_{2}$$

(30)

$${L}_{1}\epsilon =-{P}_{{{{{{{{\rm{harm}}}}}}}}}^{\perp }{({N}_{1}^{* }{V}_{2})}^{-}{{{{{{{\bf{1}}}}}}}}{\alpha }_{2}.$$

(31)

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 L1 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]

$${\dot{\theta }}_{g}:={P}_{{{{{{{{\rm{grad}}}}}}}}}{\dot{\theta }}^{(1)}= -{N}_{0}\sin ({N}_{0}^{* }{\theta }_{g})\\ -{P}_{{{{{{{{\rm{grad}}}}}}}}}{({N}_{1}^{* }{V}_{2})}^{-}\sin ({V}_{2}{N}_{1}{\theta }_{c}+{\alpha }_{2})$$

(32)

$${\dot{\theta }}_{c}:={P}_{{{{{{{{\rm{curl}}}}}}}}}{\dot{\theta }}^{(1)}=-{P}_{{{{{{{{\rm{curl}}}}}}}}}{({N}_{1}^{* }{V}_{2})}^{-}\sin ({V}_{2}{N}_{1}{\theta }_{c}+{\alpha }_{2}),$$

(33)

where θ(1) = θg + θc + θh from the Hodge decomposition. Notice that the Hodge decomposition directly implies that N1θ(1) = N1θc as θg is in the range of N0, and N1N0 = 0, thus only the curl component of θ, θc, can survive in the sine terms containing N1, a similar argument applies for θg. These two equations are coupled by the term $${P}_{{{{{{{{\rm{grad}}}}}}}}}{({N}_{1}^{* }{V}_{2})}^{-}\sin ({V}_{2}{N}_{1}\theta )$$ which vanishes if α2 = 0 because $${P}_{{{{{{{{\rm{grad}}}}}}}}}{({N}_{1}^{* }{V}_{2})}^{-}\sin ({V}_{2}{N}_{1}\theta )={P}_{{{{{{{{\rm{grad}}}}}}}}}{N}_{1}^{* }\sin ({N}_{1}\theta )=0$$, since the range of N1 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 in24 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., $${\dot{\theta }}_{c,\infty }=0$$. The coupling term in the grad projection is then a constant, which we denote as δ, and the Kuramoto dynamics reduces to

$${\dot{\theta }}_{g}=-\delta -{N}_{0}\sin ({N}_{0}^{* }{\theta }_{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 $$\frac{\pi }{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 exponent36,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.

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 states15. 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 $${{{{{{{\rm{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.