We study a system of N identical interacting particles moving on the unit sphere in d-dimensional space. The particles are self-propelled and coupled all to all, and their motion is heavily overdamped. For d=2, the system reduces to the classic Kuramoto model of coupled oscillators; for d=3, it has been proposed to describe the orientation dynamics of swarms of drones or other entities moving about in three-dimensional space. Here, we use group theory to explain the recent discovery that the model shows low-dimensional dynamics for all N3 and to clarify why it admits the analog of the Ott–Antonsen ansatz in the continuum limit N. The underlying reason is that the system is intimately connected to the natural hyperbolic geometry on the unit ball Bd. In this geometry, the isometries form a Lie group consisting of higher-dimensional generalizations of the Möbius transformations used in complex analysis. Once these connections are realized, the reduced dynamics and the generalized Ott–Antonsen ansatz follow immediately. This framework also reveals the seamless connection between the finite and infinite-N cases. Finally, we show that special forms of coupling yield gradient dynamics with respect to the hyperbolic metric and use that fact to obtain global stability results about convergence to the synchronized state.

Exactly solvable models have long played a central role in nonlinear dynamics, from Newton’s work on the gravitational two-body problem to breakthroughs in understanding solitons in the 1970s. Often, the solvability of a model reflects an underlying symmetry or other special structure in its governing equations. In this paper, we discuss a many-body system of current interest, known as the Kuramoto model on a sphere, whose unexpectedly low-dimensional dynamics call out for explanation. The model consists of N identical overdamped particles moving on a (d1)-dimensional sphere in d-dimensional Euclidean space, yielding a state space of dimension N(d1). Yet, despite the presence of damping, the model exhibits enormously many constants of motion. Here, we show that its trajectories are confined to invariant manifolds of dimension d(d+1)/2 for all N3 and trace the origin of this low-dimensional behavior to an underlying group-theoretic structure in the system. Specifically, the Kuramoto model on a sphere turns out to be the flow induced by the action of the group of Möbius transformations on the d-dimensional ball, and its invariant manifolds are the associated group orbits. For certain forms of coupling, the model acquires a further structure (hyperbolic gradient dynamics) that forces almost all solutions to converge to the perfectly synchronized state.

In 1975, Kuramoto introduced a model for a large population of coupled oscillators with randomly distributed natural frequencies.1 Kuramoto’s model displayed many remarkable features: It was exactly solvable (at least in some sense) despite being nonlinear and infinite-dimensional.2 Its solution shed analytical light on a phase transition to mutual synchronization that Winfree had previously discovered in a similar but less convenient system of oscillators.3,4 Since then, the Kuramoto model has been an object of fascination for nonlinear dynamicists, as well as a simplified model for many real-world instances of coupled oscillators in physics, biology, chemistry, and engineering.5–12 

From a mathematical standpoint, one of the most intriguing problems has been to explain the tractability of the Kuramoto model. What symmetry or other hidden structure accounts for its solvability?

The first clues came from work on an adjacent topic: series arrays of N identical overdamped Josephson junctions. The governing equations for those superconducting oscillators are closely related to the equations of the Kuramoto model13,14 and themselves displayed remarkable dynamical features, such as surprisingly low-dimensional invariant tori15,16 and ubiquitous neutral stability of splay states17 despite the presence of damping and driving in the governing equations. These features were explained in 1993 by the discovery of a certain change of variables, now called the Watanabe–Strogatz transformation,18,19 which showed that the governing equations have N3 constants of motion for all N3. Goebel20 then pointed out that the Watanabe–Strogatz transformation could be viewed as a time-dependent version of a linear fractional transformation, a standard tool in complex analysis. For more than a decade, however, these results did not attract much attention perhaps because they were assumed to be restricted to problems about Josephson junctions and within that specialized setting, they were even further restricted to junctions that were strictly identical.

A breakthrough occurred in 2008 with the work of Ott and Antonsen.21,22 They found an astonishing way to capture the macroscopic dynamics of the infinite-N Kuramoto model even when the oscillators’ frequencies were non-identical and randomly distributed. First, they wrote down an ansatz—seemingly pulled out of thin air—for the density function ρ(θ,ω,t) of oscillators having phase θ and intrinsic frequency ω at time t. Their ansatz had the form of a time-dependent Poisson density (a density better known for its role in the study of partial differential equations, specifically for the solution of Laplace’s equation on a disk, given the values of the unknown function on the bounding circle). By making this ansatz of a Poisson density, Ott and Antonsen reduced the infinite-N Kuramoto model, an integro-partial differential equation, to an infinite set of coupled ordinary differential equations. Then, by further assuming that the intrinsic frequencies of the oscillators were randomly distributed according to a Lorentzian (Cauchy) distribution, Ott and Antonsen showed that the order parameter dynamics of the Kuramoto model could be reduced tremendously, all the way down to an ordinary differential equation for a single scalar variable, the amplitude of the order parameter.21 With this discovery, the floodgates were now open. Almost immediately, the Ott–Antonsen ansatz was used to solve many longstanding problems about the Kuramoto model and its variants, as well as to generate and solve many new problems.10 

Still, a lot of old questions hung in the air. Both the Watanabe–Strogatz transformation and the Ott–Antonsen ansatz appeared somewhat unmotivated and almost miraculous. Where did they come from, and why did they work? It was also not clear whether they were connected or perhaps even equivalent. There were reasons to doubt that they were linked: the Watanabe–Strogatz transformation could be used for any finite N3 but seemed restricted to identical oscillators, whereas the Ott–Antonsen ansatz allowed for non-identical oscillators but seemed restricted to the continuum limit of infinite N. Also, why were linear fractional transformations and Poisson densities—tools from other branches of mathematics—popping up in these studies of dynamical systems?

Later work made sense of all of this. The Josephson arrays and the Kuramoto model both turned out to have deep mathematical ties to group theory, hyperbolic geometry, and projective geometry, and both the Watanabe–Strogatz transformation and the Ott–Antonsen ansatz were tapping into these structures.10,23–27 For the Josephson arrays, the governing equations turned out to be generated by a group action, specifically the action of the Möbius group of linear fractional transformations of the unit disk to itself. Seen in this light, the constants of motion for the Josephson arrays were cross-ratios, and the invariant tori were group orbits. The same group-theoretic structure was found to underlie the Kuramoto model (in the special case where all the oscillator frequencies are identical) as well as other sinusoidally coupled systems of identical phase oscillators.24,26

In the past few years, several researchers wondered how far this story could be pushed. Are there quantum or higher-dimensional extensions of the Kuramoto model that might show similar reducibility? A number of results along these lines have now been found.28–45 In particular, several researchers have explored a generalization of the Kuramoto model in which the oscillators move on spheres instead of the unit circle; the spheres could be either the ordinary two-dimensional sphere or higher-dimensional spheres. In particular, a system of particles moving on the two-sphere has been used to model the orientation dynamics of swarms of drones flying around in three dimensions.46 As we shall demonstrate, these higher-dimensional oscillator models exhibit a dynamical reduction generalizing the reduction for the Kuramoto model. Consequently, the computational effort required to simulate these systems can be substantially reduced. A counterpart of the Ott–Antonsen ansatz has also been discovered for the continuum version of the Kuramoto model for identical oscillators on the d-dimensional sphere and used to reduce its infinite-dimensional dynamics to a lower dimensional set of ordinary differential equations (ODEs).38 However, as before, some of the results appear disconnected and a bit miraculous.

Our goal in this paper is to show that hyperbolic geometry and group theory can unify and clarify our understanding of the Kuramoto model on a sphere and make all the latest results seem natural, just as they did before for the traditional Kuramoto model. We focus exclusively on the case of identical oscillators, though we expect that our methods will extend to systems with multiple populations of identical oscillators or a continuum of oscillator families with a distribution of natural frequencies. Our approach explains the model’s reducibility for any finite number of oscillators, as well as for the continuum limit, and it reveals why Poisson densities arise again in this setting. There is a close connection to Laplace’s equation and harmonic analysis, as we will see in Sec. V. We also find that complex analysis is not really essential, which is just as well, since it does not generalize to the higher-dimensional spheres being considered here. On the contrary, the proper mathematical setting is harmonic analysis and hyperbolic geometry on higher-dimensional balls. Our work also allows us to go beyond merely unifying existing results. For instance, by establishing that linearly coupled systems of identical Kuramoto oscillators on a sphere have a hyperbolic gradient structure, we can prove new global stability results about convergence to the synchronized state, as described in Sec. VIII.

In a pioneering paper, Lohe28 observed that there are at least two natural generalizations of the Kuramoto model to higher dimensions. One of them replaces the phases θj of the original Kuramoto model1,2 with complex numbers exp(iθj) on the unit circle and then views those as equivalent to 2×2 rotation matrices parameterized by a rotation angle θj. From there, it is a natural step to consider other Lie groups of matrices, many of which are non-Abelian.

Our concern in this paper, however, is with a different generalization of the Kuramoto model. Instead of regarding oscillators as particles moving on the unit circle, we think of them as particles moving on the unit sphere. The sphere could be the surface of the ordinary unit ball in three dimensions or some higher-dimensional sphere Sd1 in Rd. When d=2, the sphere reduces to the unit circle in the plane, and the model reduces to the original Kuramoto model.

The governing equations for the Kuramoto model on a sphere are

x˙i=Aixi+ZZ,xixi,i=1,,N,
(1)

where xi is a point on the unit sphere Sd1Rd, each Ai is an antisymmetric d×d matrix, and ZRd is a d-dimensional vector analogous to the complex order parameter for the classic Kuramoto model. In Eq. (1), the matrix Ai and the vector Z are functions of the configuration (x1,,xN) of points on the sphere. Note that Z does not depend on i; like the usual Kuramoto order parameter, it plays the role of a mean-field quantity that couples all the “oscillators” xi together. The antisymmetric matrix Ai is the higher-dimensional counterpart of an intrinsic frequency ωi in the original Kuramoto model.

A straightforward computation shows that the dot product between an oscillator’s instantaneous position and instantaneous velocity satisfies xi,xi˙=0, which proves that oscillators that start on the unit sphere stay on it forever. The state space for this system is the N-fold product X=(Sd1)N, which has dimension N(d1). Later, we will also consider the natural infinite-N analog of (1), where a state is a probability measure on Sd1.

In what follows, we allow Z to be any smooth function on the state space X, though in examples, we usually restrict to fairly simple functions, such as a linear combination of the form

Z=i=1Naixi,

where the ai are real constants.

There is a general technique for dimensional reduction of systems such as (1) that we pause to describe. Suppose we have a smooth manifold X, which we think of as a state space, and a group G acting on X, where G is also a smooth manifold (in other words, G is a Lie group). Then, the group action induces a space of vector fields on X, the so-called infinitesimal generators of the action.

To construct these generators, let γ(t) be a smooth curve in G with γ(0)=e, the identity element in G. Then, the derivative γ˙(0)=v, where v is a vector in the tangent space TeG of G at e. This vector v is in turn associated very naturally with a corresponding vector v~ in the tangent space of X, as follows. For each xX, tγ(t)x is a smooth curve in X, and its derivative at t=0 defines a vector v~x in the tangent space at x. The vector field v~=(v~x) is the infinitesimal generator corresponding to the element v in the tangent space of G at e. At each point xX, the infinitesimal generators vx span a linear subspace Vx of the tangent space TxX, which is exactly the set of vectors tangent to the group orbit Gx at the element x.

Now, suppose we have a vector field ξ on X, which defines a dynamical system on the state space X. If ξxVx at each point xX, then the flow corresponding to the vector field ξ will be constrained to lie on the group orbits Gx. If the dimension of G is less than that of the state space X, then this will give us a dimensional reduction of the dynamics from dimX to dimG. Now, suppose, as is the case in applications of this methodology, the correspondence GGx is one-to-one for generic xX; equivalently, the stabilizer subgroup Gx={e} for generic xX. Then, for each xX, the flow on the group orbit Gx is equivalent to a flow on G, which is a lower dimensional space than the state space X.

Here is a familiar example: Consider the orthogonal group G=SO(d) consisting of orientation-preserving linear isometries of Rd. (For an intuitive picture, think of these isometries as rotations.) Then, G acts on Rd, and the corresponding infinitesimal generators are the linear vector fields νx=Ax, where A is any skew-symmetric matrix. We can also let G act on the product space X=(Rd)N of N-tuples x=(xi), xiRd, and then the infinitesimal generators have the form (νx)i=Axi for some skew-symmetric matrix A. We could also, if we like, restrict X to N-tuples (xi) with xiSd1, the state space for (1). Now, suppose we had a dynamical system on X of the form x˙i=Axi, where the skew-symmetric matrix A is a function of the configuration x=(xi); this is just the special case of (1) with Z=0. Then, the dynamics on X reduces to dynamics on G, which has dimension d(d1)/2, and for large N, this is much smaller than the dimension of X, which is N(d1). Basically, the configuration (xi) of points on the sphere Sd1 can only move collectively by a rotation of the sphere; therefore, the dynamics reduces to a dynamical system on SO(d).

We want to apply this methodology to the system (1). However, since that system generally has Z0, we need a different group action to make this strategy work. Fortunately, vector fields of the form seen in the Kuramoto model,

x˙=Ax+ZZ,xx,
(2)

turn out to arise as the infinitesimal generators of the group action of a larger group G acting on the sphere Sd1 and its interior, the unit ball Bd. This larger group is the Möbius group of isometries of the hyperbolic geometry on Bd. It contains the orthogonal group as a proper subgroup but has bigger dimension d(d+1)/2.

Therefore, for the Kuramoto system (1), if the matrices Ai all happen to be identical, we can reduce the dynamics of the system to a much smaller system on this Möbius group G and use this reduction to understand the dynamics on the larger state space X. This is the essence of the approach we take. Ultimately, we apply it to prove a synchronization theorem for the system (1) for the special order parameter Z=i=1Naixi with ai>0. However, first, we need to show that vector fields on Sd1 of the form in (2) are indeed the infinitesimal generators of the action of the larger Möbius group.

In this paper, a Möbius transformation is a composition of Euclidean isometries and spherical inversions of Rd mapping the unit ball homeomorphically to itself and preserving orientation. This is a more restrictive definition than the commonly defined Möbius transformations, which, in general, do not need to preserve the unit ball.

As in the case d=2, flows of the form (1) are intimately related to the natural hyperbolic geometry on the unit ball Bd with boundary Sd1. This geometry has metric

ds=2|dx|1|x|2,

where |dx| is the ordinary Euclidean metric. Isometries are assumed to be with respect to this hyperbolic geometry unless otherwise qualified as Euclidean. The metric ds has a constant (sectional) curvature equal to 1, and we can describe its isometries, which generalize the Möbius transformations preserving the unit disk for d=2. For d=2, let wB2 and consider the Möbius transformation

Mw(x)=xw1w¯x,

which preserves the unit disk B2 and its boundary S1. To generalize this to higher dimensions, we need to express Mw(x) without reference to complex arithmetic operations or conjugation. This is the goal of Subsection II C 1.

1. Möbius transformations in higher dimensions

Using the identity 2w,x=w¯x+wx¯, we see that

(xw)(1wx¯)(1w¯x)(1wx¯)=xww|x|2+w2x¯12w,x+|w|2|x|2=xww|x|2+w(2w,xw¯x)12w,x+|w|2|x|2=(1|w|2)x(12w,x+|x|2)w12w,x+|w|2|x|2.

This form of Mw generalizes to higher dimensions: Let wBd and define

Mw(x)=(1|w|2)x(12w,x+|x|2)w12w,x+|w|2|x|2=(1|w|2)(x|x|2w)12w,x+|w|2|x|2w,

where xBd or Sd1. We call Mw a boost transformation.

If |x|=1, this formula simplifies to

Mw(x)=(1|w|2)(xw)|xw|2w.

Now, see that

Mw(w)=(1|w|2)w(12w,w+|w|2)w12w,w+|w|2|w|2=(1|w|2)w(1|w|2)w12w,w+|w|2|w|2=0.

Alternatively, we can use the second form to show

Mw(w)=(1|w|2)(w|w|2w)12w,w+|w|2|w|2w=(1|w|2)2w(1|w|2)2w=0.

Similar computations show that M0 is the identity, Mw1=Mw, and Mw(0)=w.

It is a standard result in hyperbolic geometry (e.g., see Theorem 3.5.1 in Beardon47) that any orientation-preserving isometry of Bd can be expressed uniquely as the product of a boost and a rotation (an orientation-preserving orthogonal transformation), and these operations can be done in either order. In other words, any such isometry can be written uniquely in the form

g(x)=ζMw(x)

and also uniquely in the form

g(x)=Mz(ξx)

for some vectors w,zBd and rotations ζ,ξSO(d), where SO(d) denotes the group of orientation-preserving orthogonal linear transformations on Rd. Therefore, counting the extra d dimensions that we get from the vector w or z, we see that the Möbius group has dimension d+d(d1)/2=d(d+1)/2.

Depending on the situation, one of these two forms might be more useful than the other, even though they are equivalent. In the interest of flexibility, it is useful to find how the parameter pairs w,ζ and z,ξ are related. We can do this by comparing the linearizations of the two formulas for g(x) at x=0,

g(x)ζ(w+(1|w|2)x)z+(1|z|2)ξx.

Equating coefficients implies z=ζw (hence, |z|=|w|) and ξ=ζ.

2. Infinitesimal generators

Having parameterized the Möbius transformations, we are now ready to derive the associated infinitesimal generators of the Möbius group action on the ball Bd. We will show that they correspond to flows of the form

y˙=AyZ,yy+12(1+|y|2)Z,
(3)

where A is an antisymmetric d×d matrix and ZRd is a vector. Note that, as advertised, this flow extends to a flow on Sd1 of the Kuramoto form in (2), as we can see by restricting (3) to vectors y on the unit sphere where |y|=1.

To derive (3), we work separately with the boost and rotation components. Let us start with the boost component. Replace w by tw and expand Mtw(x) to first order in t,

Mtw(x)x|x|2tw12tw,xtwx+t(2w,xx(1+|x|2)w).

The derivative of this expression at t=0 (which is just the coefficient of t) gives us the infinitesimal generator: it is an “infinitesimal boost” of the form (3) with Z=2w and A=0. Next, recall that the infinitesimal generators corresponding to the rotation components are flows of the form x˙=Ax for an antisymmetric matrix A. Together with the infinitesimal boosts, we then get all flows of the form (3). The group G acts on the space X in the natural way (component by component), and the infinitesimal generators of this group action on X are flows of the form (1) with all Ai identical. Therefore, by the general philosophy discussed earlier, the evolution of any initial point pX under the system (1) with all Ai=A lies in the group orbit Gp.

The given Kuramoto system has N(d1) degrees of freedom for some large N. However, since the flow of the system is determined via an action of the d(d+1)/2 dimensional Lie group G, we can alternatively study the auxiliary dynamical system on G, which we call the reduced equations. By ignoring rotations, we can further restrict our attention to a system on the d-dimensional quotient G/SO(d)Bd. The dimensional reduction not only makes the reduced equations easier to analyze than the original Kuramoto system, but the reduced equations require fewer computational resources to numerically integrate.

As before, assume all the terms Ai in (1) are equal so that Ai=A for some skew-symmetric matrix A. Fix a base point p=(p1,,pN)X. Then, if the points pi are in sufficiently general position, every element in the G-orbit of p can be expressed uniquely as gp for some gG, with parameters w,z, and ζ. We wish to derive the corresponding evolution equations for w,z, and ζ. Let (xi(t)) be any solution to (1) in the group orbit Gp; we do not require that the initial point (xi(0))=p. Then, for i=1,,N, we have xi(t)=gt(pi) for a unique gtG, which determines the parameters w,z,ζ as functions of t. Now, consider Eq. (3), with coefficients A and Z evaluated at (xi(t)). This is a non-autonomous ODE on Bd¯, and its time-t flow must be given by some g~tG. This ODE has solutions xi(t)=gt(pi)=gt(g01(xi(0))), which implies that g~t=gtg01.

Therefore, for any y0Bd¯,

y(t)=gt(g01(g0(y0)))=gt(y0)=ζMw(y0)=Mz(ζy0)

must satisfy the ODE (3) with A and Z evaluated at (xi(t)) at time t. In particular, if we let y0=0, then y(t)=ζw=z; therefore, z satisfies the ODE (3).

Now, expand y=ζMw(y0)=Mz(ζy0) to first order in y0, using the variables z and ζ,

yz+(1|z|2)ζy0;

therefore,

y˙z˙2z˙,zζy0+(1|z|2)ζ˙y0.

On the other hand, (3) gives

y˙=Ay+12(1+|y|2)ZZ,yyAz+12(1+|z|2)ZZ,zz+(1|z|2)(Aζy0+z,ζy0ZZ,zζy0Z,ζy0z).

Setting y0=0 gives the z˙ equation

z˙=Az+12(1+|z|2)ZZ,zz
(4)

as expected, and since Az,z=0, this in turn implies that

z˙,z=12(1|z|2)Z,z.

Equating the y0 terms, factoring out 1|z|2, and canceling the common term Z,zζy0 gives

ζ˙y0=Aζx0+z,ζy0ZZ,ζy0z.

Together, the last two terms above define a special type of antisymmetric operator of ζy0: Given any y1,y2Rd, define the antisymmetric operator α as

α(y1,y2)y=y1,yy2y2,yy1;

this operator has range = span(y1,y2) provided that y1 and y2 are linearly independent; otherwise, α(y1,y2)=0. Then, for all y0Rd,

ζ˙y0=Aζy0+α(z,Z)ζy0;

therefore,

ζ˙=(A+α(z,Z))ζ.

Differentiating z=ζw gives

Az+12(1+|z|2)ZZ,zz=ζw˙ζ˙w;

therefore,

ζw˙=(A+α(z,Z))zAz12(1+|z|2)Z+Z,zz=Az+|z|2ZZ,zzAz12(1+|z|2)Z+Z,zz=12(1|z|2)Z;

hence,

w˙=12(1|w|2)ζ1Z.
(5)

Summing up, the evolution equations for the (z,ζ) coordinate system on Gp are

z˙=Az+12(1+|z|2)ZZ,zz,
(6a)
ζ˙=(A+α(z,Z))ζ,
(6b)

with A and Z evaluated at Mz(ζp) and for the (w,ζ) coordinate system on Gp are

w˙=12(1|w|2)ζ1Z,
(7a)
ζ˙=(Aα(ζw,Z))ζ,
(7b)

with A and Z evaluated at ζMw(p). Note that these equations generalize the evolution equations for the parameters w and ζ given in Chen et al.26 for the classic case d=2.

The z˙ equation (4) is an extension of the system equation (1) on Sd1. However, for finite N, the z˙ equation does not uncouple from ζ since Z is evaluated at Mz(ζp). The exception to this is in the infinite-N limit: if the base point p is now the uniform density on Sd1, then ζp=p (the uniform density is invariant under rotations) and the density Mz(p) is a hyperbolic Poisson density on Sd1 whose centroid is a function of z. In the case d=2, this Poisson density has centroid z. Unfortunately, this simple relationship is false for d3 (we will give more details on this in Sec. V).

The advantage of the w˙ equation (5) is that for an order parameter function of the form

Z=i=1Naixi,

with aiR, ζ drops out of the w˙ equation, and we get the reduced equation

w˙=12(1|w|2)Z(Mw(p)).

The parameter w essentially defines the “phase relations” among xi; two configurations have the same w if and only if they are related by a rotation. Therefore, w is the key parameter that determines whether the system is approaching synchrony or incoherence.

The w variable also has a nice invariance under change of base points. Suppose p=M(p)Gp, then, we have coordinates w,ζ associated with the base point p. Any qGp has two expressions

q=ζMw(p)=ζMwp=ζMw(M(p)).

Assuming the coordinates of p are in sufficiently general position, this implies ζMw=ζMw°M, and hence,

0=ζMw(w)=ζMw(M(w)).

Therefore, the unique solution to Mw(y)=0 is w, and hence, w=M(w). In other words, the coordinates w and w transform exactly as the base points p and p.

Next, we consider the dynamics of the Kuramoto model (1) in the limit N. We assume that the rotation terms Ai=A are constant across the population, corresponding to identical “oscillators.”

Let us also assume that the order parameter Z is proportional to the centroid of the population,

Z=KNi=1Nxi.

In the continuum limit, a state of the system is a probability measure ρ on Sd1, and the order parameter becomes

Z=KSd1xdρ(x).

The measure ρ evolves according to the continuity equation (also known as the noiseless Fokker–Planck equation) associated with the flow in (1). Naturally, this flow must preserve group orbits under the action of G. Recall that if MG, then the measure Mρ is defined by the adjunction formula

Sd1f(x)d(Mρ)(x)=Sd1f(M(x))dρ(x).

In particular, we can consider the G-orbit of the uniform probability measure σ on Sd1. This orbit is special; whereas a typical group orbit Gρ has dimension equal to the dimension of G, namely, d(d+1)/2, the orbit Gσ has dimension only d. This is because the stabilizer of σ is SO(d); any rotation fixes σ, whereas the boosts deform σ. Hence, the orbit Gσ has dimension d. Any element in Gσ can be written as (Mz)σ, with zBd. The evolution equation for z is (4), with

Z(z)=KSd1xd(Mz)σ(x)=KSd1Mz(x)dσ(x).
(8)

In the case d=2 with x=ζS1, we have

dσ(ζ)=12πidζζ;

therefore, the integral

Z(z)=K2πiS1ζ+z1+z¯ζdζζ=Kζ+z1+z¯ζ|ζ=0=Kz

by the Cauchy integral formula. Therefore, (4) simplifies to the equation

z˙=iωz+K2(1|z|2)z

when d=2. Unfortunately, the formula Z(z)=Kz is not correct for d3; though as we shall see later, this formula is correct in higher dimensions for the complex hyperbolic model in even dimensions, which we discuss in Sec. VI. For d=2, the two geometries agree, which explains the coincidence for d=2.

Any Riemannian manifold X has a Laplace–Beltrami operator Δ associated with its metric; functions f on X satisfying the equation Δf=0 are called harmonic. For functions on the ball Bd with the hyperbolic metric, this operator is

Δhyp=(1|x|2)2Δeuc+2(d2)(1|x|2)i=1dxixi,

where

Δeuc=i=1d2xi2

is the standard Laplace operator (see Stoll,48 Chap. 3). We will call solutions to the equation Δhypf=0hyperbolic harmonic functions; for d=2, these coincide with ordinary (Euclidean) harmonic functions. We can consider the hyperbolic analog of the classical Dirichlet problem: given a continuous function f on Sd1, extend f to a hyperbolic harmonic function f~ on Bd. Assuming that this problem has a unique solution, then for any rotation ζSO(d), we must have f°ζ~=f~°ζ since rotations preserve the hyperbolic metric. If we average f°ζ on Sd1 over all rotations ζSO(d), we get the constant function

fave=Sd1f(x)dσ(x)

on Sd1, and any constant is hyperbolic harmonic on Bd. Therefore, the average on Bd of f°ζ~=f~°ζ over all ζSO(d) must be the constant fave. However, f~(ζ(0))=f~(0) for all ζ; therefore, we must have

f~(0)=Sd1f(x)dσ(x).

Now, let zBd; since Mz preserves the hyperbolic metric, we must have f°M~z=f~°Mz, which implies

f~(z)=f°M~z(0)=Sd1f(Mz(x))dσ(x)=Sd1f(x)d((Mz)σ)(x).

As shown in Chap. 5 of Stoll,48 the measure (Mz)σ is given by the formula

d((Mz)σ)(x)=Phyp(z,x)dσ(x),

with hyperbolic Poisson kernel function

Phyp(z,x)=(1|z|2|zx|2)d1.
(9)

Thus, the solution to the hyperbolic Dirichlet problem with boundary function f on Sd1 is given by the hyperbolic Poisson integral

f~(z)=Sd1Phyp(z,x)f(x)dσ(x),zBd.

The orbit Gσ consists of all hyperbolic Poisson measures P(z,x)dσ(x), parameterized by zBd. By contrast, the Euclidean Poisson kernel function is

Peuc(z,x)=1|z|2|zx|d;

therefore, the hyperbolic Poisson measures agree with the Euclidean Poisson measures only if d=2.

Now, we can calculate the expression Z(z) in the general case d2. We see from (8) that Z(z) is the hyperbolic Poisson integral of the function Kx on Sd1. The function Kx is (Euclidean) harmonic and homogeneous of degree 1 on Rd; following the recipe in Chap. 5 of Stoll,48 we see that its extension from Sd1 to a hyperbolic harmonic function on Bd is given by

Z(z)=KF(1,1d/2;1+d/2;|z|2)F(1,1d/2;1+d/2;1)z,
(10)

where F is the hypergeometric function

F(a,b;c;t)=k=0(a)k(b)k(c)ktkk!,

with (a)0=1 and (a)k=a(a+1)(a+k1) for k1. Notice that if a or b=0, then F(a,b;c;t)=1; this gives Z(z)=Kz for d=2, as expected.

There is an alternative generalization of Kuramoto networks to higher-dimensional oscillators when d=2m is even. Then, Rd=Cm, and we can study systems of the form

x˙j=Ajx+Zxj,Zxj,i=1,,N,
(11)

where now xi is a point on the unit sphere S2m1Cm, Ai is an anti-Hermitian m×m complex matrix, ZCm, and , denotes the complex-valued Hermitian inner product. These systems are the same as the real case when d=2,m=1 but are different for m2. To see this, suppose

Ax+Yx,YRx=Bx+Zx,ZCx

for all xS2m1Cm=Rd, where A is antisymmetric, B is anti-Hermitian, Y,ZCm, and we use the subscripts R and C to distinguish the real and complex inner products. Then,

(AB)x=ZY+(x,YRx,ZC)x,

and therefore, (AB)(x)=(AB)x for all xS2m1, which implies A=B. This implies

YZ=(x,YRx,ZC)x

for all xS2m1; hence, YZspanC(x) for all xCm; if m2, this implies Y=Z. However, then we have

x,YR=x,YC

for all xCm, which can only hold if Y=0. Hence, for m2, the only flows simultaneously of the form (1) and (11) have Z=0 and A anti-Hermitian.

Flows of the form (11) are related to the complex hyperbolic geometry on the complex unit ball Bm with the Bergman metric (see Rudin,49 Chap. 1). The orientation-preserving isometries of this metric are generated by unitary transformations ζU(m) and boost transformations of the form

Mw(x)=1|w|2x+(x,w1+1|w|21)w1x,w=xw+x,ww|w|2x1+1|w|21x,w.

Notice that when m=1, this reduces to the standard complex Möbius map Mw. As in the real case, M0 is the identity, Mw1=Mw, Mw(w)=0, and Mw(0)=w. Any orientation-preserving isometry of Bd can be expressed uniquely in the form

g(x)=ζMw(x)=Mz(ξx),

where w,zBm, but now ζ,ξU(m), the complex unitary group. Linearizing at x=0 gives

g(x)ζ(wx,ww+1|w|2x+x,ww1+1|w|2)ζ(w+1|w|2x1|w|2x,ww1+1|w|2)z+1|z|2ξx1|z|2ξx,zz1+1|z|2,

which implies z=ζw (hence |z|=|w|) and ξ=ζ, as before. The corresponding infinitesimal transformations are given by flows on the complex unit ball Bm of the form

y˙=Ay+Zy,Zy,
(12)

with A being anti-Hermitian m×m and ZCm. This flow extends to a flow on S2m1 of the form in (11). Note the absence of the quadratic term |y|2Z here. To derive these infinitesimal transformations, we can apply our prior power series expansion, noting that |x|=1 to obtain

Mtw(x)xtw1tx,wx+2(x,wxw)t.

Therefore, the infinitesimal generator is an “infinitesimal boost” of the form (12) with Z=12w and A=0. The infinitesimal generators corresponding to the rotation components are flows of the form x˙=Ax with A being anti-Hermitian; together with the infinitesimal boosts, we get all flows of the form (12).

We mention in passing that the complex model (11) has a natural quantum network analog, studied by Lohe,29 in which the complex vector xi is replaced by a normalized wave function |ψi. We expect that our reduction techniques extend to this infinite-dimensional quantum model.

Many of the results above can be found in some form in the work of earlier authors.28,30,31,34–39,43,44 Three papers, in particular, overlap considerably with the present work.

Tanaka30 demonstrates in his 2014 paper that the dynamics of (1) can be reduced using Möbius transformations that fix the unit ball, similar to what Marvel et al. found24 for the traditional Kuramoto model. Tanaka writes his Möbius transformations differently from ours, but he uses the same group of transformations and he also gets reduced equations for his Möbius parameters. Tanaka’s equation (10b) looks similar to our z˙ equation (4), except without the |z|2 term, which is puzzling. He does not mention the reduction down to dimension d in the finite-N case that we get with the w˙ equation (5). Tanaka also notes the complex case when d=2m is different and generalizes the Ott–Antonsen residue calculation to this case, which is the highlight of his paper. In the real case, Tanaka’s equation (15) is similar to our equation (10), though we were not able to show that the two expressions are equivalent. Finally, Tanaka also presents a generalization of the Ott–Antonsen reduction21 for the complex version of the system.

Lohe35 also looks at the same system as (1) [see his Eq. (22)], and he derives a similar reduction as ours by using Möbius transformations for the finite-N model. His transformation (30) on Sd1 is our Mw (with v=w) and his Eq. (31) is the same as our z˙ equation (4). He also has something that looks like the w˙ equation (5), which he says is independent of the rest of the reduced system for (in our notation) an order parameter function of the form

Z=1Ni=1NλiQixi,

where QiO(d) and λiR. However, such a Z does not satisfy the identity ζZ(p)=Z(ζp) for all rotations ζ unless Qi=±I; therefore, we do not see how the ζ term cancels.

Lohe’s map M in his Eq. (55) (ignoring the R factor) agrees with our map Mv on the sphere Sd1, but not on the ball Bd. Therefore, it is not a Möbius transformation of the type we are using. For example, M(v)=v, whereas Mv(v)=0. We are not sure why Lohe35 prefers these maps over the boosts; he claims that M preserves cross-ratios, but we do not see why this is advantageous. His map F in Eq. (63) (again ignoring the R factor) is exactly our Mv.

Chandra et al.38 concentrate on the infinite-N or continuum limit system and derive a dynamical reduction for a special class of probability densities on Sd1, generalizing the Poisson densities used in the Ott–Antonsen reduction. They proceed directly to the infinite-N version of (1). They make a very clever guess [their Eq. (7)] of the form of the special densities that generalize the Poisson densities for d=2 and then calculate the exponent in the denominator of their expression, getting exactly the hyperbolic Poisson kernel densities in (9) above. Their Eq. (15) is exactly the same as our z˙ equation (4) in the infinite-N limit. The integral in their Eq. (19) can be evaluated, as shown above in (10).

We conclude with an analysis of the system (1) with a weighted order parameter

Z=i=1Naixi,
(13)

where the ai are real constants.

We have been discussing a system of particles on a sphere that coalesce (in other words, they spontaneously synchronize) when certain conditions are met. To visualize this coalescence, we implemented the Runge–Kutta algorithm to numerically solve the Kuramoto model on the two-dimensional sphere in a three-dimensional space corresponding to the case d=3 in (1). For simplicity, we simulated N=100 particles with equal weights (ai=1/N in the order parameter Z) and set the rotation term A to zero for all the particles, a choice that is tantamount to ignoring the rotational influence, or equivalently, rotating the frame of reference along with the entire system as it evolves.

In the simulation shown in Fig. 1, randomly chosen points on the sphere were used as initial conditions. As time increases, one can see that the particles coalesce to a limit point, mimicking the spontaneous synchronization that is well known for the traditional Kuramoto model (d=2) when the oscillators are identical.

FIG. 1.

A first-order linear Kuramoto system on the two-dimensional sphere S2 with equal weights ai=1/N and randomly chosen initial conditions. The states shown are at t=0,t=10, and t=40, respectively. This simulation was written in Python and visualized with Plotly.

FIG. 1.

A first-order linear Kuramoto system on the two-dimensional sphere S2 with equal weights ai=1/N and randomly chosen initial conditions. The states shown are at t=0,t=10, and t=40, respectively. This simulation was written in Python and visualized with Plotly.

Close modal

Later in this section, we will prove that this synchronization behavior holds more generally for Kuramoto models on the sphere having weighted order parameters of the form (13), provided that the weights ai are all positive and sum to 1 and no individual weight exceeds 1/2. A partial result in this direction was obtained by Choi and Ha.32 These authors prove that in the case where all ai are equal and positive, initial conditions satisfying xi(0),Z(0)>0 for all i will synchronize. Geometrically, this condition is equivalent to requiring that the oscillators all lie on the hemisphere given by xi,y>0 for some vector y0. More generally, similar partial synchronization results are obtained by Ha et al.43 for the case

Z=(aI+W)i=1Nxi,

with W skew-symmetric and by Ha and Park41 for the complex system (11) with Z=xi.

Figure 2 shows a simulation in which we weighted each particle according to the terms in a Riemann sum approximating the integral of the normal probability distribution. The pink particles, which have higher weights, exert greater influence over the final synchronization location of the particles, but there is still synchronization.

FIG. 2.

A first-order linear Kuramoto system on S2 with weights distributed according to a Riemann sum, which approximates the integral of a normal distribution and randomly chosen initial conditions. Pink particles contribute to the order parameter with greater weights than the blue particles do. The states shown are at t=0,t=10, and t=40, respectively.

FIG. 2.

A first-order linear Kuramoto system on S2 with weights distributed according to a Riemann sum, which approximates the integral of a normal distribution and randomly chosen initial conditions. Pink particles contribute to the order parameter with greater weights than the blue particles do. The states shown are at t=0,t=10, and t=40, respectively.

Close modal

When time runs backward, almost all initial conditions of the particles tend toward a limiting configuration where their centroid is at the origin. The exception is when we have a majority cluster, depicted in Fig. 3, where one particle has a weight that exceeds the weight of all other particles. When this condition holds, it is impossible to arrange the particles so that their weighted centroid is at the origin, so the backward-time limit will tend toward an antipodal configuration, where all particles not in the majority cluster will coalesce around the antipode of the cluster. This is the configuration that minimizes the magnitude of the weighted centroid. We do not include the proof that the backward-time limit is antipodal in this case, but it is a straightforward generalization of the result that we do prove below.

FIG. 3.

A first-order linear Kuramoto system on S2 with a majority cluster, where one particle is chosen to have a weight that exceeds the combined weights of all other particles (or equivalently, where all the particles have equal weight, but a majority of them cluster into a single point and, therefore, act if they were a single giant particle, hence the name “majority cluster”). The states shown are at t=0,t=10, and t=40, respectively; we have chosen to depict time running backward to highlight that the backward-time limit tends toward an antipodal configuration. In this simulation, one particle, depicted in pink, was chosen to have a weight of 0.6, and the remaining 99 particles, depicted as blue, were chosen to have equal weights of 0.4/99.

FIG. 3.

A first-order linear Kuramoto system on S2 with a majority cluster, where one particle is chosen to have a weight that exceeds the combined weights of all other particles (or equivalently, where all the particles have equal weight, but a majority of them cluster into a single point and, therefore, act if they were a single giant particle, hence the name “majority cluster”). The states shown are at t=0,t=10, and t=40, respectively; we have chosen to depict time running backward to highlight that the backward-time limit tends toward an antipodal configuration. In this simulation, one particle, depicted in pink, was chosen to have a weight of 0.6, and the remaining 99 particles, depicted as blue, were chosen to have equal weights of 0.4/99.

Close modal

As mentioned above, if Z has the form in (13), then the w˙ equation in (5) reduces to

w˙=12(1|w|2)Z(Mw(p))
(14)

independent of the parameter ζ. We will show that this is a gradient flow on the unit ball Bd with respect to the hyperbolic metric. In the presence of a Riemannian metric, we can associate a 1-form to any vector field, and the vector field is gradient if and only if the associated 1-form is exact; since the unit ball is simply connected, this holds if and only if the associated 1-form is closed. For the Euclidean metric on Bd (or any open subset of Rd) and standard coordinates w1,,wd, the 1-form associated with the vector field with components f1,,fd is

ω=f1dw1++fddwd.

If we scale the Euclidean metric by a positive smooth function ϕ, then the associated 1-form with respect to the metric ds=ϕ|dw| is

ω=ϕ2(f1dw1++fddwd).

Therefore, the gradient of a function Φ with respect to this scaled metric is given by

Φ=ϕ2eucΦ,

where euc denotes the ordinary Euclidean gradient operator. We have ϕ(w)=2(1|w|2)1 for the hyperbolic metric; therefore, the hyperbolic gradient operator on Bd is given by

hypΦ(w)=14(1|w|2)2eucΦ(w).

Now, let us consider the vector field V defined by (14). By linearity, it suffices to treat the case Z=xi, and we can take i=1 without loss of generality. Then, the associated 1-form is

ω=4(1|w|2)2(12(1|w|2))j=1d((1|w|2)(p1,jwj)|p1w|2wj)dwj=2j=1d(p1,jwj|p1w|2wj1|w|2)dwj,

where p1,j denotes the jth component of the point p1Sd1. Let Ej denote the coefficient of dwj in parentheses above; then,

dω=2j,k=1dEjwkdwkdwj.

Applying the chain and quotient rules gives

Ejwk=2(p1,jwj)(p1,kwk)|p1w|4+2wjwk(1|w|2)2

for jk, which is symmetric in j and k; hence, the sum above for dω simplifies to dω=0. Thus, ω is closed, and we see that the flow (14) is gradient for any order parameter function of the form (13).

Next, we show that the hyperbolic potential for V, up to an additive constant, is given by

Φ(w)=i=1Nailog1|w|2|wpi|2=1d1i=1NailogPhyp(w,pi).
(15)

Here, we follow the convention that the potential decreases along trajectories; therefore, we are asserting that hypΦ=V. To derive this, we use the identity euc|ww0|2=2(ww0) for any constant vector w0Rd. Then,

eucΦ(w)=i=1Nai(2w1|w|22(wpi)|wpi|2)=21|w|2i=1Nai((1|w|2)(piw)|wpi|2w)=21|w|2i=1NaiMw(pi)=21|w|2Z(Mw(p)).

Hence, we see that

hypΦ(w)=12(1|w|2)Z(Mw(p))=V(w),

as desired.

We can use the existence of the potential Φ(w) for the flow on Bd to prove a global synchrony result for the system (1) when the coefficients ai in the order parameter Z are all positive. Specifically, we assume that 0<ai<1/2 for all i, and i=1Nai=1. We also assume N3, and all the rotation terms Ai in (1) are equal. Under these conditions, almost all trajectories for (1) converge in forward time to the (d1)-dimensional diagonal manifold ΔX as t, meaning that the system self-synchronizes. In contrast, in backward time, the system tends to an incoherent state having zero order parameter: as t, almost all trajectories for (1) converge to the codimension-d subspace ΣX consisting of states with Z(p)=0.

The proof is modeled after Theorem 1 in Chen et al.27 and will be based on two preliminary lemmas. In each of these lemmas, we assume the conditions on the ai above and that the base point p=(pi) for the flow (14) has all distinct coordinates.

We begin with a general observation about gradient flows in the ball Bd: if w0Bd is any initial condition and wBd is in the forward limit set Ω+(w0), then w is a fixed point for the flow. To see this, let Φ be a potential for the flow, and suppose w(tn)wBd for some sequence tn. Since the potential decreases along trajectories,

limtΦ(w(t))=limnΦ(w(tn))=Φ(w).

Let Ft denote the time-t flow map. If w is not a fixed point, then for any s>0,

limtΦ(w(t))=limnΦ(w(tn+s))=limnΦ(Fs(w(tn)))=Φ(Fs(w))<Φ(w),

which is a contradiction; therefore, w must be a fixed point. (Compact limit sets are connected; therefore, Ω+(w0) cannot consist of two or more but finitely many fixed points; however, it is possible that forward or backward limit sets for gradient flows consist of a continuum of fixed points. We will see that this is not the case for our system on Bd.

Lemma 1

Any fixed point for the flow (14) in Bd is repelling.

Proof.
Suppose wBd is a fixed point for (14). As discussed above, an advantage of using the w-parameter is the equivariance with respect to change of base point p. Consequently, we can assume w=0 without loss of generality; therefore, Z(p)=i=1Naipi=0. To first order in w,
Mw(pi)=piw12w,piw=(piw)(1+2w,pi)w=pi2w+2w,pipi.
The linearization of (14) at the fixed point w=0 is
w˙=12i=1Nai(pi2w+2w,pipi)=wi=1Naiw,pipi.
We claim that the linear map
Tw=i=1Naiw,pipi
has ||T||<1; to see this, suppose |w|=1. Then, |w,pipi|1, and Tw is a convex combination of the vectors w,pipi. We can only obtain |Tw|=1 if all terms w,pipi=u with |u|=1, which implies all pi=±u, and this cannot happen if at least three of the pi are distinct. Hence, ||T||<1, and therefore, the eigenvalues μi of T satisfy |μi|<1. The eigenvalues for the w˙ linearization are λi=1μi; therefore, we see that Reλi>0 for all i, establishing that the fixed point w is repelling.
Lemma 2

lim|w|1Φ(w)=.

Proof.
It suffices to show that
limnΦ(wn)=
for any sequence wnBd with wnxSd1. The result is clear if xpi: as n, the terms |wnpi| in the potential (15) are bounded away from 0 and 1|wn|20. Therefore, let us say that wnp1. We rewrite Φ(wn) as
Φ(wn)=log(1|wn|2)2a1log|wnp1|2i=2Nailog|wnpi|=log(1|wn|)2a1log|wnp1|+log(1+|wn|)2i=2Nailog|wnpi|.
The latter two terms above have finite limit as n; therefore, we focus on the first two terms. We have 1|wn||wnp1|; therefore,
log(1|wn|)2a1log|wnp1|(12a1)log|wnp1|
as n, which proves our result. Notice that we need the assumption ai<1/2 for this argument.
Theorem

Under the conditions above, almost all trajectories for (1) converge to Δ as t and to Σ as t.

Proof.

Let p=(p1,,pN)X be any point with all distinct coordinates. The points on Gp are parameterized by wBd and ζSO(d), and the dynamics for these parameters are given by (5). We begin with the dynamics as t. Let w(t) be a trajectory for (14) with initial condition w0Bd and consider the backward-time limit set Ω(w0); this is a nonempty, compact, connected subset of Bd¯. The potential Φ is decreasing along all trajectories w(t), hence bounded below as t; therefore, Lemma 2 implies that the limit set Ω(w0) must be contained in the interior Bd. We know that any wΩ(w0) is a fixed point for the flow. By Lemma 1, w is repelling, and therefore, any trajectory w(t) that comes sufficiently close to w must have w(t)w as t; therefore, Ω(w0)={w}. This proves the existence of fixed points for (14) and that every trajectory w(t) converges to a fixed point as t. If the flow had multiple fixed points, we would obtain a partition of Bd into the disjoint open basins of repulsion of the fixed points, violating connectedness of the ball. Therefore, (14) has a unique fixed point w, and w(t)w as t for all trajectories. The fixed point w has Z(Mw(p))=0; therefore, all trajectories in Gp converge to Σ as t.

In forward time, the limit set Ω+(w0) for any w0w must be completely contained in the boundary Sd1 since the unique fixed point wBd is repelling. Suppose we remove the factor (1/2)(1|w|2) in the flow (14); the scaled vector field on Bd given by
w˙=i=1NaiMw(pi)=wi=1Nai((1|w|2)(piw)|piw|2)
(16)
has the same trajectories as the original flow, just with different time parameterizations. Observe that this scaled vector field extends smoothly to Rd{pi} and coincides with the radial vector field x at any xSd1 with xpi. Therefore, there is a unique trajectory passing through each point xSd1, flowing from the interior to the exterior of the sphere, as long as xpi. Consequently, the original flow (14) has a unique trajectory w(t) in Bd with w(t)x as t, as long as xpi. This also shows that there is a neighborhood U of Sd1{pi} such that if w(t0)U for some t0, then w(t)xpi for some xSd1. Therefore, if Ω+(w0) contains some xpi, then the trajectory w(t) of w0 must enter the neighborhood U, and therefore, w(t)xSd1 as t.
Since limit sets are connected, the only other possibility is Ω+(w0)={pi} for some i; equivalently, w(t)pi. We will show that there is a unique trajectory with this behavior for each pi. Assuming this, we see that with N+1 exceptions, any trajectory w(t) converges to a point xSd1 with xpi (the exceptions are the N trajectories converging to the base point coordinates pi and the fixed point trajectory w). The corresponding trajectory in Gp has coordinates
ζ(t)Mw(t)(pi)=ζ(t)((1|w(t)|2)(piw(t))|piw(t)|2w(t)).
We have |w(t)|1 and |piw(t)| is bounded away from 0 as t; as a result, Mw(t)(pi)x for each i; therefore, the trajectory ζ(t)Mw(t)(p) in Gp converges to Δ as t.
This analysis breaks down at x=pi because the scaled vector field above does not have a unique limit as wpi; alternatively, its limit depends on the direction of the approach. To see this, write w=p1ru, where 0<r<1 and |u|=1 (with this convention, u=p1 corresponds to w approaching p1 radially). Then, |p1w|=r and
|w|2=12rp1,u+r2;
therefore,
(1|w|2)(p1w)|p1w|2=(2rp1,ur2)rur2=(2p1,ur)u.
As r0, the magnitude of this term is 2p1,u, which depends on the angle of approach given by u (note that p1,u>0 because u points outward at p1).

To complete the proof, we will examine the scaled system (16) using the polar representation (r,u) and show that the polar system has the unique fixed point r=0,u=p1, which has a unique attracting trajectory because it is a saddle with a (d1)-dimensional unstable manifold.

We see that the scaled system has
w˙=p1rua1(2p1,ur)u+O(r)=p12a1p1,uu+O(r),
where the O(r) term is a smooth function of r and u for |r|<ε=min|pip1|, i2. This condition ensures that |piw||pip1||r|>0; therefore, the i2 terms in the scaled w˙ equation are all smooth functions of r and u. Also, we can allow r<0 here, even though it is not relevant to the w˙ system. Now, r2=|wp1|2; therefore,
rr˙=wp1,w˙=ru,w˙,
which gives
r˙=(12a1)p1,u+O(r).
Differentiating ru=p1w gives
ru˙=r˙uw˙=(12a1)p1,uu(p12a1p1,uu)+O(r)=p1,uup1+O(r).
Hence, the scaled system in polar form can be written as
rr˙=(12a1)rp1,u+O(r2),ru˙=p1,uup1+O(r).
We emphasize that the O(r) and O(r2) terms are smooth functions of r,u as long as |r|<ε. We consider the “semi-scaled” polar system
r˙=(12a1)rp1,u+O(r2),
(17a)
u˙=p1,uup1+O(r),
(17b)
which has the same trajectories as the original system, just with different time parameterizations. The advantage of this modified system is that the equations are smooth on (ε,ε)×Sd1.
Observe that the system (17) has {0}×Sd1 invariant and has a fixed point (r,u)=(0,p1). The fixed point p1 is repelling on the invariant manifold {0}×Sd1; to see this, observe that
p1,u˙=p1,u21
when r=0. In fact, if we assign the coordinate θ on any great circle joining p1 and p1 on Sd1 so that u=eiθ and p1=1, then the system reduces to θ˙=sinθ. We also see that the r˙ equation linearized at (0,p1) is r˙=(12a1)r; therefore, the linearization of (17) has the single negative eigenvalue (12a1) and d1 positive eigenvalues +1. Therefore, (0,p1) is a saddle with a one-dimensional stable manifold and hence has a unique trajectory (r(t),u(t))(0,p1) with r(t)>0.
Now, suppose we have a trajectory w(t)p1 in our original system (14). The corresponding trajectory for (17) will have r(t)0; we cannot achieve r(t)=0 in finite time because the manifold {r=0} is invariant for (17). We must prove that u(t)p1 so that this trajectory is, in fact, the saddle stable manifold. Observe that
p1,u˙=p1,u21+O(r).
Also, note that p1,u(t)>0 since |w(t)|<1. Let 0<c<1; then, for some T0, tT implies O(r(t))(1c2)/2. Now, suppose 0<p1,u(t0)<c for some t0T; then, 0<p1,u(t)c for all tt0. This is because the function tp1,u(t) is decreasing if 0<p1,u(t)<c,
p1,u(t)˙c21+12(1c2)=12(1c2)
as long as 0<p1,u(t)<c. However, this also implies that eventually p1,u(t)<0, which is a contradiction. Hence, we must have p1,u(t)c for all tT, which proves that u(t)p1.

The natural hyperbolic geometry on the unit ball, with isometries consisting of the higher-dimensional Möbius group, is key to understanding the dynamics of the Kuramoto model on a sphere. Using this framework, we see that dynamical trajectories of (1) are constrained to lie on group orbits, and we can explicitly give the equations for the reduced dynamics on the group orbits. For the special class of linear order parameters, the dynamics can be further reduced to a flow on the unit ball Bd, which is gradient with respect to the hyperbolic metric. We analyze this flow and prove a global synchronization result for the system (1) for linear order parameters with positive weights and no weight greater than half the total. This illustrates the power of the geometric/group-theoretic approach.

We conclude with some directions for future research. The case of linear order parameters with both positive and negative weights can, in principle, be explored using similar methods; it will also reduce to a hyperbolic gradient system on the ball Bd. In particular, the case when the sum of the weights is 0 should have some intriguing dynamics. For the original (d=2) Kuramoto model, the dynamics are Hamiltonian and equivalent to the vector field on the unit disk given by placing a collection of point charges on the unit circle, with total charge 0, as shown by Chen et al.27 Of course, this result cannot generalize to all higher dimensions d; Hamiltonian dynamics is only possible in even dimensions.

Another possible direction to explore is the case of systems where the oscillator population is divided into two or more families with different intrinsic rotational terms Ai, and the coupling across families differs from the coupling within families. For the case d=2, these systems often support “chimera states,” in which one or more families synchronize while others tend to a partially disordered configuration. These states can be dynamically stable within their Möbius group orbits. Our framework enables the dynamical reduction of multi-family networks of higher-dimensional oscillators and makes possible the study of the dynamics of these networks without necessarily passing to the continuum limit, as is often done as a simplifying step in the analysis of Kuramoto networks.

This research was supported in part by the National Science Foundation (NSF) Research Training Group Grant: Dynamics, Probability, and PDEs in Pure and Applied Mathematics (Grant No. DMS-1645643) and by the NSF (Grant No. DMS-1910303). We thank Vladimir Jaćimović and Max Lohe for helpful comments on a preprint of this paper.

The data that support the findings of this study are available within the article.

1.
Y.
Kuramoto
, “Self-entrainment of a population of coupled non-linear oscillators,” in International Symposium on Mathematical Problems in Theoretical Physics (Springer, 1975), pp. 420–422.
2.
Y.
Kuramoto
,
Chemical Oscillations, Waves, and Turbulence
(
Springer
,
1984
).
3.
A. T.
Winfree
, “
Biological rhythms and the behavior of populations of coupled oscillators
,”
J. Theor. Biol.
16
,
15
42
(
1967
).
4.
A. T.
Winfree
,
The Geometry of Biological Time
(
Springer
,
1980
).
5.
S. H.
Strogatz
, “
From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators
,”
Phys. D
143
,
1
20
(
2000
).
6.
A.
Pikovsky
,
M.
Rosenblum
, and
J.
Kurths
,
Synchronization: A Universal Concept in Nonlinear Sciences
(
Cambridge University Press
,
2003
), Vol. 12.
7.
S.
Strogatz
,
Sync
(
Hyperion
,
2003
).
8.
J. A.
Acebrón
,
L. L.
Bonilla
,
C. J. P.
Vicente
,
F.
Ritort
, and
R.
Spigler
, “
The Kuramoto model: A simple paradigm for synchronization phenomena
,”
Rev. Mod. Phys.
77
,
137
(
2005
).
9.
F.
Dörfler
and
F.
Bullo
, “
Synchronization in complex networks of phase oscillators: A survey
,”
Automatica
50
,
1539
1564
(
2014
).
10.
A.
Pikovsky
and
M.
Rosenblum
, “
Dynamics of globally coupled oscillators: Progress and perspectives
,”
Chaos
25
,
097616
(
2015
).
11.
F. A.
Rodrigues
,
T. K. D.
Peron
,
P.
Ji
, and
J.
Kurths
, “
The Kuramoto model in complex networks
,”
Phys. Rep.
610
,
1
98
(
2016
).
12.
C.
Bick
,
M.
Goodfellow
,
C. R.
Laing
, and
E. A.
Martens
, “
Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: A review
,”
J. Math. Neurosci.
10
,
9
(
2020
).
13.
K.
Wiesenfeld
,
P.
Colet
, and
S. H.
Strogatz
, “
Synchronization transitions in a disordered Josephson series array
,”
Phys. Rev. Lett.
76
,
404
(
1996
).
14.
K.
Wiesenfeld
,
P.
Colet
, and
S. H.
Strogatz
, “
Frequency locking in Josephson arrays: Connection with the Kuramoto model
,”
Phys. Rev. E
57
,
1563
(
1998
).
15.
K. Y.
Tsang
,
R. E.
Mirollo
,
S. H.
Strogatz
, and
K.
Wiesenfeld
, “
Dynamics of a globally coupled oscillator array
,”
Phys. D
48
,
102
112
(
1991
).
16.
J. W.
Swift
,
S. H.
Strogatz
, and
K.
Wiesenfeld
, “
Averaging of globally coupled oscillators
,”
Phys. D
55
,
239
250
(
1992
).
17.
S.
Nichols
and
K.
Wiesenfeld
, “
Ubiquitous neutral stability of splay-phase states
,”
Phys. Rev. A
45
,
8430
(
1992
).
18.
S.
Watanabe
and
S. H.
Strogatz
, “
Integrability of a globally coupled oscillator array
,”
Phys. Rev. Lett.
70
,
2391
(
1993
).
19.
S.
Watanabe
and
S. H.
Strogatz
, “
Constants of motion for superconducting Josephson arrays
,”
Phys. D
74
,
197
253
(
1994
).
20.
C. J.
Goebel
, “
Comment on ‘Constants of motion for superconductor arrays
,’”
Phys. D
80
,
18
20
(
1995
).
21.
E.
Ott
and
T. M.
Antonsen
, “
Low dimensional behavior of large systems of globally coupled oscillators
,”
Chaos
18
,
037113
(
2008
).
22.
E.
Ott
and
T. M.
Antonsen
, “
Long time evolution of phase oscillator systems
,”
Chaos
19
,
023117
(
2009
).
23.
A.
Pikovsky
and
M.
Rosenblum
, “
Partially integrable dynamics of hierarchical populations of coupled oscillators
,”
Phys. Rev. Lett.
101
,
264103
(
2008
).
24.
S. A.
Marvel
,
R. E.
Mirollo
, and
S. H.
Strogatz
, “
Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action
,”
Chaos
19
,
043104
(
2009
).
25.
I.
Stewart
, “
Phase oscillators with sinusoidal coupling interpreted in terms of projective geometry
,”
Int. J. Bifurcation Chaos
21
,
1795
1804
(
2011
).
26.
B.
Chen
,
J. R.
Engelbrecht
, and
R.
Mirollo
, “
Hyperbolic geometry of Kuramoto oscillator networks
,”
J. Phys. A: Math. Theor.
50
,
355101
(
2017
).
27.
B.
Chen
,
J. R.
Engelbrecht
, and
R.
Mirollo
, “
Dynamics of the Kuramoto-Sakaguchi oscillator network with asymmetric order parameter
,”
Chaos
29
,
013126
(
2019
).
28.
M.
Lohe
, “
Non-Abelian Kuramoto models and synchronization
,”
J. Phys. A: Math. Theor.
42
,
395101
(
2009
).
29.
M.
Lohe
, “
Quantum synchronization over quantum networks
,”
J. Phys. A: Math. Theor.
43
,
465301
(
2010
).
30.
T.
Tanaka
, “
Solvable model of the collective motion of heterogeneous particles interacting on a sphere
,”
New J. Phys.
16
,
023016
(
2014
).
31.
D.
Chi
,
S.-H.
Choi
, and
S.-Y.
Ha
, “
Emergent behaviors of a holonomic particle system on a sphere
,”
J. Math. Phys.
55
,
052703
(
2014
).
32.
S.-H.
Choi
and
S.-Y.
Ha
, “
Complete entrainment of Lohe oscillators under attractive and repulsive couplings
,”
SIAM J. Appl. Dyn. Syst.
13
,
1417
1441
(
2014
).
33.
S.-Y.
Ha
,
D.
Ko
,
J.
Park
, and
X.
Zhang
, “
Collective synchronization of classical and quantum oscillators
,”
EMS Surv. Math. Sci.
3
,
209
267
(
2016
).
34.
S.-Y.
Ha
,
D.
Ko
, and
S. W.
Ryoo
, “
On the relaxation dynamics of Lohe oscillators on some Riemannian manifolds
,”
J. Stat. Phys.
172
,
1427
1478
(
2018
).
35.
M.
Lohe
, “
Higher-dimensional generalizations of the Watanabe–Strogatz transform for vector models of synchronization
,”
J. Phys. A: Math. Theor.
51
,
225101
(
2018
).
36.
V.
Jaćimović
and
A.
Crnkić
, “
Low-dimensional dynamics in non-Abelian Kuramoto model on the 3-sphere
,”
Chaos
28
,
083105
(
2018
).
37.
S.
Chandra
,
M.
Girvan
, and
E.
Ott
, “
Continuous versus discontinuous transitions in the D-dimensional generalized Kuramoto model: Odd D is different
,”
Phys. Rev. X
9
,
011002
(
2019
).
38.
S.
Chandra
,
M.
Girvan
, and
E.
Ott
, “
Complexity reduction ansatz for systems of interacting orientable agents: Beyond the Kuramoto model
,”
Chaos
29
,
053107
(
2019
).
39.
M.
Lohe
, “
Systems of matrix Riccati equations, linear fractional transformations, partial integrability and synchronization
,”
J. Math. Phys.
60
,
072701
(
2019
).
40.
L.
DeVille
, “
Synchronization and stability for quantum Kuramoto
,”
J. Stat. Phys.
174
,
160
187
(
2019
).
41.
S.-Y.
Ha
and
H.
Park
, “
From the Lohe tensor model to the Lohe Hermitian sphere model and emergent dynamics
,”
SIAM J. Appl. Dyn. Syst.
19
,
1312
1342
(
2020
).
42.
M.
Lohe
, “
On the double sphere model of synchronization
,”
Phys. D
412
,
132642
(
2020
).
43.
S.-Y.
Ha
,
D.
Kim
,
H.
Park
, and
S. W.
Ryoo
, “
Constants of motion for the finite-dimensional Lohe type models with frustration and applications to emergent dynamics
,”
Phys. D
416
,
132781
(
2021
).
44.
V.
Jaćimović
and
A.
Crnkić
, “
On reversibility of macroscopic and microscopic dynamics in the Kuramoto model
,”
Phys. D
415
,
132762
(
2021
).
45.
X.
Dai
,
K.
Kovalenko
,
M.
Molodyk
,
Z.
Wang
,
X.
Li
,
D.
Musatov
,
A.
Raigorodskii
,
K.
Alfaro-Bittner
,
G.
Cooper
,
G.
Bianconi
et al., “
D-dimensional oscillators in simplicial structures: Odd and even dimensions display different synchronization scenarios
,”
Chaos, Solitons Fractals
146
,
110888
(
2021
).
46.
R.
Olfati-Saber
, “Swarms on sphere: A programmable swarm with synchronous behaviors like oscillator networks,” in Proceedings of the 45th IEEE Conference on Decision and Control (IEEE, 2006), pp. 5060–5066.
47.
A.
Beardon
,
The Geometry of Discrete Groups
(
Springer
,
1983
).
48.
M.
Stoll
, “Harmonic function theory on real hyperbolic space”; see https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.561.4447&rep=rep1&type=pdf.
49.
W.
Rudin
, Function Theory in the Unit Ball of Cn (
Springer
,
1980
).