Optimizing higher-order network topology for synchronization of coupled phase oscillators

Networks in nature have complex interactions among agents. One significant phenomenon induced by interactions is synchronization of coupled agents, and the interactive network topology can be tuned to optimize synchronization. Previous studies showed that the optimized conventional network with pairwise interactions favors a homogeneous degree distribution of nodes for undirected interactions, and is always structurally asymmetric for directed interactions. However, the optimal control on synchronization for prevailing higher-order interactions is less explored. Here, by considering the higher-order interactions in a hypergraph and the Kuramoto model with 2-hyperlink interactions, we find that the network topology with optimized synchronizability may have distinct properties. For undirected interactions, optimized networks with 2-hyperlink interactions by simulated annealing tend to become homogeneous in the nodes’ generalized degree. We further rigorously demonstrate that for directed interactions, the structural symmetry can be preserved in the optimally synchronizable network with 2-hyperlink interactions. The results suggest that controlling the network topology of higher-order interactions leads to synchronization phenomena beyond pairwise interactions. Synchronization is a widespread emergent feature of complex systems. Here, the authors investigate the optimization of synchronization in phase oscillators with higher-order interactions, and find that optimized networks are more homogeneous in the nodes’ degree for undirected interactions, and for directed interactions they are generally structurally asymmetric, but can be symmetric, which differs from the pairwise case.

C omplex interactions are ubiquitous in physical 1,2 , biological 3,4 , and social systems 5,6 . The interactions form a complex network of the coupled agents. While systems with coupled agents can be modeled by networks with pairwise interactions 7,8 , where nodes of the network are connected by links, higher-order interactions are prevailing in various systems, including the network of neurons 9,10 , the contagion network 11,12 , and social networks 13 . An emerging direction in network science started to uncover the significance of higher-order interactions [14][15][16][17][18][19][20] , which induce diverse phenomena beyond pairwise interactions.
Synchronization is one of the remarkable behaviors for coupled agents [21][22][23][24] . Previous studies revealed that synchronization depends on the topology of the network with pairwise interactions. For the undirected interaction, the more synchronizable networks of identical coupled agents tend to be homogeneous in the nodes' degree distribution 25,26 . For the directed interaction, except for the fully-connected network, the optimal network in synchronizability is always structurally asymmetric 27 . Different from the pairwise interaction, higher-order interactions induce intriguing effects on synchronization [28][29][30][31] . However, except for the specific network structure, e.g., star-clique topology 28 , how the network topology of higher-order interactions affects synchronization has been seldom explored. Then, the question raises: for optimal synchronization, do the conclusions on the network topology with pairwise interaction [25][26][27] still hold when higherorder interactions are present?
In this paper, we investigate the effect of higher-order network topology on the phase synchronization of coupled oscillators on a type of hypergraph. By treating cycles on conventional networks as hyperlinks, a corresponding hypergraph can be obtained 14,32 . Therein, higher-order interactions can be formulated as higher-order hyperlinks 15,[29][30][31] , such as 1-hyperlink of two nodes (first-order interaction, pairwise interaction), 2-hyperlink of three nodes (second-order interaction), etc (Fig. 1a). The higher-order interactions from hyperlinks considered here are similar to the simplex interactions in 28 , and are different from simplicial complexes 33 or the multilayer network 34 . To analyze the network topology for optimal synchronization, we consider the Kuramoto-type coupling function for identical phase oscillators, focus on 2-hyperlink interactions, and search for the optimal network in synchronizability. Through analytical treatments and numerical estimations, we find that 2-hyperlink interactions can lead to distinct properties of the optimized networks compared with 1-hyperlink interactions.
For the undirected interaction, we rewire 2-hyperlink interactions and use simulated annealing to optimize synchronizability by minimizing the eigenratios of the generalized Laplacian matrices 28,30 . Similar to the conclusion for 1-hyperlink interactions 25,26 , the optimized networks with 2-hyperlink interactions become homogeneous in the generalized nodes' degree. For the directed interaction, we provide an example of optimally synchronizable network with directed 2-hyperlink interactions that preserves structural symmetry (in the sense that each node has the same number and same type of higher-order interactions). We rigorously demonstrate that the optimally synchronizable network with higher-order interactions can be symmetric, which is different from the result for 1-hyperlink interactions 27 . Still, the optimally synchronizable directed networks found are typically asymmetric by further numerical optimizations. Overall, the present result uncovers that the properties of synchronizable networks with pairwise interactions may or may not hold for higher-order interactions, indicating that novel behaviors can emerge in higher-order networks.

Results
Synchronizability of coupled phase oscillators with higherorder interactions. In this section, we first present the Kuramoto model with higher-order interactions. The generalized Laplacians are introduced by linearizing the Kuramoto model to study synchronization. We then focus on the case with 2-hyperlink interactions to demonstrate the effect of higher-order interactions on synchronization.
The model and generalized Laplacians for higher-order interactions. We first present the formulation of the coupled oscillator system to study synchronization with higher-order interactions and the generalized higher-order Laplacians. We consider the Kuramoto-type model with the following set of ordinary differential equations for N interacting phase oscillators (N is also the network size): where θ i ∈ [0, 2π) is the one-dimensional state variable (the phase) for the i-th oscillator, f describes the local dynamics, K d (d = 1, 2, …, D; D ≤ N − 1) are the coupling constants. Synchronization of the oscillators' phases is under consideration, which can be extended to state's synchronization by the master stability analysis 30 . For each order d, a ðdÞ il 1 l d are adjacency tensors. For example, the first-order interaction (1-hyperlink) has the conventional adjacency matrix: a il 1 ¼ 1 if the oscillators (i, l 1 ) have a pairwise interaction and 0 otherwise; the second-order interaction (2hyperlink) has a il 1 l 2 ¼ 1 if the oscillators (i, l 1 , l 2 ) have a 2-hyperlink interaction and 0 otherwise; etc. The interactions are undirected if the adjacency tensors are invariant under all permutation of indices 28 , and correspondingly they are directed if such invariance does not hold, i.e., the adjacency tensors are variant under some permutation of indices.
The functions g (d) are coupling functions for synchronization, which is assumed to be non-invasive (g (d) (θ, θ, …, θ) = 0 for ∀ d) 31 . The Kuramoto type of coupling functions 21,22,35 have g ð1Þ ðθ i ; θ l 1 Þ ¼ sinðθ l 1 À θ i Þ, g ð2Þ ðθ i ; θ l 1 ; θ l 2 Þ ¼ sinðθ l 1 þ θ l 2 À 2θ i Þ, …, g ðDÞ ðθ i ; θ l 1 ; ; θ l D Þ ¼ sinð∑ D d¼1 θ l d À Dθ i Þ. Then, the master stability equation 30,36 only depends on the adjacency tensors, i.e., generalized Laplacians, as the Jacobian terms in the master stability equation are constant. Besides, the present coupling constant K d can be denoted by the coefficients in 28 Based on the linearized equation in the master stability analysis (Methods), the generalized Laplacian matrix of the d-order interaction can be defined as 28 : where k ðdÞ ij is the generalized d-order degree between the nodes i, j, i.e., the number of d-hyperlink between i, j, and k ðdÞ i is the generalized d-order degree of node i. Note that for higher-order cases (d > 2), the Laplacian here is different from another definition of the Laplacian Eq. (38). A detailed comparison on the various definitions of the generalized Laplacians 28,30,31 is in Methods. Further, we use the adjacency tensors to represent higher-order interactions, and employ higher-order Laplacians defined from adjacency tensors. Alternatively, higher-order interactions can be introduced by the boundary matrix acting on simplicial complexes 33,37 , which leads to a different way to define higher-order Laplacians.
To study the dependence of synchronizability on the network structure with higher-order interactions, we focus on the Kuramoto type of coupling function for the coupled oscillator. Then, the master stability equation, Eq. (34), belongs to the case of Eq. (15) in 30 . As noted below Eq. (15) there, the situation is conceptually equivalent to synchronization in networks with only pairwise interactions. The summation of higher-order Laplacians now plays the same role of the conventional Laplacian from pairwise interactions. Thus, synchronizability depends on generalized Laplacian matrices 38,39 , and can be characterized by the eigenvalues of Laplacian matrices 25,27,40 .
The case with 2-hyperlink interactions. We demonstrate the higher-order effect by studying 2-hyperlink interactions, i.e, the interaction among three nodes (d = 2). We further consider the identical oscillators, i.e., each oscillator has an identical frequency 27 , where the function f(θ) = ω in Eq. (5), with ω denoting the natural frequency of the oscillators. Then,the dynamical equation becomes (Methods): Its linearized synchronization dynamics is: Synchronizability is determined by the second-order Laplacian matrix: With 2-hyperlink interactions only, all the definitions of the generalized Laplacians 28,30,31 are the same (Methods). Then, Eq. (6) can be rewritten as: In the next sections, we will separately study the optimized undirected and directed interactions for synchronizability.
Optimized synchronizable networks with undirected 2-hyperlink interactions Optimizing synchronizability by the eigenratio of the generalized Laplacian. In this subsection, we present the framework to Fig. 1 Schematic on optimizing synchronizability of the network with undirected higher-order interactions. The present higher-order interactions are from hyperlinks, similar as the simplex interactions 28 . a Various orders of interactions, such as 1-hyperlink, 2-hyperlink, etc. The 2-hyperlink is an interaction among three nodes on hypergraphs, which is the present major focus. b An illustration on optimizing synchronizability for a network with undirected 2-hyperlink interactions. An initialized network with 6 nodes (numbered from 1 to 6) and 8 2-hyperlinks is rewired to optimize its synchronizability, by minimizing the eigenratio of the generalized Laplacian. The colored triangles denote 2-hyperlink interactions, and are put on top of each other with shifting their positions for visualization. The rewiring process keeps the number of 2-hyperlinks constant, by deleting a few 2-hyperlinks and randomly adding the same number of 2-hyperlinks, as specified by the triplets and the arrow between them, e.g., (1,4,5) → (1, 3, 5). At each step, the deleted triangles have dashed lines, and the added triangles have thicker solid lines. The procedure finds the network topology with smaller eigenratios and better synchronizability. optimize synchronization of the network with undirected interactions. For the synchronized state of the system Eq. (5), we consider its bounded and connected stability region, where synchronizability of coupled oscillators can be quantified in terms of the eigenratio of Laplacian matrices Eq. (2), which was used mainly for networks with first-order interaction 25,27,40 . We optimize the networks with 2-hyperlink interactions for the linearized system Eq. (6) to determine synchronizability.
Specifically, we calculate the eigenvalues of higher-order Laplacian matrices Eq. (2), which can be arranged as 0 = λ 1 < λ 2 ≤ ⋯ ≤ λ N . Note that the eigenvalues are all real, as generalized Laplacians are symmetric. The smallest nonzero eigenvalue λ 2 is known as the spectral gap. The eigenratio quantifies synchronizability 30 . By diagonalizing the Laplacian matrix Eq. (2), we can get its eigenvalues and the eigenratio.
We make remarks about the search for the optimally synchronizable networks. First, we focus on using 2-hyperlink interactions to exemplify the effect of higher-order interactions in the optimized synchronizable network topology. The implementation can be extended to cases with higher-order interactions by a similar procedure, including the case with multiorder interactions by using the sum of higher-order Laplacians 28 . Second, since we investigate the optimal network topology for synchronizability, we have considered an identical frequency for the oscillators. Under this case, the system displays global synchronization instead of cluster synchronization 41,42 , and synchronizability is determined by the eigenratio.
Next, we provide the numerical protocol of optimizing networks with undirected 2-hyperlink interactions. We start with various randomly initialized networks, rewire second-order interactions with keeping a fixed number of 2-hyperlinks, and numerically search for optimal networks by simulated annealing 27 to minimize the eigenratio, Eq. (11) of the Laplacian matrix Eq. (7).
For the initialization, we randomly generate different networks with certain numbers of 2-hyperlinks. For N-node networks, there are C 3 N ¼ NðN À 1ÞðN À 2Þ=6 combinations of three nodes, i.e., the number of 2-hyperlinks. To demonstrate the optimization procedure, we initialize the network by first adding the 2-hyperlink interactions for the nodes i, i + 1, i + 2(i = 1, …, N). It ensures that each node has at least one 2-hyperlink interaction, such that the network does not have isolated nodes. Then, we randomly add N 2-hyperlink interactions to the network, such that each realization of optimization starts with these N 2-hyperlinks generated differently. Note that the 5-node case by such an initialization is a fully-connected network and is already optimal in synchronization. Therefore, when rewiring undirected interactions, we have chosen the minimal network size to be 6.
When rewiring the network, we delete randomly a fixed proportion of 2-hyperlink interactions, such as 20% of the existing 2-hyperlinks, and add the same number of 2-hyperlink interactions to three randomly chosen nodes which did not have an interaction before. For the 2-hyperlink interaction of the nodes l 1 , l 2 , l 3 to be deleted, a ð2Þ l i ;l j ;l k ¼ 0 with i, j, k being all the permutations of 1, 2, 3. After each rewiring step, we calculate the eigenratio of the rewired network, and a Metropolis acceptreject step is used for the rewired network with smaller eigenratio 27 . The chosen proportion of rewiring at each step would not dramatically affect the final optimization result, when the ratio of the rewiring was not too larger (e.g., ≥40%) and a sufficient number of optimization steps were conducted (Supplementary Fig. 1), because the eigenratios converge to a stable range during the optimization. However, the global minimum of eigenratios is not always guaranteed, and the numerical minimization may gradually lead to the global minimum with an expense of longer computation.
The optimization procedure is repeated until the eigenratio is smaller than a chosen target value, which is initially set as the smallest possible eigenratio 1. The target eigenratio cannot be too small, because otherwise it may not be reached due to the sparsity of 2-hyperlink interactions (2N) when the number of nodes increases. We thus increase the target eigenratio by 10% if the number of rewiring steps runs over 100 times without reaching the current target eigenratio. This procedure automatically increases the target eigenratio, to reduce the computational time of searching for too small eigenratio that may not be achieved for the sparse large networks. It still ensures to optimize synchronizability by minimizing the eigenratio.
The above completes one realization of optimization, and 1000 realizations are conducted with various configurations of initialized networks. In total, we have two hyperparameters about the iteration. The first is the maximum number of iterations to reduce the eigenratio to the target value (100 times) before increasing the target value. The second is the number of realizations (1000), i.e., the number of different initialized networks. The computational time increases with these two hyperparameters, and also increases with the number of nodes. When operating on a personal desktop computer, the computation time can be hours when the network size exceeds 100. Larger values of these two hyperparameters can be used to enable the search for smaller eigenratios, with a cost of more computational resources.
The network is rewired instead of simply deleting the 2-hyperlink interactions, because after deleting the 2-hyperlinks the eigenratio may be similar while all eigenvalues continue to become smaller. It eventually leads to an optimal network with much fewer 2-hyperlink interactions than the initial network. The optimization with only adding 2-hyperlink also gives less constrained network structures. Thus, we choose to rewire the network by adding the same number of 2-hyperlink interactions after the deletion, which preserves the total number of 2-hyperlink interactions. Different types of constraints can be employed in the optimization procedure, to search for the optimized network with desired properties.
Optimizing synchronizability of a 6-node network. In this subsection, we give an example with its second-order adjacency tensor and generalized Laplacian to exemplify the network with 2-hyperlink interactions. We further demonstrate the optimization procedure by this example.
Specifically, we consider the network with 6 nodes and 8 2hyperlinks, and randomly initialize the 2-hyperlink interactions. It has in total 8 2-hyperlinks connecting the nodes: Its generalized Laplacian matrix L (2) by Eq. (7) is: The eigenvalues of the generalized Laplacian matrix in ascending order are 0, 6.171, 8.167, 10.000, 10.549, 13.111, and the eigenratio is 2.125.
We next optimize synchronizability of this example. In each step of the optimization, one or two 2-hyperlinks may be rewired to generate a network with a smaller eigenratio. An illustration is given in Fig. 1(b). After conducting the numerical optimization with preserving 6 nodes and 8 2-hyperlinks, the resultant optimized network has the following 8 2-hyperlinks connecting the nodes: The generalized Laplacian matrix L (2) by Eq. (7) is: The eigenvalues in ascending order are 0, 8,9,9,11,11, and the eigenratio is 1.375. We note that the optimized network may not be the best in synchronization until sufficient numerical search is conducted. However, with the present numerical optimization, this network is at least near to the optimal network in synchronizability as its eigenratio is close to 1.
The optimized networks with various sizes. In this subsection, we provide the optimization result for undirected networks with various sizes. Examples of the initial and rewired network are given in Fig. 2a, b, with the number of nodes N = 7, 8,9,10. It demonstrates that the optimization procedure rewires 2-hyperlink interactions and reduces eigenratios (Fig. 2c). For illustration, we have shown networks with a small number of nodes. In Fig. 3a, the eigenratios before and after the optimization are provided, where the numbers of nodes include 6 and those from 10 to 100 with a step size 5. The optimized networks have smaller eigenratios after the optimization, showing better synchronizability. Indeed, directly simulating the Kuramoto dynamics in Eq. (5) validate that the optimized network has better synchronizability than the initial network ( Supplementary  Fig. 2).
With a fixed number of 2-hyperlinks, the specific configuration of the optimized network does not dramatically affect the final eigenratio when conducting multiple numerical replicates, such as 10 times, on the network structures, as shown by the errorbars in Fig. 3b. The number of 2-hyperlinks is a crucial factor in determining synchronizability. The numbers of initialized 2-hyperlinks are 2N, which gives a sparse network and only allows a relatively large eigenratio after the optimization. The eigenratios are smaller when the number of 2-hyperlinks increases to improve synchronizability, i.e., larger network densities lead to relatively smaller eigenratios ( Supplementary Fig. 3). Specifically, besides the density of 2-hyperlinks 2N in Fig. 3, we used a different number of triangles, e.g., 3N, 4N, 5N in Supplementary Fig. 3. The eigenratios of initialized networks become smaller with the increasing density, because more links improve the synchronizability of networks. At the same time, there is less room to rewire the network for optimizing synchronizability if the number of 2-hyperlinks becomes abundant. Under the chosen densities, our numerical algorithm can still find better synchronizability and reduce the eigenratio by rewiring the networks. In addition, the optimized networks under various densities also tend to have more homogeneous nodes' degrees.
Different ways of initialization can be implemented to investigate the dependence on the initialized networks. For example, one may generate random hypergraphs with 2N 2-hyperlinks and only use those connected hypergraphs to do the optimization. The results in Supplementary Figs. 4, 5 show a same qualitative conclusion as Figs. 2, 3. This indicates that the specific way of initialization does not affect the optimization result, once sufficient rewiring steps are conducted during the optimization.
We have further calculated the generalized degree k ð2Þ i in Eq. (8) of each node, which quantifies the number of 2-hyperlinks participated by each node. The distribution of the generalized degree for the optimized networks with various sizes is in Fig. 4a-f. When considering a identical frequency distribution of oscillators, the optimal network tends to be more homogeneous in the nodes' degree for the network with firstorder interactions 25,26 . Similarly, the optimized network with 2-hyperlink interactions also becomes more homogeneous, as the nodes' degree distribution concentrates to fewer values of degree in Fig. 4. We further calculate the homogeneity parameter [43][44][45][46] : where κ is the mean degree of the network defined as κ N À1 ∑ N i k ð2Þ i for the 2-hyperlink interactions and σ 2 is the variance of the nodes' degrees. Note that g becomes a delta distribution for fully homogeneous network, and tends to 0 for networks with more heterogeneous degree distributions. The homogeneity parameter for rewired network is small compared with the initial network (Fig. 4), showing that the rewired networks are more homogeneous. Other measures such as generalizing the clustering coefficient 47 , degree-degree correlations 43 , characteristic path length and heterogeneity measures 48 to high-order networks can help quantify and control synchronizability of the network.
Optimized synchronizable networks with directed 2-hyperlink interactions. In this section, we study directed 2-hyperlink Fig. 3 Eigenratios of networks with 2-hyperlink interactions before and after the optimization on synchronizability. a The eigenratio variation before and after the optimization, which shows that the optimizing procedure decreases eigenratios and enhances synchronizability. The network sizes include 6 and those from 10 to 100 with a step size 5, denoted by color. b The eigenratios of the optimized network with 2-hyperlink interactions for various size 6, 10, 20, 50, 80, 100. The errorbar denotes the standard deviation of 10 numerical replicates, and each replicate is the optimal from 1000 different randomly initialized networks. The generalized nodes' degree distribution shows that the more synchronizable network tends to be homogeneous. Correspondingly, the homogeneity parameter 43-46 g given by Eq. (16) is bigger for the rewired network, indicating more homogeneous nodes' degrees. The property of the more synchronizable 2-hyperlink interactions is consistent with 1-hyperlink interactions, where nodes also tend to be homogeneous in degree 26 . interactions to demonstrate the effect of directed higher-order interactions. We focus on the relation between synchronizability and structural symmetry of the network 27 . A network is regarded as structurally symmetric when each node has the same number and same type of (higher-order) interactions.
The optimized directed network with only first-order interactions has been studied 27 : it has been proved that the optimally synchronizable directed network is always structurally asymmetric (except for the fully-connected network). The proof was done by establishing a contradiction showing that the structurally symmetric network cannot be optimal in synchronizability. However, whether this conclusion holds for networks with higher-order interactions remains unknown.
For a higher-order directed network, we find that the symmetry may hold for the network with directed higherorder interactions, i.e., there are structurally symmetric networks which are optimal in synchronizability. Intuitively, higher-order interactions provide more capacity to reach optimal network design, enabling optimally synchronizably directed networks to preserve structural symmetry.
Before presenting the result, we recapitulate the definition of the directed hyperlinks. The directed 2-hyperlinks are defined such that each permutation of three nodes leads to a distinct 2hyperlink's direction: a 2-hyperlink interaction is "directed" once the tensor a ð2Þ i;j;k , with i, j, k being all the permutations of 1, 2, 3, has one of its six elements to be nonidentical. Then, each a ð2Þ i;j;k with i, j, k having a specific order can be regarded as a directed hyperedge, i.e., an ordered pair of disjoint subsets of vertices 49 . The ordering i, j, k specifies the direction, as illustrated in Fig. 5a. For example, with a ð2Þ i;j;k ¼ 1, i is a source node while j, k are its target nodes, and j is a source node while k is its target node. By using the adjacency tensors for directed interactions, the linearized synchronization dynamics is also given by Eq. (6), and synchronizability depends on eigenvalues of the generalized Laplacian matrix in Eq. (7).
An example of optimally synchronizable network with structural symmetry. First, we demonstrate that a network with N = 7 nodes and a symmetric structure can be optimal in synchronizability. As illustrated in Fig. 5(b), the network has the 3 × 2 × 7 directed 2-hyperlinks given by the following nonzero elements in the adjacency tensor: Other elements in the tensor a (2) are zero. From this adjacency tensor, its generalized Laplacian matrix by Eq. (7) is: L ð2Þ ¼ 6 À1 À1 À1 À1 À1 À1 À1 6 À1 À1 À1 À1 À1 À1 À1 6 À1 À1 À1 À1 À1 À1 À1 6 À1 À1 À1 À1 À1 À1 À1 6 À1 À1 À1 À1 À1 À1 À1 6 À1 À1 À1 À1 À1 À1 À1 6 The Laplacian has all nonzero eigenvalues identical as 7 and eigenratio 1. This network with 2-hyperlink interactions has the structural symmetry, as each node has the same number and type of 2-hyperlink interactions: each line of the 2-hyperlinks in Eq. (17) belongs to the same type of 2-hyperlink for the first node (source node). Therefore, it is a counterexample of only the structurally asymmetric network being optimal in synchronizability. Higherorder interactions make it possible to have an optimally synchronizable network with symmetry. We expect that larger networks can have more symmetric optimal structure.
The directed 2-hyperlink interactions enable structural symmetry. The unweighted higher-order network is most synchronizable when the real eigenratio (Eq. (21)) is the smallest. That is, the nonzero eigenvalues of the Laplacian matrix satisfy λ 2 = λ 3 = ⋯ = λ N when the network is most synchronizable 40 . i;j;k has its subscript indices separately as source, 1st target, 2nd target nodes as in 56 . Two directed 2-hyperlinks corresponding to a ð2Þ 1;3;5 ¼ 1, a ð2Þ 1;5;3 ¼ 1 are shown. b The 7-node network with directed 2-hyperlink interactions. The colored triangles denote various directed 2-hyperlink interactions, as given by Eq. (17). For example, two directed 2-hyperlinks a ð2Þ 1;3;5 ¼ 1, a ð2Þ 1;5;3 ¼ 1 in total have the arrows 1 → 3, 1 → 5, 3 → 5, 5 → 3 by the rule defined in panel a for setting the directed 2-hyperlink. This 7-node network supports that the structurally symmetric network can be optimal in synchronization for 2-hyperlink interactions. Differently, the optimal directed networks with pairwise interactions are always structurally asymmetric 27  For pairwise interactions, this condition implies that the eigenvalues of the optimally synchronizable network are all real, such that they can be ordered, even though the directed networks generally have complex eigenvalues. This property can be extend to the higher-order cases with using generalized Laplacian 30 . By a similar procedure in the Supplementary of 40 , we found that the identical eigenvalues are integers as follows. Specifically, we first define λ ∑ N i¼2 λ i =ðN À 1Þ, and λ ¼ λ 2 ¼ λ 3 ¼ Á Á Á ¼ λ N > 0 when the network is most synchronizable. The characteristic polynomial of the generalized Laplacian is: detðL ð2Þ À λIÞ ¼ Àλðλ À λÞ where I is the identity matrix. As L (2) has all integer entries, the coefficients of the characteristic polynomial should be all integers, and thus C λ NÀ1 is an integer. According to the definition of the generalized Laplacian L (2) in Eq. (7), trðL ð2Þ Þ ¼ l ¼ ðN À 1Þ λ, where l denotes the number of elements 1 in the adjacency tensor a (2) .
where integers s and t do not have common factors. For Ct N−1 = s N−1 , any prime factor p of C must be a factor of s N−1 and consequently a factor of s. Since there are no common factors in s and t, p N−1 needs to be a factor of k. Therefore, any factor of C has multiplicity N − 1, giving C = q N−1 with an integer q. This leads to λ NÀ1 ¼ q NÀ1 and that λ is an integer. By using the above properties of the eigenvalues to establish a contradiction, the authors 27 found that an optimally synchronizable network with first-order interaction (except for the fullyconnected network) must be structurally asymmetric. Below, we provide an attempt to establish such a contraction for the network with directed 2-hyperlink interactions, and find that the contraction can no longer be established.
On the one hand, in a symmetric network, the nodes are structurally identical, i.e., each node has the same number and type of higher-order interactions. It implies that the in-degrees and out-degrees from the 2-hyperlink interactions of all nodes must be equal. Thus, l needs to be divisible by N if the network is symmetric. As an example, for the network with 4 nodes and 2-hyperlink interactions, the condition of structural identity requires each node to have the same degree and the same number of 2-hyperlinks. Then, this 4-node case needs to have a fullyconnected network, with all the 4 2-hyperlinks and l = 12 divisible by 4.
On the other hand, when the network is optimal in synchronizability, the eigenvalues are integers and equal: λ 2 ¼ λ 3 ¼ Á Á Á ¼ λ N ¼ λ, and trðL ð2Þ Þ ¼ l ¼ ðN À 1Þ λ. These properties imply that λ ¼ l=ðN À 1Þ and that l must be divisible by N − 1 if the network is most synchronizable. Therefore, l must be divisible by (N − 1)N and then For the network with 2-hyperlink interactions, The two conditions can be satisfied simultaneously when (N − 1) N ≤ l ≤ (N − 2)(N − 1)N. It no longer constrains the network to be fully-connected as the case of the network with only first-order interactions 27 . Thus, the contradiction that the structurally symmetric network is not optimal in synchronizability cannot be reached for the network with 2-hyperlink interactions. We remark that our proof can be extended to higher-order interactions. For example, the 3-hyperlink would have l ≤ (N − 3) (N − 2)(N − 1)N in Eq. (20). Then, following the same procedure of proof, the 3-hyperlink also allows the structurally symmetric network to be optimal in synchronizability, different from the first-order interaction.
Optimizing synchronizability by the real eigenratio of the generalized Laplacian. Generalized Laplacians for higher-order interactions are typically not structurally symmetric, and lead to complex eigenvalues. The eigenvalues of the second-order Laplacian L (2) can be listed in ascending order of their real parts: λ 1 , λ 2 , λ 3 , …, λ N . In the strong coupling regime, since the stability region of the fully-synchronized state is bounded and connected, synchronizability can be quantified by an eigenratio of the real part 27 : where ℜ denotes the real part of the complex eigenvalues. The network will be more synchronizable if this eigenratio is smaller. Note that here we can extend this property to the higherorder case, because generalized Laplacians play the same role on quantifying synchronizability when Eq. (6) has identical oscillators and specific coupling functions (see the discussion under Eq. (15) of 30 ). We numerically optimize synchronizability of the directed network by minimizing Eq. (21). We next provide an example of optimized networks by the numerical optimization. The initial network has 6 nodes and the following 2 × 8 undirected 2-hyperlinks: The other elements in the adjacency tensor are zero. By deleting directed 2-hyperlinks, various configurations of the optimal network with R = 1 can be reached by the numerical optimization. For example, one optimized directed network from the numerical optimization has the nonzero elements in the adjacency tensor: with other elements of a (2) zero. Its generalized Laplacian matrix by Eq. (7) is: The eigenvalues in ascending order of real parts are 0, 2, 2, 2, 2, 2, and the eigenratio is R = 1.
Optimized directed networks are structurally asymmetric in general. Though the counterexample shows that the optimal network can be structurally symmetric when higher-order interaction presents, it does not mean that optimal directed networks with higher-order interactions generally tend to be symmetric over asymmetric. We use simulated annealing to search for the optimized synchronizable network numerically. For the directed interaction, we separately delete and rewire the directed 2-hyperlink interaction randomly to achieve better synchronizability. After removing or rewiring directed higher-order interactions to optimize synchronizability, the network typically becomes structurally asymmetric. We first present the optimization by removing directed 2-hyperlink interactions, such that the network can become directed and asymmetric. By removing directed 2-hyperlinks, such as setting a 1,2,3 , a 2,3,1 or a 3,1,2 , etc to zero, the network can have directed 2-hyperlink interactions instead of undirected interactions. We employ simulated annealing for the optimization, where the input is a three-dimensional tensors for a network with 2-hyperlink interactions. We then calculate the generalized Laplacian matrix L ð2Þ ij for a given tensor a l 1 ;l 2 ;l 3 by Eq. (7), modify the network and calculate Eq. (21) after each step of the modification. We have used the same procedure to set the target value in the previous section, where the target value is increased by 10% if the network modification runs over 100 steps. Then, 1000 realizations are conducted with different configurations of initialized networks.
We considered network sizes 6, 10, 20, 50, 80, 100. After optimizing synchronizability by modifying the network, we use two quantities to measures the asymmetry of the optimized network. First, we count the number of directed 2-hyperlink interactions of each node, which measures the structural asymmetry from each 2-hyperlink interaction for each node: l 1 ;l 2 ;l 3 À a ð2Þ where A 1 denotes the first asymmetry measure. Second, we measure the asymmetry on the directed in-and-out interaction for three nodes of each 2-hyperlink interaction: ½a ð2Þ l 2 ;l 1 ;l 3 þ a ð2Þ l 2 ;l 3 ;l 1 þ ½a ð2Þ l 3 ;l 1 ;l 2 þ a ð2Þ l 3 ;l 2 ;l 1 À 2½a ð2Þ l 1 ;l 2 ;l 3 þ a ð2Þ where A 2 denotes the second asymmetry measure. We estimate these asymmetry measures for each node, and then average them over all the nodes. The two asymmetric measures A 1 , A 2 of 30 numerical replicates are plotted as violin distributions in Fig. 6a, b. They show that the optimized networks are generally asymmetric, because most numerical replicates generate the optimized network with nonzero asymmetric measures. Note that these measures may not quantify structural asymmetry of having 2-hyperlinks symmetrically in the network. We further evaluate the network density for the directed network in Fig. 6c, by the average value of the generalized nodes' degree Eq. (8), showing the network densities for the optimized directed networks.
The phenomenon that asymmetry enhances synchronization is for the directed networks, which is not contradictory to optimal synchronization of fully-homogeneous undirected networks 14,25 .
Besides, after deleting 2-hyperlink interactions, the eigenvalues overall become smaller. For the optimized network obtained above, we have chosen the minimum number of deletions to reach the target eigenratio.
We have further rewired the directed 2-hyperlink interactions without deletion (Supplementary Fig. 6). When the 2-hyperlink interactions are rewired without being deleted, the asymmetric measures are larger as more interactions are kept. Regardless of the exact values of the asymmetric measures, the optimized directed networks in general are also structurally asymmetric.

Conclusion
In summary, our result demonstrates that higher-order interactions can lead to distinct properties of optimized synchronizability compared with the conventional network with pairwise interactions. For the undirected interaction, the more synchronizable networks tend to be homogeneous, consistent with pairwise interactions 25,26 . For the directed interaction, the optimized synchronizable network is structurally asymmetric in general but can be symmetric, beyond the first-order case 27 . The optimization on synchronization of higher-order interactions may find uses in real networks, such as controlling the asynchronous state 50 of brain networks with high-order structures 9 .
Recent studies 13,15 revealed significant roles of higher-order interactions in network science. In light of these studies, we investigate how higher-order interactions affect synchronization by optimizing network topology. Specifically, we have demonstrated the effect of higher-order interactions by focusing on the network topology of 2-hyperlinks, considering the phase synchronization of coupled oscillators, and employing a Kuramoto type of coupling function. The cases with general coupling function can be studied by the master stability analysis 30,31,36 , which can determine the stability of the synchronized state. Under that case, an interplay between the coupling function and the network topology needs to be analyzed. However, for a large class of coupled oscillator systems 27 , where the stability region is bounded and connected, synchronizability can be quantified by the eigenvalues of the generalized Laplacian.
When calculating eigenvalues of the generalized Laplacian matrices, the analytical solvable cases are restricted to special networks 28 , and numerical estimations are generally required. In the numerical implementation, the initialization on the network topology in Results ensures a sufficient number of 2-hyperlinks to be optimized. On the other hand, abundant 2-hyperlinks may cause an ineffectiveness of the optimization, because the network is already near to optimal synchronizability when the number of 2-hyperlinks is abundant. Different types of initialized networks can be used to further improve the search of the optimal network, with a cost of longer computational time. Similarly, though the chosen network sizes are sufficient to show the distinct properties of higher-order networks compared with the pairwise interactions, applying our numerical package to larger networks is useful when more computational resources are available.
When oscillatory frequencies are heterogeneous, synchronizability depends on both network structure and oscillators' frequencies to be optimized. For pairwise interactions 26 , synchronization can be enhanced by a match between the heterogeneity of frequencies and network structure. For higher-order networks, one also needs to optimize the frequencies and the alignment function simultaneously. Here, we have focused on the network topology and considered identical oscillators. Besides, to study the network topology, we have treated the network as a single cluster, rather than a network with a sub-cluster coupled to other nodes 27 .
We have used adjacency tensors to encode higher-order interactions, which can be formulated as simplicial complexes or hypergraphs 51 . We consider the case with pure 2-hyperlink interactions 28 . Instead, simplicial complexes require collections of simplices 15 . Extending the present result to simplicial complexes needs to include the various orders of simplices. To include multiorder interactions simultaneously depends on the coupling function. For general cases, there is an issue of diagonalizing the multiorder Laplacian matrices simultaneously 31 , because the multiorder Laplacian matrices cannot be directly added up due to the coupling function. For the Kuramoto type of coupling function, the Laplacian matrices can be added up 28 , and then the multiorder eigenvalues determine the stability and quantify synchronizability. The numerical implementation is extendable to higher-order interactions by using generalized Laplacian matrices in Eq. (2). Future work includes to investigate the optimal network design for the case with multiorder interactions 52,53,54 and with attractive and repulsive terms 28 .
To define the direction in higher-order networks, we have extended the definition of the directed network with pairwise interactions 27 to the 2-hyperlink case: the interactions are directed if the adjacency tensors are variant under the permutation of the indices. The direction of 2-hyperlink is assigned in a same way as that of the triangle in directed simplicial complexes 55 . Moreover, the direction of hyperedge in general hypergraphs needs to be carefully defined 56 , and one needs to control hyperedge including their directions in the hypergraph. To extend the present result of 2-hyperlink interactions to general directed hypergraphs is another interesting topic, which starts to be uncovered 57 .

Methods
The master stability equation for the system of coupled oscillators. In this subsection, we present the method of master stability function 36 , which can be used to obtain the synchronization phase diagram for general coupled-oscillator systems Eq. (1). We consider the system with general coupling functions for the first and second-order interactions, as higher-order interactions can be treated similarly 30 .
In order to study the stability of the synchronization solution, we consider small perturbations around the synchronous state, i.e., δθ i ¼ θ i À θ s i , where θ s denotes the steady state. We then perform a linear stability analysis on the vector δθ = (δθ 1 , δθ 2 , …, δθ N ). The master stability equation 36 can be obtained as Eq. (11) of 30 : where ⨂ is the matrix Kronecker product, I N is the n-dimensional identity matrix and the Jacobian terms at the synchronized states are: After diagonalizing the first-order Laplacian matrix L (1) by a linear coordinate transformation δθ → δη, the master stability equation becomes Eq. (13) of 30 : where i = 2, …, N denotes the different modes transverse to the synchronization manifold, λ i are the eigenvalues of the Laplacian L (1) , andL ð2Þ is the transformed second-order Laplacian L (2) . With the master stability equation Eq. (31), the master stability function Λ max , i.e., the maximum transverse Lyapunov exponents (the modes except for the first stable mode), can be obtained by numerically tracking the norm of: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ∑ N i¼2 ∑ m j ½η ðjÞ i 2 q , with η i ðη ð1Þ i ; η ð2Þ i ; ; η ðmÞ i Þ solved from Eq. (31) under given parameters. It is beyond solely using the eigenvalue of the Laplacian matrices for the Kuramoto-type coupling function.
When considering the Kuramoto-type coupling function and the linearized coupling terms 25 , the interaction term becomes variable-independent and eigenvalues of Laplacian matrix can quantify synchronization. Under this case, one could analyze eigenvalues of the generalized Laplacian matrices Eq. (2) to search for the optimally synchronizable network. Without the linearization, the coupling function is variable-dependent. Then, we need to solve Eq. (31) to quantify synchronizability, which prohibits to efficiently search for optimal network structure. We thus consider the Kuramoto-type coupling function 28 in the main text. As for the linearization, in the strong coupling regime, the system Eq. (5) will converge to a synchronized state with θ i ≈ θ j for any i, j. The linearized synchronization dynamics for the synchronized state is: The linearized equation reduces to the same form of Eq. (2) in 28 . We further use the rotating reference frame by θ i = θ i − ωt, where the synchronized solution is θ i (t) = 0, (i = 1, …, N). Then, the master stability equation for Eq. (5) is: a ð2Þ il 1 l 2 ðδθ l 1 þ δθ l 2 À 2δθ i Þ; ð34Þ which determines the linear stability of the synchronized state. This master stability equation is a case of Eq. (15) in 30 . It leads to Eq. (5) in the main text, where only the 2-hyperlink interaction is kept.
Generalized Laplacians for higher-order interaction. In this section, we compare the different definitions on generalized Laplacians for higher-order interactions. In short, though their definitions for higher-order interactions in 28,30,31 may differ, the first-order and second-order Laplacians are the same. As this paper focuses on the topology of 2-hyperlink interactions, the obtained results are valid for using any one of the generalized Laplacians defined in 28,30,31 . For convenience, we have used Eq. (35) as in 28 when presenting the generalized higher-order Laplacians in the main text.

•
The Laplacian in 28 is given by their Eqs. (3,4,5): a ðdÞ i;l 1 ; ;l d ; ð36Þ a ðdÞ i;j;l 2 ; ;l d ; which is used in Eq. (2) of the main text. We list it here for the comparison with other definitions of Laplacians.
• The generalized Laplacians for higher-order interaction terms defined in 30 is by their Eq. (28): where k ðdÞ ij also denotes the generalized d-order degree of the nodes i, j and k ðdÞ i denotes the generalized d-order degree of node i. Therein, their Eqs. (24,26) give: a ðdÞ i;l 1 ; ;l d ; ð39Þ a ðdÞ i;j;l 1 ; ;l dÀ1 : Here, d specifies the order, which gives a 3-dimensional tensor when d = 2. Note that k ðdÞ i , k ðdÞ ij here are the same as those in Eq. (35). However, the definitions of Laplacians are different: Eq. (35) multiplied by (d − 1)! becomes the same as Eq. (38). Still, for the first and second orders (d = 1, 2), the two definitions of Laplacians reduce to the same.
• The generalized Laplacian for higher-order interactions used in 31 are: to retain the zero row-sum property. It is consistent with that in 30 for the second-order.
With using generalized Laplacians for the first-order and second-order cases, we obtain the optimized networks with 2-hyperlink interactions.

Data availability
The authors declare that the data supporting the findings of this study are available within the paper.