Optimal global synchronization of partially forced Kuramoto oscillators

TL;DR

This paper develops an effective low-dimensional model to analyze and optimize the global synchronization of large networks of Kuramoto oscillators with partial external forcing, providing analytical predictions and numerical validation.

Contribution

It introduces a dimensional reduction approach and an optimization scheme for maximizing synchronization in partially forced Kuramoto networks.

Findings

Optimal forced oscillators are those with natural frequencies farthest from the forcing frequency.

The model accurately predicts critical parameters for synchronization onset.

Numerical simulations confirm the effectiveness of the optimization strategy.

Abstract

We consider the problem of global synchronization in a large random network of Kuramoto oscillators where some of them are subject to an external periodically driven force. We explore a recently proposed dimensional reduction approach and introduce an effective two-dimensional description for the problem. From the dimensionally reduced model, we obtain analytical predictions for some critical parameters necessary for the onset of a globally synchronized state in the system. Moreover, the low dimensional model also allows us to introduce an optimization scheme for the problem. Our main conclusion, which has been corroborated by exhaustive numerical simulations, is that for a given large random network of Kuramoto oscillators, with random natural frequencies $\omega_i$, such that a fraction of them is subject to an external periodic force with frequency $\Omega$, the best global synchronization properties correspond to the case where the fraction of the forced oscillators is chosen to be those ones such that $|\omega_i-\Omega|$ is maximal. Our results might shed some light on the structure and evolution of natural systems for which the presence or the absence of global synchronization are desired properties. Some properties of the optimal forced networks and its relation to recent results in the literature are also discussed.

Equations

$$\frac{d{\theta }_i}{dt}={\omega }_i+1_C(i)F\sin \left(\Omega t-{\theta }_i\right)+\lambda \sum _{j=1}^N{A}_{ij}\sin ({\theta }_j-{\theta }_i),$$

$$r(t)e^{i\psi (t)}=\frac{1}{N}\sum _{j=1}^N{e}^{i{\theta }_j(t)},$$

$${\theta }_i(t)=\alpha (t)+\frac{{\omega }_i}{\Omega }\beta (t)+\Omega t$$

$$r(\beta )=⟨\cos \frac{\beta \omega }{\Omega }⟩=\int d\omega g(\omega )\cos \frac{\beta \omega }{\Omega }$$

$$\psi =\alpha +\Omega t,$$

$$⟨h⟩=\frac{1}{N}\sum _{k=1}^N{h}_k,$$

$$r\approx 1-\frac{⟨{\omega }^2⟩}{2{\Omega }^2}{\beta }^2+\frac{⟨{\omega }^4⟩}{4!{\Omega }^4}{\beta }^4+\cdots ,$$

$$\dot{\alpha }+\frac{{\omega }_i}{\Omega }\dot{\beta }$$

$$\dot{\alpha }=-\Omega -F{\mathcal{I}}_1(\alpha ,\beta ),$$

$${\mathcal{I}}_1(\alpha ,\beta )=\frac{1}{N}\sum _{i=1}^N{1}_C(i)\sin \left(\alpha +\frac{{\omega }_i}{\Omega }\beta \right).$$

$$\dot{\beta }=\Omega -\frac{F{\Omega }^2}{⟨{\omega }^2⟩}{\mathcal{I}}_2(\alpha ,\beta )+\frac{\lambda {\Omega }^2}{⟨{\omega }^2⟩}{\mathcal{I}}_3(\beta ),$$

$${\mathcal{I}}_2(\alpha ,\beta )=\frac{1}{N}\sum _{i=1}^N{1}_C(i)\frac{{\omega }_i}{\Omega }\sin \left(\alpha +\frac{{\omega }_i}{\Omega }\beta \right)$$

$${\mathcal{I}}_3(\beta )=\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^N{A}_{ij}\frac{{\omega }_i}{\Omega }\sin \frac{\beta }{\Omega }({\omega }_j-{\omega }_i)$$

$${\mathcal{I}}_1(\alpha ,\beta )$$

$${\mathcal{I}}_2(\alpha ,\beta )$$

$${\mathcal{I}}_3(\beta )$$

$$\dot{\alpha }$$

$$\dot{\beta }$$

$$r_{\mathrm{n}}=\exp \left(-\frac{⟨{\omega }^2⟩{\beta }^2}{2{\Omega }^2}\right)$$

$$r_{\mathrm{u}}=\frac{\sin \sqrt{3⟨{\omega }^2⟩}\frac{\beta }{\Omega }}{\sqrt{3⟨{\omega }^2⟩}\frac{\beta }{\Omega }},$$

$$\left|\frac{fF}{\Omega }\right|<1.$$

$$F>F_{\mathrm{c}}=\frac{\Omega }{f}.$$

$$\frac{\partial(\dot{\alpha},\dot{\beta})}{\partial(\alpha,\beta)}=-\left(\begin{array}[]{cc}fF\cos\alpha&0\\ 0&{\lambda\left\langle d\right\rangle}+fF\cos\alpha\end{array}\right),$$

$$\frac{\lambda ⟨d⟩}{\Omega }+\frac{fF}{\Omega }<\zeta ,$$

$$\zeta =-\frac{⟨{\omega }^2⟩}{{\Omega }^2\min r^{\prime }(\beta )}>0.$$

$$\dot{\alpha }$$

$$\dot{\beta }$$

$$\mathcal{L}=\frac{{\sum }_{e(i,j)}({\omega }_i-{\omega }_j{)}^2}{{\sum }_{k=1}^N{\omega }_k^2},$$

$$\frac{{\beta }_{\ast }}{\Omega }=\frac{1+\frac{⟨\omega {⟩}_C}{⟨{\omega }^2⟩}\Omega }{fF\left(\frac{⟨{\omega }^2{⟩}_C}{⟨{\omega }^2⟩}-\frac{⟨\omega {⟩}_C}{\Omega }\right)\cos {\alpha }_{\ast }+\lambda \mathcal{L}},$$

$${\theta }_i^{\prime }=\frac{{\omega }_i}{\Omega }+1_C(i)\frac{F}{\Omega }\sin \left(\tau -{\theta }_i\right)+\frac{\lambda }{\Omega }\sum _{j=1}^N{A}_{ij}\sin ({\theta }_j-{\theta }_i),$$

$$M\left(\begin{array}[]{c}\dot{\alpha}\\ \dot{\beta}\end{array}\right)=\left(\begin{array}[]{c}\langle\omega\rangle-\Omega-F{\cal I}_{1}\\ \frac{\langle\omega^{2}\rangle}{\Omega}-\langle\omega\rangle-F{\cal I}_{2}+\lambda{\cal I}_{3}\end{array}\right),$$

$$M=\left(\begin{array}[]{cc}1&\frac{\langle\omega\rangle}{\Omega}\\ \frac{\langle\omega\rangle}{\Omega}&\frac{\langle\omega^{2}\rangle}{\Omega^{2}}\end{array}\right).$$

Full text

We consider here the dynamics of a large number of interacting Kuramoto oscillators. Each oscillator has its own natural frequency $\omega_{i}$, which is assumed to be a random variable, and the interactions among them are associated with the edges of a random network. Despite of being essentially a random system, there are plenty of robust results obtained from this kind of model which have proven to be relevant in many different areas. This is the case, for instance, of synchronization phenomena, see Review1, and Review2, for comprehensive reviews on the subject. Here, we are concerned with one of the main variations of the network of Kuramoto oscillators, the case where some of the oscillators are subject to an external periodically driven force with frequency $\Omega$. By exploring a recently proposed analytical approach, we introduce an optimization scheme for the onset of the so-called global synchronization in the system, a regime where all oscillators rigidly rotate, forming a compact swarm, in the same pace of the external force, with frequency $\Omega$. We show that the best global synchronization properties correspond to the case where the set of forced oscillators is chosen to be those ones such that the value of $|\omega_{i}-\Omega|$ is maximal. Our results may help to understand the evolution and structure of natural systems with many interacting agents for which global synchronization plays an important dynamical role.

I Introduction

Synchronization in complex networks of oscillators is a paradigmatic dynamical problem which has received huge attention recently. Its intrinsically rich dynamics and vast applicability range have motivated a myriad of studies in very different areas, see Refs. Review1, and Review2, for comprehensive reviews on the subject. The mostly used model for synchronization studies is still the well-known Kuramoto oscillator, introduced more than forty years ago kuramoto1975 in the context of chemical oscillations. The list of situations involving interacting individual agents which can be well described dynamically by a network of Kuramoto oscillators is vast, ranging from Biology to mass communication technologiesReview2 . In many cases, the real underlying interaction network is so large and entangled that any intention to model it will necessary require some kind of random description, typically with certain prescribed statistical properties for the nodes and interactions among themNewman . Roughly, this is namely the key concept of a complex network for our purposes: a random graph with certain prescribed statistical properties for its vertices and edges attributes. The connection pattern and the frequencies distribution in a network of Kuramoto oscillators are known to affect its synchronization properties, and since in many applications synchronization is a dynamical regime expected to be favored (or suppressed), many recent works have been devoted to the study of the optimization of synchronization in complex networkspecora1988 ; barahona2002 ; m1 ; m2 ; m3 ; donetti2005 ; dwivedi2015 ; PintoSaa ; optimal ; AISync . Such optimization schemes are useful not only to enhance the synchronization capability of a given complex network by adjusting its connection topology, but also to understand the evolution and structure of natural systems with many interacting agents for which synchronization plays a fundamental role.

One of the main variations of the classical synchronization problem is the case of a network of Kuramoto oscillators subject to external periodically driven forcesp1 ; p2 ; p3 ; p4 ; p5 ; p6 ; p7 . The dynamics become much richer in such a case, with distinct notions of synchronization available and the appearance of possible bifurcations among them. Here, we consider another variant of the classical problem: the case where only a fraction of the oscillators is subject to the external periodical force. There are plenty of motivations in natural systems and many recent examples of application involving networks of partially forced Kuramoto oscillators. In Biology, for instance, there are several situations where the response of specific groups of cells are triggered by external stimuli acting on other cells. Typical examples include the receiving and processing of visual, olfactory, or haptic signals, but the list is certainly much larger. The dynamical models for these cases consist commonly in a complex network where the cells, which pertinent dynamical states are described by Kuramoto oscillators, are located at the vertices, with the edges corresponding to the interaction among them. Recent applications of partially forced Kuramoto oscillators include the study of the response to external stimuli of the C. elegans neural networkmodular ; p8 (see book1, for a review on complex network models for the response of neural networks to external stimuli), self-organization and pattern formation in the growing and developing of vertebratesp9 , auditory signals in amphibiansp10 , and several issues on the circadian rhythmp11 ; p12 ; p13 .

The problem we are concerned here corresponds to a connected network with $N$ Kuramoto oscillators governed by following dynamical equations

[TABLE]

where $\theta_{i}$ stands for the dynamical state (phase) of the Kuramoto oscillator with natural frequency $\omega_{i}$, which is assumed to be located at the $i^{\rm th}$-node of the underlying network. The connections among the oscillators are represented by the usual undirected adjacency matrix with components $A_{ij}$, and $\lambda$ stands for the uniform coupling strength among the oscillators. We will consider here the so-called attractive case, for which $\lambda>0$. The subset of $N_{C}<N$ forced nodes is denoted by $C$, and $\mathbf{1}_{C}(i)$ is its indicator function, i.e., $\mathbf{1}_{C}(i)=1$ if $i\in C$, or zero otherwise. The external driven force is also assumed to have uniform intensity $F$ and frequency $\Omega$ for all nodes $i\in C$. Without loss of generality, we can assume $F\geq 0$. In fact, if $\theta_{i}(t)$ is a solution of (1) for some external force $F$, $\theta_{i}(t)+\pi$ will be a solution for the case corresponding to $-F$, and both cases have the same asymptotic dynamical properties. We will also assume $\Omega>0$ in our analysis, but we will see that the case $\Omega<0$ can be treated analogously.

Our analysis will be restricted to the situations where the underlying network is an undirected and unweighted large random network. We assume also that $\omega_{i}$ is a random variable over the network, with null average, symmetric distribution $g(\omega)$, and all moments finite. Notice that we could also have written the system (1) in an autonomous manner by introducing an extra node $(N+1)$ with natural frequency $\Omega$ and connected in a directional way to the forced nodes $i\in C$. For a recent discussion on the dynamics of directed and weighted networks of Kuramoto oscillators, see p14, . However, since our proposed approach is indeed more appropriate to the case of undirected and unweighted large random networks, we will deal with the original non-autonomous formulation of the problem. There are several dynamical behaviors for the system (1) which would deserve to be named synchronization. Since the external driven force has its own frequency $\Omega$, we will focus here on the network dynamical states for which the oscillators rigidly rotate in the same pace of the external force, with frequency $\Omega$, and in a phase-locked way. This situation is also called in the literature as the global synchronization for (1). Of course, it is quite natural to expect the existence of a certain threshold $F_{\rm c}>0$ for the external driven force such that, for $F<F_{\rm c}$, the system (1) would be rather insensitive to the external excitation. In fact, this is a well-known property of the system (1) for the case of a fully connected network topology, see p3, for further details.

We will employ the dimensional reduction approach recently proposed by Gottwald gottwald2015 , which in turn is based in the introduction of some collective coordinates for the Kuramoto oscillators in the same spirit of the Ott-Antonsen ansatz p2 ; OA1 ; OA2 , to investigate the global synchronization properties of (1) for the case of large random networks. The Gottwald approach was recently extended for the case of Stuart-Landau complex oscillators in AISync, to incorporate also amplitude effects in the dynamics. From the dimensionally reduced system, we will derive some conditions on the parameter space of the problem which will allow (or prevent) the onset of a globally synchronized state for the system. Moreover, in the same line of the rewiring algorithm discussed in optimal, ; AISync, , we will explore the dimensional reduction approach to propose an optimization scheme for global synchronization by selecting judiciously the forced nodes subset $C$ for a given random network. The scheme has been exhaustively tested with numerical simulations and it has always resulted in an enhancement in the global synchronization properties of the system. Our main conclusion is that the optimal subset $C$, in the sense that will be properly define in Section III, should consist in the nodes $i$ such that $|\Omega-\omega_{i}|$ is maximal, a result which is indeed in agreement with some recent related works optimal ; AISync ; SciRep , and that might shed some light on the evolution and structure of natural systems for which global synchronization is a desired property. Despite our analysis is explicitly done for the $\Omega>0$ case with symmetric distributions $g(\omega)$, we will show that the optimization criterion of maximizing $|\Omega-\omega_{i}|$ for choosing the subset $C$ of forced nodes is valid for more general cases as well.

The present paper is organized as follows. In the next section, the Gottwald dimensional reduction approachgottwald2015 is adapted for the global synchronization problem of partially forced Kuramoto oscillators governed by (1). A mean field analysis is performed, and we derive some conditions on the parameter space of the problem for the onset of a globally synchronized state. In Section III, we introduce our optimization scheme and present the results of our numerical simulations for large random networks. The last section is devoted to some concluding remarks on the implications of our results and on the role played by possible network symmetries on the global synchronization problem of forced Kuramoto oscillators, in the context of the recently introduced asymmetry-induced synchronization (AISync) scenarioNishMotter ; ZhangNishMotter ; ZhangMotter .

II The dimensional reduction approach

We will describe the global state of the system (1) by using the standard order parameters $r$ and $\psi$ defined as

[TABLE]

whose respective behaviors are well known: $r\approx N^{-1/2}$ for incoherent motion, whereas $r\approx 1$ for a fully synchronized state. Our global synchronization regime corresponds to a phase locked configuration with a rigid rotation $\Omega$, and hence to the case $r\approx 1$ and $\dot{\psi}\approx\Omega$. We will look for globally synchronized states of (1) by exploring a two dimensional model based on the Gottwald dimensional reduction approach introduced in gottwald2015, . Gottwald ansatz is based on the empirical observation that the set of synchronized Kuramoto oscillators in a complex network typically have their phases obeying $\theta_{i}\propto\omega_{i}$, up to a rigid rotation with frequency $\langle\omega\rangle$. By considering the collective ansatz $\theta_{i}(t)=\mu(t)\omega_{i}+\langle\omega\rangle$, we see that synchronization corresponds to the existence of a fixed point for $\mu(t)$ for large $t$, see optimal, for further details. In the present case, the Gottwald approach can be adapted as the collective ansatz

[TABLE]

for the entire network, which has been indeed extensively tested in our numerical simulations, as we will discuss in the next section. Notice that, by construction, both reduced dynamical variables $\alpha$ and $\beta$ are dimensionless. Using such ansatz and recalling that $\omega_{i}$ is assumed to be a random variable with null average and symmetric distribution $g(\omega)$, we will have for large random networks

[TABLE]

and

[TABLE]

from where the dynamical interpretation of $\alpha$ and $\beta$ will become rather clear. We denote the simple average of a given variable $h_{k}$ over the entire network as

[TABLE]

and the continuous (mean field) approximation consisting basically in the substitution of the sum with the integral with the pertinent probability distribution, as we have done in (4). Notice that synchronized states with $r\approx 1$ will demand $\beta\approx 0$. In fact, the expression (4) can be expanded as

[TABLE]

for any distribution $g(\omega)$ with finite moments. It is clear that asymptotic globally synchronized states will correspond to solutions such that $\beta\to 0$ and $\alpha\to$ constant for $t\to\infty$. Notice that the ansatz (3) can also accommodate the description of a synchronized regime rigidly rotating with frequency $\Omega^{\prime}\neq\Omega$. This case would correspond to a solution with $\beta\to 0$ and $\alpha\approx(\Omega^{\prime}-\Omega)t$ for large $t$.

In order to obtain the two-dimensional reduced system in the new variable $\alpha$ and $\beta$, let us substitute (3) in (1), which will result in

[TABLE]

Averaging both sides according to (6) leads to

[TABLE]

where the dimensionless function ${\cal I}_{1}(\alpha,\beta)$ is given by

[TABLE]

On the other hand, multiplying Eq. (II) by $\frac{\omega_{i}}{\Omega}$ and averaging both sides again, we get

[TABLE]

where

[TABLE]

and

[TABLE]

are also dimensionless functions. The continuous (mean field) approximation in the present case corresponds to

[TABLE]

where $f=N_{C}/N$ is the fraction of forced nodes, $\left\langle d\right\rangle$ stands for the average degree of the network and the prime denotes derivative. In the derivation of such approximations, we employ the crucial hypothesis that the $N_{C}$ forced nodes were also randomly chosen and, in particular, that their natural frequencies has the same distribution $g(\omega)$. For further details, see the Appendix A. Equations (9) and (11) in this case read simply

[TABLE]

It is important to stress that the quantities (14)-(16), which give origin to the dynamical two-dimensional system (17)-(18), are completely determined from the network and oscillators data. For instance, for normal and homogeneous distributions $g(\omega)$, the most commonly used ones in the literature, we have the following continuous approximations for $r$ from (4)

[TABLE]

and

[TABLE]

respectively, from where all the quantities (14)-(16) can be derived. The mean-field approximations are accurate for large random networks, see Fig. 1

for a typical example.

For $r(\beta)=1$, Eq. (17) reduces to the well-known Adler equation, exactly as in the full-connected topology case, see p3, and p8, for further details. Since $0<|r(\beta)|\leq 1$, equation (17) will never admit fixed points if

[TABLE]

We have just established the critical value of the external force $F_{\rm c}$ for the onset of the global synchronization, which will always require

[TABLE]

Of course, our main interest here is the stable fixed points $(\alpha_{*},\beta_{*})$ of the system (17) - (18), which Jacobian matrix for $\beta\approx 0$ reads

[TABLE]

from where we have that any stable fixed point corresponding to a synchronized state should have $fF\cos\alpha_{*}>0$. One can explore the right-handed side of equation (18) to derive some critical value for $\lambda$ in order to assure a globally synchronized state in the same way we did for $F$ in (22), but it turns out that the situation is very similar to the case without the external force ($F=0$), which was treated in detail in some previous works optimal ; AISync . In particular, one should not expect any stable fixed point for (17) - (18) if

[TABLE]

where $\zeta$ is the constant

[TABLE]

One can evaluate the constant $\zeta$ for the normal and the homogeneous distributions $g(\omega)$. Taking into account (19) and (20), we have $\zeta={\sqrt{e}}\approx 1.649$ for the uniform distribution, while for the uniformly distributed case we can numerically determine that $\zeta\approx 2.293$. However, in contrast with the critical force prediction (22), which has proved to be indeed quite accurate in our numerical simulations, the condition (24) for a possible $\lambda_{\rm c}$ seems to be a rather conservative one. Typically, one will need larger values for $\lambda$ in order to have robust synchronized states. This situation is completely analogous to the $F=0$ case discussed previously in optimal, ; AISync, .

III The optimization scheme

Our main purpose here is to present a prescription to select judiciously the subset $C$ of forced nodes in order to assure a better synchronization capability for the network. We will consider that a network has a better capability for a globally synchronized state if one can attain states with $r\approx 1$ and $\dot{\psi}\approx\Omega$ with smaller values for the external force intensity $F$. Since the globally synchronized states require $\beta\approx 0$, let us linearize (9) and (11) around $\beta=0$ and abandon the hypothesis that $C$ is a random subset of nodes of the network. We will have in this case

[TABLE]

where $\langle\ \ \rangle_{C}$ denotes the average in the subset $C$ and

[TABLE]

where the sum in the numerator is performed over the edges $e(i,j)$ of the network. The quadratic quantity $\mathcal{L}$ is known to play a crucial role in the usual ($F=0$) synchronization problem optimal ; AISync . In particular, we have that the larger the value of $\mathcal{L}$, the better the synchronization properties of the underlying network, see AISync, for further details. The increasing of $\cal L$ by means of some rewiring operations in the network is the central point of the optimization algorithm introduced in optimal, , which consists basically in changing the edges $e(i,j)$ of the network in order to connect oscillators such that $|\omega_{i}-\omega_{j}|$ is maximal.

The fixed points of (26) and (III) are such that

[TABLE]

from where we can see that one effectively attains better values for $r$ (close to 1, requiring smaller $\beta_{*}$) if $\langle\omega^{2}\rangle_{C}$ is maximal and $\langle\omega\rangle_{C}$ is minimal. This is equivalent to select the subset $C$ of forced nodes as those ones with the minimal values of their frequencies $\omega_{i}$ in the network. Since we are assuming $\Omega>0$, one can say that the optimal subset $C$ is formed by the oscillator such that $|\omega_{i}-\Omega|$ is maximal. As we will see, this criterion is also valid for the $\Omega<0$ case and even for more general distributions $g(\omega)$. Notice that, from (29), we have that $\lambda\cal L$ also plays an important role here. As in the $F=0$ case, the larger the value of $\lambda\mathcal{L}$, the larger the value of $r$ (smaller $\beta_{*}$). On the other hand, the system is rather insensitive to the precise value of $\alpha_{*}$. Our predictions obtained from the dimensionally reduced model, including the optimization criterion, were exhaustively tested in numerical simulations, whose main details will be presented in the following subsection.

III.1 Numerical results

For our simulations, we have numerically solved the equations (1) for large random networks and construct several synchronization diagrams. We have made extensive use of the 3.7 Python packages NetworkXNetworkX , which allows us to build many types of random networks with some prescribed topological and statistical properties, and SciPySciPy , particularly its integrate.solve_ivp function. For our purposes here, it is more convenient to introduce the dimensionless evolution parameter $\tau=\Omega t$ in (1), leading to the following system of ordinary differential equations

[TABLE]

where the prime denotes differentiation with respect to $\tau$. As we can see, for a given random network with adjacency matrix $A_{ij}$ and oscillator frequencies $\omega_{i}$, our parameter space is effectively bi-dimensional and spanned by $(F/\Omega,\lambda/\Omega)$. In this case, $\Omega^{-1}$ simply plays the role of the unit of time for the problem. The equations (30) are the base of all our numerical analysis. We perform the simulations for different network topologies and frequencies distributions $g(\omega)$ with null average, and we have not detected any appreciable dependence of our results on the network topology or the employed distribution $g(\omega)$. All results presented here correspond to the case of an Erdos-Renyi random network with 1500 nodes, average degree $\langle d\rangle\approx 7.7$ (the maximal node degree is 18, the minimal is 1), and with 151 nodes subjected to the external force. We employ a normal distribution $g(\omega)$ with $\sigma_{\omega}/\Omega=1$ for the oscillators natural frequencies. None of these parameters has demonstrated to have strong influence on the results. The initial conditions $\theta_{i}(0)$ for the numerical solutions of (30) were randomly chosen, with uniform distribution, in the circle $[0,2\pi)$.

Our first results, depicted in Fig. 2,

correspond to the verification that the ansatz (3) is indeed valid for the study of global synchronization in the system (30).

The figure shows the $(\omega_{i},\theta_{i}(\tau)\mod 2\pi)$ graphics, for a fixed $\tau$ in a globally synchronized regime with $\lambda/\Omega=2$ and $fF/\Omega=3$. Each point in the graphics corresponds to a Kuramoto oscillator in the network. The depicted line is the ordinary linear regression of the data. The linear correlation between $\omega_{i}$ and $\theta_{i}$, as incorporated in the ansatz (3), is evident. Moreover, we can clearly identify two displaced oscillator populations with similar linear regression slopes. We will return to this interesting dynamical behavior in the last section.

In order to test our optimization scheme, we have considered several synchronization diagrams of the type $r(\tau)$ and $\psi(\tau)$ versus $\lambda/\Omega$ and $F/\Omega$, for $\tau$ sufficiently large to assure the relaxation of any transient regime. The order parameters $r(\tau)$ and $\psi(\tau)$ are evaluated accordingly their definition (2) for the numerical solutions of (30), taking into account the usual statistical precautions, see PintoSaa, ; optimal, for further details. Our results are summarized in Fig. 3, from where one can appreciate that the selection of the subset of forced nodes according to our optimization scheme effectively gives origin to networks with better global synchronization capabilities. The diagrams $r\times fF/\Omega$ and $r\times\lambda/\Omega$ have evident meaning and can be easily understood. The case for the order parameter $\psi$ is more involved. Its asymptotic behavior for large values of $\tau$ in a synchronized regime is expected to be $\psi(\tau)\sim a\tau.$ For a globally synchronized state, we will have $a=1$, or $\psi\sim\Omega t$, in terms of the original time variable $t$ of (1). However, different values for $a$ are also associated with synchronized states, but for which the oscillators do not follow the same pace of the external force. In particular, the usual synchronization state for the Kuramoto network with $F=0$ is known to be such that $a=\langle\omega\rangle$. We have explicitly compared three cases for the subset $C$, namely the optimal, the random, and the worst cases. The optimal and the worst cases correspond, respectively, to selection of the forced nodes as those ones with maximal and minimal values of $|\omega_{i}-\Omega|$. The results are in complete agreement with the expectations for our optimization scheme: the optimal subset always implies better global synchronization capabilities than the random one, while the worst subset always exhibits worse capabilities when compared with the random case.

IV Final Remarks

Although all our analyses have effectively been done for the case with $\Omega>0$ and for symmetric distributions $g(\omega)$ with null average, the optimization criterion of selecting the subset $C$ as those oscillators with maximal value of $|\omega_{i}-\Omega|$ is also valid for the general case. Let us consider, first, the case with $\Omega<0$. From (1), supposing that $\theta_{i}(t)$ is a solution for the case with $\Omega>0$ for a network of Kuramoto oscillators with natural frequencies $\omega_{i}$, we have that $-\theta_{i}(t)$ will be the solution for the case corresponding to $\Omega<0$, but for a network with natural frequencies $-\omega_{i}$. Both situations will have the same fixed points for the equations equivalent to (26) and (III), but now the criterion of minimal $\langle\omega\rangle_{C}$ for $\Omega<0$ will effectively select those nodes with larger natural frequencies, which indeed corresponds to maximize $|\omega_{i}-\Omega|$. The case of non-symmetric distributions is a bit more involved since one needs to keep track of the terms proportional to $\langle\omega\rangle$ in the derivation of (9) and (11) and all the subsequent manipulations. For this case, Eq. (9) and (11) will read

[TABLE]

where

[TABLE]

Since $\det M=\frac{\sigma^{2}_{\omega}}{\Omega^{2}}$, we see that, unless all $\omega_{i}$ be equal, a situation where our optimization procedure obviously does not apply, the fixed points $(\alpha_{*},\beta_{*})$ in this case correspond to the zeros of the right-handed side of (31). However, the linearization of the right-handed side of (31) around $\beta=0$ gives origin to exactly the same condition (29), and hence the same analysis of the symmetric case does apply here. Hence, our optimization criterion is indeed valid for the case of more general $\Omega$ and $g(\omega)$.

The role played by possible symmetries of a given network in its synchronization capability was recently discussed in AISync, , for the case of multilayer networks of Stuart-Landau complex oscillators, in the context of the so-called asymmetry-induced synchronization (AISync) scenarioNishMotter ; ZhangNishMotter ; ZhangMotter . The main conclusion was that the presence of certain regularities in the interlayer connection pattern tends to diminish the synchronization capability of the network or, in other words, asymmetries in the network tend to enhance its synchronization properties. The key point here is the quadratic quantity $\cal L$ given by (28), which always decreases if interlayer symmetries are present, see AISync, for the details. It is important to stress that the same conclusions hold here: the presence of possible symmetries as those ones discussed in AISync, tends to diminish the global synchronization capabilities of the forced network. It is also interesting to notice that, in the autonomous formulation of the problem (1), in which an extra extra $i=N+1$ node with natural frequency $\Omega$ is connected in a directional way to the forced nodes $i\in C$, one can interpret our optimization scheme in the same way of the optimal synchronization for the $F=0$ caseoptimal , in particular that anti-correlation between the frequencies of neighbor nodes always favors synchronization, see also skardal2014, . The investigation of such phenomena might shed some light on the evolution and structure of natural systems for which global synchronization is a desired property. Many recent studies have been devoted for theses issues as, for instance, the study of the response to external stimuli of the C. elegans neural networkmodular ; p8 , self-organization and pattern formation in the growing and developing of vertebratesp9 , auditory signals in amphibiansp10 , and several topics on the circadian rhythmp11 ; p12 ; p13 .

We finish recalling the dynamical behavior involving the two oscillator populations we could identify in the globally synchronized regime, see Fig. 2. Typically, the displaced oscillators are the forced ones, giving origin to a curious dynamical configuration. When the global synchronization is completely attained, the forced nodes clump together and move ahead of the remaining swarm of oscillators. The dynamics in this case consist in a large group of oscillators (the free ones) chasing the smaller group of the forced ones. The existence of the two disjoint populations prevents the globally synchronized state to have an order parameter $r$ very close to 1, which can be clearly seen in the top left panel in Fig. 3, where one can observe a sudden decrease in $r$ when global synchronization is attained. It is certainly worth to incorporate such two populations behavior in a second order approximation for the collective ansatz (3). These points are now under investigation.

Appendix A Mean-field approximations for the random network

The mean-field expression (14) is obtained directly from the evaluation of (10). Notice that

[TABLE]

where $f=N_{C}/N$ is the fraction of forced nodes and the sum is performed on the forced nodes. With the hypothesis that $C$ is a random subset of the oscillators and that the natural frequencies of the elements in $C$ have the same distribution $g(\omega)$, one can write

[TABLE]

and then (14) follows from the definition (4) and that $g(\omega)$ is symmetric and has null average. The same hypothesis of a random subset $C$ allows us to write (12) as

[TABLE]

where the condition of symmetric $g(\omega)$ with null average was used again, and hence (15) follows directly from (4) as well. The evaluation of (16) a little bit more involved. First, notice that

[TABLE]

where

[TABLE]

with $d_{i}=\sum_{j}A_{ij}$ standing for the degree of the $i^{\rm th}$-node, which is a positive integer since our networks are always assumed to be connected. The sum in (37) is performed over the set $D_{i}$ of nodes connected to the $i^{\rm th}$-node in the network. Since we are dealing with random networks with connection pattern independent of the nodes natural frequencies, one has

[TABLE]

where we have used again that $g(\omega)$ is symmetric and has null average. Inserting (38) in (36) leads to

[TABLE]

where $\langle d\rangle$ is the average degree of the network and we have used the condition that $d_{i}$ and $\omega_{i}$ are independent random variables. Eq. (16) follows now directly from the definition of $r(\beta)$ given by (4).

Acknowledgment

The authors acknowledge the financial support of CNPq and FAPESP (Grant 2013/09357-9). They also wish to thank M.A.M de Aguiar and C.A. Moreira for enlightening discussions and the anonymous referees for valuable suggestions.