How to grow an oscillators’ network with enhanced synchronization

Many natural and artificial oscillator systems must be in a synchronized state to carry out their functions properly, including power-grid networks Filatrella, Nielsen, and Pedersen (2008), heart cells, neural networks, biological circadian clocks Acebrón et al. (2005), etc. For stable operation, maintaining the synchronized state is important for such networks even in the presence of newcomers. However, few studies have been conducted on how an added oscillator affects the coherent behavior of the system Wang et al. (2011). In this context, it is worth studying the optimal condition for a newly added oscillator to enhance the synchronization of a growing network.

The Kuramoto model is one of the popular models that describe the synchronization phenomenon between coupled oscillators. The original Kuramoto model describes a system of free oscillators at their natural frequency with all-to-all pairwise interaction between them Kuramoto (2003); Acebrón et al. (2005). In many real systems, though, the interaction is inhomogeneous. In this case, the system can be described as Kuramoto oscillators in a network, where the nodes and links correspond to the oscillators and the interactions, respectively. One of the significant questions is how the topology of the interaction influences the degree of synchronization Huang et al. (2006); Gómez-Gardenes, Moreno, and Arenas (2007); Gómez-Gardeñes, Moreno, and Arenas (2007); Arenas et al. (2008); Gómez-Gardenes et al. (2011); Skardal et al. (2013).

In attempts to answer this question, several studies have investigated the effect of the interaction topology on synchronization in various network structures such as small-world networks Niebur et al. (1991); Hong, Choi, and Kim (2002); Barahona and Pecora (2002), scale-free networks Moreno and Pacheco (2004); McGraw and Menzinger (2007); Hong, Park, and Tang (2007), regular lattices Strogatz and Mirollo (1988); Hong, Park, and Choi (2004, 2005); Hong et al. (2007), and motifs Vega, Vázquez-Prada, and Pacheco (2004). These studies have shown that the topological properties of interactions significantly influence the nature of the synchronization phase transition. To understand this correlation between network structure and dynamics analytically, various mean-field approaches have been developed Ichinomiya (2004); Restrepo, Ott, and Hunt (2005); Ichinomiya (2005); Lee (2005); Restrepo, Ott, and Hunt (2006).

In particular, previous studies Skardal, Taylor, and Sun (2014); Pinto and Saa (2015); Lei et al. (2022) have suggested ways to find the optimal structure or the natural frequency distribution by introducing objective functions obtained from macroscopic order parameters in the strong coupling limit. They have shown that by searching the structure or frequency distribution that lowers the objective functions, one can approach the optimal condition. However, these methods were only applied to systems with a fixed size.

In this study, we demonstrate that one can utilize the same approach as in Refs. Skardal, Taylor, and Sun, 2014; Lei et al., 2022 to find an effective condition for enhancing synchronization in growing systems. To be precise, we find the optimal natural frequency of a newly added oscillator as the network evolves though a predetermined process. This paper is organized as follows. In Sec. II, we explain the model and our main goal, and in Sec. III we present our main results. The results are applied to simple tractable cases in Sec. IV and are numerically verified in more complicated growing networks in Sec. V. In Sec. VI, we give the conclusions of this paper.

We consider partially coupled Kuramoto oscillators. The interaction structure is revealed through the network topology, where dynamic nodes correspond to oscillators and links indicate pairs of the interacting nodes. To take into account system growth, we consider the interaction network to grow under a predetermined rule. Here, we suppose that only a single node is added at each growing step [see Fig. 1 (a)]. In this setup, our purpose is to find the optimal natural frequency of the new node that leads to the highest degree of synchronization.

After a single growing step, the system’s equation of motion is given by the Kuramoto equation Kuramoto (2003),

where $\theta_{i}(t)$ and $\omega_{i}$ are the phase and the natural frequency of the $i$-th oscillator, respectively, for $i=0,1,...,N$ with the initial system size $N$, $K$ is the coupling strength, and $a_{ij}$ is the adjacency matrix characterizing the interaction network. Here, the index $i=0$ is used to indicate the newly added node.

To clarify our purpose, we need to consider the macroscopic measure of the degree of synchronization, namely the Kuramoto order parameter defined as

With this, we express our main goal as finding the $\omega_{0}$ with which the order parameter $r$ takes its maximum.

In the next section, we find the optimal $\omega_{0}$ in the strong coupling regime $K\gg|\omega_{i}|$, which gives an approximate solution. The previous study Lei et al. (2022) showed that one may use alternative order parameters to yield better approximations. Here, we adopt an additional macroscopic measure, referred to as the total local order parameter [see Fig. 1 (c)], given by

The local order parameter $\tilde{r}_{i}$ was originally introduced to understand the effect of the network topology on the formation of synchronized clusters in the relaxation process Arenas, Diaz-Guilera, and Pérez-Vicente (2006a, b); Almendral and Díaz-Guilera (2007).

We want to remark that the conventional order parameter $r$ only takes into account the phase information of the nodes, and thus the information of the network topology is implicitly contained in the steady-state values of the phases [see Fig. 1 (b)]. In contrast, the total local order parameter, being represented as the sum of the differences in pairs of phases over all links in the network [see Fig. 1 (c)], contains the topological information directly. Later in the paper, we show that, interestingly, the $\tilde{r}$-based results outperform $r$-based results in many cases.

Our first task is to express the order parameters as a function of $\omega_{0}$. For this purpose, we have to find the relation between $\theta_{i}(t)$ and $\omega_{i}$. By plugging the stationary state condition $\theta_{i}(t)=\theta_{i}^{\ast}+\bar{\omega}t$ into the equation of motion in Eq. (1), we obtain the relation

with the average frequency $\bar{\omega}=\sum_{i=0}^{N}\omega_{i}/(N+1)$. By substituting the solution for $\theta_{i}^{\ast}$ into Eqs. (2) and (3), one may find the desired forms. Since it is impossible to solve this relation, we take the strong coupling limit. Without loss of generality, we set our reference frame such that $\sum_{i=0}^{N}\theta^{\ast}=0$, which leads to $|\theta_{i}^{\ast}|\ll 1$ in the strong coupling limit. The relation is then linearized as

where $L_{ij}=k_{i}\delta_{ij}-a_{ij}$ represent the elements of the Laplacian matrix $\mathbf{L}$ of the network with the degree $k_{i}=\sum_{j=0}^{N}a_{ij}$ of node $i$ and Kronecker delta $\delta_{ij}$.

This system of linear equations is undetermined because $\mathbf{L}$ has a null eigenvalue. Thus, we have to utilize a generalized inverse, called the Moore–Penrose inverse $\mathbf{L}^{\dagger}=\sum_{\nu=1}^{N}\lambda_{\nu}^{-1}\bm{v}_{\nu}\bm{v}_{\nu}^{\rm T}$, where the superscript ${\rm T}$ stands for the transpose and $\lambda_{\nu}$ and $\bm{v}_{\nu}$ are the eigenvalues in ascending order and the corresponding eigenvectors, respectively, such that $\mathbf{L}\bm{v}_{\nu}=\lambda_{\nu}\bm{v}_{\nu}$ for $\nu=0,1,...,N$. By denoting the vector of the natural frequencies and the steady-state phases by $\bm{\omega}=(\omega_{0},\omega_{1},...,\omega_{N})^{\rm T}$ and $\bm{\theta}^{\ast}=(\theta_{0}^{\ast},\theta_{1}^{\ast},...,\theta_{N}^{\ast})^{\rm T}$, we obtain

In the strong coupling regime, the order parameters up to the leading order are given by

By substituting Eq. (6), we obtain

with the objective functions in the same form as Refs. Skardal, Taylor, and Sun, 2014; Lei et al., 2022,

Here, we used the relation $\mathbf{L}^{\dagger}\mathbf{L}\mathbf{L}^{\dagger}=\mathbf{L}^{\dagger}$. The objective functions have been coined as the adjusted Lyapunov function Lei et al. (2022) (ALF) for $p=1$ and the synchrony alignment function Skardal, Taylor, and Sun (2014) (SAF) for $p=2$. In the strong coupling regime, the optimal synchronization is attained when the objective functions are minimized. As the minimum point $\omega_{0}^{\ast}={\rm argmin}_{\omega_{0}}[O_{p}(\omega_{0})]$ deviates from the optimal value when the coupling strength is not sufficiently large, we regard $\omega_{0}^{\ast}$ as an approximate solution.

In order to find an available form of $\omega_{0}^{\ast}$, it is convenient to represent the matrices and vectors in block form,

with $\tilde{A}_{ij}=a_{0i}\delta_{ij}$, $\tilde{\bm{\omega}}=(\omega_{1},\omega_{2},...,\omega_{N})^{\rm T}$, and the Laplacian matrix $\tilde{\mathbf{L}}$ of the prior network before the growing step. Here, the tilde is used to indicate vectors and square matrices with the dimension reduced by $1$.

While the inverse of any invertible block matrix can be easily obtained from the block inversion formula, it is complicated to find a general formula for the Moore–Penrose inverse of any block matrix. Thus, we utilize specific properties of $\mathbf{L}$ to find $\mathbf{L}^{\dagger}$: (i) its lowest eigenvalue is $\lambda_{0}=0$, (ii) the corresponding null vector is $(1/\sqrt{N+1})\bm{1}$, and (iii) all the other eigenvalues are positive, $\lambda_{i}>0$ for $i>0$, unless the network is separated. These properties show that we can get rid of the singularity of $\mathbf{L}$ by adding $\varepsilon\mathbf{I}$ with any positive number $\varepsilon$ and the identity matrix $\mathbf{I}$. Thus, by subtracting the lowest eigenvalue term from $(\mathbf{L}+\varepsilon\mathbf{I})^{-1}$ and then taking the limit as $\varepsilon\rightarrow 0$, we can find the Moore–Penrose inverse

As derived in Appendix A, the formula gives

where the auxiliary matrices are defined by

for any integer $p$ with a vector of zeros $\tilde{\bm{0}}=(0,0,...,0)^{\rm T}$ and a matrix of ones $J_{ij}=1$.

This result yields the explicit form of the objective functions in terms of $\omega_{i}$. For convenience, we represent the objective functions as a function of $\bar{\omega}$ rather than $\omega_{0}$. Then they are written as quadratic functions of $\bar{\omega}$,

Their minimum is taken when

This optimal condition is our main result. We note that the coefficients are normalized, $\tilde{\bm{c}}_{p}^{\rm T}\tilde{\bm{1}}=1$. The condition shows that the optimal $\omega_{0}$ is simply given by the sum of all contributions from each node, which are factorized into the product of structure-dependent coefficients and individual natural frequencies. This reveals that $\tilde{\bm{c}}_{p}$ reflects the strength and the tendency of each prior node’s contribution to the optimal value.

For $p=1$, we can understand the meaning of the coefficients $\tilde{\bm{c}}_{1}$ by considering an auxiliary dynamics where all the oscillators are identical with the same natural frequency denoted by $f_{i}=f$ except for the added node having a different frequency, $f_{0}=f+\Delta f$. We denote the phase of these identical oscillators by $\phi_{i}$ to highlight that they are not of the original system. We remark that this setup resembles systems with an extraordinary oscillator, referred to as a pacemaker Kori and Mikhailov (2004). In such a system, however, the interaction between the pacemaker and the other nodes is non-reciprocal.

At the steady-state, Eq. (5) gives

with $\bar{f}=f+\Delta f/(N+1)$. By using the relation $\tilde{\mathbf{A}}\tilde{\bm{1}}=(\tilde{\mathbf{L}}+\tilde{\mathbf{A}})\tilde{\bm{1}}$, one can recast the lower block as

with the relative phase $\Delta\tilde{\bm{\phi}}^{*}=\tilde{\bm{\phi}}^{*}-\phi_{0}^{*}\tilde{\bm{1}}$ with respect to the distinctive oscillator’s phase. We note that the relative phase $\phi^{*}_{i}-\phi^{*}_{0}$ is a topological node property independent of any dynamical properties. Small $\Delta\tilde{\bm{\phi}}^{*}$ means that the $i$-th node is prone to synchronize with the new node, intrinsically. In this sense, one can regard $\Delta\phi^{*}_{i}$, or equivalently the coefficients $c_{1,i}$, as a measure of the synchronization anti-affinity between node $i$ and $0$.

To demonstrate the validity of our result, we apply it to simple solvable examples: identical oscillators in an arbitrary network, one-to-all connection in an arbitrary network, and two cliques connected by a single link.

As a next step, we examine the validity of our methods for two popular growing networks: a random growing network and the Barabasi–Albert (BA) network. We construct a random growing network by connecting a new node to the others with probability $P=0.1$. If the new node is connected to none of the existing nodes, we reject this step and try again until the new node is coupled to at least one of the existing nodes. This constraint leads the growing random network to have different characteristics from the Erdos–Reny network as reported in Ref. Callaway et al., 2001. Meanwhile, we construct a BA network by using the well-known preferential attachment algorithm Barabási and Albert (1999).

To see how the system size influences the performance of the SAF and ALF methods, we increase the system size iteratively starting from an initial seed network of size $N=30$. We build the seed network by growing the system from a single node in accordance with the growing rule, and then assigning random numbers from the standard Gaussian distribution $\mathcal{N}(0,1)$ as the natural frequencies of the initial nodes. After that, the natural frequencies of the added nodes are taken to be the optimal $\omega_{0}^{\ast}$ obtained from Eq. (16).

To demonstrate the performance of our results, we compare them to other plausible strategies: global average (GA), local average (LA), and random assignment (RA). The added node’s natural frequency is chosen to be $\sum_{i=1}^{N}\omega_{i}/N$ for GA, $\sum_{i=1}^{N}a_{0i}\omega_{i}/k_{0}$ for LA, and a random number sampled from the same distribution as the initial natural frequency distribution for RA. We remark that the GA method corresponds to the optimal choice obtained from the objective function $O_{p}$ with $p=0$.

Figure 3 shows the system-size dependence of the order parameter $r$ evaluated numerically in the successive growing process. Each data point is obtained from the average over 300 seed networks. We used the fourth-order Runge-–Kutta method (RK4) with a time interval of $dt=0.001$. The relaxation time is set to be $5\times 10^{5}$ steps prior to growth and $2\times 10^{5}$ steps during the growing process.

The upper and lower panels in Fig. 3 exhibit the order parameter changes in the random growing and BA network, respectively. Each plot corresponds to different values of coupling strengths such that the initial $r$ has different values of approximately $0.4$, $0.65$, and $0.9$ (left to right). Although in principle our method may show poor performance away from the strong coupling regime, the ALF method outperforms the others even in cases with the lowest coupling strength, whereas the SAF method underperforms the GA method in the random growing network. This implies that the performance of the SAF method, as in the previous two-clique network example, is more sensitive to the interaction structure. The numerical results confirm that the ALF method provides an effective approximation in a wide range of coupling strengths.

To demonstrate the robustness of the ALF method against the coupling strength, we plot the order parameter at the terminal size as a function of the coupling strength $K$ in Fig. 4. Here, the columns represent different initial and final system sizes. For comparison, we add black dashed curves to indicate the initial order parameter. The left panels show that the GA method achieves the highest performance when the system size and the coupling strength are small. But in most cases, the average-based methods are less effective, and even have a harmful influence on synchronization in the BA networks. In contrast, the ALF and SAF methods display strong performances in both network types for different coupling constants and seed network sizes. This certifies the effectiveness and reliability of our suggested approaches.

In the present paper, we investigated the effect of system growth on synchronization. We have found that the addition of a new oscillator can lead to a synchronization enhancement when we assign the new oscillator with a proper natural frequency. Furthermore, we have suggested an effective and useful approximation of the optimal natural frequency by using a linearized dynamics in the strong coupling constant regime. We have applied our approximations to various examples and found that the proposed methods outperform other plausible methods based on the statistics of the existing nodes’ natural frequencies. Between the two approaches, the ALF method showed robust performance even far from the strong coupling regime.

We also observed that the performance of the ALF method is the most insensitive to the interaction structure. In particular, with stronger network heterogeneity, the other methods show relatively poorer performance. Thus, we can conclude that the ALF method is the best strategy in the sense that it leads to a notable amount of synchronization enhancement regardless of the interaction topology in a wide range of coupling strengths.

Moreover, our results revealed that the approximate optimal value is represented as the linear combination of the natural frequencies of the existing nodes. From this, we argued that the linear coefficients can be regarded as an intrinsic node property indicating a certain type of synchronization anti-affinity with respect to the added node. It would be an interesting future work to compare this quantity with the dynamical connectivity matrix Arenas, Diaz-Guilera, and Pérez-Vicente (2006a, b).

It would also be an interesting future work to apply our findings to real systems described by swing equations, such as power-grid systems. In this case, we expect that it is straightforward to extend our method to second-order Kuramoto models. We believe that our work will pave the way for research into the synchronization of growing systems.

Finally, we want to remark that our approaches can also be utilized to study the synchronization of fixed-size systems. For example, by regarding each underlying node as a newcomer, one may examine whether its natural frequency is near the optimal value. Through this approach, we would learn how to adjust the natural frequency to enhance synchronization.

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