Global synchronization on time-varying higher-order structures

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 [1].This phenomenon, observed in a wide range of human-made and natural systems [2], 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 [3].However, the classical static network representation has its limitation in modeling many empirical systems, such as social networks [4], brain networks [5,6], 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 [7,8], 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 [9][10][11].
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 [12][13][14][15][16].By considering the simultaneous interactions of many agents, higher-order structures, namely hypergraphs [17] and simplicial complexes [18], 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 [19,20], random walks [21,22], pattern formation [12,23,24], synchronization [12,[25][26][27][28][29], social contagion and epidemics [30,31].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 [32].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 [20] and to predict the onset of endemic state in epidemic spreading [31].
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 [33], that investigates synchronization in time-varying simplicial complexes, limited only to fast switching [34,35] 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 [36] 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 where x ∈ R m denotes the state vector and f : R m → 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 where for d = 1, . . ., D, q d > 0 denotes the coupling strength, g (d) : R (d+1)m → R m the nonlinear coupling function and A (d) (t) the tensor encoding which units are interacting together.More precisely A (d) ij1...j d (t) = 1 if the units i, j 1 , . . ., 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 π(ij1...j d ) (t) = 1 for any permutation π of the indexes i, j 1 , . . ., 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, s(t), of system (1).We are interested in determining the conditions under which the orbit ( s(t), . . ., s(t)) ⊤ 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., g (d) ( s, . . ., s) = 0, for d = 1, . . ., 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 g (d) to be diffusive-like, namely for each d there exists a function h (d) : R dm → R m such that (3) In this way we can straightforwardly ensure that the coupling term in Eq. ( 3) vanishes once evaluated on the orbit ( s(t), . . ., s(t)) ⊤ , 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 ( s(t), . . ., s(t)) ⊤ with a spatially inhomogeneous term, meaning that ∀i ∈ {1, . . ., n} we define x i = s + δ x i .Substituting the latter into Eq.( 2) and expanding up to the first order, we obtain where being δ ij1j2...jD the generalized multi-indexes Kroneckerδ, and the (time-varying) d-degree of node i is given by which represents the number of hyperedges of order d incident to node i at time t.Observe that if i (t) counts both the number and the weight, it is thus the generalization of the strength of a node.Let us now define namely the number of hyperedges of order d containing both nodes i and j at time t.Again, once A (d) is weighted, then k ji (t).Finally, we define the generalized time-varying higher-order Laplacian matrix for the interaction of order d as Observe that such a matrix is symmetric because of the assumption of the tensors 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 where we used the fact the ∂ h (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 we can rewrite Eq. ( 8) in compact form This is a non-autonomous linear differential equation determining the stability of the perturbation δ 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 h (d) satisfy the condition of natural coupling (see Section II A) or the higher-order structures exhibit regular topologies (see Section II B).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.

A. Natural coupling
Let us assume the functions h (d) to satisfy the condition of natural coupling, namely and it allows to eventually rewrite Eq. ( 9) as follows where Let us observe that the matrix M(t) is a Laplace matrix; it is non-positive definite (as each one of the L (d) (t) matrices does for any d = 1, . . ., D and any t > 0, and q d > 0), it admits µ (1) = 0 as eigenvalue associated to the eigenvector φ (1) = (1, . . ., 1) ⊤ and it is symmetric.So there exists an orthonormal time-varying eigenbasis, φ (α) (t), α = 1, . . ., n, for M(t) with associated eigenvalues µ (α) ≤ 0. Let us define [10] the n × n time dependent matrix c(t) that quantifies the projections of the time derivatives of the eigenvectors onto the independent eigendirections, namely By recalling the orthonormality condition we can straightforwardly conclude that c is a real skewsymmetric matrix with a null first row and first column, i.e., c αβ + c βα = 0 and c 1α = 0.
To make one step further, we consider Eq. ( 11), and we project it onto the eigendirections, namely we introduce and recalling the definition of c we obtain Let us observe that the latter formula and the following analysis differ from the one presented in [37] 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 c = 0.The same assumption is also at the root of the results by [38]; 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.

B. Regular topologies
An alternative approach to study Eq. ( 9) is to assume regular topologies [23], namely hypergraphs such that L (d) (t) = α d L (1) (t), for d = 1, . . ., D, with α 1 = 1 and α d ∈ R + .Indeed we can use this assumption to obtain from Eq. ( 9) where that results in a sort of weighted nonlinear coupling term.We can now make use of the existence of a timevarying orthonormal basis of L (1) (t), namely ψ (α) (t), α = 2, . . ., n, associated to eigenvalues Λ (α) < 0, ψ (1) (t) = (1, . . ., 1) ⊤ and Λ (1) = 0, to project δ x i onto the n eigendirections, i .Because the latter vary in time we need to define a second n × n time dependent matrix b(t) given by that it is again real, skew-symmetric, with a null first row and first column, i.e., b αβ + b βα = 0 and b 1α = 0, because of the orthonormality condition of eigenvectors.By projecting Eq. ( 15) onto ψ (α) (t), we get 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 c(t) or b(t) and the eigenvectors µ (α) (t) or Λ (α) (t), while the dynamics (resp.the coupling) in the Jacobian J f (resp.J h (1) or J ĥ).
It is important to notice that as the eigenvalues µ (1) = 0, Λ (1) = 0 and the skew-symmetric matrices c(t), b(t) have null first row and column, in analogy with the MSF approaches carried over static networks [36] and higherorder structures [39], 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 c(t) and 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 [10].
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 c(t) and 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.
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 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.
For the sake of definitiveness, let us fix let us observe that the latter functions do not satisfy the condition for natural coupling, indeed h (1) (w) = w = w 2 = h (2) (w, w).
Let us assume to deal with regular topology, namely L (2) = α 2 L (1) .Hence following Eq.( 16) we can define J ĥ = q 1 J h (1) + q 2 α 2 J h (2) .Let us perturb the limit cycle solution w LC (t) = σ ℜ /β ℜ e iωt by defining w j = W LC (1 + ρ j )e iθ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 d dt where ω = σ ℑ − β ℑ σ ℜ /β ℜ is the frequency of the limit cycle solution.
By exploiting the eigenvectors ψ (α) (t) and eigenvalues Λ (α) (t) of L (1) (t) to project the perturbation ρ j and θ j we obtain: where the matrix 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 b, for a given Ω ≥ 0 One can show (see Appendix B and [10]) that those assumptions on the hypergraph correspond to two eigenvectors rotating in a plane orthogonal to the constant eigenvector ψ (1) ∼ (1, . . ., 1) ⊤ with frequency Ω > 0.
The case Ω = 0 corresponds thus to a static higher-order network structure.

FIG. 1.
Synchronization on time-varying regular higher-order network of coupled SL oscillators.We report the MSF as a function of q 1,ℑ and q 2,ℑ for two different values of Ω, Ω = 0 (panel (a)) and Ω = 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 α2 = 2, σ = 1.0 + 4.3i, β = 1.0 + 1.1i, q 1,ℜ = 0.1, q 2,ℜ = 0.1, Λ (2) = −1, and Λ 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,ℑ and q 2,ℑ 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,ℑ and q 2,ℑ .The corresponding results are shown in Fig. 1 for Ω = 0 (panel (a)) and Ω = 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 Ω > 0 than for Ω = 0. To study the combined effect of both coupling strengths q 1 and q 2 , we set q 1 = ǫ 1 q 1,0 and q 2 = ǫ 2 q 2,0 , and we compute the MSF as a function of ǫ 1 and ǫ 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 (Ω = 0, panel (a)) and time-varying (Ω = 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.
To support our analysis, we performed numerical simulations of the SL defined on the simple 3 nodes timevarying hypergraph.We selected (ǫ 1 , ǫ 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 Ω = 0 and negative if Ω = 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 Ω = 2 (right panel), while this is not clearly the case for Ω = 0 (left panel).

B. 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 For the sake of definitiveness, let us fix Consider again to perturb the limit cycle solution w LC (t) = σ ℜ /β ℜ e iωt by defining w j = W LC (1+ρ j )e iθ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, where ω = σ ℑ − β ℑ σ ℜ /β ℜ is the frequency of the limit cycle solution and M is the matrix q 1 L (1) (t) + q 2 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 φ (α) (t) and eigenvalues µ (α) (t) of M(t), and the matrix c (see Eq. ( 13)), we can project the perturbation ρ j and θ j on the eigenbasis and thus rewrite the time variation of the perturbation as follows FIG. 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 µ (2) and µ (3)  for some Ω ≥ 0. The eigenvalue µ (1) = 0 of M determines the dynamics parallel to the synchronous manifold.
On the other hand, the equations obtained for µ (2) and µ (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 µ (2) and µ (3) while keeping the remaining parameters in Eq. ( 26) fixed at generic nominal values.
In panel (a), we consider a static hypergraph, i.e., Ω = 0, while in panel (b) a time-varying hypergraph, i.e., Ω = 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 µ (2)  and µ (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 ) where the system parameters are kept fixed at a 1 = 10, a 2 = 8 3 , a 3 = 28 for which individual nodes exhibits chaotic trajectory.The pairwise and higher-order structures are related to each other by L (2) = α 2 L (1) .We assume the eigenvalues of the Laplacian L (1) to be constant and the matrix b to be given by b Let us thus select as reference solution 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, ǫ 1 and ǫ 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., Ω = 0, while the panel (b) to a time-varying one, i.e., Ω = 3, one can thus appreciate that the latter setting allows a negative MSF for a larger range of parameters ǫ 1 and ǫ 2 and hence we can conclude that time-varying hypergraph enhance synchronization also in the case of chaotic oscillators.

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 2).We can observe that in both panels the critical value of coupling strengths ǫj (Ω) to achieve synchronization is smaller for Ω > 0 than for Ω = 0.This implies that synchronization can occur more easily on a time-varying higher-order structure than on a static one.
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.Figure 8 represent the result for the non-invasive coupling assumption.Here, we consider the non-invasive function so that ϕ ′ (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 µ (2) and µ (3) while keeping the remaining parameters in Eq. (A12) fixed at generic nominal values.In panel (a), we consider a static hypergraph, i.e., Ω = 0, while in the (b) panel, a time-varying hypergraph, i.e., Ω = 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 µ (2) and µ (3) and thus the SL system can achieve synchronization more easily.
generalizes the link strength.Let us observe that because of the invariance of A (d) under index permutation, we can conclude that k
FIG.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 µ(2) and µ(3) for two different choices of Ω, Ω = 0 (panel (a)) and Ω = 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 σ = 1.0 + 4.3i, β = 1.0 + 1.1i.

FIG. 6 .
FIG. 6.Synchronization on time-varying regular higher-order network of coupled Lorenz oscillators.We report the MSF as a function of the coupling strengths, ǫ1 and ǫ2, for two different values of Ω, Ω = 0 (panel (a)) and Ω = 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 a1 = 10, a2 =8  3 , a3 = 28, and α2 = 2.

FIG. 7 .
FIG.7.We show the MSF in the plane (Ω, ǫ1) (panel (a)) for ǫ2 = 0.01) and in the plane (Ω, ǫ2) (panel (b)) for ǫ1 = 0.2).We can observe that in both panels the critical value of coupling strengths ǫj (Ω) to achieve synchronization is smaller for Ω > 0 than for Ω = 0.This implies that synchronization can occur more easily on a time-varying higher-order structure than on a static one.

FIG. 8 .
FIG.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 µ(2) and µ(3) for two different values of Ω (a) Ω = 0, (b) Ω = 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 σ = 1.0 + 4.3i, β = 1.0 + 1.1i, ϕ ′ (0) = 1.