Complexity reduction in the 3D Kuramoto model
CComplexity reduction in the 3D Kuramoto model
Ana Elisa D. Barioni and Marcus A. M. de Aguiar
Instituto de F´ısica ‘Gleb Wataghin’, Universidade Estadual de Campinas,Unicamp 13083-970, Campinas, SP, Brazil
The dynamics of large systems of coupled oscillators is a subject of increasing im-portance with prominent applications in several areas such as physics and biology.The Kuramoto model, where a set of oscillators move around a circle representingtheir phases, is a paradigm in this field, exhibiting a continuous transition betweendisordered and synchronous motion. Reinterpreting the oscillators as rotating unitvectors, the model was extended to allow vectors to move on the surface of D-dimensional spheres, with D = 2 corresponding to the original model. It was shownthat the transition to synchronous dynamics was discontinuous for odd D, raising alot of interest. Inspired by results in 2D, Ott et al [1] proposed an ansatz for densityfunction describing the oscillators and derived equations for the ansatz parameters,effectively reducing the dimensionality of the system. Here we take a different ap-proach for the 3D system and construct an ansatz based on spherical harmonicsdecomposition of the distribution function. Our result differs significantly from thatproposed in [1] and leads to similar but simpler equations determining the dynam-ics of the order parameter. We derive the phase diagram of equilibrium solutionsfor several distributions of natural frequencies and find excellent agreement withsimulations. We also compare the dynamics of the order parameter with numericalsimulations and with the previously derived equations, finding good agreement in allcases. We believe our approach can be generalized to higher dimensions and help toachieve complexity reduction in other systems of equations. a r X i v : . [ n li n . AO ] F e b I. INTRODUCTION
The study of synchronization between elements of a system has a long history that datesback to Huygens [2]. Over the past years several examples of synchronization have beenfound in different physical systems, from coupled metronomes [3] and electochemical oscil-lators [4] to neural networks [5], motivating the development of new mathematical methodsto understand their global behavior. One aspect of great importance is the transition fromdisordered to synchronized motion as the coupling intensity between the elements increases.In 1975 Yoshiki Kuramoto proposed a simple model of N coupled oscillators that couldbe solved analytically in the limit where N goes to infinity [6, 7]. The oscillators (or par-ticles) were described by phases θ i and natural frequencies ω i , selected from a symmetricdistribution g ( ω ). The oscillators interact according to the equations˙ φ i = ω i + KN N (cid:88) j =1 sin ( φ j − φ i ) . (1)where K is the coupling strength and i = 1 , ..., N . The complex order parameter z = pe iψ ≡ N N (cid:88) i =1 e iθ i (2)measures the degree of phase synchronization of the particles: disordered motion implies p ≈ p ≈
1. Kuramoto showed that the onset of sychronization couldbe described, in equilibrium, as a second order phase transition, where p remains very smallfor 0 < K < K c and increases as p = c √ K − K c for K > K c . These results opened theway to large number of modifications and generalizations of his model, including differenttypes of coupling functions [8–10], introduction of networks of connections (so that notall oscillators are connected to each other) [11, 12], different distributions of the oscillator’snatural frequencies (including frequencies proportional to the number of connections, leadingto explosive synchronization) [13, 14], inertial terms [15–17] and external periodic drivingforces [18–20]. More recently, interest has been shifted to understand oscillations in largerdimensions. It has been shown, in particular, that the natural extension of the Kuramotomodel to more dimensions exhibit first order phase transitions in odd dimensions and secondorder transitions in even dimensions [21].In terms of mathematical methods, a breakthrough insight was given by Ott and Antonsen[22], who considered the continuity equation satisfied by the oscillators and proposed anansatz for their distribution on the unit circle involving only two time dependent parameters.In the special case where the natural frequencies of the oscillators were given by a Lorentziandistribution, the differential equations for the ansatz parameters could be translated intoequations for the order parameter z , allowing the system to be studied through two simpleequations. This dimensional reduction permitted, in particular, the construction of completebifurcation diagrams when the oscillators are driven by periodic forces [18]. External drivingforces acting only on subsets of oscillators [19] were also used to probe modular structuresin neural networks [20].The possibility of describing a large system of N coupled particles by a small set oftime-dependent parameters was pointed out by Watanabe and Strogatz [23, 24] and furtherdeveloped in [25]. They showed that the trajectory of the i-th Kuramoto oscillator can bedescribed by the Moebius transformation z i ( t ) = e iσ ( t ) a ( t ) − z i (0)1 − a ∗ ( t ) z i (0) (3)where z i = e iφ i , σ is real and a complex. The Moebius transformation was generalized byTanaka [26] who extended it to multidimensional real of complex variables. For complexvariables, attractive coupling between particles (similar to the Kuramoto model) and mul-tivariate Lorentzian distribution of natural frequencies, equations of motion for the orderparameter similar to the Stuart-Landau equation could be derived. For real variables, cor-responding, for instance, to the Kuramoto model in 3 dimensions, an equation for the orderparameter could only be obtained for the case of identical oscillators.In this work we consider the 3-dimensional Kuramoto model proposed in [21]. In thismodel the phases of the 2D oscillators are re-interpreted as unit vectors rotating on thesurface of a circle and extended to unit vectors moving on the surface of the sphere. Chandraet al [21] showed that for N odd the onset of synchronization happens through a first orderphase transition. They also proposed an extension of the Ott-Antonsen ansatz to the N-dimensional system [1] by modifying slightly the original function proposed in [22]. Herewe take a different approach to construct an ansatz for the system: instead of guessing theform of the distribution of oscillators on the sphere, we expand the distribution in sphericalharmonics and make an ansatz for the expansion coefficients, similar to the constructionof the 2D case. We obtain a function that differs from the one derived in [21] and leadsto a simpler relation between the ansatz and the order parameter. We derive a differentialequation for the ansatz parameters and the phase transition diagram for three types ofdistributions of natural frequencies. Our results agree well with numerical simulations using N = 10000 oscillators.This paper is organized as follows: in section II we re-derive the reduced equations for the2D Kuramoto model using unit vectors. We show that Ott and Antonsen’s ansatz describesthe phase transition for any distribution of frequencies, not only the Lorentzian. In sectionIII we briefly review the 3D Kuramoto model, propose our ansatz and derive the differentialequations for the ansatz parameters using the continuity equation. We then consider thecase of identical oscillators and then the equilibrium solutions for general distributions ofnatural frequencies to study the behavior of the order parameter in terms of the couplingintensity, where we find a first order phase transition in dimension 3. In section IV wediscuss the differences and similarities between our reduced system and that proposed in [1]. II. 2D KURAMOTO MODEL
In this section we rewrite the original Kuramoto model in vector form and show how theOtt-Antonsen ansatz can be used to characterize the behavior of the order parameter as afunction of the coupling constant in equilibrium.
A. Vector formulation
Following Chandra et al, we describe the oscillators by unit vectors (cid:126)σ i = (cos φ i , sin φ i ) . It is easy to show that if φ i satisfies Kuramoto’s equation (1) then d (cid:126)σ i dt = W i (cid:126)σ i + KN (cid:88) j [ (cid:126)σ j − ( (cid:126)σ i · (cid:126)σ j ) (cid:126)σ i ] (4)where W i is an anti-symmetric matrix containing the natural frequency ω i : W i = − ω i ω i . (5)The complex order parameter z , Eq.(2), is replaced by the vector (cid:126)p = 1 N (cid:88) i (cid:126)σ i (6)describing the center of mass of the system.We also consider an extension of the model where the coupling constant K is replacedbe a 2 × K i with elements K ij = K ij ( φ i ) that might depend on φ i . It is convenientto define an auxiliary vector (cid:126)q i = K i (cid:126)p and rewrite the generalization of the dynamicalequations as d (cid:126)σ i dt = W i (cid:126)σ i + [ (cid:126)q i − ( (cid:126)σ i · (cid:126)q i ) (cid:126)σ i ] . (7)Norm conservation, | (cid:126)σ i | = 1, is guaranteed for any set of regular matrices K i (as can beseen by taking the scalar product of Eq.(8) with (cid:126)σ i ). Moreover, the equations of motion areinvariant under the transformation (cid:126)q i → (cid:126)q i + a i (cid:126)σ i , for any function a i ( φ i ). It is also convenientto add a third (artificial) dimension and write (cid:126)σ i = (cos φ i , sin φ i ,
0) and (cid:126)ω i = (0 , , ω i ) . Withthese definitions we obtain d (cid:126)σ i dt = (cid:126)ω i × (cid:126)σ i + [ (cid:126)q i − ( (cid:126)σ i · (cid:126)q i ) (cid:126)σ i ] . (8) B. Continuity equation
In the limit of infinitely many oscillators we define f ( ω, φ, t ) as the density of oscillatorswith natural frequency ω at position φ in time t . It satisfies (cid:90) π dφf ( ω, φ, t ) = g ( ω ) , (9)where g ( ω ) is the distribution of natural frequencies and (cid:90) π dφ (cid:90) dωf ( ω, φ, t ) = (cid:90) g ( ω ) dω = 1 . (10)Conservation of oscillator number implies the continuity equation ∂f∂t + ∂ ( f v φ ) ∂φ = 0 (11)with velocity field (cid:126)v = ω × ˆ r + (cid:126)q − (ˆ r · (cid:126)q )ˆ r = ( ω + q φ ) ˆ φ ≡ v φ ˆ φ. (12)Equation (11) can be written explicitly as ∂f∂t + ( ω + q φ ) ∂f∂φ − f q r = 0 (13)where q r = (cid:126)q · ˆ r , q φ = (cid:126)q · ˆ φ , ˆ r = (cos φ, sin φ ) and ˆ φ = ( − sin φ, cos φ ) (we shall omit theartificial third dimension whenever possible). We also note that ∂q φ ∂φ = − q r . Eq.(6) for theorder parameter becomes (cid:126)p ( t ) = (cid:90) ˆ r ( φ ) f ( ω, φ, t ) dφ dω. (14) C. Ansatz for density function
The density of oscillators is a periodic function of φ and, therefore, can be expanded inFourier series: f ( ω, φ, t ) = g ( ω )2 π ∞ (cid:88) m = −∞ a m ( ω, t ) e imφ . (15)Ott & Antonsen [22] showed that the ansatz a m = ρ | m | e − im Φ is self-consistent, in the sensethat it preserves this form at all times, remaining in this restricted subset of density functions.The ansatz is, therefore, parametrized by ρ ( ω, t ) and Φ( ω, t ): f ( ω, φ, t ) = g ( ω )2 π ∞ (cid:88) m = −∞ ρ | m | e im ( φ − Φ) . (16)Summing the geometric series we obtain f ( ω, φ, t ) = g ( ω )2 π (1 − ρ )1 + ρ − ρ cos ξ (17)where ξ = φ − Φ. It is convenient to define the vector (cid:126)ρ = ρ (cos Φ , sin Φ ,
0) so that cos ξ =cos( φ − Φ) = ˆ r · ˆ ρ .Substituting Eq.(16) into (14) we see that only the terms with m = 1 and m = − (cid:126)p ( t ) = (cid:90) g ( ω )2 π ˆ r ( φ ) ρ (cid:0) e iξ + e − iξ (cid:1) dφ dω = (cid:90) g ( ω ) π (cid:20)(cid:90) π dφ ˆ r ( φ )ˆ r T ( φ ) (cid:21) (cid:126)ρ dω where the dyadic product ˆ r ( φ )ˆ r T ( φ ) is the projector in the ˆ r direction. It is easy to checkthat the integral between brackets results π and (cid:126)p ( t ) = (cid:90) (cid:126)ρ ( ω, t ) g ( ω ) dω. (18)This is a key equation, that connects the parameters of the ansatz with the order param-eter. The specific relation depends on the form of g ( ω ). We will also use polar variables forthe order parameter, (cid:126)p = p (cos Ψ , sin Ψ ,
0) and ν = φ − Ψ. D. Equations of motion for ansatz parameters
The derivatives of f are ∂f∂φ = − ∂f∂ Φ = 2 ρ φ (1 − ρ ) D , (19) ∂f∂ρ = − ρ − (1 + ρ ) cos ξ ] D , (20)where we defined D = 1 + ρ − ρ cos ξ. We can write Eq.(13) explicitly as ∂f∂ρ ˙ ρ + ∂f∂ Φ ˙Φ + ( ω + q φ ) ∂f∂φ − f q r = 0 (21)Replacing the derivatives we find, after some simplifications presented in Appendix A, thefollowing equation for the dynamics of (cid:126)ρ ˙ (cid:126)ρ = (cid:126)ω × (cid:126)ρ + 12 (1 + ρ ) K (cid:126)p − ( K (cid:126)p · (cid:126)ρ ) (cid:126)ρ. (22)Taking the scalar product with ˆ ρ we also get˙ ρ = 12 (1 − ρ )( K (cid:126)p · ˆ ρ ) . (23) E. Equilibrium solutions
In this subsection we particularize for the case K = K . If K = 0 the only equilibriumsolution is ρ = p = 0. For K (cid:54) = 0 we can solve Eq.(22) with ˙ (cid:126)ρ = 0 for a given (cid:126)p andcheck that the solution allows the determination of (cid:126)p via Eq.(18) self-consistently. If ρ (cid:54) = 0,then either ρ = 1, and we only need to find the direction ˆ ρ , or (cid:126)ρ is perpendicular to (cid:126)p (seeEq.(23)). Fixing ˆ x as the direction of (cid:126)p , so that (cid:126)p = p ˆ x and (cid:126)p · ˆ ρ = p cos Φ, and setting ρ = 1we find for the ˆ x and ˆ y components of Eq.(22) − ω sin Φ + Kp − Kp cos Φ = 0 (24) ω cos Φ − Kp cos Φ sin Φ = 0 . (25)The solution is sin Φ = ωKp if | ω | ≤ Kp . For | ω | > Kp we set (cid:126)ρ = ρ ˆ y and get ρ = ω/Kp [1 ± (cid:112) − K p /ω ]. One can check that the stable solution is the one with the minus Figure 1. Solution (cid:126)ρ ( ω ) for ω > < ω < Kp (cid:126)ρ stays on the unit circle, starting atΦ = 0 and moving towards Φ = π/
2. For ω > Kp (cid:126)ρ points in the y-direction and decreases to zeroas ω → ∞ . For ω < sign. The final solution, therefore is (cid:126)ρ ( ω ) = ωKp (cid:113) K p ω − x + ˆ y for | ω | ≤ Kp (cid:20) − (cid:113) − K p ω (cid:21) ˆ y for | ω | > Kp. (26)Figure 1 illustrates the dependence of (cid:126)ρ on ω .The order parameter can now be computed for any distribution function g ( ω ). Assuming g ( ω ) to be symmetric around ω = 0 we obtain p = Kp (cid:90) π/ − π/ cos (Φ) g ( Kp sin Φ) d Φ (27)which is the result found by Kuramoto. Notice that this was obtained from the ansatz,not from the exact equations. The trivial solution is p = 0. The non-trivial solution isobtained canceling p on each side and solving for p . The critical value of K , where thenon-trivial solution starts, is obtained by setting p = 0 in the resulting equation, leading to K c = 2 /πg (0). Expanding g around p = 0 to second order we find p = (cid:115) πK g (cid:48)(cid:48) (0) (1 − K/K c ) , (28)which is exactly what was found by Kuramoto.The complete p × K phase diagram can be constructed for specific distributions g ( ω ) [27].Noting that the right hand side of Eq.(27) is a function F ( Kp ) the parametric plot p = p ( K ) p K ( a ) ( b ) p K Figure 2. Order parameter as a function of the coupling constant for the Gaussian (a) andLorentzian (b) distribution of natural frequencies. Lines show the analytical result and pointsthe result of simulations with N = 1000 oscillators. can be obtained as ( K, p ) = ( x/F ( x ) , F ( x )) for x ∈ (0 , ∞ ). For Gaussian distribution withunit variance we get F ( x ) = (cid:114) π xe − x / [ I ( x /
4) + I ( x / I n ( x ) is the modified Bessel function of the first kind. For the Lorentzian distributionwith unit width we get F ( x ) = 1 x (cid:104) √ x − (cid:105) . (30)In both cases K c can be obtained as lim x → x/F ( x ) which results K c = (cid:112) /π and K c = 2for the Gaussian and Lorentzian distributions respectively. Figure 2 compares these resultswith numerical simulations. We note that numerically solving the pair of Eqs.(8) is muchfaster than the original Eq.(1). III. 3D KURAMOTO MODEL
The 3D model is a direct extension of the vector 2D model and represents unit vectors (cid:126)σ i rotating on the surface of a sphere. The equations are [21] d (cid:126)σ i dt = W i (cid:126)σ i + 1 N (cid:88) j [ K i (cid:126)σ j − ( (cid:126)σ i · K i (cid:126)σ j ) (cid:126)σ i ] (31)0where W i is a 3 × W i = − ω i ω i ω i − ω i − ω i ω i (32)and K i is the 3 × G ( (cid:126)ω ). The orderparameter is defined as in the 2D case as the system’s center of mass (cid:126)p = 1 N (cid:88) i (cid:126)σ i . (33)Defining the 3D vector (cid:126)w i = ( ω i , ω i , ω i ), renaming (cid:126)σ ≡ ˆ r and defining K i (cid:126)p ≡ (cid:126)q i (34)we obtain the dynamical equations˙ˆ r i = (cid:126)w i × ˆ r i + (cid:126)q i − (ˆ r i · (cid:126)q i )ˆ r i . (35)As in the 2D case, the equations of motion are invariant under the transformation (cid:126)q i → (cid:126)q i + a i ˆ r i and we shall use this freedom later. A. Continuity equation
The continuity equation on the sphere is ∂f∂t + 1sin θ ∂∂θ ( f sin θv θ ) + 1sin θ ∂∂φ ( f v φ ) = 0 (36)with velocity field given by Eq. (35): (cid:126)v = (cid:126)ω × ˆ r + (cid:126)q − (ˆ r · (cid:126)q )ˆ r. (37)Computing the partial derivatives the continuity equation takes the form ∂f∂t − f q r + ( ω φ + q θ ) ∂f∂θ + 1sin θ ( − ω θ + q φ ) ∂f∂φ = 0 (38)and the order parameter becomes: (cid:126)p ( t ) = (cid:90) ˆ r ( θ (cid:48) , φ (cid:48) ) f ( θ (cid:48) , φ (cid:48) , (cid:126)ω, t ) d Ω (cid:48) d ω (39)1 B. Ansatz for the density function
In 3D we can expand the density function in spherical harmonics: f ( (cid:126)ω, θ, φ, t ) = G ( (cid:126)ω ) ∞ (cid:88) l =0 l (cid:88) m = − l f lm Y lm ( θ, φ ) . (40)Now we propose an ansatz analogous to the 2D case, but with three parameters ρ ,Θ, Φ, with the form f lm = ρ l Y ∗ lm (Θ , Φ). The parameters define a vector (cid:126)ρ ( (cid:126)ω, t ) = ρ (sin Θ cos Φ , sin Θ sin Φ , cos Θ). We write f ( (cid:126)ω, θ, φ, t ) = G ( (cid:126)ω )4 π (cid:34) π ∞ (cid:88) l =1 l (cid:88) m = − l ρ l Y ∗ lm (Θ , Φ) Y lm ( θ, φ ) (cid:35) . (41)Note that for ρ = 1 we obtain f ( ω, θ, φ, t ) = G ( (cid:126)ω ) δ ( θ − Θ) δ ( φ − Φ)sin θ indicating full synchrony. For ρ = 0, on the other hand, f ( ω, θ, φ, t ) = G ( (cid:126)ω )4 π and the oscillators are uniformly spread overthe sphere.To sum the series and calculate the distribution f explicitly we use the relations4 π l + 1 l (cid:88) m = − l Y lm (ˆ y ) Y ∗ lm (ˆ x ) = P l (ˆ x · ˆ y ) (42) ∞ (cid:88) l =0 y l P l ( x ) = 1 (cid:112) y − xy (43)Applying these relations to (41) we find f ( (cid:126)ω, θ, φ, t ) = G ( (cid:126)ω )(1 − ρ )4 π (1 + ρ − ρ ˆ r · ˆ ρ ) . (44)Notice that f depends on Θ and Φ only throughˆ r · ˆ ρ ≡ cos ξ = sin θ sin Θ cos( φ − Φ) + cos θ cos Θ . (45)In addition, using the relation 4 π (cid:88) m = − Y m (ˆ y ) Y ∗ m (ˆ x ) = ˆ x · ˆ y (46)and the orthogonality of the spherical harmonics we can write the order parameter as (cid:126)p = (cid:90) G ( (cid:126)ω ) ρ π (cid:20)(cid:90) ˆ r (cid:48) ˆ r (cid:48) T d Ω (cid:48) (cid:21) · ˆ ρd ω (47)2The angular integral of the dyadic matrix can be easily performed and results (4 π/ . Weobtain (cid:126)p = (cid:90) G ( (cid:126)ω ) (cid:126)ρ d ω. (48)This is formally identical to the 2D case (see Eq. (18)), and much simpler than the relationfound in [1]. Note that this simplification was only possible because we applied orthogonalityproperties of the spherical harmonics. Also, Eq.(44) should be contrasted with the expressionderived in [1], where the denominator is raised to the power 2, instead of 3 /
2, and thenumerator (1 − ρ ) to power 2, instead of 1. These differences will have consequences forthe derivation of the equations of motion. C. Continuity equation for ansatz distribution
The partial derivatives of f are ∂f∂ρ = cos ξ (3 + ρ ) − ρ (5 − ρ ) D (49)and ∂f∂x = ∂f∂ cos ξ ∂ cos ξ∂x = 3 ρ (1 − ρ ) D ∂ cos ξ∂x , (50)for x = θ, φ, Θ , Φ where D ≡ ρ − ρ cos ξ and ∂ cos ξ∂θ = ρ θ ρ ∂ cos ξ∂φ = sin θ ρ φ ρ∂ cos ξ∂ Θ = ˆ r · ˆΘ ∂ cos ξ∂ Φ = sin Θˆ r · ˆΦ (51)Writing the time derivative of f as ∂f∂t = ∂f∂ρ ˙ ρ + ∂f∂ Θ ˙Θ + ∂f∂
Φ ˙Φ (52)we obtain D ∂f∂t = (cid:2) cos ξ (3 + ρ ) − ρ (5 − ρ ) (cid:3) ˙ ρ + 3 ρ (1 − ρ )ˆ r · [ ˙Θ ˆΘ + sin Θ ˙Φ ˆΦ] . (53)Using eqs.(38) and (49)-(51) we can also write (see Appendix B) D ∇ ( f(cid:126)v ) = 3 ρ (1 − ρ )ˆ r · ( ˆ ρ × (cid:126)ω ) + (1 − ρ ) (cid:2) − q r (1 + ρ ) + 4 q r ρ r + 3 q θ ρ θ + 3 q φ ρ φ (cid:3) . (54)3 D. Compatibility conditions
At this point we need to specify the vector (cid:126)q that will describe the coupling between theoscillators. The natural choice (cid:126)q = K (cid:126)p , does not work. To see this consider, for instance K = K I or (cid:126)q = K(cid:126)p . In this case the last three terms in Eq.(54) become4 ρ r p r + 3 ρ θ p θ + 3 ρ φ p φ = 3 (cid:126)ρ · (cid:126)p + ρ r p r = 3 (cid:126)ρ · (cid:126)p + ρp cos ξ cos ν (55)and the term cos ξ cos ν has no counterpart in the continuity equation. Chandra et al [1]have a similar problem, that they solve by choosing the exponent of their tentative ansatzfunction. Here we do not have this freedom. However, there is still a way out, using theinvariance of the exact equations under changes in the radial part of (cid:126)q . We can cancel theunwanted terms without altering the exact equations of motion if we choose (cid:126)q = (cid:20) −
14 ˆ r ˆ r (cid:62) (cid:21) K (cid:126)p + β ˆ r (56)where ˆ r ˆ r (cid:62) is the projection in the radial direction and the coefficient β will be chosen laterto ensure consistency of the equations of motion. With this choice Eq.(54) becomes D ∇ ( f(cid:126)v ) = 3 ρ (1 − ρ )ˆ r · ( ˆ ρ × (cid:126)ω )+ (1 − ρ ) (cid:20) − β (1 + ρ ) + 3 (cid:126)ρ · ( K (cid:126)p ) −
32 (1 + ρ )ˆ r · ( K (cid:126)p ) + 4 β ˆ r · (cid:126)ρ (cid:21) . (57) E. Equations of motion
The continuity equation for the ansatz function, after multiplying by D / , is given by (cid:2) cos ξ (3 + ρ ) − ρ (5 − ρ ) (cid:3) ˙ ρ + 3 ρ (1 − ρ )ˆ r · [ ˙Θ ˆΘ + sin Θ ˙Φ ˆΦ] + 3 ρ (1 − ρ )ˆ r · ( ˆ ρ × (cid:126)ω ) − (1 − ρ ) (cid:20) − β (1 + ρ ) + 3 (cid:126)ρ ( K (cid:126)p ) −
32 (1 + ρ )ˆ r · ( K (cid:126)p ) + 4 β ˆ r · (cid:126)ρ (cid:21) = 0 (58)We simplify this equation in Appendix C. This leads us to choose β = (cid:126)ρ · ( K (cid:126)p )4 . (59)and the final equations of motion become˙ (cid:126)ρ = (cid:126)ω × (cid:126)ρ + 12 (1 + ρ )( K (cid:126)p ) − [ (cid:126)ρ · ( K (cid:126)p )] (cid:126)ρ. (60)4From this equation it also follows that˙ ρ = 12 (1 − ρ )[ ˆ ρ · ( K (cid:126)p )] . (61)The connection between (cid:126)p and (cid:126)ρ is given by Eq.(48). Eq.(60) is identical to the equationobtained in [1], but Eq.(48) is not. In our formulation the order parameter (cid:126)p ( t ) is the averageover the distribution of natural frequencies of the ansatz parameter (cid:126)ρ ( ω, t ), just like the 2-dimensional case, whereas the expression in [1] includes an extra integral that depends onthe dimension of system. For D = 2 both formulations recover the original result in [22]. F. Identical oscillators with symmetric coupling and external forces
If all oscillators have identical natural frequencies, G ( (cid:126)ω ) = δ ( (cid:126)ω − (cid:126)ω ), then (cid:126)p ( t ) = (cid:126)ρ ( (cid:126)ω , t ).Moreover, if the coupling matrix is diagonal, K = K , the order parameter satisfies˙ (cid:126)p = (cid:126)ω × (cid:126)p + K − p ) (cid:126)p and ˙ p = K − p ) p. (62)The equilibrium points are p = 0, which is stable for K < p = 1, stable for K > p = 1.Time dependent external forces can be included making K · (cid:126)p → K · (cid:126)p + (cid:126)F in Eqs. (56)and (59)-(61). For the case of diagonal coupling matrix this gives˙ (cid:126)p = (cid:126)ω × (cid:126)p + (cid:20) K − p ) − (cid:126)p · (cid:126)F (cid:21) (cid:126)p + 12 (1 + p ) (cid:126)F (63)and ˙ p = 12 (1 − p )( Kp + (cid:126)F · ˆ p ) . (64)5 p t ( a ) - 1 . 0 - 0 . 5 0 . 0 0 . 5 1 . 00 . 00 . 51 . 0 - 1 . 0 - 0 . 5 0 . 0 0 . 5 1 . 0- 1 . 0 - 0 . 5 0 . 0 0 . 5 1 . 00 . 00 . 51 . 0 - 1 . 0 - 0 . 5 0 . 0 0 . 5 1 . 0- 1 . 0 - 0 . 5 0 . 0 0 . 5 1 . 00 . 00 . 51 . 0 - 1 . 0 - 0 . 5 0 . 0 0 . 5 1 . 0 p z p y p x Figure 3. Time evolution of the order parameter for identical oscillators with (cid:126)ω = ˆ z and K = 1.(a) module p ( t ); (b) evolution on the p x − p y − p z space. Thick black line shows the result ofsimulations with N = 10000 oscillators. Thin red and blue line shows the ansatz solution usingour approach and that from [1], respectively. G. Equilibrium solutions for symmetric coupling
Here we consider equilibrium solutions for K = K but general distributions of naturalfrequencies. In this case Eqs. (60) and (61) simplify to:˙ (cid:126)ρ = (cid:126)ω × (cid:126)ρ + K ρ ) (cid:126)p − (cid:126)ρ · (cid:126)p ) (cid:126)ρ ] and ˙ ρ = K − ρ )( ˆ ρ · (cid:126)p ) . (65)We first consider the limit K →
0, or K = 0 + . According to Eqs. (65), for (cid:126)p (cid:54) = 0 equilibriumrequires ρ = 1 and ˆ ρ parallel to ˆ ω , i.e., (cid:126)ρ = ± ˆ ω . However, only one of these solutions isstable. To find out which we make (cid:126)ρ = ± ˆ ω + δ(cid:126)ρ , or ρ = 1 + δρ and expand the moduleEq.(65) to first order to obtain δ ˙ ρ = − K δρ ( (cid:126)p · ˆ ρ ) . (66)Therefore, for each (cid:126)ω the stable solution is the one for which (cid:126)p · ˆ ρ >
0, i.e., the one in thehemisphere defined by the direction of (cid:126)p . Setting (cid:126)p = p ˆ z we have for G ( (cid:126)ω ) = g ( ω ) / πp = (cid:90) G ( (cid:126)ω ) ρ z d ω = 14 π (cid:90) π d Φ (cid:90) π/ d Θ = 12 , (67)where the factor 2 comes from the fact that, in the upper hemisphere, the stable solution isˆ ρ = ˆ ω and in the lower hemisphere it is ˆ ρ = − ˆ ω . So when we cross from one hemisphere toanother both cos( θ ) and ˆ ρ changes signal, making the integral over 0 < θ < π symmetricalwith respect to π/ K (cid:54) = 0 we define the orthonormal set of vectorsˆ ξ = (cos θ cos φ, cos θ sin φ, − sin θ )ˆ ξ = ( − sin φ, cos φ,
0) (68)ˆ ω = (sin θ cos φ, sin θ sin φ, cos θ )and expand (cid:126)ρ = α ˆ ξ + β ˆ ξ + γ ˆ ω. For K = 0 we use α = β = 0 and γ = 1, with 0 < θ < π/
2. We again set (cid:126)p = p ˆ z andobtain (cid:126)p = − p sin θ ˆ ξ + p cos θ ˆ ω. (69)Projecting Eq.(60) with ˙ (cid:126)ρ = 0 in these directions we obtain0 = − β − ¯ K sin θ (1 + ρ ) / − ¯ Kα ( − α sin θ + γ cos θ )0 = α − ¯ Kβ ( − α sin θ + γ cos θ ) (70)0 = ¯ K cos θ (1 + ρ ) / − ¯ Kγ ( − α sin θ + γ cos θ ) . where ρ = α + β + γ and ¯ K = Kp/ω . These equations can be solved analytically: β = − (1 + ¯ K )2 ¯ K sin θ (cid:34) − (cid:115) − K sin θ (1 + ¯ K ) (cid:35) (71) α = − (cid:115) β (cid:18) K cos θ ¯ K sin θ (cid:19) (72) γ = ¯ Kβα cos θ (73)and satisfy ρ = α + β + γ = 1.For a delta function distribution G ( (cid:126)ω ) = δ ( ω − ω )4 πω we find p = (cid:90) π/ ρ z ( ¯ K, θ ) sin θdθ ≡ F ( ¯ K ) (74)where ρ z = − α sin θ + γ cos θ depends only on ¯ K = Kp/ω and θ . Once the function F ( ¯ K ) is calculated numerically, p ( K ) can be computed in parametric form for ω = 1 as( p = F ( x ) , K = x/F ( x )). The result is shown in Fig. 4. Equilibrium is found only for K > .
5. For
K < . K = 2 our solution is slightly closer to theexact, whereas Chandra’s is advanced, reaching equilibrium before the exact. In both cases7 p K ( a ) ( b ) p t ( c ) p t Figure 4. (a) Order parameter as a function of the coupling constant for the uniform angular distri-bution of natural frequencies and | (cid:126)ω i | = 1. The continuous curve shows the numerical integrationfrom Eq. (74) and the dots the value obtained by integrating the exact equations of motion (31).For K < . p ( t ) for K = 1; (c) p ( t ) for K = 2. Thin red and blue line shows the ansatz solutionusing our approach and that from [1], respectively. the solutions present fluctuations that depend on the (random) initial conditions: every timethe equations are integrated a slightly different curve is obtained.For small K we can obtain analytical expressions expanding Eqs. (73) to second order.We find α ≈ − ¯ K sin θ cos θ , β ≈ − ¯ K sin θ and γ ≈ − ( ¯ K /
2) sin θ . This gives p = 12 + K p ω ≈
12 + K ω . (75)In the limit of very large coupling ¯ K → ∞ , on the other hand, we get α ≈ − sin θ , β ≈− sin θ/ ¯ K and γ ≈ cos θ , which gives p = (cid:90) π/ ( − α sin θ + γ cos θ ) sin θdθ = (cid:90) π/ sin θdθ = 1 . (76)For the Gaussian distribution G ( (cid:126)ω ) = √ √ π ∆ π e − ω / we have p = √ √ π ∆ (cid:90) π/ dθ sin θ (cid:90) ∞ dω ω e − ω / ρ z ( Kp/ω, θ ) . (77)Changing the integration variable to Ω = ω/Kp and setting ∆ = 1 we obtain p = √ √ π K p (cid:90) π/ dθ sin θ (cid:90) ∞ d Ω Ω e − Ω p K / ρ z (1 / Ω , θ ) ≡ H ( Kp ) . (78)Once again the function H ( x ) can be computed numerically and the parametric curve isgiven by ( p = H ( x ) , K = x/H ( x )). The result is shown in Fig. 5. For small K we obtain p = 12 + K p ≈
12 + K . (79)8 p K ( a ) ( b ) p t Figure 5. (a) Order parameter as a function of the coupling constant for the Gaussian distributionof natural frequencies with ∆ = 1. The continuous curve shows the numerical integration fromEq. (78) and the dots the value obtained by integrating the exact equations of motion (31). Thedashed line is the quadratic approximation, Eq.(79). (b) p ( t ) for K = 2. Thin red and blue lineshows the ansatz solution using our approach and that from [1], respectively. IV. CONCLUSIONS
In this paper we studied the problem of dimensional reduction for the Kuramoto modelin 2 and 3 dimensions. In 2D we used the ansatz proposed by Ott and Antonsen [22] andsolved the equations for the ansatz parameters at equilibrium for arbitrary distributions ofnatural frequencies, re-deriving Kuramoto expressions near the bifurcation point.We extended the ansatz to 3D expanding the distribution function in spherical harmonics.The resulting function differs from that obtained by Chandra et al [1] and satisfies thecontinuity equation only if the coupling term is modified. Our ansatz is connected to aespecific choice of the coupling between the oscillators. This modification, however, is only inits radial component, which does not affect the exact equations of motion. The interpretationof the approximate density function is very similar to that in 2D: for ρ → ρ → , Φ). The continuity equation leads to a vector equation for the ansatz parameters thatis exactly that obtained in [1]. However, the connection between the ansatz and the orderparameter differs and it is simpler in our approach. The relationship we find is the naturalextension of the 2D case, where (cid:126)p is the integral of (cid:126)ρ over G ( (cid:126)ω ), whereas that found in [1]9is more complicated. This might facilitate the analysis of more general systems where thecoupling is a full matrix K or in the presence of external forces.As applications we consider two types of G ( (cid:126)ω ) consisting of uniform angular distributionswith delta or Gaussian distribution of frequency module. In both cases the ansatz equationsdescribe very well the full system and show the first order transition behavior as expected.Curiously, the case of identical oscillators, although agrees exactly at equilibrium, shows adelay in the dynamics, as shown in figure 2. Among the examples we explored, this was theonly case where the ansatz in [1] was more accurate than ours.The complete dimensional reduction is only achieved for the case of identical oscillators.For all other cases we have to solve Eq.(60) together with (48). However, as pointed outin [1], the latter can be solved by Monte Carlo methods, sampling the distribution G ( (cid:126)ω ),which converges fast for most distributions of interest.Finally, we remark that our ansatz points to an alternative way to treat the Kuramotomodel in higher dimensions. It allows us to easily obtain the dynamics of the order param-eter and study the behavior at equilibrium for different distributions of oscillator’s naturalfrequencies. ACKNOWLEDGMENTS
This work was partly supported by FAPESP, grants 2019/2027-5 (MAMA), 2019/24068-0(ANDB), 2016/01343-7 (ICTP-SAIFR FAPESP) and CNPq, grant 301082/2019-7 (MAMA).We would like to thank Alberto Saa and Jose A. Brum for suggestions and careful readingof this manuscript.
Appendix A: Derivation of the dynamics of ρ in 2D Rewriting Eq.(21) explicitly as ∂f∂ρ ˙ ρ + ∂f∂ Φ ˙Φ + ( ω + q φ ) ∂f∂φ − f q r = 0 (A1)and replacing the derivatives we find, after multiplying by D , − ρ − (1 + ρ ) cos ξ ] ˙ ρ − ρ φ (1 − ρ ) ˙Φ +2 ω (1 − ρ ) ρ φ + 2(1 − ρ )( ρ φ q φ + ρ r q r ) − (1 − ρ ) q r = 0 (A2)0Using (cid:126)q = K (cid:126)p and the relations ρ φ = (cid:126)ρ · ˆ φ = − ρ ˆ r · ˆΦ, ρ r = (cid:126)ρ · ˆ r and ωρ φ = − ˆ r · ( (cid:126)ω × (cid:126)ρ )we can greatly simplify Eq.(A2) to2(1 − ρ )ˆ r · (cid:20) (cid:126)ρ × (cid:126)ω + ρ ˙Φ ˆΦ + ˙ ρ ˆ ρ (1 + ρ )1 − ρ −
12 (1 + ρ ) K (cid:126)p (cid:21) − ρ ˙ ρ + 2(1 − ρ ) (cid:126)ρ · K (cid:126)p = 0 . Adding and subtracting ˙ ρ ˆ ρ inside the square brackets we obtain2(1 − ρ )ˆ r · (cid:20) (cid:126)ρ × (cid:126)ω + ˙ (cid:126)ρ + ˙ ρ ˆ ρ (2 ρ )1 − ρ −
12 (1 + ρ ) K (cid:126)p (cid:21) + 4 (cid:20)
12 (1 − ρ )( K (cid:126)p · (cid:126)ρ ) − ρ ˙ ρ (cid:21) = 0 . Since this equation must hold for all values of ˆ r each term in the brackets must be zero:˙ (cid:126)ρ = (cid:126)ω × (cid:126)ρ − ˙ ρ ˆ ρ (2 ρ )1 − ρ + 12 (1 + ρ ) K (cid:126)p (A3)˙ ρ = 12 (1 − ρ )( K (cid:126)p · ˆ ρ ) . (A4)Taking the scalar product of the first of these equations with (cid:126)ρ and noting that (cid:126)ρ · ˙ (cid:126)ρ = ρ ˙ ρ we can check that the second equation is recovered. Therefore the equations are compatibleand we can replace ˙ ρ given by the second into the first equation to finally obtain˙ (cid:126)ρ = (cid:126)ω × (cid:126)ρ + 12 (1 + ρ ) K (cid:126)p − ( K (cid:126)p · (cid:126)ρ ) (cid:126)ρ. (A5) Appendix B: Derivation of Eq.(54)
From Eq.(52) we obtain D ∂f∂t = (cid:2) cos ξ (3 + ρ ) − ρ (5 − ρ ) (cid:3) ˙ ρ + 3 ρ (1 − ρ )ˆ r · [ ˙Θ ˆΘ + sin Θ ˙Φ ˆΦ] (B1)Using eq.(38) and eqs. (49)-(51) we write D ∇ ( f(cid:126)v ) = − − ρ ) q r D + ( ω φ + q θ )3 ρ (1 − ρ ) ∂ cos ξ∂θ + 1sin θ ( − ω θ + q φ )3 ρ (1 − ρ ) ∂ cos ξ∂φ = 3 ρ (1 − ρ ) (cid:20) ω φ ∂ cos ξ∂θ − ω θ sin θ ∂ cos ξ∂φ (cid:21) + (cid:20) − − ρ ) q r D + 3 ρ (1 − ρ ) q θ ∂ cos ξ∂θ + 3 ρ (1 − ρ ) q φ sin θ ∂ cos ξ∂φ (cid:21) . (B2)The first line contains terms that are independent of the coupling vector (cid:126)q . They can besimplified as follows: ω φ ∂ cos ξ∂θ − ω θ sin θ ∂ cos ξ∂φ = 1 ρ ( ω φ ρ θ − ω θ ρ φ ) = ˆ r · ( ˆ ρ × (cid:126)ω ) . (B3)1Finally, using eqs. (49)-(51), the terms containing (cid:126)q can be written as − − ρ ) q r D + 3(1 − ρ )( q θ ρ φ + q φ ρ θ )= (1 − ρ ) (cid:2) − q r (1 + ρ ) + 4 q r ρ r + 3 q θ ρ θ + 3 q φ ρ φ (cid:3) . (B4) Appendix C: Derivation of the dynamics of ρ in 3D To simplify Eq.(58) we first add ˙ ρ ˆ ρ to the middle term of the first line to obtain 3 ρ (1 − ρ )ˆ r · ˙ (cid:126)ρ . We also subtract an equal term to keep the equation unaltered. This gives3(1 − ρ )ˆ r · [ (cid:126)ω × (cid:126)ρ − ˙ (cid:126)ρ ] + ˙ ρ (cid:2) ρ (5 − ρ ) − ρ cos ξ (cid:3) − (1 − ρ ) (cid:20) − β (1 + ρ ) + 3 (cid:126)ρ ( K · (cid:126)p ) −
32 (1 + ρ )ˆ r · ( K · (cid:126)p ) + 4 β ˆ r · (cid:126)ρ (cid:21) = 0 (C1)Next we put all terms containing the angular coordinates ˆ r together. We obtain3(1 − ρ )ˆ r · [ (cid:126)ω × (cid:126)ρ − ˙ (cid:126)ρ + 12 (1 + ρ )( K · (cid:126)p ) − β (cid:126)ρ − − ρ ) ρ ˙ ρ(cid:126)ρ ]+ ˙ ρρ (5 − ρ ) − (1 − ρ ) (cid:2) (cid:126)ρ · ( K · (cid:126)p ) − β (1 + ρ ) (cid:3) = 0 (C2)This implies that the following equations must be satisfied simultaneously:˙ (cid:126)ρ = (cid:126)ω × (cid:126)ρ + 12 (1 + ρ )( K · (cid:126)p ) − β (cid:126)ρ − − ρ ) ρ ˙ ρ(cid:126)ρ (C3)and ˙ ρρ (5 − ρ ) = (1 − ρ ) (cid:2) (cid:126)ρ · ( K · (cid:126)p ) − β (1 + ρ ) (cid:3) (C4)Taking the scalar product of (C3) with (cid:126)ρ and noting that (cid:126)ρ · ˙ (cid:126)ρ = ρ ˙ ρ we find that equations(C3) and (C4) are compatible only if β = (cid:126)ρ · ( K · (cid:126)p )4 . (C5)With this condition (C3) and (C4) become˙ (cid:126)ρ = (cid:126)ω × (cid:126)ρ + 12 (1 + ρ )( K · (cid:126)p ) − [ (cid:126)ρ · ( K · (cid:126)p )] (cid:126)ρ (C6)and ˙ ρ = 12 (1 − ρ )[ ˆ ρ · ( K · (cid:126)p )] . (C7)2 [1] Sarthak Chandra, Michelle Girvan, and Edward Ott. Complexity reduction ansatz for systemsof interacting orientable agents: Beyond the kuramoto model. Chaos: An InterdisciplinaryJournal of Nonlinear Science , 29(5):053107, 2019.[2] Henrique M Oliveira and Lu´ıs V Melo. Huygens synchronization of two pendulum clocks. arXiv preprint arXiv:1410.7926 , 2014.[3] James Pantaleone. Synchronization of metronomes.
American Journal of Physics , 70(10):992–1000, 2002. doi:10.1119/1.1501118. URL https://doi.org/10.1119/1.1501118 .[4] Istv´an Z. Kiss, Yumei Zhai, and John L. Hudson. Resonance clustering in globally cou-pled electrochemical oscillators with external forcing.
Phys. Rev. E , 77:046204, Apr 2008.doi:10.1103/PhysRevE.77.046204. URL https://link.aps.org/doi/10.1103/PhysRevE.77.046204 .[5] a. V. Novikov and E. N. Benderskaya. Oscillatory neural networks based on the Kuramotomodel for cluster analysis.
Pattern Recognition and Image Analysis , 24(3):365–371, 2014.ISSN 1054-6618. doi:10.1134/S1054661814030146. URL http://link.springer.com/10.1134/S1054661814030146 .[6] Yoshiki Kuramoto. Self-entrainment of a population of coupled non-linear oscillators. In
International Symposium on Mathematical Problems in Theoretical Physics , pages 420–422.Springer-Verlag, Berlin/Heidelberg, 1975. doi:10.1007/BFb0013365. URL .[7] Yoshiki Kuramoto. Chemical Waves. In
Chemical Oscillations, Waves, and Turbulence , pages89–110. Springer Berlin Heidelberg, 1984. doi:10.1007/978-3-642-69689-3˙6. URL .[8] Hyunsuk Hong and Steven H Strogatz. Kuramoto model of coupled oscillators with posi-tive and negative coupling parameters: an example of conformist and contrarian oscillators.
Physical Review Letters , 106(5):054102, 2011.[9] MK Stephen Yeung and Steven H Strogatz. Time delay in the kuramoto model of coupledoscillators.
Physical Review Letters , 82(3):648, 1999.[10] Michael Breakspear, Stewart Heitmann, and Andreas Daffertshofer. Generative models ofcortical oscillations: neurobiological implications of the kuramoto model.
Frontiers in human neuroscience , 4:190, 2010.[11] Francisco A. Rodrigues, Thomas K D M Peron, Peng Ji, and J??rgen Kurths. The Ku-ramoto model in complex networks. Physics Reports , 610:1–98, 2016. ISSN 03701573.doi:10.1016/j.physrep.2015.10.008. URL http://dx.doi.org/10.1016/j.physrep.2015.10.008 .[12] Joyce S. Climaco and Alberto Saa. Optimal global synchronization of partially forced ku-ramoto oscillators.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 29(7):073115,2019. doi:10.1063/1.5097847. URL https://doi.org/10.1063/1.5097847 .[13] Jesus Gomez-Gardenes, Sergio Gomez, Alex Arenas, and Yamir Moreno. Explosive synchro-nization transitions in scale-free networks.
Physical Review Letters , 106(12):1–4, 2011. ISSN00319007. doi:10.1103/PhysRevLett.106.128701.[14] Peng Ji, Thomas K Dm Peron, Peter J. Menck, Francisco A. Rodrigues, and J??rgen Kurths.Cluster explosive synchronization in complex networks.
Physical Review Letters , 110(21):1–5,2013. ISSN 00319007. doi:10.1103/PhysRevLett.110.218701.[15] Juan A. Acebr´on, L. L. Bonilla, Conrad J P´erez Vicente, F´elix Ritort, and Renato Spigler.The Kuramoto model: A simple paradigm for synchronization phenomena.
Reviews of ModernPhysics , 77(1):137–185, 2005. ISSN 00346861. doi:10.1103/RevModPhys.77.137.[16] Florian D¨orfler and Francesco Bullo. On the critical coupling for kuramoto oscillators.
SIAMJournal on Applied Dynamical Systems , 10(3):1070–1099, 2011.[17] Simona Olmi, Adrian Navas, Stefano Boccaletti, and Alessandro Torcini. Hysteretic transi-tions in the kuramoto model with inertia.
Physical Review E , 90(4):042905, 2014.[18] Lauren M. Childs and Steven H. Strogatz. Stability diagram for the forced Kuramoto model.
Chaos , 18(4):1–9, 2008. ISSN 10541500. doi:10.1063/1.3049136.[19] Carolina A Moreira and Marcus AM de Aguiar. Global synchronization of partially forcedkuramoto oscillators on networks.
Physica A: Statistical Mechanics and its Applications , 514:487–496, 2019.[20] Carolina A Moreira and Marcus AM de Aguiar. Modular structure in c. elegans neuralnetwork and its response to external localized stimuli.
Physica A: Statistical Mechanics andits Applications , 533:122051, 2019.[21] Sarthak Chandra, Michelle Girvan, and Edward Ott. Continuous versus discontinuous tran-sitions in the d-dimensional generalized kuramoto model: Odd d is different.
Physical Review X , 9(1):011002, 2019.[22] Edward Ott and Thomas M. Antonsen. Low dimensional behavior of large systems of globallycoupled oscillators. Chaos , 18(3):1–6, 2008. ISSN 10541500. doi:10.1063/1.2930766.[23] Shinya Watanabe and Steven H Strogatz. Constants of motion for superconducting josephsonarrays.
Physica D: Nonlinear Phenomena , 74(3-4):197–253, 1994.[24] Charles J. Goebel. Comment on “constants of motion for superconductor arrays”.
Physica D: Nonlinear Phenomena , 80(1):18 – 20, 1995. ISSN 0167-2789. doi:https://doi.org/10.1016/0167-2789(95)90049-7. URL .[25] Seth A. Marvel, Renato E. Mirollo, and Steven H. Strogatz. Identical phase oscillators withglobal sinusoidal coupling evolve by m¨obius group action.
Chaos: An Interdisciplinary Journalof Nonlinear Science , 19(4):043104, 2009. doi:10.1063/1.3247089. URL https://doi.org/10.1063/1.3247089 .[26] Takuma Tanaka. Solvable model of the collective motion of heterogeneous particles interactingon a sphere.
New Journal of Physics , 16, 01 2014. doi:10.1088/1367-2630/16/2/023016.[27] Xin Hu, S Boccaletti, Wenwen Huang, Xiyun Zhang, Zonghua Liu, Shuguang Guan, and Choy-Heng Lai. Exact solution for first-order synchronization transition in a generalized kuramotomodel.