Global synchronization in generalized multilayer higher-order networks

Networks incorporating higher-order interactions are increasingly recognized for their ability to introduce novel dynamics into various processes, including synchronization. Previous studies on synchronization within multilayer networks have often been limited to specific models, such as the Kuramoto model, or have focused solely on higher-order interactions within individual layers. Here, we present a comprehensive framework for investigating synchronization, particularly global synchronization, in multilayer networks with higher-order interactions. Our framework considers interactions beyond pairwise connections, both within and across layers. We demonstrate the existence of a stable global synchronous state, with a condition resembling the master stability function, contingent on the choice of coupling functions. Our theoretical findings are supported by simulations using Hindmarsh-Rose neuronal and R\"{o}ssler oscillators. These simulations illustrate how synchronization is facilitated by higher-order interactions, both within and across layers, highlighting the advantages over scenarios involving interactions within single layers.


I. INTRODUCTION
In recent decades, the study of complex networks has gained considerable traction, emerging as a prominent area of research.This surge in interest can be attributed to their remarkable capacity to model interconnected dynamical systems across various fields, including physics, biology, ecology, social sciences, and engineering [1][2][3].These networks are comprised of nodes, representing entities or elements, and links, representing connections or pairwise interactions between them.Furthermore, researchers have introduced the concept of multilayer networks to extend the traditional network framework.Many real-world systems can effectively be conceptualized as multilayer networks [4,5].Examples include transportation networks [6], neuronal networks in the brain [7,8], and various types of social networks [9].A multilayer network consists of individual networks, each with its set of nodes and links (referred to as intralayer links), interconnected through interlayer links.The representation of multilayer networks hinges on a fundamental assumption: the complex connections among individuals within and across layers are comprehensively elucidated through pairwise links.
While pairwise interactions, such as interlayer and intralayer links, are foundational and have yielded valuable insights, real-world systems often involve more intricate relationships.Indeed, systems in the real world, spanning human communications in social networks to neuronal interactions in the brain, can be accurately depicted through multilayer networks where interactions frequently occur among groups of three or more individuals simultaneously [10][11][12][13].For example, in neuronal networks, neurons are interconnected through electrical and * dibakar@isical.ac.in chemical synapses, giving rise to a multilayer structure [8,14].Moreover, recent findings emphasize the existence of group interactions among individual neurons [15][16][17][18][19]. Similarly, the process of epidemic spreading among people involves groups of three or more interacting through virtual and physical layers [20,21].Traditional pairwise interactions fall short in capturing these group dynamics, necessitating the consideration of higher-order interactions [13,22].These complexities are represented through simplicial complexes (hypergraphs), consisting of simplices (hyperedges) of varying dimensions [23,24], where a simplex (hyperedge) of order d denotes a set of (d + 1) nodes.Recently, higher-order interactions have garnered significant attention from researchers due to their influence on collective phenomena [25][26][27][28][29][30] across various fields, including epidemiology [31,32], ecological systems [33], consensus dynamics [34], and pattern formation [35].
One captivating collective phenomenon that has garnered significant attention in the realm of multilayer networks is synchronization.Synchronization refers to the remarkable ability of coupled individuals to selforganize and exhibit collective harmony in their behavior [36,37].Various types of synchronization phenomena have been identified within multilayer networks, including cluster synchronization [38,39], antiphase synchronization [40], mixed synchronization [41], explosive synchronization [42], intralayer synchronization [14,43], interlayer synchronization [44], and relay interlayer synchronization [45,46].Each type of synchronization represents unique and intriguing behaviors exhibited by multilayer networks, making them an exciting and fertile area of exploration in network research.However, most prior works on synchronization in multilayer frameworks have considered interactions between individuals within and across layers to be pairwise, represented by links.Only a few have explored the impact of higher-order interactions on synchronization in multilayer networks [14,43,[47][48][49].However, these studies have assumed that group interactions are confined to individuals within the layers of a multilayer network, or that the underlying node dynamics are governed by Kuramoto oscillators.Recently, authors in [47] investigated a multilayer framework where interactions both within and across layers are considered beyond pairwise.However, in this study, the authors assumed all-to-all connection mechanisms among individuals and focused solely on the case of the Kuramoto oscillator.Thus, there remains ample opportunity to explore synchronization phenomena in generalized higher-order multilayered systems beyond the conventional Kuramoto model and all-to-all framework.
In this study, we aim to address this gap by exploring the phenomenon of synchronization, particularly global synchronization, within a more generalized multilayer higher-order network.This network encompasses interactions beyond pairwise connections, both within and across layers.Our proposed mathematical framework surpasses the constraints of all-to-all coupling configurations and specific oscillator models.Instead, we employ a collection of identical dynamical systems distributed across different layers, engaging in interactions within and across layers in groups of two or three (up to order three for simplicity).In this broad context, we demonstrate the emergence of global synchronization within the multilayer network, contingent upon the cancellation of coupling functions.We examine these coupling functions as either linear or nonlinear diffusive couplings.Subsequently, we establish the necessary condition for the stability of the synchronization solution, resulting in a system of coupled linear equations known as the master stability equation.Notably, our stability condition bears a resemblance to the well-established master stability function (MSF) form under specific circumstances, enabling the derivation of fully decoupled master stability equations.Finally, we validate our findings using random higher-order multilayer networks comprising paradigmatic coupled chaotic systems, such as the Hindmarsh-Rose neuronal models and chaotic Rössler oscillators.
The paper is structured as follows: Section II introduces our mathematical model for describing the generalized higher-order multilayer network.The linear stability analysis of the global synchronization state is presented in Section III.Section IV presents the numerical results corresponding to the two coupled systems, viz, Hindmarsh-Rose neurons and Rössler oscillators.Finally, Section V provides a summary and conclusion.

II. MATHEMATICAL MODEL OF THE GENERALIZED MULTILAYER HIGHER-ORDER NETWORK
We start by considering a M -layered higher-order multilayer network, which is schematized for two layers (M = 2) in Fig. 1.Each layer comprises N number of nodes, and the ith node in one layer is identical with the inode in all the other layers.Here, for simplicity's sake, we consider an identical number of nodes in each layer.The state of the i-th node in the l-th layer at time t is given by the d-dimensional state vector X li (t), where l = 1, 2, • • • , M , and i = 1, 2, • • • , N .Now, we consider that the nodes in each layer are connected to one another through pairwise links (1-simplices) and groups of three coupled differential equations give the equation of motion governing the dynamics of the multilayer network: Here, each node's dynamics are identical and illustrated as f : R d → R d .ǫ ll1 1 is the pairwise coupling strength between the nodes of layers l and l 1 connected through pairwise links.ǫ ll1l2 strength among the nodes of layers l, l 1 , and l 2 interconnected via 2-simplices.Certainly, when l = l 1 , the pairwise interactions occur among the nodes of the same layer (i.e., intralayer connections); conversely, for l = l 1 , they take place across two different layers (i.e., interlayer connections).Similarly, in the context of three-body interactions, when l = l 1 = l 2 , the triadic interactions happen among the nodes of a specific layer l (i.e., intralayer threebody interactions).For l = l 1 = l 2 , two of the nodes in a group of three are from one layer, and the third one is from another layer.Whereas in the case of l = l 1 = l 2 , all the three nodes participating in a three-body interaction are from three different layers.The last two instances correspond to interlayer three-body interactions.For the sake of simplicity, we consider the pairwise coupling strengths as ǫ ll1 1 = ǫ 1 if the pairwise interaction occurs within a specific layer, and ǫ ll1 1 = αǫ 1 for interlayer connections.The same convention has been taken for the three-body coupling strengths, i.e., ǫ ll1l2 2 = ǫ 2 if the three-body interactions take place within a specific layer and ǫ ll1l2 2 = αǫ 2 when the three-body interactions are interlayer.Here, α is an arbitrary constant satisfying 0 < α < 1, indicating that for both pairwise and three-body interactions, the intralayer coupling strengths are more effective than the interlayer ones.Thus, the parameter α quantifies the strength of interlayer interactions in our study.Indeed, the described convention mirrors real-world scenarios where individuals within the same group tend to form stronger bonds among themselves compared to their interactions with individuals from different groups.A ll1 is a N × N matrix that describes how the nodes of the layers l and l 1 are connected to one another through pairwise links.The elements A l1l2 ij = 1 if the nodes i, j in the layers l 1 , l 2 , respectively are connected, and otherwise A l1l2 ij = 0. Thus, when l = l 1 , it corresponds to the intralayer adjacency matrix, and for l = l 1 corresponds to the interlayer adjacency matrix.Similarly, B ll1l2 is the adjacency tensor that provides information about which nodes participate in three-body interactions.The elements B l1l2l3 ijk = 1 if there is a 2-simplex interaction between the nodes i, j, k in the layers l 1 , l 2 , and l 3 , respectively, while B l1l2l3 ijk = 0 if there are no such three-body interactions.h (1) : are the pairwise and three-body coupling functions such that h (1) (x, x) = 0 and h (2) (x, x, x) = 0, i.e., the coupling functions cancel out when the state of all connected individuals are equal.Thus, here we consider the pairwise and three-body coupling functions to be generalized diffusive, i.e., h (1)  [50,51].

III. STABILITY ANALYSIS OF THE GLOBAL SYNCHRONIZATION STATE
Here, we investigate the stability of the global synchronization state, wherein the nodes of the multilayer structure oscillate in unison.In mathematical terms, X li = X l ′ j for all l, l ′ = 1, 2, . . ., M and i, j = 1, 2, . . ., N .The selection of diffusive pairwise and triadic coupling functions indeed ensures the existence and invariance of the global synchronization state in our multilayer system.Now, before proceeding to the stability analysis, it is beneficial to simplify the model.To do so, we set X li = X i+(l−1)N where (l = 1, 2, . . ., M ), and (i = 1, 2, . . ., N ), effectively transforming the M layered multilayer network with N nodes in each layer into a network with M N nodes.Thus, the model Eq. ( 1) eventually becomes, where A (1) ∈ R MN ×MN is the block matrix given by, A l1l2 ∈ R N ×N such that the elements for l 1 , l 2 = 1, 2, . . ., M and i, j = 1, 2, . . ., N are defined as Here i ↔ 1 j means an intralayer link (i.e., l 1 = l 2 ) exists between the nodes i and j, and i 1 j means no such link exists.i ↔ 2 j means there exists an interlayer link (i.e., l 1 = l 2 ) between the nodes i and j, and if there is no such link, then i 2 j.Similarly, the tensor The symbols △ 1 and △ 2 indicate an intralayer triangle (l 1 = l 2 ) and an interlayer triangle (l 1 = l 2 ), respectively.
If there are such triangles between the nodes i, j, k, then i, j, k ∈ △ 1 or, i, j, k ∈ △ 2 , and if there are no such triangles, then i, j, k / ∈ △ 1 and i, j, k / ∈ △ 2 .Now, suppose X s is the global synchronization state of the multilayer network.To analyze the stability of the given state, we consider small deviations around the synchronized state, denoted as X i = X s + δX i .Subsequently, we linearize Eq. ( 2) with respect to these perturbations, which yields the variational equations as, where Jf (X s ) and JH(X s ) are the Jacobian matrices of the functions f and H evaluated at the synchronous solution.J 1 G(X s , X s ) and J 2 G(X s , X s ) are the Jacobians with respect to the first and second variables of the coupling function G at the synchronization state.To proceed further, we use the zero-row sum Laplacian matrices corresponding to the pairwise and higher-order interactions [50].If L 1 and L 2 are the Laplacians corresponding to the pairwise and three-body interactions, then we have Incorporating the definitions of these Laplacian matrices, the variational equation ( 3) can be rewritten as follows: Now considering the state variable δX = (δX T 1 δX T 2 . . .δX T MN ) T and using matrix Kronecker product ⊗, the variational equation ( 4) can be written in vectorial form as where () T represents vector transpose and I is the identity matrix.
The coupled linearized variational equation ( 5) compromises two components: one governing motion along the synchronization manifold, referred to as parallel modes, and the other handling motion across the manifold, known as transverse modes.The stability of the synchronization state is contingent upon the convergence of all transverse modes to zero over time.To perform linear stability analysis, we thus need to decompose the variational equation into parallel and transverse modes and derive the condition for the latter's extinction.Now both the Laplacian matrices L (1) and L (2) are real symmetric matrices with zero row-sum.As a result, these matrices are diagonalizable by their basis of eigenvectors, possessing nonnegative real eigenvalues.Their common smallest eigenvalue is denoted as λ 1 = 0, and the corresponding eigenvectors form an orthonormal basis.
Therefore, to distinguish transverse and parallel modes in Eq. ( 5), we project the perturbed variable δX onto the basis of eigenvectors associated with the pairwise Laplacian L (1) , where T .The choice of the eigenvector basis is arbitrary.An alternative basis can be selected, and other sets of eigenvectors can be transformed into the chosen basis through a unitary matrix transformation.Thus, expressed in the introduced variables, the variational Eq. ( 5) undergoes the following transformation: where . Now, since the Laplacian L (2) is a zero-row sum matrix, L (2) must have null first row and column.The transformed variational Eq. ( 6) then can be partitioned into two components as follows: and Here, η 1 represents the parallel mode, while the transverse modes are denoted by η j for i = 2, 3, ..., M N .Consequently, the stability analysis for the synchronized state is simplified to solving the coupled linear differential equation ( 8) associated with the transverse modes to determine the maximum Lyapunov exponent (MLE).This equation is referred to as the master stability equation.The stability criterion for the synchronization state mandates that the MLE transverse to the synchronous solution must be negative for stability to be established.The coupled master stability equation ( 8) cannot be further simplified due to the following reasons: (i) the Laplacian matrices L (1) and L (2) generally do not commute, preventing simultaneous diagonalization with respect to a single eigenbasis, (ii) the Jacobians JH(X s ) and J 1 G(X s , X s ) + J 2 G(X s , X s ) are unrelated, making it impossible to combine the Laplacian matrices to construct an effective Laplacian matrix with orthonormal eigenbasis that could simplify the equation.
To simplify the stability analysis further, we proceed with the second scenario where the Jacobians JH(X s ) and J 1 G(X s , X s ) + J 2 G(X s , X s ) are related with each other.This condition can be achieved by selecting pairwise and three-body diffusive coupling functions in such a way that they satisfy an additional requirement, namely, G(X, X) = H(X), which is commonly referred to as the natural coupling condition in the literature [35,50,51].Now, if the coupling functions satisfy the natural coupling conditions, then we have (9) Thus, in this scenario, we can start our analysis from the variational equation ( 5).Substituting the above condition (9) in Eq. ( 5) yields, (10) This allows us to introduce an effective matrix L = (ǫ 1 L (1) + ǫ 2 L (2) ), which is again a real symmetric zerorow sum matrix, possessing the smallest eigenvalue γ 1 = 0. Furthermore, L is diagonalizable with respect to its orthonormal eigenbasis.Imposing the eigenbasis for L allows us to take an additional step, projecting the coupled linear system onto fully decoupled linear systems with dimensions equal to the dimension of individual node dynamics.Each of these linear systems depends on a single eigenvalue, γ i , of L, more precisely: where is the projection of δX j on the eigenvector v (i) of the effective Laplacian matrix L. This obtained master stability equation ( 11) is fully decoupled and thus exhibits a resemblance to the classical MSF approach [52], except for the fact that here the eigenvalues of the Laplacian depend on the coupling strengths.Solving Eq. ( 11) for maximum transverse Lyapunov exponent (which corresponds to η i , i = 2, 3, . . ., M N ) provides the necessary condition for the global synchronous solution to be stable.The negative values of the maximum transverse Lyapunov exponent indicate the emergence of the stable synchronous solution.

IV. RESULTS
Here, we present a series of numerical results that serve to validate our theoretical findings.We consider two different chaotic dynamical systems, namely the Hindmarsh-Rose (HR) model [53] and the Rössler model [54] as the individual dynamics of the nodes, coupled through pairwise and three-body interactions within and across the layers of a multilayer network.The HR-neurons are coupled by a linearly diffusive coupling scheme and the Rössler oscillators are coupled by a non-linear diffusive coupling scheme.For the sake of simplicity, we proceed with only a two-layered (M = 2) system, as depicted in Fig. 1.Global synchrony occurs when all the unitary components (nodes) in the multilayer network display identical oscillations over time, i.e., X ij (t) = X s (t) (i = 1, 2; j = 1, 2, . . ., N ).To quantify the global synchronization state, we introduce the instantaneous synchronization error X ik − X jl which is zero when the multilayer network exhibits a globally synchronized state and nonzero finite for asynchronous dynamics.In the following, we will typically consider the time average of the synchronization error better to estimate the transition between the synchronous and asynchronous states.The multilayer network, given by Eq. ( 1), is therefore evaluated using the fourth-order Runge-Kutta method for a period of 3 × 10 5 time steps with integration step size δt = 0.01, and the last 10 5 time units are taken for calculating the average synchronization error.We consider identical network topologies within each layer to delve into the influence of higher-order interactions, particularly three-way interactions and the parameter α.The nodes within the layers are connected pairwisely to each other randomly with a probability P 1 following the algorithm of Erdős Rényi random networks [55].Similarly, the pairwise connections spanning across the layers are also established randomly, with a probability denoted as P 2 .Three-body connections (2-simplices) are formed both within and across the layers by promoting all the triangles formed during the pairwise network construction process.Unless specified otherwise, the number of nodes in each layer is set to N = 100, with connection probabilities within the layers and across the layers being P 1 = 0.1 and P 2 = 0.05, respectively.In the case of the neuronal multilayer network, the pairwise and three-body coupling functions are considered as H(X) = (m+n)ΓX, and G(X i , X j ) = mΓX i + nΓX j , where Γ is the inner coupling matrix, characterizing through which variable the nodes interact with each other.Clearly, our chosen coupling functions satisfy the natural coupling condition.Without loss of generality, we consider m = n = 1 2 and the inner coupling matrix is structured in such a way that all elements except Γ 11 are zero, with Γ 11 being equal to 1. Thus, the nodes within and across the layers interact with one another through the first variable of the unitary node dynamics.The chosen linear diffusive cou-pling function resembles the electrical synaptic coupling for neuronal connections [56,57].It is important to note that a linear higher-order interaction can be combined as pairwise interactions, eventually providing weighted pairwise interactions.Thus, choosing a linear higher-order interaction may not provide good insight into the effect of higher-order interactions.Here, we choose a linear higher-order coupling for neuronal networks mainly for two reasons.First, the linear diffusive coupling mimics the electrical synaptic coupling in neurons.Second, our developed analytical theory for dimension reduction of stability problems is limited to the case of diffusive coupling (linear or nonlinear).One can choose a nonlinear coupling such as chemical synaptic or noninvasive couplings.However, providing analytical support in those cases needs further restrictions such as regular topologies [14], which is beyond the scope of the present study.However, to show that both the linear and nonlinear couplings result in almost similar effects on the global synchronization, we illustrate numerical results with nonlinear noninvasive coupling for HR neuronal network in Appendix A.
On the other hand, in the case of coupled Rössler oscillators, we consider the pairwise and non-pairwise coupling functions to be nonlinear diffusive, defined as h (1) (X i , X j ) = (0, y 3 j − y 3 i , 0) T , and h (2) (X i , X j , X k ) = (0, y 2 j y k − y 3 i , 0) T , respectively.Note that this nonlinear coupling scheme also fulfills the aforementioned natural coupling condition.
In the subsequent analysis, we explore the emergence of global synchronization in the considered random multilayer structure.Additionally, we delve into the impact of various coupling parameters, including α, ǫ 1 , and ǫ 2 .Results pertaining to an alternative connection topology within the layers, specifically small-world connectivity [58], are discussed in Appendix B. Apart from that, how the interlayer and intralayer synchronization emerge due to the combined effect of higher-order and multilayer structure has been discussed in Appendix C.

A. Hindmarsh-Rose model with linear diffusive coupling scheme
The Hindmarsh-Rose (HR) neuron model is a mathematical model used to describe the behavior of a simplified neuron.This model is known for its ability to exhibit complex neuronal dynamics, including periodic spiking, bursting, chaos, and other nonlinear behaviors.It is often used in computational neuroscience and mathematical biology to study neural activity.The three coupled differential equations that describe the evolution of the three variables (x, y, z describe the membrane potential, slow current for recovery variable, and adapting current,  We consider that the dynamics of each of the nodes in the multilayer network (1) is governed by the HR neuronal model where neurons are coupled through the membrane potential by gap junctions, allowing for rapid and direct communication between neurons by facilitating the exchange of ions and small molecules.To study the emergence of global synchronization phenomena in the multilayer framework (1), we start by evaluating the synchronization error E by varying the strength of synaptic coupling of the pairwise ǫ 1 between the neurons for different synaptic coupling strength of higher-order ǫ 2 and a fixed value of α = 0.3.The corresponding results are depicted in Fig. 2(a).When there are no interactions involving three neurons (i.e., ǫ 2 = 0.0), global synchronization is achieved at a critical value of pairwise synaptic coupling strength, ǫ 1 = 0.373 (magenta curve in Fig. 2(a)).We then introduce the three-neuronal interactions and observe the effect of higher-order synaptic interac-tions on achieving global synchronization by increasing the three-body coupling strength, ǫ 2 .As ǫ 2 increases (e.g., ǫ 2 = 0.05), the synchronization state is achieved at a relatively smaller coupling strength ǫ 1 = 0.31 (shown in a black curve).As ǫ 2 is tuned up further, the system achieves synchronization much earlier, eventually leading to a significant enhancement in synchrony.For ǫ 2 = 0.1 and 0.15, the global synchronization state emerges at ǫ 1 = 0.26 and 0.233, depicted by red and blue curves, respectively.
To confirm the validity of these findings, we assess the stability of the synchronous solution using the master stability function approach and analyze the maximum Lyapunov exponent (Λ max ) of the variational Eq. ( 11), transverse to the synchronous manifold.In Figure 2(b), we present a plot of Λ max as a function of pairwise synaptic coupling strength (ǫ 1 ), with ǫ 2 values as mentioned earlier.The curves of Λ max attain negative values at the same critical ǫ 1 values, where the synchronization error (E) becomes zero (indicated by dashed vertical lines).Thus, the analysis using the master stability function approach confirms our earlier observation that introducing three-neuronal interactions in the multilayer network advances the emergence of global synchrony.To study the combined effect of pairwise and higherorder synaptic coupling on the emergence of the global synchronization state in the multilayer network, we compute the maximum Lyapunov exponent Λ max (obtained from Eq. ( 11)) as a function of ǫ 1 and ǫ 2 for four different values of the parameter α.The corresponding results are represented in Fig. 3, where the variation of Λ max is delineated through the color bars.Our findings reveal that within a certain range of ǫ 1 , the higherorder coupling strength ǫ 2 does not impact the synchronization state for all values of α ∈ (0, 1).However, a significant influence of ǫ 2 becomes evident beyond this range.Specifically, the neurons in the network do not synchronize for any value of ǫ 2 when the pairwise synaptic coupling strength falls within a certain range, which depends on the specific value of α.As this range is exceeded, the critical value of ǫ 2 , or the necessary higherorder synaptic coupling strength for synchrony, gradually decreases to zero as ǫ 1 increases.We also observed a similar impact of both synaptic coupling strengths (ǫ 1 , ǫ 2 ) on synchrony by varying the value of α.The main difference observed for the four different values of α lies in the emerged region of global synchrony (cyan-colored region) or asynchrony (magenta-colored region).We consider α = 0.2 in Fig. 3(a) and observe that after a certain range (0 < ǫ 1 < 0.25) of ǫ 1 at the beginning, the required value of ǫ 2 to get the global synchrony is found to be decreasing as ǫ 1 increases up to ǫ 1 = 0.4.In 3(b), for α = 0.4, within the range 0 < ǫ 1 < 0.25, there is no effect of the group interaction between the neurons on the emergence of synchrony.Beyond ǫ 1 = 0.25, with increasing ǫ 2 , the global synchrony emerges at a relatively lower critical pairwise coupling.Furthermore, it is observable that beyond ǫ 1 = 0.3 the critical value of ǫ 2 to achieve the global synchrony is decreased to zero, i.e., beyond ǫ 1 = 0.3, the global synchronization can emerge even if there are no three-body interactions.In 3(c), and 3(d), we take α = 0.5, 0.9, respectively, and the critical values of ǫ 1 up to which the higher-order synaptic coupling strength, ǫ 2 remain inactive, are found to be at ǫ 1 = 0.19, 0.18, respectively.Also, the critical values on the ǫ 1 -axis, where the multilayer network is in the global synchronization state without group interaction, are found at ǫ 1 = 0.26, 0.22, respectively.Notably, the synchrony region increases while the desynchrony region decreases as α ranges from 0 to 1. Thus, as the strength of interlayer coupling increases, the multilayer network can achieve global synchronization more easily.Next, we study the combined effect of α and ǫ 1 on the global synchronization for different values of ǫ 2 .The variation of Λ max of the linearized Eq. ( 11) in the (α, ǫ 1 )parameter space are shown for four values of ǫ 2 in Fig. 4. In Fig. 4(a), we consider the case when the strength of the higher-order synaptic coupling is zero (ǫ 2 = 0), which means when there is no effect of higher-order neuronal interaction in the multilayer model.It is observable that as the value of the parameter α (or, ǫ 1 ) increases, the critical value of the pairwise coupling strength ǫ 1 ( α) to achieve global synchronization decreases.Also when three-neuron interactions are introduced, i.e., for nonzero values of ǫ 2 in the remaining subfigures of Fig. 4, one can observe a similar kind of effect of the parameters α, ǫ 1 on the global synchronization state and notably an enhancement of the stable region of global synchronization is noticeable at higher values of ǫ 2 .In Fig. 4(b), we set ǫ 2 = 0.1.Here, for a fixed value of the parameter α (or, ǫ 1 ), the critical value of ǫ 1 (or, α) decreases with increasing higher-order coupling strength.In Fig. 4(c) and Fig. 4(d), we increase higher-order coupling strength to ǫ 2 = 0.2, and ǫ 2 = 0.4, respectively, and a qualitatively similar behavior is observed.This observation strengthens our statement that introducing three-neuron interactions amplifies the likelihood of synchronization among the neurons within a multilayer neuronal network.
The combined effects of the parameters α and ǫ 2 for a fixed value of the pairwise synaptic coupling strength (ǫ 1 = 0.2) are shown in Fig. 5.One can observe that global synchronization is forbidden for 0 < α < 0.06, and beyond α = 0.06, the required strength of the higherorder synaptic coupling for global synchronization decreases with the increase in α.When α is tuned up further and crosses α = 0.6, the multilayer network achieves global synchrony even without higher-order interactions (i.e., ǫ 2 = 0).FIG. 5.The variation of the maximum Lyapunov exponent, Λmax(ǫ1, ǫ2, α) of the linearized Eq. ( 11) in (α, ǫ2) parameter space for a fixed ǫ1 = 0.2, with N = 100 HR-neurons in each layer of the multilayer network (1).

B. Rössler model with non-linear diffusive coupling scheme
This section aims to show that the previously presented results hold true beyond the example of the dynamical system shown above, i.e., the HR neuronal model.For the sake of definiteness, we thus used the paradigmatic chaotic Rössler system as the underlying dynamics governing individual nodes within the multilayer network.The dynamics of the Rössler model are governed by the following equations, dx dt = −y − z, where for a = 0.2, b = 0.2 and c = 5.7, the system exhibits chaotic behavior.In our previous analysis using the HR model, we assumed the underlying pairwise and higher-order interactions to be linear diffusion as it mimics the synaptic interactions among the neurons.However, in the context of the Rössler system, we explore nonlinear pairwise and higher-order coupling schemes to delve into the impact of nonlinear interactions on the emergence of global synchrony within multilayer networks.Specifically, the coupling functions h (1) (X i , X j ) and h (2) (X i , X j , X k ) are assumed to be of the functional forms (0, y 3 j − y 3 i , 0) T and (0, y 2 j y k − y 3 i , 0) T , respectively so that it satisfies the natural coupling condition.In a manner akin to the previous exploration of the HR multilayer network, the investigation initiates by assessing global synchronization error (E) while varying the coupling strength of pairwise interactions (ǫ 1 ) for various levels of coupling strength of higher-order interactions (ǫ 2 ), with a constant α value of 0.3.The obtained results are illustrated in Fig. 6(a).When higher-order interactions are absent (i.e., ǫ 2 = 0.0), global synchronization is attained at a critical value of pairwise coupling strength, ǫ 1 = 0.00046 (red curve in Fig. 6(a)).Subsequently, the introduction of higher-order interactions reveals the impact of higher-order interactions on achieving global synchronization by increasing the higher-order coupling strength (ǫ 2 ).As ǫ 2 tuned up from zero, the synchronization state is achieved at relatively lower values of ǫ 1 .For ǫ 2 = 0.0001, synchronization emerges at ǫ 1 = 0.0004, which is depicted by the blue curve.The magenta curve corresponds to ǫ 2 = 0.0005 where global synchronization occurs at ǫ 1 = 0.00032, and the black curve represents ǫ 2 = 0.001 with the critical ǫ 1 value found to be ǫ 1 = 0.0002.To validate these findings, the stability of the synchronous solution is assessed using the master stability function approach, and the maximum Lyapunov exponent (Λ max ) of the variational Eq. ( 11) transverse to the synchronous manifold is analyzed.The curves of Λ max exhibit negative values at the same critical ǫ 1 values, coinciding with zero synchronization error (E) (marked by dotted vertical lines).The validation of the master stability function approach confirms the earlier observation that introducing higher-order interactions in the multilayer network advances the emergence of global synchrony, underscoring the significant impact of changing higher-order coupling strengths on the global synchronization of the multilayer network.Then, in the parameter space (ǫ 1 , ǫ 2 ), we investigate the influence of both pairwise and higher-order coupling strengths on the global synchronization state for four different values of α.Using color bars in Fig. 7, we depict variations in the maximum Lyapunov exponents Λ max , which is obtained from the linearized Eq. (11).Our findings reveal that within a certain range of ǫ 1 , the higherorder coupling strength ǫ 2 does not affect the synchronization state for all α values within the interval (0, 1) as in the previous case of HR multilayer network.However, a noticeable influence of ǫ 2 emerges beyond this range.Specifically, nodes (or, Rössler oscillators) in the network fail to synchronize for any ǫ 2 value when ǫ 1 falls within a specific range, dependent on the precise α value.As this range is surpassed, the critical ǫ 2 value, or the requisite higher-order coupling strength for synchrony, gradu-ally diminishes to zero as ǫ 1 increases.We also observe a similar impact of both coupling strengths (ǫ 1 , ǫ 2 ) on synchrony by varying α.It is observable that the region of synchronization expands while the desynchrony region contracts as α ranges from 0 to 1, revealing that stronger interlayer connections play an important role in the emergence of global synchronization.Thereafter, we investigate the combined impact of α and ǫ 1 on global synchronization for various values of ǫ 2 .The variation of Λ max in the (α, ǫ 1 space is depicted for four distinct ǫ 2 values in Fig. 8.In Fig. 8(a), we examine the scenario where the strength of higherorder coupling is zero (ǫ 2 = 0), indicating no influence of higher-order interaction in the multilayer model.As α (or ǫ 1 ) increases, we observe a decrease in the critical value of pairwise coupling strength ǫ 1 (or α) required for global synchronization.When higher-order interaction is introduced, i.e., for non-zero values of higher-order coupling strength in the remaining subfigures of Fig. 8, a similar trend is observed for the parameters α and ǫ 1 regarding the global synchronization state.Notably, we find an expansion of the stable region of global synchronization at higher values of ǫ 2 .This observation once again reinforces our assertion that introducing higherorder interactions enhances the likelihood of synchronization among nodes within multilayer networks.
To understand the effect of pairwise and higher-order interactions better, in Fig. 9, we plot the synchronization threshold values in (α, ǫ 1 ) parameter plane by considering three different choices of higher-order coupling strength ǫ 2 , namely ǫ 2 < ǫ 1 , ǫ 2 = ǫ 1 and ǫ 2 > ǫ 1 .The threshold couplings are the values of (α, ǫ 1 ) for which the maximum Lyapunov exponent (Λ max ) becomes zero while crossing from positive to negative regime.Thus, the region to the left (right) of the threshold curve indicates the region of desynchrony (synchrony).When ǫ 2 < ǫ 1 , we consider the relation ǫ 1 = 3.0ǫ 2 , and plotted the threshold points in (α, ǫ 1 )-plane by magenta circular points.The red square shaped points are the threshold points when ǫ 1 = ǫ 2 , and the blue diamond shaped points are the threshold points when ǫ 1 < ǫ 2 satisfying ǫ 1 = 0.5ǫ 2 .It is observed that the synchronous region in the (α, ǫ 1 )-plane is maximum when ǫ 2 > ǫ 1 , i.e., when the higher-order interactions are more effective than the pairwise ones.On the other hand, when the pairwise interactions are more effective, the obtained global synchronization region is the smallest compared to the other two cases.So far, we have discussed the results associated with a fixed interlayer connection probability, P 2 , meaning that the number of links and triangles joining the two layers remains constant.However, P 2 is a crucial parameter as it controls the number of links and triangles between the layers.When P 2 = 0, the layers are disconnected from one another, while for P 2 = 1, each node in a layer is connected to all the nodes in the other layer.Therefore, with increasing P 2 , the number of pairwise and higherorder interactions increases between the layers.Thus, to examine the effect of P 2 on the emergence of global synchronization, we calculate the critical value of coupling strengths (ǫ 1 , ǫ 2 ) for different values of P 2 , specifically, P 2 = 0.03, 0.05, 0.08, and 0.2.Here, the probability for the random intralayer connection is fixed at P 1 = 0.1, i.e., the intralayer topology remains invariant.The corresponding critical curves are portrayed in Fig. 10 for fixed interlayer coupling strength α = 0.2.The regions to the left of these critical curves represent the domain of desynchronization, while the regions to the right correspond to the domain of global synchronization.It is evident that with increasing P 2 , the region of global synchrony becomes wider.This implies that a higher number of pairwise and higher-order connections between the layers makes it easier to achieve global synchronization.Thus, the global synchronous region can be expanded by introducing more links and triangles in the interlayer connection topology.Furthermore, one can observe an interesting result for P 2 = 0.2 > P 1 , i.e., when the number of pairwise and higher-order interactions between the layers is higher than within the layers.In this scenario, the global synchronization emerges even if ǫ 1 = 0.0 (black curve).This implies that the presence of sufficient triangles in the multilayer network alone can promote global synchronization, even when there are no pairwise interactions.The threshold coupling values in the (ǫ1, ǫ2)parameter space for four distinct random interlayer topologies of the multilayer network.Each layer of the multilayer network consists of N = 100 Rössler oscillators, with a random topology for both intralayer and interlayer connections.The probability for the random intralayer connections is fixed at P1 = 0.1.The points on the red (triangle) curve represent the threshold points when the probability for the random interlayer connections is P2 = 0.03, the points on the blue (diamond) curve are for P2 = 0.05, the magenta (square) one is for P2 = 0.08, and the threshold points on the black (star) curve are for P2 = 0.2.The interlayer coupling strength is set fixed at α = 0.2.

V. DISCUSSION
In summary, we introduce a generalized mathematical model aimed at capturing intricate higher-order interactions within and across layers of a multilayer framework.Our investigation focuses on global synchronization.Within the constraints of the invariance condition for the synchronous solution, we derive the necessary criteria for achieving a stable synchronization state, thereby extending the well-established master stability function approach to multilayer structures with higher-order interactions.The complexity introduced by multiple layers and higher-order interactions in our examined system is exemplified by the manifestation of the master stability formalism as a set of coupled linear differential equations rather than a singular parametric variational equation.Despite this intricacy, we demonstrate that, in a specific scenario, our formalism simplifies to yield a set of uncoupled parametric variational equations, each possessing dimensions equivalent to those of a single dynamical unit.In addition to the theoretical foundations, we incorporate a collection of numerical results that affirm the validity and applicability of our methodology.Our findings demonstrate that introducing higher-order interactions within and across layers of a multilayer network expands the parameter space wherein global synchronization can be achieved.Therefore, our in-depth analytical exploration has provided valuable insights into the impact of higher-order interactions on the emergence of global synchronization within a generalized multilayer structure.However, it is important to acknowledge the possibility of numerous unexplored avenues for further investigation.Within this context, a natural extension of our framework involves examining synchronization phenomena in generic multilayer structures featuring distinct intralayer and interlayer coupling schemes.Additionally, delving into the influence of higher-order interactions in generalized multilayer networks on the emergence of other synchronization phenomena, such as intralayer and interlayer synchronization in detail, presents itself as a potentially intriguing research direction.
Appendix A: Hindmarsh-Rose model with nonlinear coupling scheme In our proposed multilayer HR-neuronal model given in Section IV A, both pairwise and non-pairwise coupling functions were initially considered linear due to our analytical reliance on diffusive couplings and resemblance with the gap junction (electrical synaptic) couplings.However, in this section, we adopt a nonlinear, non-diffusive coupling scheme previously used by Gambuzza et al. [50].The coupling functions are defined as h (1) , 0, 0 .This coupling is noninvasive, so it guarantees the emergence of global synchronization without any further assumptions.We, therefore, investigate the impact of non-pairwise coupling strength on the emergence of global synchronization error in the multilayer network with this coupling scheme.The result corresponding to this scenario is presented in Fig. 11, where we plot global synchronization error (E) with respect to the parameter α by varying the higher-order coupling strength ǫ 2 .Fixing the pairwise coupling strength at ǫ 1 = 0.1, we consider four different values of the nonpairwise coupling strengths, namely, ǫ 2 = 0.1, 0.25, 0.5, and 0.8.One can observe that the global synchronous state emerges at relatively smaller values of the interlayer strength α as the non-pairwise coupling strength ǫ 2 is increased.Thus, analogous to the linear higher-order coupling, the nonlinear higher-order couplings within and across the layers make achieving the global synchronization state easier.

Appendix B: Small-world connectivity
Here, we broaden the scope of our investigation by examining results under an alternative connectivity topology among nodes within the layers of the multilayer structure.Specifically, we depart from random connections among nodes and, instead, adopt the Watts-Strogatz small-world algorithm [58] to establish connections to the nodes within each layer of the multilayer framework.Nevertheless, the connections across the layers remain randomized, consistent with our previous approach.The aim is to demonstrate that the findings we have previously obtained are not confined solely to the random network structures among nodes within the layer.Small-world connectivity within each layer is formed so that each node can interact with its 10 (k = 10) nearest neighbors.Additionally, a probability p sw = 0.1, allows nodes to connect with a few more distantly located within the layers.We assume the HR model governs the dynamics of each node.Additionally, we consider the same pairwise and higherorder coupling schemes as earlier.Considering the smallworld intralayer topology in each layer and taking the HR model as the dynamics of each neuron, the maximum Lyapunov exponent of the Eq. ( 11) has been depicted in Fig. 12.We take α = 0.1, 0.3, 0.6, and 0.9 in Fig. Combining these four figures, we can conclude that the qualitative combined effect of the two synaptic coupling strengths ǫ 1 and ǫ 2 is the same as the previous (Fig. 3), i.e., the region of synchronization increases with increasing higher-order coupling.However, a significant difference can be observed.Here, the enhancement rate is fast from the previous, and in this case, the impact of the higher-order synaptic coupling strength is observable even in the absence of pairwise interactions, i.e., the sole presence of higher-order interactions is enough for the emergence of global synchronization.The effect of the parameter α on the synchronous region in the (ǫ 1 , ǫ 2 )plane is completely the same as the previous case; a rise in α value results in an enhancement of the synchronous region.Throughout the main text, all the results discussed focus on the emergence of global synchronization in the higher-order multilayer framework.Here, we delve shortly into the emergence of intralayer and interlayer synchronization in our considered multilayer network given by Eq. (1).To do so, we consider the previously taken multilayer network constituting 100 Rössler oscillators in each layer with nonlinear diffusive couplings given in the Sec.IV B. To characterize the intralayer and interlayer synchronization states, we introduce two synchronization errors as follows, X 1k − X 2j .

(C1)
Here E 1 denotes the synchronization error corresponding to layer-1.Thus, when E 1 becomes zero, the nodes within the layer-1 become synchronized.Now, it is well known that in the multilayer network, each layer synchronizes simultaneously [59].Hence, whenever layer-1 synchronizes, so does the layer-2?Eventually, the zero value of E 1 indicates the emergence of intralayer synchronization.On the other hand, E 12 represents the syn-chronization error between the layers, i.e., the interlayer synchronization error.Therefore, zero value E 12 indicates the emergence of interlayer synchronization in the multilayer network.Therefore, we investigate the variation of intralayer and interlayer synchronization errors by varying the system parameters to observe the occurrence of the synchronized states.In Fig. 13, we plot the intralayer synchronization error (E 1 ), the interlayer synchronization error (E 12 ) and global synchronization error (E) with respect to the interlayer coupling parameter α.We fix the pairwise and non-pairwise coupling strength at ǫ 1 = 0.0003 and ǫ 2 = 0.0001.We choose these two strengths very small enough so that both layers are not in the synchronous state when α = 0.0, i.e., when the layers are totally disconnected from one another.Now, after gradually increasing the value of the parameter α, it is observable that the interlayer and intralayer synchronization states, and simultaneously the global synchronization state emerge at the same value of α (α ≈ 0.5).Therefore, for fixed pairwise and higherorder coupling (drawn from the desynchrony region), inter and intralayer synchrony occur at the same threshold value of the connecting/disconnecting parameter α between the layers.
FIG. 1. Schematic diagram of the multilayer higherorder network.There are two layers, Layer-1 and Layer-2, colored in red, each composed of 10 nodes.The black solid circles denote the nodes.The solid black lines represent the intralayer pairwise interactions, and the dashed black lines represent the interlayer pairwise interactions.Triangles filled in blue represent the intralayer three-body interactions.The triangles filled in grey and orange represent interlayer three-body interactions.In the grey triangle, one node from Layer-2 is connected with two nodes from Layer-1.Whereas the orange triangle indicates the three-body interactions in which one node from Layer-1 participates with two nodes from Layer-2.
respectively) of the model are given by dx dt = y − ax 3 + bx 2 − z + I, dy dt = c − dx 2 − y, dz dt = r(s(x − x 0 ) − z)).(12) The model has eight parameters: a, b, c, d, r, s, x 0 , and I. Here, the parameter I represents an external current input to the neuron.It can be used to simulate the effect of synaptic inputs or other external influences on the neuron's behavior.Other control parameters used often in the literature are a, b, c, d, r, the first four models of the working of the fast ion channels, and the last one of the slow ion channels, respectively.The control parameter x 0 delays or advances the activation of the slow current in the modeled neuron.These eight parameter values are taken as a = 1.0, b = 3.0, c = 1.0, d = 5.0, r = 0.005, s = 4.0, x 0 = −1.6, and I = 3.25 to obtain chaotic behavior.
FIG.10.The threshold coupling values in the (ǫ1, ǫ2)parameter space for four distinct random interlayer topologies of the multilayer network.Each layer of the multilayer network consists of N = 100 Rössler oscillators, with a random topology for both intralayer and interlayer connections.The probability for the random intralayer connections is fixed at P1 = 0.1.The points on the red (triangle) curve represent the threshold points when the probability for the random interlayer connections is P2 = 0.03, the points on the blue (diamond) curve are for P2 = 0.05, the magenta (square) one is for P2 = 0.08, and the threshold points on the black (star) curve are for P2 = 0.2.The interlayer coupling strength is set fixed at α = 0.2.

5 EFIG. 11 .
FIG. 11.Variation of global synchronization error as a function of α for Hindmarsh-Rose neuronal model.Each layer of the multilayer network consists of N = 100 HR-neurons, with a random topology for both intralayer and interlayer connections and nonlinear coupling functions.The pairwise coupling strength is kept fixed at ǫ1 = 0.1.The four curves colored in magenta (circle), red (square), blue (inverted triangle), and green (diamond) represent global synchronization error for four distinct values of the non-pairwise coupling strengths given by ǫ2 = 0.1, 0.25, 0.5, and 0.8, respectively.

FIG. 13 .
FIG.13.The intralayer, interlayer, and global synchronization errors as a function of α.The multilayer is consisting of N = 100 Rössler oscillators in both the layers of the network (1) with the random intralayer and interlayer network topology.Pairwise and non-pairwise coupling strengths are fixed at ǫ1 = 0.0003 and ǫ2 = 0.0001, respectively.Three curves of colored blue (inverted triangle), red (circle), and magenta (star), correspond to the intralayer, interlayer, and global synchronization errors.