Dynamics of the Kuramoto-Sakaguchi Oscillator Network with Asymmetric Order Parameter

Dynamics of the Kuramoto-Sakaguchi Oscillator Network with Asymmetric Order Parameter

Bolun Chen, Jan R. Engelbrecht and Renato Mirollo Departments of Physics and Mathematics, Boston College, Chestnut Hill, MA 02467
Abstract

Abstract: We study the dynamics of a generalized version of the famous Kuramoto-Sakaguchi coupled oscillator model. In the classic version of this system, all oscillators are governed by the same ODE, which depends on the order parameter of the oscillator configuration. The order parameter is the arithmetic mean of the configuration of complex oscillator phases, multiplied by some constant complex coupling factor. In the generalized model we consider, the order parameter is allowed to be any complex linear combination of the complex oscillator phases, so the oscillators are no longer necessarily weighted identically in the order parameter. This asymmetric version of the K-S model exhibits a much richer variety of steady-state dynamical behavior than the classic symmetric version; in addition to stable synchronized states, the system may possess multiple stable (N-1,1) states, in which all but one of the oscillators are in sync, as well as multiple families of neutrally stable asynchronous states or closed orbits, in which no two oscillators are in sync. We present an exhaustive description of the possible steady state dynamical behaviors; our classification depends on the complex coefficients that determine the order parameter. We use techniques from group theory and hyperbolic geometry to reduce the dynamic analysis to a 2D flow on the unit disc, which has geometric significance relative to the hyperbolic metric. The geometric-analytic techniques we develop can in turn be applied to study even more general versions of Kuramoto oscillator networks.

pacs:
05.45.Xt,74.81.Fa

The Kuramoto-Sakaguchi coupled oscillator model is a famous, well-studied dynamical system that models the dynamics of a network of identical coupled oscillators. In the classic formulation, all the oscillators are driven by identical coupling to the system’s order parameter, which is just the average of the complex oscillator phases, multiplied by some complex coupling constant. We study a generalized version of the Kuramoto-Sakaguchi system, in which the order parameter is now a complex-linear combination of the complex oscillator phases, but not necessarily with identical coefficients. We analyzed the dynamics of this asymmetric K-S system, and found that it supports a much richer variety of dynamical behaviors than the classic K-S model, which typically has steady state dynamics that are completely synchronized or completely asynchronous with all oscillators out of sync with each other. The asymmetric K-S model also can support multiple stable  states, in which all but one of the oscillators are in sync, as well as multiple families of neutrally stable asynchronous states or closed orbits. We introduce new group-theoretic and geometric techniques to study the asymmetric K-S model, which effectively reduce the -dimensional dynamics to a 2D flow on the unit disc, and use the natural hyperbolic geometry on the disc to study this flow. The techniques we develop lead to a complete classification of the dynamics of the asymmetric K-S model in terms of the order parameter coefficients, and can be used to study more general oscillator networks via a similar dimensional reduction. This also builds a connection between the somewhat distant fields of oscillator network dynamics and hyperbolic geometry / low-D complex dynamics.

I Introduction

Our subject is the study of networks of Kuramoto oscillators, which are dynamical systems governed by equations of the form

(1)

Here is an angular variable (i.e. an element of ) and the coefficients are smooth functions of . The state space for this system is the -fold torus . Kuramoto oscillator networks often arise as idealized models of physical dynamical systems, like Josephson junction series arrays [1; 2], and also as the result of averaging more complex dynamical systems [3]. Beginning with the original work of Kuramoto over forty years ago [4], Kuramoto networks have been a very fertile research subject in applied dynamics [5]. When the functions are symmetric in the variables , we say the system is a network of identical Kuramoto oscillators; any permutation of the components of a solution results in another solution to the system. This is the case for the famous Kuramoto-Sakaguchi (K-S) model, which has equations

(2)

In this paper we will investigate the dynamics of Kuramoto networks where the functions are not symmetric in the . As we shall see, dropping the symmetry assumption leads to a richer variety of dynamic behavior. The focus of our work is a variation of the K-S model: the asymmetric K-S network given by

(3)

The dynamics of this network are governed by its order parameter, which we can express in complex form, with and , as

(4)

It turns out, broadly speaking, that the dynamics depend largely on the sum of the , which we denote by

As we discussed in [6], the system (3) is invariant under the phase shift for any constant . Hence we can identify states which are equal up to a phase shift, and reduce the dynamics to an -dimensional state space, which is the torus . In this reduced state space, there is a unique state with all equal, which we refer to as the sync state or just sync.

Our goal is to understand the generic long-term behavior of trajectories in the reduced state space, in forward and backward time. Some terminology: an asynchronous state has all distinct; an  state has all but one equal. The first result, Theorem 1, is that if then almost all trajectories in the reduced state space converge in forward time to sync, and in backward time to an asynchronous state or to one of finitely many  states. The dynamics are similar if , except reversed in time. In the course of preparing this manuscript we learned that this result was independently discovered by M. A. Lohe in [7], Lohe presents an argument that is essentially correct, but has some subtle technical gaps which we address in the discussion following the proof of Theorem 1. Our approach, which is different than in [7], is based on the correspondence between the system dynamics and a flow on the hyperbolic disc with special attention to the behavior at the boundary circle. The geometric techniques we develop in preparation for the proof of Theorem 1 also form the basis of the proofs of our subsequent Theorems 2-5.

The case is covered in Theorem 2; then the system has Hamiltonian structure, and almost all trajectories in the reduced state space are periodic or homoclinic connections to and from sync, with one exceptional case, which is if at least one of the coefficients satisfies the condition . Then in addition to the behavior above, there also exists a positive measure set of initial conditions with trajectories that converge in forward time to  states (and also a positive measure set of initial conditions with trajectories that converge in backward time to  states). We will also discuss three special cases of this system, where we can describe additional details of the dynamics: the case , where the dynamics have both gradient and Hamiltonian structure; the case of real with , where the dynamics flow along a 2D electrostatic field, and the case of real , which includes the classic Kuramoto model with all . The results are summarized in Table 1 in the Discussion section at the end of this paper.

The organization of this paper is as follows: We begin by summarizing some of our earlier work on Kuramoto networks in [6], which exploits a connection to hyperbolic geometry to simplify the analysis of the network dynamics. We will then derive some general properties of the system (3), and prove Theorems 1 and 2. After this we proceed with the analysis of the special cases of (3) mentioned above, and conclude with some suggestions for future research.

Ii Complex Formulation

It is desirable to express the system (1) in complex form, with . Let ; is a complex-valued function on which we define to be the order parameter for the system. Then using we obtain governing equations

(5)

The asymmetric Kuramoto-Sakaguchi system (3) has and order parameter given by (4) above. We will assume henceforth that the function and that the order parameter has this form.

Iii Reduction To 3D System

In their seminal paper [8], Watanabe and Strogatz demonstrated that the trajectories for any system of the form (1) are constrained to lie in submanifolds of the state space with dimension at most three. Subsequently, it was shown that these submanifolds are the group orbits under a natural action of the Möbius group on the torus [9]. Here is the 3D group of Möbius transformations which preserve the unit disc (and hence its boundary ). An element can be expressed uniquely in the form

(6)

where the parameters and satisfy and . When , we denote the above Möbius transformation by . If and then

defines the group action of on . The group orbits are the sets .

Now fix a point , and assume that at least three of the are distinct. Then any point in the group orbit can be expressed in the form for a unique choice of and . In effect, and can be thought of as coordinates on the group orbit . As derived in [6], the system (5) on is equivalent to the system in and given by

(7)

Iv Reduction To 2D System

When the order parameter function has the form (4), we can cancel the and in the equation above, which then simplifies to an equation in alone:

(8)

The variable determines a point on the group orbit up to rotation by some ; in effect, determines the phase relations among the coordinates .

More formally, if we identify and for any , then the full state space for this reduced model is the -dimensional torus ; the group orbits under this identification give us reduced group orbits , which are invariant under the reduced dynamics. In this reduced state space, sync is the unique fully synchronized state represented by any . For a base point with at least three distinct coordinates, its reduced -orbit is parametrized by , and equation (8) gives the dynamics on the reduced orbit. Fixed points in the reduced system correspond to either fixed points or uniformly rotating solutions with constant phases in the original -dimensional system. Since we are primarily interested in how the phases among the coordinates evolve, we will henceforth make this reduction and work with the 2D dynamical system given by (8).

Observe that changing the signs of all the in 8) is equivalent to reversing the direction of time for the system; this time-reversal property will be used frequently in the sequel.

V Boundary Correspondence

Fix a base point with all distinct; then parametrizes all possible phase configurations in the reduced group orbit . We wish to describe what happens to these phase configurations as approaches the boundary of the disc. Suppose a sequence converges to some , and for all . Then

which corresponds to the sync state in the reduced state space . However, if for some , then need not exist. To see this, suppose for example . Write , with , , so is the angle at which is approaching the boundary circle. Then

Therefore exists iff exists, and in this case

Thus we see that if for some , and the approach angle is asymptotically , then the limiting configuration of in the reduced state space is the  state in which the th oscillator has phase relative to the others in sync.

If all the coordinates are distinct, then this analysis shows that the boundary of the reduced -orbit consists of all  states together with the sync state. This is topologically a union of circles, all meeting at the sync state. (If does not have all distinct , then the boundary of will consist of circles meeting at sync, where is the number of distinct , provided . If then the boundary is sync.) In the next section we describe the dynamics on these boundary components, which are invariant for any system with order parameter given by (4).

Figure 1: Points in correspond to asynchronous states in . As , , the corresponding states in  approach sync; this also holds if tangentially. As at a fixed non-tangential angle , the corresponding states in  approach the  state with out of phase by .

The boundary correspondence is illustrated in Figure 1. A trajectory approaching a point on the boundary of the disc distinct from any (or approaching a tangent to the boundary circle) corresponds to a trajectory approaching sync in . A trajectory approaching a base point coordinate at a fixed angle corresponds a trajectory to approaching the  state with oscillator phase shifted by relative to the synchronized oscillators. We have shown previously [10] that any attracting or repelling states in  of a system of the form (5) must be sync or  states, which are the the common boundary of asynchronous -orbits. So trajectories approaching the boundary circle are particular particularly important to understanding the dynamics of (5).

Vi Dynamics

The boundary dynamics for the system (3) correspond to the dynamics for the system 5) with , which is

with and for constants . Since we are mainly interested in the phase between and , we let and analyze the dynamics of . Note that if , then and depend only on the constant , and must be equal; otherwise we can see from the above equations that the phase between and would change in time. So solutions to correspond to solutions which are rotating at the same constant angular velocity. The evolution equation for is

We see that , which corresponds to sync, is always a fixed point. If we express , then , and we obtain the equivalent equation

(9)

This flow always has fixed point , corresponding to sync, with eigenvalue . If , then there is an additional fixed point with opposite eigenvalue ; is stable if , unstable if , and has the opposite stability. If but is pure imaginary, then is the only fixed point, and is attracting globally but not locally near (the flow has the same direction everywhere on the circle). The flow is identically iff . We mention in passing that the uniformly rotating solutions corresponding to fixed states usually do not have angular velocity equal to . For example, the sync solution has angular velocity , which is not equal to unless .

For any partition of into two disjoint nonempty sets, the set of states where for is a 1D manifold in the reduced state space , invariant under the dynamics for any system (3), and these 1D manifolds all meet at the sync state. The dynamics on these two-cluster manifolds are given by the polar equation above, where measures the phase difference between the two clusters; the appropriate values of the coefficients are found by summing the over each of the two clusters. The tangent directions to these 1D invariant manifolds at are eigenvectors for the linearization of the system at sync, and they span the full tangent space at sync. Therefore we see from the polar equation above that the unique eigenvalue for the linearized dynamics at sync is . Consequently we see that the sync state in is linearly stable for , unstable for and linearly neutral for . This is of course consistent with the much stronger result of Theorem 1, that sync is globally stable when and globally unstable when .

Vii Fixed Point Analysis

In this section we study the fixed points for the flow (8). We begin with a lemma which will be crucial to the proofs of all our theorems.

Lemma 1. Assume that and at least three . Then the flow on the disc given by (8) has at most fixed points.

Proof. The details are easier to carry out if we transform the system to an equivalent system on the upper half plane, via the Möbius transformations

which give a correspondence between the upper half plane and the disc . (Notice that correspond to respectively.) Then the flow transforms to

(10)

So we must prove that the equation

has at most solutions . Let us assume without loss of generality that all , and . Observe that

Let

which is a rational function in with degree . The equation for fixed points is equivalent to . Express

and is a polynomial in with degree at most . Since has poles at , and have no common factors. Note that

which shows that the degree of is exactly , unless and then the degree is at most . Therefore we are done if ; the equation can have at most roots in .

If , then WLOG we can set and rewrite the fixed point equation for in the form . Observe that

both and are degree monic polynomials, since . Therefore has degree at most . We also see that and have no common factors, since any common factor would also divide . Therefore is a rational function with degree . Let

which is a rational function in . Our fixed point equation is equivalent to for . Any fixed point for the map is also a fixed point for the second iterate , which is a rational function in with degree , and therefore has at most fixed points (a rational function of degree has at most fixed points as a map of the extended complex plane ). The map has fixed points at , since has a pole at . Therefore the number of fixed points is at most .

We remark that the bound is far from sharp; using the Lefschetz fixed point theorem and some other results from the theory of iterated rational maps [11], we can improve the bound to (which coincidentally agrees with the previous bound for ). Since we only need the finiteness of the number of fixed points, we omit the proof, which would take us somewhat far afield.

The condition that at least three is necessary to insure finitely many fixed points; to see this, suppose only and all other . Then the fixed points for (8) in the disc are given by an equation of the form

for some nonzero . Since for any and , we see that we must have to have any solutions. We also must have to have solutions, since the map is one-to-one. For , the equation above is the equation of a circular arc in joining the points and . When , this arc is the unique geodesic joining and for the hyperbolic metric on the disc, which we will discuss later. For other , these circular arcs form the family of curves of constant curvature joining and , which are called hypercircles in hyperbolic geometry.

Next, we study the linear stability of fixed points for the system (8). The choice of coordinate depends on the base point , which without loss of generality we can choose so that is a fixed point for the flow in the reduced phase space . With this choice the flow (8) has fixed point at , and so . To first order in , the flow is given by

If we write , and as usual , then the 2D linear system for has matrix

which has , , and eigenvalues

(11)

Now suppose ; then the fixed point must have at least one eigenvalue with , and hence is a repelling node or spiral, a saddle, or a non-hyperbolic fixed point with one positive and one zero eigenvalue. In the first two cases there are respectively or trajectories that converge to the fixed point as . In the non-hyperbolic case with one , if the fixed point at is isolated (which is the case if at least three ) then there are at most two trajectories as (see [12] Section 2.11, Theorem 1). Therefore assuming at least three , we can conclude that in all cases there are at most two trajectories that converge to the fixed point as .

When the eigenvalues have the form with either real or pure imaginary, so fixed points can never be attracting in this case, and as we shall prove in the discussion preceding Theorem 2, can only attract finitely many trajectories. The case is equivalent to the case with time reversed. The eigenvalues at the fixed point are completely determined by the quantities and . More generally, if is a fixed point for the flow (8), the eigenvalues at are given by (11) with

Now we consider the question of constraints on the number and type of fixed points for (8). The equations for to be a fixed point for (8) with prescribed value are

Suppose we fix distinct and for , and also fix . We wish to find coefficients such that and (8) has fixed points at with . This is a system of linear equations in the coefficients :

We claim this system has solutions if . To prove this, consider the associated homogeneous system with . If we transform to the upper half plane, as we did in the proof of Lemma 2, with , then the homogeneous system in the is equivalent to

Using the identity

together with , we see that the homogeneous system is equivalent to the system

This implies that the rational function in Lemma 2 has for , which implies that its numerator has double roots at the . If is not identically , then we must have , which is a contradiction. Hence is identically , which means that all . Therefore the inhomogeneous system has a unique solution for , and infinitely many solutions for . In other words, we can find systems of the form (8) with as many fixed points as we desire, and can even prescribe the eigenvalues at the fixed points as we like, within the constraints imposed by the form of the eigenvalue equation (11).

Viii Boundary Flow Analysis

To fully understand the dynamics of (8), we need to analyze the flow near the boundary of the disc . Consider first any point which is distinct from any of the . Near , the trajectories are the same curves as for the modified flow without the factor , given by

This modified flow extends to a smooth flow on any open subset of which excludes the . If , then

So we see that if , then this flow points outwards and crosses the circle at every point . Hence the original flow will have a unique trajectory converging in forward time to each . If is a point on this trajectory, then the forward limit set must be the single point . Similarly if , then there is a unique trajectory converging in backward time to each . We also see that if but , then no trajectory can converge in forward or backward time to any , since the modified flow has trajectories along the arcs of the circle obtained by removing the .

Next, we analyze the flow near the . As in the proof of Lemma 1, it is easier to transform the system to the upper half plane and study the equation (10). Assume WLOG that , so and the remaining . We wish to analyze the flow near , and to do this we will employ the polar representation , where . We see that

(12)

where represents sums of terms of the form with . Using the polar conversions

we get the equivalent polar system

(13)

When , the equation for can be written in the form

which is equivalent to (9) if we let and . This makes sense, because the points with correspond to the points on the  boundary component with all oscillators in sync except the first. This equation has a unique fixed point , defined by

provided . (There are no solutions in if , unless we also have ; in this case, for all .) As we saw earlier, when , the fixed point has eigenvalue and the fixed point has eigenvalue , in the direction along the interval with . The linearization of the equation at is

(14)

which is independent of .

Now suppose . Then the fixed point has positive eigenvalue , so must be a repelling node, a saddle, or a non-hyperbolic fixed point with exactly one non-zero eigenvalue. If we assume that at least three , then this fixed point is isolated. If it is a repelling node then no trajectory converges to the fixed point in forward time. A saddle has two attracting trajectories, but they must be on opposite sides of the repelling trajectories (unstable manifolds) along the axis, so there is one attracting trajectory with . If the fixed point has a zero eigenvalue, it must be a repelling node, topological saddle or saddle-node (see [12] Section 2.11, Theorem 1), and in each case has at most one attracting trajectory with . Therefore in all cases there is at most one trajectory converging to this fixed point in forward time. In backward time, a set of initial conditions with positive measure will converge to this fixed point if it is a repelling node; no trajectories converge to this fixed point if it is a saddle or topological saddle (the unstable manifolds have ), and at most a single trajectory with can converge to the fixed point if it is a saddle-node. These results will be crucial for the proof of Theorem 1.

Ix Hyperbolic Geometry and the Gradient and Hamiltonian Conditions

The factor in the equation (8) suggests that this flow has connections to hyperbolic geometry. The Poincare model for hyperbolic geometry on the unit disc has metric

This metric is conformal with the Euclidean metric (i.e. angle measures agree), has constant negative curvature , and its geodesics are lines or arcs of circles which meet the boundary in angles. Since the reduced -orbits are in one-to-one correspondence with via the coordinate , we can transfer this metric to the reduced -orbits. In fact, as shown in [6], this metric on the reduced -orbits is independent of the choice of base point.

In 2D Riemannian geometry the simplest flows are given by gradient and Hamiltonian vector fields. A gradient vector field has the form for some smooth real function , where the gradient is defined in terms of the Riemannian metric. A Hamiltonian flow is the rotation of a gradient field , and therefore has as a conserved quantity. The hyperbolic gradient of a function in complex form is given by

where is the hyperbolic metric factor, is the ordinary Euclidean gradient and

In Ref [6] we derived criteria for the flow to be gradient or Hamiltonian: define the differential operator on the torus with coordinates by

Then the flow is gradient for the hyperbolic metric on all reduced orbits iff everywhere on , and similarly is Hamiltonian iff everywhere on . For the asymmetric Kuramoto-Sakaguchi model with order parameter (4), these conditions reduce to

In particular, we see that the symmetric K-S model, which has all , is gradient iff or , and Hamiltonian iff (as first pointed out in Ref. [8]).

We showed in [6] that the flow for the symmetric K-S model with is the hyperbolic gradient flow for the function

It is not hard to modify this function to find the corresponding potential for the asymmetric case, assuming is real. Observe first that

We will need the identity

which follows from

Using this, we see that