# Small-world networks of Kuramoto oscillators

## Abstract

The Kuramoto model of coupled phase oscillators on small-world (SW) graphs is analyzed in this work. When the number of oscillators in the network goes to infinity, the model acquires a family of steady state solutions of degree , called -twisted states. We show that this class of solutions plays an important role in the formation of spatial patterns in the Kuramoto model on SW graphs. In particular, the analysis of -twisted states elucidates the role of long-range random connections in shaping the attractors in this model.

We develop two complementary approaches for studying -twisted states in the coupled oscillator model on SW graphs: linear stability analysis and numerical continuation. The former approach shows that long-range random connections in the SW graphs promote synchronization and yields the estimate of the synchronization rate as a function of the SW randomization parameter. The continuation shows that the increase of the long-range connections results in patterns consisting of one or several plateaus separated by sharp interfaces.

These results elucidate the pattern formation mechanisms in nonlocally coupled dynamical systems on random graphs.

## 1 Introduction

Coupled dynamical systems on graphs model the evolution of complex networks of diverse nature. Examples include neuronal networks in biology [6]; Josephson junctions and coupled lasers in physics [14, 24]; communication, sensory, and power networks in technology [8], to name a few. In these models, each node of the graph supports a local dynamical system, whose dynamics is typically well understood. Local systems at the adjacent nodes interact with each other, thus, producing collective dynamics of the coupled system. The spatial structure of the network may vary from regular as in lattice dynamical systems to completely random. Thus, it is not surprising that dynamical systems on graphs exhibit a wide range of spatio-temporal dynamics.

The main challenge in studies of dynamical networks has been relating the structure of the graph of the network to its dynamics. In view of the enormous complexity of this class of systems, our understanding of the link between the structure and dynamics has been guided by the analysis of certain successful models illuminating various aspects of network dynamics. In this respect the Kuramoto model of coupled oscillators stands out as one of the most influential in the field [13]. Derived as the phase description of large ensembles of limit cycle oscillators, the Kuramoto model has been used by physicists and biologists to gain insights into many important effects in large networks such as synchronization, multistability, phase transitions, wave propagation, and spatio-temporal chaos, to name a few [11]. Recent experiments with the Kuramoto model with nonlocal coupling led to the discovery of chimera states, spectacular patterns combining coherent and irregular turbulence-like behavior [12].

The Kuramoto model describes dynamics of coupled phase oscillators on a graph. Thus, before we introduce the model, we need to review some graph-theoretic notation. Let be an undirected graph, where stands for the set of nodes of , and the edge set contains the unordered pairs of connected nodes. Furthermore, is called a simple graph, if it does not contain loops or multiple edges.

Let be a sequence of simple graphs. If then graphs are called dense.

The following graph sequence will appear frequently throughout this paper. For sufficiently large and such that , we define a sequence of -nearest-neighbor (-NN) graphs as follows

where . Let then is a dense graph sequence.

We now introduce the Kuramoto model of coupled phase oscillators on

(1.1) |

where models the evolution of the phase oscillator , and is its intrinsic frequency. The scaled sum on the right hand side of (1.2) describes the interactions between oscillator and the oscillators connected with it. In the classical Kuramoto model, . This is so-called attractive coupling case, because when , synchrony is stable. Recently, there has been an interest in studying the Kuramoto model with , i.e., the repulsive coupling case [10]. Existence of steady state solutions clearly does not depend on the sign of . However, the stability of these solutions does. In this paper, we are interested in the case of dense graphs , which corresponds to the nonlocal coupling in the Kuramoto model (1.1).

If all intrinsic frequencies are equal, as we assume for the remainder of this paper, they can be eliminated by recasting (1.1) in the rotating frame of coordinates. Thus, we consider

(1.2) |

Our work was inspired by the analysis of the Kuramoto model on the sequence of -NN graphs in [26]. In this paper, Wiley, Strogatz, and Girvan identified an important class of steady state solutions of degree ,

(1.3) |

These solutions are called -twisted states. The -twisted states correspond to spatially uniform (synchronous) solutions. Steady states with form linear patterns that wind around times (see Fig. 1). It was shown in [26] that the stability of -twisted states for the Kuramoto model depends on the range of nonlocal coupling . This provides a nice case study illustrating how the network topology is shaping the structure of the phase space of the coupled system. The follow-up study of -twisted states in the Kuramoto model with repulsive coupling suggests their relevance to the chimera states and, more generally, to the coherence-incoherence transition in coupled systems [10].

Interestingly, the original motivation for the analysis in [26] was the dynamics on SW networks. SW graphs, proposed by Watts and Strogatz, interpolate between graphs with regular local connections and completely random Erdős-Rényi graphs [25]. In the classical Watts-Strogatz SW graph, one starts with the -NN graph on a circle (). Then each edge in is replaced with probability to a randomly chosen one (see Fig. 2). We will refer to as the randomization parameter. If then is transformed into the Erdős-Rényi graph. For analytical convenience Newman and Watts and independently Monasson used another variant of SW connectivity, where the new random edges are added to the graph, while no edges are removed [20, 19].

Like many other random graphs, SW graphs have small average path length, but they also have large clustering coefficients inherited from the regular graphs. The combination of these two properties has been shown in many real-life networks, including neuronal networks [25, 3], the power grid of the western United States, and professional networks [25]. Thus, studying dynamical systems on SW graphs has the potential of elucidating principles governing dynamics in diverse realistic networks. A specific question posed in [26] is how adding a small number of random long range connections to otherwise regular network affects stability of the synchronous state. The synchronizing effect of long-range connections was suggested by the previous studies of related models [21, 2]. However, the rigorous understanding of this effect in SW networks was lacking.

In this paper, we study the role of the long-range random connections in shaping stationary solutions of the Kuramoto model on SW graphs. For this, we use two complementary approaches: the linear stability analysis of -twisted states and numerical continuation. Neither of this approaches can be applied to the problem as it stands, and, therefore, some preliminary work is required. We comment on the stability analysis first. In [26], Wiley, Strogatz, and Girvan obtained a sufficient condition for stability of -twisted states for the -NN coupled Kuramoto model. Their analysis employed the continuum limit of the Kuramoto model. Unlike the regular -NN coupling case, on SW graph when the number of oscillators is finite, -twisted states are no longer solutions of the Kuramoto equation. We show that they become steady state solutions again in the appropriate probabilistic sense in the limit as . Thus, as before the continuum limit provides a natural setting for the stability analysis. Our previous work [18] shows how to derive continuum limits for the Kuramoto model on W-random graphs including the SW graphs. This allows the stability analysis of the -twisted states for the SW model to be performed in analogy to the analysis of the deterministic -NN networks in [26].

The stability analysis of -twisted states in the continuous setting clearly shows the respective contributions of the regular local and random long-range connections of the SW network topology to the spectrum of the linearized problem. It shows that adding new long-range random connections enhances stability of the synchronous regime, while inhibiting other -twisted states. Furthermore, the estimate of the leading eigenvalue of the linearized problem yields the synchronization rate as a function of the randomization parameter . The synchronization rate grows monotonically with . Thus, our quantitative estimates of the eigenvalues of the linearized problem show that adding new random connections in the SW topology strengthens the synchrony.

As a complementary approach aimed at elucidating the impact of long-range connections on the solutions of the Kuramoto model and on the spatial patterns that it forms, we developed a continuation method that reveals the transformations of the -twisted states under the variation of the randomization parameter . The major obstacle in implementing the continuation procedure was the fact that in general the right hand side of the Kuramoto model is not continuous in . This is easy to understand by noting that the adjacency matrices of two independent copies of SW graphs, and , do not have to be close when is small or even when it is equal to . To overcome this problem we constructed a special family of SW graphs , whose members depend continuously on . Here, the continuity is meant in the appropriate probabilistic sense (see Section 5 for details). We call such family of SW graphs consistent. With the consistent family of SW graphs at hand, we set up a numerical continuation algorithm based on the Newton’s method and computed branches of steady state solutions of the Kuramoto model on SW graphs for by starting from -twisted states for . The numerical continuation revealed an unexpected effect. It turns out that for increasing values of , the -twisted states are transformed into patterns composed of plateaus separated by interfaces. We observed this effect in all random realizations of the consistent families of SW graphs that we used in our experiments. The formation of plateaus in solutions continued from -twisted states is another manifestation of the synchronizing effect of the long-range connections.

This paper is organized as follows. In Section 2, we redefine SW graphs using the framework of W-random graphs [18] and derive the continuum limit of the Kuramoto model on SW graphs. In Section 3, we show that in the limit as , -twisted states become stationary solutions of the Kuramoto model on SW graphs. This is followed by the linear stability analysis of -twisted states. The implications of the stability analysis for synchronization are discussed in Section 4. The numerical continuation from the -twisted states is developed in Section 5. The final Section 6 contains concluding remarks.

## 2 The continuum limit of the Kuramoto model

The stability analysis of -twisted states in Section 3 uses the continuum (thermodynamic) limit of the Kuramoto equation (1.2) as . The thermodynamic limit has proved to be a useful tool for elucidating various aspects of the dynamics of large networks (see, e.g., [12, 1, 26, 10, 22]). In [17], we provided a rigorous justification of the continuum limit for the models of Kuramoto type on dense graphs, using the ideas from the theory of graph limits [15, 16, 5]. This work was further extended to cover dynamical systems on random graphs in [18]. In particular, the analysis in [18] shows how to set up and to justify the continuum limit for the Kuramoto model on the SW graphs. Before we apply the continuum limit to study the stability of -twisted states, we review a few facts illuminating the basis for the limiting dynamics of (1.2) as .

As the first step toward the continuum limit, we redefine SW graphs as W-random graphs [15, 16]. The latter provide a flexible framework for constructing random graphs, which is especially convenient for studying coupled dynamical systems. Below, we review the construction of W-random graphs and use it to define SW graphs (cf. [18]).

Let be a class of symmetric measurable functions on with values in . Here and throughout this paper, stands for the unit interval .

###### Definition 2.1.

[16] Let and

be a set of distinct points from . Define such that for each

where stands for the probability of an event. Graph is called a -random graph and denoted by .

In the remainder of this paper, we use

(2.1) |

and

(2.2) |

where is a fixed parameter, and

(2.3) |

With these choices of and , is the -NN graph, which is used as a starting point in the construction of SW graphs [25, 20].

We are now in a position to interpret SW networks as W-random graphs (cf. [18]). To this end, let

(2.4) |

Let us call a local edge if , and a long-range one otherwise. The construction of can be interpreted as follows: is obtained from by letting the latter to shed each local connection with probability and to add each long-range one with the same probability (see Fig. 2b). In particular, coincides with , and is the Erdős-Rényi graph with the edge density . The construction of resembles that of the Watts-Strogats SW graphs [25], but differs from it in some details. For other -random graphs, which match the classical Watts-Strogatz and Manasson-Newman-Watts models of SW graphs [25, 20, 19] more closely, we refer an interested reader to [18, Section 5]. In the remainder of this paper, (2.4) will be used as a model of SW connectivity.

Now we turn to the main point of this section: the derivation of the continuum limit for the Kuramoto model on the sequence of SW graphs . To explain the intuition behind the limiting procedure, let us rewrite (1.2) as

(2.5) |

where stands for the adjacency matrix of . Our goal is to derive the infinite-dimensional analog of (2.5) as . To this end, we define a step-function such that

(2.6) |

where is a solution of the discrete model (2.5). We expect that as , the sequence of step-functions approaches (in the appropriate norm). The latter function is interpreted as the solution for the continuum model.

Next, we have to decide how to interpret the adjacency matrix as . The answer to this question is suggested by the theory of graph limits, which for every sequence of convergent dense graphs, like , defines the limiting object - a symmetric measurable function on the unit square. The limit of a convergent graph sequence is called a graphon [15]. From the theory of graph limits, we know that , as a sequence of -random graphs, converges to the underlying graphon (cf. [16]). Below, we provide an informal intuitive interpretation of the graph limit , which is meant to explain its appearance in the continuum limit of (2.5). A detailed exposition of the theory of graph limits can be found in the monograph by Lovász [15], or in the following papers [5, 16, 23].

To see the meaning of the graph limit , it is instructive to recall the geometric representation of the adjacency matrices of SW graphs in Fig. 2b,c. The matrices depicted in these plots correspond to and . These graphs belong to the same sequence of SW graphs generated by graphon with (see (2.4)). Note that for large , as in Fig. 2c, black pixels fill different regions of the unit square uniformly with two different densities and , just as in the definition of the function (see (2.2) and (2.4)). In fact, the graphon represents the density of connections between different subsets of nodes of the SW graph for large .

We have identified and as continuous counterparts of and in the large limit. Furthermore, the right-hand side in (2.5) has the form of a Riemann sum. Therefore, by sending in (2.5), we formally arrive at the following evolution equation

(2.7) |

The continuum limit (2.7) was rigorously justified in [18]. Specifically, Theorem 3.4 of [18] implies that for any , (see (2.6)) converge to in the norm in probability as , provided the initial conditions and are sufficiently close.

## 3 The stability analysis of q-twisted states

The goal of this section is to elucidate the role of -twisted states (1.3) in the long-time asymptotic behavior of solutions of the Kuramoto model (1.2) on SW graphs (2.4) for large . First, we show that -twisted states become steady state solutions of (1.2), (2.4) as with probability (see Lemma 3.1). Next, we use the continuum limit (2.7) to study the linear stability of the continuous -twisted states

(3.1) |

(see Theorem 3.3). The results of this section will be used in Section 4 to study synchronization in SW networks.

Note that while the continuous -twisted states (3.1) solve (2.7) for all values of , their discrete counterparts generically do not solve (1.2) for . They nonetheless remain relevant for the following reasons. First, for large , (1.2) has solutions that remain close to on finite time intervals (cf. [18, Theorem 3.4]). Furthermore, the following lemma shows directly (i.e., without invoking the limiting equation (2.7)) that become the steady-state solutions of (1.2) for in the limit as .

###### Lemma 3.1.

Proof. Assume first that . Let be arbitrary but fixed and let . To simplify notation, throughout this proof we suppress the dependence of various sets and variables on . Denote

and Note that for , tend to infinity with

Let be two arrays of independent identically distributed (IID) random variables (RVs) such that

(3.4) |

Here, stands for the probability distribution of a RV, and denotes the binomial distribution with parameter . Define

(3.5) |

From this point, the proof follows the lines of the proof of the strong law of large numbers [4, Theorem 6.1]. Specifically, we estimate

(3.11) | |||||

By Markov inequality, using (3.10) and (3.11), for any we have

(3.12) |

which in turn implies via the first Borel-Cantelli lemma [4, Theorem 4.3]

By Theorem 5.2(i) in [4], the last statement is equivalent to convergence of to almost surely. Therefore, the combination of (3.6),

yields (3.3) for .
If then , and (3.3) follows by repeating the above argument
for only.

Having shown, that -twisted states solve the Kuramoto equation on in the large limit, we turn to study the stability of -twisted states, as steady state solutions of the continuum equation. For concreteness, in the remainder of this paper we restrict to the attractive coupling case, i.e., in (1.2). The repulsive coupling case can be treated similarly. Thus, we rewrite (2.7) as follows

(3.13) |

where

(3.14) |

and was defined in (2.3).

###### Remark 3.2.

The results of this section hold for a more general class of kernels . It is sufficient to assume that is a even bounded function.

We are now in a position to formulate a sufficient condition for stability of -twisted states (see (3.1)).

###### Theorem 3.3.

For and , the -twisted state is a steady state of (3.13). Moreover, it is linearly stable with respect to perturbations from provided

(3.15) |

where

(3.16) |

###### Remark 3.4.

We precede the proof of Theorem 3.3 with the following observation, which follows from the odd-symmetry of the product .

###### Lemma 3.5.

Let be a solution of (3.13). Then

(3.17) |

Proof. By integrating both parts of (3.13) over and using Fubuni’s theorem, we have

(3.18) |

where we also interchanged differentiation with integration on the left-hand side,
because and is bounded on
(cf. [17, Theorem 3.1]). The integral on the right-hand side
is zero, because the domain of integration is symmetric
with respect to the line and the integrand is an odd symmetric function.

Proof. (Theorem 3.3) Consider as a subspace of . Then

where

By Lemma 3.5, and are invariant under the flow (2.7). It is easy to see that -twisted state solutions of (2.7) remain stable under spatially homogeneous perturbations from . Therefore, it is sufficient to study the stability of -twisted states with respect to the perturbations from .

To this end, we linearize (3.13) about the -twisted state , i.e., we plug into (3.13) and retain only the linear terms on the right-hand side of the resultant equation

(3.19) |

which we supply with the initial condition

(3.20) |

To simplify notation we rescale the time variable in (3.19) to absorb the factor on the right-hand side. Thus, (3.19) can be rewritten in the following form

(3.22) |

where denotes the convolution and

By applying the Fourier transform to both sides of (3.22), we have

(3.23) |

Thus, using Parseval’s identity, we arrive at the sufficient condition for linear stability of the -twisted state

(3.24) |

Further, it is easy to verify that

Thus, (3.24) can be rewritten as follows

(3.25) |

where

Since is an even function, it is sufficient to restrict (3.25)
for .
This concludes the proof of Theorem 3.3.

## 4 Synchronization

In this section, we apply Theorem 3.3 to study synchronization in the Kuramoto model on SW graphs (3.13). Throughout this section, we continue to use the rescaled linearized equation (3.22).

By applying Theorem 3.3 to (3.13) with (3.14), we obtain the following condition for stability of the -twisted states

(4.2) |

where stands for the second-order difference

(4.3) |

Using (4.1), we rewrite (4.2) as follows

(4.4) |

where

and the Kronecker delta if and is equal to otherwise.

Equation (4.4) yields the following sufficient condition for synchronization in SW networks

(4.5) |

where

(4.6) |

Conditions (4.4) and (4.5) show how the random long-range connections in the SW network affect the stability of the -twisted states. For the -NN coupled networks, the linear stability of the -twisted states was studied in detail by Wiley, Strogatz, and Girvan [26]. They found that the -twisted state is stable provided , where is a solution of a certain transcendental equation (see [26, Equation (18)]). For these values of , we have

Thus, the synchronous solution () is stable for the full range of . In addition, for small , there may be many other stable -twisted states. For such values of , these -twisted states remain stable on SW graphs with

(4.7) |

From Equation (4.5) we find that the synchronous state remains stable for all . Moreover, for , Equation (4.5) provides a relaxed condition compared to the corresponding conditions for . In contrast, from Equation (4.4) we find that for increasing values of the conditions for stability become more stringent. In particular, (4.7) shows that a -twisted state, which is stable for , may lose stability at some positive value of . From this, we conclude that long-range random connections promote synchrony, while inhibiting other -twisted states.

Next, we discuss the stability of the synchronous regime in SW networks in a little more detail. To this end, we employ the leading exponent

(4.8) |

which reflects the rate of convergence to the synchronous state. Larger values of indicate stronger synchronization. For small values of the coupling range , from (4.5), we estimate

(4.9) |

Equation (4.9) shows that larger values of result in faster synchronization. Furthermore, from (4.5) and (4.9), we estimate the leading exponent for SW connectivity

(4.10) |

Since for small , , we have

i.e., for increasing values of the leading exponent grows. This once again shows that adding the long-range random connections helps synchronization in the Kuramoto model on SW graphs.

## 5 Continuation

In this section, we develop a numerical continuation method, which will be used to study the impact of the long-range random connections on the spatial structure of the stationary solutions of the Kuramoto model (1.2).

### 5.1 Preliminaries

First, we adjust our notation to the purpose of this section. To this end, we rewrite (1.2) as follows

(5.1) |

where and . Matrix is an adjacency matrix of the SW graph with the randomization parameter . Below we will use the -norm of matrix

The vector form of (5.1) is given by

(5.2) |

where . Recall that for any , -twisted state

is a stationary solution of (5.2) with , i.e.,

(5.3) |

Thus, one can try to continue from for near to see how the addition of new random edges affects the structure of the steady state solutions of (5.2). However, one immediately runs into the following problem: in general, is not (stochastically) continuous in , as the following lemma shows.

###### Lemma 5.1.

Let and be two independent copies of SW graphs with randomization parameters . Then

(5.4) |

where

and and stand for the adjacency matrices of and respectively.

Proof. Let

be the set of indices corresponding to the edges of -NN graph, and let be its complement.

Denote RVs

We show that are IID binomial RVs. Indeed, let and . Then

Likewise, for and we have

Thus, and

###### Remark 5.2.

Note that for , stays bounded away from as .

### 5.2 Consistent family of SW graphs

Lemma 5.1 shows that is not stochastically continuous in . Therefore, nonlinear equation

(5.7) |

as it stands can not be used for continuation. To rectify this problem, we introduce a special family of SW graphs, along which is stochastically continuous in .

###### Definition 5.4.

Let and . Consider a family of SW graphs , where . Graphs are said to form a consistent family of SW graphs if