Necessary condition for the synchronization of dynamical systems with higher-order interactions

The object of our study is the most general simplicial complex of N coupled dynamical units. This means that the different dynamical units are subject not only to pairwise interactions, but also to three-body interactions, four-body interactions and so on. The precise microscopic structure of the interactions is described by underlying simplicial complex, which can have any dimension D ≥ 1 (for all details and notations on simplicial complexes, see the “Methods”). In the particular case of D = 1, our system coincides with the standard case of a complex network of N coupled dynamical units. We assume that the equations of motion governing the dynamics of our D-dimensional simplicial complex can be written as:

As the notation may result somehow cumbersome, for the sake of clarity in what follows we illustrate the case of D = 2, so that a reader will be able to appreciate each and every conceptual action we are making. At the end, we will then summarize the steps one has to do in order to extrapolate the results to all values of D.

Let us then consider the following set of coupled differential equations

Existence and invariance of the synchronized solution x^s(t) = x₁(t) = … = x_N(t) are guaranteed by the noninvasiveness of the coupling functions. In order to study the stability of the synchronization solution, one considers small perturbations around the synchronous state, i.e., δx_i = x_i − x^s, and perform a linear stability analysis of Eq. (2). To do this, one can perform the following transformation of the variables. Consider $$ \delta \mathbf{x}={\left[\delta {\mathbf{x}}_1^T,\delta {\mathbf{x}}_2^T,\dots ,\delta {\mathbf{x}}_N^T\right]}^T$$, and let us take, as a reference basis of $$ {\mathbb{R}}^{}$$, the one made by the eigenvectors v₁, v₂, …, v_N of the classic Laplacian matrix $$ {\mathcal{L}}^{\left(1\right)}$$. This allows to define new variables η = (V⁻¹ ⊗ I_m)δx, where V = [v₁, v₂, …, v_N]. To express the dynamics of the system in terms of the new variables η, one needs to extend the notion of classical Laplacian matrix, which accounts for pairwise interactions, to a set of generalized Laplacian matrices, where the generic matrix of order d, indicated as $$ {\mathcal{L}}^{\left(d\right)}$$, accounts for (d + 1)-body interactions (for a formal definition see “Methods”). In the specific case of D = 2, we will therefore describe the systems with two matrices, $$ {\mathcal{L}}^{\left(1\right)}$$ and $$ {\mathcal{L}}^{\left(2\right)}$$, respectively.

Through a series of three conceptual steps detailed in the “Methods”, the following equations can be derived

The problem of stability is then reduced to: (i) simulating a single, uncoupled, nonlinear system; (ii) using the obtained trajectory to feed up the elements of the Jacobians JG⁽¹⁾ and JG⁽²⁾; (iii) simulating the dynamics of a system of N − 1 coupled linear equations, and tracking the behavior of the norm $$ \sqrt{{\sum }_{i=2}^N{\sum }_{j=1}^m{\left({\eta }_i^{\left(j\right)}\right)}^2}$$ for the calculation of the maximum Lyapunov exponent (being $$ {\eta }_i\equiv \left({\eta }_i^{\left(1\right)},{\eta }_i^{\left(2\right)},...,{\eta }_i^{\left(m\right)}\right)$$).

Stability of the synchronous solution requires as a necessary condition that Λ_max, the maximum among the (conditional) Lyapunov exponents associated to all transverse modes, be negative. Given the node dynamics and the coupling functions, Λ_max is in general function of the topology of the two-body interactions, the topology of the three body interactions, and the two coupling strengths σ₁ and σ₂, i.e., Λ_max = $$ {\mathrm{\Lambda }}_{\max }\left({\sigma }_1,{\sigma }_2,{\mathcal{L}}^{\left(1\right)},{\mathcal{L}}^{\left(2\right)}\right)$$.

It is important to notice that, in analogy with the classical MSF approach, also in the case of simplicial complexes one is, therefore, able to separate the motion along the synchronization manifold and that transverse to it, and such a crucial separation ultimately enables the study of stability of the synchronous manifold. For simplicial complexes, however, the higher complexity in the structure of the interactions yields a formalism requiring the analysis of a set of coupled differential equations, rather than of a single parametric variational equation (as in the case of the MSF). In other words, in the fully general case the set of equations describing the motion transverse to the synchronous manifold cannot be further decomposed into independent, decoupled modes, as it happens in the network case; however, the analysis of stability still requires a straightforward computation of a single quantity, i.e., the maximum Lyapunov exponent, which has to be performed on such a set of coupled, linear equations. In the more general case, the transverse modes are intertwined, such that stability has to be analyzed without reduction in dimensionality. However, we will momentarily show that there are relevant instances where such an expression can be simplified, up to recover a formalism that is identical to the classical MSF, allowing separation of the modes and reduction of the dimensionality of the problem to a single parameteric variational equation.

In analogy with the classification of systems made for synchronization of complex networks (Chapter 5 in ref. ¹), one immediately realizes that, once specified the dynamical system taking place in each node (i.e., the function f), the various coupling functions g⁽¹⁾ and g⁽²⁾, and the structure of the simplicial complex (i.e., $$ {\mathcal{L}}^{\left(1\right)}$$ and $$ {\mathcal{L}}^{\left(2\right)}$$), all possible cases can be divided in three classes: (i) class I problems, where Λ_max is positive in all the half plane (σ₁ ≥ 0, σ₂ ≥ 0), and therefore synchronization is never stable; (ii) class II problems, for which Λ_max is negative within a unbounded area of the half plane (σ₁ ≥ 0, σ₂ ≥ 0); and (iii) class III problems, for which the area of the half plane (σ₁ ≥ 0, σ₂ ≥ 0) in which Λ_max is negative is instead bounded, and therefore additional instabilities of the synchronous motion may occur at larger values of the coupling strengths. While class I problems are trivial (in that synchronization is never observed), examples of class II and class III problems are shown in Fig. 1 for simplicial complexes of Rössler oscillators⁶², and one easily sees that the predictions made by solving Eq. (3) are indeed fully confirmed by the simulations of the original system in Eq. (2).

Fig. 1

Synchronization in simplicial complexes of Rössler oscillators.

Contour plots of the time averaged (over an observation time T = 500) synchronization error E (see “Methods” for definition and the vertical bars of each panel for the color code) in the plane (σ₁, σ₂) for some examples of simplicial complexes (whose sketches are reported in the top left of each panel). Simulations refer to coupled Rössler oscillators (x = (x, y, z)^T and f = (−y−z, x+ay, b+z(x−c))^T) with parameters fixed in the chaotic regime (a = b = 0.2, c = 9). In a–d, $$ {\mathbf{g}}^{\left(1\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)={\left[x_j-x_i,0,0\right]}^T$$, while in (e) $$ {\mathbf{g}}^{\left(1\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)={\left[0,y_j-y_i,0\right]}^T$$. As for the other coupling function, one has $$ {\mathbf{g}}^{\left(2\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,y_j^2{y}_k-y_i^3,0\right]}^T$$ in (d) and $$ {\mathbf{g}}^{\left(2\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[x_j^2{x}_k-x_i^3,0,0\right]}^T$$ in all other panels. The blue continuous lines are the theoretical predictions of the synchronization thresholds obtained from Eq. (3). a, b, and c are examples of class III problems, whereas panels d and e are examples of class II problems.

Far from being limited to the case of D = 2, our approach can be extended straightforwardly to simplicial complexes of any order D. Each term on the right hand side of Eq. (1) can, indeed, be manipulated following exactly the same three conceptual steps described in the “Methods”. Once again, one is entitled to select the eigenvector set which diagonalizes $$ {\mathcal{L}}^{\left(1\right)}$$, to introduce the new variables η = (V⁻¹ ⊗ I_m)δx. Following the very same steps which led us to write Eq. (3), one then obtains

The MSF for synchronization in simplicial complexes

Our results can be greatly simplified in a series of relevant cases in which either the topology of the connectivity structure, or the coupling functions, allow to formulate our approach in terms of MSF. Once again, for the sake of illustration, we will start considering first the case of D = 2, and then the extension to any order D.

The first case is an all-to-all coupling, for which every two and three-body interaction is active. In this case, the classical Laplacian matrix is.

Consequently, one has that

Starting from Eq. (2), applying the steps detailed in the Methods and noticing that in the all-to-all configuration λ₂ = …λ_N = N, for each η_i (with $$ i\in \left\{2,\dots ,N\right\}$$), one obtains

In other words, in the all-to-all case, the variables η_i come out to be all uncoupled to each other, so that Λ_max uniquely depends on σ₁, σ₂ and N, i.e., Λ_max = Λ_max(σ₁, σ₂, N).

In the more general case of a D-dimensional simplicial complex, it is easy to write the generalized Laplacian of order d as a function of the classical Laplacian matrix. In fact, the number of d-simplices having node i as a vertex and the number of d-simplices formed by the link (i, j) are respectively

Given the definition of the generalized Laplacian, we find that

Once again, one can derive a parametric equation analogous to Eq. (8), with a MSF (once fixed both the node dynamics and the coupling functions) which solely depends on the coupling coefficients and the number of nodes, i.e. Λ_max = Λ_max(σ₁, σ₂, …, σ_D, N)

Another interesting case is that of generalized diffusion interactions with natural coupling functions. This amounts to consider diffusive coupling functions, given by

where $$ {\mathbf{h}}^{\left(1\right)}:{\mathbb{R}}^m\to {\mathbb{R}}^m$$ and $$ {\mathbf{h}}^{\left(2\right)}:{\mathbb{R}}^{2m}\to {\mathbb{R}}^m$$. In addition, a condition of natural coupling is considered:

Equation (14) expresses, indeed, the fact that the coupling to node i from two-body and three-body interactions is essentially similar, in that a three-body interaction where two nodes are on the same state is equivalent to a two-body interaction. Here, our approach takes the form of a MSF with a particularly convenient expression, as it can be written as a function of a single parameter. In fact, in this case, the transverse modes can be fully decoupled (see the “Methods” for the full derivation) and a single parameter MSF can be defined, starting from the following m-dimensional linear parametric variational equation

from which the maximum Lyapunov exponent is calculated: Λ_max = Λ_max(α) with $$ \alpha =\lambda \left({\sigma }_1{\mathcal{L}}^{\left(1\right)}+{\sigma }_2{\mathcal{L}}^{\left(2\right)}\right)$$ or $$ \alpha ={\sigma }_1\lambda \left({\mathcal{L}}^{\left(1\right)}+r{\mathcal{L}}^{\left(2\right)}\right)={\sigma }_1\lambda \left(\mathcal{M}\right)$$, where $$ \mathcal{M}$$ is given by $$ {\mathcal{L}}^{\left(1\right)}+r{\mathcal{L}}^{\left(2\right)}$$ with $$ r=\frac{{\sigma }_2}{{\sigma }_1}$$. The situation, is therefore, conceptually equivalent to that of synchronization in complex networks, with the effective matrix $$ \mathcal{M}$$ playing the same role of the classical Laplacian: given the dynamical system f, the coupling functions h⁽¹⁾ and h⁽²⁾, and the structure of connection of the simplicial complex (i.e., $$ {\mathcal{L}}^{\left(1\right)}$$ and $$ {\mathcal{L}}^{\left(2\right)}$$) one can define three possible classes of problems:

(i)class I problems, for which the curve Λ_max = Λ_max(α) does not intercept the abscissa and it is always positive. In this case synchronization is always forbidden, no matter which simplicial complex is used for connecting the dynamical systems;

(ii)class II problems, for which the curve Λ_max = Λ_max(α) intercepts the abscissa only once at α_c, and for which, therefore, the synchronization threshold is given by the self consistent equation $$ {\sigma }_1^{\mathrm{critical}}={\alpha }_c/{\lambda }_2\left[\mathcal{M}\left({\sigma }_1^{\mathrm{critical}},{\sigma }_2^{\mathrm{critical}}\right)\right]$$, i.e. it scales with the inverse of the second smallest eigenvalue of the effective matrix;

(iii)class III problems, for which the curve Λ_max = Λ_max(α) intercepts the abscissa twice at α₁ and α₂ > α₁. In this case, synchronization can be observed only if the entire eigenvalue spectrum of the effective matrix is such that $$ {\sigma }_1{\lambda }_2\left(\mathcal{M}\right)>{\alpha }_1$$ and, at the same time, $$ {\sigma }_1{\lambda }_N\left(\mathcal{M}\right)<{\alpha }_2$$. In this case, the parameter $$ \frac{{\lambda }_2\left(\mathcal{M}\right)}{{\lambda }_N\left(\mathcal{M}\right)}$$ can be considered as a proxy measure of synchronizability of the simplicial complex, in that the closer is such a parameter to unity (the more compact is the spectrum of eigenvalue of $$ \mathcal{M}$$) the larger can be the range of coupling strengths for which the two above synchronization conditions can be satisfied.

We have so far considered the case of D = 2. In the fully general scenario, the condition for natural coupling is given by

The equation for the MSF is formally analogous to Eq. (15), where now $$ \alpha ={\sigma }_1\lambda \left({\mathcal{M}}^{\left(D\right)}\right)$$ parameterizes the eigenvalues of the effective matrix of order D

In summary, we have shown that, while in the general case the transverse modes are intertwined, in the case of all-to-all coupling or of natural coupling functions a significant dimensionality reduction of the stability analysis problem is obtained, through the formulation of a MSF (Eq. (12) for all-to-all coupling and Eq. (15) for natural coupling).

Synchronization in simplicial complexes of chaotic systems

Following is a series of results confirming the validity and wide applicability of our approach. We focus on two paradigmatic three-dimensional ($$ \mathbf{x}={\left(x,y,z\right)}^T\in {\mathbb{R}}^3$$) chaotic systems, namely the Rössler oscillator⁶² and the Lorenz system⁶³, and, as a real-world example of neuron dynamics, on the Hindmarsh-Rose (HR) model⁶⁴ (see the “Methods” for the equations describing the three systems, as well as for the setting of parameters and of stipulations for the numerical simulations). In particular, we start with considering the more general case with diffusive coupling, then we discuss our results on neuron dynamics and on the MSF cases of all-to-all and natural coupling, where we also show an analysis carried out on a real-world structure. Finally, we move away from the study of complete synchronization and illustrate an example of cluster synchronization in simplicial complexes.

The general case

Our discussion begins with going back to Fig. 1, where we have considered a few elementary configurations of simplicial complexes, chosen in order to illustrate the classes of problems that one can deal with even when the structures involve only a small number of nodes. In particular, Fig. 1 reveals that synchronization in the general case crucially depends on the topology and the coupling functions: the same configuration can in fact feature different dynamics when diverse mechanisms regulate the coupling and, conversely, the same coupling functions may lead to different behaviors when the topology of interactions changes.

As an example, let us consider the full dynamical equations of coupled Rössler oscillators, when the coupling functions are chosen as $$ {\mathbf{g}}^{\left(1\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)={\left[x_j-x_i,0,0\right]}^T$$ and $$ {\mathbf{g}}^{\left(2\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[x_j^2{x}_k-x_i^3,0,0\right]}^T$$. They read

In each of the configurations considered, the state of the system is monitored by the average synchronization error E defined in the “Methods”. Figure 1 reports E(σ₁, σ₂) for different simplicial complexes (shown as insets in the panels) and coupling functions, along with the theoretical predictions provided by Eq. (3) (the blue, continuous, lines superimposed to the diagrams of the synchronization error). In all the cases, the numerical simulations are in very good agreement with the theoretical predictions for the synchronization thresholds.

The results of Fig. 1 suggest several interesting considerations. Indeed, in the cases reported in panels a and b of Fig. 1 synchronization may be achieved using either two-body or three-body interactions only (for very small σ₁ indeed there is a range of values of σ₂ leading to synchronization, and viceversa), while in the case of panel c synchronization is forbidden for very small values of σ₁. In the last case, in fact, the two triangles do not have a common edge as in Fig. 1a, nor a common node as in Fig. 1b, and therefore interactions through links becomes essential for synchronization. Finally, one notice that there are scenarios, as in panels d and e, where the synchronization region is unbounded. As already mentioned, Fig. 1 provides examples of two of the three possible classes of behavior, with class III behavior in Fig. 1a–c, and class II in Fig. 1d (where synchronization exists in an unbounded region of the coupling coefficient regulating pairwise interactions, i.e., σ₁) and in Fig. 1e (where synchronization exists in an unbounded region of the coupling coefficient regulating three-body interactions, i.e., σ₂).

Applications to neuron dynamics

We now discuss the applicability of our framework to the study of neuron dynamics. Synchronization in neural activity is of utmost importance. On the one hand, synchronized network oscillations are known to play a role in establishing ensembles of neurons in a task-dependent, flexible manner⁴⁰. On the other hand, synchronization of neural activity is associated to epileptic seizures^41–43. Recent evidences in neuroscience have pointed out the existence of higher-order interactions between neurons⁴⁴. In particular, astrocytes and other glial cells are considered a plausible biological source of high-order interactions^45,46, as they make contact with thousands of synapses and actively modulate their function⁶⁵. Although how to account for these interactions in nonlinear models of neuron dynamics is still an open problem, here we discuss an example showing the suitability of our framework to study neuronal synchronization in the presence of high-order coupling. As a pratical case study, we here consider an ensemble of HR neurons, subject not only to pairwise coupling but also to three-body interactions. The system is described by

Fig. 2

Synchronization in simplicial complexeses of Hindmarsh–Rose neurons.

Contour plots of the time averaged (over an observation time T = 500) synchronization error E (see “Methods” for definition and the vertical bars of each panel for the color code) in the plane (σ₁, σ₂) for simplicial complexes of HR neurons coupled as in Eq. (19). Parameters are fixed in the chaotic regime (r = 0.006, s = 4, I = 3.2). a–c refer to three different simplicial complexes corresponding to the structures considered in Fig. 1a–c. The blue continuous lines are the theoretical predictions of the synchronization thresholds obtained from Eq. (3).

MSF cases

Let us now move to discuss other results, which refer to the cases where our approach yields a MSF. We start with the all-to-all coupling case where, according to Eq. (8), one obtains a MSF that is function of N, σ₁ and σ₂. We then consider a simplicial complex of Rössler oscillators with all-to-all coupling, described by

The results are shown in Fig. 3 for three values of N (N = 10, N = 50, and N = 100): the synchronous manifold is stable in a bounded region of the semiplane (σ₁ > 0, σ₂ > 0) delimited by blue (N = 10), red (N = 50) and black (N = 100) lines. One immediately sees that such a stability region moves toward the origin when N is increased. Hence, increasing N reduces the lower and upper thresholds for achieving synchronization.

Fig. 3

Synchronization in a simplicial complex of Rössler oscillators with all-to-all coupling.

Lower and upper boundary curves for the region where synchronization is stable, at different values of N. The color codes for the different curves is reported at the top of the panel.

Finally, we consider the case of natural coupling. Here, in full analogy with what occurs for networks, the MSF is a function of a single parameter, i.e., Λ_max = Λ_max(α) with $$ \alpha =\lambda \left({\sigma }_1{\mathcal{L}}^{\left(1\right)}+{\sigma }_2{\mathcal{L}}^{\left(2\right)}\right)$$ or $$ \alpha ={\sigma }_1\lambda \left({\mathcal{L}}^{\left(1\right)}+r{\mathcal{L}}^{\left(2\right)}\right)={\sigma }_1\lambda \left(\mathcal{M}\right)$$. This enables the study of synchronization stability into two steps, one pertaining only to the node dynamics and coupling functions, providing Λ_max = Λ_max(α), and a second step, where the condition Λ_max(α) < 0 is checked at the points $$ \alpha =\left\{{\sigma }_1{\lambda }_2\left(\mathcal{M}\right),\dots ,{\sigma }_1{\lambda }_N\left(\mathcal{M}\right)\right\}$$.

We calculated the MSF for the Rössler oscillator and the Lorenz system with several choices of the coupling functions: $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[x_j^3,0,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[x_j^2{x}_k,0,0\right]}^T$$; $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,x_j^3,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,x_j^2{x}_k,0\right]}^T$$; $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,0,x_j^3\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,0,x_j^2{x}_k\right]}^T$$; $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[y_j^3,0,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[y_j^2{y}_k,0,0\right]}^T$$... $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,0,z_j^3\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,0,z_j^2{z}_k\right]}^T$$.

The results are shown in Fig. 4 for the Rössler oscillator and in Fig. 5 for the Lorenz system. Both cases exhibit a variety of behaviors that actually encompass all possible classes of MSF. In the case of Rössler oscillator we have one class III example (Fig. 4a), one class II example (Fig. 4e), while all remaining cases do correspond to class I. In the case of the Lorenz system we have several examples of class I behavior (Fig. 5c, f–h); three class II examples (Fig. 4a, d, e), and one class III example with a very narrow region for synchronization (Fig. 4b). Moreover, in Fig. 4i the MSF assumes negative values in two different intervals of α; overall, this represents a further example of class III behavior, providing however the extra scenario where increasing the coupling strength one can achieve alternating regions of synchronization and desynchronization.

Fig. 4

Synchronization in simplicial complexes of Rössler oscillators, in the case of natural coupling.

The Master Stability Function is here calculated taking into account several coupling functions. a $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[x_j^3,0,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[x_j^2{x}_k,0,0\right]}^T$$, b $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[y_j^3,0,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[y_j^2{y}_k,0,0\right]}^T$$, c $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[z_j^3,0,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[z_j^2{z}_k,0,0\right]}^T$$, d $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,x_j^3,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,x_j^2{x}_k,0\right]}^T$$, e $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,y_j^3,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,y_j^2{y}_k,0\right]}^T$$, f $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,z_j^3,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,z_j^2{z}_k,0\right]}^T$$, g $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,0,x_j^3\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,0,x_j^2{x}_k\right]}^T$$, h $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,0,y_j^3\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,0,y_j^2{y}_k\right]}^T$$, i $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,0,z_j^3\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,0,z_j^2{z}_k\right]}^T$$.

Fig. 5

Synchronization in simplicial complexes of Lorenz systems, in the case of natural coupling.

The Master Stability Function is here calculated taking into account several coupling functions. a $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[x_j^3,0,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[x_j^2{x}_k,0,0\right]}^T$$, b $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[y_j^3,0,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[y_j^2{y}_k,0,0\right]}^T$$, c $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[z_j^3,0,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[z_j^2{z}_k,0,0\right]}^T$$, d $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,x_j^3,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,x_j^2{x}_k,0\right]}^T$$, e $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,y_j^3,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,y_j^2{y}_k,0\right]}^T$$, f $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,z_j^3,0\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,z_j^2{z}_k,0\right]}^T$$, g $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,0,x_j^3\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,0,x_j^2{x}_k\right]}^T$$, h $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,0,y_j^3\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,0,y_j^2{y}_k\right]}^T$$, i $$ {\mathbf{h}}^{\left(1\right)}\left({\mathbf{x}}_j\right)={\left[0,0,z_j^3\right]}^T$$ and $$ {\mathbf{h}}^{\left(2\right)}\left({\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[0,0,z_j^2{z}_k\right]}^T$$.

Real-world structures

As an example of a real-world structure, we apply our method to a social system modeling the interactions between the members of a university sport club, the so-called Zachary karate club data set⁶⁶. The original social system is described in terms of a network consisting of N = 34 nodes and 78 links. Since the links form 45 triangles, several simplicial complexes can be constructed from this network, depending on which and how many nodes forming a triangle are effectively taken into consideration as being part of a 2-simplex or, on the contrary, as only connected by three pairwise interactions. In this way, we will be able to investigate the relevance of three-body interactions in mechanisms of collective behavior, such as consensus and synchronization, in social systems. It is indeed well known that pairwise interactions are not always enough to capture the complex behavior of many systems, including social systems⁴⁴. For instance, processes of social contagion can occur in different ways, either through pairwise interactions (the links of a network), or in groups of three or more individuals (higher-order simplices), and it has been shown that models of diffusion on simplicial complexes can reproduce well the complex mechanisms of influence and reinforcement that are at work in the formation of opinions and in the adoption of novelties^67,68. At first, let us consider the case where all triangles are considered as 2-simplexes. In this way, the members of the Zachary karate club may have both pairwise, when they are connected by a network link, and three-body interactions, when they belong to the same triangle (2-simplex). The presence of a link indicates a social interaction among the two nodes of the link, whereas a 2-simplex can be interpreted as a social interaction involving three members of the club, such as a discussion to which all of them simultaneously participated. Oscillators have usually been used to describe the units of a coupled dynamical system when modeling opinion formation in social systems^50,69. As dynamical units we have decided to use chaotic oscillators, as it can be relevant to study synchronization in the more general scenario in which the opinions do not necessarily converge to a fixed stationary state^70,71. In particular, we associate to each node a Rössler oscillator and focus on the class III case, selecting the coupling functions as $$ {\mathbf{g}}^{\left(1\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)={\left[x_j^3-x_i^3,0,0\right]}^T$$ and $$ {\mathbf{g}}^{\left(2\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j,{\mathbf{x}}_k\right)={\left[x_j^2{x}_k-x_i^3,0,0\right]}^T$$. With these assumptions, the dynamics of each node i is described by

Equations (21) are then simulated for different values of σ₁ and σ₂. The average synchronization error and the predictions provided by the MSF (15) are illustrated in Fig. 6a that shows the crucial role played by the pairwise links, as synchronization turns out to be impossible when only three-body interactions are considered, i.e., when σ₁ = 0.

Fig. 6

Synchronization in Zachary karate club structure.

Synchronization is studied in simplicial complexes extracted from the interactions characterizing the Zachary karate club network. a Synchronization error (color code reported in the bar at the right of the panel) vs. σ₁ and σ₂ for the simplicial complex obtained when all the triangles are considered as being 2-simplexes. The red line delimits the area of stability of the synchronous solution predicted by the MSF. b λ₂ vs. the percentage of 2-simplexes in the structure, p_2s (see text for definition); c λ₂/λ_N vs. p_2s. In b and c three different values of r are considered, with the color code for the plotted curves being reported in the corresponding insets.

Next, we take the original network, and build different simplicial complexes by considering an increasing percentage (labeled as p_2s) of triangles in the original structure as true 2-simplexes, that is, an increasing percentage of social interactions taking place among groups of three members of the club. For each of these structures, we determine the effective matrix $$ \mathcal{M}$$ in (41), and calculate its spectrum of eigenvalues, and in particular we calculate the quantities $$ {\lambda }_2\left(\mathcal{M}\right)$$ and $$ {\lambda }_2\left(\mathcal{M}\right)/{\lambda }_N\left(\mathcal{M}\right)$$ (to simplify the notation here we shortly refer to these quantities as λ₂ and λ₂/λ_N). The former quantity provides the scaling of synchronization for class II systems, while the latter quantity (λ₂/λ_N) is a proxy of synchronizability for class III systems. The larger are the two quantities, the easier is to obtain synchronization. Fig. 6b, c illustrates the results at three values of $$ r=\frac{{\sigma }_2}{{\sigma }_1}$$. One finds that increasing p_2s has the effect of increasing λ₂ (thus it facilitates synchronization in class II systems), but simultaneously dwindles λ₂/λ_N (thus hindering synchronization in class III). Furthermore, Fig. 6 reveals that a larger value of $$ r=\frac{{\sigma }_2}{{\sigma }_1}$$ leads to larger values of λ₂, but smaller values of λ₂/λ_N, thus suggesting a beneficial impact of stronger three-body interactions for class II systems and an opposite effect on class III systems.

Cluster synchronization in simplicial complexes

In complex networks, symmetries may induce cluster synchronization, a regime where nodes group into clusters of units synchronized to each other³⁰. Network symmetries are permutations of the nodes preserving the connectivity pattern; they form a mathematical group, where each element may be represented by a permutation matrix R with elements r_ij = 1 if nodes i and j permute, and r_ij = 0 otherwise. The relevant property is that the orbits of the symmetry group associated to the network represent a partition into clusters that contain nodes that may synchronize each other. Group-theoretical considerations determine the exact composition and stability of clusters defined by the symmetry group³⁰.

Extending the notion of cluster synchronization to simplicial complexes requires a formal definition of symmetries in the general framework of multi-body interactions. Although this goes beyond the purpose of this paper, here, we illustrate an example of a simplicial complex displaying cluster synchronization. The core idea is that symmetry-related nodes must be flow invariant. Since the flow now comprises pairwise as well as higher-order interactions, symmetries must preserve the invariance for all the interactions taking place in the simplicial complex. Note that the same general principle is at the basis of the onset of cluster synchronization in multilayer networks where the symmetries guarantee that synchronized nodes have equal dynamical variables when inter-layer coupling is also included⁷².

In simplicial complexes, the flow invariance of symmetry-related nodes is obtained when the same symmetries hold for all the Laplacians involved in the dynamical equations of the nodes. Indicating with R_i with i = 1, …, n_p the n_p permutation matrices describing representing the symmetry group associated to the cluster synchronization state, this requires that the Laplacians satisfy the following Lyapunov equations^30,73:

To illustrate our results, we consider the two simplicial complexes shown in Fig. 7a, b, that have the same set of links, but different 2-simplices. The symmetries existing for the common network backbone induce the following partition of the nodes: V₁ = {1, 2}, V₂ = {7, 8}, V₃ = {9, 10}, V₄ = {11, 12}, V₅ = {3}, V₆ = {4}, V₇ = {5}, V₈ = {6}. Hence, the two simplicial complexes satisfy Eq. (22) for this symmetry group. However, for the simplicial complex in Fig. 7a, Eq. (22) do not hold, while the simplicial complex in Fig. 7b satisfies them.

Fig. 7

Cluster synchronization in simplicial complexes of Rössler oscillators.

a A simplicial complex where the symmetries of $$ {\mathcal{L}}^{\left(2\right)}$$ do not match those of $$ {\mathcal{L}}^{\left(1\right)}$$. b A simplicial complex where the symmetries of $$ {\mathcal{L}}^{\left(2\right)}$$ match those of $$ {\mathcal{L}}^{\left(1\right)}$$. c Synchronization error as a function of σ₁ for the simplicial complex ia. d Synchronization error as a function of σ₁ for the simplicial complex in panel b. In both cases, σ₂ has been set to σ₂ = 0.2.

We numerically study cluster synchronization monitoring the average synchronization error in each non-trivial cluster $$ e_{V_h}={⟨\sqrt{{\sum }_{i,j\in V_h}\parallel {\mathbf{x}}_i\left(t\right)-{\mathbf{x}}_j\left(t\right){\parallel }^2}⟩}_T$$ for h = 1, …, 4 and the average overall synchronization error E as in Eq. (45) to measure the onset of complete synchronization. Figure 7c, d shows the results for σ₁ ∈ [0, 1] and σ₂ = 0.2. As it can be observed, for the simplicial complex of Fig. 7b there is a interval of values of σ₁, i.e., σ₁ ∈ [0.1, 0.62], where the units in the four clusters are synchronized in the absence of complete synchronization, whereas for the simplicial complex of Fig. 7a V₁ is the only cluster that synchronizes before complete synchronization is achieved. These results show that higher-order interactions can modulate the pattern of synchronization emerging in the simplicial complex, as a diverse arrangement of the same number of 2-simplices (two in our example) led to different synchronous clusters.

Networks and higher-order structures of interactions

A network is a collection of nodes and of edges connecting pairs of nodes. Mathematically, it is represented by a graph $$ \mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$$, which consists of a set $$ \mathcal{V}$$ with $$ N=\mid \mathcal{V}\mid $$ elements called vertices (or nodes), and a set $$ \mathcal{E}$$ whose K elements, called edges or links, are pairs of nodes (i, j) (i, j = 1, 2, …, N and i ≠ j). As graphs explicitly refer to pairwise interactions, networks have been very successful in capturing the properties of coupled dynamical systems in all such cases in which the interactions can be expressed (or approximated) as a sum of two-body terms⁷⁵. Conversely, their limits emerge when it comes to model higher-order interactions. In fact, the presence of a triangle of three nodes i, j, k in a network, e.g., the presence of the three links (i, j), (i, k), (j, k) in the corresponding graph, is not able to capture the difference between a three-body interaction of the three individuals, from the sum of three pairwise interactions. Notice that these are two completely different situations, with completely different social mechanisms and dynamics at work⁶⁷.

Simplicial complexes are instead the proper mathematical structures for describing high-order interactions. A simplicial complex is an aggregate of simplices, objects that generalize links and can in general be of different dimension. A d-simplex, or simplex of dimension d, σ is, in its simplest definition, a collection of d + 1 nodes. In this way, a 0-simplex is a node, a 1-simplex is a link, a 2-simplex (i, j, k) is a two-dimensional object made by three nodes, usually called a (full) triangle, a 3-simplex is a tetrahedron, i.e., a three-dimensional object and so on. It is now possible to differentiate between a three-body interaction, and three bodies in pairwise interactions: the first case will be represented by a complete triangle, a two-dimensional simplex, while the second case will consist of three one-dimensional objects. Hence, in the following of this paper, simplices of dimension d will be used to describe the structure of (d + 1)-body interactions.

Finally, a simplicial complex $$ \mathcal{S}$$ on a given set of nodes $$ \mathcal{V}$$, with $$ \mid \mathcal{V}\mid =N$$, is a collection of M simplices, $$ \mathcal{S}=\left\{{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_M\right\}$$, with the extra requirement that, for any simplex $$ \sigma \in \mathcal{S}$$, all the simplices $$ {\sigma }^{\prime }$$ with $$ {\sigma }^{\prime }\subset \sigma $$, i.e., all the simplices built from subsets of σ, are also contained in $$ \mathcal{S}$$. Due to this requirement, simplicial complexes are a very particular type of hypergraphs⁷⁶. Simplicial complexes have shown to be appropriate in the context of social systems^67,77,78 and they will turn very useful to study coupled dynamical systems. We indicate as M_d, d = 1, 2, …, D the number of d-simplices present in $$ \mathcal{S}$$ (where D, the order of the simplicial complex, is the dimension of the largest simplex in $$ \mathcal{S}$$), with the constraint that $$ {\sum }_{d=1}^D{M}_d=M$$.

As a mathematical representation of simplicial complexes, we use here a formalism which generalizes directly the concept of adjacency matrix for a network. The adjacency matrix A of a graph $$ \mathcal{G}$$ is a N × N matrix, such that entry a_ij is 1 when edge $$ \left(i,j\right)\in \mathcal{E}$$, and 0 otherwise. The idea can be extended to simplicial complexes by considering tensors instead of matrices. In fact, for each dimension d, we can define the $$ {\underbrace{N\times N\times \dots \times N}}_{d+1}$$ adjacency tensor A^(d), whose entry $$ a_{i_1,\dots ,i_{d+1}}^{\left(d\right)}$$ is equal to 1 if the d-simplex (i₁, …, i_d+1) belongs to the simplex $$ \mathcal{S}$$, and is 0 otherwise¹⁷. Notice that each tensor is symmetric with respect to its d + 1 indices, which means that the value of a given entry $$ a_{i_1,\dots ,i_{d+1}}^{\left(d\right)}$$ is equal to the value of the entries corresponding to any permutation of the indices.

With the definition above, A⁽¹⁾ coincides with the standard adjacency matrix A, while the N  ×  N  ×  N adjacency tensor A⁽²⁾ characterizes two-dimensional objects: one has $$ a_{ijk}^{\left(2\right)}=1$$ if the three nodes i, j, k form a full triangle, and otherwise $$ a_{ijk}^{\left(2\right)}=0$$. As a conclusion, it is possible to map completely the connectivity structure of a simplicial complex $$ \mathcal{S}$$ into the entire set of D adjacency tensors A^(d), d = 1, 2, …, D.

A node i of a simplicial complex $$ \mathcal{S}$$ cannot be, therefore, characterized only by giving its degree $$ k_i={\sum }_{j=1}^N{a}_{ij}^{\left(1\right)}$$, but one needs instead to account for the number of simplices of any dimension, incident in i. It is therefore extremely useful to define the generalized d-degree, $$ k_i^{\left(d\right)}$$, of a node i as

Analogously, we can also define the generalized d-degree $$ k_{ij}^{\left(d\right)}$$ of a link (i, j) as the number of d-simplices to which link (i, j) is part of. We can write its expression in terms of the adjacency tensor A^(d) of dimension d, with d = 1, 2, …, D, as¹⁷

The Laplacian is a matrix that is of particular importance in many linear processes such as diffusion in graphs, but also turns useful in the linearization of nonlinear systems, for instance when we study the stability of a synchronized state in a networked dynamical system. The Laplacian matrix L = {l_ij} of a graph can be defined as L = K − A, where K is the diagonal matrix having the node degrees as diagonal elements. We give here a definition of generalized Laplacian describing the case of systems with high-order interactions. The generalized Laplacian of order d, with d = 1, 2, …, D, is a matrix $$ {\mathcal{L}}^{\left(d\right)}$$ whose elements are defined as

Notice that $$ {\mathcal{L}}^{\left(1\right)}$$ recovers exactly the classical Laplacian matrix. This definition of generalized Laplacian will turn useful in our study.

Derivation of Eq. (3)

To derive the conditions for the stability of the synchronization solution x^s, one first considers small perturbations around the synchronous state, i.e., δx_i = x_i − x^s, and performs a linear stability analysis of Eq. (2). One has

Then, one can factor out the terms $$ \frac{\partial {\mathbf{g}}^{\left(1\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)}{\partial {\mathbf{x}}_i}{\mid }_{\left({\mathbf{x}}^s,{\mathbf{x}}^s\right)}\delta {\mathbf{x}}_i$$ and $$ \frac{\partial {\mathbf{g}}^{\left(2\right)}\left({\mathbf{x}}_i,{\mathbf{x}}_j,{\mathbf{x}}_k\right)}{\partial {\mathbf{x}}_i}{\mid }_{\left({\mathbf{x}}^s,{\mathbf{x}}^s,{\mathbf{x}}^s\right)}\delta {\mathbf{x}}_i$$ in the summations (both of them, indeed, do not depend on the indices of the summations). Furthermore, one has that $$ {\sum }_{j=1}^N{a}_{ij}^{\left(1\right)}=k_i^{\left(1\right)}$$ and $$ \sum _{j=1}^N{\sum }_{k=1}^N{a}_{ijk}^{\left(2\right)}=2k_i^{\left(2\right)}$$. Plugging back the resulting terms inside the summations, and using Eq. (31), one eventually obtains

Already at this stage, it is fundamental to remark that our approach does not require a diffusive functional form for the interplay among the network nodes, and therefore we are actually encompassing an extremely broad class of coupling functions. For instance, our approach allows the formal treatment of the Kuramoto model⁵⁰, where m = 1, each network unit i is identified by the instantaneous phase θ_i of an oscillator, and the coupling between nodes i and j is given by the function $$ \sin \left({\theta }_j-{\theta }_i\right)$$, which is not diffusive.

Let us now make our second, conceptual, step, which will allow us to greatly simplify the last term on the right hand side of Eq. (32). Such a term refers to three-body interactions, and we now show how to map it into a single summation involving the generalized Laplacian matrix. This is done by remarking that the two Jacobian matrices J₁g⁽²⁾(x^s, x^s, x^s) and J₂g⁽²⁾(x^s, x^s, x^s) are both independent on k and j. Accordingly, Eq. (32) becomes

Then, using the symmetric property of T, namely $$ {\sum }_k{\tau }_{ijk}={\sum }_k{\tau }_{ikj}$$, we have

Let us now rewrite Eq. (35) in block form by introducing the stack vector $$ \delta \mathbf{x}={\left[\delta {\mathbf{x}}_1^T,\delta {\mathbf{x}}_2^T,\dots ,\delta {\mathbf{x}}_N^T\right]}^T$$ and denoting by JF = Jf(x^s), JG⁽¹⁾ = Jg⁽¹⁾(x^s, x^s) and JG⁽²⁾ = J₁g⁽²⁾(x^s, x^s, x^s) + J₂g⁽²⁾(x^s, x^s, x^s). One obtains

The third, and final, conceptual step is to remark that all generalized Laplacians $$ {\mathcal{L}}^{\left(d\right)}$$ are symmetric real-valued zero-row-sum matrices. Therefore: (i) they are all diagonalizable; (ii) for each one of them the set of eigenvalues is made of real non-negative numbers, and the corresponding set of eigenvectors constitutes a orthonormal basis of $$ {\mathbb{R}}^N$$; (iii) they all share, as the smallest of their eigenvalues, λ₁ ≡ 0, whose associated eigenvector $$ \frac{1}{\sqrt{N}}{\left(1,1,1,...,1\right)}^T$$ is aligned along the synchronization manifold; (iv) as in general they do not commute, the sets of eigenvectors corresponding to all others of their eigenvalues are different from one another, and yet any perturbation to the synchronization manifold (which, by definition, lies in the tangent space) can be expanded as linear combination of one whatever of such eigenvector sets (the relevant consequence is that one can arbitrarily select any of the generalized Laplacians as the reference for the choice of the basis of the transverse space, and all other eigenvector sets will map to such a basis by means of unitary matrix transformations).

We are then fully entitled to take, as reference basis, the one constituted by the eigenvectors of the classic Laplacian $$ {\mathcal{L}}^{\left(1\right)}$$ (V = [v₁, v₂, …, v_N]), and consider new variables η = (V⁻¹ ⊗ I_m)δx. We get

As $$ {\mathcal{L}}^{\left(2\right)}$$ is zero-row sum (i.e. $$ {\mathcal{L}}^{\left(2\right)}{\mathbf{v}}_1=0$$), Eq. (3) are finally obtained from Eq. (38).

Derivation of the MSF in the case of natural coupling

In the case of natural coupling as in Eq. (14), one has that J₁h⁽²⁾(x^s, x^s) + J₂h⁽²⁾(x^s, x^s) = Jh⁽¹⁾(x^s). The consequence is that the equations of the linearized dynamics in Eq. (35) can be rewritten as follows

Alternatively, one can consider the zero-row sum, symmetric, effective matrix $$ \mathcal{M}$$, given by

Equation (41) allows to establish a formal full analogy between the case of a simplicial complex and that of a network with weights given by the coefficients of the effective matrix $$ \mathcal{M}$$. In this case, by diagonalizing the effective matrix $$ \mathcal{M}$$, the transverse modes can be fully decoupled such that Eq. (15) is obtained, which prompts for the definition of a single-parameter MSF.

Numerical simulations

In our numerical simulations, we used two paradigmatic chaotic systems for the study of synchronization in systems of coupled units. The isolated dynamics of the Rössler oscillator is described by

In both cases, the parameters are fixed so as the resulting dynamics is chaotic. Namely, for the Rössler oscillator we selected a = b = 0.2, c = 9, and for the Lorenz system σ = 10, ρ = 28, and β = 8/3.

Furthermore, as a real-world example, we have considered the HR model for the neuron, whose isolated dynamics is described by the set of equations

In all cases, the state of the system is monitored by the average synchronization error defined as

Numerical integrations of the simplicial complexes of chaotic units are performed by means of an Euler algorithm, with integration step δt = 10⁻⁴, in a windows of time equal to 2T with T = 500.

For the calculation of the maximum Lyapunov exponent of the transverse modes in Eq. (3) we used the algorithm reported in ref. ⁷⁹ (pp. 116–117) with the following parameters: integration step size δt = 10⁻³, number of iterations per cycle I = 10000, number of cycles C = 5.

For the calculation of the MSF (15) we made use of the algorithm for the computation of the entire spectrum of Lyapunov exponents in ref. ⁸⁰ (with parameters: integration step size of the Euler algorithm δt = 10⁻⁵, length of the simulation L = 2500, windows of averaging T = 0.9L).