Higher-order synchronization on the sphere

We construct a system of $N$ interacting particles on the unit sphere $S^{d-1}$ in $d$-dimensional space, which has $d$-body interactions only. The equations have a gradient formulation derived from a rotationally-invariant potential of a determinantal form summed over all nodes, with antisymmetric coefficients. For $d=3$, for example, all trajectories lie on the 2-sphere and the potential is constructed from the triple scalar product summed over all oriented 2-simplices. We investigate the cases $d=3,4,5$ in detail, and find that the system synchronizes from generic initial values, for both positive and negative coupling coefficients, to a static final configuration in which the particles lie equally spaced on $S^{d-1}$. Completely synchronized configurations also exist, but are unstable under the $d$-body interactions. We compare the relative effect of 2-body and $d$-body forces by adding the well-studied 2-body interactions to the potential, and find that higher-order interactions enhance the synchronization of the system, specifically, synchronization to a final configuration consisting of equally spaced particles occurs for all $d$-body and 2-body coupling constants of any sign, unless the attractive 2-body forces are sufficiently strong relative to the $d$-body forces. In this case the system completely synchronizes as the 2-body coupling constant increases through a positive critical value, with either a continuous transition for $d=3$, or discontinuously for $d=5$. Synchronization also occurs if the nodes have distributed natural frequencies of oscillation, provided that the frequencies are not too large in amplitude, even in the presence of repulsive 2-body interactions which by themselves would result in asynchronous behaviour.


Introduction
Emerging phenomena in complex systems have generally been modelled by means of networks with interacting pairs of nodes, for which the prevailing paradigm is the Kuramoto model [1] and its many generalizations. There has been an increasing realization, however, that pairwise interactions and associated network structures are not adequate to describe the group interactions which occur in, and possibly dominate, the behaviour of many complex systems of interest. Higher-order interactions are known to occur in neuroscience, ecology, social systems and many others [2,3]. We refer to the review [4] for a detailed discussion of higher-order interactions and their applications with an extensive list of references, including the various frameworks which have been used to model higher-order systems, as well as the recent overview [5] of higher-order networks.
We describe here a hierarchy of models with higher-order interactions which are similar to the well-known Kuramoto systems, and focus firstly on the properties of exclusive d-body interactions, where in our formulation d is also the dimension of the system. Kuramoto models on the unit sphere with conventional 2-body interactions have been well-studied and generalize the simple Kuramoto model which corresponds to the case d = 2. In this approach we associate to each of the N nodes a unit vector x i of length d, which lies therefore on S d−1 , and evolves according to first-order equations in the gradient formẋ i = ∇ i V d , where the bilinear potential V d = V d (x 1 , . . . , x N ) is invariant under rotations in SO(d). Because the potential is bilinear, only 2-body interactions occur.
We can model higher-order interactions, however, by choosing instead a multilinear determinantal form for V d which, by antisymmetry, allows only d-body interactions. For d = 3, for example, for which all trajectories x i (t) lie on the 2-sphere, V 3 is constructed from the triple scalar product which, being antisymmetric, allows only 3-body interactions, i.e. each node i interacts with the other nodes only through the oriented 2-simplices containing the i th node. We find that synchronization occurs for all values of the coupling constant, whether positive or negative, with the system evolving from all initial values to the same final static configuration, up to an arbitrary rotation. The equations for larger values of d are of formidable complexity due to the nonlinearities of degree d + 1, and so we consider only d 5, although some general results may be derived. We restrict our investigation to the simplest possible models with d-body interactions by choosing the coupling coefficients to be signature functions, which are antisymmetric and of magnitude either zero or one, but we also briefly describe the effect of oscillations by including distributed frequency matrices.
Besides the modelling of exclusive d-body forces, we can also include pairwise interactions in any dimension d by adding the usual bilinear terms to the potential, and so we are able to investigate the relative effect of combined 2-body and d-body forces, in principle for any d, but here in detail only for d 5. For d = 3, for example, we find that the 3-body forces enhance synchronization, and so if the coupling constant for the 2-body interactions is negative, resulting in repulsive 2-body forces that would otherwise prevent the system from synchronizing, the combined system nevertheless synchronizes due to the 3-body interactions. This is consistent with previous observations that "higher-order interactions can stabilize strongly synchronized states even when the pairwise coupling is repulsive" [3,6].
For a general description of synchronization with respect to higher-order interactions we refer to [4] (Section 6) also [5] (Section 4) and [7], and note that extensions of the Kuramoto model to higher-order networks have been extensively investigated, but generally with trajectories restricted to the unit circle S 1 [8,9,10,11,12,13]. For a description and properties of simplicial complexes, and the relation to algebraic topology we refer to [4,12,14], however our focus here is on the properties of specific dynamical systems defined on the unit sphere which allow higherorder interactions, and their synchronization behaviour. The effect of 2-body forces in combination with 3-and 4-body forces has previously been considered, but only in the context of phase oscillators on S 1 [6,8]. Properties such as multistability, which refers to the coexistence of multiple steady states, has been observed in these and other models [8,9,14], and it is generally thought "that higher-order interactions favor multistability" [4], but it remains an open question as to whether such properties are model-specific, or are general consequences of higher-order interactions. Multistability, for example, does not appear in the models that we consider, although to any stable steady state solution there corresponds, by rotational covariance, a family of rotated steady states. A similar observation can be made with respect to "explosive synchronization" which has been observed in a higher-order Kuramoto model [13], although similar properties are wellknown to occur in models with only 2-body interactions [15,16,17]. This behaviour is also absent in the higher-order models that we consider, similarly for the property of "abrupt transitions" as appears in various models [4,6,9]. We do find, however, that there is a discontinuous transition for d = 5 at which the 2-body forces suddenly become sufficiently strong to overcome the 5-body forces, as discussed in Section 6. One difference in our approach, which possibly eliminates discontinuous behaviours, is that our equations take a gradient form which can be used to deduce the stability of steady states of the model, as explained in [14].

Outline and summary of results
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 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 κ 2 , κ 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 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 κ 2 /κ 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 S 2 . These states exist and are stable, unless the ratio of coupling constants κ 2 /κ 3 becomes sufficiently large, at which point the system transitions continuously to a completely synchronized formation. We can state precisely the value of κ 2 /κ 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 κ 2 /κ 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 κ 2 /κ 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 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.

Synchronization models on the unit sphere
Synchronization models in which the N trajectories of the complex system are confined to the unit sphere have been extensively investigated, particularly for the special case of the Kuramoto model. In general, a unit d-vector x i is associated with the i th node of the network and evolves according to the nonlinear equations obtained by minimizing a potential V d (x 1 , . . . , x N ), where the constraint x i x i = 1 for every i = 1, . . . N is enforced by means of Lagrange multipliers.
The first-order evolution equations have the gradient formẋ i = ∇ i V d , to which one can add further terms such as Ω i x i , where Ω i is a d × d antisymmetric frequency matrix, in which case each node has one or more natural modes of oscillation as determined by the eigenvalues of Ω i . Properties of the model depend on the choice of V d , which is invariant under transformations of the rotation group SO(d). There are two possible such choices for V d , either the bilinear inner product, or the triple scalar product for d = 3 and its generalization to any d. This leads to two types of models, those with pairwise interactions, and those with d-body interactions. We can also additively combine these two possibilities and hence investigate the relative effect of d-body and 2-body forces within the system.

Pairwise interactions on the unit sphere
The choice for V d which has been widely investigated is where the coefficients a ij are symmetric, and in this case the evolution equations areẋ i = −λ i x i + j a ij x j where λ i are Lagrange multipliers, supplemented by the constraint x i x i = 1. Hence λ i = j a ij x i x j and so we obtain the system of N equationsẋ where we have included d × d antisymmetric frequency matrices Ω i , as well as a normalized coupling coefficient κ 2 which measures the relative strength of the natural oscillations compared to the nonlinear interactions. Since x i ẋ i = 0, the unit length of x i is maintained as the system evolves. The Kuramoto model corresponds to the case d = 2 with the parametrization x i = (cos θ i , sin θ i ) and Ω i = 0 −ω i ω i 0 , for which (2) reduces toθ i = ω i + κ 2 N N j=1 a ij sin(θ j − θ i ). In this case the potential (1) reduces to V 2 = i,j a ij cos(θ j − θ i ) which is well-known as a Lyapunov function for the Kuramoto model [18,19]. It is convenient to write (2) in the forṁ where X i = κ 2 N j=1 a ij 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 X av defined by then we can write (2) as: The system (2) has been intensively investigated for general d, usually with global coupling and for the homogeneous case Ω 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 Ω i = Ω, it is known that for κ 2 > 0 the order parameter r = X av evolves exponentially quickly to the value r ∞ = 1, in which case all nodes are co-located to form a completely synchronized state, or for κ 2 < 0 to a state with r ∞ = 0. Indeed, it is clear from (5) that for Ω i = 0 either x i = u = X av for some fixed unit vector u, or X av = 0, each provide a steady state solution. For nonidentical matrices Ω i with a sufficiently large positive coupling κ 2 , the system undergoes "practical synchronization" [22] as κ 2 → ∞, meaning that the configuration becomes asymptotically close to a completely synchronized state [38]. If κ 2 < 0, with distributed frequency matrices Ω 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 κ 2 > κ c for some critical coupling κ c > 0 [40,41,42]. The d = 4 case, for which trajectories lie on S 3 , equivalently on the SU 2 group manifold, is similar to d = 2, as numerical studies show [20].

Higher-order interactions on the unit sphere
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 where x i is a unit column vector of length d, and i 1 , i 2 , . . . i d ∈ N N . Define the potential where the summation is over all i 1 , i 2 , . . . i d ∈ N N and the coefficients a [i 1 ,i 2 ,...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 . . , x i d ) remains unchanged because det R = 1. If, however, R ∈ O(d) with det R = −1, then V d changes sign, for example if d is odd, then V d changes sign under parity inversion x i → −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 ∈ O(d) with det R = −1, and the sign of V d is reversed. The potential (6) leads to a hierarchy of synchronization models on the unit sphere 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 x i , and hence to the dimension of the unit sphere S d−1 . For d = 2 we write x i = (cos θ i , sin θ i ) and so V 2 = N i,j=1 a [ij] sin(θ j − θ i ), which we consider in Section 7, and for d = 3 we obtain to be analyzed in Section 3. The coefficients a [i 1 ,i 2 ,...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 ,...i d ] to be the signature, denoted ε i 1 i 2 ...i d , of the permutation P = {i 1 , i 2 , . . . 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, ε ij = 1 if j > i, ε ij = −1 if j < i, and ε ij = 0 if j = i for any i, j ∈ N N . Similarly, ε ijk = 1 if i < j < k, or for any cyclic permutation of i, j, k, ε 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 ε i 1 i 2 ...i d , with a sign depending on the orientation. There are N!/(N − d)! nonzero coefficients ε i 1 i 2 ...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.

Nonlinear d-body equations on the sphere
In order to calculate the gradient of V d as defined in (6), we write the determinant in terms of the Levi-Civita symbol ε a 1 a 2 ...a d and the components x a i of x i , as follows: For d = 3, for example, we have det( then v i 2 ...i d behaves as a vector under rotations, i.e. if , for any fixed vector u, provides a convenient method for calculating v i 2 ...i d . In particular, we have Vectors such as v i 2 ...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 v i 2 ...i d either directly from (9) or by means of u v i 2 ...i d = det(u, x i 2 , . . . , x i d ), which holds for any fixed d-vector 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 V d defined by (8) is given for unconstrained vectors x i by: and so after including Lagrange multipliers in order to enforce the constraints, the equations of motion are: where we have included a normalized coupling constant κ 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 κ d by this factor. We can add oscillator terms Ω i x i to the right-hand side, as for the system (2), but if the matrices Ω i are identical, Ω i = Ω, then we can replace x i → e Ωt x i in the usual way and so because e Ωt is a rotation matrix, we can in effect set Ω = 0. In this case the coefficient κ d /N d−1 in (10) can be set to ±1 by rescaling the time variable. For every solution of (10) with κ d > 0 there corresponds another solution with κ d < 0, obtained by transforming If d is odd, for example, parity inversion x i → −x i in effect changes the sign of κ d . The properties of the system (10) are therefore independent of the sign of κ 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 which necessarily depends on i, and so the coupling takes place through all possible (d − 1)-simplices containing the i th node.
We wish to determine the behaviour of the system (10) for generic initial values x i (0) = x 0 i , 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 X av defined in (4), from which we calculate the rotationally invariant order parameter r = X av . Synchronization occurs when r achieves a constant asymptotic value r ∞ , where 0 r ∞ 1. If r ∞ = 1 all particles are co-located at a common point, usually referred to as a completely synchronized configuration, and if r ∞ = 0 then X av = 0 as t → ∞, 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].

Steady state solutions
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 x i x j .
We observe firstly that (10) has the fixed point solution x i = ±u, where u is any constant unit vector, with a sign that can depend on i, because by antisymmetry v i 2 ...i d is zero if the vectors 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 where λ is independent of i, and therefore takes the value It follows from (11) that the right-hand side of (10) is zero, indeed this is true even if λ depends on i, however we find numerically that static solutions satisfy the special form (11) with λ 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 where κ 2 denotes the strength of the 2-body forces. We can satisfy the static equations in this case by solving for constants λ 1 , λ 2 independent of i, since then we have and hence the right-hand side of (12) is zero, provided that More directly, define then we can rewrite (12) equivalently asẋ i = w i −x i (x i w i ), and so we see immediately thatẋ i = 0 for w i = 0. We note here that the signs of λ 1 , λ 2 in (13) are each reversed under the transformation which is equivalent to reversing the sign of κ d , and so the signs of λ 1 , λ 2 are correlated with the sign of κ 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 λ 2 as a function of these parameters, then (14) expresses these parameters explicitly in terms of the ratio κ 2 /κ d . Let us proceed now to the case d = 3, for which we can find exact steady state solutions and so verify (13).

3-body interactions on the 2-sphere
For d = 3 the N equations (10) for the unit 3-vectors x i ∈ S 2 are: where we have included frequency 3-vectors ω i = (ω 1 i , ω 2 i , ω 3 i ). For the homogeneous case ω i = 0, the right-hand side of (16) has a gradient form derived from the potential V 3 = N i,j,k=1 ε ijk x i x j × x k , and the coupling constant κ 3 /N 2 can be scaled to unity, with a sign that can be reversed by means of the parity transformation x i → −x i for all i. The system has the fixed point x i = u, where u is any constant unit vector but, numerically, we find that such completely synchronized configurations are unstable. Instead, the system synchronizes from generic initial values to a static configuration, with an orientation that depends on the initial values, of the form shown in Fig. 1(a). We refer to this as ring synchronization, since all nodes lie equally spaced on a circle arising from the intersection of a plane with the unit sphere. Fig. 1(a) shows the configuration for N = 40 nodes (red), and the unit normal n (green), where n = X av / X av with X av defined in (4). The asymptotic configuration in all cases satisfies n We can write down an explicit expression for ring-synchronized configurations such as in Fig. 1(a), up to an arbitrary SO(3) rotation, and show that these are exact fixed point solutions of (16). The stability of these configurations is a numerical observation, however, although we prove this to be the case for N = 3 in Section 3.3 by solving exactly for the rotational invariants.

Steady state solutions on S 2
For any static synchronized configuration as shown in Fig. 1(a) we can point the normal n along the Z-axis, i.e. we can set n = (0, 0, 1) by means of an SO(3) rotation. Such rotations can be performed algebraically or numerically as explained in Appendix D. Then, since all nodes lie in the plane z = 1 √ 3 , the steady state solution is given by: for i = 1, . . . N, where we have also performed an SO(2) rotation about the Z-axis so that x N = 1 √ 3 √ 2, 0, 1 is aligned along the X-axis. The sign of x i varies according to the sign of κ 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 These invariants are related by the identity 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 x i − x i+1 is independent of i. In addition, this distance is equal to x 1 − x N , showing that the nodes form a closed loop, as follows from the symmetry 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 x 0 = 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 S 3 , see Section 5. If we define a unit 2-vector u i from (17) according to ) then i u i = 0 and so 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 Appendix A, see in particular (A.4), we obtain X av = 1 √ 3 (0, 0, 1). The vectors (17) also satisfy, as proved in Appendix B: for all i, which verifies the relation (11) (17) is an exact static solution of (16), which numerically we find to be stable for κ 3 > 0, while for κ 3 < 0 its negative is stable.

3-body systems with distributed frequencies
Consider next the system (16) for distributed frequency vectors ω i . It is well-known for the 2-body system (2) that in this case phase-locked synchronization, in the sense that r(t) = X av (t) approaches a constant value, does not occur exactly for d = 3. Rather, r varies between narrow limits which decrease as κ 2 → ∞, 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 κ 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 |κ 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 1 √ 3 . The frequency vectors ω i in this example have random entries with approximate unit length, and κ 3 = 20. Following the initial transient, the configuration rotates on S 2 with small variations in the relative positions of the nodes. For small values of |κ 3 | the particle motion is asynchronous, suggesting that there is a critical value κ c such that practical synchronization occurs only for |κ 3 | > κ c > 0.

An exact solution for N = 3
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 ω i = 0 and κ 3 /N 2 = 1, (16) reduces to: together with the cyclic permutations. Denote Consider now the fixed points of the system (21). These satisfy x 1 × x 2 = 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, x i ×x j = 0, or x 123 = 0 which implies that the three vectors x i /x 123 form an orthonormal set. In the first case, x 1 , x 2 , 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 x 0 i are consistent with this, i.e. only if |x 0 12 | = |x 0 13 | = |x 0 23 | = 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 = 0, we must have |x 123 | = 1 and so {x 1 , x 2 , 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 2 123 in terms of x 12 , x 23 , x 13 by means of a well-known identity obtained as follows: define the 3 then x 2 123 = p(u). Fromu = −4u x 123 we obtainu 2 = 16u 2 x 2 123 = 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 equationu 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ü = V ′ (u) with u(0) = u 0 ,u(0) = −4u 0 x 0 123 . Let us determine the properties of p and hence of V . We have p(0) = 1, p(1) = −(c 1 − c 2 ) 2 and p(−1) = −(c 1 + c 2 ) 2 , hence p has a real root in (0, 1], denoted r + , and a second real root in [−1, 0), denoted r − . The third root r 3 = −1/(2c 1 c 2 r − r + ) is therefore also real, and is positive or negative according to the sign of c 1 c 2 , and lies outside the interval (−1, 1). The possible values of u(t) for all t > 0 are restricted by the condition p(u) = x 2 123 0, which implies that −1 r − u(t) r + 1 for all t > 0. We factorize p according to p(u) = 2c we can obtain u as an explicit elliptic integral of the third kind, see for example the expression in Gradshteyn and Ryzhik [46], Section 3.137, item (3.) with r = 0.
We can determine how the solution behaves, however, without an explicit expression for u. There is a mechanical analogy with a particle moving under the influence of the potential V with the Lagrangian L = T − V = 1 2u 2 − V (u), with the corresponding equation of motionü = V ′ (u), but with boundary conditions such that the first integral u 2 = 2V (u) holds. The potential V (u) = 8u 2 p(u) has zeroes in [−1, 1] at u = 0 and at u = r ± , and for all initial values, except u 0 = r ± , the particle settles into the stable minimum of V at u = 0. An example of V is shown in Fig. 2(a) for the parameters c 1 = − 3 4 , c 2 = 1 3 , corresponding to a specific choice for the initial values x 0 1 , x 0 2 , x 0 3 . The motion of a particle located at u is restricted to lie between the limits r ± , the roots of p as shown in Fig. 2(a), and are given explicitly for this example by r − = −0.905, r + = 0.703. These points are static solutions ofu 2 = 16u 2 p(u), but do not correspond to any fixed point solutions of the full system (21), unless |c 1 | = |c 2 | = 1. They are, in any case, unstable fixed points of the combined equations for u and x 123 , for any c 1 , c 2 . From x 2 123 = p(u) andu = −4u x 123 we obtainẋ 123 = −2u p ′ (u), the right-hand side of which is non-negative for all u, and so x 123 is an increasing function of t. In particular,ẋ 123 is nonzero at either endpoint r ± , and so these endpoints are unstable static solutions oḟ u 2 = 16u 2 p(u). If a particle initially moves to one of these endpoints, it simply bounces off the barrier imposed by the condition r − u r + , and then approaches the stable equilibrium point at u = 0. For the potential V in Fig. 2(a), we have plotted the explicit solutions u(t), x 123 (t) for the initial values u 0 = 1 2 , x 0 123 = − p(u 0 ) in Fig. 2(b), showing x 123 (t) (in red) as an increasing function of t with x 123 (t) → 1, and u bouncing off the point at u = r + and then approaching the asymptotic limit of zero, which is the stable minimum of V . The fact that x 123 is an increasing function of time is a general property of the gradient formulation of the system, since by regarding the potential V d as a function of time, we haveV We deduce that for all initial values x 0 1 , x 0 2 , x 0 3 , except for the unstable fixed points with |c 1 | = |c 2 | = 1, the system synchronizes to x 12 = x 13 = x 23 = 0. The vectors x 1 , x 2 , x 3 therefore form an orthonormal set which by means of a rotation we can align with the X, Y, Z axes (up to a left-right orientation), and so the nodes lie in the plane x + y + z = 1 with the unit normal n = 1 √ 3 (1, 1, 1). As before, we obtain r ∞ = 1 √ 3 . The relation (20) is satisfied with λ = 2N √ 3 cot π N = 2. Although we have determined the rotational invariants x 12 , x 13 , x 23 , x 123 as functions of t, we omit the further step of solving (21) explicitly for x 1 (t), x 2 (t), x 3 (t), which requires us to choose a suitable basis, since synchronization of (21) follows from the properties of the rotational invariants and the fixed points.

Combined 2and 3-body interactions on the 2-sphere
Of particular interest with all higher-order interactions is how they compare with the well-studied 2-body interactions, and whether one is dominant in some sense. The 2body and 3-body systems for d = 3, given by (2) and (16) respectively, are ideal for such an investigation since, separately, they each synchronize although under different conditions, and generally with incompatible properties. The 2-body system (2), with identical frequency matrices Ω i = Ω and global coupling a ij = 1, completely synchronizes for κ 2 > 0, with the order parameter r evolving exponentially quickly to the asymptotic value r ∞ = 1 [22,21]. Evidently all nodes strongly attract each other from all initial values. For κ 2 < 0, however, the system evolves to a final state with r ∞ = 0, in which the particles are distributed over S 2 [21], and so for this case the interactions are repulsive, a property well-known also for the Kuramoto model with d = 2. By contrast, the 3-body system (16) with identical frequencies ω i = ω, synchronizes for any κ 3 from generic initial values to a final configuration with r ∞ = 1 √ 3 , in which the particles are equally separated as shown in Figure 1 (left), which we refer to as ring synchronization. For combined 2-body and 3-body interactions the question arises, therefore, if the initial particles are clustered together in a tight bunch, do they synchronize to a point under the influence of the attractive 2-body interactions, or do they separate according to the 3-body forces to form a ring on S 2 ? As will be shown, the system synchronizes for any coupling coefficients, whether positive or negative, and the 3-body forces prevail unless the 2-body coupling is positive and sufficiently strong by comparison to the 3-body coupling. The specific system we consider, the d = 3 case of (12), is given by: and properties of the system depend on the coupling constants only through the ratio κ 2 /κ 3 . We find that (23) achieves ring synchronization for all positive and negative values of κ 2 , κ 3 , unless κ 2 /κ 3 exceeds a critical value, specifically, (23) completely synchronizes only for κ 2 /|κ 3 | > 2 N cot π N , otherwise the 3-body forces prevail. The system ring synchronizes for all κ 2 < 0 for which the 2-body forces are repulsive, as we show next.

Steady state solutions on S 2
Numerically we find that solutions of (23), for generic initial values, approach a static configuration. The vectors x i = ± u are fixed points of (23), where the sign can depend on i, and indeed for 2-body interactions with κ 3 = 0 and κ 2 > 0 we obtain complete synchronization with plus signs for all i [21], and then x i = u is a stable fixed point. We can also obtain bipolar synchronization, in which the signs vary, if we include multiplicative factors of any sign in the 2-body interactions [48,38].
The following configuration, which generalizes (17), is a steady state solution of the combined system (23): where r ∞ , α, to be determined, are parameters satisfying r 2 ∞ (1 + α 2 ) = 1, which ensures that x i is a unit vector. As before, we have used the rotational covariance of (23) to rotate any planar static configuration to lie in the plane z = r ∞ , together with an SO(2) rotation about the Z-axis so that x N = r ∞ (α, 0, 1) is aligned along the X-axis. The rotational invariant corresponding to (24) is and so again the nodes are equally spaced. We determine r ∞ , and hence α, as a function of κ 2 /κ 3 by requiring that the steady state (24) should satisfy (23). Let us firstly verify that the parameter r ∞ in (24) corresponds to X av . We have where we have used (A.4), and hence X av = r ∞ as required. Also, from (24): and By means of the summation formulas in Appendix B we obtain: for all i. Evidently (27) corresponds to the general formula (13) with λ 1 = 2Nr ∞ cot π N and λ 2 = N cot π N (1−3r 2 ∞ )/r ∞ . Hence (23) is satisfied provided that (14) holds, i.e. provided that κ 2 r ∞ + κ 3 (1 −3r 2 ∞ ) 1 N cot π N = 0. We have established therefore that (24) is a steady state solution of (23), which we find numerically to be stable for κ 3 > 0, under conditions to be determined. For κ 3 < 0 the stable solution is found by means of the parity change x i → −x i in (23), by reversing the sign of κ 3 . The condition therefore that (24), or its negative, be a stable steady state for (23) is: which leads to the following explicit expression for r ∞ as a function of κ 2 /κ 3 : Evidently r ∞ is positive for any sign combination of κ 2 , κ 3 , and for κ 2 = 0 we regain the value r ∞ = 1 √ 3 as discussed in Section 3. If κ 3 = 0 the vector (24) is not a fixed point of the 2-body system, as (26) shows directly, unless r ∞ = 1, α = 0.

Transition to complete synchronization
The formula (30) for r ∞ violates the bound r ∞ 1 if κ 2 /|κ 3 | is too large. The critical ratio is found by setting r ∞ = 1 in (29), hence we require in order for the static solution (24) of the combined system (22) to exist. We can now determine the behaviour of the system as κ 2 /κ 3 varies. As κ 2 /|κ 3 | approaches the critical ratio from below, it follows from (30) that r ∞ → 1, which means that the particles become more tightly bunched, and so the synchronized ring shrinks in diameter and eventually reduces to a point at which equality holds in (31). Complete synchronization occurs for all larger values of κ 2 /|κ 3 |, i.e. the 2-body forces prevail through a continuous transition. Otherwise, whenever κ 2 /|κ 3 | is less than 2 N cot π N , ring synchronization prevails, in particular as κ 2 decreases to large negative values, we have r ∞ → 0 as is evident from the formula (30), and so the ring expands to a great circle on S 2 . In this case the 3-body forces overcome the repulsive 2-body forces, however large, to maintain ring synchronization.
Numerical calculations confirm these properties. Fig. 3(a) shows ring synchronization for the system (23) with N = 40, κ 3 = 2, κ 2 = 1, for which κ 2 /|κ 3 | is less than the critical ratio 0.635, and numerically we find r ∞ = 0.896, in agreement with (30). Evidently, the attractive 2-body forces shrink the ring of synchronized nodes down to a smaller diameter, compared to that in Fig. 1(a). If we include frequency terms ω i ×x i on the right-hand side of (23), then exact synchronization is replaced by practical synchronization, in which the final configuration, now time-dependent, resembles either the ring or the completely synchronized configuration that occurs for ω i = 0, depending on κ 2 /κ 3 . As an example, in Fig. 3(b) we show the effect of nonzero frequency vectors ω i such that ω i < 1/20, for the same parameters as in Fig. 3(a), where evidently the synchronized ring is slightly distorted by the small nonzero oscillations at each node.
Practical synchronization occurs even for κ 2 < 0, when the 2-body forces are repulsive, provided that |κ 3 | is sufficiently large to overcome the natural oscillations. If we fix κ 2 = −1 and choose the same frequencies ω i as in Fig. 3(b), then the system synchronizes with r ∞ ≈ 0.366 for κ 3 = 2, and for any larger values, but not for κ 3 = 1, which suggests that there exists a critical value for |κ 3 |, in order for synchronization to occur. Again, the 3-body forces enhance synchronization, since without 3-body interactions of sufficient strength the system would evolve asynchronously for any κ 2 < 0.

4-body interactions on the 3-sphere
The 4-body forces have properties which differ from those for d = 3, consistent with previous observations that features of even-dimensional systems on the sphere differ from those for odd dimensions [20,36]. Whereas the d = 3 system always synchronizes with the same value r ∞ = 1 √ 3 , now r ∞ depends on N, and although the nodes arrange themselves to lie equally spaced, in this case on a torus imbedded in S 3 , the nodes do not form a closed loop. This also occurs for d = 2, see for example Fig. 5(a).

The homogeneous system
Consider firstly the homogeneous 4-body systeṁ where the antisymmetric vectors v jkl are defined by (9). We can reverse the sign of κ 4 by means of an element R ∈ O(4) such that det R = −1, and by rescaling we set κ 4 /N 3 = 1.
Numerically, we find that the system (32) synchronizes from generic initial values to a final static configuration in which the scalars n x i = n i , where n = X av / X av , depend on i, always with n i > 0. Hence the final configuration is confined to one hemisphere of S 3 but, unlike the d = 3 case, does not lie in a hyperplane which intersects S 3 . The following rotationally invariant formula for the static final configuration holds for all i, j: which implies that the nodes are equally spaced, i.e. x i − x i+1 is independent of i for i = 1, . . . N − 1. The symmetry x i x N = −x N −i x N also follows. We deduce the following expression for x i , up to a constant rotation: These N nodes do not form a closed loop since for i = 0 we have x 0 = (1, 0, 1, 0)/ √ 2 which equals −x N , rather than x N . The configuration (34) consists of N equally spaced nodes lying on a torus imbedded in S 3 , with the first two components restricted to a semicircle, while the last two components trace out a full circle plus a semicircle for i = 1, . . . N.
The expression (34) satisfies (32) exactly, as a consequence of the following relation, which is the d = 4 case of (11): which we prove in Appendix C. From (34), with the help of (A.5), we find: from which we obtain Evidently, r ∞ is a function of N, decreasing from its maximum value r ∞ = 1 2 at N = 4, to a minimum which is attained as N → ∞ for which r ∞ → 2 √ 5/(3π) ≈ 0.4745. Although the nodes (34) lie on S 3 , we can visualize the configuration by omitting one of the four components, then normalizing the remaining components to form a unit 3-vector on S 2 . We plot such a configuration in Fig. 4(a) for N = 80, in which the second component has been dropped, showing equal spacing with endpoints that are diametrically opposite on S 2 . Another way of visualizing the evolving system is by means of the Hopf fibration in which S 3 is mapped to S 2 according to (x, y, z, w) → 2(xz − yw), 2(xw + yz), In general, this maps circles on S 3 to points on S 2 , but here distinct points on S 3 are mapped one-to-one or two-to-one to points on S 2 , which enables us to visualize the evolving system as it synchronizes. The Hopf fibration and therefore the final configuration, however, does not respect the SO(4) covariance of the system (32) and so the final nodes are generally not arranged in a symmetrical configuration. Equal spacing of the nodes, for example, which follows from the rotationally invariant expression (33), is not evident. However, if we rotate the final configuration to obtain (34), then under the Hopf fibration (37) we find that x i → cos 4πi N , sin 4πi N , 0 , and so for odd N the nodes are mapped to N equally spaced distinct points on the equator of S 2 , and to N/2 distinct points for even N.

Combined pairwise and 4-body interactions on the 3-sphere
The combined 2-and 4-body system is given by (12), namely: Numerically this system synchronizes in a way similar to the d = 3 system (23), in that the 4-body forces dominate so that the final configuration is similar to that with κ 2 = 0, unless the ratio κ 2 /κ 4 is sufficiently large, in which case complete synchronization occurs. We find that the expression (33) generalizes to: where α, β, θ are unknown parameters. This relation can be verified numerically to high accuracy for suitably fitted parameters α, β, θ, indicating that the form (39) is exact, however there is no simple functional relation between α and β. For the homogeneous case (34) we have θ = π 4 and α = β = 1, but otherwise α, β satisfy nontrivial trigonometric relations that depend on N. Consistent with (39), we have and the relations between α, β, θ are determined by requiring that (13) be exactly satisfied for constants λ 1 , λ 2 which depend on α, β, θ. Numerically, we find indeed that (13) is satisfied in all cases, for constants λ 1 , λ 2 which are independent of i for all initial values. It is possible in principle to derive exact expressions for all quantities of interest by means of (40), as we discuss briefly in Appendix C, however the resulting equations are of such complexity, due to the fact that α, β are not integers, that we confine ourselves here mainly to numerical observations. One way of visualizing the configuration (40) is by means of the Hopf fibration (37), for which showing that the nodes form a ring of equally spaced points around the 2-sphere in the plane z = cos 2θ. The nodes become more tightly spaced as α, β decrease, in which case r ∞ increases. We can determine the following expression for r 2 ∞ from either (39) or (40): which reduces to (36) for θ = π 4 , α = β = 1. We wish to determine the properties of the system as the 2-body forces increase in strength relative to the 4-body forces. For d = 3 we found that a continuous transition occurs at a critical ratio, from ring synchronization to complete synchronization, but here for d = 4 we find that this transition is discontinuous, i.e. as κ 2 /κ 4 increases the system suddenly completely synchronizes at, and beyond, a critical ratio. The order parameter r 2 ∞ as given in (41) attains the value unity, signifying complete synchronization, only if both α, β → 0, as can be deduced from (41), but numerically we find that neither α nor β approach zero as the transition point is approached, which is indicative of a discontinuity. Other properties of the combined system (38) are similar to the d = 3 case, for example the system synchronizes to a configuration satisfying (39) for all negative values of κ 2 . If we fix κ 4 and allow κ 2 to decrease to large negative values, r ∞ decreases accordingly, with the equal spacing of nodes in all cases maintained under the repulsive 2-body forces. Again, the 4-body forces prevail over the repulsive 2-body forces so as to enforce synchronization.

Combined 5-body and 2-body interactions on S 4
The largest value of d that we consider is d = 5 for the combined system (12), in order to confirm that it has properties similar to the d = 3 case. The main difference compared with d = 3 is that as the 2-body forces increase in relative strength, the transition to a completely synchronized configuration is discontinuous. We have the equations: and again find numerically that the system synchronizes for all generic initial values, either to a configuration with equally spaced nodes, even for negative κ 2 , or else to a point if κ 2 /κ 5 is sufficiently large. This holds for either sign of κ 5 , which can be reversed by replacing x i → −x i .

The homogeneous 5-body system
For κ 2 = 0 the system synchronizes to a static configuration in which all nodes lie in a hyperplane defined by n x i = r ∞ = 1 √ 5 where n = X av /r ∞ , as occurs for d = 3. These configurations satisfy from which we deduce that: for which X av = (0, 0, 0, 0, 1 √ 5 ). The equally spaced nodes form a closed sequence since x 0 = x N . We can write these vectors as x i = 1 √ 5 (2u i , 1) where u i ∈ S 3 satisfies i u i = 0, corresponding to balanced spacing. Specifically: These unit vectors are similar to those encountered for d = 4 in (34), which lie on a torus in S 3 , except that here the nodes form a closed sequence. We can visualize u i by deleting one component, and normalizing the remaining components to unity, and then plot the resulting configuration on S 2 , as for d = 4. As an example, we have plotted u i in Fig. 4(b) for N = 80, in which the last component has been dropped, showing a closed sequence of points in S 2 , in contrast to the d = 4 case in Fig. 4(a). Under the Hopf fibration (37) we have u i → cos 6πi N , sin 6πi N , 0 and so all nodes are mapped to points on the equator.

The 5-body system with 2-body forces
Of particular interest is the combined system (42) with κ 2 = 0, by means of which we can investigate the relative effect of 5-body and 2-body forces. We find, similar to the d = 3, 4 cases, that 5-body forces prevail over the 2-body forces unless κ 2 is positive and sufficiently large compared to |κ 5 |, in which case the system completely synchronizes. The transition to complete synchronization is discontinuous, as for d = 4, except that here we are able to determine the precise conditions under which this occurs.
The synchronized configuration takes the form where the parameters r ∞ , α satisfy r 2 ∞ (1 + α 2 ) = 1. This expression reduces to (44) for r ∞ = 1 √ 5 and α = 2. It follows that which shows that the nodes are equally spaced, and since x 0 = x N the nodes form a closed sequence. We require the configuration (46) to satisfy the static equations (13), namely: N j,k,l,m=1 which determines r ∞ , and hence α, as functions of κ 2 /κ 5 . We can verify that this equation is indeed satisfied, and calculate λ 1 , λ 2 , see Appendix C, to obtain (50) According to (14), the static equations (42) are satisfied by (48) provided that λ 2 = κ 2 N 4 /κ 5 , hence we obtain the following equation which fixes r ∞ as a function of κ 2 /κ 5 : where we have included the absolute value |κ 5 | so as to allow for negative κ 5 , noting that the signs of κ 5 and λ 1 , λ 2 , change under x i → −x i . The configuration (46) exists therefore only if (51) holds. The case κ 2 = 0 corresponds to r ∞ = 1 √ 5 , whereas for r ∞ = 1 we have α = λ 1 = λ 2 = 0 and (46) becomes an unstable fixed point. The relation (51) can be satisfied for any negative value of κ 2 , however large, corresponding to arbitrarily small values for r ∞ , and numerically we find indeed that the system synchronizes such that (47) and (51) are satisfied in all cases.
We conclude that 5-body forces enhance synchronization in a way similar to d = 3, 4, in that synchronization occurs for all coupling constants κ 2 , κ 5 in any sign combination.

The symmetric-antisymmetric Kuramoto models
The simplest of the hierarchy of models (10) is the case d = 2 which can be termed the antisymmetric Kuramoto model which, it appears, has previously escaped attention, but is useful as a guide to possible behaviours for d > 2. We wish to compare the properties for antisymmetric couplings ε ij to those of the standard Kuramoto model with symmetric couplings a ij = 1, even though in both cases we have 2-body interactions.
We parametrize x i = (cos θ i , sin θ i ), then as defined in (9), v i = (sin θ i , − cos θ i ), and the potential (6) is given by V 2 = N i,j=1 ε ij sin(θ j − θ i ). The equations of motion arė where we have included distributed frequencies ω i , and have denoted the corresponding coupling constant by κ a which, by comparison with the general system (10), is equal to −κ d for d = 2. By writing cos(θ j −θ i ) = sin(θ j −θ i + π 2 ) we see that this model is similar to the well-known Sakaguchi-Kuramoto model [45] with a frustration angle α = π 2 , except that here the couplings are antisymmetric. The condition |α| < π 2 is imposed in [45] in order to ensure that the system synchronizes, however the antisymmetric couplings in (52) allow the system to synchronize for both positive and negative values for κ a , for any frequencies ω i and for any frustration angles, provided that |κ a | exceeds a critical value κ c . The transformation θ i → −θ i together with ω i → −ω i in (52) is equivalent to κ a → −κ a , hence the behaviour of the system (52) is indifferent to the sign of κ a . We will see that the combination of symmetric and antisymmetric couplings enhances the synchronizability of the system.

Synchronization for identical frequencies
Taking firstly the case ω i = 0 for all i, then we find numerically that the system (52) synchronizes from random initial values θ i (0) for any positive or negative κ a to a static configuration in which the nodes are equally spaced on the half circle, as shown for example in Fig. 5(a) for N = 40. Such configurations are similar to splay states [47,44], except that here they are restricted to the half-circle. A steady state solution of (52) with ω i = 0 is given by where θ 0 is a constant angle, as we now show. Let ζ = exp iπ N , then by means of the geometric series we obtain hence N j=1 ε ij cos (j−i)π N = 0 and so (53) is a steady state solution of (52) with ω i = 0. Let us also verify that the relation (11) holds. We have, choosing θ 0 = 0 in (53): then from (54) we deduce that N j=1 ε ij v j = cot π 2N x i , in accordance with (11). The value of r ∞ is given by (55) Numerical results show that the steady state (53), or its negative depending on the sign of κ a , is stable since it is attained from generic (random) initial values. The fixed point θ i = θ 0 corresponds to a completely synchronized configuration but is evidently unstable, consistent with previous observations for d > 2. There is some similarity with the case d = 4 which we considered in Section 5, since in both cases r ∞ depends on N, see (36), and the synchronized nodes form a closed loop only with the anti-periodic extension in which we include the nodes −x i . For distributed nonzero frequencies ω i the system again synchronizes provided that |κ a | > κ c for some critical value κ c , where κ a can be positive or negative, with a phaselocked frequency Ω = 1 N j ω j . The particles in such a synchronized configuration are no longer exactly spaced apart but, as κ increases, the configuration approaches the steady state shown in (53).

Combined symmetric-antisymmetric models
Of particular interest is the effect of the antisymmetric couplings in (52) compared to the standard symmetric interactions in the Kuramoto model, and whether the antisymmetric couplings enhance synchronization, as occurs for 3-body interactions. Indeed, we find that the antisymmetric interactions dominate the usual symmetric interactions, for example whereas the Kuramoto model does not synchronize for negative coupling constants, this is over-ridden by antisymmetric couplings of any sign, for which the system always synchronizes to a final configuration with the nodes distributed around the circle. We consider therefore the systeṁ with independent coupling constants κ s , κ a of any sign, which govern the relative strengths of the symmetric and antisymmetric interactions respectively. Consider firstly the case ω i = 0 for all i then, as is well-known, if κ a = 0 the system completely synchronizes for κ s > 0, i.e. all particles are co-located at a single point with r ∞ = 1, and if κ s < 0 the particles are distributed around the circle such that r ∞ = 0. However for general κ s , κ a , whether positive or negative in any combination, the system always synchronizes with the final state resembling that for κ s = 0, i.e. the nodes are equally spaced over an arc of the unit circle. The steady state solution in this case is of the form where θ 0 is a constant angle and α is an unknown parameter which is determined by substituting θ i into (56). Let us verify in this case that the relation (13) is satisfied, and hence find the condition which fixes α. We have firstly, with the help of (A.2,A.3): then (13) reads N j=1 ε ij v j = λ 1 x i − λ 2 X av , which is satisfied for all i with Hence (57) is a steady state solution of (56) provided that (14) holds, namely: where we have replaced the symmetric coupling coefficient κ 2 → κ s in (14), together with κ d → −κ a , consistent with the sign of κ a in (56). The value of r ∞ , calculated from (58) or with the help of (A.6), is: which generalizes the formula (55) for which α = ±1.
The relation (59) shows how the system behaves as the coupling ratio κ s /κ a varies, in particular how the spacing between nodes changes, as measured by α. If κ s > 0 we set α = − 2 π arctan κa κs , but if κ s < 0 numerical results show that we must choose either α = −2− 2 π arctan κa κs for κ a > 0, or α =2− 2 π arctan κa κs for κ a < 0. The case κ s = 0 corresponds to α = ∓1 according to the sign of κ a , since either (53) or its negative is a stable fixed point.
If κ s is large and positive, for either positive or negative κ a , (59) shows that α is small, with r ∞ close to unity since from (60) we have r ∞ → 1 as α → 0. Hence the particles are tightly bunched together, but are always equally spaced. Complete synchronization never occurs, even for large values of κ s , in contrast to the d = 3 case, see Section 4.1. If κ s is small, whether positive or negative, the configuration is similar to that when κ s = 0, since |α| is close to unity. If κ s is large and negative, |α| increases, but with |α| < 2, and so the nodes spread out to occupy almost the full circle. As an example, we have plotted such a configuration for N = 40 in Fig. 5(b) with κ a = κ s = −1, for which α = 3 2 . The κ s coupling therefore either increases the length of the arc on which the synchronized points lie, corresponding to a repulsive effect for κ s < 0, or binds the nodes together more tightly for κ s > 0.
Consider finally the case of distributed frequencies ω i , then synchronization occurs with a phase-locked frequency Ω = 1 N j ω j , but only if at least one of κ a , κ s is sufficiently large in magnitude. The Kuramoto model does not synchronize if the coupling coefficient is negative, but now we find that synchronization does indeed occur if |κ a | is sufficiently large, i.e. the antisymmetric couplings override the repulsion of nodes arising from the negative symmetric couplings. These synchronized configurations are similar to those encountered in models with distributed phase lag, which can be understood if we write (56) in the forṁ The main difference between this system and the Sakaguchi-Kuramoto model with distributed phase-lag is that the angles are antisymmetric, a consequence of which is that (61) has a gradient formulation, unlike the Sakaguchi-Kuramoto model for which symmetric angles are usually assumed.
As an example of a synchronized configuration for the system (61), equivalently (56), we show in Fig. 5(c) a final synchronized state with N = 40 for κ a = κ s = −5, with distributed frequencies such that |ω i | < 1. Since κ s < 0 this system does not synchronize unless we include a nonzero κ a coupling of sufficiently large magnitude, and then synchronization occurs even if κ s is very large and negative.

Conclusion
We have constructed a hiercharchy of models on the unit sphere in d dimensions which have either exclusive d-body interactions or combined d-and 2-body interactions. We have shown that such systems synchronize under very general conditions to a steady state with equally spaced particles on S d−1 , unless the 2-body forces are sufficiently strong in which case the system completely synchronizes. We have derived exact expressions for the synchronized configurations for d 5, and have calculated for d = 3, 5 the exact ratio of coupling constants at which the system transitions to a completely synchronized configuration. Our findings support previous observations that higher-order interactions enhance the synchronization of the combined systems even when the pairwise couplings are repulsive, but we have not observed other reported phenomena such as multistability or explosive synchronization.
Some of the properties that we have derived for d 5 generalize without difficulty to any d, for example in the simplest case in which N = d, an exact static solution of the d-body system (10)  Our aim has been to determine the dynamical characteristics of higher-order interactions with antisymmetric coupling coefficients, and it remains to establish properties of more general couplings with their associated network structures. There are also many mathematical questions to be finalized, such as the stability of the asymptotic configurations, as well as possible applications in this vast area of research.

Appendix A.
We derive here several well-known trigonometric summation formulas. Let ζ = exp iπ N then ζ N = −1, and from the geometric series, q k=p x k = (x p − x q+1 )/(1 − x), for any parameter α: and so From (A.1) we also obtain: which, since the right-hand side is real, is equal to N i,j=1 cos απ(j−i) N .

Appendix B.
Here we prove that (20) and (27) are satisfied by the vectors (17). Firstly, it follows from (17) that and we wish to evaluate N j,k=1 ε ijk (x j × x k ). We first evaluate N k=1 ε ijk ζ −2k where ζ = exp iπ N . If i < j then ε ijk = 1 for 1 k < i, ε ijk = −1 for i < k < j and ε ijk = 1 for j < k N. By means of the geometric series we obtain for i < j: which is antisymmetric in i, j, and therefore holds for all i, j. By summing over j and using N j=1 ζ −2j = 0, as follows from (A.1) with α = −2, we obtain: Hence: where we swapped the j, k summations, as well as Also, from (B.2): again using N j=1 ζ 2j = 0, and so 1 N N j,k=1 ε ijk sin 2(j − k)π N = − cot π N .
By application of the above formulas we evaluate N j,k=1 ε ijk (x j × x k ) using (B.1) to obtain (20) and (27).

Appendix C.
We prove here (35) for d = 4 and similar relations for d = 5. We first calculate the components of v jkl from the identity u v jkl = det(u, x j , x k , x l ) for an arbitrary vector u and find for the first component, with the help of a computer algebra system: 4 √ 2 (v jkl ) 1 = cos (3j − k − 3l)π N − cos (3j + k − 3l)π N + cos (j + 3k − 3l)π N − cos (3j − 3k − l)π N + cos (3j − 3k + l)π N − cos (j − 3k + 3l)π N , (C. 1) and similarly for the other three components. We wish to evaluate N j,k,l=1 ε ijkl v jkl . For the first component we obtain from (C.1): where we have used the antisymmetric properties of ε ijkl to combine the six terms on the right-hand side of (C.1) into a single term. In order to evaluate this sum, let ζ = exp iπ N then with methods similar to those in Appendix B we derive N k,l=1 which is valid for all i, j and all odd integers α, β such that α + β = 0. For α = −1, β = −3 we obtain N j,k,l=1 The real part of (C.4) now reads: which establishes the first component of (35). The second component is verified by taking the imaginary part of (C.4), and a similar calculation holds for the two remaining components.
For the more general expression (40) for x i , which includes the effect of 2-body forces, we calculate the components of v jkl as before, and again find that each component is a sum of six terms, generalizing (C.1). By means of the antisymmetry properties of ε ijkl we obtain an expression which reduces to (C.2) for θ = π 4 , α = β = 1, namely: We can evaluate this expression in principle by summing N j,k,l=1 ε ijkl ζ 3βj−3βl−αk , and then extracting the real part, however because α, β are not integers, the result does not simplify to an expression such as (C.4). Due to the resulting complexity we resort to numerical observations for this case.
For d = 5 we wish to prove that (46) satisfies (48), and find explicit expressions for λ 1 , λ 2 . We firstly evaluate v jklm by means of u v jklm = det(u, x j , x k , x l , x m ) again with the help of a computer algebra system, to obtain an expression for each component, similar to (C.1) for d = 4, except that there are now 24 terms, instead of 6, on the righthand side. When summed over the signature, however, these 24 terms can be reduced to a single term, similar to (C.2) and (C.6) for d = 4, by using the antisymmetric properties of the signature function, hence we obtain N j,k,l,m=1 It is sufficient to consider only the first component of (48) in order to find λ 1 , since we have X av = (0, 0, 0, 0, r ∞ ). Having found λ 1 we then determine λ 2 by considering the 5th component of (48). It is easiest to first find λ 1 by setting i = 0 in (C.7), then by means of sums such as (C.3) we obtain: N j,k,l,m=1 from which we can identify λ 1 as given in (49), and similarly obtain λ 2 . Expressions such as (48) can also be verified exactly for specific, but small, N by means of a computer algebra system, as well by high-accuracy numerical computations.

Appendix D.
We show here how to rotate static configurations, such as that shown in Fig. 1(a), to a fixed orientation so as to directly compare final states for different initial values.
By equating the right-hand side of this expression at t = 1 with the computed normal n, we can identify ω, and therefore obtain the rotation matrix R with the required property Rn = m, where R is given in explicit form by R = e −Ω . Hence we can rotate a configuration x i at any fixed time for all i according to x i → Rx i with the property that X av aligns with m.