Synchronization of Kuromoto Oscillators on Simplicial Complexes: Hysteresis, Cluster Formation and Partial Synchronization

TL;DR

This paper studies how oscillators synchronize on complex structures, showing how higher-order interactions and network geometry influence synchronization patterns.

Contribution

The study introduces new insights into synchronization dynamics on simplicial complexes, highlighting the role of topology and higher-order interactions.

Findings

Synchronization transitions can be abrupt or smooth depending on the structure and oscillator properties.

Higher-order interactions significantly affect synchronization patterns and can lead to cluster formation.

Topological effects and shared face geometry influence synchronization in systems like nanomaterials and the brain.

Abstract

The analysis of the synchronization of oscillator systems based on simplicial complexes presents some interesting features. The transition to synchronization can be abrupt or smooth depending on the substrate, the frequency distribution of the oscillators and the initial distribution of the phase angles. Both partial and complete synchronization can be seen as quantified by the order parameter. The addition of interactions of a higher order than the usual pairwise ones can modify these features further, especially when the interactions tend to have the opposite signs. Cluster synchronization is seen on sparse lattices and depends on the spectral dimension and whether the networks are mixed, sparse or compact. Topological effects and the geometry of shared faces are important and affect the synchronization patterns. We identify and analyze factors, such as frustration, that lead to these effects. We note that these features can be observed in realistic systems such as nanomaterials and the brain connectome.

Figures

Full text

1. Introduction

Complex networks can be seen in a huge variety of natural and engineered contexts. Such complex networks are characterized by quantities such as path lengths, clustering coefficients and betweenness centralities. The simplicial characterization of such networks has also revealed the existence of a hidden geometry in their structure and interactions. These networks include complex systems from the brain [1,2], designed materials [3,4,5,6] and physics problems [7,8], as well as the structures emerging from online social endeavors and large-scale social networks [9]. Such simplicial structures naturally include many-body interactions and can incorporate interactions of a higher order than the usual pairwise ones [10,11]. In the case of the dynamics of complex networks, the interplay between the dynamical evolution, the higher-order interactions and the structure of the substrate give rise to novel effects, which are a topic of recent investigations [12,13].

There are several different directions in the current research on the simplicial structure of networks. One direction is the examination of empirical data to identify the underlying simplicial structure and identify dynamical patterns via methods such as Q-analysis, persistent homology, topological information and data analysis [14,15,16,17,18,19,20]. Another direction is to investigate the dynamic and stochastic processes of these networks, to uncover the impact of higher-order structures and interactions and to identify the elements of new phenomena. To achieve this end, it is useful to model assemblies of simplicial complexes with specific structural characteristics to examine the functioning of the dynamical processes they can support. The self-assembly of growing structures using cliques is one way of constructing such structures [7,21,22]. Some striking examples of simplicial networks in the real world include (multi)-brain networks originating from human connectomes [23,24,25,26,27,28,29,30] and experimental signals [31,32]. Higher-order interactions based on simplices can appear as constitutive characteristics of complex systems in a wide variety of systems across many contexts, as mentioned above.

Recent studies have focussed on generative models of high-dimensional simplicial complexes of a controlled structure. Self-assembly models of cliques of different orders were developed in Refs. [6,22], where cliques of different sizes can be involved and the assembly is controlled by two factors: geometric compatibility and chemical affinity. The chemical affinity factor controls whether the growing structures will be dense or sparse, leading to a variety of substrate networks and important dynamical effects. The geometric compatibility factor is used to identify the common faces of the existing structure where new simplices can be attached.

2. The Self-Assembly of Simplices

Here, we outline a model for the self-assembly of simplices The generative models employ a bottom-up approach to grow a network with simplicial complexes with a network structure determined by controlling parameters. Here, the self-assembly of cliques is constrained by the geometric compatibility and chemical affinity factors [22]. More precisely, a simplex (clique) of some size n between 2 to 10 is added to an existing structure by sharing one of its faces, i.e., sub-cliques of the size $[eqn]$ $[eqn]$ , with an already existing clique in the growing structure. Here, faces of the order $[eqn]$ correspond to a single node, link, triangle, tetrahedron, etc. The order of the largest sub-clique is $[eqn]$ . At a given point in time, t, a new clique is attached in two steps. First, the order of the face q to be shared with a similar face on the existing structure is determined. Then, the remaining number of new nodes $[eqn]$ are added and all the selected nodes are joined together to form the clique.

The face that is to be shared is identified based on the geometric compatibility of a clique face with the existing structure at the given time. There is also another parameter, the chemical affinity parameter $[eqn]$ , which controls the number $[eqn]$ of new nodes that are added. The probability of a clique of the order $[eqn]$ attaching along the $[eqn]$ face is given by

$$[eqn]$$

where $[eqn]$ is the chemical affinity parameter [22] and $[eqn]$ is the number of geometrically compatible locations available for attaching a simplex of order q on the growing structure at time t. The chemical affinity factor modifies the geometric compatibility via its sign. For a large $[eqn]$ , the system favors the addition of nodes. Thus, the cliques preferably share a single node. For the opposite limit, with increasing $[eqn]$ , the cliques preferably share their larger faces, leading to compact structures. When $[eqn]$ , the process is governed by geometric compatibility alone. Consequently, each face can be shared with a finite probability, depending on its size and the actual network structure. An online demo can be seen in [33]. We note that $[eqn]$ leads to compact structures, $[eqn]$ leads to mixed structures and $[eqn]$ leads to sparse structures. The three types of networks can be seen in Figure 1. Since these have distinctly different substrates, the hysteresis effects observed for each type are different, even when other parameters are identical.

3. Simplicial Structure of Interactions: Higher-Order Interactions

The simplicial structures which are achieved after self-assembly provide useful and natural substrates to support higher-order interactions. Such effects have been analyzed in a variety of contexts, including spin systems and oscillator systems.

Here, we discuss continuous oscillator systems based on networks grown by the generative algorithm described above, using 5-simplex structures as the base geometry. Here, we discuss pairwise and higher interactions. We note here that the nature of the substrate, as well as the chemical affinity parameter, has its own effect on the hysteresis behavior, the transition, in this case the synchronization transition, and the patterns which appear in the system.

3.1. Phase Synchronization with Competing Interactions

We consider an ensemble of N-coupled Kuramoto oscillators associated with the nodes of a given simplicial complex [12,21]. The equation governing the evolution of the phase of i-th oscillator is given by [12]

$$[eqn]$$

where $[eqn]$ s represent the intrinsic frequencies of the phase oscillators. The second and third terms in Equation (2) represent 1-simplex and 2- simplex interactions, respectively. Here, $[eqn]$ is an element of the 1-simplex adjacency matrix and $[eqn]$ if nodes i and j are connected by a link, and $[eqn]$ is an element of the 2-simplex adjacency tensor, with $[eqn]$ , if nodes i, j and l belong to the same triangle, and 0 otherwise. The normalization factors $[eqn]$ and $[eqn]$ indicate the number of links and triangles of the node i.

Thus, equal weight is given to all the terms contributing to the sum of each interaction term. In Equation (2), $[eqn]$ is the intrinsic frequency of the i-th oscillator, which governs its dynamics when there is no interaction with other oscillators in the network. The pairwise interactions seek to reduce the difference between the phase of the i-th oscillator and each of its neighbouring oscillators when $[eqn]$ . In contrast, the oscillators tend towards opposite phases when $[eqn]$ . The third term, representing three-node interactions of the i-th oscillator based on each 2-simplex incident on node i, is a natural generalization of the pairwise interaction term. It is important to note that the interactions between these three nodes occur over the faces of the simplicial complex and not over any given three nodes. Furthermore, this interaction term is symmetric in i in that it is unaffected by permutations in the other two indices.

Here, through synchronization, by which we mean perfect synchronization, all the phase angles at all the considered nodes evolve with identical values. We can also use cluster synchronization, where some of the nodes that are connected to each other remain in step with each other, whereas other nodes are not in step. There can be one or more synchronized clusters in the network.

The degree of synchronization of the whole network is quantified by the usual Kuramoto-order parameter, defined by

$$[eqn]$$

Here, the angular brackets indicate the time average. Hence, $[eqn]$ represents the case of perfect synchronization, i.e., all phases are equal, and $[eqn]$ corresponds to the disordered phase. Meanwhile, the stable states with $[eqn]$ indicate the presence of more complex synchronization patterns, which arise due to the competition between interactions and the effect of shared faces.

3.2. Effects of Interplay

In this context, we examine the emergence and disappearance of collective dynamic behavior, as measured by the order parameter on the hysteresis loop when the relative strengths of the various interactions embedded in the geometry vary. For example, the presence of triangle-based interactions is understood to disrupt the promoted order by increasing the strength of pairwise couplings and can cause abrupt desynchronization [7,34,35]. Furthermore, the occurrence of partially synchronized phases with negative pair interactions is another striking feature of geometric interactions on simplicial complexes; see, for example, [12] and references therein.

The influence of topology on diffusive processes has been examined, e.g., via a spectral analysis of the network [36,37,38,39,40,41,42,43]. One possible measure is the spectral dimension as derived from the eigenvalue spectrum of the Laplace operator associated with the network adjacency matrix. In the context of phase synchronization, it has been observed that networks with $[eqn]$ are required to enable stable global synchronization, while such states cannot be reached when $[eqn]$ . There has been some analysis of the conditions for global synchronization [8,44]. However, the origin of partial synchronization on simplicial complexes and the nature of the underlying dynamical states needs further analysis. In addition to the topological dimension of a simplicial complex, the role of its architecture in the synchronization processes remains an open problem, especially in relation to by the presence of higher-order interactions and the internal inhomogeneity of the nodes’ dynamics.

An important physical consequence of the manner in which the extent to which these factors influence the transition to synchronization and desynchronization can be inferred is the behavior of the hysteresis loop. As seen above, $[eqn]$ and $[eqn]$ denote the strength of the pair and triangle interactions. The hysteresis is studied via the behavior of the order parameter, which is a measure of the level of synchronization in the system of increasing strengths in the pairwise coupling (forward sweep), as well as decreasing ones (reverse sweep) for fixed values of the triplet coupling. We examined both positive and negative values of $[eqn]$ .

We examined the response of the system via the hysteresis loop of the system for different simplicial geometries, and for varying contributions from pairwise and triangular interactions. These factors have significant effects. The structure of the networks grown here is controlled via the growing algorithm, which incorporates the geometric compatibility and the chemical affinity, which essentially control the number of neighbors of each node, leading to varying effects. These include partial or frustrated synchronization and the presence of fluctuations and meta-stable states. While many aspects of the observed patterns are under analysis, we discuss the current status in the next section.

3.3. Hysteresis Effects

We start with the simplest case: the synchronization process using only the pairwise interactions. These interactions are enabled by the edges of the simplicial complex underlying the network via the connectivity matrix. Hence, different network features from the local to global levels are expected to play a role in the cooperative behaviors, depending on the interaction strength $[eqn]$ .

For the simulations here, the initial condition for $[eqn]$ was chosen from the uniform distribution $[eqn]$ and the frequencies $[eqn]$ at each node i were chosen from a the Gaussian distribution with a zero mean, respectively. The numerical solution of the set of equations was obtained using the Python library SciPy (version 1.0 (2020)), with 50,000 iteration steps, always considering the previous state of each dynamical variable (parallel updates).

For the four-dimensional simplicial complex, as $[eqn]$ increased from zero to $[eqn]$ , we observed a continuous transition from a desynchronized state where $[eqn]$ to a completely synchronized state $[eqn]$ ; see Figure 2. On the other hand, as $[eqn]$ decreased from zero to $[eqn]$ , a partially synchronized state ( $[eqn]$ ) was reached. The same behavior was seen for both increasing and decreasing increments in $[eqn]$ . Thus, there was no hysteresis. Randomized versions of the network show a fully desynchronized regime for $[eqn]$ .

The introduction of two-simplex interactions changes the situation drastically. A hysteresis loop now opens up due to the presence of triangle-based interactions. The hysteresis loop grows in size as $[eqn]$ is increased, affecting both the forward sweep on the positive $[eqn]$ side and the size of the first-order jump. Similar behavior is seen for a wider range of values of the 2-simplex interaction strength as long as the large positive pairwise interactions are sufficient to maintain complete synchrony. However, for larger values of $[eqn]$ , the fully synchronized state is no longer accessible; instead, partial synchronization is reached in the presence of the competing interactions. The detailed behavior that is observed depends on the sign and strength of the pairwise interactions. These features can be seen in Figure 2 and Figure 3.

The partial synchronization seen in the system can be attributed to frustration effects [43]. It is necessary to analyze small lattices to identify the exact frustration effects and patterns. However, cluster synchronization can clearly be discerned via the distribution of phases for these large lattices. We note that such effects have been observed in other systems and analyzed via eigenvalue analysis [45].

It should also be noted that if all oscillators rotate with a single common frequency $[eqn]$ , then discontinuous transitions and hysteresis behavior can be seen, even for the pure pairwise interaction case (where $[eqn]$ ) [21]. Thus, higher-order triangle-based interactions are not necessary to see a first-order transition for this case (see Figure 4).

4. Patterns and Dynamics of Individual Nodes

The hysteretic behavior discussed above was explored via the global-order parameter r. It is also instructive to examine the phase angle behavior at individual nodes, as well as the time evolution of the average angle, $[eqn]$ , and the distribution of phase angles at a fixed time t. These can be seen for specific values of $[eqn]$ in Figure 5, which were selected from the forward and backward branch of the hysteresis loops for the $[eqn]$ values $[eqn]$ , $[eqn]$ and $[eqn]$ .

1. $[eqn]$ , compact networks: Here, we choose the values $[eqn]$ on the forward branch with a pure pairwise interaction, f and $[eqn]$ , on the backward branch. The behavior of the three quantities of interest, $[eqn]$ , $[eqn]$ and the distribution of $[eqn]$ , can be seen in panels $[eqn]$ of Figure 5. It can be seen from panels a1 and a3 that the $[eqn]$ has one prominent peak value in the distribution, with a range of $[eqn]$ values in the rest of the distribution. The average angle in the time evolution from 490 to 500 steps shows one clear oscillation and then flattens off.At $[eqn]$ , $[eqn]$ , we examine the first 500 steps and can hence see the transient before the iterates settle. Panel $[eqn]$ shows that the $[eqn]$ tends to synchronize. The average angle $[eqn]$ shows small transient values due to the spread of phase angles and then shows a clear oscillation. The distribution of phase angles peaks around a single value.At $[eqn]$ , $[eqn]$ , the $[eqn]$ -s shows a wide spread (see panels $[eqn]$ and $[eqn]$ ). The average angle shows a reasonably regular oscillation between time steps 490 and 1000 (panel $[eqn]$ ).2. $[eqn]$ , sparse networks: The behavior of the sparse network is shown in the $[eqn]$ panels. The d and e panels show the behavior of the pairwise interactions with $[eqn]$ on the forward and backward branch alone for a negative value of K. The f panels show the effects of the higher-order term alone, i.e., $[eqn]$ , $[eqn]$ , forward branch. The pairwise backward branch shows a banded structure in the $[eqn]$ plot, which translates into oscillations in the $[eqn]$ , as seen in $[eqn]$ . The forward branch shows a more or less constant value of the $[eqn]$ . The higher-order interactions also show clear oscillations in panel $[eqn]$ . Hence, peaks are seen in the distribution in $[eqn]$ and $[eqn]$ . However the higher-order interactions show a more pronounced peak.We note that the clustering is much more clear in the compact structure. The sparse structure tends to distribute the phase angles more uniformly across the available range, although the higher-order interactions promote some clustering. Thus, both the kind of interaction and the substrate influence the synchronization and clustering seen in this system.

The networks tend to not synchronize if only a single node is shared between the simplices and show clustered behavior as higher q structures, such as links, triangles and tetrahedra are shared between the structures. At $[eqn]$ , the structures show full synchronization, as seen from the hysteresis loops. Switching on the $[eqn]$ interactions shows that $[eqn]$ , $[eqn]$ and $[eqn]$ have similar behavior, with one node shared between the simplices. Sharing a triangle leads to two clusters for $[eqn]$ and both $[eqn]$ = −1.0 and 0.0. For $[eqn]$ and $[eqn]$ and $[eqn]$ , with a common link, show clustered behavior again, with the number of nodes that are synchronized increasing with an increase in the strength of pairwise interactions. Common tetrahedra with $[eqn]$ and $[eqn]$ and $[eqn]$ show two clustered behaviors.

Partial Synchrony and Patterns

The partial synchrony seen here is a consequence of cluster formation, as can be confirmed through the behavior of larger lattices. The number of clusters that is formed decreases with the increase in the number of nodes in the common faces, as might be expected. The more compact structures have a greater tendency to synchronize globally, in contrast with the sparse structures. Negative pairwise interactions tend to inhibit or ‘frustrate’ the synchronization by pushing the phase angles apart. We will now discuss this in further detail.

We note that some studies of Kuramoto oscillators on regular two-dimensional lattices showed frustration effects [46] for the case of repulsive pairwise coupling. Here, at higher coupling parameter values, the oscillators relaxed into a state that was $[eqn]$ out of phase with its neighbor. Such states are called splay states in the literature, and were observed in the context of defects [46] in oscillator systems, as well as in coupled Sine circle map lattices [47]. We note that such states lead to partial synchronization in the system.

5. Small Systems: Two Simplex Systems

Small simplex systems constitute the simplest case in which the patterns can be analyzed. Here, we analyze the cases of two 5- simplex cliques for the cases where the two simplices share a node, a link, a triangle and a tetrahedron. The hysteretic behavior seen here is shown in Figure 6. The intrinsic frequency was $[eqn]$ at all nodes. We noted hysteretic behavior in all cases and also stable behavior once full synchronization was reached, despite the decreasing coupling, regardless of whether pairwise or triplet interactions occurred. A partially synchronized case was observed, even in the single simplex case, for $[eqn]$ . For the two simplex cases, the nine node cases (one node shared) showed partial synchronization everywhere. Desynchronization was only observed at two points in the hysteretic loops for $[eqn]$ . The triangle shared case (eight nodes) showed partial synchronization in all cases, with a stable step for $[eqn]$ . The case with shared tetrahedra showed a gentler increase to full synchronization.

The detailed evolution at each node for the two simplex case is shown in Figure 7, Figure 8 and Figure 9 for the $[eqn]$ and $[eqn]$ values indicated by crosses ( $[eqn]$ ), stars and diamonds in Figure 6:

• 1.Crosses ( $[eqn]$ ): Here the node with two simplices shows maximum desynchronization. The link that shared two simplices shows four clusters; the triangle that shared two simplices reduced to two clustes; the tetrahedron case showed four clusters again. The ( $[eqn]$ ) case shows full synchronization.
• 2.The cases denoted by stars ( $[eqn]$ and $[eqn]$ ) show several clusters for the node shared case and two clusters can be observed for the shared $[eqn]$ cases denoted by a triangle.
• 3.Diamonds: ( $[eqn]$ , and $[eqn]$ case, link shared). The first two cases are multiclustered, while the last case show two clustered solutions. The shared tetrahedra cases marked by diamonds ( $[eqn]$ and $[eqn]$ cases show three oscillating cluster solutions.

We note that the backward sweep is very stable, and once synchronization is achieved, the two clusters are very stable. It may be interesting to see what happens if the frequencies are no longer uniform for all nodes and are locked to specific values.

Patterns of the Two Simplex Structure

The patterns shown in the two-simplex case can be deciphered more easily via a space–time plot. These were plotted between iterates 400 and 500 after discarding a transient of 400 iterates, where the patterns had stabilized. Following the same sequence as above, the following can be obtained:

• 1.Crosses ( $[eqn]$ ): The patterns here can be seen in the first row of Figure 10. The figure shows the behavior at each node for two simplex situations sharing nodes, links, triangles and tetrahedra. The lead and lag patterns can be seen in the space–time plot. The phase angles at nodes labeled 1 to 4 follow node 4 which leads, while node 5 lags and node 6 is synchronized with node 4, which is ahead of nodes 7 to 10. For the shared link case, nodes 1 and 5 are synchronized and the remaining nodes are desynchronized. The triangle shared case shows three synchronized clusters of two nodes and a single node, none of which are in step with each other. The shared tetrahedra can be fully synchronized ( $[eqn]$ ) or fully desynchronized ( $[eqn]$ ).
• 2.Stars ( $[eqn]$ , and $[eqn]$ ): The shared node cases can be seen in the second row of the space–time plot. As mentioned above, the $[eqn]$ node linked case is multi-clustered, with nodes 4 and 7 being synchronized with each other. The remaining nodes show no systematic behavior. The $[eqn]$ case shows three clusters: nodes $[eqn]$ are synchronized nodes, 1 and 7 are synchronized and nodes 9 and 8 are synchronized. The $[eqn]$ , $[eqn]$ node linked case is similar. The two shared cases denoted by triangles are seen in the same row. Here, of the $[eqn]$ nodes $[eqn]$ are in step, 1 and 4 are in step and node 5 evolves independently. For $[eqn]$ , no clear pattern is seen though there are two clusters.
• 3.Diamonds: ( $[eqn]$ , and $[eqn]$ case, link shared). This is the last row in the space–time plot. Here, there are three link shared cases and two tetrahedra shared cases. For $[eqn]$ =0.0, $[eqn]$ , (link shared) node 5 leads all the other nodes. For $[eqn]$ and $[eqn]$ , nodes $[eqn]$ are in step, as are nodes $[eqn]$ and $[eqn]$ . For links $[eqn]$ , $[eqn]$ , there are two clearly synchronized clusters of nodes $[eqn]$ and nodes $[eqn]$ . For the shared tetrahedra, we can observed a synchronized cluster where $[eqn]$ nodes are in step, $[eqn]$ are in step and $[eqn]$ are in step (three clusters). For the last case, again, there are three clusters in step: $[eqn]$ , $[eqn]$ and $[eqn]$ .

The observed patterns are due to cluster synchronization and lead to partially synchronized values for the ordered parameter. These appear to be stable, at least during the period of analysis. Further eigenvalue analysis needs to be carried out to identify the bifurcation behavior that leads to these.

6. Summary and Conclusions

We explored two aspects of the dynamics of simplicial structures: the role of the simplicial structure of the substrate and the role of the higher-order interactions. Some recent studies on higher-order interactions include [48,49]. We note that both features make significant contributions to the qualitative features of the dynamics and the transition to synchronization. The substrate structure was controlled via a growing algorithm which incorporates geometric compatibility and chemical structure, leading to compact, mixed and sparse structures. We note that compact structures have a greater tendency to synchronize and to show abrupt synchronization. This also depends on the distribution of the frequencies of the oscillators. We also note that the paths to synchronization are achieved via cluster synchronization, with the exact paths being dependent on the structure of the lattice. Some structures show partial synchronization due to the presence of specific dynamical states. This can depend on the eigenvalue distribution. Similar behavior can be seen in spin systems.

The hysteresis phenomenon allows for a very stable response to synchronization in the reverse phase. In the context of phase locked systems, this could permit very stable synchronization and frequency-locking on the reverse phase, even if individual oscillators change frequencies.

In the context of small systems, we note the phase shifts between individual oscillators and that the system geometries control how shared nodes lead to synchronization between subsets of oscillators with different node labels. This could be put to use for message encryptation purposes.

Noise effects can be very important in these systems, and have not been widely explored. Individual nodes can be more or less stable to noise, depending on the number of shared neighbours. Intermittency effects are possible in such systems.

Two-cluster and single-cluster configurations can also be utilized to switch between frequencies, contributing to resonance effects. Further investigation could lead to quite useful applications. Our other results may have relevance to studies of nanostructures and the human connectome.