The Kuramoto model on a sphere: Explaining its low-dimensional dynamics with group theory and hyperbolic geometry

TL;DR

This paper explains the low-dimensional dynamics of a generalized Kuramoto model on a sphere using group theory and hyperbolic geometry, unifying finite and infinite particle cases and analyzing stability.

Contribution

It introduces a group-theoretic and geometric framework to understand the low-dimensional behavior of the sphere-based Kuramoto model, including the continuum limit and stability analysis.

Findings

Low-dimensional dynamics explained via hyperbolic geometry and group theory.

Unified understanding of finite and infinite particle limits.

Global stability results for synchronized states under certain couplings.

Abstract

We study a system of $N$ 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 \ge 3$, and to clarify why it admits the analog of the Ott-Antonsen ansatz in the continuum limit $N \rightarrow \infty$. The underlying reason is that the system is intimately connected to the natural hyperbolic geometry on the unit ball $B^d$. In this geometry, the isometries form a Lie group consisting of higher-dimensional generalizations of the M\"obius 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.

Equations

$${\dot{x}}_i=A_i{x}_i+Z-⟨Z,x_i⟩x_i,i=1,\dots ,N,$$

$$Z=\sum _{i=1}^N{a}_i{x}_i$$

$$\dot{x}=Ax+Z-⟨Z,x⟩x,$$

$$ds=\frac{2|dx|}{1-|x{|}^2},$$

$$M_w(x)=\frac{x-w}{1-\overline{w}x},$$

$$\frac{(x-w)(1-w\overline{x})}{(1-\overline{w}x)(1-w\overline{x})}$$

$$=\frac{x-w-w|x{|}^2+w(2⟨w,x⟩-\overline{w}x)}{1-2⟨w,x⟩+|w{|}^2|x{|}^2}$$

$$=\frac{(1-|w{|}^2)x-(1-2⟨w,x⟩+|x{|}^2)w}{1-2⟨w,x⟩+|w{|}^2|x{|}^2}.$$

$$M_w(x)$$

$$=\frac{(1-|w{|}^2)(x-|x{|}^{2}w)}{1-2⟨w,x⟩+|w{|}^2|x{|}^2}-w,$$

$$M_w(x)=\frac{(1-|w{|}^2)(x-w)}{|x-w{|}^2}-w.$$

$$M_w(w)$$

$$=\frac{(1-|w{|}^2)w-(1-|w{|}^2)w}{1-2⟨w,w⟩+|w{|}^2|w{|}^2}=0.$$

$$M_w(w)=\frac{(1-|w{|}^2)(w-|w{|}^{2}w)}{1-2⟨w,w⟩+|w{|}^2|w{|}^2}-w=\frac{(1-|w{|}^2{)}^{2}w}{(1-|w{|}^2{)}^2}-w=0.$$

$$g(x)=\zeta M_w(x)$$

$$g(x)=M_{-z}(\xi x),$$

$$g(x)\approx \zeta \left(-w+(1-|w{|}^2)x\right)\approx z+(1-|z{|}^2)\xi x.$$

$$\dot{y}=Ay-⟨Z,y⟩y+\frac{1}{2}(1+|y{|}^2)Z,$$

$$M_{tw}(x)\approx \frac{x-|x{|}^{2}tw}{1-2t⟨w,x⟩}-tw\approx x+t\left(2⟨w,x⟩x-(1+|x{|}^2)w\right).$$

$$y(t)=g_t(g_0^{-1}(g_0(y_0)))=g_t(y_0)=\zeta M_w(y_0)=M_{-z}(\zeta y_0)$$

$$y\approx z+(1-|z{|}^2)\zeta y_0,$$

$$\dot{y}\approx \dot{z}-2⟨\dot{z},z⟩\zeta y_0+(1-|z{|}^2)\dot{\zeta }{y}_0.$$

$$\dot{y}$$

$$\approx Az+\frac{1}{2}(1+|z{|}^2)Z-⟨Z,z⟩z$$

$$\displaystyle+(1-|z|^{2})\Bigl{(}A\zeta y_{0}+\langle z,\zeta y_{0}\rangle Z-\langle Z,z\rangle\zeta y_{0}-\langle Z,\zeta y_{0}\rangle z\Bigr{)}.$$

$$\dot{z}=Az+\frac{1}{2}(1+|z{|}^2)Z-⟨Z,z⟩z$$

$$⟨\dot{z},z⟩=\frac{1}{2}(1-|z{|}^2)⟨Z,z⟩.$$

$$\dot{\zeta }{y}_0=A\zeta x_0+⟨z,\zeta y_0⟩Z-⟨Z,\zeta y_0⟩z.$$

$$\alpha (y_1,y_2)y=⟨y_1,y)y_2-⟨y_2,y)y_1;$$

$$\dot{\zeta }{y}_0=A\zeta y_0+\alpha (z,Z)\zeta y_0$$

$$\dot{\zeta }=(A+\alpha (z,Z))\zeta .$$

$$Az+\frac{1}{2}(1+|z{|}^2)Z-⟨Z,z⟩z=-\zeta \dot{w}-\dot{\zeta }w$$

$$\zeta \dot{w}$$

Full text

††preprint: AIP/12https://www.overleaf.com/project/607b700b5b6aed13ffadacf13-QED

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-1)$-dimensional sphere in $d$-dimensional Euclidean space, yielding a state space of dimension $N(d-1)$. 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\geq 3$ 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 further structure (hyperbolic gradient dynamics) that force 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. Kuramoto (1975) Kuramoto’s model displayed many remarkable features: It was exactly solvable (at least in some sense), despite being nonlinear and infinite-dimensional. Kuramoto (1984) 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. Winfree (1967, 1980) 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. Strogatz (2000); Pikovsky, Rosenblum, and Kurths (2003); Strogatz (2003); Acebrón et al. (2005); Dörfler and Bullo (2014); Pikovsky and Rosenblum (2015); Rodrigues et al. (2016); Bick et al. (2020)

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 Wiesenfeld, Colet, and Strogatz (1996, 1998), and themselves displayed remarkable dynamical features, such as surprisingly low-dimensional invariant toriTsang et al. (1991); Swift, Strogatz, and Wiesenfeld (1992) and ubiquitous neutral stability of splay statesNichols and Wiesenfeld (1992), 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 Watanabe and Strogatz (1993, 1994), which showed that the governing equations have $N-3$ constants of motion for all $N\geq 3$. Goebel Goebel (1995) 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, even further restricted to junctions that were strictly identical.

A breakthrough occurred in 2008 with the work of Ott and Antonsen. Ott and Antonsen (2008, 2009) 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. Ott and Antonsen (2008) 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. Pikovsky and Rosenblum (2015)

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 also was not clear whether they were connected or perhaps even equivalent. There were reasons to doubt they were linked: the Watanabe-Strogatz transformation could be used for any finite $N\geq 3$, 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. Pikovsky and Rosenblum (2015, 2008); Marvel, Mirollo, and Strogatz (2009); Stewart (2011); Chen, Engelbrecht, and Mirollo (2017, 2019) 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.Marvel, Mirollo, and Strogatz (2009); Chen, Engelbrecht, and Mirollo (2017)

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. Lohe (2009); Tanaka (2014); Chi, Choi, and Ha (2014); Ha et al. (2016); Ha, Ko, and Ryoo (2018); Lohe (2018); Jaćimović and Crnkić (2018); Chandra, Girvan, and Ott (2019a, b); Lohe (2019); DeVille (2019); Ha et al. (2021); Jaćimović and Crnkić (2021); Dai et al. (2021) In particular, several researchers have explored a generalization of the Kuramoto model in which the oscillators move on spheres instead of the unit circle. These spheres could be either the ordinary two-dimensional sphere or higher-dimensional spheres. A counterpart of the Ott-Antonsen ansatz has been discovered for the continuum version of the Kuramoto model on the $d$-dimensional sphere and used to reduce its infinite-dimensional dynamics to a lower dimensional set of ordinary differential equations (ODEs). Chandra, Girvan, and Ott (2019b) But 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. 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 Section 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. Instead, 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 Section VIII.

II Preliminaries

II.1 The Kuramoto Model on a Sphere

In a pioneering paper, Lohe Lohe (2009) 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 modelKuramoto (1975, 1984) with complex numbers $\exp(i\theta_{j})$ on the unit circle and then views those as equivalent to $2\times 2$ rotation matrices parametrized 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 $S^{d-1}$ in $\mathbb{R}^{d}$. 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

[TABLE]

where $x_{i}$ is a point on the unit sphere $S^{d-1}\subset{\mathbb{R}}^{d}$, each $A_{i}$ is an antisymmetric $d\times d$ matrix, and $Z\in{\mathbb{R}}^{d}$ is a $d$-dimensional vector analogous to the complex order parameter for the classic Kuramoto model. In Eq. (1), the matrix $A_{i}$ and the vector $Z$ are functions of the configuration $(x_{1},\dots,x_{N})$ 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” $x_{i}$ together. The antisymmetric matrix $A_{i}$ 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 $\langle x_{i},\dot{x_{i}}\rangle=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=(S^{d-1})^{N}$, which has dimension $N(d-1)$. Later we will also consider the natural infinite-$N$ analogue of (1), where a state is a probability measure on $S^{d-1}$.

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, like a linear combination of the form

[TABLE]

where the $a_{i}$ are real constants.

II.2 General philosophy: Lie groups and reducible systems

There is a general technique for dimensional reduction of systems like (1) which 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 $\dot{\gamma}(0)=v$, where $v$ is a vector in the tangent space $T_{e}G$ of $G$ at $e$. This vector $v$ is in turn associated very naturally with a corresponding vector $\tilde{v}$ in the tangent space of $X$, as follows. For each $x\in X$, $t\mapsto\gamma(t)x$ is a smooth curve in $X$, and its derivative at $t=0$ defines a vector $\tilde{v}_{x}$ in the tangent space at $x$. The vector field $\tilde{v}=(\tilde{v}_{x})$ is the infinitesimal generator corresponding to the element $v$ in the tangent space of $G$ at $e$. At each point $x\in X$, the infinitesimal generators $v_{x}$ span a linear subspace $V_{x}$ of the tangent space $T_{x}X$, 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}\in V_{x}$ at each point $x\in X$, 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 the state space $X$, then this will give us a dimensional reduction of the dynamics from $\dim X$ to $\dim G$. Now suppose, as is the case in applications of this methodology, the correspondence $G\to Gx$ is one-to-one for generic $x\in X$; equivalently, the stabilizer subgroup $G_{x}=\{e\}$ for generic $x\in X$. Then for each $x\in X$, 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 ${\mathbb{R}}^{d}$. (For an intuitive picture, think of these isometries as rotations.) Then $G$ acts on ${\mathbb{R}}^{d}$, 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=({\mathbb{R}^{d}})^{N}$ of $N$-tuples $x=(x_{i})$, $x_{i}\in{\mathbb{R}}^{d}$, and then the infinitesimal generators have the form $(\nu_{x})_{i}=Ax_{i}$ for some skew-symmetric matrix $A$. We could also, if we like, restrict $X$ to $N$-tuples $(x_{i})$ with $x_{i}\in S^{d-1}$, the state space for (1). Now suppose we had a dynamical system on $X$ of the form $\dot{x}_{i}=Ax_{i}$, where the skew-symmetric matrix $A$ is a function of the configuration $x=(x_{i})$; 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-1)/2$, and for large $N$ this is much smaller than the dimension of $X$, which is $N(d-1)$. Basically the configuration $(x_{i})$ of points on the sphere $S^{d-1}$ can only move collectively by a rotation of the sphere, so the dynamics reduces to a dynamical system on $SO(d)$.

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

[TABLE]

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

So for the Kuramoto system (1), if the matrices $A_{i}$ 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=\sum_{i=1}^{N}a_{i}x_{i}$ with $a_{i}>0$. But first we need to show that vector fields on $S^{d-1}$ of the form in (2) are indeed the infinitesimal generators of the action of the larger Möbius group.

II.3 Hyperbolic geometry and Möbius transformations

In this paper, a Möbius transformation is a composition of Euclidean isometries and spherical inversions of $\mathbb{R}^{d}$ 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 $B^{d}$ with boundary $S^{d-1}$. This geometry has metric

[TABLE]

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 constant (sectional) curvature equal to $-1$, and we can describe its isometries, which generalize the Möbius transformations preserving the unit disc for $d=2$. For $d=2$, let $w\in B^{2}$ and consider the Möbius transformation

[TABLE]

which preserves the unit disc $B^{2}$ and its boundary $S^{1}$. To generalize this to higher dimensions, we need to express $M_{w}(x)$ without reference to complex arithmetic operations or conjugation. This is the goal of the next subsection.

II.3.1 Möbius transformations in higher dimensions

Using the identity $2\langle w,x\rangle=\overline{w}x+w\overline{x}$, we see that

[TABLE]

This form of $M_{w}$ generalizes to higher dimensions: Let $w\in B^{d}$ and define

[TABLE]

where $x\in B^{d}$ or $S^{d-1}$. We call $M_{w}$ a boost transformation.

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

[TABLE]

Now see that

[TABLE]

Alternatively, we can use the second form to show

[TABLE]

Similar computations show that $M_{0}$ is the identity, $M_{w}^{-1}=M_{-w},$ and $M_{w}(0)=-w$.

It is a standard result in hyperbolic geometry (e.g., see Theorem 3.5.1 in Beardon Beardon (1983)) that any orientation-preserving isometry of $B^{d}$ 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

[TABLE]

and also uniquely in the form

[TABLE]

for some vectors $w,z\in B^{d}$ and rotations $\zeta,\,\xi\in SO(d)$, where $SO(d)$ denotes the group of orientation-preserving orthogonal linear transformations on ${\mathbb{R}}^{d}$. So 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-1)/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$:

[TABLE]

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

II.3.2 Infinitesimal generators

Having parametrized the Möbius transformations, we are now ready to derive the associated infinitesimal generators of the Möbius group action on the ball $B^{d}$. We will show that they correspond to flows of the form

[TABLE]

where $A$ is an antisymmetric $d\times d$ matrix and $Z\in{\mathbb{R}}^{d}$ is a vector. Note that, as advertised, this flow extends to a flow on $S^{d-1}$ 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 $M_{tw}(x)$ to first order in $t$:

[TABLE]

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 $\dot{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 $A_{i}$ identical. Therefore, by the general philosophy discussed earlier, the evolution of any initial point $p\in X$ under the system (1) with all $A_{i}=A$ lies in the group orbit $Gp$.

III Reduced Equations

The given Kuramoto system has $N(d-1)$ 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)\cong B^{d}$. 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.

Now suppose all the terms $A_{i}$ in (1) are equal. Fix a base point $p=(p_{1},\dots,p_{N})\in X$. Then if the points $p_{i}$ are in sufficiently general position, every element in the $G$-orbit of $p$ can be expressed uniquely as $gp$ for some $g\in G$, with parameters $w,z$ and $\zeta$. We wish to derive the corresponding evolution equations for $w,z$ and $\zeta$. Let $(x_{i}(t))$ be any solution to (1) in the group orbit $Gp$; we do not require that the initial point $(x_{i}(0))=p$. Then for $i=1,\dots,N$ we have $x_{i}(t)=g_{t}(p_{i})$ for a unique $g_{t}\in G$, which determines the parameters $w,z,\zeta$ as functions of $t$. Now consider the equation (3), with coefficients $A$ and $Z$ evaluated at $(x_{i}(t))$. This is a non-autonomous ODE on $\overline{B^{d}}$, and its time-$t$ flow must be given by some $\tilde{g}_{t}\in G$. This ODE has solutions $x_{i}(t)=g_{t}(p_{i})=g_{t}(g_{0}^{-1}(x_{i}(0)))$, which implies that $\tilde{g}_{t}=g_{t}g_{0}^{-1}$.

So for any $y_{0}\in\overline{B^{d}}$,

[TABLE]

must satisfy the ODE (3) with $A$ and $Z$ evaluated at $(x_{i}(t))$ at time $t$. In particular, if we let $y_{0}=0$, then $y(t)=-\zeta w=z$, so $z$ satisfies the ODE (3).

Now expand $y=\zeta M_{w}(y_{0})=M_{-z}(\zeta y_{0})$ to first order in $y_{0}$, using the variables $z$ and $\zeta$:

[TABLE]

so

[TABLE]

On the other hand, (3) gives

[TABLE]

Setting $y_{0}=0$ gives the $\dot{z}$ equation

[TABLE]

as expected, and since $\langle Az,z\rangle=0$, this in turn implies that

[TABLE]

Equating the $y_{0}$ terms, factoring out $1-|z|^{2}$ and canceling the common term $\langle Z,z\rangle\zeta y_{0}$ gives

[TABLE]

Together, the last two terms above define a special type of antisymmetric operator of $\zeta y_{0}$: Given any $y_{1},y_{2}\in{\mathbb{R}}^{d}$, define the antisymmetric operator $\alpha$ as

[TABLE]

this operator has range = ${\rm span}(y_{1},y_{2})$ providing $y_{1}$ and $y_{2}$ are linearly independent; otherwise $\alpha(y_{1},y_{2})=0$. Then for all $y_{0}\in{\mathbb{R}}^{d}$,

[TABLE]

and therefore

[TABLE]

Differentiating $z=-\zeta w$ gives

[TABLE]

so

[TABLE]

hence

[TABLE]

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

[TABLE]

with $A$ and $Z$ evaluated at $M_{-z}(\zeta p)$, and for the $(w,\zeta)$ coordinate system on $Gp$ are

[TABLE]

with $A$ and $Z$ evaluated at $\zeta M_{w}(p)$. Note that these equations generalize the evolution equations for the parameters $w$ and $\zeta$ given in Chen et al.Chen, Engelbrecht, and Mirollo (2017) for the classic case $d=2$.

IV Comparison of Z Versus W Coordinates

The $\dot{z}$ equation (4) is an extension of the system equation (1) on $S^{d-1}$. However, for finite $N$, the $\dot{z}$ equation does not uncouple from $\zeta$, since $Z$ is evaluated at $M_{-z}(\zeta p)$. The exception to this is in the infinite-$N$ limit: if the base point $p$ is now the uniform density on $S^{d-1}$, then $\zeta p=p$ (the uniform density is invariant under rotations) and the density $M_{-z}(p)$ is a hyperbolic Poisson density on $S^{d-1}$ 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\geq 3$ (we will give more details on this in the next section).

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

[TABLE]

with $a_{i}\in\mathbb{R}$, $\zeta$ drops out of the $\dot{w}$ equation and we get the reduced equation

[TABLE]

The parameter $w$ essentially defines the “phase relations” among the $x_{i}$; two configurations have the same $w$ if and only if they are related by a rotation. So $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^{\prime}=M(p)\in Gp$; then we have coordinates $w^{\prime},\zeta^{\prime}$ associated to the base point $p^{\prime}$. Any $q\in Gp$ has two expressions

[TABLE]

Assuming the coordinates of $p$ are in sufficiently general position, this implies $\zeta M_{w}=\zeta^{\prime}M_{w^{\prime}}\circ M$, and hence

[TABLE]

But the unique solution to $M_{w^{\prime}}(y)=0$ is $w^{\prime}$, and hence $w^{\prime}=M(w)$. In other words, the coordinates $w$ and $w^{\prime}$ transform exactly as the base points $p$ and $p^{\prime}$.

V Continuum Limit

Next, we consider the dynamics of the Kuramoto model (1) in the limit $N\to\infty$. Let us assume first that the rotation terms $A_{i}=A$ are constant across the population, corresponding to identical “oscillators.” Later we will consider the case where $A$ varies depending on some distribution.

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

[TABLE]

In the continuum limit, a state of the system is a probability measure $\rho$ on $S^{d-1}$, and the order parameter becomes

[TABLE]

The measure $\rho$ evolves according to the continuity equation (also known as the noiseless Fokker-Planck equation) associated to the flow in (1). Naturally, this flow must preserve group orbits under the action of $G$. Recall that if $M\in G$, then the measure $M_{\ast}\rho$ is defined by the adjunction formula

[TABLE]

In particular, we can consider the $G$-orbit of the uniform probability measure $\sigma$ on $S^{d-1}$. 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_{-z})_{\ast}\sigma$, with $z\in B^{d}$. The evolution equation for $z$ is (4), with

[TABLE]

In the case $d=2$ with $x=\zeta\in S^{1}$, we have

[TABLE]

so the integral

[TABLE]

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

[TABLE]

when $d=2$. Unfortunately, the formula $Z(z)=Kz$ is not correct for $d\geq 3$; 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 the next section. 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 to its metric; functions $f$ on $X$ satisfying the equation $\Delta f=0$ are called harmonic. For functions on the ball $B^{d}$ with the hyperbolic metric, this operator is

[TABLE]

where

[TABLE]

is the standard Laplace operator (see StollStoll , Chapter 3). We will call solutions to the equation $\Delta_{hyp}f=0$ hyperbolic harmonic functions; for $d=2$ these coincide with ordinary (Euclidean) harmonic functions. We can consider the hyperbolic analogue of the classical Dirichlet problem: given a continuous function $f$ on $S^{d-1}$, extend $f$ to a hyperbolic harmonic function $\tilde{f}$ on $B^{d}$. Assuming this problem has a unique solution, then for any rotation $\zeta\in SO(d)$ we must have $\widetilde{f\circ\zeta}=\tilde{f}\circ\zeta$, since rotations preserve the hyperbolic metric. If we average $f\circ\zeta$ on $S^{d-1}$ over all rotations $\zeta\in SO(d)$ we get the constant function

[TABLE]

on $S^{d-1}$, and any constant is hyperbolic harmonic on $B^{d}$. Therefore the average on $B^{d}$ of $\widetilde{f\circ\zeta}=\tilde{f}\circ\zeta$ over all $\zeta\in SO(d)$ must be the constant $f_{ave}$. But $\tilde{f}(\zeta(0))=\tilde{f}(0)$ for all $\zeta$, so we must have

[TABLE]

Now let $z\in B^{d}$; since $M_{-z}$ preserves the hyperbolic metric, we must have $\widetilde{f\circ M}_{-z}=\tilde{f}\circ M_{-z}$, which implies

[TABLE]

As shown in Chapter 5 in Stoll Stoll , the measure $(M_{-z})_{\ast}\sigma$ is given by the formula

[TABLE]

with hyperbolic Poisson kernel function

[TABLE]

Thus the solution to the hyperbolic Dirichlet problem with boundary function $f$ on $S^{d-1}$ is given by the hyperbolic Poisson integral

[TABLE]

The orbit $G\sigma$ consists of all hyperbolic Poisson measures $P(z,x)\,d\sigma(x)$, parametrized by $z\in B^{d}$. By contrast, the Euclidean Poisson kernel function is

[TABLE]

so 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\geq 2$. We see from (8) that $Z(z)$ is the hyperbolic Poisson integral of the function $Kx$ on $S^{d-1}$. The function $Kx$ is (Euclidean) harmonic and homogeneous of degree 1 on ${\mathbb{R}}^{d}$; following the recipe in Chapter 5 in Stoll Stoll , we see that its extension from $S^{d-1}$ to a hyperbolic harmonic function on $B^{d}$ is given by

[TABLE]

where $F$ is the hypergeometric function

[TABLE]

with $(a)_{0}=1$ and $(a)_{k}=a(a+1)\cdots(a+k-1)$ for $k\geq 1$. 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 ${\mathbb{R}}^{d}={\mathbb{C}}^{m}$, and we can study systems of the form

[TABLE]

where now $x_{i}$ is a point on the unit sphere $S^{2m-1}\subset{\mathbb{C}}^{m}$, $A_{i}$ is an anti-Hermitian $m\times m$ complex matrix, $Z\in{\mathbb{C}}^{m}$ and $\langle,\rangle$ 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\geq 2$. To see this, suppose

[TABLE]

for all $x\in S^{2m-1}\subset{\mathbb{C}}^{m}={\mathbb{R}}^{d}$, where $A$ is antisymmetric, $B$ is anti-Hermitian, $Y,Z\in{\mathbb{C}}^{m}$ and we use the subscripts $\mathbb{R}$ and $\mathbb{C}$ to distinguish the real and complex inner products. Then

[TABLE]

and so $(A-B)(-x)=(A-B)x$ for all $x\in S^{2m-1}$, which implies $A=B$. This implies

[TABLE]

for all $x\in S^{2m-1}$, hence $Y-Z\in{\rm span}_{\mathbb{C}}(x)$ for all $x\in{\mathbb{C}}^{m}$; if $m\geq 2$, this implies $Y=Z$. But then we have

[TABLE]

for all $x\in{\mathbb{C}}^{m}$, which can only hold if $Y=0$. Hence for $m\geq 2$, 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 $B^{m}$ with the Bergman metric (see RudinRudin (1980), Chapter 1). The orientation-preserving isometries of this metric are generated by unitary transformations $\zeta\in U(m)$ and boost transformations of the form

[TABLE]

Notice that when $m=1$, this reduces to the standard complex Möbius map $M_{w}$. As in the real case $M_{0}$ is the identity, $M_{w}^{-1}=M_{-w}$, $M_{w}(w)=0$ and $M_{w}(0)=-w$. Any orientation-preserving isometry of $B^{d}$ can be expressed uniquely in the form

[TABLE]

where $w,z\in B^{m}$ but now $\zeta,\,\xi\in U(m)$, the complex unitary group. Linearizing at $x=0$ gives

[TABLE]

which implies $z=-\zeta w$ (hence $|z|=|w|$) and $\xi=\zeta$, as before.

The corresponding infinitesimal transformations are given by flows on the complex unit ball $B^{m}$ of the form

[TABLE]

with $A$ anti-Hermitian $m\times m$ and $Z\in{\mathbb{C}}^{m}$. This flow extends to a flow on $S^{2m-1}$ of the form in (11). Note the absence of the quadratic term $|y|^{2}Z$ here. To derive these infinitesimal transformations, we can apply our prior power series expansion, noting that $|x|=1$ to obtain

[TABLE]

So the infinitesimal generator is an “infinitesimal boost” of the form (12) with $Z=-\frac{1}{2}w$ and $A=0$. The infinitesimal generators corresponding to the rotation components are flows of the form $\dot{x}=Ax$ with $A$ 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. Lohe (2009); Tanaka (2014); Chi, Choi, and Ha (2014); Ha, Ko, and Ryoo (2018); Lohe (2018); Jaćimović and Crnkić (2018); Chandra, Girvan, and Ott (2019a, b); Lohe (2019); Ha et al. (2021); Jaćimović and Crnkić (2021) Three papers in particular overlap considerably with the present work.

Tanaka Tanaka (2014) 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, Mirollo, and Strogatz foundMarvel, Mirollo, and Strogatz (2009) 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 $\dot{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 $\dot{w}$ equation (5). Tanaka also notes that 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 reductionOtt and Antonsen (2008) for the complex version of the system.

Lohe Lohe (2018) also looks at the same system as (1) (see his equation (22)) and he derives a similar reduction as ours by using Möbius transformations for the finite-$N$ model. His transformation (30) on $S^{d-1}$ is our $M_{w}$ (with $v=w$) and his equation (31) is the same as our $\dot{z}$ equation (4). He also has something that looks like the $\dot{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

[TABLE]

where $Q_{i}\in O(d)$ and $\lambda_{i}\in\mathbb{R}$. But such a $Z$ does not satisfy the identity $\zeta Z(p)=Z(\zeta p)$ for all rotations $\zeta$, unless $Q_{i}=\pm I$, so we do not see how the $\zeta$ term cancels.

Lohe’s map $M$ in his equation (55) (ignoring the $R$ factor) agrees with our map $M_{-v}$ on the sphere $S^{d-1}$, but not on the ball $B^{d}$. So it is not a Möbius transformation of the type we are using. For example, $M(-v)=v$ whereas $M_{-v}(-v)=0$. We are not sure why Lohe Lohe (2018) 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 equation (63) (again ignoring the $R$ factor) is exactly our $M_{-v}$.

Chandra, Girvan, and Ott Chandra, Girvan, and Ott (2019b) concentrate on the infinite-$N$ or continuum limit system, and derive a dynamical reduction for a special class of probability densities on $S^{d-1}$, 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 equation (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 equation (15) is exactly the same as our $\dot{z}$ equation (4) in the infinite-$N$ limit. The integral in their equation (19) can be evaluated, as shown above in (10).

VIII An Example: First-Order Linear Order Parameter Gives Gradient System

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

[TABLE]

where the $a_{i}$ are real constants.

VIII.1 A Computer Visualization

Underlying all of the mathematics presented is a system of particles on a sphere which coalesce when certain conditions are met. To visualize this, we implemented the Runge-Kutta algorithm to numerically solve the Kuramoto model on the 2-dimensional sphere in 3-space, corresponding to the case $d=3$ in (1). For simplicity, we simulated $N=100$ particles with equal weights ($a_{i}=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 this synchronization behavior holds more generally for Kuramoto models on the sphere having weighted order parameters of the form (13), provided the weights $a_{i}$ are all positive and satisfy an upper bound. Specifically, the forward-limit synchronization behavior remains regardless of how the particles are weighted, provided the $a_{i}$ sum to $1$ and no individual weight exceeds $1/2$, as we later prove.

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 backwards, almost all initial conditions of the particles will tend towards 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 which exceeds the weight of all other particles. When this occurs, it is impossible to arrange the particles so their weighted centroid is at the origin, so the backwards time limit will tend towards an antipodal configuration, where all particles not in the majority cluster will coalesce around the antipode of the cluster. This is the configuration which minimizes the magnitude of the weighted centroid. We do not include the proof that the backwards-time limit is antipodal in this case, but it is a straightforward generalization of the result which we do prove below.

VIII.2 Existence of Hyperbolic Gradient

As mentioned above, if $Z$ has the form in (13), then the $\dot{w}$ equation in (5) reduces to

[TABLE]

independent of the parameter $\zeta$. We will show that this is a gradient flow on the unit ball $B^{d}$ 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 $B^{d}$ (or any open subset of ${\mathbb{R}}^{d}$) and standard coordinates $w_{1},\dots,w_{d}$, the $1$-form associated to the vector field with components $f_{1},\dots,f_{d}$ is

[TABLE]

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

[TABLE]

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

[TABLE]

where $\nabla_{euc}$ denotes the ordinary Euclidean gradient operator. We have $\phi(w)=2(1-|w|^{2})^{-1}$ for the hyperbolic metric, so the hyperbolic gradient operator on $B^{d}$ is given by

[TABLE]

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

[TABLE]

where $p_{1,\,j}$ denotes the $j$th component of the point $p_{1}\in S^{d-1}$. Let $E_{j}$ denote the coefficient of $dw_{j}$ in parentheses above; then

[TABLE]

Applying the chain and quotient rules gives

[TABLE]

for $j\neq k$, 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

[TABLE]

Here we follow the convention that the potential decreases along trajectories, so we are asserting that $\nabla_{hyp}\Phi=-V$. To derive this, use the identity $\nabla_{euc}|w-w_{0}|^{2}=2(w-w_{0})$, for any constant vector $w_{0}\in{\mathbb{R}}^{d}$. Then

[TABLE]

Hence we see that

[TABLE]

as desired.

VIII.3 Analysis of Dynamics

We can use the existence of the potential $\Phi(w)$ for the flow on $B^{d}$ to prove a global synchrony result for the system (1) when the coefficients $a_{i}$ in the order parameter $Z$ are all positive. Specifically, we assume that $0<a_{i}<1/2$ for all $i$, and $\sum_{i=1}^{N}a_{i}=1$. We also assume $N\geq 3$ and all the rotation terms $A_{i}$ in (1) are equal. Under these conditions, almost all trajectories for (1) converge in forward time to the $(d-1)$-dimensional diagonal manifold $\Delta\subset X$ as $t\to\infty$, meaning that the system self-synchronizes. In contrast, in backwards time the system tends to an incoherent state having zero order parameter: as $t\to-\infty$ almost all trajectories for (1) converge to the codimension-$d$ subspace $\Sigma\subset X$ consisting of states with $Z(p)=0$.

The proof is modeled after Theorem 1 in Chen et al. Chen, Engelbrecht, and Mirollo (2019) and will be based on two preliminary lemmas. In each of these lemmas we assume the conditions on the $a_{i}$ above, and that the base point $p=(p_{i})$ for the flow (14) has all distinct coordinates.

We begin with a general observation about gradient flows in the ball $B^{d}$: if $w_{0}\in B^{d}$ is any initial condition and $w^{\ast}\in B^{d}$ is in the forward limit set $\Omega_{+}(w_{0})$, then $w^{\ast}$ is a fixed point for the flow. To see this, let $\Phi$ be a potential for the flow, and suppose $w(t_{n})\to w^{\ast}\in B^{d}$ for some sequence $t_{n}\to\infty$. Since the potential decreases along trajectories,

[TABLE]

Let $F_{t}$ denote the time-$t$ flow map. If $w^{\ast}$ is not a fixed point, then for any $s>0,$

[TABLE]

which is a contradiction, so $w^{\ast}$ must be a fixed point. (Compact limit sets are connected, so $\Omega_{+}(w_{0})$ cannot consist of two or more but finitely many fixed points; however it is possible that forward or backward limits sets for gradient flows consist of a continuum of fixed points. We will see that this is not the case for our system on $B^{d}$.)

Lemma 1**.**

Any fixed point for the flow (14) in $B^{d}$ is repelling.

Proof.

Suppose $w^{\ast}\in B^{d}$ 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^{\ast}=0$ without loss of generality, so $Z(p)=\sum\limits_{i=1}^{N}a_{i}p_{i}=0$. To first order in $w$,

[TABLE]

The linearization of (14) at the fixed point $w^{\ast}=0$ is

[TABLE]

We claim that the linear map

[TABLE]

has $||T||<1$; to see this, suppose $|w|=1$. Then $|\langle w,p_{i}\rangle p_{i}|\leq 1$ and $Tw$ is a convex combination of the vectors $\langle w,p_{i}\rangle p_{i}$. We can only obtain $|Tw|=1$ if all terms $\langle w,p_{i}\rangle p_{i}=u$ with $|u|=1$, which implies all $p_{i}=\pm u$, and this cannot happen if at least three of the $p_{i}$ are distinct. Hence $||T||<1$ and so the eigenvalues $\mu_{i}$ of $T$ satisfy $|\mu_{i}|<1$. The eigenvalues for the $\dot{w}$ linearization are $\lambda_{i}=1-\mu_{i}$, so wee see that ${\rm Re}\,\lambda_{i}>0$ for all $i$, establishing that the fixed point $w^{\ast}$ is repelling. ∎

Lemma 2**.**

$\displaystyle\lim_{|w|\to 1}\Phi(w)=-\infty.$**

Proof.

It suffices to show that

[TABLE]

for any sequence $w_{n}\in B^{d}$ with $w_{n}\to x\in S^{d-1}$. The result is clear if $x\neq p_{i}$: as $n\to\infty$ the terms $|w_{n}-p_{i}|$ in the potential (15) are bounded away from [math], and $1-|w_{n}|^{2}\to 0$. So let’s say that $w_{n}\to p_{1}$. We rewrite $\Phi(w_{n})$ as

[TABLE]

The latter two terms above have finite limit as $n\to\infty$, so we focus on the first two terms. We have $1-|w_{n}|\leq|w_{n}-p_{1}|$, so

[TABLE]

as $n\to\infty$, which proves our result. Notice that we need the assumption $a_{i}<1/2$ for this argument. ∎

Theorem. Under the conditions above, almost all trajectories for (1) converge to $\Delta$ as $t\to\infty$ and to $\Sigma$ as $t\to-\infty$.

Proof.

Let $p=(p_{1},\dots,p_{N})\in X$ be any point with all distinct coordinates. The points on $Gp$ are parametrized by $w\in B^{d}$ and $\zeta\in SO(d)$, and the dynamics for these parameters are given by (5). We begin with the dynamics as $t\to-\infty$. Let $w(t)$ be a trajectory for (14) with initial condition $w_{0}\in B^{d}$, and consider the backward time limit set $\Omega_{-}(w_{0})$; this is a nonempty, compact, connected subset of $\overline{B^{d}}$. The potential $\Phi$ is decreasing along all trajectories $w(t)$, hence bounded below as $t\to-\infty$, so Lemma 2 implies that the limit set $\Omega_{-}(w_{0})$ must be contained in the interior $B^{d}$. We know that any $w^{\ast}\in\Omega_{-}(w_{0})$ is a fixed point for the flow. By Lemma 1, $w^{\ast}$ is repelling and so any trajectory $w(t)$ which comes sufficiently close to $w^{\ast}$ must have $w(t)\to w^{\ast}$ as $t\to-\infty$; therefore $\Omega_{-}(w_{0})=\{w^{\ast}\}$. This proves the existence of fixed points for (14), and that every trajectory $w(t)$ converges to a fixed point as $t\to-\infty$. If the flow had multiple fixed points, we would obtain a partition of $B^{d}$ into the disjoint open basins of repulsion of the fixed points, violating connectedness of the ball. Therefore (14) has a unique fixed point $w^{\ast}$, and $w(t)\to w^{\ast}$ as $t\to-\infty$ for all trajectories. The fixed point $w^{\ast}$ has $Z(M_{w^{\ast}}(p))=0$, so all trajectories in $Gp$ converge to $\Sigma$ as $t\to-\infty$.

In forward time, the limit set $\Omega_{+}(w_{0})$ for any $w_{0}\neq w^{\ast}$ must be completely contained in the boundary $S^{d-1}$, since the unique fixed point $w^{\ast}\in B^{d}$ is repelling. Suppose we remove the factor $(1/2)(1-|w|^{2})$ in the flow (14); the scaled vector field on $B^{d}$ given by

[TABLE]

has the same trajectories as the original flow, just with different time parametrizations. Observe that this scaled vector field extends smoothly to ${\mathbb{R}}^{d}-\{p_{i}\}$, and coincides with the radial vector field $x$ at any $x\in S^{d-1}$ with $x\neq p_{i}$. Therefore there is a unique trajectory passing through each point $x\in S^{d-1}$, flowing from the interior to the exterior of the sphere, as long as $x\neq p_{i}$. Consequently the original flow (14) has a unique trajectory $w(t)$ in $B^{d}$ with $w(t)\to x$ as $t\to\infty$, as long as $x\neq p_{i}$. This also shows that there is a neighborhood $U$ of $S^{d-1}-\{p_{i}\}$ such that if $w(t_{0})\in U$ for some $t_{0}$, then $w(t)\to x\neq p_{i}$ for some $x\in S^{d-1}$. So if $\Omega_{+}(w_{0})$ contains some $x\neq p_{i}$, then the trajectory $w(t)$ of $w_{0}$ must enter the neighborhood $U$, and therefore $w(t)\to x\in S^{d-1}$ as $t\to\infty$.

Since limit sets are connected, the only other possibility is $\Omega_{+}(w_{0})=\{p_{i}\}$ for some $i$; equivalently, $w(t)\to p_{i}$. We will show that there is a unique trajectory with this behavior for each $p_{i}$. Assuming this, we see that with $N+1$ exceptions, any trajectory $w(t)$ converges to a point $x\in S^{d-1}$ with $x\neq p_{i}$ (the exceptions are the $N$ trajectories converging to the base point coordinates $p_{i}$, and the fixed point trajectory $w^{\ast}$). The corresponding trajectory in $Gp$ has coordinates

[TABLE]

We have $|w(t)|\to 1$ and $|p_{i}-w(t)|$ is bounded away from [math] as $t\to\infty$, so $M_{w(t)}(p_{i})\to-x$ for each $i$ and therefore the trajectory $\zeta(t)M_{w(t)}(p)$ in $Gp$ converges to $\Delta$ as $t\to\infty$.

This analysis breaks down at $x=p_{i}$ because the scaled vector field above does not have a unique limit as $w\to p_{i}$; rather, its limit depends on the direction of the approach. To see this, write $w=p_{1}-ru$, where $0<r<1$ and $|u|=1$ (with this convention, $u=p_{1}$ corresponds to $w$ approaching $p_{1}$ radially). Then $|p_{1}-w|=r$ and

[TABLE]

so

[TABLE]

As $r\to 0$, the magnitude of this term is $2\langle p_{1},u\rangle$, which depends on the angle of approach given by $u$ (note that $\langle p_{1},u\rangle>0$ because $u$ points outwards at $p_{1}$).

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^{\ast}=0,u^{\ast}=p_{1}$, which has a unique attracting trajectory because it is a saddle with a $(d-1)$-dimensional unstable manifold.

We see that the scaled system has

[TABLE]

where the $O(r)$ term is a smooth function of $r$ and $u$ for $|r|<\epsilon=\min|p_{i}-p_{1}|$, $i\geq 2$. This condition insures that $|p_{i}-w|\geq|p_{i}-p_{1}|-|r|>0$, so the $i\geq 2$ terms in the scaled $\dot{w}$ equation are all smooth functions of $r$ and $u$. And we can allow $r<0$ here, even though it is not relevant to the $\dot{w}$ system. Now $r^{2}=|w-p_{1}|^{2}$, so

[TABLE]

which gives

[TABLE]

Differentiating $ru=p_{1}-w$ gives

[TABLE]

Hence the scaled system in polar form can be written

[TABLE]

We emphasize that the $O(r)$ and $O(r^{2})$ terms are smooth functions of $r,u$ as long as $|r|<\epsilon$. We consider the “semi-scaled” polar system

[TABLE]

which has the same trajectories as the original system, just with different time parametrizations. The advantage of this modified system is that the equations are smooth on $(-\epsilon,\epsilon)\times S^{d-1}$.

Observe that the system (17b) has $\{0\}\times S^{d-1}$ invariant, and has fixed point $(r^{\ast},u^{\ast})=(0,p_{1})$. The fixed point $p_{1}$ is repelling on the invariant manifold $\{0\}\times S^{d-1}$; to see this, observe that

[TABLE]

when $r=0$. In fact, if we assign the coordinate $\theta$ on any great circle joining $p_{1}$ and $-p_{1}$ on $S^{d-1}$ so that $u=e^{i\theta}$ and $p_{1}=1$, then the system reduces to $\dot{\theta}=\sin\theta$. We also see that the $\dot{r}$ equation linearized at $(0,p_{1})$ is $\dot{r}=-(1-2a_{1})r$, so the linearization of (17b) has the single negative eigenvalue $-(1-2a_{1})$ and $d-1$ positive eigenvalues $+1$. Therefore $(0,p_{1})$ is a saddle with a one-dimensional stable manifold, and hence has a unique trajectory $(r(t),u(t))\to(0,p_{1})$ with $r(t)>0$.

Now suppose we have a trajectory $w(t)\to p_{1}$ in our original system (14). The corresponding trajectory for (17b) will have $r(t)\to 0$; we cannot achieve $r(t)=0$ in finite time because the manifold $\{r=0\}$ is invariant for (17b). We must prove that $u(t)\to p_{1}$, so that this trajectory is in fact the saddle stable manifold. Observe that

[TABLE]

Also note that $\langle p_{1},u(t)\rangle>0$ since $|w(t)|<1$. Let $0<c<1$; then for some $T\geq 0$, $t\geq T$ implies $O(r(t))\leq(1-c^{2})/2$. Now suppose $0<\langle p_{1},u(t_{0})\rangle<c$ for some $t_{0}\geq T$; then $0<\langle p_{1},u(t)\rangle c$ for all $t\geq t_{0}$. This is because the function $t\mapsto\langle p_{1},u(t)\rangle$ is decreasing if $0<\langle p_{1},u(t)\rangle<c$:

[TABLE]

as long as $0<\langle p_{1},u(t)\rangle<c$. But this also implies that eventually $\langle p_{1},u(t)\rangle<0$, which is a contradiction. Hence we must have $\langle p_{1},u(t)\rangle\geq c$ for all $t\geq T$, which proves that $u(t)\to p_{1}$. ∎

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 $B^{d}$, 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 $B^{d}$. In particular, the case when the sum of the weights is [math] 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 disc given by placing a collection of point charges on the unit circle, with total charge [math], as shown by Chen et al. Chen, Engelbrecht, and Mirollo (2019). 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 $A_{i}$, 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.

Acknowledgements.

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

Data Availability

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