On the number of stable solutions in the Kuramoto model

 View Online  Export Citation RESEARCH ARTICLE | SEPTEMBER 20 2023 On the number of stable solutions in the Kuramoto model Special Collection: Nonlinear dynamics, synchronization and networks: Dedicated to Jürgen Kurths' 70th birthday Alex Arenas  ; Antonio Garijo ; Sergio Gómez ; Jordi Villadelprat Chaos 33, 093127 (2023) https://doi.org/10.1063/5.0161977 27 August 2024 08:12:10 Chaos ARTICLE pubs.aip.org/aip/cha On the number of stable solutions in the Kuramoto model Cite as: Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 Submitted: 13 June 2023 ·Accepted: 25 August 2023 · Published Online: 20 September 2023 View Online Export Citation CrossMark Alex Arenas,1,2,a)Antonio Garijo,1Sergio Gómez,1and Jordi Villadelprat1 AFFILIATIONS 1Departament d’Enginyeria Informàtica i Matemàtiques, Universitat Rovira i Virgili, 43007 Tarragona, Spain 2Pacific Northwest National Laboratory, 902 Battelle Blvd, Richland, Washington 99354, USA Note: This paper is part of the Focus Issue on Nonlinear dynamics, synchronization and networks: Dedicated to Juergen Kurths’ 70th birthday. a)Author to whom correspondence should be addressed: [email protected] ABSTRACT We consider a system of ncoupled oscillators described by the Kuramoto model with the dynamics given by ˙ θ=ω+Kf(θ). In this system, an equilibrium solution θ∗is considered stable when ω+Kf(θ∗)=0, and the Jacobian matrix Df(θ∗)has a simple eigenvalue of zero, indicating the presence of a direction in which the oscillators can adjust their phases. Additionally, the remaining eigenvalues of Df(θ∗)are negative, indicating stability in orthogonal directions. A crucial constraint imposed on the equilibrium solution is that |0(θ∗)| ≤ π, where |0(θ∗)| represents the length of the shortest arc on the unit circle that contains the equilibrium solution θ∗. We provide a proof that there exists a unique solution satisfying the aforementioned stability criteria. This analysis enhances our understanding of the stability and uniqueness of these solutions, offering valuable insights into the dynamics of coupled oscillators in this system. © 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0161977 Synchronization in ensembles of network-coupled heterogeneous oscillators is crucial in various natural and engineered phenomena, ranging from cell cycles to robust power systems. One of the most prominent and elegant models for studying synchronization is the Kuramoto model,1which provides a mathematically tractable description of this phenomenon. Kuramoto recognized the mean-field approach as the most suitable method for analytical treatment and introduced an all-to-all purely sinusoidal coupling scheme, deriving the governing equations for each oscillator in the system. The Kuramoto model is a mathematically tractable description of synchronization in network-coupled heterogeneous oscillators. We investigate the conditions for stable equilibrium solutions in this model. Our main finding is that a stable equilibrium solution θ∗exists when ω+Kf(θ∗)=0, and the Jacobian matrix Df(θ∗)has zero as a simple eigenvalue and negative eigenvalues in orthogonal directions. We also establish the constraint |0(θ∗)| ≤ π, indicating that the equilibrium lies within the shortest arc on the unit circle containing it. This analysis contributes to understanding the dynamics and stability of the Kuramoto model. I. INTRODUCTION The Kuramoto model is a widely studied mathematical model that describes the synchronization behavior in a system of ncoupled oscillators. It was introduced by Kuramoto in his seminal works1,2 and has since become a fundamental framework for understanding synchronization phenomena in various fields. The model assumes n connected oscillators, each characterized by a phase variable θiand a natural frequency ωi. The dynamics of each i-oscillator, coupled with the rest, is described by the equation ˙ θi=ωi+K n n X j=1 sin(θj−θi), (1) where ˙ θidenotes the time derivative of θiand Krepresents the coupling strength among the oscillators. The Kuramoto model has been extensively studied due to its ability to capture and explain synchronization phenomena in a wide range of systems, including biological, physical, and social systems.3–6It has been applied to understand phenomena, such Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 33, 093127-1 © Author(s) 2023 27 August 2024 08:12:10 Chaos ARTICLE pubs.aip.org/aip/cha as neuronal synchronization,7power grid synchronization,8and opinion formation in social networks.9 In this paper, we focus on investigating the stability and properties of equilibrium solutions in the Kuramoto model for a finite system of noscillators. We analyze the conditions under which stable synchronization emerges and explore the dynamics of the system as the coupling strength Kvaries. Our findings contribute to a deeper understanding of synchronization mechanisms in finite oscillator networks. II. PRELIMINARIES AND MAIN RESULT In this work, we focus on the finite Kuramoto model, which consists of noscillators arranged on the unit circle. The phase space for this model is the n-torus, defined as Tn=S1×...×S1. (2) We use the usual model in each circle S1=R/2πZ, identifying two angles θ1and θ2in Rif and only if θ1−θ2=2πkfor some index k∈Z. Thus, given any point θ∈Tn, we can always express θ=(θ1,...,θn), where θi∈(−π,π]. This phase space serves as the domain for a system of first-order ordinary differential equations, initially introduced by Kuramoto,1          ˙ θ1=ω1+Kf1(θ1,...,θn), ˙ θ2=ω2+Kf2(θ1,...,θn), . . .. . . ˙ θn=ωn+Kfn(θ1,...,θn), (3) where fi(θ1,...,θn)= n X j=1,j6=i aij sin (θj−θi)for 1 ≤i≤n. (4) The network’s topology is represented by an adjacency matrix Aof size n×n. Specifically, aii =0, and aij =aji =1 if nodes i and jare connected, while aij =aji =0 otherwise. Furthermore, we require the network to be connected. Systems (3) and (4) described above are known as the “Kuramoto model” of synchronization associated with the network defined by the adjacency matrix A. However, we refer to the specific case where aij =1 for all i6= jand aii =0 as the “classic Kuramoto model” since it corresponds to the original formulation. Using the compact notation ω=(ω1,...,ωn),θ=(θ1,...,θn), and f=(f1,...,fn). The Kuramoto model can be expressed as ˙ θ=ω+Kf(θ), (5) where K≥0 is a real parameter. For any fixed initial condition θ0 on the torus, there exists a unique solution ϕ: [0, ∞)→Tnof (5), defined for all t≥0, passing through θ0at t=0. An θ∗is referred to as an equilibrium point of (5) if ω+Kf(θ∗)=0. In this case, the solution ϕ(t)≡θ∗, defined for all t≥0, satisfies the ordinary differential equation (ODE), and the local behavior of this solution is determined by the spectrum of Df(θ∗), i.e., the eigenvalues of the differential matrix. We denote the differential matrix as H:=Df. Thus, H(θ)=     h1,1(θ)h1,2(θ) . . . h1,n(θ) h1,2(θ)h2,2(θ) . . . h2,n(θ) . . .. . ..... . . h1,n(θ) . . . hn−1,n(θ)hn,n(θ)      , (6) where hi,i(θ)=∂fi ∂θi (θ)= − n X j=1,j6=i aij cos (θj−θi), hi,j(θ)=∂fi ∂θj (θ)=aij cos (θj−θi)i6= j. (7) We observe that the matrix H(θ)is symmetric, which implies that its eigenvalues are real. This property allows us to compute the derivative of the eigenvalues of H(θ)and control how they evolve (as demonstrated in Proposition III.3). Notably, the spectrum of H(θ) remains unaffected by variations in ωor K. Before delving into the analysis of the spectrum of H(θ), we present an alternative approach to this problem. The symmetry of the matrix Henables us to find a primitive of f. More precisely, we can consider the scalar map V: [−π,π]n⊂Rn→R given by V(θ)=C−ω·θ−K n X i=1 n X j=i+1 aij cos (θi−θj), (8) where the expression a·bdenotes the inner product of two vectors aand bin Rnand Cis any constant value. It is important to note that the function Vis defined within the hypercube [−π,π]nbut not on the torus Tn, except in the case when ω=0. This is due to the fact that the inner product ω·θis not a 2π-periodic map. To illustrate this claim, let us consider the case where ω6= 0, and without loss of generality, assume that ω16= 0. We observe that the inner product ω·(π, 0, ..., 0)=πω1, whereas ω·(−π, 0, ..., 0)= −πω1. This demonstrates that the inner product does not exhibit 2πperiodicity, leading to the conclusion that V is not well-defined on the torus Tnin the presence of non-zero ω. Using the expression for Vgiven above, it is evident that the dynamics of the Kuramoto model (5) can be represented as the potential flow generated by the scalar map V(8). This can be seen by observing that ˙ θ=ω+Kf(θ)= −∇V(θ). (9) Therefore, locally, the flow of the Kuramoto model can be interpreted as the potential flow induced by the scalar map V. The Hessian matrix of V, denoted as HV, and the differential of f, denoted as Df, are related by HV(θ)= −KDf(θ)= −KH(θ). (10) It is worth noting that the symmetry of H(θ)is not surprising, as it is a multiple of the Hessian matrix of V. This relationship establishes the connection between the symmetry of H(θ)and the potential flow representation of the Kuramoto model. Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 33, 093127-2 © Author(s) 2023 27 August 2024 08:12:10 Chaos ARTICLE pubs.aip.org/aip/cha The Kuramoto model exhibits a finite number of equilibrium points. This can be demonstrated by rewriting the nonlinear system ω+Kf(θ)=0as a quadratic system and applying Bézout’s theorem (for more details, refer to Ref. 10). In previous discussions, we have defined an equilibrium solution of the Kuramoto system (5). Now, we introduce the concept of a stable equilibrium solution, as defined in Ref. 11. Definition II.1: We say that θ∗=θ∗ 1,...,θ∗ nis a stable solution of (5) if and only if (a) ω+Kf(θ∗)=0, (b) H(θ∗)is negative semi-definite, and (c) dim(Ker(H(θ∗))) =1. One important remark about the above definition is that the requirement for negative definiteness of H(θ∗)is not feasible in the Kuramoto model. Due to the nature of the equation driving system (5), stable solutions are not isolated. Specifically, if θ∗is an equilibrium point, then for any angle α∈S1,θ∗+α=(θ∗ 1+α,...,θ∗ n +α) is also an equilibrium point. We denote all these solutions as [θ∗]=θ∗+α,α∈S1. This fact is evident in the eigenvalues of H(θ), as 1=(1, ..., 1)is an eigenvector of H(θ)with eigenvalue λ=0 for all θ∈Tn. Thus, λ=0 always appears in the spectrum of H(θ). The third condition requires that the remaining eigenvalues of H(θ∗)are strictly negative. A second important remark concerns the stability of the equilibrium point θ∗. On one hand, the existence of a strictly positive eigenvalue of Df(θ∗)implies that the unstable manifold of θ∗has a dimension greater than or equal to one. On the other hand, the stability of θ∗can also be linked to the local behavior of Vnear θ∗. Thus, we impose that Vexhibits a local minimum at θ∗. Consequently, the Hessian map HV(θ∗)is positive semi-definite, and, therefore, H(θ∗)is negative semi-definite since HV(θ∗)= −KH(θ∗) [see Eq. (10)]. The following definition will play a fundamental role in the classification of stable solutions of the Kuramoto model. Definition II.2: For any point θ=(θ1,...,θn)on Tn, we denote 0(θ)as a closed shortest arc on S1containing all angles. The length of this arc is denoted by |0(θ)|(see Fig. 1). The literature on the Kuramoto model is extensive, and it has served as the foundation for the study of various synchronization phenomena. Providing a comprehensive list of all contributions would be impractical. However, we can outline some key findings in the finite Kuramoto model based on different scenarios involving the frequency vector ω. FIG. 1. Three different examples of θ=(θ1,...,θ7)∈T7. We show the closed shortest arc 0(θ) containing all angles and their length |0(θ)|. The first case corresponds to ω=0. Taylor12 demonstrated that the origin θ∗=(0, ..., 0)is the only stable solution in the classic Kuramoto model. Subsequently, several authors12–14 showed that for sufficiently dense networks, the origin is the unique stable solution. Network density is measured using a parameter µ, which indicates that each oscillator has at least µ(n−1)connections with other oscillators. The classic Kuramoto model corresponds to µ=1. Taylor12 proved that the origin is the unique stable solution for networks with density parameter µ≥0.9395. Ling et al.13 established the same result for µ≥(3−√2)/2 ∼0.7929 and Kassabov et al.14 for µ≥0.75. Therefore, it is possible to have multiple stable solutions for small values of µ.15 When ω=0, the Kuramoto model can be analyzed using Morse theory developed by Milnor.16 In this framework, the system of ODEs (5) represents the downhill flow of the map V:Tn→R(8). However, the global behavior becomes more complex, as Morse theory reveals the existence of multiple unstable solutions. The number of these unstable solutions is related to the Betti numbers of Tn.10 For instance, since V:Tn→Ris a continuous map defined on a compact manifold, it must reach at least one local maximum. Consequently, there always exist initial conditions that do not converge to the stable solution θ∗=0. Additionally, it should be noted that in the case of ω=0, the Kuramoto model is independent of parameter K. In the case where ω6= 0, parameter Kplays a crucial role in the dynamics of the Kuramoto model (5). Numerical simulations reveal the existence of a critical parameter Kc>0, such that for 0<K<Kc, the Kuramoto model does not possess any stable solutions. This fact can be easily demonstrated. Let ω=(ω1,...,ωn) 6= 0, and assume without loss of generality that ω16= 0. Then, there exists 0such that for 0 <K< 0, we have ω1+Kf1(θ)6= 0, since f1is a bounded map. Therefore, the Kuramoto model (5) does not exhibit any stable solutions for 0 <K< 0. Numerous results have been obtained regarding the estimation of the critical parameter Kc (see Refs. 3and 17, and references therein). On the other hand, when parameter Kis sufficiently large, the Kuramoto model becomes similar to the case when ω=0. Thus, for large enough K, the Kuramoto model possesses a stable equilibrium θKthat converges to 0as Ktends to infinity. However, the number of stable equilibrium points in the Kuramoto model is still unknown. The stability of solutions in the Kuramoto model has been extensively studied by various authors, including Refs. 11,18, and 19. The objective of this work is to investigate the number of stable solutions in the Kuramoto model. Our main result can be stated as follows: Theorem A: Let θ∗and η∗be two stable solutions of the Kuramoto model (5) satisfying |0(θ∗)| ≤ πand |0(η∗)| ≤ π. Then, [θ∗]=[η∗]. It is important to note that the above result does not impose any explicit assumptions on the system variables, such as the frequencies ω, the parameter K, or the adjacency matrix A. The Kuramoto model is known to exhibit the possibility of multiple stable solutions (see Ref. 15, and references therein). In Sec. III, we provide a concrete example where the Kuramoto model demonstrates two stable solutions. Furthermore, we discuss how this example relates to the tools utilized in the proof of Theorem A. Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 33, 093127-3 © Author(s) 2023 27 August 2024 08:12:10 Chaos ARTICLE pubs.aip.org/aip/cha III. PROOF OF THEOREM A AND EXAMPLES The following lemma is well known and can be proven using the Gershgorin circle theorem (see Chap. 6 of Ref. 20). We include it here for completeness. Lemma III.1: Let M =(mi,j)be an n ×n real, symmetric and diagonal dominant matrix, i.e., |mi,i| ≥ n X j=1, j6=i|mi,j|for 1≤i≤n. The following conditions hold: •Assume that all the entries in the diagonal verify mi,i≥0. Then, all the eigenvalues of M are non-negative, and, therefore, M is positive semi-definite and xTMx≥0for all vectors x∈Rn. •Assume that all the entries in the diagonal verify mi,i≤0. Then, all eigenvalues of M are non-positive, and, therefore, M is negative semi-definite and xTMx≤0for all vectors x∈Rn. In our work, the derivative of an eigenvalue with respect to a real parameter will play a fundamental role. The following result can be used to compute the derivative of a simple or multiple eigenvalue with respect to a real parameter. It is a combination of two theorems: Theorem 5 of Ref. 21 for the derivative of a simple eigenvalue and Theorem 2.3 of Ref. 22 for the derivative of a multiple eigenvalue depending on a single real parameter. Additional references related to this result include Refs. 22 and 23. Lemma III.2: Let A(t)be a n ×n real and symmetric matrix and such that t 7→ A(t)is a real analytic function of all t ∈(a,b) ⊂R. Suppose that λ(t0)is an eigenvalue of A(t0)with multiplicity r≥1, where t0∈(a,b). Then, there exists  > 0, and real analytic functions λ1(t),...,λr(t)and v1(t),...,vr(t), such that A(t)vs(t)=λs(t)vs(t),∀t∈(t0−,t0+), λs(t0)=λ(t0),s=1, ...,r, vs(t)Tvs(t)=1, (11) and we have that λ0 s(t0)=dλs dt (t0)=vs(t0)TA0(t0)vs(t0),s=1, ...,r. Using the above result, we can prove that, under certain conditions, the derivative of the eigenvalues of the matrix H(tθ)is not negative. In the following lemma, we collect this result. Proposition III.3: Let θ=(θ1,...,θn)be any point in Tn(2). For any t ≥0, we define the n ×n matrix A(t)=H(tθ). Let λ(t0)be an eigenvalue of A(t0). Then, λ0(t0)≥0for all 0<t0≤1/2. Moreover, assuming that θi∈[0, π]for all 1≤i≤n, then λ0(t0)≥0for all 0<t0≤1. Proof. We consider the n×nsymmetric matrix given by A(t)=H(tθ). From the definition of H(θ)[see Eqs. (6) and (7)], we have that A(t)=(hi,j(tθ)) for 1 ≤i,j≤n, where the functions hi,j(tθ)are given by hi,i(tθ)= − n X j=1,j6=i aij cos tθj−θi, hi,j(tθ)=aij cos tθj−θi. Computing the derivative with respect to the real parameter t, we obtain dhi,i(tθ) dt = n X j=1,j6=i aij(θj−θi)sin (t(θj−θi)), dhi,j(tθ) dt = −aij(θj−θi)sin (t(θj−θi)). (12) Thus, we have obtained an explicit expression of the derivative A0(t)given by A0(t)=             dh1,1(tθ) dt dh1,2(tθ) dt ... dh1,n(tθ) dt dh2,1(tθ) dt dh2,2(tθ) dt ... dh2,n(tθ) dt . . .. . ..... . . dhn,1(tθ) dt dhn,2(tθ) dt ... dhn,n(tθ) dt             . (13) We introduce the auxiliary function gt(x)=xsin(tx)for t>0. We left to the reader to check the following properties of the map gt (see Fig. 2). For any value of t>0, the function gtis an even map, i.e., gt(x)=gt(−x)for all x∈Rand gt(x)≥0 for all x∈[−π/t,π/t]. Moreover, if we pick any value of t∈(0, 1/2], then gt(x)≥0 for all x∈[−2π, 2π]. For any t>0, the matrix A0(t)is given by A0(t)=       Pn j=1,j6=1a1jgt(θj−θ1)−a12gt(θ2−θ1) . . . −a1ngt(θn−θ1) −a21gt(θ1−θ2)Pn j=1,j6=2a2jgt(θj−θ2) . . . −a2ngt(θn−θ2) . . .. . ..... . . −an1gt(θ1−θn)−an2gt(θ2−θn) . . . Pn j=1,j6=nanjgt(θj−θn)        . (14) We claim that A0(t)is a symmetric, diagonally dominant, and positive semi-definite matrix for all t∈(0, 1/2]. The symmetry of A0(t) follows from the evenness of the function gt(x)and the symmetry of the adjacency matrix A=(aij). By hypothesis, θ=(θ1,...,θn)is a point in the n-torus Tn, where −π < θi≤πfor all 1 ≤i≤n. Consequently, −2π≤θi −θj≤2πfor all 1 ≤i,j≤n. Therefore, for any t∈(0, 1/2], we have gt(θi−θj)≥0 for all 1 ≤i,j≤n. It follows that for any Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 33, 093127-4 © Author(s) 2023 27 August 2024 08:12:10 Chaos ARTICLE pubs.aip.org/aip/cha FIG. 2. Graph of gt(x)=xsin(tx)for t=1/2 (blue) and t=1 (red). 1≤k≤n,  n X j=1,j6=k akjgt(θj−θk)= n X j=1,j6=k akjgt(θj−θk) = n X j=1,j6=k−akjgt(θj−θk). This proves that matrix A0(t)is diagonally dominant for any t∈(0, 1/2]. Finally, matrix A0(t)is positive semi-definite for any t∈(0, 1/2] since it is diagonally dominant and all entries in its diagonal are non-negative (Lemma III.1). We select 0 <t0≤1/2 and consider λ(t0)as a real eigenvalue of matrix A(t0). Since A(t0)is a symmetric matrix, all its eigenvalues are real. Applying Lemma III.2, we can compute its derivative as λ0(t0)=v(t0)TA0(t0)v(t0), where v(t0)is a normalized eigenvector of A(t0)corresponding to the eigenvalue λ(t0). By our previous argument that A0(t)is positive semi-definite for all t∈(0, 1/2], we conclude that ztA0(t)z≥0 for all z∈Rn. Now, let us consider the case where θ=(θ1,...,θn)satisfies 0≤θi≤πfor all 1 ≤i≤n. From this assumption, we can easily conclude that −π≤θi−θj≤πfor all 1 ≤i,j≤n. As a result, gt(θi−θj)≥0 for all t∈(0, 1] where gt(x)=xsin(tx)is the auxiliary function defined earlier (see Fig. 2). Consequently, matrix A0(t) is symmetric, diagonally dominant, and positive semi-definite for all t∈(0, 1]. Finally, let λ(t0)be any eigenvalue of A(t0). By applying Lemma III.2, we have λ0(t0)=v(t0)T,A0(t0),v(t0)≥0, since A0(t0) is a positive semi-definite matrix for any t0∈(0, 1].  The previous proposition can be interpreted geometrically in the following sense. At t=0, and independently on θ, the n×n matrix A(0)=H(0)is given by A(0)=       −Pn j=1,j6=1a1ja12 ... a1n a21 −Pn j=1,j6=2a2j... a2n . . .. . ..... . . an1... an2−Pn j=1,j6=nanj        . (15) It is easy to see that A(0)is a symmetric, diagonally dominant, and semi-definite negative matrix. So, all its eigenvalues are nonpositive (see Lemma III.1). As we travel through the ray tθfor t∈(0, 1/2], the spectrum of A(tθ)moves to the right with respect to the spectrum of A(0). Remark 1: The Laplacian matrix is a matrix representation of a network. In particular, the rank of the Laplacian matrix is related to the number of connected components of the network (see Chap. 13 of Ref. 24 for details). Let L be the Laplacian matrix associated to the network of oscillators. From the expression of A(0)(15), we just observe that L = −A(0). Moreover, it is well known that λ=0is an eigenvalue of L whose multiplicity coincides with the number of connected components of the graph (Lemma 13.1.1 of Ref. 24). In our case, we have assumed that our network of oscillators form a connected graph. So, we conclude that λ=0is a simple eigenvalue of A(0)and the rest of its eigenvalues are strictly negative real numbers. Proof of Theorem A We define set C, related with the set of points where function V is convex, C= {θ∈[−π,π]n|H(θ)is semi-definite negative and dim(Ker(H(θ))) =1}. (16) We start showing that Cis an open and nonempty set. We first prove that 0=(0, ..., 0)belongs to C. Taking t=0 the matrix H(0)=A(0)(15) is a symmetric, diagonally dominant, and semidefinite negative matrix. Moreover, λ=0 is a simple eigenvalue of H(0)since the network formed by all the oscillators is connected (see Remark 1). We claim that Sis an open set. Let θ0be a point in C. We denote by pθ0(x)the characteristic polynomial of H(θ0). By hypothesis pθ0(x)=x·qθ0(x)with qθ0(0)6= 0. Moreover, all the roots of qθ0(x)are real and strictly negative numbers. Hence, in a sufficiently small neighborhood of θ0, the roots of qθ(x)are still strictly negative, showing that Cis an open set. We denote by ∂Cand Cthe boundary and the closure of C, respectively. Set ∂Ccontains all the θ’s such that H(θ)is negative semi-definite and λ=0 is a multiple eigenvalue. Furthermore, the set Ccoincides with the set of points where the map Vis a convex map. Finally, the open set [−π,π]n\Ccontains all the θ’s such that H(θ)has at least one strictly positive eigenvalue. As we mention before C6= ∅ is an open set. Hence, we can decompose Cinto its disjoint connected components, and there are at most a countably many connected components of C. So, C=∞ [ i=1 Ciwith Ci∩Cj= ∅for i6= j. Without loss of generality, we can assume that C1is the connected component of Ccontaining the origin 0.A priori we do not know how many of those Cifor i>1 are different from the Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 33, 093127-5 © Author(s) 2023 27 August 2024 08:12:10 Chaos ARTICLE pubs.aip.org/aip/cha empty set. We assume that θ∗=(θ1,...,θn)and η∗=(η1,...,ηn) are two stable solutions of the Kuramoto model (5) with |0(θ∗)| ≤πand |0(η∗)| ≤ π. As we mention before θ∗+(α,...,α) and η∗+(β,...,β) are also stable solutions for all α∈(−π,π] and β∈(−π,π]. We select α0and β0such that 0 ≤θi+α0≤π and 0 ≤ηi+β0≤πfor all 1 ≤i≤n. So, we just choose a stable solution in the upper half part of the n-torus. Thus, renaming θ∗=(θ1,...,θn)and η∗=(η1,...,ηn), if necessary, we can assume without loss of generality that θi,ηi∈[0, π] for all 1 ≤i≤n. We first assume that θ∗and η∗belong to two different connected components. Thus, one of them, for example, η∗belongs to Ci∗with i∗>1. We consider the ray tη∗for t∈[0, 1]. This ray crosses (at least) two different connected components C1and Ci∗. The first one since the origin 0is contained in C1and the second one since η∗belongs to Ci∗. Thus, we can assume that tη∗∈C1 for t∈[0, t0)and tη∗∈Ci∗for t∈(t1, 1]. Moreover, t0η∗∈∂C1and t1η∗∈∂Ci∗. We have proved that the derivative of any eigenvalue λ(t)of H(tη∗)verifies λ0(t)≥0 (see Proposition III.3). Now, suppose that the ray tη∗exits the set Cand enters [−π,π]n\Cfor some t∈(t0,t1). This is a contradiction with the fact that λ0(t)≥0 for t∈(0, 1], since at least one eigenvalue needs to be non-negative in C1, then positive in the complement of Cand then again non-negative in Ci∗. On the other hand, suppose that the ray tη∗does not exist C. This could be the case, for example, if t0=t1. We observe that in this case when t→t− 0all the eigenvalues of H(tη∗)increase and (at least) one of them collides to λ=0 since for t=t0and for t→t− 1, this eigenvalue needs to come back, so the derivative at this eigenvalue needs to be strictly negative. We further assume that θ∗and η∗belong to the same connected component. In this scenario, both minima correspond to the same point, as it is implausible to have two distinct local minima in a region where the map is convex, unless [θ∗]=[η∗]. We conclude this section with two particular examples involving n=5 oscillators as related to Theorem A. The first example illustrates a Kuramoto model with a unique stable solution θ1such that |0(θ1)|< π. The second example demonstrates a Kuramoto model with two stable solutions, θ1and θ2, satisfying |0(θ1)|< π and |0(θ2)|> π, respectively. We first consider a Kuramoto model with five oscillators, where only the central oscillator is connected to the rest of the elements in the network (see Fig. 3 on the left). The corresponding adjacency matrix is given by A=     01111 10000 10000 10000 10000      . (17) For any 0 < ε ≤1, we consider the Kuramoto model with adjacent matrix Aand the following choice of parameters: ωε= 0, √3 2, sin π 2(1−ε),−√3 2,−sin π 2(1−ε)! and K=1. FIG. 3. Two stable configurations in the Kuramoto model. On the left hand side with adjacent matrix A(17) and on the right hand side with adjacent matrix A(18). We also show the length of the shortest arc containing the five oscillators. We claim that the point θ∗ ε=0, π 3,π 2(1−ε),−π 3,−π 2(1−ε) is the unique stable equilibrium point of the Kuramoto model for all 0< ε ≤1. To verify this claim, we must ensure that ωε+f(θ∗ ε)=0 and that the matrix H(θ∗ ε)has four negative eigenvalues for all 0< ε ≤1. In this specific case, we can explicitly compute the five eigenvalues of H(θ∗ ε), denoted by λi(ε) for i=1, ..., 5. Simple computations reveal that λ1(ε) =0, λ2(ε) = −1 2,λ3(ε) = −τ(ε), λ4(ε) =1 4−3−6τ(ε) −q9−4τ(ε) +36τ(ε)2, λ5(ε) =1 4−3−6τ(ε) +q9−4τ(ε) +36τ(ε)2, where τ() =cos π 2(1−ε). Thus, there are always four negative eigenvalues for any 0 < ε ≤1. Moreover, |0(θ∗ ε)| = π(1−ε) < π, which tends to πas ε→0. We next consider a specific example where the Kuramoto model exhibits more than one stable solution, as demonstrated in Ref. 15 and other references. We study a Kuramoto model with five oscillators in which each oscillator is connected to its two nearest neighbors (see Fig. 3 on the right). The corresponding adjacency matrix is presented as follows: A=     01001 10100 01010 00101 10010      . (18) In this example, we assume that the vector of frequencies ωis zero and K=1. The Kuramoto model has two equilibrium points located at 0=(0, 0, 0, 0, 0)and θ∗=0, 2π 5,4π 5,6π 5,8π 5. These equilibrium points are stable since the matrices H(0)and H(θ∗)have 4 strictly negative eigenvalues, as shown in Ref. 15. Although this case is not covered by Theorem A since |0(θ∗)| = 8π 5> π (see Fig. 3 right), we can apply Proposition III.3 to this case since it is applicable in a general context. We consider the matrix A(t)=H(tθ∗)given by Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 33, 093127-6 © Author(s) 2023 27 August 2024 08:12:10 Chaos ARTICLE pubs.aip.org/aip/cha A(t)=                   −cos 8πt 5−cos 2πt 5cos 2πt 50 0 cos 8πt 5 cos 2πt 5−2 cos 2πt 5cos 2πt 50 0 0 cos 2πt 5−2 cos 2πt 5cos 2πt 50 0 0 cos 2πt 5−2 cos 2πt 5cos 2πt 5 cos 8πt 50 0 cos 2πt 5−cos 2πt 5−cos 8πt 5                   . (19) In this particular example, it is possible to obtain the explicit expression of the eigenvalues λk(t)for k=1, ..., 5 of A(t). More precisely, we have that λ1(t)=0, λ2(t)= −1 2√5+5cos 2πt 5, λ3(t)=1 2√5−5cos 2πt 5, λ4(t)= −3 2cos 2πt 5−cos 8πt 5+√2 4s9+5 cos 4πt 5−4 cos 6πt 5−4 cos(2πt)+4 cos 16πt 5, λ5(t)= −3 2cos 2πt 5−cos 8πt 5−√2 4s9+5 cos 4πt 5−4 cos 6πt 5−4 cos(2πt)+4 cos 16πt 5. In Fig. 4, we have to plot the graph of the five eigenvalues of A(t)=H(tθ∗)for 0 ≤t≤1. Thus, for t=0, the eigenvalues of H(0)are given by λ1(0)=0, λ2(0)=λ5(0)= −1 2√5+5 ∼ −3.618 034, and λ3(0)=λ4(0)=1 2√5−5∼ −1.381 966, proving, thus, that 0is a stable equilibrium point. Similarly, for t=1, the eigenvalues of H(θ∗)are given by λ1(1)=0, λ2(1) =λ5(1)= −1 2√5+5cos 2π 5∼ −1.118 034 and λ3(1)=λ4(1) FIG. 4. Graph of the eigenvalues λ1(t)(red), λ2(t)(dark lilac), λ3(t)(green), λ4(t)(black), and λ5(t)(dark blue) of A(t)(19) for 0 ≤t≤1. It is also shown the vertical line t=1/2. =1 2√5−5cos 2π 5∼ −0.427 051 showing that θ∗is also a stable equilibrium solution of the Kuramoto model. In Proposition III.3, we have proved that all the eigenvalues verify λ0 k(t)≥0 for 0 <t≤1/2. In Fig. 4, we also have to plot the vertical line t=1/2 and we can see that the functions λk(t)are increasing for all k=1, ..., 5 as it is proved in Proposition III.3. Moreover, it is possible to prove that λ0 4(t) < 0 and λ0 5(t) < 0 for some value of 1/2<t≤1. Although Theorem A does not apply to this particular example (since |0(θ∗)|> π), the ideas used in the proof of Theorem A can used to understand the existence of these two stable equilibrium solutions. Thus, these two stable equilibrium solutions 0and θ∗belong to two different connected components of set C[see (16)]. We recall that Cis the set of points θwhere four of the eigenvalues of H(θ)are strictly negative. We denote by C1the connected component of Ccontaining 0and C2the connected component of Ccontaining θ∗. In Fig. 4, it is shown the evolution of all the eigenvalues from t=0 to t=1. We observe that λ4is the only eigenvalue that changes their sign. More precisely, λ4(t)≥0 for t∈[t0,t1] and λ4(t)≤0 for t∈[0, t0]∪[t1, 1] (see Fig. 4). Hence, at t=t0, point t0θ∗belongs to the boundary of C1, and for t∈(t0,t1), point tθ∗is the complement of Cand at t=t1, point t1θ∗belongs to the boundary of C2and for t1<t≤1, point tθ∗belongs to C2. Thus, eigenvalue λ4(t)needs to increase to exit C1and decrease to enter into C2or in other words their derivative change their sign. Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 33, 093127-7 © Author(s) 2023 27 August 2024 08:12:10 Chaos ARTICLE pubs.aip.org/aip/cha IV. DISCUSSION In this paper, we have demonstrated that the Kuramoto model yields a unique stable solution [θ∗] satisfying |0(θ∗)| ≤ π. This implies that half of the unit circle can be selected to include all the oscillators. Such solutions can be perceived as a type of “entrained” solution, characterizing what is commonly seen as a cluster of entrained oscillators around a mutual phase. This is distinct from other solution types, such as the splay state solution depicted in Fig. 3. The stability of any equilibrium solution of (5) is captured within the symmetric matrix H(θ∗)as detailed in Eq. (6). Our proof of this principal finding hinges on controlling the derivative of the eigenvalues of the function t7→ λ(t), where λ(t)denotes an eigenvalue of the matrix A(t)=H(tθ). ACKNOWLEDGMENTS A.A. and S.G. also acknowledge support from Spanish Ministerio de Ciencia e Innovacion (No. PID2021-128005NB-C21), Generalitat de Catalunya (Nos. 2021SGR-633 and PDAD14/20/00001), and Universitat Rovira i Virgili (No. 2019PFR-URV-B2-41). A.G. and J.V. also acknowledge support from the Ministry of Science, Innovation and Universities of Spain through Grant No. MTM 2020-118281GB-C33. A.A. also acknowledges support from ICREA Academia, and the James S. McDonnell Foundation (220020325), the Joint Appointment Program at Pacific Northwest National Laboratory (PNNL). PNNL is a multi-program national laboratory operated for the U.S. Department of Energy (DOE) by Battelle Memorial Institute under Contract No. DE-AC05-76RL01830, and the European Union’s Horizon Europe Programme under the CREXDATA project, Grant Agreement No. 101092749. AUTHOR DECLARATIONS Conflict of Interest The authors have no conflicts to disclose. Author Contributions Alex Arenas: Conceptualization (lead); Formal analysis (equal); Supervision (lead); Writing – review & editing (equal). Antonio Garijo: Formal analysis (lead); Supervision (equal); Writing – original draft (lead). Sergio Gómez: Formal analysis (equal); Writing – review & editing (equal). Jordi Villadelprat: Formal analysis (equal); Supervision (equal). DATA AVAILABILITY Data sharing is not applicable to this article as no new data were created or analyzed in this study. REFERENCES 1Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, Kyoto University, Kyoto, 1975, Lecture Notes in Physics Vol. 39 (Springer, Berlin, 1975), pp. 420–422. 2Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, New York, NY, 1984). 3S. H. Strogatz, Physica D 143, 1 (2000). 4A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge University Press, 2003), p. 12. 5J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005). 6A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. 469, 93 (2008). 7M. Breakspear, S. Heitmann, and A. Daffertshofer, Front. Hum. Neurosci. 4, 190 (2010). 8G. Filatrella, A. H. Nielsen, and N. F. Pedersen, Eur. Phys. J. B 61, 485 (2008). 9H. Hong and S. H. Strogatz, Phys. Rev. E 84, 046202 (2011). 10J. Baillieul and C. I. Byrnes, IEEE Trans. Circuits Syst. 29, 724 (1982). 11J. C. Bronski, L. DeVille, and M. J. Park, Chaos 22, 033133 (2012). 12R. Taylor, J. Phys. A 45, 055102 (2012). 13S. Ling, R. Xu, and A. S. Bandeira, SIAM J. Optim. 29, 1879 (2019). 14M. Kassabov, S. H. Strogatz, and A. Townsend, Chaos 31, 073135 (2021). 15A. Townsend, M. Stillman, and S. H. Strogatz, Chaos 30, 083142 (2020). 16J. Milnor, Morse theory, Annals of Mathematics Studies No. 51 (Princeton University Press, Princeton, NJ, 1963), pp. vi+153, based on lecture notes by M. Spivak and R. Wells. 17F. Dörfler and F. Bullo, Automatica 50, 1539 (2014). 18R. E. Mirollo and S. H. Strogatz, Physica D 205, 249 (2005). 19T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti, IEEE Trans. Control Netw. Syst. 7, 302 (2020). 20R. Horn and C. Johnson, Matrix Analysis, 2nd ed. (Cambridge University Press, 2013). 21P. Lancaster, Numer. Math. 10, 377 (1964). 22J. G. Sun, Linear Algebra and Its Applications 137–138, 183–211 (1990). 23F. Rellich, Perturbation Theory of Eigenvalue Problems (Gordon and Breach Science Publishers, New York, 1969), pp. x+127, assisted by J. Berkowitz, with a preface by Jacob T. Schwartz. 24C. Godsil and G. Royle, Algebraic Graph Theory, Graduate Texts in Mathematics (Springer, 2001), p. 207. Chaos 33, 093127 (2023); doi: 10.1063/5.0161977 33, 093127-8 © Author(s) 2023 27 August 2024 08:12:10