# Global synchronization on time-varying higher-order structures

Md Sayeed Anwar  Dibakar Ghosh  Physics and Applied Mathematics Unit, Indian
Statistical Institute, 203 B. T. Road, Kolkata 700108, India  Timoteo Carletti
Department of Mathematics and Namur Institute for Complex Systems, naXys,
University of Namur, 2 rue Grafé, Namur B5000, Belgium

###### Abstract

Synchronization has received a lot of attention from the scientific community
for systems evolving on static networks or higher-order structures, such as
hypergraphs and simplicial complexes. In many relevant real world
applications, the latter are not static but do evolve in time, in this paper
we thus discuss the impact of the time-varying nature of high-order structures
in the emergence of global synchronization. To achieve this goal we extend the
master stability formalism to account, in a general way, for the additional
contributions arising from the time evolution of the higher-order structure
supporting the dynamical systems. The theory is successfully challenged
against two illustrative examples, the Stuart-Landau nonlinear oscillator and
the Lorenz chaotic oscillator.

##  I Introduction

In the realm of complex systems, synchronization refers to the intriguing
ability of coupled nonlinear oscillators to self-organize and exhibit a
collective unison behavior without the need for a central controller Arenas
_et al._ (2008). This phenomenon, observed in a wide range of human-made and
natural systems Boccaletti  _et al._ (2018), continues to inspire scientists
seeking to unravel its underlying mechanisms.

To study synchronization, network science has proved to be a powerful and
effective framework. Here, the interconnected nonlinear oscillators are
represented as nodes, while their interactions are depicted as links Barabási
(2016). However, the classical static network representation has its
limitation in modeling many empirical systems, such as social networks
Wasserman  _et al._ (1994), brain networks Valencia  _et al._ (2008); Bassett
_et al._ (2011), where the connections among individual basic units are
adaptable enough to be considered to evolve through time. Therefore, the
framework of networks has been generalized as to include time-varying networks
Holme and Saramäki (2012); Masuda and Lambiotte (2016), whose connections vary
with time. The results presented in this framework support the claim that
synchronization is enhanced by the dynamics of the supporting medium Ghosh
_et al._ (2022); Carletti and Fanelli (2022); Anwar  _et al._ (2022).

Another intrinsic limitation of networks is due to their capability to only
model pairwise interactions. To go beyond this issue, scholars have brought to
the fore the relevance of higher-order structures, which surpass the
traditional network setting that models the interactions between individual
basic units only through pairwise links Carletti  _et al._ (2020a); Battiston
_et al._ (2020, 2021); Majhi  _et al._ (2022); Boccaletti  _et al._ (2023). By
considering the simultaneous interactions of many agents, higher-order
structures, namely hypergraphs Berge (1973) and simplicial complexes Bianconi
(2021), offer a more comprehensive understanding of complex systems. These
higher-order structures have been proven to produce novel features in various
dynamical processes, including consensus Neuhäuser  _et al._ (2020, 2021),
random walks Carletti  _et al._ (2020b); Schaub  _et al._ (2020), pattern
formation Carletti  _et al._ (2020a); Muolo  _et al._ (2023); Gao  _et al._
(2023), synchronization Carletti  _et al._ (2020a); Skardal and Arenas (2020,
2019); Carletti  _et al._ (2023); Anwar and Ghosh (2022a, b), social contagion
and epidemics Iacopini  _et al._ (2019); Chowdhary  _et al._ (2021).
Nevertheless, the suggested framework is not sufficiently general for
describing systems with many-body interactions that vary with time. As an
example, group interactions in social systems have time-varying nature as the
interactions among groups of individuals are not always active but rather
change throughout time Cencetti  _et al._ (2021). Some early works have begun
to investigate the time-varying aspect of many-body interactions in various
dynamical processes. For instance, time-varying group interactions have been
demonstrated to influence the convergence period of consensus dynamics
Neuhäuser  _et al._ (2021) and to predict the onset of endemic state in
epidemic spreading Chowdhary  _et al._ (2021).

The present work is motivated by these recent research directions, and it aims
to take one step further by considering the impact of time-varying higher-
order structures in the synchronization of nonlinear oscillators. In this
context, a preliminary effort has been reported in Anwar and Ghosh (2022c),
that investigates synchronization in time-varying simplicial complexes,
limited only to fast switching Stilwell  _et al._ (2006); Petit  _et al._
(2017) among distinct static simplicial configurations, implying that the time
scale of the simplicial evolution is exceedingly fast compared to that of the
underlying dynamical system. In contrast, in the present work, we allow the
higher-order structures to evolve freely with time, thus removing any
limitations on the imposed time evolution of the higher-order structure. We
present the results in the framework of hypergraphs, but they hold true also
for simplicial complexes. Under such broad circumstances, we develop a theory
to determine the conditions ensuring the stability of a globally synchronized
state that generalizes the Master Stability Equation Pecora and Carroll (1998)
to a setting where the time evolution of underlying higher-order structures is
explicitly considered. The generalized framework we discuss here assumes that
the coupling functions cancel out when the dynamics of individual oscillators
are identical, which is a necessary condition that must be met for the
extended system to have a synchronous solution and it has been frequently used
in the literature across various domains. The developed theory reveals that
the consideration of temporality in group interactions can induce
synchronization more easily than static group interactions, tested on higher-
order structures of coupled Stuart Landau oscillators and paradigmatic Lorenz
systems.

##  II The model

To start with, let us consider a $m$-dimensional dynamical system whose time
evolution is described by the following ordinary differential equation

$\frac{d\vec{x}}{dt}=\vec{f}(\vec{x})\,,$ (1)

where $\vec{x}\in\mathbb{R}^{m}$ denotes the state vector and
$\vec{f}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ some smooth nonlinear
function; let us assume moreover that system (1) exhibits an oscillatory
behavior, being the latter periodic or irregular; we are thus considering the
framework of generic nonlinear oscillators. Let us now consider $n$ identical
copies of system (1) coupled by a symmetric higher-order structure; namely, we
allow the nonlinear oscillators to interact in couples, as well as in
triplets, quadruplets, and so on, up to interactions among $D+1$ units. We can
thus describe the time evolution of the state vector of the $i$-th unit by

$\dot{\vec{x}}_{i}=\vec{f}(\vec{x_{i}})+\sum\limits_{d=1}^{D}q_{d}\sum\limits_{j_{1},\dots,j_{d}=1}^{n}A_{ij_{1}\dots
j_{d}}^{(d)}(t)\vec{g}^{(d)}(\vec{x}_{i},\vec{x}_{j_{1}},\dots,\vec{x}_{j_{d}})\,,$
(2)

where for $d=1,\dots,D$, $q_{d}>0$ denotes the coupling strength,
$\vec{g}^{(d)}:\mathbb{R}^{(d+1)m}\rightarrow\mathbb{R}^{m}$ the nonlinear
coupling function and $\mathbf{A}^{(d)}(t)$ the tensor encoding which units
are interacting together. More precisely ${A}^{(d)}_{ij_{1}\dots j_{d}}(t)=1$
if the units $i,j_{1},\dots,j_{d}$ do interact at time $t$, observe indeed
that such tensor depends on time, namely the intensity of the coupling as well
which units are coupled, do change in time. Finally, we assume the time-
varying interaction to be symmetric, namely if ${A}^{(d)}_{ij_{1}\dots
j_{d}}(t)=1$, then ${A}^{(d)}_{\pi(ij_{1}\dots j_{d})}(t)=1$ for any
permutation $\pi$ of the indexes $i,j_{1},\dots,j_{d}$. Let us emphasize that
we consider the number of nodes to be fixed, only the interactions change in
time; one could relax this assumption by considering to have a sufficiently
large reservoir of nodes, from which the core of the system can recruit new
nodes or deposit unused nodes.

Let us fix a periodic reference solution, $\vec{s}(t)$, of system (1). We are
interested in determining the conditions under which the orbit
$(\vec{s}(t),\dots,\vec{s}(t))^{\top}$ is a solution of the coupled system
(2), and moreover it is stable, namely the $n$ units globally synchronize and
behave at unison. A necessary condition is that the coupling functions vanish
once evaluated on such orbit, i.e., $\vec{g}^{(d)}(\vec{s},\dots,\vec{s})=0$,
for $d=1,\dots,D$. This assumption is known in the literature as non-invasive
condition.

For the sake of pedagogy, we will hereby consider a particular case of non-
invasive couplings and we will refer the interested reader to Appendix A for a
general discussion. We are thus assuming the coupling functions
$\vec{g}^{(d)}$ to be diffusive-like, namely for each $d$ there exists a
function $\vec{h}^{(d)}:\mathbb{R}^{dm}\rightarrow\mathbb{R}^{m}$ such that

$\vec{g}^{(d)}(\vec{x}_{i},\vec{x}_{j_{1}},\dots,\vec{x}_{j_{d}})=\vec{h}^{(d)}(\vec{x}_{j_{1}},\dots,\vec{x}_{j_{d}})-\vec{h}^{(d)}(\vec{x}_{i},\dots,\vec{x}_{i})\,.$
(3)

In this way we can straightforwardly ensure that the coupling term in Eq. (3)
vanishes once evaluated on the orbit $(\vec{s}(t),\dots,\vec{s}(t))^{\top}$,
allowing thus to conclude that the latter is also a solution of the coupled
system.

To study the stability of the reference solution, let us now perturb the
synchronous solution $(\vec{s}(t),\dots,\vec{s}(t))^{\top}$ with a spatially
inhomogeneous term, meaning that $\forall i\in\\{1,\dots,n\\}$ we define
$\vec{x}_{i}=\vec{s}+\delta\vec{x}_{i}$. Substituting the latter into Eq. (2)
and expanding up to the first order, we obtain

$\delta\dot{\vec{x}}_{i}=\frac{\partial\vec{f}}{\partial\vec{x}_{i}}\Big{\rvert}_{\vec{s}}\delta\vec{x}_{i}+\sum_{d=1}^{D}q_{d}\sum_{j_{1},\dots,j_{d}=1}^{n}B_{ij_{1}\dots
j_{d}}(t)\sum_{\ell=1}^{d}\frac{\partial\vec{h}^{(d)}}{\partial\vec{x}_{j_{\ell}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}\delta\vec{x}_{j_{\ell}}\,,$
(4)

where

$\displaystyle B_{ij_{1}}(t)$ $\displaystyle=$ $\displaystyle
A_{ij_{1}}^{(1)}(t)-k^{(1)}_{i}(t)\delta_{ij_{1}}\,,$ $\displaystyle
B_{ij_{1}j_{2}}(t)$ $\displaystyle=$ $\displaystyle
A_{ij_{1}j_{2}}^{(2)}(t)-2k_{i}^{(2)}(t)\delta_{ij_{1}j_{2}}\,,\dots$
$\displaystyle B_{ij_{1}j_{2}\dots j_{D}}(t)$ $\displaystyle=$ $\displaystyle
A_{ij_{1}j_{2}\dots j_{D}}^{(D)}(t)-D!k_{i}^{(D)}(t)\delta_{ij_{1}j_{2}\dots
j_{D}}\,,$

being $\delta_{ij_{1}j_{2}\dots j_{D}}$ the generalized multi-indexes
Kronecker-$\delta$, and the (time-varying) $d$-degree of node $i$ is given by

$k_{i}^{(d)}(t)=\frac{1}{d!}\sum_{j_{1},..,j_{d}=1}^{n}A_{ij_{1}\dots
j_{d}}^{(d)}(t)\,,$ (5)

which represents the number of hyperedges of order $d$ incident to node $i$ at
time $t$. Observe that if $\mathbf{A}^{(d)}$ is weighted, then
$k_{i}^{(d)}(t)$ counts both the number and the weight, it is thus the
generalization of the strength of a node. Let us now define

$k_{ij}^{(d)}(t)=\frac{1}{(d-1)!}\sum_{j_{1},...,j_{d-1}}^{n}A_{ijj_{1}\dots
j_{d-1}}^{(d)}(t)\,,$ (6)

namely the number of hyperedges of order $d$ containing both nodes $i$ and $j$
at time $t$. Again, once $\mathbf{A}^{(d)}$ is weighted, then
$k_{ij}^{(d)}(t)$ generalizes the link strength. Let us observe that because
of the invariance of $\mathbf{A}^{(d)}$ under index permutation, we can
conclude that $k_{ij}^{(d)}(t)=k_{ji}^{(d)}(t)$. Finally, we define the
generalized time-varying higher-order Laplacian matrix for the interaction of
order $d$ as

$L_{ij}^{(d)}(t)=\begin{cases}-d!k_{i}^{(d)}(t)&\text{if $i=j$}\\\
(d-1)!k_{ij}^{(d)}(t)&\text{if $i\neq j$}\end{cases}\,.$ (7)

Observe that such a matrix is symmetric because of the assumption of the
tensors $\mathbf{A}^{(d)}$. Let us also notice the difference in sign with
respect to other notations used in the literature.

We can then rewrite Eq. (4) as follows

$\displaystyle\delta\dot{\vec{x}}_{i}$ $\displaystyle=$
$\displaystyle\frac{\partial\vec{f}}{\partial\vec{x}_{i}}\Big{\rvert}_{\vec{s}}\delta\vec{x}_{i}+\sum_{d=1}^{D}q_{d}\left[\sum_{j_{1}=1}^{n}\frac{\partial\vec{h}^{(d)}}{\partial\vec{x}_{j_{1}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}\delta\vec{x}_{j_{1}}\sum_{j_{2},\dots,j_{d}=1}^{n}B_{ij_{1}\dots
j_{d}}(t)+\dots+\sum_{j_{d}=1}^{n}\frac{\partial\vec{h}^{(d)}}{\partial\vec{x}_{j_{d}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}\delta\vec{x}_{j_{d}}\sum_{j_{1},\dots,j_{d-1}=1}^{n}B_{ij_{1}\dots
j_{d}}(t)\right]$ (8) $\displaystyle=$
$\displaystyle\frac{\partial\vec{f}}{\partial\vec{x}_{i}}\Big{\rvert}_{\vec{s}}\delta\vec{x}_{i}+\sum_{d=1}^{D}q_{d}\sum_{j=1}^{n}L^{(d)}_{ij}(t)\left[\frac{\partial\vec{h}^{(d)}}{\partial\vec{x}_{j_{1}}}+\dots+\frac{\partial\vec{h}^{(d)}}{\partial\vec{x}_{j_{d}}}\right]_{(\vec{s},\dots,\vec{s})}\delta\vec{x}_{j}\,,$

where we used the fact the
$\frac{\partial\vec{h}^{(d)}}{\partial\vec{x}_{j_{1}}}+\dots+\frac{\partial\vec{h}^{(d)}}{\partial\vec{x}_{j_{d}}}$
is independent from the indexes being the latter just place holders to
identify the variable with respect to the derivative has to be done. Finally,
by defining

$\displaystyle\mathbf{J}_{f}$ $\displaystyle:=$
$\displaystyle\frac{\partial\vec{f}}{\partial\vec{x}_{i}}\Big{\rvert}_{\vec{s}(t)}\text{
and }$ $\displaystyle\mathbf{J}_{h^{(d)}}$ $\displaystyle:=$
$\displaystyle\sum_{\ell=1}^{d}\frac{\partial\vec{h}^{(d)}}{\partial\vec{x}_{j_{\ell}}}\Big{\rvert}_{(\vec{s}(t),\dots,\vec{s}(t))}\forall
d\in\\{1,\dots,D\\}\,,$

we can rewrite Eq. (8) in compact form

$\delta\dot{\vec{x}}_{i}=\mathbf{J}_{f}\delta\vec{x}_{i}+\sum_{d=1}^{D}q_{d}\sum_{j=1}^{n}L^{(d)}_{ij}(t)\mathbf{J}_{h^{(d)}}\delta\vec{x}_{j}\,.$
(9)

This is a non-autonomous linear differential equation determining the
stability of the perturbation $\delta\vec{x}_{i}$, for instance, by computing
the largest Lyapunov exponent. To make some analytical progress in the study
of Eq. (9), we will consider two main directions: the functions
$\vec{h}^{(d)}$ satisfy the condition of natural coupling (see Section II.1)
or the higher-order structures exhibit regular topologies (see Section II.2).
The aim of each assumption is to disentangle the dependence of the nonlinear
coupling functions from the higher-order Laplace matrices and thus achieve a
better understanding of the problem under study.

###  II.1 Natural coupling

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

$\vec{h}^{(D)}(\vec{x},\dots,\vec{x})=\dots=\vec{h}^{(2)}(\vec{x},\vec{x})=\vec{h}^{(1)}(\vec{x})\,,$
(10)

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

$\delta\dot{\vec{x}}_{i}=\mathbf{J}_{f}\delta\vec{x}_{i}+\sum_{j=1}^{n}M_{ij}(t)\mathbf{J}_{h^{(1)}}\delta\vec{x}_{j}\,,$
(11)

where

$M_{ij}(t):=\sum_{d=1}^{D}q_{d}L^{(d)}_{ij}(t)\quad\forall i,j=1,\dots n\,.$
(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}>0$), it admits $\mu^{(1)}=0$ as
eigenvalue associated to the eigenvector $\phi^{(1)}=(1,\dots,1)^{\top}$ and
it is symmetric. So there exists an orthonormal time-varying eigenbasis,
$\phi^{(\alpha)}(t)$, $\alpha=1,\dots,n$, for $\mathbf{M}(t)$ with associated
eigenvalues $\mu^{(\alpha)}\leq 0$. Let us define Carletti and Fanelli (2022)
the $n\times n$ time dependent matrix $\mathbf{c}(t)$ that quantifies the
projections of the time derivatives of the eigenvectors onto the independent
eigendirections, namely

$\frac{d\vec{\phi}^{(\alpha)}(t)}{dt}=\sum_{\beta}c_{\alpha\beta}(t)\vec{\phi}^{(\beta)}(t)\quad\forall\alpha=1,\dots,n\,.$
(13)

By recalling the orthonormality condition

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

we can straightforwardly conclude that $\mathbf{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$.

To make one step further, we consider Eq. (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 $\mathbf{c}$ we obtain

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

Let us observe that the latter formula and the following analysis differ from
the one presented in Van Gorder (2021) 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}=\mathbf{0}$. The
same assumption is also at the root of the results by Zhang and Strogatz
(2021); indeed, commuting time-varying networks implies to deal with a
constant eigenbasis. In conclusion, Eq. (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.

###  II.2 Regular topologies

An alternative approach to study Eq. (9) is to assume regular topologies Muolo
_et al._ (2023), namely hypergraphs such that
$\mathbf{L}^{(d)}(t)=\alpha_{d}\mathbf{L}^{(1)}(t)$, for $d=1,\dots,D$, with
$\alpha_{1}=1$ and $\alpha_{d}\in\mathbb{R}_{+}$. Indeed we can use this
assumption to obtain from Eq. (9)

$\delta\dot{\vec{x}}_{i}=\mathbf{J}_{f}\delta\vec{x}_{i}+\sum_{j=1}^{n}L^{(1)}_{ij}(t)\mathbf{J}_{\hat{h}}\delta\vec{x}_{j}\,,$
(15)

where

$\mathbf{J}_{\hat{h}}:=\sum_{d=1}^{D}q_{d}\alpha_{d}\mathbf{J}_{h^{(d)}}\,,$
(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)}<0$,
$\psi^{(1)}(t)=(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 need to define a second $n\times n$ time
dependent matrix $\mathbf{b}(t)$ given by

$\frac{d\vec{\psi}^{(\alpha)}(t)}{dt}=\sum_{\beta}b_{\alpha\beta}(t)\vec{\psi}^{(\beta)}(t)\quad\forall\alpha=1,\dots,n\,,$
(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 Eq. (15) onto
$\psi^{(\alpha)}(t)$, we get

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

Let us conclude by observing that the latter equation has the same structure
of (14). Those equations determine the generalization of the Master Stability
Equation 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 $\mathbf{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 Pecora and Carroll (1998) and
higher-order structures Gambuzza  _et al._ (2021), also in the case of time-
varying higher-order structures, we can decouple the Master Stability Equation
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 Maximum
Lyapunov Exponent (MLE) associated with the transverse modes measures the
exponential growth rate of a tiny perturbation in the transverse subspace. It
serves as an enhanced form of Master Stability Function (MSF) and provides
valuable insights into the stability of the reference orbit. 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 Carletti and Fanelli (2022).

To demonstrate the above introduced theory and emphasize the outcomes arising
from the modified Master Stability Equations (14) and (18), we will present
two key examples in the following sections. Indeed, we will utilize the
Stuart-Landau limit cycle oscillator and the chaotic Lorenz system as
prototype dynamical systems anchored to each individual nodes. To simplify the
calculations, we assume that the hypergraph consists of only three nodes,
three links and one triangle (face), whose weights change in time.
Additionally, the eigenvector projection matrices $\mathbf{c}(t)$ and
$\mathbf{b}(t)$ do not vary in time; this assumption results from a suitable
choice of the Laplace eigenbasis as explained later in Appendix B. Finally, to
simplify the analysis we also assume the Laplace eigenvalues to be constant in
time. Let us stress that despite such assumptions, the proposed framework is
very general and can be applied to any time varying hypergraphs.

##  III Synchronization of Stuart-Landau oscillators coupled via time-varying
higher-order networks

The aim of this section is to present an application of the theory above
introduced. We decided to use the Stuart-Landau (SL) model as a prototype
example for two reasons; first, it provides the normal form for a generic
system close to a supercritical Hopf-bifurcation, second, because of its
structure, the Jacobian of the reaction part becomes constant once evaluated
on the reference orbit and this simplifies the presentation of the results.

A SL oscillator can be described by a complex amplitude $w$ that evolves in
time according to $\dot{w}=\sigma w-\beta|w|^{2}w$, where
$\sigma=\sigma_{\Re}+i\sigma_{\Im}$ and $\beta=\beta_{\Re}+i\beta_{\Im}$ are
complex model parameters. The system admits a limit cycle solution
$w_{LC}(t)=\sqrt{\sigma_{\Re}/\beta_{\Re}}e^{i\omega t}$, where
$\omega=\sigma_{\Im}-\beta_{\Im}\sigma_{\Re}/\beta_{\Re}$, that is stable
provided $\sigma_{\Re}>0$ and $\beta_{\Re}>0$, conditions that we hereby
assume.

To proceed in the analysis, we couple together $n$ identical SL oscillators,
each described by a complex amplitude $w_{j}$, with $j=1,...,n$, anchored to
the nodes of a time-varying hypergraph as prescribed in the previous section,
namely

$\frac{dw_{j}}{dt}=\sigma w_{j}-\beta
w_{j}|w_{j}|^{2}+\sum\limits_{d=1}^{D}q_{d}\sum\limits_{j_{1},\dots,j_{d}=1}^{n}A_{jj_{1}\dots
j_{d}}^{(d)}(t)\vec{g}^{(d)}(w_{j},w_{j_{1}},\dots,w_{j_{d}})\,.$ (19)

For the sake of simplicity, we restrict our analysis to pairwise and three-
body interactions, namely $D=2$ in Eq. (19). We hereby present and discuss the
SL synchronization under the diffusive-like coupling hypothesis and by using
two different assumptions: regular topology and natural coupling. The case of
non-invasive coupling will be presented in Appendix A.1.

###  III.1 Diffusive-like and regular topology

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{array}[]{l}g^{(1)}(w_{j},w_{j_{1}})=h^{(1)}(w_{j_{1}})-h^{(1)}(w_{j})\text{
and }\\\ \\\
g^{(2)}(w_{j},w_{j_{1}},w_{j_{2}})=h^{(2)}(w_{j_{1}},w_{j_{2}})-h^{(2)}(w_{j},w_{j})\,.\end{array}$

For the sake of definitiveness, let us fix

$h^{(1)}(w)=w\text{ and }h^{(2)}(w_{1},w_{2})=w_{1}w_{2}\,,$ (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 Eq. (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_{LC}(t)=\sqrt{\sigma_{\Re}/\beta_{\Re}}e^{i\omega t}$ by defining
$w_{j}=W_{LC}(1+\rho_{j})e^{i\theta_{j}}$, where $\rho_{j}$ and $\theta_{j}$
are real and small functions for all $j$. A straightforward computation allows
to write the time evolution of $\rho_{j}$ and $\theta_{j}$

$\dfrac{d}{dt}\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)+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)\,,$ (21)

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

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

$\dfrac{d}{dt}\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)+\Bigl{\\{}\left(\begin{matrix}-2\sigma_{\Re}&0\\\
-2\beta_{\Im}\frac{\sigma_{\Re}}{\beta_{\Re}}&0\end{matrix}\right)+\Lambda^{(\beta)}\left[\left(\begin{matrix}q_{1,\Re}&-q_{1,\Im}\\\
q_{1,\Im}&q_{1,\Re}\end{matrix}\right)\\\
+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)\right]\Bigr{\\}}\left(\begin{matrix}{\rho_{\beta}}\\\
{\theta_{\beta}}\end{matrix}\right)\,,$ (22)

where the matrix $\mathbf{b}$ has been defined in Eq. (17).

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=3$ nodes, three links and one triangle. Moreover, the eigenvalues are
assumed to be constant and the time-derivative of the associated eigenvectors
projected on the eigenbasis to return a constant matrix $\mathbf{b}$, for a
given $\Omega\geq 0$

$\mathbf{b}=\begin{pmatrix}0&0&0\\\ 0&0&\Omega\\\ 0&-\Omega&0\end{pmatrix}\,.$
(23)

One can show (see Appendix B and Carletti and Fanelli (2022)) that those
assumptions on the hypergraph correspond to two eigenvectors rotating in a
plane orthogonal to the constant eigenvector
$\psi^{(1)}\sim(1,\dots,1)^{\top}$ with frequency $\Omega>0$. The case
$\Omega=0$ corresponds thus to a static higher-order network structure.

Figure 1: Synchronization on time-varying regular higher-order network of
coupled SL oscillators. We report the MSF as a function of $q_{1,\Im}$ and
$q_{2,\Im}$ for two different values of $\Omega$, $\Omega=0$ (panel (a)) and
$\Omega=2$ (panel (b)), by using a color code, we determine the region of
stability (black) and the region of instability (yellow). The remaining
parameters have been fixed at the values $\alpha_{2}=2$, $\sigma=1.0+4.3i$,
$\beta=1.0+1.1i$, $q_{1,\Re}=0.1$, $q_{2,\Re}=0.1$, $\Lambda^{(2)}=-1$, and
$\Lambda^{(3)}=-2$.

Under those assumptions, Eq. (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 $[-5,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 Master Stability Function (MSF) of Eq. (22), as a
function of $q_{1,\Im}$ and $q_{2,\Im}$. The corresponding results are shown
in Fig. 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 MSF (black region), i.e., to the stability of the SL limit
cycle and thus to global synchronization, is larger for $\Omega>0$ than for
$\Omega=0$.

Figure 2: Synchronization on time-varying regular higher-order network of
coupled SL oscillators. The MSF is reported as a function of $\epsilon_{1}$
and $\epsilon_{2}$ for two different values of $\Omega$, $\Omega=0$ (panel
(a)) and $\Omega=2$ (panel (b)). The color code represents the values of the
MSF, negative values (black) while positive values (yellow). The remaining
parameters have been fixed at the values $\alpha_{2}=2$, $\sigma=1.0+4.3i$,
$\beta=1.0+1.1i$, $q_{1,0}=0.1-0.5i$, $q_{2,0}=0.1+0.5i$, $\Lambda^{(2)}=-1$,
and $\Lambda^{(3)}=-2$.

To study the combined effect of both coupling strengths $q_{1}$ and $q_{2}$,
we set $q_{1}=\epsilon_{1}q_{1,0}$ and $q_{2}=\epsilon_{2}q_{2,0}$, and we
compute the MSF as a function of $\epsilon_{1}$ and $\epsilon_{2}$, 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 Fig. 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.

Figure 3: Synchronization domains. We show the MSF in the plane
$(\Omega,\epsilon_{1})$ (panel (a)) for $\epsilon_{2}=0.02$ and in the plane
$(\Omega,\epsilon_{2})$ (panel (b)) for $\epsilon_{1}=0.02$. We can observe
that in both panels, the critical value of coupling strengths
$\hat{\epsilon}_{j}(\Omega)$ to achieve synchronization is smaller for
$\Omega>0$ than for $\Omega=0$. Furthermore, in panel (a) existence of an
interval $\mathcal{I}_{1}=[\Omega_{1},\Omega_{2}]$ can be observed such that
for all $\Omega\in\mathcal{I}_{1}$, there exist three different values of
critical coupling $\hat{\epsilon}_{1}$ for the occurrence of synchronization.
In panel (b), we can observe the existence of two intervals
$\mathcal{I}_{2}=[\Omega_{3},\Omega_{4}]$ and
$\mathcal{I}_{3}=[\Omega_{5},\Omega_{6}]$ such that for all
$\Omega\in\mathcal{I}_{2}$ there exist two critical values of
$\hat{\epsilon}_{2}$ and for all $\Omega\in\mathcal{I}_{3}$ there exist three
critical values of $\hat{\epsilon}_{2}$ for the emergence of synchronization.
The remaining parameters are kept fixed at the values $\alpha_{2}=2$,
$\sigma=1.0+4.3i$, $\beta=1.0+1.1i$, $q_{1,0}=0.1-0.5i$, $q_{2,0}=0.1+0.5i$,
$\Lambda^{(2)}=-1$, and $\Lambda^{(3)}=-2$.

Our last analysis concerns the relation between the frequency $\Omega$ and the
size of the coupling parameters $\epsilon_{1}$, $\epsilon_{2}$, still assuming
$q_{1}=\epsilon_{1}q_{1,0}$ and $q_{2}=\epsilon_{2}q_{2,0}$, on the onset of
synchronization. In Fig. 3 we report the MSF in the plane
$(\Omega,\epsilon_{1})$ for a fixed value of $\epsilon_{2}$ (panel (a)), and
in the plane $(\Omega,\epsilon_{2})$ for a fixed value of $\epsilon_{1}$
(panel (b)). Let us observe that the synchronization can be easier achieved
the smaller the value $\epsilon_{j}$, $j=1,2$, for which the MSF is negative,
having fixed $\Omega$. Let us thus define
$\hat{\epsilon}_{1}(\Omega)=\min\\{\epsilon>0:\mathrm{MSF}(\epsilon,\epsilon_{2},\Omega)<0\\}$,
for fixed $\epsilon_{2}$, and similarly $\hat{\epsilon}_{2}(\Omega)$. The
results of Fig. 3 clearly show that
$\hat{\epsilon}_{1}(\Omega)<\hat{\epsilon}_{1}(0)\sim 3.5$ and
$\hat{\epsilon}_{2}(\Omega)<\hat{\epsilon}_{2}(0)\sim 4.2$ and thus support
our claim that time-varying structures allow to achieve synchronization
easier.

To support our analysis, we performed numerical simulations of the SL defined
on the simple $3$ nodes time-varying hypergraph. We selected
$(\epsilon_{1},\epsilon_{2})=(2.5,0.5)$ and the remaining parameters values as
in Fig. 2. By observing the latter figure, we 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 Fig. 4 confirm
these conclusions; indeed, we can observe that (real part of) the complex
state variable is in phase for all $i$ in the case $\Omega=2$ (right panel),
while this is not clearly the case for $\Omega=0$ (left panel).

Figure 4: Temporal evolution of $\Re w_{i}$ for a suitable choice of coupling
parameters $\epsilon_{1}=2.5$ and $\epsilon_{2}=0.5$. In the top panel, we set
$\Omega=0$ while $\Omega=2$ in the bottom panel. The other parameters values
are the same as in Fig. 2, i.e., $\alpha_{2}=2$, $\sigma=1.0+4.3i$,
$\beta=1.0+1.1i$, $q_{1,0}=0.1-0.5i$, $q_{2,0}=0.1+0.5i$, $\Lambda^{(2)}=-1$,
and $\Lambda^{(3)}=-2$.

###  III.2 Diffusive-like and natural coupling

The aim of this section is to replace the condition of regular topology with a
condition of natural coupling and consider thus 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

$h^{(1)}(w)=h^{(2)}(w,w)\,.$

For the sake of definitiveness, let us fix

$h^{(1)}(w)=w^{3}\text{ and }h^{(2)}(w_{1},w_{2})=w_{1}(w_{2})^{2}\,.$ (24)

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

$\dfrac{d}{dt}\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(2\omega
t)&-\sin(2\omega t)\\\ \sin(2\omega t)&\cos(2\omega
t)\end{matrix}\right)\left(\begin{matrix}{\rho_{l}}\\\
{\theta_{l}}\end{matrix}\right)\,,$ (25)

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

By invoking the eigenvectors $\phi^{(\alpha)}(t)$ and eigenvalues
$\mu^{(\alpha)}(t)$ of $\mathbf{M}(t)$, and the matrix $\mathbf{c}$ (see Eq.
(13)), we can project the perturbation $\rho_{j}$ and $\theta_{j}$ on the
eigenbasis and thus rewrite the time variation of the perturbation as follows

$\dfrac{d}{dt}\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)+\biggl{[}\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^{(\beta)}\left(\begin{matrix}\cos(2\omega
t)&-\sin(2\omega t)\\\ \sin(2\omega t)&\cos(2\omega
t)\end{matrix}\right)\biggr{]}\left(\begin{matrix}{\rho_{\beta}}\\\
{\theta_{\beta}}\end{matrix}\right)\,.$ (26)

Figure 5: Synchronization on time-varying higher-order network of coupled SL
oscillators with diffusive-like natural coupling. We report the MSF as a
function of the eigenvalues $\mu^{(2)}$ and $\mu^{(3)}$ for two different
choices of $\Omega$, $\Omega=0$ (panel (a)) and $\Omega=2$ (panel (b)) by
using a color code, black is associated to negative values while positive ones
are shown in yellow. We characterize the range of the axes by considering the
absolute values of the eigenvalues. The remaining parameters are kept fixed at
$\sigma=1.0+4.3i$, $\beta=1.0+1.1i$.

Let us assume again to deal with an hypergraph made by $3$ nodes and consider
a time-independent matrix $\mathbf{c}$

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

for some $\Omega\geq 0$. The eigenvalue $\mu^{(1)}=0$ of $\mathbf{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 the latter equations and provide the conditions for a
global stable synchronous solution to exist. In Fig. 5, 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 Eq. (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 easier
synchronize.

##  IV Synchronization of Lorenz systems nonlinearly coupled via time-varying
higher-order networks

The aim of this section is to show that our results hold true beyond the
example of the dynamical system shown above, i.e., the Stuart-Landau. We thus
decide to present an application of synchronization for chaotic systems on a
time-varying higher-order network. For the sake of definitiveness, we used the
paradigmatic chaotic Lorenz model for the evolution of individual nonlinear
oscillators.

We consider again the scenario of regular topology with the toy model
hypergraph structure composed of $n=3$ nodes described previously, the whole
system can thus be described by

$\begin{cases}\dot{x}_{i}&=a_{1}(y_{i}-x_{i})+\epsilon_{2}\sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N}A^{(2)}_{ijk}(x_{j}^{2}x_{k}-x_{i}^{3})\\\
\dot{y}_{i}&=x_{i}(a_{3}-z_{i})-y_{i}+\epsilon_{1}\sum\limits_{j=1}^{N}A^{(1)}_{ij}(y_{j}-y_{i})\\\
\dot{z}_{i}&=x_{i}y_{i}-a_{2}z_{i}\end{cases}\,,$ (27)

where the system parameters are kept fixed at $a_{1}=10$, $a_{2}=\frac{8}{3}$,
$a_{3}=28$ for which individual nodes exhibits chaotic trajectory. The
pairwise and higher-order structures are related to each other by
$\mathbf{L}^{(2)}=\alpha_{2}\mathbf{L}^{(1)}$. We assume the eigenvalues of
the Laplacian $\mathbf{L}^{(1)}$ to be constant and the matrix $\mathbf{b}$ to
be given by

$\mathbf{b}=\begin{pmatrix}0&0&0\\\ 0&0&\Omega\\\
0&-\Omega&0\end{pmatrix}\quad\text{for some $\Omega\geq 0$.}$

Let us thus select as reference solution $\vec{s}(t)$ a chaotic orbit of the
isolated Lorenz model and consider as done previously the time evolution of a
perturbation about such trajectory. Computations similar to those reported
above, allow to obtain a linear non-autonomous system ruling the evolution of
the perturbation, whose stability can be numerically inferred by computing the
largest Lyapunov exponent, i.e., the MSF. We first considered the impact of
the coupling strength, $\epsilon_{1}$ and $\epsilon_{2}$ on synchronization;
results are reported in Fig. 6 where we present the level sets of the MSF as a
function of the above parameters by using a color code: black dots refer to
negative MSF while yellow dots to positive MSF. The panel (a), refers to a
static hypergraph, i.e., $\Omega=0$, while the panel (b) to a time-varying
one, i.e., $\Omega=3$, one can thus appreciate that the latter setting allows
a negative MSF for a larger range of parameters $\epsilon_{1}$ and
$\epsilon_{2}$ and hence we can conclude that time-varying hypergraph enhance
synchronization also in the case of chaotic oscillators.

Figure 6: Synchronization on time-varying regular higher-order network of
coupled Lorenz oscillators. We report the MSF as a function of the coupling
strengths, $\epsilon_{1}$ and $\epsilon_{2}$, for two different values of
$\Omega$, $\Omega=0$ (panel (a)) and $\Omega=3$ (panel (b)), by using a color
code, where black dots stand for a negative MSF, i.e., global synchronization,
while yellow dots for a positive MSF. The remaining parameters are kept fixed
at $a_{1}=10$, $a_{2}=\frac{8}{3}$, $a_{3}=28$, and $\alpha_{2}=2$.

We conclude this analysis by studying again the relation between the frequency
$\Omega$ and the size of the coupling parameters $\epsilon_{1}$,
$\epsilon_{2}$ on the onset of synchronization. In Fig. 7 we show the MSF in
the plane $(\Omega,\epsilon_{1})$ for a fixed value of $\epsilon_{2}=0.01$
(panel (a)), and in the plane $(\Omega,\epsilon_{2})$ for a fixed value of
$\epsilon_{1}=0.2$ (panel (b)). By using again
$\hat{\epsilon}_{1}(\Omega)=\min\\{\epsilon>0:\mathrm{MSF}(\epsilon,\epsilon_{2},\Omega)<0\\}$,
for fixed $\epsilon_{2}$, and similarly $\hat{\epsilon}_{2}(\Omega)$, we can
conclude that $\hat{\epsilon}_{1}(\Omega)<\hat{\epsilon}_{1}(0)\sim 1.4$ and
$\hat{\epsilon}_{2}(\Omega)<\hat{\epsilon}_{2}(0)\sim 0.04$ and thus
supporting again our claim that time-varying structures allow to achieve
synchronization easier.

Figure 7: We show the MSF in the plane $(\Omega,\epsilon_{1})$ (panel (a)) for
$\epsilon_{2}=0.01$) and in the plane $(\Omega,\epsilon_{2})$ (panel (b)) for
$\epsilon_{1}=0.2$). We can observe that in both panels the critical value of
coupling strengths $\hat{\epsilon}_{j}(\Omega)$ to achieve synchronization is
smaller for $\Omega>0$ than for $\Omega=0$. This implies that synchronization
can occur more easily on a time-varying higher-order structure than on a
static one.

##  V Conclusions

To sum up we have here introduced and studied a generalized framework for the
emergence of global synchronization on time-varying higher-order networks and
developed a theory for its stability without imposing strong restrictions on
the functional time evolution of the higher-order structure. We have
demonstrated that the latter can be examined by extending the Master Stability
Function technique to the novel framework for specific cases based either on
the inter-node coupling scheme or the topology of the higher-order structure.
Our findings reveal that the behavior of the higher-order network is
represented by a matrix that changes over time and possesses skew symmetry.
This matrix is derived from the time-dependent evolution of the eigenvectors
of the higher-order Laplacian. Additionally, the eigenvalues associated with
these eigenvectors can also vary over time and have an impact on shaping the
evolution of the introduced disturbance. We have validated the proposed theory
on time-varying hypergraphs of coupled Stuart-Landau oscillators and chaotic
Lorenz systems, and the results obtained indicate that incorporating temporal
aspects into group interactions can facilitate synchronization in higher-order
networks compared to static ones.

The framework and concepts presented in this study create opportunities for
future research on the impact of temporality in systems where time-varying
group interactions have been observed but not yet thoroughly explored due to
the absence of a suitable mathematical setting. Importantly, the fact that our
theory does not require any restrictions on the time evolution of the
underline structure could offer the possibility to apply it for a diverse
range of applications other than synchronization.

## References

  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Physics Reports 469, 93 (2008).
  * Boccaletti  _et al._ (2018) S. Boccaletti, A. N. Pisarchik, C. I. Del Genio, and A. Amann,  _Synchronization: from coupled systems to complex networks_ (Cambridge University Press, 2018).
  * Barabási (2016) A.-L. Barabási,  _Network science_ (Cambridge University Press, 2016).
  * Wasserman  _et al._ (1994) S. Wasserman, K. Faust, _et al._ ,  _Social Network Analysis: Methods and Applications_ (Cambridge University Press, 1994).
  * Valencia  _et al._ (2008) M. Valencia, J. Martinerie, S. Dupont, and M. Chavez, Physical Review E 77, 050905 (2008).
  * Bassett  _et al._ (2011) D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proceedings of the National Academy of Sciences 108, 7641 (2011).
  * Holme and Saramäki (2012) P. Holme and J. Saramäki, Physics Reports 519, 97 (2012).
  * Masuda and Lambiotte (2016) N. Masuda and R. Lambiotte,  _A guide to temporal networks_ (World Scientific, 2016).
  * Ghosh  _et al._ (2022) D. Ghosh, M. Frasca, A. Rizzo, S. Majhi, S. Rakshit, K. Alfaro-Bittner, and S. Boccaletti, Physics Reports 949, 1 (2022).
  * Carletti and Fanelli (2022) T. Carletti and D. Fanelli, Chaos, Solitons & Fractals 159, 112180 (2022).
  * Anwar  _et al._ (2022) M. S. Anwar, S. Rakshit, D. Ghosh, and E. M. Bollt, Physical Review E 105, 024303 (2022).
  * Carletti  _et al._ (2020a) T. Carletti, D. Fanelli, and S. Nicoletti, Journal of Physics Complexity 1, 035006 (2020a).
  * Battiston  _et al._ (2020) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Physics Reports 874, 1 (2020).
  * Battiston  _et al._ (2021) F. Battiston, E. Amico, A. Barrat, G. Bianconi, G. Ferraz de Arruda, B. Franceschiello, I. Iacopini, S. Kéfi, V. Latora, Y. Moreno, _et al._ , Nature Physics 17, 1093 (2021).
  * Majhi  _et al._ (2022) S. Majhi, M. Perc, and D. Ghosh, Journal of the Royal Society Interface 19, 20220043 (2022).
  * Boccaletti  _et al._ (2023) S. Boccaletti, P. De Lellis, C. del Genio, K. Alfaro-Bittner, R. Criado, S. Jalan, and M. Romance, Physics Reports 1018, 1 (2023).
  * Berge (1973) C. Berge,  _Graphs and hypergraphs_ , North-Holland Pub. Co. (American Elsevier Pub. Co, 1973).
  * Bianconi (2021) G. Bianconi,  _Higher-order networks: An introduction to simplicial compelxes_ (Cambridge University Press, 2021).
  * Neuhäuser  _et al._ (2020) L. Neuhäuser, A. Mellor, and R. Lambiotte, Physical Review E 101, 032310 (2020).
  * Neuhäuser  _et al._ (2021) L. Neuhäuser, R. Lambiotte, and M. T. Schaub, Physical Review E 104, 064305 (2021).
  * Carletti  _et al._ (2020b) T. Carletti, F. Battiston, G. Cencetti, and D. Fanelli, Phys. Rev. E 101, 022308 (2020b).
  * Schaub  _et al._ (2020) M. T. Schaub, A. R. Benson, P. Horn, G. Lippner, and A. Jadbabaie, SIAM Review 62, 353 (2020).
  * Muolo  _et al._ (2023) R. Muolo, L. Gallo, V. Latora, M. Frasca, and T. Carletti, Chaos, Solitons & Fractals 166, 112912 (2023).
  * Gao  _et al._ (2023) S. Gao, L. Chang, M. Perc, and Z. Wang, Physical Review E 107, 014216 (2023).
  * Skardal and Arenas (2020) P. S. Skardal and A. Arenas, Communications Physics 3, 1 (2020).
  * Skardal and Arenas (2019) P. S. Skardal and A. Arenas, Physical Review Letters 122, 248301 (2019).
  * Carletti  _et al._ (2023) T. Carletti, L. Giambagli, and G. Bianconi, Physical Review Letters 130, 187401 (2023).
  * Anwar and Ghosh (2022a) M. S. Anwar and D. Ghosh, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 033125 (2022a).
  * Anwar and Ghosh (2022b) M. S. Anwar and D. Ghosh, Physical Review E 106, 034314 (2022b).
  * Iacopini  _et al._ (2019) I. Iacopini, G. Petri, A. Barrat, and V. Latora, Nature Communications 10, 2485 (2019).
  * Chowdhary  _et al._ (2021) S. Chowdhary, A. Kumar, G. Cencetti, I. Iacopini, and F. Battiston, Journal of Physics: Complexity 2, 035019 (2021).
  * Cencetti  _et al._ (2021) G. Cencetti, F. Battiston, B. Lepri, and M. Karsai, Scientific Reports 11, 7028 (2021).
  * Anwar and Ghosh (2022c) M. S. Anwar and D. Ghosh, arXiv preprint arXiv:2212.01081 (2022c).
  * Stilwell  _et al._ (2006) D. J. Stilwell, E. M. Bollt, and D. G. Roberson, SIAM Journal on Applied Dynamical Systems 5, 140 (2006).
  * Petit  _et al._ (2017) J. Petit, B. Lauwens, D. Fanelli, and T. Carletti, Physical Review Letters 119, 148301 (2017).
  * Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, Physical Review Letters 80, 2109 (1998).
  * Van Gorder (2021) R. A. Van Gorder, Proceedings of the Royal Society A 477, 20200753 (2021).
  * Zhang and Strogatz (2021) Y. Zhang and S. H. Strogatz, Nature Communications , 3273 (2021).
  * Gambuzza  _et al._ (2021) L. V. Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora, and S. Boccaletti, Nature Communications 12, 1255 (2021).

##  Appendix A Non-invasive couplings

Here we will discuss the results corresponding to a slightly more general
hypothesis for $\vec{g}^{(d)}$, namely to be non-invasive, i.e.,

$\vec{g}^{(d)}(\vec{s},\dots,\vec{s})=0\quad\forall d=1,\dots,D\,,$ (28)

whose goal is again to guarantee that the coupling term in Eq. (3) vanishes
once evaluated on the orbit $(\vec{s}(t),\dots,\vec{s}(t))^{\top}$. Indeed by
using again $\vec{x}_{i}=\vec{s}+\delta\vec{x}_{i}$ and expanding Eq. (3) up
to the first order we get

$\delta\dot{\vec{x}}_{i}=\mathbf{J}_{f}\delta\vec{x}_{i}+\sum_{d=1}^{D}q_{d}\sum_{j_{1},\dots,j_{d}=1}^{n}B_{ij_{1}\dots
j_{d}}(t)\left[\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{i}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}\delta\vec{x}_{i}+\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{1}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}\delta\vec{x}_{j_{1}}+\dots+\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{d}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}\delta\vec{x}_{j_{d}}\right]\,;$
(29)

from Eq. (28) we can obtain

$\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{i}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}+\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{1}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}+\dots+\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{d}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}=0\,,$

and thus rewrite (29) as follows

$\delta\dot{\vec{x}}_{i}=\mathbf{J}_{f}\delta\vec{x}_{i}+\sum_{d=1}^{D}q_{d}\sum_{j_{1},\dots,j_{d}=1}^{n}B_{ij_{1}\dots
j_{d}}(t)\left[\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{1}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}(\delta\vec{x}_{j_{1}}-\delta\vec{x}_{i})+\dots+\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{d}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}(\delta\vec{x}_{j_{d}}-\delta\vec{x}_{i})\right]\,.$
(30)

Recalling the definition of $k^{(d)}_{ij}$ given in Eq. (6) we get

$\delta\dot{\vec{x}}_{i}=\mathbf{J}_{f}\delta\vec{x}_{i}+\sum_{d=1}^{D}q_{d}(d-1)!\left[\sum_{j_{1}=1}^{n}k^{(d)}_{ij_{1}}(t)\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{1}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}(\delta\vec{x}_{j_{1}}-\delta\vec{x}_{i})+\dots+\sum_{j_{l}=1}^{n}k^{(d)}_{ij_{d}}(t)\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{d}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}(\delta\vec{x}_{j_{d}}-\delta\vec{x}_{i})\right]\,.$
(31)

By using the definition of the higher-order Laplace matrix (7) we eventually
obtain

$\delta\dot{\vec{x}}_{i}=\mathbf{J}_{f}\delta\vec{x}_{i}-\sum_{d=1}^{D}q_{d}\sum_{j=1}^{n}L^{(d)}_{ij}(t)\left[\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{1}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}+\dots+\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{d}}}\Big{\rvert}_{(\vec{s},\dots,\vec{s})}\right]\delta\vec{x}_{j}\,.$
(32)

Let us consider now a particular case of non-invasive function, we assume thus
there exists a function
$\vec{\varphi}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$, such that
$\vec{\varphi}(0)=0$ and define

$g^{(d)}(\vec{x}_{i},\vec{x}_{j_{1}},\dots,\vec{x}_{j_{d}})=\sum_{\ell=1}^{d}\vec{\varphi}(\vec{x}_{i}-\vec{x}_{j_{\ell}})\,,$
(33)

then

$\frac{\partial\vec{g}^{(d)}}{\partial\vec{x}_{j_{\ell}}}=-\mathbf{J}_{\varphi}(0)\,,$

where $\mathbf{J}_{\varphi}(0)$ is the Jacobian of the function
$\vec{\varphi}$ evaluated at $0$. In conclusion (32) rewrites as follows

$\delta\dot{\vec{x}}_{i}=\mathbf{J}_{f}\delta\vec{x}_{i}-\sum_{d=1}^{D}q_{d}\sum_{j=1}^{n}L^{(d)}_{ij}(t)(-d)\mathbf{J}_{\varphi}(0)\delta\vec{x}_{j}=\mathbf{J}_{f}\delta\vec{x}_{i}+\sum_{j=1}^{n}G_{ij}(t)\mathbf{J}_{\varphi}(0)\delta\vec{x}_{j}\,,$
(34)

where $\mathbf{G}(t)=\sum_{d=1}^{D}dq_{d}\mathbf{L}^{(d)}(t)$ can be
considered as an effective time-varying simplicial complex or hypergraph.

Let us now observe that the effective matrix $\mathbf{G}(t)$ is a Laplace
matrix; it is non-positive definite (as each one of the $\mathbf{L}^{(d)}(t)$
does for any $d=1,\dots,D$ and any $t>0$), it admits $\mu^{(1)}=0$ as
eigenvalue associated to the eigenvector $\phi^{(1)}=(1,\dots,1)^{\top}$ and
it is symmetric. So there exist a orthonormal time-varying eigenbasis,
$\phi^{(\alpha)}(t)$, $\alpha=1,\dots,n$, for $\mathbf{G}(t)$ with associated
eigenvalues $\mu^{(\alpha)}\leq 0$. Similar to before, we define the $n\times
n$ time dependent matrix $\mathbf{c}(t)$ that quantifies the projections of
the time derivatives of the eigenvectors onto the independent eigendirections,
namely

$\frac{d\vec{\phi}^{(\alpha)}}{dt}(t)=\sum_{\beta}c_{\alpha\beta}(t)\vec{\phi}^{(\beta)}(t)\quad\forall\alpha=1,\dots,n\,.$
(35)

By recalling the orthonormality condition
$\left(\vec{\phi}^{(\alpha)}(t)\right)^{\top}\cdot\vec{\phi}^{(\beta)}(t)=\delta_{\alpha\beta}$
we can again straightforwardly conclude that $\mathbf{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$.

Thereafter, we consider Eq. (34), 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 $\mathbf{c}$ we obtain

$\frac{d\delta\hat{\vec{x}}_{\beta}}{dt}=\sum_{\alpha}c_{\beta\alpha}(t)\delta\hat{\vec{x}}_{\alpha}+\left[\mathbf{J}_{f}+\mu^{(\beta)}(t)\mathbf{J}_{\varphi}(0)\right]\delta\hat{\vec{x}}_{\beta}\,.$
(36)

This is the required Master Stability Equation, solving which for the
calculation of maximum Lyapunov exponents provide the condition for stability
of the synchronous solution.

###  A.1 Synchronization of Stuart-Landau oscillators with non-invasive
coupling assumption

Figure 8: Synchronization on time-varying higher-order network of coupled SL
oscillators with non-invasive coupling configuration. Region of synchrony and
desynchrony are depicted by simultaneously varying $\mu^{(2)}$ and $\mu^{(3)}$
for two different values of $\Omega$ (a) $\Omega=0$, (b) $\Omega=2$, where the
domain in black indicates the area of the stable synchronous solution. The
range of the axes is characterized by considering the absolute values of the
eigenvalues. All the other values are kept fixed at $\sigma=1.0+4.3i$,
$\beta=1.0+1.1i$, $\varphi^{\prime}(0)=1$.

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 $\varphi$ such that $\varphi(0)=0$, $\varphi^{\prime}(0)\neq 0$
and

$\begin{array}[]{l}g^{(1)}(w_{1},w_{2})=\varphi(w_{1}-w_{2})\,,\text{ and }\\\
g^{(2)}(w_{1},w_{2},w_{3})=\varphi(w_{1}-w_{2})+\varphi(w_{1}-w_{3})\,.\end{array}$
(37)

By reasoning as before, we get

$\begin{array}[]{l}\dfrac{d}{dt}\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}(0)\sum_{\ell}\left(q_{1}L^{(1)}_{j\ell}+q_{2}L^{(2)}_{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}$ (38)

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

$\begin{array}[]{l}\dfrac{d}{dt}\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)+\biggl{[}\left(\begin{matrix}-2\sigma_{\Re}&0\\\
-2\beta_{\Im}\frac{\sigma_{\Re}}{\beta_{\Re}}&0\end{matrix}\right)+\varphi^{\prime}(0)\mu^{(\beta)}\left(\begin{matrix}1&0\\\
0&-1\end{matrix}\right)\biggr{]}\left(\begin{matrix}{\rho_{\beta}}\\\
{\theta_{\beta}}\end{matrix}\right).\end{array}$ (39)

Figure 8 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 $\mathbf{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 Eq. (39) 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.

##  Appendix B Structure of the small hypergraph

The goal of this section is to provide more details about the construction of
the simple time-varying hypergraph used as support for the numerical
simulations in the main text. To start with we need to obtain the time-
evolution of eigenvectors $\vec{\psi}^{(\alpha)}(t)$, which follows the
equation

$\begin{array}[]{l}\dfrac{d\vec{\psi}^{(\alpha)}}{dt}=\sum\limits_{\alpha}b_{\beta\alpha}\vec{\psi}^{(\alpha)}\,,\end{array}$
(40)

where the matrix $\mathbf{b}$ has been given in Eq. (23). The eigenvector
associated with the least eigenvalue $\Lambda^{(1)}=0$ is constant and is
given by $\vec{\psi}^{(1)}=\frac{1}{\sqrt{3}}(1,1,1)^{\top}$. The other two
eigenvectors are obtained by solving the previous equation and are represented
as $\vec{\psi}^{(2)}(t)=\vec{v}_{1}\cos(\Omega t)+\vec{v}_{2}\sin(\Omega t)$
and $\vec{\psi}^{(3)}(t)=-\vec{v}_{1}\sin(\Omega t)+\vec{v}_{2}\cos(\Omega
t)$, where $\vec{v}_{1}$, $\vec{v}_{2}$ are the unknown vectors that should be
determined using the constraints to have orthonormal eigenbasis for every $t$.
Following a few steps of calculation, we can obtain the other two eigenvectors
as follows

$\begin{array}[]{l}\vec{\psi}^{(2)}(t)=\dfrac{1}{\sqrt{6}}\begin{pmatrix}1\\\
-2\\\ 1\end{pmatrix}\cos(\Omega t)+\dfrac{1}{\sqrt{2}}\begin{pmatrix}-1\\\
0\\\ 1\end{pmatrix}\sin(\Omega t)\;\;\mbox{and},\\\
\vec{\psi}^{(3)}(t)=-\dfrac{1}{\sqrt{6}}\begin{pmatrix}1\\\ -2\\\
1\end{pmatrix}\sin(\Omega t)+\dfrac{1}{\sqrt{2}}\begin{pmatrix}-1\\\ 0\\\
1\end{pmatrix}\cos(\Omega t).\end{array}$ (41)

Now recalling our assumption about constant eigenvalues and using the relation
$\mathbf{L}^{(1)}_{ij}(t)=\sum\limits_{\alpha}\Lambda^{(\alpha)}\vec{\psi}^{(\alpha)}_{i}(t)\vec{\psi}^{(\alpha)}_{j}(t)$,
we can obtain the entries of the pairwise Laplace matrix as

$\begin{array}[]{l}{L}^{(1)}_{ij}(t)=\Lambda^{(2)}\vec{\psi}^{(2)}_{i}(t)\vec{\psi}^{(2)}_{j}(t)+\Lambda^{(3)}\vec{\psi}^{(3)}_{i}(t)\vec{\psi}^{(3)}_{j}(t),\end{array}$
(42)

where we use the fact that $\Lambda^{(1)}=0$ for all time $t$. Finally by
using the relation between pairwise adjacency and Laplace matrices
$L^{(1)}_{ij}(t)=A^{(1)}_{ij}(t)$, for $i\neq j$, we obtain the temporal
evolution of the links as

$\begin{array}[]{l}A^{(1)}_{12}(t)=\dfrac{1}{2}-\dfrac{1}{3}\cos(\frac{\pi}{3}+2\Omega
t),\\\ \\\ A^{(1)}_{13}(t)=\dfrac{1}{2}+\dfrac{1}{3}\cos(2\Omega t),\\\ \\\
A^{(1)}_{23}(t)=\dfrac{1}{2}-\dfrac{1}{3}\cos(\frac{\pi}{3}-2\Omega
t),\end{array}$ (43)

where we have used the fact that the non-zero eigenvalues are given by
$\Lambda^{(2)}=-1$ and $\Lambda^{(3)}=-2$.

Again from the regular structure of the hypergraph, we have
$\mathbf{L}^{(2)}(t)=\alpha_{2}\mathbf{L}^{(1)}(t)$, for all $t$. Therefore,
following the relation (42), entries of the $2nd$-order Laplacian
$\mathbf{L}^{(2)}$ can be represented as,

$\begin{array}[]{l}L^{(2)}_{ij}(t)=\alpha_{2}[\Lambda^{(2)}\vec{\psi}^{(2)}_{i}(t)\vec{\psi}^{(2)}_{j}(t)+\Lambda^{(3)}\vec{\psi}^{(3)}_{i}(t)\vec{\psi}^{(3)}_{j}(t)].\end{array}$
(44)

Now, the definition of higher-order Laplacian implies that,
$L^{(2)}_{ij}(t)=\sum\limits_{k}A^{(2)}_{ijk}(t)$, $i\neq j$. Hence, using the
above relation and Eq. (44), we can obtain the temporal evolution of the
$3$-hyperedge as

$\begin{array}[]{l}A^{(2)}_{123}(t)=1-\dfrac{2}{3}\cos(\frac{\pi}{3}+2\Omega
t),\end{array}$ (45)

where we have again used the fact that the non-zero eigenvalues are
$\Lambda^{(2)}=-1$, and $\Lambda^{(3)}=-2$, and the value of the parameter
$\alpha_{2}$ has been set $\alpha_{2}=2$. Due to the assumption of undirected
hypergraph, we also trivially have,
$A^{(2)}_{123}(t)=A^{(2)}_{\pi{(123)}}(t)$, where $\pi{(123)}$ indicates any
permutation of $(123)$. Fig. 9 portrays the temporal evolution of the links
and $3$-hyperedge 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)}_{123}(t)$.

Figure 9: Temporal evolution of edges and $3$-hyperedge obtained from Eqs.
(43) and (45) for a particular value of $\Omega=2$.