Metamorphoses and explosively remote synchronization in dynamical networks

In physics, the existence of a symmetry implies the conservation of a physical quantity and a great deal can be learned about the system without analyzing the intermediary details of the system. The principle of symmetry often leads to the discoveries of unexpected and surprising phenomena. In complex networks, intricate dynamical phenomena such as cluster synchronization are the results of symmetry.^1–5 Here, we report such a phenomenon: synchronization can switch abruptly from one group of nodes to another as a bifurcation parameter passes through a critical point, leading to explosively remote synchronization, signifying a synchronization metamorphosis in nonlinear dynamical networks.

Synchronization has been an extremely active research topic in nonlinear and complex dynamical systems,^6,7 which is manifested as the emergence of coherent motion among coupled dynamical units when the interaction or coupling is sufficiently strong. Some background of synchronization in nonlinear and complex systems is as follows. In a typical study, the setting is coupled nonlinear oscillators, where the bifurcation or control parameter is the coupling strength among the oscillators. Historically, a central task is to identify the critical point at which a transition from desynchronization to synchronization occurs.^6,7 Depending on the dynamics of the oscillators and the coupling function, the system can have a sequence of transitions, giving rise to distinct synchronization regimes in the parameter space. For example, for a system of coupled identical nonlinear oscillators, complete synchronization can arise when the coupling exceeds a critical strength that can be determined by the master stability function.^27,28 Systems of nonlinearly coupled phase oscillators, e.g., those described by the classic Kuramoto model⁶ can host phase synchronization and the critical coupling strength required for the onset of this type of “weak” synchronization can be determined by the mean-field theory.^29,30 Synchronization in coupled physical oscillators was experimentally studied.^31,32 A counterintuitive phenomenon is that adding connections can hinder network synchronization among time-delayed oscillators.³³ Some previous key studies of synchronization in complex networks are as follows. It was found that small-world networks, due to their small network diameters, are more synchronizable than regular networks of comparable sizes,³⁴ but heterogeneity in the network structure presents an obstacle to synchronization.³⁵ Subsequently, it was found that heterogeneous networks with weighted links can be more synchronizable than small-world and random networks.³⁶ Synchronization in complex clustered networks³⁷ and the onset of chaotic phase synchronization in complex networks of coupled heterogeneous oscillators were also studied.³⁸ The interplay between network symmetry and synchronization was uncovered^39,40 and understood.^2,3,17 In particular, a symmetry group can be generated by the possible symmetries of the network and the orbits of the group determine the partition of the synchronous clusters. In general, the phase space of the whole networked dynamical system can be decomposed into the synchronization subspace and the transverse subspace through a transformation matrix generated by the symmetry group, which determines the stability of the cluster synchronization patterns.^2,3

Recent years have witnessed the discovery of two remarkable synchronization phenomena: explosive synchronization^8–15 and remote synchronization.^16–21 In explosive synchronization, as a system parameter changes through a critical point, synchronization emerges suddenly and discontinuously (in the sense that some measure of the synchronization abruptly increases from zero to a finite value), signifying a first-order phase transition. In remote synchronization, certain non-directly coupled dynamical units are synchronized but those that are directly coupled are not. Previously, the two phenomena were regarded as unrelated and studied separately. In this paper, we demonstrate that the two phenomena can occur simultaneously in the same network system through a synchronization metamorphosis.

A mechanism for explosive synchronization is that, at the abrupt transition point, a number of small sized synchronous clusters have already existed, and the transition is essentially a percolation process leading to a giant connected component of these clusters.²² Remote synchronization is a manifestation of the fundamental symmetry in the network. Our idea is that, in a network with certain symmetry, the nodes are organized as symmetric clusters: within each cluster, nodes are locally connected but there is no coupling among nodes in different clusters. Each node in a cluster has symmetric counterparts in other clusters. The sets of symmetric nodes in all the clusters constitute various layers, where there is no direct coupling in any of the layers, in which the nodes appear completely isolated from each other. The layers are thus “virtual.” Nodes in each cluster are connected locally and can reach synchronization readily, e.g., even in the weakly coupling regime, but there is no inter-cluster synchronization. As the coupling becomes stronger, an explosive percolation²³ among the symmetric clusters can occur, leading to a sudden onset of synchronization among the nodes in each virtual layer. Because of the symmetry, the remote synchronization tends to be perfect. However, this does not guarantee inter-layer synchronization because of the lack of any symmetry among the layers. In fact, synchronization among the nodes in a cluster is fragile and can be weakened when the synchronization among the remote, symmetric nodes takes over. At the explosive transition, there is then a “transfer” of synchronization from the clusters to the virtual layers—a metamorphosis. From the point of view of the whole network, global synchronization can arise explosively at the transition but it is only partial. Complete synchronization among all nodes can eventually occur when the coupling is sufficiently strong. To demonstrate these phenomena in a concrete setting, we consider a multichain network with a star symmetry and implement the Kuramoto phase-coupled dynamics on the network—a typical setting for studying explosive synchronization. We develop a theoretical understanding of explosively remote synchronization based on symmetry considerations, and the robustness of this phenomenon and synchronization metamorphosis is established.

To demonstrate the two phenomena in a concrete setting, we consider the classic Kuramoto model,⁶ 

where $ωi$ and $θi$ are the frequency and phase of node $i$⁠, respectively, and $Aij$ is the $ij$th element of the network adjacency matrix $A$⁠: $Aij=1$ if there is a direct link between nodes $i$ and $j$ and $Aij=0$ otherwise. We exploit the star network structure to search for synchronization metamorphosis and explosively remote synchronization, as it was a previously established setting to demonstrate explosive synchronization.^13,16 Given a simple star network of $n$ peripheral nodes [Fig. 1(a)], we extend each peripheral node into a linear chain to generate a concentric chain network, where there are $n$ such chains radiating from the original central (hub) node with the same length $m$⁠, i.e., there are $m$ nodes along each linear chain. Nodes on different chains but with an equal distance from the hub node are circularly symmetrical with respect to each other, and they form a “virtual layer” because they are not directly coupled with each other. The network, thus, has $m$ virtual layers—each with $n$ nodes, as shown in Fig. 1(b). The network has $N=nm+1$ nodes and we label the central hub node with the index one.

To implement the Kuramoto dynamics, we set the initial frequency of the central hub node to be $ω1=4.0$ and those of all other nodes to be $ωi∈(1,1.01)$ for $i=2,3,…,N$⁠. The initial phases are $θi∈(−π,π)$ for $i=1,2,…,N$⁠. These initial frequencies and phases are randomly distributed in their respective ranges with a uniform probability density function. To characterize the synchronous behaviors on the network, we use three types of order parameters: (1) the global order parameter defined as

(2) the layer order parameter

where $sl$ (⁠$l=1,2,…,m$⁠) denotes the nodal set of each layer, and (3) the chain order parameter

where $sc$ (⁠$c=1,2,…,n$⁠) is the nodal set of each chain (excluding the central hub node). The order parameters $r$⁠, $rl$⁠, and $rc$⁠, thus, characterize the synchronization among all, intra-layer, and intra-chain nodes, respectively, where $r,rl,rc∈[0,1]$⁠. In particular, null values of these parameters signify that the nodal dynamics are completely incoherent and there is lack of any synchronization in the network. In the opposite extreme, $r=1$ indicates that all nodes in the network are synchronized, while $rl=1$ (⁠$rc=1$⁠) means that all nodes in a virtual layer (all nodes in a linear chain) are completely synchronized.

Figure 1(c) shows, for a star-chain network with $m=12$ layers and each with $n=10$ nodes, the chain order parameter $rc$ for an arbitrary chain vs the coupling parameter. In the forward direction, perfect synchronization is achieved as $ε$ becomes nonzero. The $rc$ value remains at approximately a constant value until $εc(f)≈0.79$⁠, at which $rc$ drops discontinuously from a near unity value to a lower value, signifying a replacement of complete synchronization along the chain by partial synchronization and certain loss of synchronization. In the backward direction, i.e., as $ε$ decreases from a relatively large value (e.g., two), $rc$ decreases smoothly from a near unity value until $εc(b)≈0.31$⁠, at which $rc$ increases suddenly to a near unity value so that complete synchronization along the chain is restored. The difference between the values of $εc(f)$ and $εc(b)$ signifies a hysteresis loop.

The phenomenon of explosively remote synchronization is demonstrated in Figs. 1(d) and 1(e), where the layer order parameter $rl$ for $l=1$ and $l=12$ (the last layer) vs $ε$ is shown. In the forward direction, a transition to explosive synchronization occurs at $εc(f)$⁠, where $rl$ immediately reaches the maximal value one as $ε$ increases through $εc(f)$⁠. This is quite striking as perfect synchronization emerges among nodes that are not directly coupled. Figures 1(d) and 1(e) also show that, in the backward direction, $rl$ drops suddenly from one to some near-zero value as $ε$ decreases through $εc(b)$⁠. Compared with the behavior of the chain order parameter $rc$ in Fig. 1(c), the hysteresis loops in Figs. 1(d) and 1(e) are nearly perfect in the sense of sudden unity changes in the order parameter.

Figures 1(c)–1(e) present evidence for the phenomenon of synchronization metamorphosis. In particular, when explosively remote synchronization sets in at $εc(f)$⁠, as shown in Figs. 1(d) and 1(e), $rl$ increases from a near -zero value to a near unity value, while the value of $rc$ decreases from one to a value about 0.8 [Fig. 1(d)], indicating a weakening of the synchronous behavior along the chain. That is, remote synchronization along the circumferential direction is achieved upon an infinitesimal change in $ε$ at the loss of certain degree of synchronization in the radial direction. At $εc(f)$⁠, there is then a sudden jump in the group of synchronized nodes, signifying a synchronization metamorphosis. As shown in Figs. 1(c)–1(e), metamorphosis in the backward direction is more pronounced because, as $ε$ decreases through $εc(b)$⁠, remote synchronization is lost almost completely but synchronization along the chain becomes near perfect, both occurring in an abrupt and explosive manner. Because of the circular symmetry of the network structure, the phenomenon of synchronization metamorphosis occurs with respect to any circular layer and any linear chain in the radial direction.

Figure 1(f) shows the global order parameter $r$ vs $ε$⁠. Globally, in the forward direction, there is a transition to explosive synchronization as $ε$ increases through $εc(f)$⁠, where $r≳0$ for $ε<εc(f)$ and $r$ increases abruptly to a large value of about 0.8 for $ε=εc(f)+0$⁠. As $ε$ increases further from $εc(f)$⁠, $r$ approaches the maximally possible value one. Note that, global synchronization as characterized by the order parameter $r$ can never be perfect in the sense that the value of $r$ can never reach one for any finite coupling, but both local cluster synchronization along a chain and remote synchronization in a virtual layer can be perfect as the values of $rc$ and $rl$ can reach unity. The hysteresis loop exhibited in Figs. 1(c)–1(f) is typical of explosive synchronization.^8–15 

To visualize the actual pattern of remote synchronization among the phase oscillators, we display in Figs. 2(a) and 2(b) all the phase variables of the network for $ε=0.79$ (right after the onset of explosively remote synchronization in the forward direction) and $ε=0.31$ (right before the destruction of the remote synchronization in the backward direction), respectively. Various plateaus represent the nearly constant phases in different virtual layers, i.e., remote synchronization.

The results in Fig. 1 demonstrate that a synchronization metamorphosis occurs in both the forward and backward directions but at a different critical point, leading to a hysteresis loop in the metamorphosis. The parameter interval of the hysteresis loop is given by $Δε=εc(f)−εc(b)$⁠. How does $Δε$ depend on the network structural parameters $m$ and $n$? Figure 3(a) shows, for fixed $m=12$ circular layers, $Δε$ vs $n$⁠, the number of linear chains. It can be seen that, while the values of $εc(f)$ and $εc(b)$ depend on $n$⁠, the size of the hysteresis interval $Δε$ is relatively large and hardly changes with $n$⁠, suggesting that the metamorphosis hysteresis is robust. Similarly, for a fixed value of $n$⁠, e.g., $n=10$⁠, $Δε$ maintains at a large value, regardless of how many circular layers hosting remote synchronization are in the network, as shown in the inset in Fig. 3(a). Globally, for the whole network, the metamorphosis leads to explosive synchronization with a hysteresis loop in the global order parameter, which is robust with respect to variations in the number of chains and virtual layers, as exemplified in Fig. 3(b) for fixed $m=12$ and in Fig. 3(c) for $n=10$⁠, respectively.

We demonstrate that remote synchronization can occur in networks with a symmetry. Specifically, we consider a simple network of seven nodes and six edges with the standard Kuramoto phase dynamics, as illustrated in Fig. 4(a), where there are three pairs of symmetric nodes: {2,3}, {4,5}, {6,7}. We set the initial frequencies of the seven nodes as $ω={2.305,1.102,1.102,0.803,0.803,0.506,0.506}$ and choose their initial phase from the interval $(−π,π)$⁠. Figure 4(b) shows, for $ε=0.35$⁠, the phase evolution of the seven nodes. For this relatively weak coupling, the network as a whole is not synchronized. However, the two nodes in each symmetric pair are synchronized, which is in fact remote synchronization, as the nodes are not directly connected. For this coupling, the nodes that are directly connected are not synchronized. For slightly stronger coupling, e.g., $ε=0.8$⁠, the network achieves global phase synchronization, in addition to cluster synchronization of each symmetric pair, as shown in Fig. 4(c). For larger coupling, e.g., $ε=5.0$⁠, there is complete synchronization of all nodes in the network, as shown in Fig. 4(d).

To demonstrate the generality of explosively remote synchronization in networks with a symmetry, we study the synchronization dynamics of the Cayley tree, as shown in Fig. 5(a). A node is in the virtual layer $m$ when the path length from this node to the central node is $m$⁠. For a tree network with four virtual layers, the numbers of nodes in the layers are 3, 6, 12, and 24, respectively. The initial frequency of the central node is set to be $ωH=0.47$ and the initial frequencies of the nodes in the first to the fourth virtual layers are $ω1∈[1.3,1.31]$⁠, $ω2∈[1.6,1.61]$⁠, $ω3∈[2.0,2.01]$⁠, and $ω4∈[3.9,3.91]$⁠, respectively. Figure 5(b) shows the global order parameter $r$ vs the coupling parameter $ε$⁠, which exhibits an explosive transition to partial global synchronization in the forward direction and a hysteresis. Figures 5(c)–5(e) show the layer order parameter $rl$ vs $ε$ for the virtual layers $m=2,3$ and 4, respectively, all exhibiting an explosive transition to perfect remote synchronization among the nodes in the layer that are not directly connected. While synchronization within each layer is perfect in the sense of $rl=1$⁠, there is no complete interlayer synchronization, so global synchronization in the entire tree network is only partial, i.e., $r<1$ at the onset of explosively remote synchronization.

We exploit the similarity between explosive percolation²³ and explosive transition to synchronization.²² In the forward direction, for $ε≲εc(f)$⁠, nodes within each chain have already achieved near perfect synchronization but there is no interchain synchronization, leading to cluster synchronization. At $εc(f)$⁠, the dynamics of these symmetric clusters merge to form a giant synchronous component in the dynamical sense. Globally, the value of $εc(f)$ is far from being sufficiently large, so global synchronization in the giant component is not perfect. If remote synchronization among the symmetric nodes in different chains is perfect, then synchronization in the chair direction, i.e., within each cluster, must be weakened, leading to a metamorphosis. The question is why the remote synchronization can be perfect. The answer lies in the symmetry of the network, which guarantees robust synchronization between symmetric nodes, even when there is no direct coupling among them, at the explosive onset.

As the coupling parameter approaches the critical value $εc(f)$⁠, the phase difference between any pair of oscillators is small, so the Kuramoto system can be linearized as

Consider a pair of symmetric nodes with the same degree: $i$ and $m$⁠, which are not directly connected. The evolution of the phase variable $θm$ is governed by a similar equation,

At the onset of synchronization, we have $θi˙−θ˙m=0$⁠, so

where $ki$ is the degree of node $i$ and $m$⁠, $ni≤ki$ is their degree after removing the common neighboring nodes, and $ij$ and $mj$ denote the neighboring sets of nodes $i$ and $m$⁠, respectively. Since nodes $i$ and $m$ are symmetric, there is an automorphic permutation $π$⁠: $m=π(i)$⁠, with $P$ as the corresponding automorphic matrix that satisfies $P⋅A=A⋅P$⁠. Let $e$ and $f$ be another pair of symmetric nodes: $e=π(f)$⁠. Permuting all symmetric nodes leads to

and

The identity $P⋅A=A⋅P$ stipulates that if $f$ is a neighbor of $i$⁠, then $e$ must be a neighbor of $m$⁠. That is, the neighbors of symmetric nodes are also symmetric to each other. Denoting the neighbors of $i$ and $m$ as $ij$ and $mj$ (⁠$j=1,…,ni$⁠) and using $θ˙ij−θ˙mj=0$⁠, we get

for $j=1,…,ni$⁠, where nodes $(ij,q)$ and $(mj,q)$ for $q=1,…,nij$ denote the sets of neighboring nodes of $ij$ and $mj$⁠, respectively. The Kuramoto equations for those neighboring nodes and in fact the equations for all nodes in the whole symmetric motif can be written down in a similar way. At the onset of synchronization, the terms on the left hand side of Eq. (7) are zero. Say there are $p$ pairs of symmetric nodes in this symmetric motif that constitute two symmetric groups of nodes, denoted as $G1$ and $G2$⁠, respectively. Converting the Kuramoto equations of all the symmetric nodes in $G1$ and $G2$ into the form of Eq. (7) and combining them lead to

where $L$ is the $p×p$ Laplacian matrix of $G1$ or $G2$ and $B$ is a diagonal matrix of elements $yl$ (⁠$l=1,…,p$⁠) with $yl$ being the number of common neighbors between this node and its symmetric counterpart, and $X$ is a $p×1$ vector of the phase differences between the symmetric nodes. If the determinant of $(L+B)$ is not zero, the only solution will be $X=0$⁠, i.e., any pair of symmetric nodes are perfectly synchronized in spite of the absence of direct links between them.

To prove $|L+B|≠0$⁠, we begin by performing elementary row and column transformations of the matrix $(L+B)$ and changing the rows and columns so that $yl≠0$ appears in the first row and the first column of the matrix. If there are an additional row and column with $yl≠0$⁠, we place them as the second row and the second column, and so on, to obtain the following matrix:

where $rij=0$ or $−1$⁠. Consider the case where $y1′≠0$ is the only nonzero value among all the $yl′$ (⁠$l=1,…,p$⁠), we can convert $R$ into a Laplacian matrix through row and column additions.

We first treat the case of $y1′=1$⁠. Adding a new node to the graph $G1$ or $G2$⁠, we obtain a new network with the following Laplacian matrix:

We can express the matrix in (10) in the form

where $α=[1,−1,0,0,…,0]T$⁠. Because the new network so obtained is connected, this matrix has only one trivial eigenvalue and its rank is $p$⁠. We transform the matrix $Q$ by adding all rows to the first row and all columns to the first column to obtain a matrix in the form $[000R]$⁠. The rank of this matrix being $p$ implies $|R|≠0$⁠.

For the more general case of $y1′=q$⁠, we expand the matrix in the same way to obtain the new Laplacian matrix as

where $0=[0,0,0,…,0]T$⁠, $β=[r12,r13,…,r1p]T$⁠, and the matrix $Qp+q$ has rank $p+q−1$⁠. Adding all the rows in the matrix (11) as the first row and designating the sum of all the columns as the first column lead to a matrix $[000Qp+q−1]$ of rank $p+q−1$⁠. We then remove the first $q−1$ rows of matrix $Qp+q−1$ to get a matrix with rank $p$⁠,

Consider the linear equation $M⋅X=0$⁠, where $X$ is a $p+q−1$ dimensional column vector with $p+q−1$ unknowns, which can be rewritten as

For fixed $x1,x2,…,xq−1$⁠, Eq. (13) has no free variables, giving a unique solution so that the matrix

has the full rank $p$⁠.

The above proof indicates that in both cases (⁠$y1′=1$ or $y1′=q$⁠), the matrix $R$ has full rank. It, thus, suffices to consider the case of $yl′≠0$⁠; then, $yl′=1$⁠. Let $yl′=1$ for $l=1,2,…,K$⁠. We expand the matrix $R$ to have the following Laplacian form:

where $rj=[rjK+1,rjK+2,…,rjp]T$⁠. Performing the same operations on matrix (14) as did matrix (10), we have that $|R|≠0$⁠.

Taken together, through elementary column and row transforms, we have proved that $|R|≡|L+B|≠0$⁠, so the only solution of Eq. (7) is zero. Physically, this means that, when the whole Kuramoto network reaches partial global synchronization, the phases of the symmetric nodes that are not directly connected reach the same value abruptly, leading to the numerically observed explosive onset of perfect remote synchronization with unity layer order parameter. In the  Appendix, we present two concrete examples to demonstrate the matrix operations.

In nonlinear dynamical systems, other metamorphic phenomena can arise such as basin boundary metamorphoses.^24–26 The main contribution of this work is the discovery of the phenomenon of synchronization metamorphosis in nonlinear dynamical networks, where an explosive transition to synchrony and remote synchronization, two previously separately studied phenomena, can occur simultaneously in a metamorphic manner. The key is symmetry, where its interplay with nonlinear dynamics leads to the observed explosive onset of remote synchronization. Since perfect global synchronization among all nodes can be achieved only for infinite coupling, the sudden emergence of perfect remote synchronization among nodes belonging to a symmetry group must occur at the expense of the deterioration of previously perfect synchrony among a distinct set of nodes, leading to a synchronization metamorphosis.

To demonstrate synchronization metamorphosis in a concrete setting, we have focused on a prototype of the dynamical network: a multichain network with a star symmetry. The processes on the network are assumed to be the Kuramoto phase-coupled dynamics. A synchronization metamorphosis occurs at the critical coupling values in both the forward and backward directions. Specifically, a perfect remote synchronization in every virtual layer is achieved at a limited loss of synchronization in the chain direction in the forward direction, and the remote synchronization in every layer is lost but the synchronization along the chain becomes perfect in the backward direction, both occurring in an explosive manner. We have also gained a theoretical understanding of explosively remote synchronization based on symmetry considerations and analyzed the robustness of synchronization metamorphosis.

The work at the East China Normal University was supported by the National Natural Science Foundation of China (NNSFC) under Grant Nos. 82161148012, 11975099, and 11835003 and the Science and Technology Commission of Shanghai Municipality under Grant No. 14DZ2260800. The work at the Arizona State University was supported by the AFOSR under Grant No. FA9550-21-1-0438.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We present two examples to demonstrate the matrix operations in the mathematical proof in Sec. III. In Fig. 6(a), the Kuramoto equations of nodes 2 and 3 are

Setting $θ˙2−θ˙3=0$ leads to

Similarly, we get

At the onset of synchronization, the terms on the left hand sides of Eqs. (A3)–(A5) are zero. Combining Eqs. (A3)–(A5), we obtain the following matrix equation:

To prove $|R|≠0$⁠, we expand $R$ into the following new matrix:

where $Q$ is the Laplacian matrix of the new network in Fig. 6(b). Since the new network is connected, $Q$ has only one trivial eigenvalue and its rank is three. We transform the matrix $Q$ by adding all rows to the first row and all columns to the first column to obtain a matrix in the form $[000R]$⁠. The rank of this matrix being three implies $|R|≠0$⁠, that is, $θ2=θ3,θ4=θ5, and θ6=θ7$⁠.

For the network in Fig. 6(c), using the Kuramoto equations of symmetric nodes, we get the matrix in a similar fashion,

We expand $R$ as

where $Q6×6$ is the Laplacian matrix of the new network in Fig. 6(d), so the rank of $Q6×6$ is five. Adding all rows in matrix (A9) as the first row and designating the sum of all the columns as the first column lead to the matrix $[000Q5×5]$ of rank five. We then remove the first two rows of $Q5×5$ to get the following matrix of rank three:

The linear equation $M⋅X=0$ with $X≡[x1,x2,x3,x4,x5]T$ can be rewritten as

For fixed $x1$ and $x2$⁠, Eq. (A11) has no free variables, so there is a unique solution and the matrix $R$ has full rank: $θ2=θ3,θ4=θ5, and θ6=θ7$⁠.