Higher-order synchronization on the sphere

We begin in Section 2 by discussing models of synchronization on the unit sphere which form a natural generalization of the widely-studied Kuramoto model. We formulate the equations as a gradient system, choosing firstly the familiar bilinear potential in Section 2.1 which leads to pairwise interactions, followed in Section 2.2 by a multilinear potential with antisymmetric couplings, which leads to $d$-body interactions on $\mathbb{S}^{d-1}$. Whereas the $2$-body equations have cubic nonlinearities with $N-1$ summations at each node, corresponding to $2$-body interactions with the other $N-1$ nodes, the $d$-body systems have nonlinearities of order $d+1$, with $(N-1)!/(N-d)!$ summations at each node, which interacts with the other nodes by means of a connected $(d-1)$-simplex. It is straightforward to combine the $2$-body and $d$-body systems, and so we investigate the relative effect of the two couplings with corresponding strengths $\kappa_{2},\kappa_{d}$. We discuss steady state solutions and their stability in Section 2.4 for both the $d$-body system and the combined $2$-body and $d$-body systems. The stable steady state configurations for $d$-body systems consist of equally spaced nodes on $\mathbb{S}^{d-1}$, rather than completely synchronized configurations as occur for $2$-body couplings. For combined systems, the $d$-body steady states prevail, unless the ratio of coupling constants $\kappa_{2}/\kappa_{d}$ is sufficiently large to favour the $2$-body forces.

In Section 3 we focus on the case $d=3$, this being the simplest of the higher-order models, although $d=2$ is also of interest but is considered later in Section 7 because it has only $2$-body interactions. We show that the $d=3$ system does not completely synchronize, as occurs for the $2$-body models on the sphere, rather for generic initial values the nodes arrange themselves asymptotically in a ring formation. Exactly why ring synchronization, rather than complete synchronization, occurs is discussed for the special case $N=3$ in Section 3.3, for which we can solve the nonlinear equations exactly. We consider combined $2$- and 3-body interactions in Section 4, obtaining an exact expression for the steady states which again form an equally spaced ring of nodes on $\mathbb{S}^{2}$. These states exist and are stable, unless the ratio of coupling constants $\kappa_{2}/\kappa_{3}$ becomes sufficiently large, at which point the system transitions continuously to a completely synchronized formation. We can state precisely the value of $\kappa_{2}/\kappa_{3}$ at which this occurs, see Section 4.2.

Properties of the general system depend on whether $d$ is even or odd, and so we have selected two even cases $d=2,4$ and two odd cases $d=3,5$ for detailed investigation. For even $d$, the equally spaced synchronized nodes do not form a closed sequence, as occurs for odd $d$, as we see in Section 5 for $d=4$. The asymptotic configurations for the combined $2$- and $4$-body systems, discussed in Section 5.2, are of a form that is not easily analyzed, but general arguments show that again synchronization is controlled by the 4-body interactions, unless the ratio $\kappa_{2}/\kappa_{4}$ becomes sufficiently large, then the attractive forces due to the $2$-body interactions force the system to completely synchronize. For the 5-body system, which is similar to the case $d=3$, we obtain exact steady state expressions which extend to the combined system with $2$-body interactions, see Section 6. In this case there is a discontinuous transition to a completely synchronized system at a critical ratio $\kappa_{2}/\kappa_{5}$.

Finally we consider the $d=2$ case in Section 7, which we have left to last because only $2$-body interactions are involved, however the system has properties which differ from those of the standard Kuramoto model, due to the antisymmetric coupling constants. There is merit in investigating this as the simplest of the $\mathbb{S}^{d-1}$ models and also as a guide to the behaviour of the systems for $d>2$. Again, in the combined system the antisymmetric couplings enhance the synchronization of the Kuramoto model.

The choice for $\mathfrak{V}_{d}$ which has been widely investigated is

$$\mathfrak{V}_{d}=\sum_{i,j=1}^{N}a_{ij}\,\bm{x}_{i}\centerdot\bm{x}_{j},$$

(1)

where the coefficients $a_{ij}$ are symmetric, and in this case the evolution equations are $\dot{\bm{x}_{i}}=-\lambda_{i}\bm{x}_{i}+\sum_{j}a_{ij}\bm{x}_{j}$ where $\lambda_{i}$ are Lagrange multipliers, supplemented by the constraint $\bm{x}_{i}\centerdot\bm{x}_{i}=1$. Hence $\lambda_{i}=\sum_{j}a_{ij}\bm{x}_{i}\centerdot\bm{x}_{j}$ and so we obtain the system of $N$ equations

$$\dot{\bm{x}_{i}}=\Omega_{i}\bm{x}_{i}+\frac{\kappa_{2}}{N}\sum_{j=1}^{N}a_{ij}\left[\bm{x}_{j}-(\bm{x}_{i}\centerdot\bm{x}_{j})\;\bm{x}_{i}\right],$$

(2)

where we have included $d\times d$ antisymmetric frequency matrices $\Omega_{i}$, as well as a normalized coupling coefficient $\kappa_{2}$ which measures the relative strength of the natural oscillations compared to the nonlinear interactions. Since $\bm{x}_{i}\centerdot\dot{\bm{x}_{i}}=0$, the unit length of $\bm{x}_{i}$ is maintained as the system evolves. The Kuramoto model corresponds to the case $d=2$ with the parametrization $\bm{x}_{i}=(\cos\theta_{i},\sin\theta_{i})$ and $\Omega_{i}=\pmatrix{0&-\omega_{i}\cr\omega_{i}&0}$, for which (2) reduces to $\dot{\theta_{i}}=\omega_{i}+\frac{\kappa_{2}}{N}\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i})$. In this case the potential (1) reduces to $\mathfrak{V}_{2}=\sum_{i,j}a_{ij}\cos(\theta_{j}-\theta_{i})$ which is well-known as a Lyapunov function for the Kuramoto model [18, 19]. It is convenient to write (2) in the form

$$\dot{\bm{x}_{i}}=\Omega_{i}\bm{x}_{i}+\bm{X}_{i}-\bm{x}_{i}\;(\bm{x}_{i}\centerdot\bm{X}_{i}),$$

(3)

where $\bm{X}_{i}=\kappa_{2}\sum_{j=1}^{N}a_{ij}\bm{x}_{j}/N$, a form which is maintained also for the $d$-body system. For global coupling, $a_{ij}=1$ for all $i,j$, every node connects to the average position $\bm{X}_{\mathrm{av}}$ defined by

$$\bm{X}_{\mathrm{av}}=\frac{1}{N}\sum_{j=1}^{N}\bm{x}_{j},$$

(4)

then we can write (2) as:

$$\dot{\bm{x}_{i}}=\Omega_{i}\bm{x}_{i}+\kappa_{2}\bm{X}_{\mathrm{av}}-\kappa_{2}\,\bm{x}_{i}\,(\bm{x}_{i}\centerdot\bm{X}_{\mathrm{av}}).$$

(5)

The system (2) has been intensively investigated for general $d$, usually with global coupling and for the homogeneous case $\Omega_{i}=0$, whether as a model of synchronization [20, 21, 22, 23, 24, 25, 26], or opinion formation and consensus studies on the unit sphere [27, 28, 29, 30, 31, 32, 33], or for the modelling of swarming behaviour [34, 35, 36]. Many synchronization properties have been established for any $d$ [21, 22, 37, 38, 39], in particular for the case of identical frequencies $\Omega_{i}=\Omega$, it is known that for $\kappa_{2}>0$ the order parameter $r=\|\bm{X}_{\mathrm{av}}\|$ evolves exponentially quickly to the value $r_{\infty}=1$, in which case all nodes are co-located to form a completely synchronized state, or for $\kappa_{2}<0$ to a state with $r_{\infty}=0$. Indeed, it is clear from (5) that for $\Omega_{i}=0$ either $\bm{x}_{i}=\bm{u}=\bm{X}_{\mathrm{av}}$ for some fixed unit vector $\bm{u}$, or $\bm{X}_{\mathrm{av}}=0$, each provide a steady state solution. For non-identical matrices $\Omega_{i}$ with a sufficiently large positive coupling $\kappa_{2}$, the system undergoes “practical synchronization” [22] as $\kappa_{2}\to\infty$, meaning that the configuration becomes asymptotically close to a completely synchronized state [38]. If $\kappa_{2}<0$, with distributed frequency matrices $\Omega_{i}$, the system remains asynchronous, a well-known property of the Kuramoto model. For specific values of $d$ more precise properties can be proved, for example for the Kuramoto model phase-locked synchronization occurs for any $\kappa_{2}>\kappa_{c}$ for some critical coupling $\kappa_{c}>0$ [40, 41, 42]. The $d=4$ case, for which trajectories lie on $\mathbb{S}^{3}$, equivalently on the $\mathrm{SU}_{2}$ group manifold, is similar to $d=2$, as numerical studies show [20].

The unit sphere models which follow from the potential (1) have only pairwise interactions, and so we now consider models with higher-order interactions. Let $\mathbb{N}_{N}=\{1,2,\dots N\}$, with $N\geqslant d$, then $(\bm{x}_{i_{1}},\bm{x}_{i_{2}},\dots,\bm{x}_{i_{d}})$ is a $d\times d$ matrix, where $\bm{x}_{i}$ is a unit column vector of length $d$, and $i_{1},i_{2},\dots i_{d}\in\mathbb{N}_{N}$. Define the potential

$$\mathfrak{V}_{d}=\sum_{i_{1},i_{2},\dots i_{d}=1}^{N}a_{[i_{1},i_{2},\dots i_{d}]}\,\det(\bm{x}_{i_{1}},\bm{x}_{i_{2}},\dots,\bm{x}_{i_{d}}),$$

(6)

where the summation is over all $i_{1},i_{2},\dots i_{d}\in\mathbb{N}_{N}$ and the coefficients $a_{[i_{1},i_{2},\dots i_{d}]}$ are antisymmetrized over these indices, corresponding to the antisymmetry of the determinant. For $d=3$, for example, we antisymmetrize any set of coefficients $a_{ijk}$ according to

$$a_{[i,j,k]}=\frac{1}{3!}(a_{ijk}+a_{kij}+a_{jki}-a_{jik}-a_{ikj}-a_{kji}).$$

$\mathfrak{V}_{d}$ as defined by (6) is rotationally invariant, i.e. if $\bm{x}_{i}\to R\bm{x}_{i}$ where $R\in\mathrm{SO}(d)$, then $\det(\bm{x}_{i_{1}},\bm{x}_{i_{2}},\dots,\bm{x}_{i_{d}})$ remains unchanged because $\det R=1$. If, however, $R\in\mathrm{O}(d)$ with $\det R=-1$, then $\mathfrak{V}_{d}$ changes sign, for example if $d$ is odd, then $\mathfrak{V}_{d}$ changes sign under parity inversion $\bm{x}_{i}\to-\bm{x}_{i}$. As a more general example, for any $d$, we can choose $R$ to be the identity matrix but with the sign of any one diagonal element reversed, then $R\in\mathrm{O}(d)$ with $\det R=-1$, and the sign of $\mathfrak{V}_{d}$ is reversed.

The potential (6) leads to a hierarchy of synchronization models on the unit sphere $\mathbb{S}^{d-1}$ with $d$-body couplings, since interactions take place between nodes only as constituents of oriented $(d-1)$-simplices. The order of the interaction is evidently tied to the dimension $d$ of the vector $\bm{x}_{i}$, and hence to the dimension of the unit sphere $\mathbb{S}^{d-1}$. For $d=2$ we write $\bm{x}_{i}=(\cos\theta_{i},\sin\theta_{i})$ and so $\mathfrak{V}_{2}=\sum_{i,j=1}^{N}a_{[ij]}\sin(\theta_{j}-\theta_{i})$, which we consider in Section 7, and for $d=3$ we obtain

$$\mathfrak{V}_{3}=\sum_{i,j,k=1}^{N}a_{[ijk]}\det(\bm{x}_{i},\bm{x}_{j},\bm{x}_{k})=\sum_{i,j,k=1}^{N}a_{[ijk]}\;\bm{x}_{i}\centerdot\bm{x}_{j}\times\bm{x}_{k},$$

(7)

to be analyzed in Section 3.

The coefficients $a_{[i_{1},i_{2},\dots i_{d}]}$ in (6) are antisymmetric but are otherwise unrestricted and so, in order to investigate the simplest possible case, and by analogy with the symmetric coupling $a_{ij}=1$ which we imposed for the pairwise interactions in the system (2), we choose the coefficients $a_{[i_{1},i_{2},\dots i_{d}]}$ to be the signature, denoted $\varepsilon_{i_{1}i_{2}\dots i_{d}}$, of the permutation $P=\{i_{1},i_{2},\dots i_{d}\}$. The signature is defined as the parity $(-1)^{n}$ of $P$, where $n$ is the number of transpositions of pairs of elements that must be composed to construct the permutation. This definition extends to indices which are not permutations by setting the value to be zero. We have, for example, $\varepsilon_{ij}=1$ if $j>i$, $\varepsilon_{ij}=-1$ if $j<i$, and $\varepsilon_{ij}=0$ if $j=i$ for any $i,j\in\mathbb{N}_{N}$. Similarly, $\varepsilon_{ijk}=1$ if $i<j<k$, or for any cyclic permutation of $i,j,k$, $\varepsilon_{ijk}=-1$ for any anticyclic permutation, and zero if any two indices are equal. In these higher-order models, therefore, each node interacts with all the other nodes in a connected $(d-1)$-simplex, with equal strength in magnitude, by means of the coupling coefficients $\varepsilon_{i_{1}i_{2}\dots i_{d}}$, with a sign depending on the orientation. There are $N!/(N-d)!$ nonzero coefficients $\varepsilon_{i_{1}i_{2}\dots i_{d}}$, and hence there are $N!/(N-d)!$ independent $d$-body couplings between the $N$ nodes, which take place through the $(d-1)$-simplices.

In order to calculate the gradient of $\mathfrak{V}_{d}$ as defined in (6), we write the determinant in terms of the Levi-Civita symbol $\varepsilon^{a_{1}a_{2}\dots a_{d}}$ and the components $x_{i}^{a}$ of $\bm{x}_{i}$, as follows:

$$\mathfrak{V}_{d}=\sum_{i_{1},i_{2},\dots i_{d}=1}^{N}\sum_{a_{1},a_{2},\dots a_{d}=1}^{d}\varepsilon_{i_{1}i_{2}\dots i_{d}}\;\varepsilon^{a_{1}a_{2}\dots a_{d}}\,x_{i_{1}}^{a_{1}}x_{i_{2}}^{a_{2}}\dots x_{i_{d}}^{a_{d}}.$$

(8)

For $d=3$, for example, we have $\det(\bm{x}_{i},\bm{x}_{j},\bm{x}_{k})=\bm{x}_{i}\centerdot\bm{x}_{j}\times\bm{x}_{k}=\sum_{a,b,c=1}^{3}\varepsilon^{abc}x_{i}^{a}x_{j}^{b}x_{k}^{c}$. For each set of indices $i_{2}\dots i_{d}$ define the vector $\bm{v}_{i_{2}\dots i_{d}}$ of length $d$ with components

$$(\bm{v}_{i_{2}\dots i_{d}})^{a}=\sum_{a_{2}\dots a_{d}=1}^{d}\varepsilon^{a\,a_{2}\dots a_{d}}\,x_{i_{2}}^{a_{2}}\dots x_{i_{d}}^{a_{d}},$$

(9)

then $\bm{v}_{i_{2}\dots i_{d}}$ behaves as a vector under rotations, i.e. if $\bm{x}_{i}\to R\bm{x}_{i}$ where $R\in\mathrm{SO}(d)$, then $\bm{v}_{i_{2}\dots i_{d}}\to R\bm{v}_{i_{2}\dots i_{d}}$. If $R\in\mathrm{O}(d)$ with $\det R=-1$, however, then $\bm{v}_{i_{2}\dots i_{d}}\to-R\bm{v}_{i_{2}\dots i_{d}}$. For $d=3$ we have $\bm{v}_{ij}=\bm{x}_{i}\times\bm{x}_{j}$. The identity $\bm{u}\centerdot\bm{v}_{i_{2}\dots i_{d}}=\det(\bm{u},\bm{x}_{i_{2}},\dots,\bm{x}_{i_{d}})$, for any fixed vector $\bm{u}$, provides a convenient method for calculating $\bm{v}_{i_{2}\dots i_{d}}$. In particular, we have $\bm{x}_{i}\centerdot\bm{v}_{i_{2}\dots i_{d}}=\det(\bm{x}_{i},\bm{x}_{i_{2}},\dots,\bm{x}_{i_{d}})$.

Vectors such as $\bm{v}_{i_{2}\dots i_{d}}$ can be viewed as the dual of elements belonging to an exterior algebra on which is defined a wedge or exterior product. The dual vector, generally referred to as the Hodge dual, has components formed using the Levi-Civita symbol as shown in (9). In three dimensions, for example, the vector product of two vectors can be viewed as the dual of the wedge product of these vectors. For our purposes, particularly for numerical evaluation or with a computer algebra system, it is sufficient to determine properties and explicit expressions for $\bm{v}_{i_{2}\dots i_{d}}$ either directly from (9) or by means of $\bm{u}\centerdot\bm{v}_{i_{2}\dots i_{d}}=\det(\bm{u},\bm{x}_{i_{2}},\dots,\bm{x}_{i_{d}})$, which holds for any fixed $d$-vector $\bm{u}$. For a detailed description of exterior algebras, properties of the antisymmetrization operator, and the definition of the Hodge dual, we refer to [43], Chapter 8.

The gradient of the potential $\mathfrak{V}_{d}$ defined by (8) is given for unconstrained vectors $\bm{x}_{i}$ by:

$$\nabla_{i}\mathfrak{V}_{d}=d\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\;\bm{v}_{i_{2}\dots i_{d}},$$

and so after including Lagrange multipliers in order to enforce the constraints, the equations of motion are:

$$\dot{\bm{x}_{i}}=\frac{\kappa_{d}}{N^{d-1}}\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\left[\bm{v}_{i_{2}\dots i_{d}}-\bm{x}_{i}\,(\bm{x}_{i}\centerdot\bm{v}_{i_{2}\dots i_{d}})\right],$$

(10)

where we have included a normalized coupling constant $\kappa_{d}$. For fixed $i$ there are $(N-1)!/(N-d)!$ nonzero terms under the summation, which for large $N$ approaches $N^{d-1}$, and so we have normalized $\kappa_{d}$ by this factor. We can add oscillator terms $\Omega_{i}\bm{x}_{i}$ to the right-hand side, as for the system (2), but if the matrices $\Omega_{i}$ are identical, $\Omega_{i}=\Omega$, then we can replace $\bm{x}_{i}\to\rme^{\Omega t}\bm{x}_{i}$ in the usual way and so because $\rme^{\Omega t}$ is a rotation matrix, we can in effect set $\Omega=0$. In this case the coefficient $\kappa_{d}/N^{d-1}$ in (10) can be set to $\pm 1$ by rescaling the time variable. For every solution of (10) with $\kappa_{d}>0$ there corresponds another solution with $\kappa_{d}<0$, obtained by transforming $\bm{x}_{i}\to R\bm{x}_{i}$, where $R\in\mathrm{O}(d)$ with $\det R=-1$. If $d$ is odd, for example, parity inversion $\bm{x}_{i}\to-\bm{x}_{i}$ in effect changes the sign of $\kappa_{d}$. The properties of the system (10) are therefore independent of the sign of $\kappa_{d}$, in contrast to the Kuramoto model, or more generally the $2$-body system (2). The form of (10) is the same as shown in (3), except that now each vector $\bm{x}_{i}$ couples to $\bm{X}_{i}=\kappa_{d}\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\bm{v}_{i_{2}\dots i_{d}}/N^{d-1}$, which necessarily depends on $i$, and so the coupling takes place through all possible $(d-1)$-simplices containing the $i^{\mathrm{th}}$ node.

We wish to determine the behaviour of the system (10) for generic initial values $\bm{x}_{i}(0)=\bm{x}_{i}^{0}$, where “generic” means we exclude unstable fixed points, and look for synchronized final configurations. A useful measure of synchronization, whether for $2$-body or $d$-body interactions, is by means of the average position $\bm{X}_{\mathrm{av}}$ defined in (4), from which we calculate the rotationally invariant order parameter $r=\|\bm{X}_{\mathrm{av}}\|$. Synchronization occurs when $r$ achieves a constant asymptotic value $r_{\infty}$, where $0\leqslant r_{\infty}\leqslant 1$. If $r_{\infty}=1$ all particles are co-located at a common point, usually referred to as a completely synchronized configuration, and if $r_{\infty}=0$ then $\bm{X}_{\mathrm{av}}=0$ as $t\to\infty$, hence the particles are distributed over the sphere with an average position of zero, sometimes referred to as a balanced configuration. Both cases occur for the 2-body system (2) with identical frequencies, depending on the sign of the coupling [21].

For the values of $d$ under consideration we find numerically that the system (10) always attains a final static configuration for all generic initial values, and so we wish to determine all stable steady state solutions. This appears to be a formidable task, considering that the nonlinearities on the right-hand side of (10) are of order $d+1$, however we show that the static equations can be solved by a simpler reduced system. Due to the rotational covariance of (10), any static configuration can be rotated to an arbitrary orientation, and so for numerical comparison of final configurations, which always depend on the initial values, we calculate rotational invariants such as $\bm{x}_{i}\centerdot\bm{x}_{j}$.

We observe firstly that (10) has the fixed point solution $\bm{x}_{i}=\pm\bm{u}$, where $\bm{u}$ is any constant unit vector, with a sign that can depend on $i$, because by antisymmetry $\bm{v}_{i_{2}\dots i_{d}}$ is zero if the vectors $\bm{x}_{i}$ are either parallel or antiparallel. There is the possibility therefore that the system completely synchronizes, as occurs for the $2$-body system (2), however numerically we find that such configurations are unstable. Instead, we look for static solutions satisfying

$$\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\bm{v}_{i_{2}\dots i_{d}}=\lambda\,\bm{x}_{i},$$

(11)

where $\lambda$ is independent of $i$, and therefore takes the value

$$\lambda=\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\bm{x}_{i}\centerdot\bm{v}_{i_{2}\dots i_{d}}=\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\det(\bm{x}_{i},\bm{x}_{i_{2}},\dots,\bm{x}_{i_{d}}).$$

It follows from (11) that the right-hand side of (10) is zero, indeed this is true even if $\lambda$ depends on $i$, however we find numerically that static solutions satisfy the special form (11) with $\lambda$ independent of $i$, and we show by direct calculation that this holds exactly for $d=3$ see (20), for $d=4$ see (35) and for $d=5$ see (49).

We are also interested in combining $2$-body and $d$-body systems in order to evaluate the relative effect of the corresponding couplings, hence we combine (5) and (10) to obtain

$$\dot{\bm{x}}_{i}=\kappa_{2}\bm{X}_{\mathrm{av}}-\kappa_{2}\,\bm{x}_{i}(\bm{x}_{i}\centerdot\bm{X}_{\mathrm{av}})+\frac{\kappa_{d}}{N^{d-1}}\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\left[\bm{v}_{i_{2}\dots i_{d}}-\bm{x}_{i}\,(\bm{x}_{i}\centerdot\bm{v}_{i_{2}\dots i_{d}})\right],$$

(12)

where $\kappa_{2}$ denotes the strength of the $2$-body forces. We can satisfy the static equations in this case by solving

$$\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\bm{v}_{i_{2}\dots i_{d}}=\lambda_{1}\bm{x}_{i}-\lambda_{2}\bm{X}_{\mathrm{av}},$$

(13)

for constants $\lambda_{1},\lambda_{2}$ independent of $i$, since then we have

$$\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\bm{x}_{i}\centerdot\bm{v}_{i_{2}\dots i_{d}}=\lambda_{1}-\lambda_{2}\,(\bm{x}_{i}\centerdot\bm{X}_{\mathrm{av}}),$$

and hence the right-hand side of (12) is zero, provided that

$$\lambda_{2}=\frac{\kappa_{2}N^{d-1}}{\kappa_{d}}.$$

(14)

More directly, define

$$\bm{w}_{i}=\frac{\kappa_{d}}{N^{d-1}}\sum_{i_{2},\dots i_{d}=1}^{N}\varepsilon_{i,i_{2}\dots i_{d}}\bm{v}_{i_{2}\dots i_{d}}-\frac{\lambda_{1}\kappa_{d}}{N^{d-1}}\bm{x}_{i}+\kappa_{2}\bm{X}_{\mathrm{av}},$$

(15)

then we can rewrite (12) equivalently as $\dot{\bm{x}}_{i}=\bm{w}_{i}-\bm{x}_{i}\,(\bm{x}_{i}\centerdot\bm{w}_{i})$, and so we see immediately that $\dot{\bm{x}}_{i}=0$ for $\bm{w}_{i}=0$. We note here that the signs of $\lambda_{1},\lambda_{2}$ in (13) are each reversed under the transformation $\bm{x}_{i}\to R\bm{x}_{i}$, where $R\in\mathrm{O}(d)$ with $\det R=-1$, which is equivalent to reversing the sign of $\kappa_{d}$, and so the signs of $\lambda_{1},\lambda_{2}$ are correlated with the sign of $\kappa_{d}$.

It is not at all evident that stable static solutions should satisfy (13) for all $i$, but numerically we find this to be true in all cases, and prove for $d=3,5$ that it holds exactly. For static solutions that depend on one or more unknown parameters, we can in principle determine $\lambda_{2}$ as a function of these parameters, then (14) expresses these parameters explicitly in terms of the ratio $\kappa_{2}/\kappa_{d}$. Let us proceed now to the case $d=3$, for which we can find exact steady state solutions and so verify (13).

For any static synchronized configuration as shown in Fig. 1(a) we can point the normal $\bm{n}$ along the $Z$-axis, i.e. we can set $\bm{n}=(0,0,1)$ by means of an $\mathrm{SO}(3)$ rotation. Such rotations can be performed algebraically or numerically as explained in D. Then, since all nodes lie in the plane $z=\frac{1}{\sqrt{3}}$, the steady state solution is given by:

$$\bm{x}_{i}=\frac{1}{\sqrt{3}}\left(\sqrt{2}\cos\frac{2i\pi}{N},\sqrt{2}\sin\frac{2i\pi}{N},1\right),$$

(17)

for $i=1,\dots N$, where we have also performed an $\mathrm{SO}(2)$ rotation about the $Z$-axis so that $\bm{x}_{N}=\frac{1}{\sqrt{3}}\left(\sqrt{2},0,1\right)$ is aligned along the $X$-axis. The sign of $\bm{x}_{i}$ varies according to the sign of $\kappa_{3}$, i.e. either (17) or its negative is a stable fixed point.

From (17) we can compute the following rotational invariants, where we denote $x_{ij}=\bm{x}_{i}\centerdot\bm{x}_{j}$ and $x_{ijk}=\bm{x}_{i}\centerdot\bm{x}_{j}\times\bm{x}_{k}$:

	$\displaystyle x_{ij}$	$\displaystyle=$	$\displaystyle\frac{1}{3}\left[1+2\cos\frac{2(i-j)\pi}{N}\right],$		(18)
	$\displaystyle x_{ijk}$	$\displaystyle=$	$\displaystyle\frac{8}{3\sqrt{3}}\sin\frac{(i-j)\pi}{N}\sin\frac{(k-i)\pi}{N}\sin\frac{(j-k)\pi}{N}.$		(19)

These invariants are related by the identity $x_{ijk}^{2}=1-x_{ij}^{2}-x_{ik}^{2}-x_{jk}^{2}+2x_{ij}x_{ik}x_{jk}$ as we discuss for $N=3$ in Section 3.3. Since these expressions are independent of the orientation of the system, they provide a convenient means for the numerical verification of all exact expressions, for any initial values, and hence for all final configurations.

We see from (18) that $x_{ij}$ depends only on the difference $|i-j|$, one consequence of which is that adjoining nodes are equally spaced, since $\|\bm{x}_{i}-\bm{x}_{i+1}\|$ is independent of $i$. In addition, this distance is equal to $\|\bm{x}_{1}-\bm{x}_{N}\|$, showing that the nodes form a closed loop, as follows from the symmetry $\bm{x}_{i}\centerdot\bm{x}_{N}=\bm{x}_{N-i}\centerdot\bm{x}_{N}$ for $i=1,2\dots N-1$. Although this property is evident from Fig. 1(a), it does not hold for even values of $d$, see Section 5 for $d=4$. It is convenient to extend the formula (17) to the value $i=0$, then we have $\bm{x}_{0}=\bm{x}_{N}$ which again shows that the nodes form a closed sequence. Configurations with equally spaced nodes are sometimes referred to as splay states, although for $d>3$ the asymptotic configurations are not restricted to a plane as occurs here; for $d=4$, for example, we obtain equally spaced nodes lying on a torus embedded in $\mathbb{S}^{3}$, see Section 5. If we define a unit $2$-vector $\bm{u}_{i}$ from (17) according to $\bm{x}_{i}=\frac{1}{\sqrt{3}}(\sqrt{2}\bm{u}_{i},1)$ then $\sum_{i}\bm{u}_{i}=0$ and so $\bm{u}_{i}$ can be regarded as a splay state as described in [44].

Directly from (17) and the definition (4), together with well-known trigonometric sums derived in A, see in particular (65), we obtain $\bm{X}_{\mathrm{av}}=\frac{1}{\sqrt{3}}(0,0,1)$. The vectors (17) also satisfy, as proved in B:

$$\frac{1}{N^{2}}\sum_{j,k=1}^{N}\varepsilon_{ijk}(\bm{x}_{j}\times\bm{x}_{k})=\frac{2}{N\sqrt{3}}\,\cot\frac{\pi}{N}\;\bm{x}_{i},$$

(20)

for all $i$, which verifies the relation (11) for $d=3$, with $\lambda=\frac{2N}{\sqrt{3}}\cot\frac{\pi}{N}$. From (20) we obtain

$$\frac{1}{N^{2}}\sum_{j,k=1}^{N}\varepsilon_{ijk}(\bm{x}_{i}\centerdot\bm{x}_{j}\times\bm{x}_{k})=\frac{2}{N\sqrt{3}}\,\cot\frac{\pi}{N},$$

from which $\sum_{j,k=1}^{N}\varepsilon_{ijk}\,\left[\bm{x}_{j}\times\bm{x}_{k}-\bm{x}_{i}\,(\bm{x}_{i}\centerdot\bm{x}_{j}\times\bm{x}_{k})\right]=0$. Hence (17) is an exact static solution of (16), which numerically we find to be stable for $\kappa_{3}>0$, while for $\kappa_{3}<0$ its negative is stable.

Consider next the system (16) for distributed frequency vectors $\bm{\omega}_{i}$. It is well-known for the $2$-body system (2) that in this case phase-locked synchronization, in the sense that $r(t)=\|\bm{X}_{\mathrm{av}}(t)\|$ approaches a constant value, does not occur exactly for $d=3$. Rather, $r$ varies between narrow limits which decrease as $\kappa_{2}\to\infty$, which has been termed practical synchronization [21, 22], i.e. the particle positions do not shrink to a single point, but are confined to a small region with a diameter inversely proportional to $\kappa_{2}$ [21]. This phenomenon of practical synchronization also appears for (16), since we find numerically that $r(t)$ again varies asymptotically between narrow limits which decrease as $|\kappa_{3}|$ increases. Ring synchronization, however, is still clearly evident, as in the example in Fig. 1(b), which shows an asymptotic configuration for $N=40$ for which the average value for $r$ is close to $\frac{1}{\sqrt{3}}$. The frequency vectors $\bm{\omega}_{i}$ in this example have random entries with approximate unit length, and $\kappa_{3}=20$. Following the initial transient, the configuration rotates on $\mathbb{S}^{2}$ with small variations in the relative positions of the nodes. For small values of $|\kappa_{3}|$ the particle motion is asynchronous, suggesting that there is a critical value $\kappa_{c}$ such that practical synchronization occurs only for $|\kappa_{3}|>\kappa_{c}>0$.

The special case $N=3$ of (16) may be solved exactly for the rotational invariants in terms of elliptic functions, and provides insight into the rate of synchronization and the stability of the synchronized configurations. In particular, we show why the fixed point corresponding to complete synchronization is unstable, whereas the ring synchronized configuration, corresponding to the $N=3$ case of (18), is stable. With $\bm{\omega}_{i}=0$ and $\kappa_{3}/N^{2}=1$, (16) reduces to:

	$\displaystyle\dot{\bm{x}_{1}}$	$\displaystyle=$	$\displaystyle 2\left[\bm{x}_{2}\times\bm{x}_{3}-\bm{x}_{1}\;(\bm{x}_{1}\centerdot\bm{x}_{2}\times\bm{x}_{3})\right]$		(21)
	$\displaystyle\dot{\bm{x}_{2}}$	$\displaystyle=$	$\displaystyle 2\left[\bm{x}_{3}\times\bm{x}_{1}-\bm{x}_{2}\;(\bm{x}_{2}\centerdot\bm{x}_{3}\times\bm{x}_{1})\right]$
	$\displaystyle\dot{\bm{x}_{3}}$	$\displaystyle=$	$\displaystyle 2\left[\bm{x}_{1}\times\bm{x}_{2}-\bm{x}_{3}\;(\bm{x}_{3}\centerdot\bm{x}_{1}\times\bm{x}_{2})\right],$

hence

$$\frac{d}{dt}(\bm{x}_{1}\centerdot\bm{x}_{2})=\bm{x}_{1}\centerdot\dot{\bm{x}_{2}}+\dot{\bm{x}_{1}}\centerdot\bm{x}_{2}=-4\bm{x}_{1}\centerdot\bm{x}_{2}\;(\bm{x}_{1}\centerdot\bm{x}_{2}\times\bm{x}_{3}),$$

together with the cyclic permutations. Denote $x_{ij}=\bm{x}_{i}\centerdot\bm{x}_{j},x_{ijk}=\bm{x}_{i}\centerdot\bm{x}_{j}\times\bm{x}_{k}$, then we have $\dot{x}_{12}=-4x_{12}x_{123},\dot{x}_{23}=-4x_{23}x_{123}$ and $\dot{x}_{13}=-4x_{13}x_{123}$. Hence the ratios $x_{23}/x_{12}$ and $x_{13}/x_{12}$ are constants of the motion. Let us choose $u=x_{12}$ to be the independent variable, then $x_{23}=c_{1}u,x_{13}=c_{2}u$ for constants $c_{1},c_{2}$ which are fixed by the initial values $\bm{x}_{i}^{0}$, i.e. we have $c_{1}=x_{23}^{0}/x_{12}^{0},c_{2}=x_{13}^{0}/x_{12}^{0}$, assuming that $u_{0}=x_{12}^{0}=\bm{x}_{1}^{0}\centerdot\bm{x}_{2}^{0}$ is nonzero. The constants $c_{1},c_{2}$ can take any value, positive or negative.

Consider now the fixed points of the system (21). These satisfy $\bm{x}_{1}\times\bm{x}_{2}=\bm{x}_{3}\,x_{123}$ and its cyclic permutations, which implies $x_{12}x_{123}=x_{13}x_{123}=x_{23}x_{123}=0$. Either $x_{123}=0$, in which case all cross products are zero, $\bm{x}_{i}\times\bm{x}_{j}=0$, or $x_{123}\neq 0$ which implies that the three vectors $\bm{x}_{i}/x_{123}$ form an orthonormal set. In the first case, $\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}$ are either parallel or antiparallel, with four possible choices of relative sign, and so the three nodes are located at either a single point, or at two opposite points on the unit circle. In either case we have $|c_{1}|=|c_{2}|=1$, and since $c_{1},c_{2}$ are constants of the motion, these fixed points exist for the system (21) only if the initial values $\bm{x}_{i}^{0}$ are consistent with this, i.e. only if $|x_{12}^{0}|=|x_{13}^{0}|=|x_{23}^{0}|=1$. These fixed points are therefore unstable, since under any perturbation of the initial values these steady state solutions cannot be attained as the system evolves.

For the second case, in which $x_{123}\neq 0$, we must have $|x_{123}|=1$ and so $\{\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}\}$ forms an orthonormal set, with an arbitrary orientation with respect to the axes, satisfying $x_{12}=x_{13}=x_{23}=0$. Let us show that (21) synchronizes to these points from all initial values, except for the unstable fixed points. Firstly, we express $x_{123}^{2}$ in terms of $x_{12},x_{23},x_{13}$ by means of a well-known identity obtained as follows: define the $3\times 3$ matrix $M=(\bm{x}_{1},\bm{x}_{2},\bm{x}_{3})$ then $(M^{\mathsf{T}}M)_{ij}=x_{ij}$, and so we obtain $(\det M)^{2}=\det M^{\mathsf{T}}\det M=\det(x_{ij})$, which is the Gram determinant, hence $x_{123}^{2}=1-x_{12}^{2}-x_{13}^{2}-x_{23}^{2}+2x_{12}x_{13}x_{23}$. Define the cubic polynomial

$$p(u)=1-(1+c_{1}^{2}+c_{2}^{2})u^{2}+2c_{1}c_{2}u^{3},$$

(22)

then $x_{123}^{2}=p(u)$. From $\dot{u}=-4u\,x_{123}$ we obtain $\dot{u}^{2}=16u^{2}x_{123}^{2}=16u^{2}p(u)=2V(u)$, where we have defined the potential $V(u)=8u^{2}p(u)$ as a fifth order polynomial in $u$. The equation $\dot{u}^{2}=2V(u)$ can be solved for $u$ with the initial value $u(0)=u_{0}$, however, in order to avoid taking the square root with an ambiguous sign, it is preferable from a numerical perspective to solve the equivalent second order equation $\ddot{u}=V^{\prime}(u)$ with $u(0)=u_{0},\dot{u}(0)=-4u_{0}x_{123}^{0}$.