Condensation of Synchronous Orbits on Complex Networks

In this paper we study frequency synchronization of Kuramoto oscillators. We ﬁnd a typical phenomenon of condensed synchronous orbits on single-layer or duplex networks through statistical mechanics analysis and numerical simulations, where the distribution of synchronous orbits is in a bell-shaped form. Further, we investigate phase synchronization on single-layer and duplex networks with diﬀerent distributions of inherent frequencies. We ﬁnd that normally distributed inherent frequencies with low variances are more beneﬁcial for phase synchronization, and separately distributed inherent frequencies can slow down the synchronization process. In the end, we investigate the inﬂuence of one layer’s inherent frequencies on the other layer’s phase synchronization through inter-layer couplings. Interestingly, we ﬁnd that one layer’s inherent frequencies with a highly condensed distribution can greatly improve phase synchronization on the other layer. The results shed new lights to our understanding of the nature of synchronization on single-layer as well as multilayer complex networks of coupled Kuramoto oscillators.


Introduction
Over the last two decades, network science has developed to become an important subject of scientific studies in the world, including many complex biological, social, technological systems [1,2,3,4,5]. Knowing the importance of interactions among network's agents, current research has been focusing on the dynamical properties of networks [6,7,8,9,10]. On the other hand, the complexity of real-world systems calls for more comprehensive and practical models, such as multilayer networks, to study the intrinsic dynamical properties of networks from different perspectives [11,12,13,14,15,16]. In [2,6,7,9,10,11,14], where synchronization on singlelayer or multilayer networks are studied for various tasks. Besides, from flocking formation to phase synchronization of oscillators [17,18], numerous issues about synchronization on complex networks have been investigated [19].
Speaking of phase synchronization of oscillators, the Kuramoto model has been extensively used [18,20]. Watts and Strogatz used the Kuramoto model to reveal the small-world properties of networks [2]. By employing the Kuramoto model on BA scale-free networks, Moreno and Pacheco found that the nodes with large degrees play an important role in determining network synchronization [21]. In 2018, we studied the Kuramoto model on multiplex spacial networks and found that inter-layer connections play a vital role in global synchronization [22]. Many significant works have been done during the past decade, including the explosive synchronization discovered by Gómez et al. [23] and the Chimera phenomenon studied in [24]. However, the above research are all based on the first-order Kuramoto model. In order to investigate the Kuramoto model with inertia and phase-frequency correlation, a second-order Kuramoto model was used [20,25]. To date, frequency synchronization of second-order Kuramoto oscillators on both single and duplex networks is rarely investigated, which is the topic for study in the present paper.
Generally, the dynamics of systems on a single-layer network can be described by where ( ) = ( ( )), ( ) ∈ R , describes the nodal dynamics for the -th node, is the coupling strength, coupling matrix L = ( ) × is the Laplacian matrix, and : R → R is the inner coupling function. One important method for studying synchronization of this model is the master stability function (MSF) method proposed by Pecora and Carroll [26]. Extensive works have been done in the area of synchronization, structure identification and control on both single-layer and multilayer networks [27,28,29,30]. Due to the complicated nodal dynamics and the innercoupling function, the synchronous orbits and the interplay between the network topology and dynamics can be very complicated. As an alternative, a simpler model is widely used to study the topological effect of a network on its synchronization [31,32], which is given by where ( ) ∈ R describes the nodal states for the -th node, is the coupling strength, and A = ( ) × is the adjacency matrix. The network model (2) is also used for studying consensus or diffusion [33,34,31]. Based on this model, Gómez et al. [11] isolated the topological factor from dynamics and found a method to optimize the structures of multilayer networks for diffusion like synchronization. In this paper, we propose an undamped second-order Kuramoto model, which is similar to Eq. (2), but the couplings between oscillators are varying according to oscillators' phases by following a cosine function. The paper is orgenized as follows. In Section II, we present the undamped second-order Kuramoto model and study the synchronous orbits of oscillators' frequencies. In order to investigate the distribution of the synchronous orbits, in Section III we introduce several methods based on statistical mechanics to analyze the steady state of frequency synchronization. Section IV presents the numerical results of frequency synchronization on single-layer networks, in which the distribution of synchronous orbits follows a bell-shaped form. In Section V, we further explore the synchronous orbits on duplex networks, and also find condensation phenomena we predicted. Finally, we study the influence of inherent frequency distributions on phase synchronization and obtain similar results as that in Section IV. Section VI presents some conclusions.

Undamped Second-order Kuramoto Model
First of all, consider a network governed by the first-order Kuramoto equation where ( ) is oscillator 's phase, 0 is the initial instantaneous or inherent frequency of oscillator , = ( ) × is the adjacency matrix of the network, and is the coupling strength. The following typical order parameter is used to quantify phase synchronization: where ( ) is the average phase of oscillators, and ( ) ∈ [0, 1]. When the coupled Kuramoto oscillators' phases are synchronized, the order parameter will reach the value of = 1, while for the completely incoherent situation, = 0. In order to investigate the instantaneous frequency variation, introduce an undamped second-order Kuramoto model, as follows: which is transformed from Eq. (3) by taking a derivative of ( ) with respect to time, where ( ) and ( ) represent the instantaneous frequencies of nodes and , respectively.
Interestingly, it has the same form as the classic consensus or diffusion problem, except that the Laplacian matrix in Eq. (6) not only is symmetric but also changes with time by correlating with oscillators' instantaneous phases. Because Eq. (6) satisfies the dissipative coupling condition, the instantaneous frequencies of the oscillators will eventually reach consensus or synchronization, as shown in Fig. 1. Furthermore, since the synchronous orbit is sensitive to oscillators' inherent frequencies (initial instantaneous frequencies), the objective here is to find a statistical pattern about the distribution of the synchronous orbits under different inherent values, but with the same distribution.

Statistical Properties of Synchronous Orbits
Complex networks are systems consisting of many coupled nodes. So, the dynamical properties of the networks must follow statistical rules of the nodes in some way. From a thermal dynamics point of view, one may consider a Kuramoto network as a closed ensemble, that is, a canonical ensemble, in which there is no matter but only energy exchange with the environment. This means that the set of nodes obey the law of conservation during the evolution. Consequently, the probability that one system occupies state can be written as where is the energy of a system when the system occupies state , = e − is the partition function, = 1 , is the thermodynamic temperature, and is the Boltzmann constant.
Assume that each synchronous orbit represents a possible state. Since synchronous orbits are represented by synchronized frequencies, the energy of state is related with the corresponding frequency. In general, a system's energy can be described by where is the synchronized frequency of the network, 0 is the average of inherent frequency of the oscillator and is a parameter which is set as 1. In this part, frequency synchronization of the Kuramoto network can be considered as an annealing process, and the synchronous orbit 0 has the lowest energy comparing with other states.
According to Eqs. (7) and (8), one has which is a Gaussian distribution of the synchronous frequency orbits . As can be seen, the expectation of the Gaussian distribution is determined by the oscillator's inherent frequency, and the variance is related with the temperature .

Condensation of Synchronous Orbits on Single-layer Networks
In this section, numerical analysis is performed to investigate frequency synchronization on single-layer networks, which evolve with Eqs. (3) and (5), and to study the influence of inherent frequency distributions of oscillators on phase synchronization. Note that, unless otherwise specified, the inherent frequencies of oscillators' phases are randomly chosen in [0, 2 ] throughout this paper. Figure 1 shows a typical example of synchronous frequency orbits of an Erdös-Rényi (ER) random network [1], composed of 100 nodes, where the inherent frequencies are randomly chosen in [0, 1]. The expectation of inherent frequencies is 0.53268, while the synchronous orbit is located at 0.51967, which means that the synchronous frequency is near the expectation of inherent frequencies rather than at its exact value. According to the analysis in Section III, the locations of synchronous orbits under the same inherent frequency distribution but different specific values may follow a Gaussian distribution.

Distribution of synchronous orbits
Employing the Runge-Kutta algorithm, the frequency synchronization is simulated under different topologies with network size = 1000 for 10000 times, and 10000 synchronous orbits are obtained, where the inherent frequency values are randomly chosen in [0, 1] and the coupling strength is = 1. Figure 2 shows the distributions of these synchronous orbits on ER random, K-regular, Watts-Strogatz (WS) small-world and Barabasi-Albert (BA) scale-free networks [3,2], which are all clearly in the bell-shaped form. The light grey columns show that there is a much lager probability that a synchronous orbit locates on the average value of the inherent frequencies, which is = 0.5. The red curves are the Gaussian fits of columns. Since the probability density diagram of synchronous frequencies is symmetric and condensed, this kind of phenomena are referred to as condensation of synchronous orbits. Obviously, the variance of the Gaussian envelope changes with network topologies. For relatively regular topologies, such as K-regular or WS networks (rewiring probability is 0.05), the variances are very small, while for more random topology, such as ER or BA networks, the variances are much lager. Interestingly, it is found that the emergence of bell-shaped distribution about synchronous orbits is uncorrelated with the distribution of the inherent frequencies. Furthermore, the variance can be influenced by the network topology, coupling strength and network size. Because the variance of the Gaussian envelope is related with the temperature as mentioned in Section III, there could be other unrevealed truths behind this phenomenon.

Impact of inherent frequency distribution on phase synchronization
Because the distribution of synchronous orbits is of a lowvariant bell-shaped form, the condensation patterns of inherent frequencies may have some positive or negative impact on phase synchronization over the network. Therefore phase synchronization on ER networks with different inherent frequency distributions is further investigated, where the expectations are set as 0.5. As shown in Fig. 3, a more condensed distribution of inherent frequencies (smaller variance ) is more beneficial for the oscillators to achieve phase synchronization, compared with uniformly distributed inherent frequencies (the royal blue circle line). Surprisingly, if half the inherent frequencies are 0.1 and the other half are 0.9, which is more like a two-peak distribution (the oranges quare line), the phase synchronization is greatly suppressed.
According to the results shown in Fig. 3, the critical values of coupling strength thresholding non-synchronization and synchronization, denoted as , are plotted in Fig. 4, with respect to different variances of inherent frequency distributions. It is obvious that is linearly correlated with . The result implies that more condensed inherent frequency distributions can better facilitate phase synchronization.

Condensation of Synchronous Orbits on Duplex Networks
In this section, the frequency synchronization is investigated on duplex networks, with dynamicsdescribed by where ( ) and ( ) are the phase and frequency of node located on layer ( ), ( )0 is the inherent frequency of oscillator on layer ( ), intra is the intra-layer coupling strength of layers and , inter is the inter-layer coupling strength between layers, and ( ) = ( ( ) ) × is the adjacency matrix of layer ( ).

Distribution of synchronous orbits
Consider a duplex ER network, where each layer is randomly generated according to the ER algorithm [1,4]. In Fig. 5, it shows that the statistical distribution of 10000 frequency synchronous orbits on duplex ER networks also follows the Gaussian distribution. Here, the network size of each layer is = 1000, and the coupling strength intra = inter = 1. The light grey columns indicate that the condensation phenomenon of synchronous orbits, as presented in Section IV, can also emerge on duplex networks. The red curve shows the Gaussian fitting of columns, where = 0.5 due to the uniformly distributed inherent values of ( )0 in [0, 1].
According to Eqs. (11a), (11b) and [11], frequency synchronization on strongly connected duplex networks will always happen due to the existence of inter-layer coupling and the zero-sum rows in the supra-Laplacian of a duplex network. Similarly, the bell-shaped form of distributed synchronous orbits does not change by varing distribution functions of inherent frequencies, network topology, coupling strength or network size. In particular, these factors can only change the expectations and the variances of the distributions.

Impact of inherent frequency distribution on phase synchronization
Consider the duplex ER network as a whole. Figure 6 shows that different inherent frequency distributions have different influences on phase synchronization on duplex ER networks, where = 100. Similarly to the observations from Fig. 3, more condensed distributions of inherent frequencies, i.e., Gaussian distributions with smaller variances, are beneficial for networks to achieve global phase synchronization, compared with uniformly distributed inherent frequencies (shown by royal blue line-circle). Meanwhile, letting half of the nodes have inherent values 0.1 and the other half 0.9, like a two-peak distribution, one obtains the result shown by the orange square lines. It is obvious that the global phase synchronization is suppressed by these two-peak-like inherent values, which is in accordance with the observation in Section IV. One can further observe from Fig. 6 that duplex ER networks are more beneficial to phase synchronization than single-layer ER networks, in that their critical coupling strengths are smaller than those shown in Fig. 3. In Fig. 7, one can observe that, for Gaussian distributed inherent frequencies with low variances, the critical intralayer coupling strength intra , above which synchronization can be achieved, is proportional to the variance . The red line represents the linear fit of data points, whose slope is approximately half of that in Fig. 4. Additional simulations show that, if a multiplex network has layers, the slope value is 1/ of that for a single-layer network. This interesting observation is worth further investigation.
Furthermore, because global phase synchronization can also be influenced by inter-layer coupling, it is meaningful to investigate the critical intra-layer coupling strength with different inter-layer couplings, with resuts shown in Fig. 8. It can be seen that the critical intra-layer coupling strength is an exponential function of the inter-layer coupling strength inter . The red curve is the exponential fit of the data. When inter > 0.4, intra will reach its lower bound, implying that the synchronous ability cannot be further improved by increasing the inter-layer coupling. Thus, the inter-layer coupling strength is set to inter = 1, to exclude its interference on the following analysis.

State control by altering adjacency layer's inherent frequency distribution
Sometimes, in the real world, a system's state is required to be confined in a particular regime, which can not be manipulated directly. A feasible way is to manipulate an accessible system that is coupled with the underlying system in need of control, as does in [28]. It is to investigate how in this subsection one layer's inherent frequencies affect the other layer's phase synchronization through inter-layer coupling, which has rarely been studied before. Here, a new way is suggested for state control of multilayer networks. Specifically, in a duplex network, the inherent frequencies are uniformly distributed on the target layer, while for the interacting layer (the layer that is coupled with and therefore intervene with the target layer), the inherent frequencies are controlled by probability functions. The objective is to influence phase synchronization of the target layer by altering the inherent frequency distribution of the interacting layer. As shown in Fig. 9, normally distributed inherent frequencies with small variances on the interacting layer can improve phase synchronization on the target layer, compared with uniformly distributed inherent frequencies (shown by royal blue line-circle). Similarly to previously shown results, when half of the nodes' inherent frequencies are set as 0.1 and the other half as 0.9, like a two-peak distribution (shown by orange line-square), the phase synchronization on the target layer is suppressed. Figure 10 shows that, for normally distributed inherent frequencies with low variances, the critical intra-layer coupling strength intra keeps almost unchanged with varying variance . This means that, to improve phase synchronization on the target layer, there is no need to set too condensed inherent frequencies on the interacting layer. Using an appropriate Gaussian distribution can achieve the same goal.
It is natural to ask if one can use a small network to manipulate the target layer to reach the same goal under a restriction of total cost. Here, use = ′ / to represent the restriction of total cost, where ′ and are the sizes of the interacting layer and the target layer, respectively. Set = 1000. Figure  11 displays the order parameter versus the intra-layer coupling strength intra for different sizes of interacting layers. Apparently, a stronger interacting layer (higher cost) is more beneficial for phase synchronization on the target layer. In all the simulations, the inherent frequencies on the interacting layer and the target layer are Gaussian distributed with variance = 0.1 and uniformly distributed, respectively. In this way, phase synchronization can be achieved on one layer by changing the inherent frequencies of its adjacency layer under certain cost restriction.

Concluding remarks
In this paper, we have introduced a modified second-order Kuramoto model to investigate the synchronous orbits on networks from a statistical mechanics perspective, and have predicted that the distribution of the synchronous orbits follows a bell-shaped form. Our numerical results on singlelayer and duplex networks have verified this prediction. We call this interesting phenomenon as condensation of synchronous orbits. Furthermore, we have investigated the influence of inherent frequencies on the phase synchronization on both single-layer and duplex networks. We found that a more condensed distribution of inherent frequencies is more beneficial for phase synchronization, and more dispersedly dis-tributed inherent frequencies can suppress the synchronization process. Finally, we have investigated how one layer's inherent frequencies affect the other layer's phase synchronization through inter-layer couplings. We found that a more condensed distribution of inherent frequencies can greatly improve phase synchronization on the other layer. Our investigation offers a different approach to studying coherent behaviors on complex networks and presents a different way to investigate network control problems.

Code availability
Not applicable