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 $N\u22653$ and to clarify why it admits the analog of the Ott–Antonsen ansatz in the continuum limit $N\u2192\u221e$. 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 $(d\u22121)$-dimensional sphere in $d$-dimensional Euclidean space, yielding a state space of dimension $N(d\u22121)$. 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 $N\u22653$ 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.

## I. INTRODUCTION

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 model^{13,14} and themselves displayed remarkable dynamical features, such as surprisingly low-dimensional invariant tori^{15,16} and ubiquitous neutral stability of splay states^{17} 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 $N\u22123$ constants of motion for all $N\u22653$. Goebel^{20} 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 $\rho (\theta ,\omega ,t)$ of oscillators having phase $\theta $ and intrinsic frequency $\omega $ 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 $N\u22653$ 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.

## II. PRELIMINARIES

### A. The Kuramoto model on a sphere

In a pioneering paper, Lohe^{28} observed that there are at least two natural generalizations of the Kuramoto model to higher dimensions. One of them replaces the phases $\theta j$ of the original Kuramoto model^{1,2} with complex numbers $exp\u2061(i\theta j)$ on the unit circle and then views those as equivalent to $2\xd72$ rotation matrices parameterized by a rotation angle $\theta 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 $Sd\u22121$ 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

where $xi$ is a point on the unit sphere $Sd\u22121\u2282Rd$, each $Ai$ is an antisymmetric $d\xd7d$ matrix, and $Z\u2208Rd$ 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,\u2026,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 $\omega i$ in the original Kuramoto model.

A straightforward computation shows that the dot product between an oscillator’s instantaneous position and instantaneous velocity satisfies $\u27e8xi,xi\u02d9\u27e9=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=(Sd\u22121)N$, which has dimension $N(d\u22121)$. Later, we will also consider the natural infinite-$N$ analog of (1), where a state is a probability measure on $Sd\u22121$.

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

where the $ai$ are real constants.

### B. General philosophy: Lie groups and reducible systems

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 $\gamma (t)$ be a smooth curve in $G$ with $\gamma (0)=e$, the identity element in $G$. Then, the derivative $\gamma \u02d9(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 $x\u2208X$, $t\u21a6\gamma (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 $x\u2208X$, 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 $\xi $ on $X$, which defines a dynamical system on the state space $X$. If $\xi x\u2208Vx$ at each point $x\u2208X$, then the flow corresponding to the vector field $\xi $ 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 $dim\u2061X$ to $dim\u2061G$. Now, suppose, as is the case in applications of this methodology, the correspondence $G\u2192Gx$ is one-to-one for generic $x\u2208X$; equivalently, the stabilizer subgroup $Gx={e}$ for generic $x\u2208X$. Then, for each $x\u2208X$, 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 $\nu 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)$, $xi\u2208Rd$, and then the infinitesimal generators have the form $(\nu x)i=Axi$ for some skew-symmetric matrix $A$. We could also, if we like, restrict $X$ to $N$-tuples $(xi)$ with $xi\u2208Sd\u22121$, the state space for (1). Now, suppose we had a dynamical system on $X$ of the form $x\u02d9i=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(d\u22121)/2$, and for large $N$, this is much smaller than the dimension of $X$, which is $N(d\u22121)$. Basically, the configuration $(xi)$ of points on the sphere $Sd\u22121$ 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 $Z\u22600$, we need a different group action to make this strategy work. Fortunately, vector fields of the form seen in the Kuramoto model,

turn out to arise as the infinitesimal generators of the group action of a larger group $G$ acting on the sphere $Sd\u22121$ 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=\u2211i=1Naixi$ with $ai>0$. However, first, we need to show that vector fields on $Sd\u22121$ of the form in (2) are indeed the infinitesimal generators of the action of the larger Möbius group.

### C. Hyperbolic geometry and Möbius transformations

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 $Sd\u22121$. This geometry has metric

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 $\u22121$, and we can describe its isometries, which generalize the Möbius transformations preserving the unit disk for $d=2$. For $d=2$, let $w\u2208B2$ and consider the Möbius transformation

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 $2\u27e8w,x\u27e9=w\xafx+wx\xaf$, we see that

This form of $Mw$ generalizes to higher dimensions: Let $w\u2208Bd$ and define

where $x\u2208Bd$ or $Sd\u22121$. We call $Mw$ a *boost transformation*.

If $|x|=1$, this formula simplifies to

Now, see that

Alternatively, we can use the second form to show

Similar computations show that $M0$ is the identity, $Mw\u22121=M\u2212w,$ and $Mw(0)=\u2212w$.

It is a standard result in hyperbolic geometry (e.g., see Theorem 3.5.1 in Beardon^{47}) 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

and also uniquely in the form

for some vectors $w,z\u2208Bd$ and rotations $\zeta ,\xi \u2208SO(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(d\u22121)/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,\zeta $ and $z,\xi $ are related. We can do this by comparing the linearizations of the two formulas for $g(x)$ at $x=0$,

Equating coefficients implies $z=\u2212\zeta w$ (hence, $|z|=|w|$) and $\xi =\zeta $.

#### 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

where $A$ is an antisymmetric $d\xd7d$ matrix and $Z\u2208Rd$ is a vector. Note that, as advertised, this flow extends to a flow on $Sd\u22121$ 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$,

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=\u22122w$ and $A=0$. Next, recall that the infinitesimal generators corresponding to the rotation components are flows of the form $x\u02d9=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 $p\u2208X$ under the system (1) with all $Ai=A$ lies in the group orbit $Gp$.

## III. REDUCED EQUATIONS

The given Kuramoto system has $N(d\u22121)$ 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)\u2245Bd$. 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,\u2026,pN)\u2208X$. 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 $g\u2208G$, with parameters $w,z$, and $\zeta $. We wish to derive the corresponding evolution equations for $w,z$, and $\zeta $. 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,\u2026,N$, we have $xi(t)=gt(pi)$ for a unique $gt\u2208G$, which determines the parameters $w,z,\zeta $ 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\xaf$, and its time-$t$ flow must be given by some $g~t\u2208G$. This ODE has solutions $xi(t)=gt(pi)=gt(g0\u22121(xi(0)))$, which implies that $g~t=gtg0\u22121$.

Therefore, for any $y0\u2208Bd\xaf$,

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)=\u2212\zeta w=z$; therefore, $z$ satisfies the ODE (3).

Now, expand $y=\zeta Mw(y0)=M\u2212z(\zeta y0)$ to first order in $y0$, using the variables $z$ and $\zeta $,

therefore,

On the other hand, (3) gives

Setting $y0=0$ gives the $z\u02d9$ equation

as expected, and since $\u27e8Az,z\u27e9=0$, this in turn implies that

Equating the $y0$ terms, factoring out $1\u2212|z|2$, and canceling the common term $\u27e8Z,z\u27e9\zeta y0$ gives

Together, the last two terms above define a special type of antisymmetric operator of $\zeta y0$: Given any $y1,y2\u2208Rd$, define the antisymmetric operator $\alpha $ as

this operator has range = $span(y1,y2)$ provided that $y1$ and $y2$ are linearly independent; otherwise, $\alpha (y1,y2)=0$. Then, for all $y0\u2208Rd$,

therefore,

Differentiating $z=\u2212\zeta w$ gives

therefore,

hence,

Summing up, the evolution equations for the $(z,\zeta )$ coordinate system on $Gp$ are

with $A$ and $Z$ evaluated at $M\u2212z(\zeta p)$ and for the $(w,\zeta )$ coordinate system on $Gp$ are

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

## IV. COMPARISON OF Z VS W COORDINATES

The $z\u02d9$ equation (4) is an extension of the system equation (1) on $Sd\u22121$. However, for finite $N$, the $z\u02d9$ equation does not uncouple from $\zeta $ since $Z$ is evaluated at $M\u2212z(\zeta p)$. The exception to this is in the infinite-$N$ limit: if the base point $p$ is now the uniform density on $Sd\u22121$, then $\zeta p=p$ (the uniform density is invariant under rotations) and the density $M\u2212z(p)$ is a hyperbolic Poisson density on $Sd\u22121$ 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 $d\u22653$ (we will give more details on this in Sec. V).

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

with $ai\u2208R$, $\zeta $ drops out of the $w\u02d9$ equation, and we get the reduced equation

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\u2032=M(p)\u2208Gp$, then, we have coordinates $w\u2032,\zeta \u2032$ associated with the base point $p\u2032$. Any $q\u2208Gp$ has two expressions

Assuming the coordinates of $p$ are in sufficiently general position, this implies $\zeta Mw=\zeta \u2032Mw\u2032\xb0M$, and hence,

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

## V. CONTINUUM LIMIT

Next, we consider the dynamics of the Kuramoto model (1) in the limit $N\u2192\u221e$. 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,

In the continuum limit, a state of the system is a probability measure $\rho $ on $Sd\u22121$, and the order parameter becomes

The measure $\rho $ 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 $M\u2208G$, then the measure $M\u2217\rho $ is defined by the adjunction formula

In particular, we can consider the $G$-orbit of the uniform probability measure $\sigma $ on $Sd\u22121$. This orbit is special; whereas a typical group orbit $G\rho $ has dimension equal to the dimension of $G$, namely, $d(d+1)/2$, the orbit $G\sigma $ has dimension only $d$. This is because the stabilizer of $\sigma $ is $SO(d)$; any rotation fixes $\sigma $, whereas the boosts deform $\sigma $. Hence, the orbit $G\sigma $ has dimension $d$. Any element in $G\sigma $ can be written as $(M\u2212z)\u2217\sigma $, with $z\u2208Bd$. The evolution equation for $z$ is (4), with

In the case $d=2$ with $x=\zeta \u2208S1$, we have

therefore, the integral

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

when $d=2$. Unfortunately, the formula $Z(z)=Kz$ is not correct for $d\u22653$; 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 $\Delta $ associated with its metric; functions $f$ on $X$ satisfying the equation $\Delta f=0$ are called *harmonic*. For functions on the ball $Bd$ with the hyperbolic metric, this operator is

where

is the standard Laplace operator (see Stoll,^{48} Chap. 3). We will call solutions to the equation $\Delta hypf=0$*hyperbolic 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 $Sd\u22121$, extend $f$ to a hyperbolic harmonic function $f~$ on $Bd$. Assuming that this problem has a unique solution, then for any rotation $\zeta \u2208SO(d)$, we must have $f\xb0\zeta ~=f~\xb0\zeta $ since rotations preserve the hyperbolic metric. If we average $f\xb0\zeta $ on $Sd\u22121$ over all rotations $\zeta \u2208SO(d)$, we get the constant function

on $Sd\u22121$, and any constant is hyperbolic harmonic on $Bd$. Therefore, the average on $Bd$ of $f\xb0\zeta ~=f~\xb0\zeta $ over all $\zeta \u2208SO(d)$ must be the constant $fave$. However, $f~(\zeta (0))=f~(0)$ for all $\zeta $; therefore, we must have

Now, let $z\u2208Bd$; since $M\u2212z$ preserves the hyperbolic metric, we must have $f\xb0M~\u2212z=f~\xb0M\u2212z$, which implies

As shown in Chap. 5 of Stoll,^{48} the measure $(M\u2212z)\u2217\sigma $ is given by the formula

with hyperbolic Poisson kernel function

Thus, the solution to the hyperbolic Dirichlet problem with boundary function $f$ on $Sd\u22121$ is given by the hyperbolic Poisson integral

The orbit $G\sigma $ consists of all *hyperbolic* Poisson measures $P(z,x)d\sigma (x)$, parameterized by $z\u2208Bd$. By contrast, the Euclidean Poisson kernel function is

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 $d\u22652$. We see from (8) that $Z(z)$ is the hyperbolic Poisson integral of the function $Kx$ on $Sd\u22121$. 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 $Sd\u22121$ to a hyperbolic harmonic function on $Bd$ is given by

where $F$ is the hypergeometric function

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

## VI. COMPLEX CASE

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

where now $xi$ is a point on the unit sphere $S2m\u22121\u2282Cm$, $Ai$ is an anti-Hermitian $m\xd7m$ complex matrix, $Z\u2208Cm$, and $\u27e8,\u27e9$ 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 $m\u22652$. To see this, suppose

for all $x\u2208S2m\u22121\u2282Cm=Rd$, where $A$ is antisymmetric, $B$ is anti-Hermitian, $Y,Z\u2208Cm$, and we use the subscripts $R$ and $C$ to distinguish the real and complex inner products. Then,

and therefore, $(A\u2212B)(\u2212x)=(A\u2212B)x$ for all $x\u2208S2m\u22121$, which implies $A=B$. This implies

for all $x\u2208S2m\u22121$; hence, $Y\u2212Z\u2208spanC(x)$ for all $x\u2208Cm$; if $m\u22652$, this implies $Y=Z$. However, then we have

for all $x\u2208Cm$, which can only hold if $Y=0$. Hence, for $m\u22652$, 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 $\zeta \u2208U(m)$ and boost transformations of the form

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

where $w,z\u2208Bm$, but now $\zeta ,\xi \u2208U(m)$, the complex unitary group. Linearizing at $x=0$ gives

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

with $A$ being anti-Hermitian $m\xd7m$ and $Z\u2208Cm$. This flow extends to a flow on $S2m\u22121$ 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

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

## VII. RELATION TO PREVIOUS RESEARCH

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.

Tanaka^{30} 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*. found^{24} 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\u02d9$ 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\u02d9$ 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 reduction^{21} for the complex version of the system.

Lohe^{35} 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 $Sd\u22121$ is our $Mw$ (with $v=w$) and his Eq. (31) is the same as our $z\u02d9$ equation (4). He also has something that looks like the $w\u02d9$ 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

where $Qi\u2208O(d)$ and $\lambda i\u2208R$. However, such a $Z$ does not satisfy the identity $\zeta Z(p)=Z(\zeta p)$ for all rotations $\zeta $ unless $Qi=\xb1I$; therefore, we do not see how the $\zeta $ term cancels.

Lohe’s map $M$ in his Eq. (55) (ignoring the $R$ factor) agrees with our map $M\u2212v$ on the sphere $Sd\u22121$, but not on the ball $Bd$. Therefore, it is not a Möbius transformation of the type we are using. For example, $M(\u2212v)=v$, whereas $M\u2212v(\u2212v)=0$. We are not sure why Lohe^{35} 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 $M\u2212v$.

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 $Sd\u22121$, 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\u02d9$ equation (4) in the infinite-$N$ limit. The integral in their Eq. (19) can be evaluated, as shown above in (10).

## VIII. AN EXAMPLE: A FIRST-ORDER LINEAR ORDER PARAMETER GIVES A GRADIENT SYSTEM

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

where the $ai$ are real constants.

### A. Computer visualization

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.

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 $\u27e8xi(0),Z(0)\u27e9>0$ for all $i$ will synchronize. Geometrically, this condition is equivalent to requiring that the oscillators all lie on the hemisphere given by $\u27e8xi,y\u27e9>0$ for some vector $y\u22600$. More generally, similar partial synchronization results are obtained by Ha *et al.*^{43} for the case

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.

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.

### B. Existence of a hyperbolic gradient

independent of the parameter $\zeta $. 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,\u2026,wd$, the $1$-form associated with the vector field with components $f1,\u2026,fd$ is

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

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

where $\u2207euc$ denotes the ordinary Euclidean gradient operator. We have $\varphi (w)=2(1\u2212|w|2)\u22121$ for the hyperbolic metric; therefore, the hyperbolic gradient operator on $Bd$ is given by

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

where $p1,j$ denotes the $j$th component of the point $p1\u2208Sd\u22121$. Let $Ej$ denote the coefficient of $dwj$ in parentheses above; then,

Applying the chain and quotient rules gives

for $j\u2260k$, which is symmetric in $j$ and $k$; hence, the sum above for $d\omega $ simplifies to $d\omega =0$. Thus, $\omega $ 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

Here, we follow the convention that the potential *decreases* along trajectories; therefore, we are asserting that $\u2207hyp\Phi =\u2212V$. To derive this, we use the identity $\u2207euc|w\u2212w0|2=2(w\u2212w0)$ for any constant vector $w0\u2208Rd$. Then,

Hence, we see that

as desired.

### C. Analysis of dynamics

We can use the existence of the potential $\Phi (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 $\u2211i=1Nai=1$. We also assume $N\u22653$, and all the rotation terms $Ai$ in (1) are equal. Under these conditions, almost all trajectories for (1) converge in forward time to the $(d\u22121)$-dimensional diagonal manifold $\Delta \u2282X$ as $t\u2192\u221e$, meaning that the system self-synchronizes. In contrast, in backward time, the system tends to an incoherent state having zero order parameter: as $t\u2192\u2212\u221e$, almost all trajectories for (1) converge to the codimension-$d$ subspace $\Sigma \u2282X$ 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 $w0\u2208Bd$ is any initial condition and $w\u2217\u2208Bd$ is in the forward limit set $\Omega +(w0)$, then $w\u2217$ is a fixed point for the flow. To see this, let $\Phi $ be a potential for the flow, and suppose $w(tn)\u2192w\u2217\u2208Bd$ for some sequence $tn\u2192\u221e$. Since the potential decreases along trajectories,

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

which is a contradiction; therefore, $w\u2217$ must be a fixed point. (Compact limit sets are connected; therefore, $\Omega +(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$.

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

$lim|w|\u21921\Phi (w)=\u2212\u221e.$

Under the conditions above, almost all trajectories for (1) converge to $\Delta $ as $t\u2192\u221e$ and to $\Sigma $ as $t\u2192\u2212\u221e$.

Let $p=(p1,\u2026,pN)\u2208X$ be any point with all distinct coordinates. The points on $Gp$ are parameterized by $w\u2208Bd$ and $\zeta \u2208SO(d)$, and the dynamics for these parameters are given by (5). We begin with the dynamics as $t\u2192\u2212\u221e$. Let $w(t)$ be a trajectory for (14) with initial condition $w0\u2208Bd$ and consider the backward-time limit set $\Omega \u2212(w0)$; this is a nonempty, compact, connected subset of $Bd\xaf$. The potential $\Phi $ is decreasing along all trajectories $w(t)$, hence bounded below as $t\u2192\u2212\u221e$; therefore, Lemma 2 implies that the limit set $\Omega \u2212(w0)$ must be contained in the interior $Bd$. We know that any $w\u2217\u2208\Omega \u2212(w0)$ is a fixed point for the flow. By Lemma 1, $w\u2217$ is repelling, and therefore, any trajectory $w(t)$ that comes sufficiently close to $w\u2217$ must have $w(t)\u2192w\u2217$ as $t\u2192\u2212\u221e$; therefore, $\Omega \u2212(w0)={w\u2217}$. This proves the existence of fixed points for (14) and that every trajectory $w(t)$ converges to a fixed point as $t\u2192\u2212\u221e$. 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\u2217$, and $w(t)\u2192w\u2217$ as $t\u2192\u2212\u221e$ for all trajectories. The fixed point $w\u2217$ has $Z(Mw\u2217(p))=0$; therefore, all trajectories in $Gp$ converge to $\Sigma $ as $t\u2192\u2212\u221e$.

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\u2217=0,u\u2217=p1$, which has a unique attracting trajectory because it is a saddle with a $(d\u22121)$-dimensional unstable manifold.

## IX. SUMMARY AND DISCUSSION

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.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

## REFERENCES

*International Symposium on Mathematical Problems in Theoretical Physics*(Springer, 1975), pp. 420–422.

*et al.*, “

*Proceedings of the 45th IEEE Conference on Decision and Control*(IEEE, 2006), pp. 5060–5066.