Global synchronization on time-varying higher-order structures

Let us assume the functions $\vec{h}^{(d)}$ to satisfy the condition of natural coupling, namely

$$\begin{align} \vec{h}^{\left(D\right)}\left(\vec{x},\dots,\vec{x}\right) = \dots = \vec{h}^{\left(2\right)}\left(\vec{x},\vec{x}\right) = \vec{h}^{\left(1\right)}\left(\vec{x}\right), \end{align} \tag{ 10 }$$

that implies $\mathbf{J}_{h^{(1)}} = \mathbf{J}_{h^{(2)}} = \dots = \mathbf{J}_{h^{(D)}}$ and it allows to eventually rewrite equation (9) as follows

$$\begin{align} \delta\dot{\vec{x}}_i = \mathbf{J}_{f}\delta\vec{x}_i+\sum_{j = 1}^n M_{ij}\left(t\right)\mathbf{J}_{h^{\left(1\right)}}\delta\vec{x}_{j}, \end{align} \tag{ 11 }$$

where

$$\begin{align} M_{ij}\left(t\right) : = \sum_{d = 1}^D q_d L^{\left(d\right)}_{ij}\left(t\right)\quad \forall i,j = 1,\dots n .\end{align} \tag{ 12 }$$

Let us observe that the matrix $\mathbf{M}(t)$ is a Laplace matrix; it is non-positive definite (as each one of the $\mathbf{L}^{(d)}(t)$ matrices does for any $d = 1,\dots, D$ and any t > 0, and $q_d\gt0$), it admits $\mu^{(1)} = 0$ as eigenvalue associated to the eigenvector $\phi^{(1)} = \frac{1}{\sqrt{N}}(1,\dots,1)^\top$ because each $\mathbf{L}^{(d)}(t)$ does the same and it is symmetric, being $\mathbf{L}^{(d)}(t)$ also symmetric for all $d = 1,\dots,D$. So there exists an orthonormal time-varying eigenbasis, $\phi^{(\alpha)}(t)$, $\alpha = 1,\dots,n$, for $\mathbf{M}(t)$ with associated eigenvalues $\mu^{(\alpha)} \lt 0$ for $\alpha = 2,\dots,n$. The goal is to use those eigenvectors to project (11) and thus to obtain the MSF; however, as we observed in the introduction, because the eigenvectors vary in time we have to take into account their time dependence. To achieve this goal we thus define [12] the n × n time dependent matrix $\mathbf{c}(t)$ that quantifies the projections of the time derivatives of the eigenvectors onto the independent eigendirections, namely

$$\begin{align} \frac{\mathrm{d} \vec{\phi}^{\left(\alpha\right)}\left(t\right)}{\mathrm{d}t} = \sum_{\beta}c_{\alpha\beta}\left(t\right)\vec{\phi}^{\left(\beta\right)}\left(t\right)\quad\forall \alpha = 1,\dots, n. \end{align} \tag{ 13 }$$

By recalling the orthonormality condition

$$\begin{align*} \left(\vec{\phi}^{\left(\alpha\right)}\left(t\right)\right)^\top\cdot \vec{\phi}^{\left(\beta\right)}\left(t\right) = \delta_{\alpha \beta} , \end{align*}$$

we can straightforwardly conclude that c is a real skew-symmetric matrix with a null first row and first column, i.e. $c_{\alpha\beta}+c_{\beta\alpha} = 0$ and $c_{1\alpha} = 0$.

Based on the above reasoning and to make one step further, we consider equation (11), and we project it onto the eigendirections, namely we introduce $\delta\vec{x}_i = \sum_\alpha \delta\hat{\vec{x}}_{\alpha}\phi^{(\alpha)}_i$ and recalling the definition of c we obtain

$$\begin{align} \frac{\mathrm{d}\delta\hat{\vec{x}}_{\beta}}{\mathrm{d}t} = \sum_\alpha c_{\beta\alpha}\left(t\right)\delta\hat{\vec{x}}_{\alpha}+\left[\mathbf{J}_{f}+ \mu^{\left(\beta\right)}\left(t\right)\mathbf{J}_{h^{\left(1\right)}}\right]\delta\hat{\vec{x}}_{\beta}. \end{align} \tag{ 14 }$$

Let us observe that the latter formula and the following analysis differ from the one presented in [40] where the perturbation is assumed to align onto a single mode, a hypothesis that ultimately translates in the stationary of the Laplace eigenvectors that is $\mathbf{c} = \boldsymbol{0}$. The same assumption is also at the root of the results presented in [41]; indeed, commuting time-varying networks implies to deal with a constant eigenbasis. In conclusion, equation (14) returns the more general description for the projection of the linearized dynamics on a generic time-varying Laplace eigenbasis, and thus allowing us to draw general conclusions without unnecessary simplifying assumptions. Let us observe that we did not introduce any assumption of the time-varying hypergraph and in particular we did not impose that nodes belonging to a k-hyperedge necessarily belong also to all k'-hyperedges, $1\unicode{x2A7D} k^{\prime} \lt k$ contained in the former, as it should be for a simplicial complex. Namely the previous analysis applies to both hypergraphs and simplicial complexes. Stated differently the adjacency tensors $\mathbf{A}^{(d)}(t)$, with $d = 1,\dots, D$ are independent each other, condition that will be no longer true in the case of regular topologies as we will show in the next section.

An alternative approach to study equation (9) is to assume regular topologies [25], namely hypergraphs such that $\mathbf{L}^{(d)}(t) = \alpha_d \mathbf{L}^{(1)}(t)$, for $d = 1,\dots,D$, with α₁ = 1 and $\alpha_d\in\mathbb{R}_+$. The latter have been introduced in [25] to propose simple structures allowing for the emergence of Turing patterns; observe that tetrahedra or icosahedra are examples of regular topologies once we assume unweighted, i.e. binary, and static adjacency tensors, the same holds true for the triangular lattice with periodic boundary conditions, i.e. a two-torus paved with triangles. Let us also emphasize that such topologies should not be confused with the concept of regular network or regular hypergraph, where, i.e. all nodes have the same degree. In the case of higher-order structures, the adjective regular refers to the fact that higher-order Laplace matrices differ only by a multiplicative factor and thus they are simultaneously diagonalizable.

By the very first property of regular topologies, equation (9) can be rewritten as

$$\begin{align} \delta\dot{\vec{x}}_i = \mathbf{J}_{f}\delta\vec{x}_i+\sum_{j = 1}^n L^{\left(1\right)}_{ij}\left(t\right)\mathbf{J}_{\hat{h}}\delta\vec{x}_{j}, \end{align} \tag{ 15 }$$

where

$$\begin{align} \mathbf{J}_{\hat{h}} : = \sum_{d = 1}^D q_d \alpha_d \mathbf{J}_{h^{\left(d\right)}}, \end{align} \tag{ 16 }$$

that results in a sort of weighted nonlinear coupling term. We can now make use of the existence of a time-varying orthonormal basis of $\mathbf{L}^{(1)}(t)$, namely $\psi^{(\alpha)}(t)$, $\alpha = 2,\dots,n$, associated to eigenvalues $\Lambda^{(\alpha)} \lt0$, $\psi^{(1)}(t) = \frac{1}{\sqrt{N}}(1,\dots,1)^\top$ and $\Lambda^{(1)} = 0$, to project $\delta\vec{x}_i$ onto the n eigendirections, $\delta\vec{x}_i = \sum_\alpha \delta\tilde{\vec{x}}_{\alpha}\psi^{(\alpha)}_i$. Because the latter vary in time we have again an extra term in the MSF and we thus need to define a second n × n time dependent matrix $\mathbf{b}(t)$ given by

$$\begin{align} \frac{\mathrm{d} \vec{\psi}^{\left(\alpha\right)}\left(t\right)}{\mathrm{d}t} = \sum_{\beta}b_{\alpha\beta}\left(t\right)\vec{\psi}^{\left(\beta\right)}\left(t\right)\quad\forall \alpha = 1,\dots, n, \end{align} \tag{ 17 }$$

that it is again real, skew-symmetric, with a null first row and first column, i.e. $b_{\alpha\beta}+b_{\beta\alpha} = 0$ and $b_{1\alpha} = 0$, because of the orthonormality condition of eigenvectors. By projecting equation (15) onto $\psi^{(\alpha)}(t)$, we get

$$\begin{align} \frac{\mathrm{d}\delta\tilde{\vec{x}}_{\beta}}{\mathrm{d}t} = \sum_\alpha b_{\beta\alpha}\left(t\right)\delta\tilde{\vec{x}}_{\alpha}+\left[\mathbf{J}_{f}+ \Lambda^{\left(\beta\right)}\left(t\right)\mathbf{J}_{\hat{h}}\right]\delta\tilde{\vec{x}}_{\beta} .\end{align} \tag{ 18 }$$

Let us conclude by observing that the latter equation has the same structure of (14). Those equations determine the generalization of the MSF to the case of time-varying higher-order structures. The time variation signature of the topology is captured by the matrices $\mathbf{c}(t)$ or $\mathbf{b}(t)$ and the eigenvectors $\mu^{(\alpha)}(t)$ or $\Lambda^{(\alpha)}(t)$, while the dynamics (resp. the coupling) in the Jacobian J_f (resp. $\mathbf{J}_{h^{(1)}}$ or $\mathbf{J}_{\hat{h}}$).

It is important to notice that as the eigenvalues $\mu^{(1)} = 0$, $\Lambda^{(1)} = 0$ and the skew-symmetric matrices $\mathbf{c}(t), \mathbf{b}(t)$ have null first row and column, in analogy with the MSF approaches carried over static networks [39] and higher-order structures [42], hence also in the case of time-varying higher-order structures, we can decouple the MSF into two components. One component describes the movement along the synchronous manifold, while the other component represents the evolution of different modes that are transverse to the synchronous manifold. The MSF ultimately relies on the computation of the maximum Lyapunov exponent (MLE) associated with the transverse modes and measures the exponential growth rate of a tiny perturbation in the transverse subspace. For the synchronous orbit to be stable, the MLE associated to all transverse modes must be negative. Moreover, the MSF approaches applied to static networks and higher-order structures can be simplified by examining the evolution of the perturbation along each independent eigendirection associated with distinct eigenvalues of the Laplacian matrix. Let us observe that this is not possible in the present because the matrices $\mathbf{c}(t)$ and $\mathbf{b}(t)$ mix the different modes and introduce a complex interdependence among them, making it challenging to disentangle their individual contributions. For this reason, one has to address numerically the problem [12].

To demonstrate the above introduced theory and emphasize the outcomes arising from the modified MSF (14) and (18), we will present two key examples in the following sections. Indeed, we will utilize the SL limit cycle oscillator and the chaotic Lorenz system as prototype dynamical systems anchored to each individual node. To simplify the calculations, we assume that the hypergraph consists of only four nodes, five links and two three-hyperedges, i.e. triangles, whose weights change in time (see appendix B for a detailed presentation of this structure). Let us stress that despite the simple hypergraph, the proposed framework is very general and can be applied to general time-varying hypergraphs. Observe on the other hand that the regular topology assumption induces constraints among the adjacency tensors, $\mathbf{A}^{(d)}(t)$, with $d = 1,\dots, D$, namely $\mathbf{A}^{(d)}(t) = \alpha_d \mathbf{A}^{(1)}(t)$. The latter ultimately implies dealing with a simplicial complex because of the presence of all k'-order connections, $1\unicode{x2A7D} k^{\prime}\lt k$, among nodes once they belong to a given k-hypergraph.

Let us thus assume the existence of two functions $h^{(1)}(w)$ and $h^{(2)}(w_1,w_2)$ such that $g^{(1)}$ and $g^{(2)}$ do satisfy the diffusive-like assumption, namely

$$\begin{align*} \begin{array}{l} g^{\left(1\right)}\left(w_j,w_{j_1}\right) = h^{\left(1\right)}\left(w_{j_1}\right)-h^{\left(1\right)}\left(w_j\right)\text{ and }\\\\g^{\left(2\right)}\left(w_j,w_{j_1},w_{j_2}\right) = h^{\left(2\right)}\left(w_{j_1},w_{j_2}\right)-h^{\left(2\right)}\left(w_j,w_j\right). \end{array} \end{align*}$$

For the sake of definitiveness, let us fix

$$\begin{align} h^{\left(1\right)}\left(w\right) = w \text{ and } h^{\left(2\right)}\left(w_1,w_2\right) = w_1w_2\, , \end{align} \tag{ 20 }$$

let us observe that the latter functions do not satisfy the condition for natural coupling, indeed $h^{(1)}(w) = w\neq w^2 = h^{(2)}(w,w)$.

Let us assume to deal with regular topology, namely $\mathbf{L}^{(2)} = \alpha_2\mathbf{L}^{(1)}$. Hence following equation (16) we can define $\mathbf{J}_{\hat{h}} = q_1 \mathbf{J}_{h^{(1)}}+q_2 \alpha_2 \mathbf{J}_{h^{(2)}}$. Let us perturb the limit cycle solution $w_\mathrm{LC}(t) = \sqrt{\sigma_\Re/\beta_\Re}\mathrm{e}^{\mathrm{i}\omega t}$ by defining $w_j = W_\mathrm{LC}(1+\rho_j)\mathrm{e}^{\mathrm{i}\theta_j}$, where ρ_j and θ_j are real and small functions for all j. A straightforward computation allows to write the time evolution of ρ_j and θ_j once we restrict ourselves to consider only linear terms

$$\begin{align} \dfrac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} {\rho_j} \\{\theta_j}\end{matrix}\right) & = \left(\begin{matrix} -2\sigma_\Re & 0\\-2\beta_\Im \frac{\sigma_\Re}{\beta_\Re} & 0 \end{matrix}\right)\left(\begin{matrix} {\rho_j} \\{\theta_j}\end{matrix}\right)+\sum_\ell L_{j\ell}^{(1)} \biggl[\left(\begin{matrix} q_{1,\Re} & - q_{1,\Im}\\ q_{1,\Im} & q_{1,\Re} \end{matrix}\right) \nonumber\\ & \quad + 2\alpha_2 \sqrt{\frac{\sigma_\Re}{\beta_\Re}} \left(\begin{matrix} \cos (\omega t) & - \sin (\omega t)\\ \sin (\omega t) & \cos (\omega t) \end{matrix}\right)\left(\begin{matrix} q_{2,\Re} & - q_{2,\Im}\\ q_{2,\Im} & q_{2,\Re} \end{matrix}\right)\biggr]\left(\begin{matrix} {\rho_\ell} \\{\theta_\ell}\end{matrix}\right) , \end{align} \tag{ 21 }$$

where $\omega = \sigma_\Im-\beta_\Im \sigma_\Re/\beta_\Re$ is the frequency of the reference limit cycle solution.

By exploiting the eigenvectors $\psi^{(\alpha)}(t)$ and eigenvalues $\Lambda^{(\alpha)}(t)$ of $\mathbf{L}^{(1)}(t)$ to project the perturbation ρ_j and θ_j we obtain:

$$\begin{align} \dfrac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} {\rho_\beta} \\ {\theta_\beta}\end{matrix}\right) & = \sum_\alpha b_{\beta\alpha}\left(\begin{matrix} {\rho_\alpha} \\{\theta_\alpha}\end{matrix}\right)+\bigg\{ \left(\begin{matrix} -2\sigma_\Re & 0\\-2\beta_\Im \frac{\sigma_\Re}{\beta_\Re} & 0 \end{matrix}\right) +\Lambda^{(\beta)} \bigg[\left(\begin{matrix} q_{1,\Re} & - q_{1,\Im}\\ q_{1,\Im} & q_{1,\Re} \end{matrix}\right) \nonumber \\ & \quad + 2\alpha_2 \sqrt{\frac{\sigma_\Re}{\beta_\Re}} \left(\begin{matrix} \cos (\omega t) & - \sin (\omega t)\\ \sin (\omega t) & \cos (\omega t) \end{matrix}\right)\left(\begin{matrix} q_{2,\Re} & - q_{2,\Im}\\ q_{2,\Im} & q_{2,\Re} \end{matrix}\right)\bigg]\bigg\}\left(\begin{matrix} {\rho_\beta} \\{\theta_\beta}\end{matrix}\right) , \end{align} \tag{ 22 }$$

where the matrix b has been defined in equation (17), $\rho_j = \sum_\beta \rho_\beta \psi^{(\beta)}_j(t)$ is the sough projection and similarly for θ_β .

As already mentioned, for the sake of definiteness and to focus on the impact of the time-varying topology, we hereby consider a simple higher-order network structure composed of n = 4 nodes, five links and two triangles with one link shared between both triangles. The number of nodes is thus fixed in time whereas the weight of the links and of the triangular face evolve in time. More precisely we consider the weights of the triangular faces to be given by $A_{\pi(123)}^{(2)}(t) = f(t)$ and $A_{\pi(234)}^{(2)}(t) = g(t)$, where f(t) and g(t) are generic positive smooth functions and $\pi(ijk)$ denotes any permutation of the indexes ijk. Whereas, the weights of the links, namely the pairwise adjacency matrix, are assumed to be given by

$$\begin{align} \begin{array}{l} A^{\left(1\right)}\left(t\right) = \begin{pmatrix} 0 & f\left(t\right) & f\left(t\right) & 0 \\ f\left(t\right) & 0 & f\left(t\right)+g\left(t\right) & g\left(t\right) \\ f\left(t\right) & f\left(t\right)+g\left(t\right) & 0 & g\left(t\right) \\ 0 & g\left(t\right) & g\left(t\right) & 0 \end{pmatrix}. \end{array} \end{align} \tag{ 23 }$$

This provides us a regular hypergraph satisfying the relation $\mathbf{L}^{(2)}(t) = \alpha_2\mathbf{L}^{(1)}(t)$ with $\alpha_{2} = 1$ (see appendix B for a detailed derivation). One can show that the Laplacian of the hypergraph possesses time-dependent eigenvalues and eigenvectors which are orthogonal each other and also to the constant eigenvector $\psi^{(1)}\sim (1,\dots,1)^\top$ corresponding to the zero eigenvalue. Consequently using the relation (17), we can obtain the skew symmetric projection matrix b such that all the entries of b are zero except $b_{34} = -b_{43}$ (we refer the interested reader to appendix B for a comprehensive derivation of the above matrix).

To realize how time-varying regular structure affects the global synchronization phenomenon, we fix, without loss of generality, the functions determining the evolution of edges and triangular faces such that $f(t) = 1+[\sin(\Omega t)]^2$ and $g(t) = 1+[\cos(\Omega t)]^2$. With this consideration, the case $\Omega = 0$ corresponds to the static higher-order network structure and $\Omega\gt0$ signifies the time-varying case. Figure B2 and the accompanying supplementary movie show the time evolution of the hypergraph structures for $\Omega = 2.0$.

Under all these assumptions, equation (22) determines a time periodic linear system whose stability can be determined by using Floquet theory. In order to illustrate our results, we let $q_{1,\Im}$ and $q_{2,\Im}$ to freely vary in the range $[-1.5,1.5]$ and $[-2.5,2.5]$, while keeping fixed to generic values the remaining parameters, and we compute the Floquet eigenvalue with the largest real part, corresponding thus to the MLE of equation (22), as a function of $q_{1,\Im}$ and $q_{2,\Im}$. The corresponding results are shown in figure 1 for $\Omega = 0$ (panel (a)) and $\Omega = 2$ (panel (b)). By a direct inspection, one can clearly conclude that the parameters region associated with a negative MLE (black region), i.e. to the stability of the SL limit cycle and thus to global synchronization, is larger for $\Omega \gt0$ than for $\Omega = 0$ in the regime of relatively smaller triadic coupling strength, i.e. $q_{2,\Im}$ small enough. However, for larger values of the latter, islands of positive MLE (yellow region) emerge for $\Omega = 2$, i.e. for the time-varying higher-order structure.

Zoom In Zoom Out Reset image size

To study the combined effect of both coupling strengths q₁ and q₂, we set $q_{1} = \epsilon_{1}q_{1,0}$ and $q_{2} = \epsilon_{2}q_{2,0}$, and we compute the MLE as a function of ε₁ and ε₂, having fixed without loss of generality $q_{1,0} = 0.1-0.5i$ and $q_{2,0} = 0.1+0.5i$. The corresponding results are presented in figure 2 for static ($\Omega = 0$, panel (a)) and time-varying ($\Omega = 2$, panel (b)) higher-order structure. We can again conclude that the region of parameters corresponding to global synchronization (black region) is larger in the case of time-varying hypergraph than in the static case when the strength of triadic interactions is small enough. Conversely, when the strength of triadic interactions increases, the region of parameters supporting global synchronization contracts in the time-varying hypergraph case as opposed to the static counterpart.

Zoom In Zoom Out Reset image size

Our last analysis concerns the relation between the frequency Ω and the size of the coupling parameters ε₁, ε₂, still assuming $q_{1} = \epsilon_{1}q_{1,0}$ and $q_{2} = \epsilon_{2}q_{2,0}$, on the onset of synchronization. In figure 3 we report the MLE in the plane $(\Omega,\epsilon_{1})$ for a fixed value of ε₂ (panel (a)), and in the plane $(\Omega,\epsilon_{2})$ for a fixed value of ε₁ (panel (b)). Let us observe that the synchronization can be easily achieved the smaller the value ε_j , $j = 1,2$, for which the MLE is negative, having fixed Ω. Let us thus define $\hat{\epsilon}_1(\Omega) = \min \{\epsilon \gt0 : \mathrm{MSF}(\epsilon,\epsilon_2,\Omega)\lt0\}$, for fixed ε₂, and similarly $\hat{\epsilon}_2(\Omega)$. The results of figure 3 clearly show that $\hat{\epsilon}_1(\Omega)\lt\hat{\epsilon}_1(0)\sim 1.2$ and $\hat{\epsilon}_2(\Omega)\lt\hat{\epsilon}_2(0)\sim 0.72$ and thus support our claim that time-varying structures allow to achieve synchronization easier. However, for larger values of Ω, one can observe that the global synchronization is forbidden when the strength of three-body coupling ε₂ is relatively larger (see figure 3(b)). Thus figure 3 illustrates that achieving synchronization is more attainable when the higher-order couplings are not overly intense. However, a noteworthy observation emerges: the interplay between a suitable level of higher-order coupling and temporality has the potential to induce desynchronization in higher-order structures.

Zoom In Zoom Out Reset image size

To complement our analysis, we performed numerical simulations of the SL defined on the simple four nodes time-varying hypergraph. We selected $(\epsilon_{1},\epsilon_{2}) = (1.3,0.5)$ and the remaining parameters values as in figure 2, we can thus conclude that for the chosen parameters, the MSF is positive if $\Omega = 0$ and negative if $\Omega = 2$, hence the SL should globally synchronize on the time-varying hypergraph while it would not achieve this state in the static case. Results of figure 4 confirm these conclusions; indeed, we can observe that (real part of) the complex variable is in phase for all i in the case $\Omega = 2$ (figure 4(b)), while this is not clearly the case for $\Omega = 0$ (figure 4(a)).

Zoom In Zoom Out Reset image size

Furthermore, to support the analytical insights suggesting that large values of Ω and ε₂ lead to desynchronization in the time-varying higher-order structure, we conducted numerical simulations of the SL model defined on the same four-node time-varying hypergraph. Specifically, we set $(\epsilon_{1},\epsilon_{2}) = (0.5,2.0)$ while keeping the other parameter values as those used in figure 2. From the latter figure one can conclude that the MSF is negative if $\Omega = 0$ and positive if $\Omega = 2$, hence the global synchronization is forbidden in the time-varying scenario. The findings presented in figure 5 validate this conclusion. Specifically, in the scenario where $\Omega = 0$ (figure 5(a)), the time series exhibit a synchronized state, as evidenced by the in-phase behavior of the real part of the complex variable across all nodes (i). On the contrary, this synchronized state is not observed in the case of $\Omega = 2$ (figure 5(b)), underscoring the influence of larger values of Ω and ε₂ in inducing desynchronization within the system.

Zoom In Zoom Out Reset image size

The results obtained so far rely on the use of a small time-varying hypergraph composed by four nodes, we claim however that their validity goes beyond this simple framework. To support this statement, we studied generic random large hypergraphs obtained by using a slightly different approach than the presented above. Indeed, we assume the Laplace eigenvalues to be constant and the time-derivative of the associated eigenvectors projected on the eigenbasis to return a constant matrix b. Namely we generated 100 random skew-symmetric matrices b (with zeros in their first row and column), we assigned random sets of eigenvalues $\Lambda^{(\beta)}\unicode{x2A7D} 0$ and we numerically solve equation (22) to compute the MSF for the corresponding structures. Then we determined the smallest nontrivial zero of the MSF, denoted as $\epsilon_{1}^{(0)}(\mathbf{b})$, while keeping the other coupling strength ε₂ fixed at a nominal value. For comparison, we computed the same quantity for a null matrix (i.e. $\mathbf{b} = \boldsymbol{0}$), representing a static hypergraph, denoted as $\epsilon_{1}^{(0)}(\boldsymbol{0})$, while keeping the same random eigenvalues.

Our claim is that synchronization can be achieved earlier in time-varying hypergraph, namely synchronization emerges for a smaller value of the parameter, ε₁ (once ε₂ has been held fixed), stated differently one should have $\Delta(\mathbf{b}): = \epsilon_{1}^{(0)}(\boldsymbol{0})-\epsilon_{1}^{(0)}(\mathbf{b})\gt0$. This is indeed what we show in the left panel of figure 6 for the case of time-varying hypergraphs made by 100 nodes and on the right panel for 200 nodes. In those figures we report the probability distribution of the values $\Delta(\mathbf{b})$ for 100 hypergraphs build as described above, and we can clearly appreciate that the probability to have $\Delta(\mathbf{b})\unicode{x2A7D} 0$ is null. Our results support thus the claim that the SL model defined on top of a temporal hypergraph consistently encompasses a broader range of coupling strengths associated to synchronization, with respect to the scenario of a static hypergraph. Furthermore, it is noteworthy that the minimum value of $\Delta(\mathbf{b})$ for hypergraphs composed of 200 nodes is comparably larger than the one for hypergraphs composed of 100 nodes, which underscores that SL oscillators coupled with large stationary hypergraphs determine a more formidable challenge for synchronization when compared to the case of time-varying hypergraphs.

Zoom In Zoom Out Reset image size

The aim of this section is to replace the condition of regular topology with a condition of natural coupling and consider again, a diffusive-like coupling. Let us thus consider now two functions $h^{(1)}(w)$ and $h^{(2)}(w_1,w_2)$ satisfying the natural coupling assumption, namely

$$\begin{align*} h^{\left(1\right)}\left(w\right) = h^{\left(2\right)}\left(w,w\right)\, . \end{align*}$$

For the sake of definitiveness, let us fix

$$\begin{align} h^{\left(1\right)}\left(w\right) = w^3 \text{ and } h^{\left(2\right)}\left(w_1,w_2\right) = w_1\left(w_2\right)^2 .\end{align} \tag{ 24 }$$

Consider again to perturb the limit cycle solution $w_\mathrm{LC}(t) = \sqrt{\sigma_\Re/\beta_\Re}\mathrm{e}^{\mathrm{i}\omega t}$ by defining $w_j = W_\mathrm{LC}(1+\rho_j)\mathrm{e}^{\mathrm{i}\theta_j}$, where ρ_j and θ_j are real and small functions for all j. A straightforward computation allows us to write the time evolution of ρ_j and θ_j as,

$$\begin{align} \dfrac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} {\rho_j} \\{\theta_j}\end{matrix}\right) = \left(\begin{matrix} -2\sigma_\Re & 0\\-2\beta_\Im \frac{\sigma_\Re}{\beta_\Re} & 0 \end{matrix}\right)\left(\begin{matrix} {\rho_j} \\{\theta_j}\end{matrix}\right) +3\frac{\sigma_\Re}{\beta_\Re}\sum_\ell M_{j\ell} \left(\begin{matrix} \cos \left(2\omega t\right) & - \sin \left(2\omega t\right)\\ \sin \left(2\omega t\right) & \cos \left(2\omega t\right) \end{matrix}\right)\left(\begin{matrix} {\rho_l} \\{\theta_l}\end{matrix}\right) , \end{align} \tag{ 25 }$$

where $\omega = \sigma_\Im-\beta_\Im \sigma_\Re/\beta_\Re$ is the frequency of the limit cycle solution and M is the matrix $q_1 \mathbf{L}^{(1)}(t)+q_2 \mathbf{L}^{(2)}(t)$ (see equation (12)). Let us observe that in this case, the coupling parameters q₁ and q₂ should be real numbers if we want to deal with real Laplace matrices, hypothesis that we hereby assume to hold true.

Proceeding along the same ideas developed in the previous section, we invoke the eigenvectors, $\phi^{(\alpha)}(t)$, eigenvalues, $\mu^{(\alpha)}(t)$, of $\mathbf{M}(t)$, and the matrix c (see equation (13)), to project the perturbation ρ_j and θ_j on the eigenbasis and thus rewrite the time variation of the perturbation as follows

$$\begin{align} \dfrac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} {\rho_\beta} \\{\theta_\beta}\end{matrix}\right) = \sum_\alpha c_{\beta\alpha}\left(\begin{matrix} {\rho_\alpha} \\{\theta_\alpha}\end{matrix}\right)+\left[ \left(\begin{matrix} -2\sigma_\Re & 0\\-2\beta_\Im \frac{\sigma_\Re}{\beta_\Re} & 0 \end{matrix}\right) +3\frac{\sigma_\Re}{\beta_\Re} \mu^{\left(\beta\right)} \left(\begin{matrix} \cos \left(2\omega t\right) & - \sin \left(2\omega t\right)\\ \sin \left(2\omega t\right) & \cos \left(2\omega t\right) \end{matrix}\right)\right]\left(\begin{matrix} {\rho_\beta} \\{\theta_\beta}\end{matrix}\right). \end{align} \tag{ 26 }$$

Let us assume to deal with a hypergraph made by three nodes, three links and one three-hyperedge, i.e. a triangular face, where links weights and face weight evolve in time. Consider then a time-independent matrix c

$$\begin{align*} \mathbf{c} = \begin{pmatrix} 0 &0 &0 \\ 0 & 0 & \Omega \\ 0 & -\Omega & 0 \end{pmatrix} , \end{align*}$$

for some $\Omega \unicode{x2A7E} 0$. The eigenvalue $\mu^{(1)} = 0$ of M determines the dynamics parallel to the synchronous manifold. On the other hand, the equations obtained for $\mu^{(2)}$ and $\mu^{(3)}$ give the dynamics of transverse modes to the synchronization manifold. Hence the MSF can be obtained by solving equation (26) and provides the conditions for a global stable synchronous solution to exist. In figure 7, we show the level sets of the MSF as a function of the eigenvalues $\mu^{(2)}$ and $\mu^{(3)}$ while keeping the remaining parameters in equation (26) fixed at generic nominal values. In panel (a), we consider a static hypergraph, i.e. $\Omega = 0$, while in panel (b) a time-varying hypergraph, i.e. $\Omega = 2$, negative values of MSF are reported in black and they correspond thus to a global synchronous state, positive values of MSF are shown in yellow. One can clearly appreciate that in the case of time-varying hypergraph, the MSF is negative for a much larger set of eigenvalues $\mu^{(2)}$ and $\mu^{(3)}$ and thus the SL system can easily synchronize.

Zoom In Zoom Out Reset image size

To validate the above results we again consider the SL oscillator with a particular case of non-invasive coupling function, namely we assume to exist a real function ϕ such that $\varphi(0) = 0$, $\varphi^{\prime}(0)\neq 0$ and

$$\begin{align} \begin{array}{l} g^{\left(1\right)}\left(w_1,w_2\right) = \varphi\left(w_1-w_2\right), \text{ and } \\ g^{\left(2\right)}\left(w_1,w_2,w_3\right) = \varphi\left(w_1-w_2\right)+\varphi\left(w_1-w_3\right). \end{array} \end{align} \tag{ A.10 }$$

By reasoning as before, we get

$$\begin{align} \begin{array}{l} \dfrac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} {\rho_j} \\{\theta_j}\end{matrix}\right) = \left(\begin{matrix} -2\sigma_\Re & 0\\-2\beta_\Im \frac{\sigma_\Re}{\beta_\Re} & 0 \end{matrix}\right)\left(\begin{matrix} {\rho_j} \\{\theta_j}\end{matrix}\right)+ \varphi^{\prime}\left(0\right)\sum_\ell \left(q_1 L^{\left(1\right)}_{j\ell} +q_2 L^{\left(2\right)}_{j\ell}\right) \left(\begin{matrix} 1 & 0\\ 0 & -1 \end{matrix}\right)\left(\begin{matrix} {\rho_l} \\{\theta_l}\end{matrix}\right). \end{array} \end{align} \tag{ A.11 }$$

By using again the eigenvectors $\phi^{(\alpha)}(t)$, eigenvalues $\mu^{(\alpha)}(t)$ of $\mathbf{G}(t)$ and the matrix c (see equation (A.8)), we can rewrite the previous formula as

$$\begin{align} \begin{array}{l} \dfrac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} {\rho_\beta} \\{\theta_\beta}\end{matrix}\right) = \sum_\alpha c_{\beta\alpha}\left(\begin{matrix} {\rho_\alpha} \\{\theta_\alpha}\end{matrix}\right)+\left[\left(\begin{matrix} -2\sigma_\Re & 0\\-2\beta_\Im \frac{\sigma_\Re}{\beta_\Re} & 0 \end{matrix}\right) + \varphi^{\prime}\left(0\right)\mu^{\left(\beta\right)} \left(\begin{matrix} 1 & 0\\ 0 & -1 \end{matrix}\right)\right]\left(\begin{matrix} {\rho_\beta} \\{\theta_\beta}\end{matrix}\right). \end{array} \end{align} \tag{ A.12 }$$

Figure A1 represent the result for the non-invasive coupling assumption. Here, we consider the non-invasive function so that $\varphi^{\prime}(0) = 1$ and the skew-symmetric projection matrix c is considered constant throughout the analysis as earlier. Here we show the level sets of the MSF as a function of the eigenvalues $\mu^{(2)}$ and $\mu^{(3)}$ while keeping the remaining parameters in equation (A.12) fixed at generic nominal values. In panel (a), we consider a static hypergraph, i.e. $\Omega = 0$, while in the (b) panel, a time-varying hypergraph, i.e. $\Omega = 2$, negative values of MSF are reported in black, and they correspond thus to a global synchronous state, positive values of MSF are shown in yellow; one can clearly appreciate that in the case of the time-varying hypergraph, the MSF is negative for a much larger set of eigenvalues $\mu^{(2)}$ and $\mu^{(3)}$ and thus the SL system can achieve synchronization more easily.

Zoom In Zoom Out Reset image size

The aim of this section is to compute explicitly the spectrum of the toy example hypergraph presented in the previous section. By using a symbolic manipulator one can easily compute the eigenvalues and eigenvectors of the three-hypergraphs and thus obtain

$$\begin{align} \Lambda^{\left(1\right)} & = 0, \Lambda^{\left(2\right)} = -3\left(f\left(t\right)+g\left(t\right)\right), \nonumber\\ \Lambda^{\left(3,4\right)} & = -\dfrac{3}{2}\left(f\left(t\right)+g\left(t\right)\right) \pm \dfrac{1}{2} \sqrt{9\left[f\left(t\right)\right]^2-14f\left(t\right)g\left(t\right)+9\left[g\left(t\right)\right]^2}, \end{align} \tag{ B.10 }$$

and

$$\begin{align} \begin{array}{l} \psi^{\left(1\right)} = \dfrac{1}{2} \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix}, \psi^{\left(2\right)} = \dfrac{1}{\sqrt{2}} \begin{pmatrix} 0 \\ -1 \\ 1 \\ 0 \end{pmatrix}, \\ \psi^{\left(3,4\right)} = \dfrac{1}{4s_{3,4}g\left(t\right)} \begin{pmatrix} 6\left(f\left(t\right)-g\left(t\right)\right) \mp 2 \sqrt{9\left[f\left(t\right)\right]^2-14f\left(t\right)g\left(t\right)+9\left[g\left(t\right)\right]^2} \\ -\left(3f\left(t\right)-g\left(t\right)\right) \pm \sqrt{9\left[f\left(t\right)\right]^2-14f\left(t\right)g\left(t\right)+9\left[g\left(t\right)\right]^2}\\ -\left(3f\left(t\right)-g\left(t\right)\right) \pm \sqrt{9\left[f\left(t\right)\right]^2-14f\left(t\right)g\left(t\right)+9\left[g(t)\right]^2} \\ 1 \end{pmatrix}. \end{array} \end{align} \tag{ B.11 }$$

Let us observe that to increase the readability of the previous equation we did not explicitly wrote the normalizing factor s₃ (resp. s₄ ) for $\psi^{(3)}$ (resp. $\psi^{(4)}$). One can easily show that eigenvalues and eigenvectors are real for all possible functions f(t) and g(t), and to prove that they are orthonormal and thus they form a basis for $\mathbb{R}^{4}$ for all t.

A lengthy but straightforward calculation allows to obtain the matrix b (see equation (17) in the main text) that describes the evolution of the derivatives of the eigenvectors on the eigenbasis. Because of the form of the latter ones, we can show that all the entries of b are zero but $b_{34} = -b_{43}$ that is given by

$$\begin{align} \begin{array}{l} b_{34} = \dfrac{N_{34}}{D_{34}}, \end{array} \end{align} \tag{ B.12 }$$

where

$$\begin{align} N_{34} & = 864g^4f^{\prime} - 1344fg^3f^{\prime} - 864fg^3g^{\prime} - 864f^3gg^{\prime} + 864f^2g^2f^{\prime} + 1344f^2g^2g^{\prime} \nonumber\\ & \quad + \sqrt{9f^2 - 14fg + 9g^2} \left(567f^3g^{\prime} - 729f^3f^{\prime}+ 791g^3f^{\prime}- 729g^3g^{\prime}- 1899fg^2f^{\prime}+ 1701f^2gf^{\prime} \right. \nonumber\\ & \quad \left.+\ 1477fg^2g^{\prime} - 1323f^2gg^{\prime}\right) + \left( 9f^2 - 14fg + 9g^2\right)^{\frac{3}{2}} \left(81ff^{\prime}- 63fg^{\prime} - 63gf^{\prime} + 81gg^{\prime}\right), \end{align} \tag{ B.13 }$$

and

$$\begin{align} D_{34} & = 2 \left(27f^2 - 42f g + 27g^2 +\left(7g - 9f\right) \sqrt{9f^2 - 14fg + 9g^2}\right)^{\frac{3}{2}}\nonumber\\ & \quad \times \sqrt{27f^2 - 42fg + 27g^2+ \left(9f - 7g\right) \sqrt{9f^2 - 14fg + 9g^2}} \sqrt{9f^2 - 14fg + 9g^2}, \end{align} \tag{ B.14 }$$

where for ease of reading we removed the explicit dependence on t and ' denotes the first derivative with respect to t.

The numerical studies performed in the main text has been realized by assuming $f(t) = 1+[\sin(\Omega t)]^2$ and $g(t) = 1+[\cos(\Omega t)]^2$, in this case the eigenvalues read

$$\begin{align} \begin{array}{l} \Lambda^{\left(1\right)} = 0, \Lambda^{\left(2\right)} = -9, \Lambda^{\left(3,4\right)} = -\dfrac{9}{2} \pm \dfrac{\sqrt{4\cos\left(4\Omega t\right)+ 13}}{2}, \end{array} \end{align} \tag{ B.15 }$$

and the non-zero element of the matrix b is

$$\begin{align} \begin{array}{l} b_{34} = \dfrac{\tilde{N}_{34}}{\tilde{D}_{34}}, \end{array} \end{align} \tag{ B.16 }$$

where $\tilde{N}_{34} = 24 \Omega \sin (2 t \Omega ) [\cos (2 t \Omega )+3]$ and $\tilde{N}_{34} = 2 \sqrt{2} [\cos (2 t \Omega )+3] [4 \cos (4 t \Omega )+13]$.

Figure B2 portrays the temporal evolution of the links and three-hyperedges weights. To better understand the evolution of the hypergraph, we provide the graphical evolution of the hypergraph in the accompanying supplementary movie, together with the time evolution of the weights of the links $A^{(1)}_{ij}(t)$ and of the hyperedge $A^{(2)}_{ijk}(t)$.

Zoom In Zoom Out Reset image size