A global synchronization theorem for oscillators on a random graph

Networks of coupled oscillators have long been studied in biology, physics, engineering, and nonlinear dynamics.^1–13 Recently, they have begun to attract the attention of other communities as well. For example, oscillator networks have been recognized as having the potential to yield “beyond-Moore’s law” computational devices¹⁴ for graph coloring,¹⁵ image segmentation,¹⁶ and approximate maximum graph cuts.¹⁷ They have also become a model problem for understanding the global convergence of gradient descent in nonlinear optimization.¹⁸ In such settings, global issues come to the fore. When performing gradient descent, for instance, one typically wants to avoid getting stuck in local minima. Conditions to enforce convergence to the global minimum then become desirable. Likewise, conditions to enforce the global synchrony of oscillators play an analogous role in dynamical systems theory.

We say that a network of oscillators globally synchronizes if it converges to a state for which all the oscillators are in phase, starting from all initial conditions except a set of measure zero. Until recently, only a few global synchronization results were known for networks of oscillators.^4,5 These results were restricted to complete graphs, in which each oscillator is coupled to all the others. In the past decade, however, several advances have been made for a wider class of network structures, starting with work by Taylor,¹⁹ who proved that if each oscillator in a network of identical Kuramoto oscillators is connected to at least 94% of the others, the network will fall into perfect synchrony for almost all initial conditions, no matter what the topology of the network is like in other respects. Taylor’s result was strengthened by Ling et al.¹⁸ in 2018, and further progress has been made since then.^20,21

Ling et al.¹⁸ also made a seminal advance in the study of random networks. They considered identical Kuramoto oscillators on an Erdős–Rényi random graph,²² in which any two oscillators are coupled with unit strength, independently and at random, with probability $0≤p≤1$⁠; otherwise, the oscillators are uncoupled. They showed that if $p≫log⁡(n)/n1/3$⁠, then with high probability, the network is globally synchronizing as $n→∞$⁠.

The open question is to find and prove the sharpest result along these lines. Intuitively, as one increases the value of $p$ from $0$ to $1$⁠, one expects to find a critical threshold above which global synchrony is extremely likely and below which it is extremely rare.

At the very least, for global synchrony to be ubiquitous, $p$ must be large enough to ensure that the random network is connected, and from the theory of random graphs,²² we know that connectedness occurs with high probability once $p>(1+ε)log⁡(n)/n$ for any $ε>0$⁠. Therefore, the critical threshold for global synchronization cannot be any smaller than this connectivity threshold and is apt to be a little above. On the basis of numerical evidence, Ling et al.¹⁸ conjectured that if $p≫log⁡(n)/n$⁠, then Erdős–Rényi graphs are globally synchronizing as $n→∞$⁠, but nobody has proven that yet. The challenge is to see how close one can get. In this paper, we come within a factor of $log⁡(n)$ and prove that $p≫log2⁡(n)/n$ is sufficient to give global synchrony with high probability.

We study the homogeneous Kuramoto model^19,23,24 in which each oscillator has the same frequency $ω$⁠. By going into a rotating frame at this frequency, we can set $ω=0$ without loss of generality. Then, phase-locked states in the original frame correspond to equilibrium states in the rotating frame. Therefore, to explore the question that concerns us, it suffices to study the following simplified system of identical Kuramoto oscillators:

where $θj(t)$ is the phase of oscillator $j$ (in the rotating frame) and the adjacency matrix $A$ is randomly generated. In particular, $Ajk=Akj=1$ with probability $p$⁠, with $Ajk=Akj=0$ otherwise, independently for $1≤j,k≤n$⁠. Thus, all interactions are assumed to be symmetric, equally attractive, and of unit strength. (The assumption that all interactions are equally attractive is made for convenience. We expect our arguments and results below can be generalized to attractive interactions of non-uniform strength.) We take the unusual convention that the network can have self-loops so that oscillator $i$ is connected to itself with probability $p$⁠; i.e., $P[Ajj=1]=p$⁠. Since $sin⁡(0)=0$⁠, this convention does not alter the dynamics of the oscillators, but it does make proof details easier to write down.

Finally, since the adjacency matrix $A$ is symmetric, we know that (1) is a gradient system.^23,24 Thus, all the attractors of (1) are equilibrium points. Therefore, we do not need to concern ourselves with the possibility of more complicated attractors, such as limit cycles, tori, chimera states, or strange attractors. However, one slight complication is that the equilibria are not isolated; each equilibrium point lies on a subspace of equilibria because of the rotational symmetry of the Kuramoto model. Specifically, if $θ=(θ1,…,θn)$ is an equilibrium, then so is $θ^=(θ1+a,…,θn+a)$ for any real number $a$⁠. Despite this non-generic feature, the long-term dynamics are simple. Since the dynamics of the oscillators is governed by a gradient system and the gradient is zero along these subspaces of equilibria, each solution of (1) will still converge to an equilibrium point (that happens to lie on a subspace of equilibria). In fact, since the time derivative of $∑j=1nθj(t)$ is zero, the average of the $θj$’s remains constant over time with the constant being determined by the initial conditions. This means that one can assume that the dynamics of the oscillators is restricted to a hyperplane on which there is a unique all-in-phase state. When we say “the all-in-phase state,” we mean the unique one in this invariant hyperplane.

We find it helpful to visualize the oscillators on the unit circle, where oscillator $j$ is positioned at the coordinate $[cos⁡(θj),sin⁡(θj)]$⁠. From this perspective, a network is globally synchronizing if starting from any initial positions (except a set of measure zero), all the oscillators eventually end up at the same point on the circle. Due to periodicity, we assume that the phases, $θj$⁠, take values in the interval $[−π,π)$⁠.

One cannot determine whether a network is globally synchronizing by numerical simulation of (1), as it is impossible to try all initial conditions. Of course, one can try millions of random initial conditions of the oscillators’ phases and then watch the dynamics of (1). However, even if all observed initial states eventually fall into the all-in-phase state, one cannot conclude that the system is globally synchronizing because other stable equilibria could still exist; their basins might be minuscule but could nevertheless have a positive measure. These small basins of attraction are very difficult to detect with random initial conditions.

With that caveat in mind, we note that such numerical experiments have been conducted, and they tentatively suggest that $p≫log⁡(n)/n$ is sufficient for global synchronization.¹⁸ In this paper, we investigate global synchrony via a theoretical study. We show that $p≫log2⁡(n)/n$ is good enough to ensure global synchronization with high probability, improving on $p≫log⁡(n)/n1/3$ proved by Ling et al.,¹⁸ and bringing us closer to the connectivity threshold of Erdős–Rényi graphs.

Although we focus on Erdős–Rényi graphs, many of the inequalities we derive hold for any random or deterministic network. To highlight this, we state many of our findings for a general graph $G$ and a general parameter $p∈R$⁠. In the end, we restrict ourselves back to Erdős–Rényi graphs and take $p$ to be the probability of a connection between any two oscillators.

Our results depend on both the adjacency matrix $A$ and the graph Laplacian matrix $L=D−A$⁠, where $D$ is a diagonal matrix and $Dii$ is the degree of vertex $i$ (counting self-loops). For any $p∈R$⁠, denote the shifted adjacency and the graph Laplacian matrix by

where $Jn$ is the $n×n$ matrix of all ones and $In$ is the $n×n$ identity matrix. It is worth noting that for Erdős–Rényi graphs, the shifts are precisely the expectation of the matrices as $E[A]=pJn$ and $E[L]=−pJn+npIn$⁠. Remarkably, we show that the global synchrony of a network can be guaranteed by ensuring that the spectral norms $‖ΔA‖$ and $‖ΔL‖$ satisfy particular inequalities. The spectral norm of a symmetric matrix is the maximum eigenvalue in absolute magnitude, and $‖ΔA‖$ and $‖ΔL‖$ are extensively studied in the random matrix literature.^25,26 We also find it appealing that the spectral norm of the graph Laplacian matrix appears naturally in our analysis, as it has been used previously to study the dynamics of networks of oscillators²⁷ as well as diffusion on graphs. From a high-level viewpoint, we prove that global synchronization of a graph is a spectral property, like many other graph properties.

While we focus in this paper on global synchronization of identical oscillators, Medvedev and his collaborators have considered other aspects of synchronization for non-identical oscillators on Erdős–Rényi graphs, as well as on other graphs, such as Cayley and Watts–Strogatz graphs, often in the continuum limit.^28–31 Their findings have a similar overarching message that synchronization occurs spontaneously above a critical threshold.

To prove that a random graph is globally synchronizing with high probability, we bridge the gap between spectral graph theory and coupled oscillator theory. The literature contains many good probabilistic estimates for the spectral norm of a random graph’s adjacency and graph Laplacian matrices, which we use to control the long-time dynamics of (1). The key to our proof is to establish conditions on these two spectral norms that force any stable equilibrium to have phases that lie within a half-circle. Confining the phases in this way then guarantees global synchronization because it is known that the only stable equilibrium of (1) with phases confined to a half-circle is the all-in-phase synchronized state.^8,23

Our first attempt at controlling the phases is the inequality we state below as (9), which can be used to guarantee that very dense Erdős–Rényi graphs are globally synchronizing with high probability. As an aside, when $p=1$ inequality (9) together with (4) provides a new proof that a complete graph of identical Kuramoto oscillators is globally synchronizing.⁵ A similar inequality to (9) was derived by Ling et al.¹⁸ to show that $p≫log⁡(n)/n1/3$ is sufficient for global synchrony with high probability.

To improve on (9), our argument becomes more intricate. We carefully examine the possible distribution of edges between oscillators whose phases lie on different arcs of the circle and show that any equilibrium is destabilized if there are too many edges between oscillators that have disparate phases (see Lemma 8 and Fig. 2).

An important quantity in the study of Kuramoto oscillators is the so-called complex order parameter, $ρ1$⁠. The magnitude of $ρ1$ is between $0$ and $1$ and measures the synchrony of the oscillators in the network. We find it useful to also look at second-order moments of the oscillator distribution for analyzing the synchrony of random networks. Higher-order moments are also called Daido order parameters and can be used to analyze oscillators with all-to-all coupling, corresponding to a complete graph.^9,32–34

For an equilibrium $θ=(θ1,…,θn)$⁠, we define the first- and second-order moments as

(For convenience, we use the notation $∑j$ to mean $∑j=1n$⁠.) Without loss of generality, we may assume that the complex order parameter $ρ1$ is real-valued and nonnegative. To see this, write $ρ1=|ρ1|eiψ$ for some $ψ$⁠. Then, $θ^=(θ1−ψ,…,θn−ψ)$ is also an equilibrium of (1) with the same stability properties as $θ$ since (1) is invariant under a global shift of all phases by $ψ$⁠. Therefore, for the rest of this paper, we assume that $ψ=0$ for any equilibrium of interest with $0≤ρ1≤1$⁠. This allows us to select a representative equilibrium point from its associated subspace of equilibria. When $ρ1=1$⁠, the oscillators are in pure synchrony, and when $ρ1≈0$⁠, the phases are scattered around the unit circle without a dominant phase.

To avoid working with complex numbers, it is convenient to consider the quantity $|ρ2|2$⁠. For $m=1,2$⁠, we have

By analyzing $ρ12$ and $|ρ2|2$⁠, one hopes to witness the rough statistics of an equilibrium to understand its potential for synchrony without concern for its precise pattern of phases.

Let $qθ=(eiθ1,…,eiθn)⊤$ and note that

where $A$ is the adjacency matrix, $qθ¯$ is the complex conjugate of $qθ$⁠, and $⊤$ denotes the transpose of a vector. Since $cos2⁡(θk−θj)=12(cos⁡(2(θk−θj))+1)$⁠, we have

where $1$ is the vector of all ones and $q2θ=(ei2θ1,…,ei2θn)⊤$⁠.

We would like to know when a stable equilibrium is close to the all-in-phase state, and we know this when $ρ1$ is close to $1$⁠. Therefore, we begin by deriving a lower bound on $ρ1$ for any stable equilibrium of (1). Similar inequalities involving $ρ12$ and $|ρ2|2$ have been used to demonstrate that sufficiently dense Kuramoto networks are globally synchronizing.^18,21

To maximize the lower bound on $ρ12$ from Lemma 1, one can optimize over $p$⁠. For a random graph where each edge has a fixed probability of being present, independently of the other edges, we usually just select $p$ to be that probability. Regardless, to make the lower bound on $ρ12$ in Lemma 1 useful, we need to find a nontrivial lower bound on $|ρ2|2$ since this quantity appears in the right hand side of (4). To obtain such a bound, we use the following technical lemma.

To see how Lemma 2 can be used to derive a lower bound on $|ρ2|2$ for any stable equilibrium state, we start by dropping the second sum in Lemma 2. Using the upper bound $‖ΔAqθ‖2≤n‖ΔA‖2$⁠, we find that $∑jsin2⁡(θj)≤‖ΔA‖2/(np2ρ12)$⁠. Since $nRe(ρ2)=∑jcos⁡(2θj)=∑j(1−2sin2⁡(θj))$⁠, we have the following lower bound on $|ρ2|$⁠:

From (4) and (7), we can deduce that when $‖ΔA‖≪np$⁠, then $ρ1$ and $|ρ2|$ must both be close to $1$⁠. Intuitively, this should mean that the corresponding stable equilibrium, $θ$⁠, must be close to the all-in-phase state with the possible exception of a small number of stray oscillators. However, our goal is to prove global synchrony, which is a more stringent condition, and we must completely rule out the existence of any stray oscillators.

To precisely control the number of stray oscillators, we define a set of indices for oscillators whose phases lie outside of a sector of half-angle $ϕ$ (centered about the all-in-phase state),

for any angle $0≤ϕ≤π$ (see Fig. 1). If we can prove that $Cπ/2$ is empty, then we know that all the phases in the equilibrium state lie strictly inside a half-circle. That would give us what we want because by a basic theorem for the homogeneous Kuramoto model, the only stable equilibrium $θ$ with all its phases confined to a half-circle is the all-in-phase state.^8,23 In terms of bounds, since $|Cπ/2|$ is integer-valued, if we can show that $|Cπ/2|<1$⁠, then we know that $Cπ/2$ is the empty set.

From Lemma 2 and $‖ΔAqθ‖2≤n‖ΔA‖2$⁠, we find that

where $δ|θj|>π/2$ is $1$ if $|θj|>π/2$ and $0$ otherwise. Therefore, by plugging in $ϕ=π/2$⁠, we see that $|Cπ/2|≤‖ΔA‖2/(np2ρ12)$⁠. Thus, if $‖ΔA‖2<np2ρ12$⁠, then $Cπ/2$ must be the empty set and the network is globally synchronizing.

Unfortunately, this kind of reasoning is not sufficient to prove global synchrony for graphs of interest to us here. For example, for an Erdős–Rényi random graph, we know that $‖ΔA‖2≈4p(1−p)n$ with high probability for large $n$ (see Sec. VI). Therefore, the upper bound on $|Cπ/2|$ is approximately $4(1−p)/(pρ12)$⁠, which for $p<1/2$ is certainly not good enough to conclude that $Cπ/2$ is empty. Instead, we must further improve our bounds on $|Cπ/2|$ by using a recursive refinement strategy that we refer to as an “amplification” argument (see Sec. V).

The precise amplification argument that we use requires bounds on the number of edges and sizes of vertex sets of a graph expressed in terms of the spectral norms of $ΔA$ and $ΔL$⁠. It is worth noting that these bounds hold for both deterministic and random graphs. For a vertex set $C$ of a graph $G$⁠, we denote the characteristic vector by $vC$⁠; i.e., $(vC)j=1$ if $j∈C$ and $(vC)j=0$ if $j∉C$⁠. We use $vC⊤$ to denote the vector transpose of $vC$⁠. We write $|C|$ to be the number of vertices in $C$ and denote the number of (directed) edges in $G$ starting at vertex set $C$ and ending at $C′$ as $EC,C′$⁠. Therefore, $EC,C$ is twice the number of (undirected) edges between vertices in $C$⁠. For Erdős–Rényi graphs, one expects to have $EC,C′≈p|C||C′|$⁠. However, for our argument to work, expectations are not adequate; instead, we need to have bounds on the difference between $EC,C′$ and $p|C||C′|$⁠. The results in this section are proved by classical techniques, and closely related bounds are well known. We give the proofs of the precise inequalities that we need so that the paper is self-contained.

We first show that for any vertex set $C$⁠, the number of edges connecting a vertex in $C$ to another vertex in $C$ deviates from $p|C|2$ by at most $‖ΔA‖|C|$ for any $0≤p≤1$⁠.

Lemma 3 controls the number of edges connecting a set of vertices. The following result controls the number of edges leaving a vertex set. We denote the vertices of $G$ that are not in $C1$ by $VG∖C1$⁠.

When a bound on $‖ΔL‖$ is available, Lemma 4 can be used to ensure that $G$ is connected. In particular, Lemma 4 tells us that a graph of size $n$ is connected if $‖ΔL‖<np$ for some $0≤p≤1$⁠.

Since $0≤EC1,C1¯≤|C1||C1¯|$⁠, we know that $|EC1,C1¯−p|C1||C1¯||≤max{p,1−p}|C1||C1¯|$⁠. Therefore, Lemma 4 is only a useful bound when $‖ΔL‖≪np$⁠. We can take Lemma 4 a step further and bound the number of edges between any two disjoint sets of vertices of $G$ using $‖ΔL‖$⁠. We denote the union of the vertex sets $C1$ and $C2$ by $C1⊔C2$⁠.

Lemmas 3 and 5 can be combined to obtain the following result. We do not regard Theorem 6 as having significant independent interest, but later, it is an important statement in the proof of Corollary 9. It only deserves the status of a theorem in the sense that it is a key technical step in our overall argument.

We are finally ready for our amplification argument, which is a way to improve the bounds on $|Cπ/2|$ from (9). We first write down a new inequality that holds for any stable equilibrium. We write it down using a kernel function $K$ that later allows us to improve our argument with the aid of a computer (see Sec. IV D). For any stable equilibrium $θ=(θ1,…,θn)$ of (1), we know that it satisfies equilibrium and stability conditions. The equilibrium conditions are given by $∑jAjksin⁡(θk−θj)=0$ for all $1≤k≤n$⁠. In addition, stability requires that the Hessian matrix [see (2.3) of Ling et al.¹⁸] associated with $θ$ is a nonnegative definite matrix; i.e., it has nonnegative eigenvalues.

The fact that the Hessian matrix has nonnegative eigenvalues can be used to form many useful inequalities. For example, if $HA(θ)$ is the Hessian matrix, then a stable equilibrium must satisfy the Rayleigh quotient $x⊤HA(θ)x≥0$ for any $n×1$ vector $x$⁠. In particular, selecting $x=ek$ (the $k$th unit vector) implies that the $k$th diagonal entry of $HA(θ)$ is nonnegative. This yields the inequalities $ek⊤HA(θ)ek=∑jAjkcos⁡(θk−θj)≥0$ for $1≤k≤n$⁠. Next, one can combine $∑jAjksin⁡(θk−θj)=0$ and $∑jAjkcos⁡(θk−θj)≥0$ in tricky ways and then use trigonometric identities to simplify the inequalities. After manipulating these inequalities for a while, we find that it is beneficial to massage them differently depending on the values of $θk$ and $θj$ that appear in them. However, then, having combined the inequalities in diverse ways, we notice that they all share a common structure. To highlight this commonality, we express the various inequalities as a unified family of the form $∑jAjkK(θj,θk)≥0$⁠, where $K$ is defined piecewise to reflect its dependence on its arguments.

In preliminary work, we have found indications that Lemma 7 can be used to prove stronger results than those we report below; see Sec. IV D for further discussion. However, we are not sure yet how to write down an argument that uses the full strength of Lemma 7 in a readily digestible fashion. So for now, we use the following simplified lemma instead. It is a key step in our amplification argument.

Lemma 7 shows us that if $θ$ is a stable equilibrium, then any oscillator with phase $>π/2$ in absolute value cannot be coupled to too many oscillators with phases $<π/2$⁠. This is because $K(α,β)=−cos⁡(α)<0$ when $|α|<π/2$ and $|β|>π/2$⁠. If there are too many negative contributions in the sum $∑jAjkK(θj,θk)$⁠, then it will end up negative, which is not allowed for a stable equilibrium. To make this precise, we can use our sets $Cβ$ in (8).