Kuramoto Phase Model with Inertia : Bifurcations Leading to the Loss of Synchrony and to the Emergence of Chaos

We consider a finite-dimensional model of phase oscillators with inertia in the case of star configuration of coupling. The system of equations is reduced to a nonlinearly coupled system of pendulum equations. We prove that the transition from synchronous to asynchronous oscillations occurs via bifurcation of saddle-node equilibrium. In this connection the asynchronous regime can be partially synchronous rotations. We find that the reverse transition from asynchronous to synchronous regime occurs via bifurcation of homoclinic orbit both of the saddle equilibrium point and of the saddle periodic orbit. In the case of homoclinic loop of the saddle point the synchrony appears only from asynchronous mode without partially synchronized rotations. In the case of the homoclinic curve of the saddle periodic orbit the system undergoes a chaotic rotation regime which results in a random return to synchrony. We establish that return transitions are hysteretic in the case of large inertia.

Synchronous behavior is one of ubiquitous collective phenomena in ensembles of oscillatory systems.The Kuramoto model is arguably the most studied model that describes synchronization [1,2,3].This model captures essential features of synchronization, observed in science and applications.Examples are arrays of coupled Josephson junctions [4], semiconductor laser arrays [5], the ensembles of the cells in the heart [6], Hodgkin-Huxley neurons [7], central pattern generator for animal locomotion [8], rhythmic applause [9], pedestrian crowd synchrony on London's Millennium bridge [10], microwave oscillator arrays [11] etc.For other examples, see [12], [13] and [14].In this paper, an analytical study of a star motif of Kuramoto oscillators is presented.The system consists of a central node, the hub, connected to an arbitrary number of peripheral nodes.The star configuration of ensemble of oscillators for the different models has been studied.In [15], starting from the analysis of the topological properties of the star configuration, some analytical considerations have been applied to derive the bifurcation diagram of the system with respect to the parameter mismatch between peripheral oscillators and hub and to the coupling strength.The analysis revealed that the system may become fully synchronized (more precisely, the peripheral oscillators are completely synchronized among each other and phase synchronized with the hub).In the case of star-coupled ensemble of phase oscillators the analytical description of the parameter regions of existence of different synchronous regimes has been obtained in [16].It was shown that peripheral oscillators compete for the synchronization with the hub and only a given number of peripheral oscillators can win this competition.

Kuramoto phase model with inertia
We consider the Kuramoto phase model with inertia of N coupled phase oscillators [17,18] where θ i is instantaneous phase, ω i is natural frequency of the i-th oscillator, K i,j are the entries of a coupling symmetric matrix K = {K ij } N N , β is a positive parameter representing an inertia of oscillators.In the case β = 0 and K ij = K, i, j = 1, 2, . . ., N the system (1) becomes the original paradigmatic Kuramoto model [1,2,3].
We consider the star configuration of coupling [19,20,21] when the matrix K has the form i.e. the first element is the hub of the configuration, and the system (1) reads We introduce new variables and new parameters where φ i is the phase difference between the hub and each another peripheral oscillator, ∆ i is the frequency mismatch of (i + 1)'s peripheral oscillator and the first hub oscillator.Using (3) we rewrite the system (2) in the form The phase synchronization of oscillators in the model ( 1) is defined as an attractor of the system (4) which trajectories (φ s i (t), φs i (t)) satisfy the conditions where • denotes the mean value, and parameter ε < π/2 is a measure of synchronization.
Respectively, the steady trajectories (φ a i (t), φa i (t)) of the system (4) defining asynchronous mode of the oscillators in (1) satisfy the condition The mutual synchronization of the peripheral oscillators is characterized by the rotation numbers , i, j = 1, 2, . . ., n.
The main candidate for synchronous state of the system (1) is the stable equilibrium point of the system (4).
From the boundness of the right parts in (4) follows the next where a = max i∈{1,2,...,n} a i , is the absorbing domain of the system (4).
Equilibria of the system (4) are the solutions of the system where Since det M = n + 1 = N the system (8) has a unique solution where ∆i are the entries of the column M −1 Γ, and (9) reads The system (9) for | ∆i | < a i has 2 n equilibria in G.The principal equilibrium corresponding to the synchronous mode is the point The stability of the equilibria is defined by the variational linear system of ODE where We seek a solution of the system (11) in the form x i = c i e pt and obtain the next characteristic equation for the system ( 11) where σ i = (βp 2 + p + 2α i )α −1 i .Hence, the equilibrium O is stable, i.e. the complete synchronization regime is stable, when the real parts of 2n roots of the equation (12) for Let us consider the particular case of the system ( 4) This condition implies the uniform coupling in the system (4): a i = a, ∆ i = ∆.In this case ∆i = ∆/N , and α i = a 2 i − (∆/N ) 2 = α, and σ i = σ.In this case the equation ( 12) takes the form Then σ − 1 = 0 gives the equation and σ + n − 1 = 0 leads to the equation From ( 15) repeated (n − 1) times and ( 16) we conclude that the real parts of all the roots of the equation ( 12) are negative and therefore the equilibrium O is asymptotically stable.In this case from equation ( 2) we obtain the following expression for frequency of complete synchronization It's easy to verify that in homogeneous case (13) all the rest equilibria are saddles with different dimensions of unstable manifolds.Hence, we proved the next statement Then the system (4) has the stable equilibrium point O, corresponding to the synchronous mode of the system (1) when the hub oscillator synchronizes the enclosing ones.
Corollary 1.The stability of the equilibrium point O is preserved for small mismatch |α i − α| < µ due to its structural stability.For large mismatch the stability conditions one can derive using (12).
Indeed, in this case the system (9) has no solution, and the synchronization loss occurs due to disappearance of the stable equilibrium O via saddle-node bifurcation.

The uniform coupling in star configuration of Kuramoto model
Consider the case of uniform coupling when a i = a = const, ∆ i = ∆ = const, but rewrite the system (4) in another form where new parameter b > 0 is not necessary equal to a.
Lemma 1.The system (19) has the invariant manifold M : The dynamics in the manifold M is determined by the pendulum equation where α = a + bn.
Indeed, each equation ( 19) after substitution u i = u becomes one and the same equation (20), and any trajectory of the system (19) The local stability of the manifold M is defined by the variational equation where ξ i = φ i − φ and φ is driven by the system (20).First we present well known [22,23,24] bifurcational diagram and qualitative phase pictures of the pendulum equation ( 20) which for new time t = α β t and new parameters λ = (αβ) −1/2 , γ = ∆/α takes the form For λ > 0 the bifurcations in this equation are saddle-node for |γ| = 1 and the homoclinic loop encircling the cylinder (φ, y = φ) at |γ| = γ h (λ), where γ h (λ) is Tricomi curve satisfying the conditions where the value λ sn ≈ 1.2 corresponds to the homoclinic loop of the saddle-node.The condition γ h (0) = 4 π one can obtain using averaging method for small parameters λ and γ in (22), and increasing property, γ h (λ) > 0, of the function γ = γ h (λ) follows from clockwise rotation of the vector field ( φ = y, ẏ) given by the equation ( 22), while parameter λ increases.
For the parameters β, α, ∆ of the equation ( 20) the bifurcations read is the homoclinic bifurcation, and is the saddle-node bifurcation.
Lemma 2. The system given by the equation ( 20) in the phase cylinder G = {φ ∈ S 1 , y ∈ R 1 } has the following phase portraits: 1) In the parameter domain the system ( 27) is globally asymptotically stable (Fig. 1,S) such that the stable focus , y f = 0) attracts the whole cylinder besides the stable separatrices of the saddle O s (φ s = π − φ f , y s = 0).
2) In the parameter domain the system ( 27) is bistable: it has the stable focus (node) and the unique stable limit cycle l c (φ = φ c (t), y = y c (t)) encircling the cylinder; the basins of the focus and the cycle are separated by the stable separatrices of the saddle (Fig. 1,B).
3) In the parameter domain d 3 : {|∆| > α} the system (27) has the unique limit cycle attracting the whole cylinder (Fig. 1,R).Consider the local stability of the trajectories l * (φ * (t), y * (t)) in the invariant manifolds M , especially of the limiting set, which consist of the stable focus (node) O f (φ, 0), the saddle O s (φ s , 0) and limit cycle l c (φ c , y c ).
In the transverse direction to the manifold (transversally to the vector (1, 1, . . ., 1)) the system (21) takes the form where , providing cos φ * (t) > 0, at least for (βa) −1/2 ≥ λ sn guarantees the local stability of this part of the manifold M .The transverse stability of the limit cycle is defined by the equation (29) for φ * = φ c (t).Since φc (t) = y c (t) > 0 the phase φ c (t) rotates and the term cos φ c (t) in (29) changes the sign thereby creating a problem of the cycle transverse stability.We solve it in the case when |∆| = (βa) −1/2 + ε, for small enough ε > 0. In this case the cycle just appearing from the homoclinic loop of the system (27) passes a small neighborhood of the saddle O s and therefore spends the most time (of order 1/ε) in the neighborhood |φ − φ s | < ε.Since cos φ s < 0, due to (29) the limit cycle l c is born being unstable.
From above reasoning we conclude Statement 2. 1) If the Lyapunov-Floquet exponents from (29) for φ * (t) = φ c (t) are negative then the asynchronous mode is such that the peripheral oscillators are synchronized with rotation numbers equaled 1.
2) The homoclinic bifurcation of the system (19) leads to an asynchronous mode of the peripheral oscillators.
For sufficiently large inertia such that (βa) −1/2 < λ sn , the transition from coherence to incoherence of oscillators is hysteretic.When the frequency difference ∆ increases the transition from the stable equilibrium O f (coherence) to the rotation mode in the solid torus G 0 (incoherence) occurs via the saddle-node bifurcation |∆| = a.Obviously, this rotation mode can be the stable cycle l c = (φ c (t), φc (t)) in the manifold M .In this case one observes the transition from complete phase synchronization to the synchronous state of the peripheral oscillators being asynchronous to the hub oscillator with mean frequency difference φc (t) .When the frequency difference ∆ decreases from large values corresponding to the rotation mode at the bifurcation of homoclinic orbit of the saddle |∆| = aγ h ((βa) −1/2 ) the reverse transition to the complete synchronization due to Statement 2 occurs only from the asynchronous mode of the peripheral oscillators.Note that this hysteretic behaviour being similar to the transitions in the Josephson junction model [25] was discussed in the recent paper [26].

Nonsymmetric coupling
We consider the general case of nonsymmetric coupling but, as an example, for three oscillators in the star configuration.Similarly to symmetrical case ( 2), ( 4) we obtain the system where β 1,2 are different inertias of peripheral oscillators, a 1,2 , b 1,2 are the coupling matrix entries, and ∆ 1,2 are frequency differences.Our goal is to introduce several parameter domains exhibiting different simple and complicated dynamics of the system (30).

Equilibria
The system (30) has four equilibria in the region giving solutions φ similarly to (9).For β 1 = β 2 = β the stability of the equilibria in this case is defined by the equation where r = a 1 α 1 +a 2 α 2 , α 1,2 = cos φ ) is stable and three other equilibria are saddles.

Comparison systems
We rewrite the equations (30) in the form of the systems Introduce two comparison 2D systems (see [27] and ref. within) for each subsystem in (34) A The vector projection of the system (35) on the cylinder (φ, y) is rotated clockwise (counterclockwise) relatively to the vector of A + i (A − i , respectively) in the half-cylinder (φ, y > 0) and vise versa in the half-cylinder (φ, y < 0).Now we depict the separatrices of the saddles and cycles simultaneously for the systems A + i and A − i .The unstable (stable) separatrices form the strips between them-sparatrix channels g u i (g s i respectively).Introduce the intersection g o i = g u i ∩g s i , called the saddle cell [27], the annulus K s i bounded by the stable cycles of the system A + i and A − i , and the absorbing domain g + i bounded by unstable separatrices of the systems A + i and A − i , and segments φ = const (see Fig. 2), i = 1, 2. We select three pairs of parameter domains standing for ∆, a and β respectively.Both systems A + i and A − i have the same qualitative phase portraits in each of these domains d ki , k = 1, 2, 3, forming mutual arrangement of the saddle channels, annulus and absorbing domain depicted in Fig. 2. Using the above geometric structures we obtain the following findings.Proof.The system (34) has no entire trajectories besides the saddles in the domains It follows from Lyapunov-Chetaev function for monotone functions sin φ 1,2 at the intervals of φ 1,2 for g o(+) 1,2 in the system (34).The domain g + = g + 1 ∪ g + 2 is the attracting domain of the trajectories of the system (34) (besides the stable manifolds of the saddles) for the parameter region d 11 ∪ d 12 due to the directing property of the comparison systems (see Fig. 2-1).The stability of the locally stable point O f in the globally attracting domain g + can be derived with the Lyapunov function using the monotonicity of sin φ 1,2 in the square |φ form the same structures of Fig. 2-3.Then the solid torus K + = K s 1 × K s 2 attracts all the trajectories of the system (34).The nonwandering set of trajectories in K + is rotating and defines the asynchronous mode of the oscillators.
This statement immediately follows from the simple structure of Fig. 2-3 forcing all trajectories of the system (34) to enter K + .
Theorem 4. Let the parameters of the system (34) be in the domain d 2 = d 21 ∪ d 22 with the same structures of comparison systems (Fig. 2-2).Then the system (34) is fourfoldstable, that is, it has four separate components of limiting set in four absorbing domains g Proof.In the parameter domain d 2 a trajectory of the system (34) given by a solution , is such that due to the comparison principal the coordinates of the first (second, respectively) subsystem remain in the first (second, respectively) absorbing domain, ( φi , ỹi ) ∈ g + i (K s i ), i = 1, 2, for any t > 0. This implies that in R 2 × T 2 the trajectory of the system (34) with any initial point in the domain g , respectively), stay in these domains for any t > 0.

Bifurcational transitions
First we note that in all cases of system the bifurcations of equilibria are simple and occur via the saddle-node when the frequency differences increase.In order to exhibit the complicated bifurcations leading to emergence of chaos we consider the reduced system (30) for b 2 = 0 corresponding to the unidirectional coupling of one of the peripheral oscillators.The second ("master") equation in (30) has the pendulum dynamics and in the simple case |∆ 2 | > a 2 which we consider has the unique rotating limit cycle φ c (t) = φ c (t + T ) (see Lemma 2).The first ("slave") equation is the pendulum one as well but it is driven by periodic force −b 1 sin φ c (t).Using the results from [25,27] we obtain the next 2) There exists an interval |∆ 1 −a 1 γ h ((β 1 a 1 ) −1/2 )| < ε corresponding to a structurally stable homoclinic orbit to the cycle l s providing a chaotic compement of the system (30) limiting set containing infinite numbers of saddle cycles.
The proof of the theorem is based [25,27] on the fact that when the parameter ∆ 1 increases from the values from the domain d 11 (corresponding to Fig. 2-1) up to the values from the domain d 12 (corresponding to Fig. 2-2) the channels g s 1 and g u 1 as well as the manifolds W s 1 and W u 1 change their mutual arrangement causing the birth, existence and death of the homoclinic orbits H 1 = W s 1 ∩ W u 1 .Omitting the details we note the important property of the hysteretic bifurcational transition in the system (30).When the parameter ∆ 1 decreases the transition from asynchronous rotation to the synchrony occurs at random values of the parameter ∆ 1 from the interval corresponding to the quasi-attractor existence.This complexity of the system (30) dynamics is similar to that of the shunted Josephson junction [27].
Hence, the example of three oscillators exhibits the complexity of the system dynamics which obviously is typical for the general system (1).

Fig. 1 .
Fig. 1.Phase portraits of the system (20) for different value of parameters

Theorem 2 .
Let the parameters of the system (34) be in the domain d 1 = d 11 ∪ d 12 when the comparison systems A +(−) 1 and A +(−) 2form the same structure of Fig.2-1.Then the equilibrium point O f (φ + 1 , φ + 2 ) is globally asymptotically stable.Herewith three oscillators are globally synchronized.

Theorem 5 . 1 ) 1 manifolds lie in the channels W s 1 ⊂ g s 1 ×K s 2 , W u 1 ⊂ g u 1
In the parameter region b 2 = 0, |∆ 2 | > a 2 , |∆ 1 | < a 1 −b 1 the system (30) in the solid torus g o 1 × K s 2 has a unique saddle cycle l s which stable W s 1 and unstable W u ×K s 2 and have mutual arrangement corresponding for d 11 to Fig.2-1 and for d 12 to Fig.2-2.