On the Critical Coupling of the Finite Kuramoto Model on Dense Networks
OOn the Critical Coupling of the FiniteKuramoto Model on Dense Networks
Shuyang Ling ∗ April 8, 2020
Abstract
Kuramoto model is one of the most prominent models for the synchronization ofcoupled oscillators. It has long been a research hotspot to understand how naturalfrequencies, the interaction between oscillators, and network topology determinethe onset of synchronization. In this paper, we investigate the critical coupling ofKuramoto oscillators on deterministic dense networks, viewed as a natural gener-alization from all-to-all networks of identical oscillators. We provide a sufficientcondition under which the Kuramoto model with non-identical oscillators has oneunique and stable equilibrium. Moreover, this equilibrium is phase cohesive andenjoys local exponential synchronization. We perform numerical simulations of theKuramoto model on random networks and circulant networks to complement ourtheoretical analysis and provide insights for future research.
The synchronization problem of coupled oscillators has a long history, dating back toChristiaan Huygens in 1665 who investigated the behavior of two pendulum clocks mountedside by side on the same support [2, 24]. Since then, the study of the synchronizationphenomenon has attracted a large amount of attention across various scientific areas in-cluding mathematics, physics, neuroscience, and engineering. Kuramoto model is oneclassical model of the synchronization of coupled oscillators [16, 17]. Kuramoto con-sidered n fully connected oscillators on a torus whose dynamics is characterized by anordinary differential equation:˙ θ i ( t ) = ω i − Kn n (cid:88) j =1 sin( θ i − θ j ) , ≤ i ≤ n (1.1)where K is called the coupling strength between these oscillators and ω i is the naturalfrequency of the i th oscillator.One central question about the Kuramoto model is to understand when the oscilla-tors { θ i ( t ) } ni =1 will synchronize, the answer of which depends on K and the strength ofnatural frequencies { ω i } ni =1 . This question has sparked extensive research. Moreover, theKuramoto model has already found numerous applications such as electric power net-works, neuroscience, chemical oscillations, spin glasses, see [1, 2, 4, 11, 12, 13, 25] andthe references therein for more details. ∗ Division of Data Science, New York University Shanghai, Shanghai, China, (Email: [email protected]) a r X i v : . [ n li n . AO ] A p r n this paper, we consider the Kuramoto model on more general networks,˙ θ i ( t ) = ω i − Kn n (cid:88) j =1 a ij sin( θ i − θ j ) , ≤ i ≤ n (1.2)where A = ( a ij ) ≤ i,j ≤ n is the adjacency matrix of the underlying network connecting theseoscillators, i.e., two oscillators i and j are linked if and only if a ij = 1. We are interestedFigure 1: The dynamics of Kuramoto model with non-identical oscillatorsin the synchronization phenomenon and the long-term behavior of this dynamical system.In particular, we will focus on answering this question: Does there exist a frequency synchronized solution? If so, is it unique? (1.3)Unlike (1.1) where the network is complete, the answer to (1.3) will depend on networktopology A besides the coupling strength K and the dissimilarities of natural frequencies. As discussed briefly, our work focuses on answering the question (1.3): the existence anduniqueness of frequency synchronized solution of finite
Kuramoto model on deterministicdense networks where each vertex has sufficiently many neighbors . As the synchronizationof coupled oscillators is such an important research topic, there inevitably exists anextensive body of literature on it. As a result, we are unable to give an exhaustiveliterature review here. Instead, we will review the related literature and point out howthese previous works inspire ours.The study of collective synchronization of coupled oscillators has found quite manyapplications in biology, physics, and engineering [4, 23, 26, 27]. Starting from the workby Winfree [35], researchers are interested in developing mathematical models and ap-proaches to analyze and explain the synchronization phenomena. One of the most suc-cessful attempts was made by Kuramoto in [16, 17], who proposed a simplified yet stillhighly non-trivial model of coupled oscillators (1.1). Kuramoto showed that if the cou-pling K is smaller than a certain threshold K c , then the oscillators are incoherent andtransiting to synchrony as K exceeds the threshold K c . This critical coupling thresholdis determined by the strength of the natural frequency ω i of each oscillator. Since then,a lot of progress has been made to advance our understanding of the Kuramoto model,see [1, 2, 27] for excellent reviews on this topic. Most of the analyses rely on tools fromstatistical physics such as calculating the thermodynamical limit of the Kuramoto model.As pointed out in [27], the Kuramoto model with finite oscillators exhibits very dif-ferent behavior from the scenario where the number of oscillators goes to infinity. It2s one major research problem to understand what the critical coupling is for the finiteKuramoto model. In [9], D¨orfler and Bullo gave a comprehensive and excellent studyon the critical coupling of the finite Kuramoto model on a complete graph. In partic-ular, [9] provided a sufficient and necessary condition on the coupling strength for theoscillators to achieve synchronization within the cohesive region where the cohesivenessis later formally defined in (2.1). Finite Kuramoto model on general complex networksare discussed in [1, 7, 10, 12, 13, 15, 25, 33]. A review of recent progresses can be foundin [1, 10, 12, 25, 28]. The work [15] by Jadbabaie, Motee, and Barahona first estimated thecritical coupling of the Kuramoto model on arbitrary complex networks for local conver-gence and showed the existence and uniqueness of stable fixed point within the cohesiveregion. Later on, various necessary and sufficient conditions are obtained [7, 10, 12, 32] forthe estimation of the critical coupling. In particular, [10, 12] established sufficient condi-tions which guarantee frequency synchronization and the existence of locally exponentialstable fixed point within the phase-cohesive region. The work [13] proposed a concise andclosed-form condition for synchronization for a large family of networks. These sufficientconditions are formulated by using the second smallest eigenvalue of graph Laplacianassociated to the underlying network and the natural frequency { ω i } ni =1 . The criticalcoupling of the Kuramoto model is also studied on complete bipartite graphs [33] andrandom graphs [6, 14, 19].It is important to note that Kuramoto model is closely related to the gradient flow of E ( θ ) := K n (cid:88) ≤ i,j ≤ n a ij (1 − cos( θ i − θ j )) − n (cid:88) i =1 ω i θ i (1.4)and thus (1.2) has the following equivalent form:d θ d t = −∇ θ E ( θ ) , θ (0) = θ where θ is the initial state. In other words, the Kuramoto oscillators move in the direc-tion which makes E ( θ ) decrease fastest. Recently, [18, 20, 29, 30] have studied the globalphase synchronization for the homogeneous case ω i = 0, which is actually a special caseof general Kuramoto model. Under ω i = 0, E ( θ ) in (1.4) is a nonnegative energy functionand also has been found related to nonconvex optimization approach applied to groupsynchronization problem on complex networks and community detection [3, 18]. Herethe main focus is to understand how the energy landscape of E ( θ ), i.e., the number ofstable fixed points of homogeneous Kuramoto dynamics, depends on the network topol-ogy. Two main families of networks are investigated: (i): deterministic dense networkand circulant networks; (ii): Erd˝os-R´enyi random graphs [31]. The work [29] by Taylorfirst proved that homogenous Kuramoto enjoys global synchronization from almost anyinitialization if the degree of each node is at least 0 . n . Later, [18] has proven that0.7929 n suffices to guarantee the uniqueness of a globally stable synchronized state, whichwas recently improved further to 0.7889 n in [20] with a more refined analysis. On theother hand, [5, 34] have constructed very interesting examples that the uniform twistedstates are possibly stable fixed points for homogeneous Kuramoto model on circulantnetworks. More precisely, there always exists a circulant network whose degree is smallerthan 15 n/ ≈ . n such that the corresponding homogeneous Kuramoto oscillatorshave a stable equilibrium which is not phase-synchronized. A recent conjecture is pro-posed in [30], stating that the critical threshold could be 0 . n . The discussion on theglobal synchronization of homogeneous Kuramoto model has been generalized to othermanifolds such as n -sphere and Stiefel manifold [21, 22].3ur work is motivated by a series of works on the critical coupling of the finite Ku-ramoto model [10, 12, 15], and the recent progress on the uniqueness of stable phase syn-chronized solution of the homogeneous Kuramoto model on dense networks [18, 20, 29, 30].Our main contribution is on establishing a sufficient condition under which the inhomo-geneous Kuramoto oscillators on dense networks have a unique and stable equilibriumwhich also is phase-cohesive and satisfies local exponential synchronization. This workis a generalization of the work [9] by D¨orfler and Bullo from all-to-all networks to densenetworks. Our work also extends the previous analysis of the energy landscape of theKuramoto model with identical oscillators [18, 20, 29] to the case with heterogeneousoscillators. Moreover, our numerical experiments provide insights for several future di-rections which are well worth exploring. We organize this paper as follows: we will discuss the basics of synchronization andpresent our main theorem in Section 2. Section 3 focuses on numerical simulations andthe proof will be given in Section 4.
We introduce notations which will be used throughout the paper. Vectors are denotedby boldface lower case letters, e.g., z . For any vector z , (cid:107) z (cid:107) := (cid:112)(cid:80) ni =1 z i denotes its (cid:96) -norm; (cid:107) z (cid:107) ∞ := max ≤ i ≤ n | z i | stands for the (cid:96) ∞ -norm. For a given matrix Z , ddiag( Z )denotes a diagonal matrix whose diagonal entries are the same as those of Z ; Z (cid:62) is thetranspose of Z . Let I n and n be the n × n identity matrix and a column vector of“1” in R n respectively. For any symmetric matrix Z ∈ R n × n , we denote Z (cid:23) Z ispositive semidefinite. For any two matrices X and Z of the same size, their inner productis denoted by (cid:104) X , Z (cid:105) := (cid:80) i,j X ij Z ij and X ◦ Z denotes their Hadamard product, i.e.,( X ◦ Z ) ij = X ij Z ij . We first review the basics of synchronization before moving to our main results. One canfind more discussion on these concepts from many sources such as [4, Chapter 16].
Definition 2.1 (Phase synchronization) . A solution θ ( t ) ∈ R n achieves phase synchro-nization if θ i ( t ) = θ j ( t ) for all t ≥ . It is well-known that phase synchronization is possible only if ω i = ω j , see [4]. Oth-erwise, { θ : θ i = θ j , i (cid:54) = j } is not even a fixed point of (1.1). For the Kuramoto modelwith non-identical oscillators, we focus on frequency synchronization. Definition 2.2 (Frequency synchronization) . A solution θ ( t ) ∈ R n achieves frequencysynchronization if ˙ θ i ( t ) = ˙ θ j ( t ) for all t ≥ . Definition 2.3 (Phase cohesiveness) . A solution θ ( t ) ∈ R n is γ -phase cohesive if thereexists a number γ ∈ [0 , π ) such that θ ( t ) ∈ ∆( γ ) for all t ≥ where ∆( γ ) := { θ : | θ i − θ j | < γ, ∀ i (cid:54) = j } . (2.1) In other words, an arc of length γ contains all phases θ i ( t ) on the unit circle for any t ≥ . ω i by ω i − n (cid:80) nj =1 ω j . In this way, we have n (cid:88) i =1 ˙ θ i ( t ) = n (cid:88) i =1 ω i = 0 . In other words, a frequency-synchronized solution θ ( t ) is a fixed point of (1.2),˙ θ i ( t ) = 0 ⇐⇒ Kn n (cid:88) j =1 a ij sin( θ i − θ j ) = ω i , ≤ i ≤ n. (2.2)The stability is determined by the spectra of Jacobian matrix. The Jacobian matrix isequal to − n − K · J ( θ ) where J ( θ ) := ddiag( AQQ (cid:62) ) − A ◦ QQ (cid:62) . (2.3)Here A = ( a ij ) ≤ i,j ≤ n is the adjacency matrix, Q ∈ R n × with its i th row q i ∈ [cos( θ i ) , sin( θ i )],and ( A ◦ QQ (cid:62) ) ij := a ij cos( θ i − θ j ) is the Hadamard product of A and QQ (cid:62) . The ( i, j )-entry of J ( θ ) is ( J ( θ )) ij = (cid:40)(cid:80) nj =1 a ij cos( θ i − θ j ) , i = j, − a ij cos( θ i − θ j ) , i (cid:54) = j, which is actually a Laplacian matrix associated with the signed weights { a ij cos( θ i − θ j ) } i,j .From the construction, we know that θ ( t ) is a stable frequency synchronized solution if (cid:107) ˙ θ (cid:107) = 0 and the second smallest eigenvalue λ ( J ( θ )) > J ( θ ) . As discussed before, a complete understanding of synchronization of the Kuramotomodel on general complex networks is still one major open problem in the field. Ourdiscussion will mainly focus on deterministic dense networks, which are defined formallyas follows.
Definition 2.4 (Deterministic dense networks) . A network with symmetric adjacencymatrix A = ( a ij ) ≤ i,j ≤ n ∈ { , } n × n is µ -dense if the degree of each node is at least µ ( n − , ≤ µ ≤ . In particular, if µ = 1 , the corresponding network is complete. The Kuramoto model on dense networks is a natural generalization of the classicalKuramoto model with all-to-all coupling. Now we return to the main questions (1.3):(a) When do the coupled oscillators on µ -dense network have a stable frequency synchro-nized solution?(b) If such a solution exists, is it unique?The answer depends on the parameter µ , critical coupling parameter K , as well as thestrength of natural frequency ω i . This will be the main focus of our paper.
Theorem 2.1.
Consider the Kuramoto model (1.2) with natural frequency { ω i } ni =1 and (cid:80) ni =1 ω i = 0 . Suppose (cid:107) ω (cid:107) ∞ := max ≤ i ≤ n | ω i | < K (cid:32)(cid:114) µ −
34 + µ − (cid:33) , (2.4) then there exists a unique frequency-synchronized solution. Moreover, this solution islocated within the phase cohesive region ∆( π/ and the coupled oscillators achieve localexponential frequency synchronization. emark 2.2. From (2.4) , we see that µ ≥ −√ is implicitly assumed since the righthand side of (2.4) is nonnegative. Therefore, this network is always connected. We briefly discuss what this theorem implies by looking at two interesting specialcases: (i) µ = 1; (ii) (cid:107) ω (cid:107) ∞ = 0, and make comparisons with the state-of-the-art results.(i) Heterogeneous oscillators over a complete graph µ = 1:The model (1.2) immediately reduces to the classical Kuramoto model (1.1). Thenthe condition (2.4) reads K > (cid:107) ω (cid:107) ∞ . Note that [9, Theorem 4.1] shows that K ≥ ω max − ω min suffices to guarantee a local exponentially stable frequency synchro-nized solution within the phase cohesive region ∆( γ ) for some γ < π/
2. Comparedwith [9], our result is slightly looser and becomes tight if ω max + ω min = 0 . (ii) Homogeneous Kuramoto model on arbitrary unweighted networks:If ω i = 0 holds for all 1 ≤ i ≤ n , then (2.4) turns into (cid:114) µ −
34 + µ − > ⇐⇒ µ > − √ ≈ . . In other words, the stable equilibrium is unique which is exactly the phase syn-chronized solution, i.e., θ i = θ j , if µ > . µ > . µ < / ≈ . µ -dense circulant networksuch that the corresponding Kuramoto model has multiple stable equilibria, i.e.,phase/frequency synchronized solutions. In other words, for µ < /
22, there al-ways exist multiple stable frequency synchronized solutions if ω i is sufficiently small.Adding tiny natural frequencies to the homogeneous Kuramoto model, regarded asa small perturbation, will not change the number of stable solutions. We provide anumerical illustration in Figure 4. One natural question is the tightness of the bound (2.4) in Theorem 2.1 as well as thepossible extension of this result to other families of networks. We address this issueby providing numerical examples. In particular, we will focus on two types of graphs:Erd˝os-R´enyi random graphs and circulant networks.
Note that the condition in (2.4) only holds for very dense networks, viewed as the worst-case scenario. It is natural to ask what the critical coupling is for Erd˝os-R´enyi (ER)random graphs. We denote ER graph as G ( n, p ) if the network is of size n and its( i, j )-entry a ij of the symmetric adjacency matrix is a Bernoulli random variable withparameter p. We want to study how the frequency synchronization depends on p and the strengthof natural frequency. We let network contain n = 100 vertices and K = 1, and sample A from G ( n, p ) with p varying from 0.05 to 1. For each random instance, we generateuniformly distributed ω i over [ − κp, κp ] for κ = 0 , . , · · · , . ,
1. In other words, therange of ω i is approximately proportional to the expected degree np of each node. For eachnetwork, we use Euler scheme to simulate the trajectory with step size τ = 2 · − . The6imulation stops if the iterate converges to an approximate stable frequency-synchronizedsolution, i.e., (cid:107) ˙ θ (cid:107) < − , λ ( J ( θ )) ≥ − , or the iteration number reaches 10 . We run 20 experiments for each pair of ( p, κ ), andcalculate the proportion of the iterate converging to a stable synchronized solution.The phase transition plot is in Figure 2: the white region stands for success and blackmeans failure. We can see if κ = (cid:107) ω (cid:107) ∞ /p < .
75 and p > .
2, the dynamics will finallysynchronize with high probability. Empirically, it seems κ does not heavily depend on p . We also compute the length of the shortest arc covering all θ i . Figure 2 indicatesthat if κ < .
7, i.e., (cid:107) ω (cid:107) ∞ ≤ . p , the oscillators will converge to frequency-synchronizedsolution in the cohesive region ∆( π/ . In fact, Figure 2 matches our previous numericalstudy in [18]. For the homogeneous Kuramoto model κ = 0 ( ω i = 0), the oscillatorsalways achieve global synchronization for p > n − log n , i.e., as long as the ER graph isconnected with high probability, there exists a synchronized solution.Compared with the bound in (2.4), it shows that our existing bound is relatively con-servative since it does not assume an additional probabilistic structure on the networks.We believe one can obtain a refined bound on the critical coupling threshold by takingthis extra prior information into account. Figure 2: Phase transition plot for frequency synchronization solutions: ER graphs andrandom initialization.
The homogeneous Kuramoto model on circulant networks exhibits very interesting be-haviors. It has been shown that the uniform twisted state θ twist := 2 πn − [1 , · · · , n ] is astable equilibrium for a large family of circulant networks such as Wiley-Strogatz-Girvan(WSG) networks [34] and Harary graphs [5]. Now we will focus on the critical couplingthresholds on WSG networks.The WSG( k ) networks are constructed by starting with an n -cycle graph and thenconnecting each node with its closest k neighbors. As a result, the degree of each nodeis 2 k. We set up numerical experiments as follows: let k = (cid:98) np/ (cid:99) , i.e., the largestinteger smaller than or equal to np/
2, for p = 0 . , . , · · · ,
1. Note that WSG( k ) isalso a µ -dense network with µ = p. We generate the natural frequencies ω i according touniform distribution [ − κp, κp ] for κ ranging from 0 , . , · · · , . For each pair of ( κ, p ),we run 20 experiments and count the number of cases in which the iterates converge to7 stable synchronized solution. In particular, to address the dependence of the dynamicson initialization, two types of initial states are chosen: (i) θ (0) is uniformly distributedover [0 , π ]; (ii) θ (0) = πn [1 , · · · , n ] is a uniform twisted state.From Figure 3, we can see that the phase transition plot looks similar to Figure 2for random initialization but the white region shrinks: on circulant networks, randominitialization will lead to a frequency-synchronized solution if p ≥ . (cid:107) ω (cid:107) ∞ ≤ . p while p ≥ . (cid:107) ω (cid:107) ∞ ≤ . p suffice for ER graphs.However, if the initialization is a uniform-twisted state , Figure 5 looks very differentfrom Figure 3. There exists a large set of ( κ, p ) in Figure 5 such that starting from theuniform twisted state, the trajectory does not converge to a cohesive solution. We presenta more concrete example in Figure 4: Figure 4 demonstrates that the oscillators reacha stable non-cohesive frequency-synchronized solution outside ∆( π/
2) for small ω i , withinitialization θ twist . However, this stable non-cohesive solution disappears if the intrinsicfrequency gets stronger. Instead, they converge to a stable solution within ∆( π/ . Figure 3: Phase transition plot for frequency synchronization solutions: WSG circulantnetworks with degree np and uniformly random initialization.Here is one explanation for Figure 5: if ω i = 0, the uniform-twisted state is alwaysa stable equilibrium (with strictly positive second smallest eigenvalue of J ( θ )) if k issmaller than 0.34 in [34]. If ω i is small compared with the second smallest eigenvalue λ ( J ( θ )), this perturbation on ω i does not change dynamics qualitatively. The uniform-twisted state still lies in the basin of attraction and the oscillators will converge to anon-cohesive frequency synchronized solution if starting from θ twist . On the other hand,if the intrinsic frequency gets larger, the iterate will escape from the basin of attractionaround the uniform twisted solution and reach another synchronization solution if thereexists one. If ω i gets too strong, a stable synchronized state no longer exists.That is why the left plot in Figure 5 shows that for fixed p < .
6, the oscillators syn-chronize to a non-cohesive solution for small ω i ; as ω i gets stronger, it starts to becomeincoherent; if the strength of ω i exceeds certain threshold, the oscillators start synchroniz-ing to a cohesive synchronized solution; but when (cid:107) ω (cid:107) ∞ is larger than 0 . p , the system nolonger synchronizes. The two plots in Figure 5 indicate that there exists a boundary for( κ, p ) which distinguishes a cohesive synchronized state from to a non-cohesive one. Weleave the explanation of this boundary for future research. On the other hand, if p > . θ (0) = θ twist yields a phase-cohesive synchronized solution for (cid:107) ω (cid:107) ∞ < . p. It is because θ twist is not a stable equilibrium for the homogeneous Kuramoto model if8 > .
68; also the interaction between oscillators gets stronger than the effect of naturalfrequency due to the higher network connectivity, which makes the global synchronizationmore likely. -1 -0.5 0 0.5 1-1-0.8-0.6-0.4-0.200.20.40.60.81 -1 -0.5 0 0.5 1-1-0.8-0.6-0.4-0.200.20.40.60.81
Figure 4: Phase cohesiveness. Left: p = 0 . (cid:107) ω (cid:107) ∞ = 0 .
1. In this case, a non-cohesivesynchronized state exists for small ω i ; Right: p = 0 . (cid:107) ω (cid:107) ∞ = 0 .
2; For large ω i ,setting θ (0) = θ twist still leads to a non-cohesive stable equilibrium. Figure 5: Phase transition plot for frequency synchronization solutions: WSG circulantnetworks with degree np and initialization θ (0) = θ twist . Define the unnormalized order parameter r ( θ ) = n (cid:88) j =1 e i θ j whose magnitude equals 1 if all oscillators are fully synchronized.9 emma 4.1. Consider the Kuramoto model (1.2) with natural frequency { ω i } ni =1 and (cid:80) ni =1 ω i = 0 . Suppose (2.4) holds, then all the stable equilibria, if existing, must be in acohesive region ∆( π/ . Proof of Lemma 4.1.
The proof generalizes the idea in [18, Theorem 3.1]: first show thatif an equilibrium is stable, its order parameter is sufficiently large, which is made possi-ble by using the second-order necessary condition; then the first-order critical conditionguarantees that all oscillators are inside the cohesive region ∆( γ ) for some γ < π/ θ is a stable equilibrium, then the Jacobian matrix is negative definite andthus J ( θ ) = diag( AQQ (cid:62) ) − A ◦ QQ (cid:62) (cid:23) J ( θ ) is defined in (2.3). As aresult, (cid:104) J ( θ ) , QQ (cid:62) (cid:105) ≥ ⇐⇒ (cid:104) A , QQ (cid:62) (cid:105) ≥ (cid:104) A , QQ (cid:62) ◦ QQ (cid:62) (cid:105) (4.1)since QQ (cid:62) (cid:23) (cid:107) r ( θ ) (cid:107) = n (cid:88) i =1 n (cid:88) j =1 cos( θ i − θ j ) = (cid:104) n (cid:62) n , QQ (cid:62) (cid:105) (4.2)= (cid:104) A , QQ (cid:62) (cid:105) + (cid:104) n (cid:62) n − A , QQ (cid:62) (cid:105)≥ (cid:104) A , QQ (cid:62) ◦ QQ (cid:62) (cid:105) + (cid:104) n (cid:62) n − A , QQ (cid:62) (cid:105) = (cid:104) n (cid:62) n , QQ (cid:62) ◦ QQ (cid:62) (cid:105) + (cid:104) n (cid:62) n − A , QQ (cid:62) − QQ (cid:62) ◦ QQ (cid:62) (cid:105) (4.3)where the first inequality follows from (4.1) and n is the n × QQ (cid:62) is rank-2, positive semidefinite, and has its trace equal to n .Therefore, (cid:104) n (cid:62) n , QQ (cid:62) ◦ QQ (cid:62) (cid:105) = (cid:107) QQ (cid:62) (cid:107) F ≥ | Tr( QQ (cid:62) ) | = n (cid:107) · (cid:107) F stands for the matrix Frobenius norm. For the second term in (4.3), wesimply use the fact that | ( QQ (cid:62) ) ij | ≤ (cid:104) A , n (cid:62) n (cid:105) is at least µn ( n −
1) + n since A is µ -dense. Then (cid:104) n (cid:62) n − A , QQ (cid:62) − QQ (cid:62) ◦ QQ (cid:62) (cid:105) ≥ − (cid:104) n (cid:62) n − A , n (cid:62) n (cid:105)≥ − n + 2 µn ( n −
1) + 2 n = − − µ ) n + 2(1 − µ ) n. Therefore, we get (cid:107) r ( θ ) (cid:107) ≥ n − − µ ) n + 2(1 − µ ) n = (cid:18) µ − (cid:19) n + 2(1 − µ ) n. (4.4)Next, we bound the phase of the fixed point (if it exists). Let’s compute the phasedifference between r ( θ ) and i th oscillator e i θ i re − i θ i = n (cid:88) j =1 e i( θ j − θ i ) = n (cid:88) j =1 (cos( θ j − θ i ) + i sin( θ j − θ i )) = (cid:107) r (cid:107) · e i( ψ − θ i ) where ψ is the phase of r ( θ ) . By taking the imaginary part, it holds that (cid:107) r (cid:107) · sin( ψ − θ i ) = n (cid:88) j =1 sin( θ j − θ i )= n (cid:88) j =1 a ij sin( θ i − θ j ) + n (cid:88) j =1 (1 − a ij ) sin( θ i − θ j )= nK − ω i + n (cid:88) j =1 (1 − a ij ) sin( θ i − θ j )10here the last inequality uses ˙ θ i = 0 and (cid:80) nj =1 a ij sin( θ i − θ j ) = nK − ω i . Therefore, weobtain (cid:107) r (cid:107) · | sin( ψ − θ i ) | ≤ nK − (cid:107) ω (cid:107) ∞ + (1 − µ )( n −
1) (4.5)where the right hand side is independent of the index i . Combining (4.5) with (4.4) leadsto an upper bound of | sin( ψ − θ i ) | and thus | ψ − θ i | , | sin( θ i − ψ ) | ≤ nK − (cid:107) ω (cid:107) ∞ + (1 − µ )( n − (cid:113)(cid:0) µ − (cid:1) n + 2(1 − µ ) n , ≤ i ≤ n. (4.6)We claim that under (2.4), { θ i } ni =1 are inside two disjoint quadrants. θ i ∈ (cid:110) θ : ψ − π < θ < ψ + π (cid:111) ∪ (cid:26) θ : ψ + 3 π < θ < ψ + 5 π (cid:27) . (4.7)To ensure this, it suffices to bound (4.6) by 1 / √ (cid:107) ω (cid:107) ∞ K ≤ (cid:114) µ −
34 + 1 − µn + µ − . In fact, as will be shown below, all { θ i } ni =1 are in the same quadrant if θ is a stableequilibrium. Otherwise, we can construct a test vector such that the quadratic formw.r.t. J ( θ ) is smaller than 0, as discussed in [29, 18]. We provide the proof here to makethe presentation self-contained. Remember the stability is directly to the second smallesteigenvalue of J ( θ ) since the eigenvalues of Jacobian matrix and J ( θ ) have completelyopposite signs. Suppose not every θ i is located in the same quadrant, then we define z = Γ − Γ c ∈ R n , z i = (cid:40) , i ∈ Γ , − , i ∈ Γ c , where Γ = { i : − π/ < θ i − ψ < π/ } and Γ c = { i : 3 π/ < θ i − ψ < π/ } . If so, z (cid:62) J z = (cid:88) ≤ i,j ≤ n a ij cos( θ i − θ j )( z i − z j ) = 8 (cid:88) i ∈ Γ ,j ∈ Γ c a ij cos( θ i − θ j ) < a ij are nonnegative with at least one a ij = 1 for some i ∈ Γ and j ∈ Γ c , andcos( θ i − θ j ) < . This implies that the smallest eigenvalue of J ( θ ) is strictly less than 0while J ( θ ) is positive semidefinite if θ is a stable equilibrium, which is a contradiction.Now, we wrap up our discussion here: we have shown under (2.4), if there exists astable equilibrium, all { θ i } ni =1 are located in the same quadrant, which means | θ i − θ j | < π , ∀ i (cid:54) = j, following from triangle inequality. Lemma 4.2. If { θ i } ni =1 are initialized in ∆( γ ) with γ < π/ , then they will stay inside ∆( γ ) if K sin( γ ) > ω max − ω min µ − , µ > . (4.8) In other words, ∆( γ ) is a basin of attraction. roof: The main idea of proof is adapted from that in [9]. We aim to show thatunder (4.8) stated in the lemma, the range of the phase oscillators will not increase ifthey are initialized in the cohesive region ∆( γ ) . First of all, we assume there exists an arc of length γ covering all { θ i } ni =1 . Let i =argmax ≤ k ≤ n θ k and j = argmin ≤ k ≤ n θ k . Note that i and j are unnecessarily unique,i.e., there could be multiple indices such that the maximum and minimum are attainedrespectively. For any ( i, j ), we have γ = θ i − θ j > γ is the range of all oscillators. We will investigate how θ i − θ j changeswith respect to the time t . From (1.2), we have˙ θ i ( t ) − ˙ θ j ( t ) = ω i − ω j − Kn n (cid:88) k =1 ( a ik sin( θ i − θ k ) + a jk sin( θ k − θ j )) . (4.9)The goal is to find an upper bound for (4.9). We claim that˙ θ i ( t ) − ˙ θ j ( t ) ≤ ω max − ω min − (2 µ − K sin( γ ) . (4.10)Note that ω i − ω j ≤ ω max − ω min holds, and thus the key is to bound the third termin (4.9). Apparently, since θ i − θ k and θ k − θ j are in [0 , γ ] ⊆ [0 , π/ µ ( n −
1) neighbors.
Case 1: If a ij = 1, then n (cid:88) k =1 ( a ik sin( θ i − θ k ) + a jk sin( θ k − θ j ))= (cid:88) k (cid:54) = i,j ( a ik sin( θ i − θ k ) + a jk sin( θ k − θ j )) + 2 sin( γ ) ≥ (cid:88) { k : a ik = a jk =1 }\{ i,j } (sin( θ i − θ k ) + sin( θ k − θ j )) + 2 sin( γ )where sin( θ i − θ j ) = γ. Note that |{ k : a ik = a jk = 1 }\{ i, j }| ≥ µ ( n − − − ( n −
2) = (2 µ − n − µ since i and j share at least (2 µ − n − µ neighbors andmin ≤ θ ≤ γ [sin( θ ) + sin( γ − θ )] ≥ sin( γ ) , γ ≤ π . (4.11)Thus n (cid:88) k =1 ( a ik sin( θ i − θ k ) + a jk sin( θ k − θ j )) ≥ ((2 µ − n − µ ) sin( γ )+2 sin( γ ) ≥ (2 µ − n sin( γ ) . Case 2: If a ij = 0, then n (cid:88) k =1 ( a ik sin( θ i − θ k ) + a jk sin( θ k − θ j )) ≥ (cid:88) { k : a ik = a jk =1 }\{ i,j } (sin( θ i − θ k ) + sin( θ k − θ j )) ≥ |{ k : a ik = a jk = 1 }\{ i, j }| · sin( γ ) ≥ (2 µ ( n − − ( n − · sin( γ ) ≥ (2 µ − n sin( γ ) . θ i ( t ) − ˙ θ j ( t ) ≤ ω i − ω j − (2 µ − K sin( γ ) ≤ ω max − ω min − (2 µ − K sin( γ )which means θ i ( t ) − θ j ( t ) has a negative derivative at t if ω max − ω min < (2 µ − K sin( γ ) ⇐⇒ K sin( γ ) > ω max − ω min µ − . This implies that if the condition above holds, the range of { θ i } ni =1 will not increase, i.e., θ stays in ∆( γ ) for γ < π/ . Lemma 4.3.
Under the assumption (4.8) , if θ (0) ∈ ∆( γ ) for γ < π , then θ ( t ) satis-fies exponential frequency synchronization and the frequency synchronized solution lies in ∆( γ ) . Proof:
First note that Lemma 4.2 implies that as long as the initialization satisfies θ (0) ∈ ∆( γ ), θ ( t ) stays in ∆( γ ). Now it is safe to only analyze the dynamics in ∆( γ ) . We consider the (cid:96) -norm of ˙ θ ( t ) with initialization θ (0) inside the cohesive region.Take the derivative w.r.t. t and we havedd t (cid:107) ˙ θ ( t ) (cid:107) = 2 ˙ θ (cid:62) ¨ θ where ¨ θ i ( t ) = − Kn n (cid:88) j =1 a ij cos( θ i ( t ) − θ j ( t )) · ( ˙ θ i ( t ) − ˙ θ j ( t )) . Writing it in matrix form gives ¨ θ ( t ) = − Kn J ( θ ( t )) ˙ θ ( t )where J ( θ ( t )) satisfies J ( θ ) = ( J ( θ )) ij , ( J ( θ )) ij = − a ij cos( θ i − θ j ) ≤ − a ij cos( γ ) ≤ θ ∈ ∆( γ ) and γ < π/ . Note that the underlying network is connected since the degreeof each node is at least n/
2. Thus the second smallest eigenvalue of J (0), the Laplacianof the corresponding adjacency matrix, is strictly positive, which is a classical result inspectral graph theory [8]. Therefore, J ( θ ) is positive semidefinite for θ ∈ ∆( γ ) and itssecond smallest eigenvalue satisfies λ ( J ( θ )) ≥ cos( γ ) · λ ( J (0)) > . Note that ˙ θ ( t ) (cid:62) n = 0, and thus it holds thatdd t (cid:107) ˙ θ ( t ) (cid:107) = 2 ˙ θ ( t ) (cid:62) ¨ θ ( t ) = − Kn ˙ θ ( t ) (cid:62) J ( θ ( t )) ˙ θ ( t ) ≤ − Kn λ ( J ( θ ( t ))) · (cid:107) ˙ θ ( t ) (cid:107) ≤ − Kn cos( γ ) · λ ( J (0)) · (cid:107) ˙ θ ( t ) (cid:107) . As a result, we have dd t (cid:16) e n − K cos( γ ) λ ( J (0)) · t (cid:107) ˙ θ ( t ) (cid:107) (cid:17) ≤ .
13y integrating the time t from 0, it holds that (cid:107) ˙ θ ( t ) (cid:107) ≤ e − n − K cos( γ ) λ ( J (0)) · t (cid:107) ˙ θ (0) (cid:107) . In other words, (cid:107) ˙ θ ( t ) (cid:107) converges to 0 exponentially fast as t → ∞ where γ < π , i.e.,converging to a frequency synchronized solution of the system. Moreover, this equilibriumis stable since it is within the cohesive region and J ( θ ) (cid:23) Proof of Theorem 2.1.
The proof of the main theorem follows from combining the afore-mentioned three results, Lemma 4.1, Lemma 4.2, and Lemma 4.3 together. First of all,Lemma 4.2 implies that if (cid:107) ω (cid:107) ∞ K ≤ µ − , then the cohesive region ∆( π/
2) is a basin of attraction. Then Lemma 4.3 shows that aslong as an initialization is chosen within ∆( π/ (cid:107) ω (cid:107) ∞ K ≤ (cid:114) µ −
34 + 1 − µn + µ − . we have shown that all the stable equilibria, if there is any, must stay within ∆( π/ . As long as 1 ≥ µ > −√ , we have (cid:114) µ −
34 + µ − ≤ µ − . This means the condition in Lemma 4.1 is stronger than those in Lemma 4.2 and 4.3. Inother words, if (2.4) holds, the whole system has a unique and stable frequency synchro-nized solution which is also phase-cohesive.
References [1] J. A. Acebr´on, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler. TheKuramoto model: A simple paradigm for synchronization phenomena.
Reviews ofModern Physics , 77(1):137, 2005.[2] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou. Synchronizationin complex networks.
Physics Reports , 469(3):93–153, 2008.[3] A. S. Bandeira, N. Boumal, and V. Voroninski. On the low-rank approach forsemidefinite programs arising in synchronization and community detection. In
Con-ference on Learning Theory , pages 361–382, 2016.[4] F. Bullo.
Lectures on Network Systems . Kindle Direct Publishing, 2019. Withcontributions by J. Cortes, F. Dorfler, and S. Martinez.[5] E. A. Canale and P. Monz´on. Exotic equilibria of Harary graphs and a new minimumdegree lower bound for synchronization.
Chaos: An Interdisciplinary Journal ofNonlinear Science , 25(2):023106, 2015.146] H. Chiba, G. S. Medvedev, and M. S. Mizuhara. Bifurcations in the Kuramoto modelon graphs.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 28(7):073109,2018.[7] N. Chopra and M. W. Spong. On exponential synchronization of kuramoto oscilla-tors.
IEEE transactions on Automatic Control , 54(2):353–357, 2009.[8] F. R. Chung.
Spectral Graph Theory , volume 92. American Mathematical Society,1997.[9] F. D¨orfler and F. Bullo. On the critical coupling for Kuramoto oscillators.
SIAMJournal on Applied Dynamical Systems , 10(3):1070–1099, 2011.[10] F. D¨orfler and F. Bullo. Exploring synchronization in complex oscillator networks. In , pages 7157–7170.IEEE, 2012.[11] F. Dorfler and F. Bullo. Synchronization and transient stability in power networksand nonuniform Kuramoto oscillators.
SIAM Journal on Control and Optimization ,50(3):1616–1642, 2012.[12] F. D¨orfler and F. Bullo. Synchronization in complex networks of phase oscillators:A survey.
Automatica , 50(6):1539–1564, 2014.[13] F. D¨orfler, M. Chertkov, and F. Bullo. Synchronization in complex oscillatornetworks and smart grids.
Proceedings of the National Academy of Sciences ,110(6):2005–2010, 2013.[14] T. Ichinomiya. Frequency synchronization in a random oscillator network.
PhysicalReview E , 70(2):026116, 2004.[15] A. Jadbabaie, N. Motee, and M. Barahona. On the stability of the Kuramoto modelof coupled nonlinear oscillators. In
Proceedings of the 2004 American Control Con-ference , volume 5, pages 4296–4301. IEEE, 2004.[16] Y. Kuramoto. Self-entrainment of a population of coupled non-linear oscillators. In
International Symposium on Mathematical Problems in Theoretical Physics , pages420–422. Springer, 1975.[17] Y. Kuramoto.
Chemical Oscillations, Waves, and Turbulence , volume 19. SpringerScience+Business Media, 1984.[18] S. Ling, R. Xu, and A. S. Bandeira. On the landscape of synchronization net-works: A perspective from nonconvex optimization.
SIAM Journal on Optimization ,29(3):1879–1907, 2019.[19] M. Lopes, E. Lopes, S. Yoon, J. Mendes, and A. Goltsev. Synchronizationin the random-field Kuramoto model on complex networks.
Physical Review E ,94(1):012308, 2016.[20] J. Lu and S. Steinerberger. Synchronization of Kuramoto oscillators in dense net-works. arXiv preprint arXiv:1911.12336 , 2019.[21] J. Markdahl, J. Thunberg, and J. Gon¸calves. Almost global consensus on the n-sphere.
IEEE Transactions on Automatic Control , 63(6):1664–1675, 2017.1522] J. Markdahl, J. Thunberg, and J. Goncalves. High-dimensional Kuramoto modelson Stiefel manifolds synchronize complex networks almost globally.
Automatica ,113:108736, 2020.[23] R. E. Mirollo and S. H. Strogatz. Synchronization of pulse-coupled biological oscil-lators.
SIAM Journal on Applied Mathematics , 50(6):1645–1662, 1990.[24] J. Pe˜na Ramirez, L. A. Olvera, H. Nijmeijer, and J. Alvarez. The sympathy of twopendulum clocks: beyond Huygens’ observations.
Scientific Reports , 6(23580), 032016.[25] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths. The Kuramoto model incomplex networks.
Physics Reports , 610:1–98, 2016.[26] S. Strogatz.
Sync: The Emerging Science of Spontaneous Order . Penguin UK, 2004.[27] S. H. Strogatz. From Kuramoto to Crawford: exploring the onset of synchronizationin populations of coupled oscillators.
Physica D: Nonlinear Phenomena , 143(1-4):1–20, 2000.[28] S. H. Strogatz. Exploring complex networks.
Nature , 410(6825):268, 2001.[29] R. Taylor. There is no non-zero stable fixed point for dense networks in the ho-mogeneous Kuramoto model.
Journal of Physics A: Mathematical and Theoretical ,45(5):055102, 2012.[30] A. Townsend, M. Stillman, and S. H. Strogatz. Circulant networks of identicalKuramoto oscillators: Seeking dense networks that do not globally synchronize andsparse ones that do. arXiv preprint arXiv:1906.10627 , 2019.[31] R. van der Hofstad.
Random Graphs and Complex Networks , volume 1 of
CambridgeSeries in Statistical and Probabilistic Mathematics . Cambridge University Press,2016.[32] M. Verwoerd and O. Mason. Global phase-locking in finite populations of phase-coupled oscillators.
SIAM Journal on Applied Dynamical Systems , 7(1):134–160,2008.[33] M. Verwoerd and O. Mason. On computing the critical coupling coefficient forthe Kuramoto model on a complete bipartite graph.
SIAM Journal on AppliedDynamical Systems , 8(1):417–453, 2009.[34] D. A. Wiley, S. H. Strogatz, and M. Girvan. The size of the sync basin.
Chaos: AnInterdisciplinary Journal of Nonlinear Science , 16(1):015103, 2006.[35] A. T. Winfree. Biological rhythms and the behavior of populations of coupled oscil-lators.