Continuous versus Discontinuous Transitions in the $D$-Dimensional Generalized Kuramoto Model: Odd $D$ is Different

The Kuramoto model, originally proposed to model the dynamics of many interacting oscillators, has been used and generalized for a wide range of applications involving the collective behavior of large heterogeneous groups of dynamical units whose states are characterized by a scalar angle variable. One such application in which we are interested is the alignment of orientation vectors among members of a swarm. Despite being commonly used for this purpose, the Kuramoto model can only describe swarms in 2 dimensions, and hence the results obtained do not apply to the often relevant situation of swarms in 3 dimensions. Partly based on this motivation, as well as on relevance to the classical, mean-field, zero-temperature Heisenberg model with quenched site disorder, in this paper we study the Kuramoto model generalized to $D$ dimensions. We show that in the important case of 3 dimensions, as well as for any odd number of dimensions, the $D$-dimensional generalized Kuramoto model for heterogeneous units has dynamics that are remarkably different from the dynamics in 2 dimensions. In particular, for odd $D$ the transition to coherence occurs discontinuously as the inter-unit coupling constant $K$ is increased through zero, as opposed to the $D=2$ case (and, as we show, also the case of even $D$) for which the transition to coherence occurs continuously as $K$ increases through a positive critical value $K_c$. We also demonstrate the qualitative applicability of our results to related models constructed specifically to capture swarming and flocking dynamics in three dimensions.


A. Background
Collective behavior in large populations of interacting elements has been a subject of intense study in physical, social, biological and technological systems [1][2][3][4][5][6][7][8][9][10]. An important, frequently encountered example is the case of interacting phase oscillators, i.e., coupling between elements whose state is characterized by a point on a unit circle. In 1967 Winfree first systematically studied the dynamics of a population of weakly coupled phase oscillators [11]. A few years later [12], Kuramoto presented a simplified version of the Winfree model which he solved in the limit of N → ∞, where N is the number of oscillators. This model, now known as the Kuramoto model, is where θ i represents the phase angle of the i th oscillator, ω i is its natural frequency of oscillation (which we will also refer to as the natural rotation), and K is the coupling strength between oscillators. Typically the ω i are chosen randomly from some unimodal distribution with a finite spread ∆, and N 1 (the case of interest in this paper) is often considered. In the N → ∞ limit * sarthakc@umd.edu Kuramoto was able to show the presence of a continuous phase transition between asynchronous and partially synchronous states of the system [12,13].
The Kuramoto model and its generalizations have since been used to study synchronization behavior in a wide variety of systems, modeling biological problems such as the behavior of cardiac pacemaker cells [14], synchronization in large groups of flashing fireflies [15,16], circadian rhythms [17,18], and neuronal synchronization [19], as well as problems in physics and engineering such as synchronization of power-grid networks [4,9], superconducting Josephson junctions [20], atomic physics [21], and neutrino oscillation [22], among others. Another class of applications of the Kuramoto model has been modeling the alignment of unit vectors representing the direction of motion of interacting members of a swarm or flock of moving agents in two dimensions [23][24][25]. Alternately, one can think of such unit vectors as characterizing the opinion of an individual in a group of interacting individuals [26]. In this later case, alignment of unit vectors can be viewed as modeling the evolution toward social consensus.
The aforementioned studies all describe the alignment of agents via a single scalar variable θ i , which characterizes the alignment state of the individual coupled units. However, for several of the above cited applications alignment in higher-dimensional spaces is important, and this is the subject of this paper. For example, the problem of alignment of velocity vectors in a flock of birds, a school of fish, or a swarm of flying drones is more realistically set up in three-dimensional space, whereas the alignment of opinion dynamics of a population could in general be multidimensional depending on the characteristics of the opinions considered. With such examples in mind, Olfati-Saber [26] introduced a higher-dimensional generalization of the Kuramoto model without the presence of any individual natural rotation [analogous to the ω i term in Eq. (1)]. (In 2013, Zhu [27] considered the equivalent case of identical natural rotations for each agent.) The choice of the generalization in Refs. [26,27] maintains the form of the coupling between two agents in all dimensions, i.e., in D dimensions the state of each agent is taken to be a Ddimensional unit vector, and the coupling between two agents is proportional to the sine of the angle between their unit vectors [28]. Network characteristics leading to complete alignment were discussed; however, no complete stability analysis of the system was presented. In our paper we consider globally coupled systems, with a spread of the individual natural rotations of each unit, which follows from the generalization of the spread in the natural frequencies of the standard Kuramoto model. These natural rotations act as constant biases to the states of the agents. In particular, for a given swarming agent, the natural rotation term can be thought of as a systematic error in the dynamics of the agent, which causes the agent to drift away from traveling in purely a straight line. We motivate the form natural rotation term in the context of flocking and swarming in D = 3 in Sec. IV. In assuming these natural rotations we set up a model more general than the one that has been studied by previous authors, leading to new and interesting results.

B. Main Result
A key point in this paper is the remarkable difference between the standard two-dimensional Kuramoto model and its generalizations to 3 dimensions (and also to odd values of D ≥ 5). A striking example of this is the nature of the transition from the incoherent state to the partially aligned state. As previously noted, the twodimensional Kuramoto model, in the limit of infinite system size, was shown by Kuramoto [12,13] to exhibit a continuous phase transition to coherence with increasing coupling strength K. This is represented by the dashed curve in Fig. 1, where the horizontal axis is the coupling strength, K, and the vertical axis represents the 'order parameter' [Sec. II, defined in Eq. (5)], which is a measure of the coherence (or degree of synchronization). The exact shape of this curve can be derived analytically [29], and it can be shown that this phase transition to synchrony is effectively a low-dimensional bifurcation [30]. The three-dimensional Kuramoto model, on the other hand, exhibits a discontinuous phase transition as we increase the coupling strength through zero (solid curve in Fig. 1): For negative values of the coupling strength (indicative of repulsive interactions between agents), the agents tend to a completely incoherent state (defined by an 'order parameter' with zero magnitude), and as we increase the coupling strength through zero, there is a discontinuous jump of the coherence as measured by the order parameter. Further, we find that this discontinuous phase transition occurs nonhysteretically.

C. Relation to Statistical Physics Models
It is interesting to note that if the time-independent frequencies ω i in Eq. (1) are replaced by independent, zero-mean, white noise of uniform strength, then the statistical equilibria and phase transitions of the Kuramoto model are the same as those of the mean-field classical XY model, which describes the interactions of classical two-dimensional spins with global coupling [31][32][33][34]. In this case, the strength of the white noise corresponds to the temperature, and the magnitude of the coherence corresponds to the magnetization. Thus the Kuramoto model can be thought of as the mean-field XY model with thermal noise replaced by quenched randomness (the randomly chosen time-independent frequencies ω i ). Specifically, the mean-field XY model and the Kuramoto model yield similar behavior [19,32] in that they both show a continuous ('second order') transition as the coupling constant increases through a critical value K c > 0 (which, for the Kuramoto model, increases with the spread ∆ in the distribution of the natural frequencies, while, for the XY model, K c increases with temperature). A surprising result of our paper is that, when these models are extended to three dimensions, the two-dimensional qualitative similarity of the behavior for the cases of the quenched randomness and thermal noise versions of the XY model no longer applies: As mentioned above, the three-dimensional Kuramoto model with quenched disorder shows a discontinuous ('first order') phase transition at a zero coupling strength. Independent of the magnitude of the spread in the rotations comprising the quenched disorder, the three-dimensional Kuramoto model always shows partial alignment for K > 0. Since in the three-dimensional model the coupling between any two agents is identical to the two-dimensional case, i.e., proportional to the sine of the angle between the unit vectors σ i of the two agents, this model also describes the interactions of classical three-dimensional spins with global coupling, i.e., the mean-field classical Heisenberg model. If this quenched disorder in terms of the spread of natural rotations were to correspond with the a temporally noisy disordered system, then allowing for larger spread would correspond to higher temperatures and larger noise. However, at finite positive temperature the classical Heisenberg model, like the XY model, has a continuous phase transition at a positive critical coupling strength K c [35]. Thus, in contrast to the twodimensional case, for these problems in three dimensions there is a qualitative difference between temperature and quenched disorder.  [29]), shown as the black dashed curve, and for the Kuramoto model generalized to three dimensions as calculated from the theory in Eq. (18), shown as the solid red curve. Note the continuous transition in the two-dimensional Kuramoto model at a critical coupling of Kc > 0, and the discontinuous transition of the three-dimensional Kuramoto model at Kc = 0. The blue dotted line represents the maximum possible value of coherence, corresponding to |ρ|= 1.

II. MODEL DESCRIPTION
In order to see how the Kuramoto model can be generalized to higher dimensions [26,27], we note that Eq. (1) for θ i can be rewritten (see Fig. 2 and its caption) in terms of the evolution of a collection of N two-dimensional unit vectors, σ i with (x, y) components (cos θ i , sin θ i ): where From this point of view the natural generalization of the Kuramoto problem, Eq. (1), to D dimensions is to consider Eq. (2), but now with σ i being a unit vector in D dimensions and W i being a real D × D antisymmetric matrix. Thus, unlike the standard Kuramoto model where the state of an agent is described by a single scalar variable θ i , the state of each agent in this generalized Kuramoto model is completely described by a D-dimensional unit vector σ i .
Each W i term can be thought of as a constant bias to the dynamics of σ i . In the uncoupled dynamics, dσ i /dt = W i σ i , each agent is acted on by a constant linear operator, which causes each agent to move along the surface of the unit sphere S. For example, in the context of swarms or flocks, it is natural to assume that each agent, in the absence of coupling (K = 0), has some imperfection that causes it to deviate away from the ideal of straight-line steady motion (dσ i /dt = 0). Antisymmetry of W i is imposed so as to ensure that the state vectors σ i are unit vectors at all times.
For example, D = 3, as discussed above, is of particular interest. In this case the term W i σ i can be represented as where ω i is referred to as the rotation vector; see Fig.  3 which schematically represents the solution of Eq. (2) for the case K = 0 and D = 3, in which σ i is shown precessing around the vectorω i = ω i /|ω i | at the rate ω i = |ω i |. (Here, and later in this paper, we use the notation |v| to represent the Euclidean norm of the vector v) Note that the dot product of the right-hand side of Eq. (2) with σ i is identically zero in all dimensions D, implying that d|σ i |/dt = 0, consistent with σ i being a unit vector.
In the context of the spin models discussed earlier in Sec. I C, for positive K, the first term in Eq. (2) corresponds to the interaction term between individual spins, with each pair of spins tending to align themselves parallel to each other. This term leads to macroscopic magnetization in the system of spins. The second term in Eq. (2), W i σ i , corresponds to the quenched disorder discussed in Sec. I C which inhibits coherence among the spins.
In the context of flocking models, each σ i is interpreted as the unit vector along the velocity vector for the i th agent. It is also assumed that the state of the agent is completely described by σ i , i.e., the agent is effectively axisymmetric about σ i . For positive K, the summation term in Eq. (2) corresponds to all-to-all communication between agents in the flock, with each agent tending to align itself with each of the other agents. This term promotes coherence within the flock. The second term, W i σ i , corresponds to a simple dispersing term causing decorrelation of the agent orientations σ i . In particular, if we wish to consider the addition of a dispersal term that maintains the norm of σ i , and for simplicity is assumed to be time independent and linear, then it must be of the form W i σ i for some antisymmetric matrix W i .
In the context of swarms and flocks of threedimensional agents, further motivation and justification for the form of the dispersing term W i σ i is presented in Sec. IV. In particular, in order to support the possible generality of our main result (exemplified in Fig. 1), in Sec. IV we consider another model, different from the generalized Kuramoto model Eq. (2), and show that our result also applies to this other model.
To better understand the dynamics of the generalized, D-dimensional Kuramoto model, we define an 'order parameter', ρ, that is analogous to the Kuramoto order parameter, N −1 j exp(iθ j ), used to analyze the system of Eq. (2) and is equivalent to it for D = 2: Like the Kuramoto order parameter, |ρ|= 1 corresponds to the system being a completely coherent state, σ i = σ j for all i, j; while |ρ|= 0 corresponds to an incoherent state. Using this order parameter, we can rewrite Eq.
(2) as It can be seen in the above equation that the dynamics of each agent is determined by the two terms on the righthand side. The first term (i.e., the term proportional to K) represents a global coupling of each agent to all the other agents through the order parameter. For positive K, this term attracts the state of each agent, σ i , towards the average orientation of the full population, characterized by the direction of ρ; whereas for negative K, this term causes dispersal of the system agents away from coherence, with each agent moving away from the average orientation of the agents. The second term also gives dispersing dynamics, with each individual agent having distinct individual dynamics when uncoupled from the other agents in the system. To completely specify the setup of the system, we need to specify the choice of the N natural rotations in Eq. (6). In the case of the standard D = 2 Kuramoto model, Eq. (1), the natural rotations are added in the form of individual distinct natural frequencies ω i for each individual agent. Assuming that the natural frequency of each agent is independently picked randomly according to a fixed unimodal distribution g(ω), the change in coordinates, θ i → θ i + ω 0 t, effectively reduces the natural frequency of each agent by any constant ω 0 . Thus the mean of the distribution g(ω) can be set to 0 without loss of generality. In the unit vector formulation of the D = 2 Kuramoto model, this is equivalent to the change of variables σ i → e W0t σ i , where The new equation after the change of variables has the rotation matrix shifted as In the case of D = 2, the matrices e W0t and W i commute, and hence the change is equivalent to the timeindependent transformation W i → W i − W 0 , allowing us to shift the mean of the distribution to any arbitrary quantity. For D > 2, however, commutation of antisymmetric matrices or rotation matrices does not generally apply (i.e., the rotation group in D > 2 is nonabelian), and hence this change of coordinates does not yield an equivalent model with a change of rotation matrices. Thus for D > 2 the mean of this distribution cannot be simply shifted as in D = 2. In general, for D dimensions, we specify the distribution over the space of antisymmetric matrices that we use to choose the individual W i for each agent i. We denote this distribution by G(W), which is analogous to the distribution g(ω) in the case of D = 2. In this paper, we restrict the choice of G(W) as follows: we choose each of the upper-triangular elements of W independently from a normal distribution with zero mean and a standard deviation of ∆. The lower-triangular elements are then set according to the constraint that W is an antisymmetric matrix. This particular choice of G(W) corresponds to an ensemble of antisymmetric matrices that has zero mean, and is invariant to rotations (choosing an anisotropic distribution, such as shifting the mean of the upper-triangular elements, or choosing the uppertriangular elements from normal distributions with unequal variance does not appear to change the qualitative results illustrated in Fig. 1). Hence, due to the rotational symmetry, |ρ|= 0 will be a solution to our system (note that this solution may be stable or unstable). Further, we also note that Eq. (2) is invariant to the transformation t → ∆ × t, K → K/∆ and W → W/∆, and hence, without loss of generality, we set ∆ to be unity for the remainder of this paper.
For future reference, it is useful to point out the following facts that apply to any real antisymmetric matrix A (such as W i ): (i) Since iA is Hermitian, the real part of all the eigenvalues of A is zero. Hence all nonzero eigenvalues will be purely imaginary or zero.
(ii) If λ is an eigenvalue of A, then so is −λ.
(iii) If D is odd, then A must have at least one zero eigenvalue [implied by (ii)]. Further, the corresponding eigenvector is real.
For example, following Eq. (2) we have noted that for D = 3 we can express W j σ j in the form ω j × σ j , with ω j = ω jωj . In terms of the above discussion,ω j is the real eigenvector corresponding to the zero eigenvalue of the 3 × 3 matrix W j (W jωj = 0), and ±iω j are the nonzero eigenvalues of W j . As discussed earlier, for D = 3 we can now represent the second term on the right-hand side of Eq. (6) as a cross product, giving Given the choice of the distribution G(W) made above, we can write the distribution of the natural rotations of individual agents as G The distribution of rotation directions, U (ω) is then isotropic, and independent of the distribution of rotation magnitudes, and the distribution of magnitudes is . This choice of the distribution G(ω) sets the mean of the distribution to always be 0. In numerically simulating this system, we observe that the order parameter always goes to a fixed point, similar to the case of the standard Kuramoto model with zero mean of the distribution of frequencies.

III. DYNAMICS AND EQUILIBRIA
To map out the interplay between the tendency to align and the natural rotation of the individual units [i.e. the two opposing tendencies represented by the two terms on the right-hand side of Eq. (2)], we plot numerically obtained phase transition diagrams for D = 2 -9 (see Fig. 4). For large N and varying values of the coupling strength K, we allow the system to reach its time asymptotic equilibrium, and then we plot the magnitude of the order parameter at equilibrium as a function of K. We note that the results in Fig. 4 apply for all the random initial realizations of the distributions of the individual states σ i that we have tested.
As would be expected from the earlier discussion, in the case of negative coupling, i.e., K < 0, the system of agents goes to a state which is incoherent, |ρ|≈ 0. For even D, as in the D = 2 Kuramoto model (Fig. 4), there exists a positive critical coupling constant K c > 0. In contrast, for odd D, coherence begins at K = 0, i.e., K c = 0. Moreover, in contrast to the even D case where the transition is continuous ('second order'), for odd D the transition is a discontinuous jump from |ρ|= 0 in K < 0 to |ρ|> 0 for K → 0 + , past which |ρ| increases continuously with increasing K, asymptoting at |ρ|= 1 as K → ∞. For example, for D = 3 we find that |ρ|= 0.5 at K = 0 + , and this result (as we shall subsequently show) is independent of the distribution g(ω). Furthermore, we find that this discontinuous transition is nonhysteretic. To better understand these observed phenomena, we now present a mathematical analysis of this system.

A. Coherent states for D = 3
We first focus on the case of a positive coupling constant K in three dimensions. We seek fixed points of the order parameter. To study these analytically we first solve for fixed points of the agents, assuming that the order parameter is at a fixed point with positive magnitude. We hence solve for σ F i .The superscript F indicates that the agent is at a fixed point. Given a spherically symmetric distribution of rotation vectors, we can choose the direction of the order parameter ρ arbitrarily. The magnitude of the order parameter must be chosen to be self consistent given the orientation of the agents, according to Eq. (5). We define a quantity µ i = ω i /(K|ρ|) to rewrite the above equation as whereρ = ρ/|ρ| is a unit vector in the direction of ρ. This vector equation can be solved (see Appendix A) to obtain where t i =ρ · σ F i , and ξ i =ρ ·ω i /ρ · σ F i . From Eq. (10) we observe that there are two fixed points for each agent, one in the same hemisphere as the order parameter vector (corresponding toρ · σ F i > 0), and the other in the opposite hemisphere (corresponding toρ · σ F i < 0). Importantly, we also observe that there is a fixed point solution σ F i for any given ω i , ρ and K. Do these solutions correspond to a stable or unstable fixed points? Given a steady-state solution with all agents at one of their fixed points, for some ρ such that |ρ|> 0,  42) and (44) for the critical coupling strength for even dimensions has been shown in correspondingly colored arrows on the x-axis in (a). For a discussion on the slight mismatch between the theory and the numerical results, see Sec. III C. We expect this mismatch to decrease with increasing N . The theoretical estimates from Eq. (13) for the magnitude of the discontinuity, i.e., |ρ| at K → 0 + are shown in correspondingly colored arrows on the y-axis in (b). Note the close match between the theoretical result for the discontinuity, and the numerical observation for |ρ| at K = 0.2.
we consider a perturbation i to the i th agent. Assuming where we have used i · σ F i = 0 so that the perturbed σ i remains a unit vector. Thus we see that the stability of the fixed point σ F i depends on the sign of ρ · σ F i , with positive (negative) ρ · σ F i implying a stable (an unstable) fixed point. Since for each agent σ i there are two solutions for σ F i with opposing signs of ρ · σ F i according to Eq. (10), each agent has a stable fixed point and an unstable fixed point. We assume that each agent will approach its stable fixed point.
This behavior is in contrast to the two-dimensional Kuramoto model, where the proportion of agents in the entrained population increases continuously from 0 as we increase K beyond K c . (This fundamental difference is due to the previously noted fact that W i for odd D always has zero as one of its eigenvalues.) To understand the presence of the discontinuous phase transition, we first look at the case of small coupling, such that 0 < K ∆. By ignoring the first term on the right-hand side of Eq. (8), or by considering the limit of µ i → ∞ in Eqs. (10) and (11), we see that σ F i = ±ω i . Since the stable fixed points corresponds to ρ · σ F i > 0, each agent will go to a stable fixed point given by [sgn(ρ ·ω i )]ω i . Note that this location of the fixed point on the unit sphere is independent of the magnitude of the agent's rotation vector, and depends only on the orientation of the rotation vector. Since the distribution of rotation vectors was chosen such that the distribution of directions U (ω) was uniform on the unit sphere, the fixed points sgn(ρ ·ω i )ω i will be a set of uniformly distributed points over the hemisphere, ρ · σ > 0, of unit radius. This is demonstrated in Fig. 5, where we illustrate the orientations of N = 5 · 10 3 agents at a fixed time for a coupling strength K = 0.1. In this plot, we have mapped the endpoints of the orientation vectors σ i on the unit sphere S to a rectangle via an area-preserving transformation (see Fig. 5 caption for details). At the initial time, corresponding to an initial uniform distribution on S, the agents are uniformly distributed on the rectangle, whereas after T = 50000 time units it can be seen that the agents are uniformly distributed over only the upper half of the rectangle, corresponding to the hemisphere ρ · σ > 0 of S.
As discussed earlier, the magnitude of the order parameter must be consistent with the orientations of the agents, according to Eq. (5). Thus, being the average of the orientations of all the agents, the order parameter will have |ρ|= 1/2, since the centroid of a hemisphere is located at a distance of half of the radius from the center of the sphere.
This result is independent of the choice of the distribution g(ω), provided the rotation vector directions are isotropically distributed. As discussed earlier, negative values of coupling result in the system going to an incoherent state, with |ρ|= 0, while here we see that for small positive coupling the order parameter attains a value of FIG. 5. Orientations of each of the N = 5 · 10 3 agents at time T = 0, T = 500 and T = 50000. We visualize the orientations by plotting the endpoints of the orientation vectors on the sphere S. The sphere is then mapped onto a rectangle using an area-preserving transformation. We choose the z-axis along ρ, and arbitrarily choose mutually orthogonal x and y axes. θ then represents the angle measured from the z-axis (cos θ =ρ · σ), and φ represents the azimuthal angle measured anti-clockwise from the x-axis. The agent state vectors are initialized with a uniform distribution on the sphere, and evolved with a coupling constant K = 0.1. Note how all of the agents tend to uniformly distribute themselves on one hemisphere. |ρ|= 0.5. This result naturally generalizes to higher odd dimensions. As for the case of D = 3, letω i be the real eigenvector corresponding to the zero eigenvalue of the D × D matrix W i . In the limit of 0 < K ∆, we can again ignore the first term on the right-hand of Eq. (2). Solving for fixed points, we set dσ i /dt = 0, and hence the fixed point solutions will be given by Following the same analysis as performed above for D = 3, we reach the conclusion that for small positive K, the agents will go to fixed points given by sgn(ρ ·ω i )ω i . Hence the magnitude of ρ at K = 0 + will be given by the position of the centroid of a uniform hemisphere in D dimensions: This matches well with numerical results shown in Fig.  4 (b), where the colored arrows indicate the theory predictions according to Eq. (13). Note the close agreement with the prediction of the magnitude of ρ indicated by these arrows at K = 0, and the K > 0 start of the phase transition curves at K = 0.2 shown by the various colored symbols. By setting up a consistency relation in a similar fashion (see below), we can calculate the magnitude of the order parameter as a function of the coupling constant for D = 3. As shown in Fig. 6, this theory (solid black curve) agrees well with results from simulations of Eq. (2) with N = 10 4 (red plus signs).
We now give our analysis for D = 3 resulting in the solid curve in Fig. 6. As earlier, we assume that ρ is in some particular fixed direction. Since the distribution of the direction of the unit vectors σ i has been taken to be isotropic, we can assume that there will be rotational symmetry of the distribution of stable fixed points of the agents about the axis along ρ. Thus Since each agent has a unique natural rotation vector, we label the agent state variables as functions of their rotation vectors, as opposed to the index label i. Since the rotation vectors are chosen from a given distribution G(ω), we can approximate the above sum as which applies in the limit N → ∞ in Eq. (2). We interpretρ · σ F (ω) as cos[θ(ω)], where θ(ω) is the angle between the direction of the order parameter, and the stable fixed point of the agent with rotation vector ω. We will later use Eq. (10) to insert the expression for cos[θ(ω)]. We write the above as |ρ|= cos[θ(ω 0 ,ω)]g(ω 0 )U (ω)dω 0 dω.
Performing a change of variables from ω 0 to µ = ω/(K|ρ|) we get |ρ|= cos[θ(µ,ω)]g(µK|ρ|)U (ω)K|ρ|dµdω, and hence This can now be numerically solved to obtain |ρ| for a given K. For example, for the particular choice of G(ω) discussed above, where the three components of each vector ω i are chosen independently from a normal distribution centered at 0 with a standard deviation of ∆, the integral in Eq. (18) overω can be split into an azimuthal integral about the axisρ, which is trivial, and an integral over the angle betweenρ andω, i.e., the ζ integral below where the integration variable ζ representsρ ·ω. Solving this integral equation numerically for |ρ| for different values of the coupling constant K we obtain the solid black curve in Fig. 6.
To complete the analysis of the coherent states for K > 0, we now discuss why the state vectors σ i approach their stable fixed points σ F i . We demonstrate this in the limit of 0 < K ∆. Under this assumption, we note that in Eq. (7), the typical magnitude of the second term on the right-hand side, O(∆), is much larger than the first term, which is O(|Kρ|). We refer to O(∆) as the fast time-scale, and O(|Kρ|) as the slow time-scale. The assumed separation of time-scales implies that, to lowest order, we can neglect the first term in Eq. (7), leading to the equation This has the solution depicted in Fig. 3, where the state vector σ i uniformly precesses rapidly aboutω i , with the quantity z i (t) = σ i (t)·ω i constant on the fast time-scale. To determine the dynamics over the slow time-scale, we consider the dot product of Eq. (7) withω i , and average both sides of the equation over the fast time scale. This gives the evolution of z i as where ρ = N −1 z iωi , with the angle brackets representing averaging over the fast time-scale. This equation has a single stable fixed point at +1 or −1 dependent on the sign of ρ ·ω. Thus starting from random initial conditions, z i (t) will move to its fixed point at sgn( ρ ·ω). This is equivalent to stating that each σ i will move to its fixed point [sgn(ρ ·ω i )]ω i . While we have only thus proved that σ i will approach σ F i in the limit of small K, we numerically observe this to be true for all K, i.e., each agent goes to its corresponding stable fixed point as discussed above.
Until now, we have restricted our discussion of coherent states to K > 0. Are there any stable coherent states for K < 0? If the answer were yes, the fixed points of the agents could be calculated as earlier resulting in Eqs. (10) and (11), and would be governed by the stability equation given in Eq. (12). Since K < 0, the stable fixed points will correspond to solutions where ρ · σ F < 0 for each of the agents. This would imply that all of the agents would point to the hemisphere that the vector ρ points away from (not toward), contradicting the definition of ρ as the average of the orientations of all the agents. Thus there cannot be any stable fixed point solutions with positive magnitude of the order parameter for negative coupling.
This, however, does not rule out the possibility of unstable coherent states with K < 0. Going back to Eq. (7), we make a few observations. First, since all natural rotations were chosen such that the distribution of rotation directions was uniform on the sphere, the transformation ω → −ω does not affect the distribution or the macroscopic dynamics of the agents. After this transformation, we note that transforming K → −K, and changing the direction of time, i.e., t → −t, leaves Eq. (7) invariant. Thus, each stable fixed point of the macroscopic order parameter, ρ, for a given value of coupling strength K > 0, is also a fixed point at a coupling strength of −K, but is unstable (since we have reversed the sign of time). Thus the curve of coherent stable states for K > 0 extends symmetrically to K < 0 representing coherent unstable states. These stable (solid black curve) and unstable states (dashed black curve) are shown in Fig. 6. We call these coherent states the 'upper' branch of the phase transition diagram.

B. Incoherent states for D = 3
When the order parameter has zero magnitude, the system is said to be incoherent. As we demonstrate, this state is stable for negative values of the coupling constant and unstable for positive vales of the coupling constant.
In order to address the incoherent state, we first consider the following question: Given a state where |ρ|= 0 for all time, what are the possible dynamics of the individual agents? Setting ρ = 0 in Eq. (2), we get dσ i /dt = W i σ i . In the case D = 3, this means that the state σ i of each agent precesses about their own rotation axes, as illustrated in Fig. 3. If each agent were randomly placed uniformly on S, then this would be consistent with |ρ|= 0, and would be a steady state. However, this is not the only such arrangement of σ i that is possible corresponding to |ρ|= 0. For example, if each agent, σ i was placed on the axis of the corresponding rotation vector, such that σ i =ω i (or σ i = −ω i ), then this would also be consistent with |ρ|= 0 [since we have assumed that U (ω) is uniform], and the agents would each be at fixed points (this will be possible whenever D is odd). In fact, the steady-state |ρ|= 0 applies for any random proportion p of agents oriented parallel to the axes of their natural rotations, and the remaining agents at uniformly distributed locations on the sphere. Thus for N → ∞ there are an infinite number of distributions of ω and σ for which |ρ|= 0 is a steady state.
To characterize these states in the limit of N → ∞, we assume that the distribution of agent orientations σ rendered onto the unit sphere S is well defined. We denote by F (σ, ω, t) the distribution of agents on S, such that F (σ, ω, t)dσdω is the fraction of agents that lie in the two-dimensional differential element on the surface S centered at σ at time t, and have a natural rotation vector within the differential element dω centered at ω. Since the natural rotations of each agent are time independent and are independent of σ, we can write where is the distribution of the antisymmetric natural rotation vectors, G(ω)dω = 1, and G(ω) = g(ω)U (ω). In terms of this distribution function F , the order parameter will be given as An example of a class of distributions in D = 3 for which |ρ|= 0 is a steady state is given by for any p ∈ [0, 1], where δ(·) represents the Dirac delta function. As we will demonstrate shortly, in the limit N → ∞, this entire class of distributions is stable to small perturbations for all K < 0, i.e., for the incoherent region demonstrated in Fig. 1. This is in sharp contrast to the case of D = 2, wherein there is a single stable incoherent steady-state distribution in the large system size limit (corresponding to f = 1/(2π)) for the incoherent region in Fig. 1.
However, we observe from numerical simulations with K < 0 (done at large, but necessarily finite N ) that, starting with an initial condition corresponding to Eq. (24) with p = 0 (i.e., with σ i distributed isotropically and independently of its corresponding ω i , for all i) we observe that σ i evolves slowly with time to either σ i = +ω i or σ i = −ω i (i.e., σ i aligns with its rotation vector), with about half of the population {σ i } going to +ω i , and half to −ω i . Furthermore, as N increases, the rate of this relaxation becomes slower and slower, approaching zero as N → ∞. In addition, the fractions of agents going to +ω i and −ω i approach 1/2 as N → ∞. Thus, taking the limit t → ∞ followed by taking the limit N → ∞, Eq. (24) with p = 1 (i.e., F (σ, ω) = g(ω)U (ω)[δ(σ −ω) + δ(σ +ω)]/2) appears to approximate the distribution of agents on S. If the order in which the limits are taken is reversed, then p = 0, its initial value (i.e., F (σ, ω) = g(ω)U (ω)U (σ)) represents the distribution of agents on S. Similar results apply for other odd values of the dimension D, whereω i is now the D-dimensional eigenvector of W i having zero eigenvalue and with magnitude one (i.e., W iωi = 0).
We illustrate these numerical results in Figs. 7 where we show the histograms of the initial (plotted in blue) and final (plotted in red) distributions ofω i · σ i over the N agents. These numerical simulations were performed with N = 1000, K = −2, ∆ = 1. In the insets we plot the time-series of z i vs time for 50 randomly chosen agents. We see that for all odd D,ω i · σ i evolves towards ±1. Note that a similar consideration of even D is inapplicable since a randomly chosen even-dimensional W i typically does not have a zero eigenvalue, and thuŝ ω i does not exist.
While macroscopically, in terms of the magnitude of the order parameter (|ρ|= 0), the N → ∞ stationary states with distributions given by Eq. (24) appear identical for all p, their stability to perturbations depends on p. To analyze the stability of this class of N → ∞ stationary states we perform a linear analysis. To do this, we first describe the dynamics of the system in terms of the distribution F . We treat Eq. (7) as a velocity field for the flow of this distribution and hence set up a continuity equation: with a velocity field v given by where ∇ S · A represents the operator for the divergence of an arbitrary vector field A, along the surface S of the unit sphere in σ-space. The order parameter, ρ is described in terms of the distribution function F according to Eq. (23). We show in Appendix B that the continuity equation Eq. (25) can be rewritten as where ∇ S Φ is the component of the gradient of a scalar field Φ that is parallel to the surface S. We consider a small perturbation, such that the distribution f (σ, ω, t) can be written as where ξ(σ, ω) is small. Inserting Eq. (28) into Eq. (27) and linearizing gives To further simplify this equation, we make a choice of basis, such that ω = ωẑ. This allows us to rewrite the above equation as sξ(σ, ω, t) + ω ∂ ∂φ ξ(σ, ω, t) = 2K(ρ · σ)f 0 (σ, ω), (30) where φ is the azimuthal coordinate around the z-axis.
In this basis, we can then write f 0 as where θ is the angle measured from the z-axis, and together θ and φ represent σ.
Inserting the form f 0 from Eq. (31) into Eq. (30), we then solve for ξ(σ, ω, t) and insert the obtained solution into Eq. (23) to obtain giving the final dispersion relation, Note that the case of p = 0 (corresponding to an initial condition with independently chosen, uniformly random σ) and the case of p = 1 (corresponding to an initial condition with each σ being eitherω or −ω) have different dispersion relations. Thus, despite having the same macroscopic characteristic of |ρ|= 0, they will have different stabilities to perturbation. In the limit of small K, s will also be small, and we can ignore the second term in the square brackets in the above expression. Thus, Note that since K is small, this represents the behavior of s for K around zero. Since s ∝ K, we see that the incoherent state, having |ρ|= 0, will be stable (s < 0) for K < 0, and unstable (s > 0) for K > 0 as has been represented in Fig. 6. We call these incoherent states the 'lower' branch of the phase transition diagram. It can be seen from Fig. 6 that the upper branch is stable whenever the lower branch is unstable (i.e., for K > 0), and the upper branch is unstable whenever the lower branch is stable (i.e., for K < 0). Thus, for no value of K are there two values of |ρ| that are stable. This lack of bistability implies that the transition from incoherence to partial coherence occurs nonhysteretically at K = 0.

C. Phase transition in even dimensions
So far our primary focus has been on the cases of odd dimensions. As discussed earlier in Sec. III and in Fig.  4, the even-dimensional cases exhibit continuous ('second order') phase transitions at positive critical coupling strength K c > 0. To better understand this, we extend the treatment of the D = 2 case (e.g., Ref. [29]) to even D > 2. Like in the case of D = 3, we assume that the system has reached an equilibrium, with the order parameter having a magnitude |ρ|. Unlike Eq. (11), wherein a fixed point for each agent exists for all values of W, this will no longer be the case for even D. Rather, only certain values of the natural rotation W will permit the existence of fixed points above a certain value of K. Similar to Ref. [29], we first determine the conditions on W that permit fixed points of the corresponding agents, and then use this to set up a consistency relation similar to Eq. (14) to determine K c . A key assumption in this approach is that for steady states with |ρ|> 0 with N → ∞, only agents for which σ is at a fixed point contribute to the sum in Eq. (5), which we prove a posteriori [see Eqs. (39) and (40) and accompanying discussion].
As in Eq. (8), we see that the fixed points of σ i must satisfy where we have dropped the index i for simplicity. Denoting the term (ρ · σ F ) as γ, we observe where 1 denotes the D-dimensional identity matrix. Since We now transform the above equation to a basis that block-diagonalizes the antisymmetric matrix W. There exists a real orthogonal matrix, R such that R T WR is a block-diagonal matrix whose j th block is the 2 × 2 matrix for all j ∈ {1, 2, . . . , D/2}. We will refer to these ω j as the Λ = D/2 natural frequencies associated with W. Further, we define ρ 2 k to be the sum of the squares of the magnitudes of the 2k − 1 th and 2k th components of Rρ. Then Eq. (36) can be simplified to Note that this change of basis does not affect the value of γ, since it is a scalar quantity. Each term in the summand of the above expression can be interpreted as being proportional to a Lorentzian function of γ centered about γ = 0, and hence has a single maximum at γ = 0. Thus, H(γ) will also have a single maximum at γ = 0, from which it follows that in order for Eq. (37) to have a real solution for γ, H(γ = 0) must be greater than or equal to 1. Hence the condition on W that will permit the existence of σ F will be For the case of the standard D = 2 Kuramoto model, the above criteria reduces to |ω|< |Kρ| (Ref. [29], Eq. (4.2)). For a given ρ, we denote the region in W-space that satisfies that the above criteria as Γ. Each W i ∈ Γ will have a corresponding fixed point for σ i and the set of such agents i will be referred to as the entrained population.
For each W j / ∈ Γ, σ j is continually in motion, and we refer to these agents as the drifting population. We now argue that the contribution to the order parameter, ρ from the drifting population will be zero, and then use the Eq. (37) to write out a consistency relation for the order parameter as calculated only from the remaining entrained population.
Assuming an equilibrium of the system, such that the order parameter is at a fixed point, the drifting agents must form a stationary distribution on S. We denote this distribution by f (σ, W), which is analogous to f (σ, ω, t) defined in Eq. (22). Since the velocity of each agent is governed by Eq. (6), stationarity of the distribution requires that f (σ, W) is inversely proportional to the magnitude of this velocity. Hence where C(W, Kρ) is a normalization constant, for each W not in Γ. Since Γ is invariant to the transformation W → −W, it follows from the definition of C(W, Kρ) that it must also be invariant to W → −W.
The contribution to the order parameter from the drifting population is then given by Applying the variable transformations of σ → −σ and W → −W we obtain ρ drift = −ρ drift , and hence |ρ drift |= 0. Thus, the only contribution to the order parameter is from the entrained population of agents. Let H(γ) = 1 give rise to some solution (ρ · σ F ) = γ ≡ γ({ω i }, {ρ i }). Then, dotting both sides of Eq. (5) with ρ in the limit of infinite system size gives As in the two-dimensional case, the critical coupling strength, K c , will be such that the magnitude of the order parameter is infinitesimally small but nonzero. We can use this to determine a value of the critical coupling as wherẽ and g(ω 1 , . . . , ω Λ ) is the joint distribution of natural frequencies associated with the distribution W (see Appendix C for details). Note that, for our particular choice of an antisymmetric matrix ensemble from which we randomly draw the W i (i.e., independently Gaussian uppertriangular matrix elements), there are known results for g andg from random matrix theory. In particular, Ref.
[36] yields[37] The predictions for the critical coupling strength, K c , made according to Eqs. (42) and (44) for D = 2, 4, 6 and 8 have been marked by vertical arrows in Fig. 4(a). We expect that with increasing N the numerically observed transitions will appear to be sharper at the marked critical coupling strength. Note that continuing the curve from large values of |ρ| to the x-axis without changing its curvature (as would be expected from the shape of the phase transition curve in D = 2; see Fig. 1) approximates the predicted values accurately.

IV. MODEL VARIANT: EXTENDED-BODY AGENTS IN THREE DIMENSIONS
From Eq. (2), the dynamics of the system of agents can be thought of resulting from the interplay of two terms, K[ρ − (ρ · σ i )σ i ], promoting coherence among agents, and W i σ i , promoting decoherence between agents. We have shown that the competition between these two opposing tendencies is resolved by a critical transition from incoherence to coherence that is qualitatively different for even and odd dimensionality (Figs. 1 and 4). In order to show that this qualitative result is not restricted to our particular assumed form of the K = 0 agent dynamics (dσ i /dt = W i σ i ), we here consider a very different model with D = 3, and show that our conclusion for the behavior shown for the solution of Eq. (2) continues to apply. Specifically, we consider a different form of the dispersal term in the context of the three-dimensional dynamics of extended objects (e.g., the fish in Fig. 8). We will also further justify the term W i σ i as a simple choice of dispersive dynamics for interacting agents.
As we discuss in Sec. II, the setup of Eq. (2) considers interactions between agents that are fully described by a single D-dimensional unit vector. For an extended object, a single unit vector does not uniquely specify the agent state. In the specific context of three-dimensional extended objects in three-dimensional space (e.g., the dynamics of flocks of birds, swarms of drones etc.), the orientation of the extended body must be specified by two unit vectors. We call such agents extended-body agents, and describe their state via the two vectors σ, which as earlier represents the direction of the velocity of the extended-body agent; and η, chosen orthogonal to σ (see Fig. 8). For simplicity, we define ν = σ × η to form the right-handed orthonormal triple {σ, η, ν}. It should be noted that extended-body agents in two dimensions are completely described by a single unit vector, σ, as in the (2), we assume that the state of an extended-body agent cannot be described by a single unit vector σ. Rather, the pair of vectors σ and η together describe the orientation and state of the agent. The direction of agent velocity is assumed to be along the direction σ as earlier. The unit vector ν is defined as ν = σ × η standard D = 2 Kuramoto model. We will first set up the dynamics of this extended-body agent when it is not coupled to other agents. Motivated by the uncoupled dynamics of this extended-body agent representing some fixed errors/biases, we assume that the uncoupled dynamics of this extended-body agent is autonomous, i.e., not explicitly dependent on time. Under this assumption, we write dσ/dt = Φ(σ, η), dη/dt = Ψ(σ, η), Further, we make the natural assumption that the dynamics do not depend on information of its orientation with respect to any fixed frame of reference, i.e., there is no 'special' direction in space that determines the dynamics of the extended-body agent. Thus, for any rotation matrix R. Since the unit vectors {σ, η, ν} form an orthonormal basis, we can write the vector field Φ in this basis, and hence the scalar a(σ, η) must be independent of σ and η, a(σ, η) = a. Similarly, b(σ, η) and c(σ, η) must also be independent of σ and η. Applying similar reasoning to all the components of Φ(σ, η), Ψ(σ, η) and Θ(σ, η), we see that they must each be linear functions of σ, η and ν. Hence [noting that Rν = (Rσ) × (Rη)] Further, since {σ, η, ν} are unit vectors forming a right-handed triple, and, using d/dt(σ · η) = d/dt(σ · ν) = d/dt(η · ν) = 0, Thus, Eqs. (52) reduce to for some scalar, extended-body-agent specific quantities α, β and γ. Having specified the uncoupled dynamics of an extended-body agent, we add the effect of inter-agent coupling, in the form of the Kuramoto-like interactions described in Sec. II. Thus, analogous to Eq. (2), we write where ρ is given by Eq. (5), which is the average of the velocity directions σ i of each extended-body agent. Note that this form of coupling treats σ as a special direction as compared to η and ν, since we assume that the goal of the swarm is to maintain coherence via coupling that aligns the velocity direction σ i of each agent i to the the motion of the swarm as a whole. We then writeη i andν i such that the constraint Eqs. (53) and (54) continue to hold for the coupled system and that K = 0 corresponds to Eqs. (45).
We perform a simulation of N = 10 4 such extendedbody agents by numerically integrating Eqs. (56), (57) and (58) for a range of values of K similar to Fig. 4. Since the dynamics captured by the Eqs. (45) represent random biases/errors, we choose the quantities α, β and γ for each agent from independent, normal distributions with zero mean and unit variance. In Fig. 9 we present the phase transition displayed by this system of evolving extended-body agents. For each value of K we numerically integrate the system until |ρ| reaches a steady-state value. Note that we continue to observe a discontinuous transition of |ρ| as K increases through 0. Further, we also numerically observed that if α, β and γ are chosen anisotropically, i.e., if they are chosen from normal distributions with zero mean but differing variance, the qualitative result shown in Fig. 9 does not change, i.e., the transition to coherence is still discontinuous at K = 0. This indicates that the phenomenon of discontinuous transitions in odd dimensions is not specific to the form of the dispersal term chosen in Eq. (2), rather, it is a more general phenomena occurring for a potentially wide range of systems of interacting agents in odd dimensions. In contrast, this model for two dimensions (β = γ = 0) is the same as the original Kuramoto model and hence has a continuous transition to coherence at a critical positive value of K.
To further examine the dynamics of the uncoupled agents Eqs. (55), we adopt notation where we represent Eqs. (55) as where U is the 3 × 3 antisymmetric matrix, and σ η ν represents a 3 × 3 matrix whose columns are the vectors σ, η and ν. We can then consider a change of basis R such that where ω 2 = α 2 + β 2 + γ 2 . Using the same convention as Eq. (59), we define the orthonormal triple of unit vectors {u, v, w} as Thus Eq. (59) becomes According to Eqs. (63), the vectors u and v, which are fixed linear combinations of σ, η and ν, undergo uniform rotation with an angular frequency of ω about the axis w. Since σ, η and ν describe the physical orientation of the uncoupled extended-body agent, the extendedbody agent will demonstrate dynamics that correspond to rotations in three dimensions, and there will exist a W such thatσ = Wσ. In particular, the axis of this rotation will be along the unit vector w given by Note that R (and hence its components) is dependent on the random biases/systematic errors present, arising from the particular form of Φ and Ψ, whereas σ(0), η(0) and ν(0) = σ(0) × η(0) depends on the initial orientation/state of the extended-body agent. Thus the axis of rotation is dependent on the initial state of the extendedbody agent, while the frequency of rotation, ω is determined solely by the random systematic errors of the extended-body agent [i.e., α, β and γ in Eq. (55)].
Thus, under the assumptions made above, the dynamics of uncoupled extended-body agents can be described asσ = Wσ for some initial-condition-dependent W (in particular, Wσ = −ωw × σ). Note however that this is not identical to the uncoupled dynamics of the agents described in Eq. (2). In particular, the axis of rotation of the extended-body agent under the dynamics described here is along the vector w, which is determined by the initial conditions of the extended-body agent state (σ(0), η(0)). However, in the uncoupled dynamics of the generalized Kuramoto agents described by Eq. (2) the axis of rotation is predetermined by the rotation matrix W i assigned to agent i and is independent of the initial condition chosen for the agent. An isotropic ensemble of rotation matrices for the generalized Kuramoto agents in the case of extended-body agents corresponds to an 'isotropic' distribution of the extended-body agent parameters (α, β, γ), as well as isotropic initial conditions of the extended-body agents.
Further, this simple interpretation of σ undergoing uniform rotation no longer holds for the case of coupled extended-body agents, and Eqs. (56), (57) and (58) cannot be simply written in the form of Eq. (2) with an initial-condition-dependent W for arbitrary K (this, however, is possible in the limit of K → 0 or |ρ|→ 0, hence our results for the stability analysis of the |ρ|= 0 state will recreate the phase transitions in higher dimensions). The qualitative dynamics of coupled extendedbody agents and coupled generalized Kuramoto agents described by Eq. (2) are in general distinct, yet our main point of discontinuous phase transitions for odd dimensions at K = 0 continues to hold.
Thus, for the case of extended-body agents, under the assumptions made in this section, rotation matrices as dispersal terms arise naturally as simple error/fixed-bias terms for the individual agents. Rather than considering the case of initial-condition-dependent rotation matrices, in Eq. (2) we have considered the simplification of choosing fixed rotation matrices W. This motivates the generalization of the Kuramoto model presented in Eq. (2) as a simple model to capture the dynamics of swarming and flocking agents. Further, we also see that the result obtained from our toy model Eq. (2) for the qualitative continuous or discontinuous behavior of the incoherentto-coherent transition continues to hold for other threedimensional agent dynamics, such as the extended-body agent dynamics described in this section.
In Sec. V we briefly describe other extensions and variants to the generalized Kuramoto model described in Sec. II.

V. DISCUSSION AND CONCLUSIONS
We have considered a generalization of the Kuramoto model to arbitrary dimensions, describing a system of interacting, orientable units, whose state is completely described by D-dimensional unit vectors. Our main result (Fig. 4) is that the macroscopic dynamics of the Kuramoto model is strongly dependent on the dimensionality of the system, with odd-dimensional systems behaving similar to one other, and likewise for even-dimensional systems. For odd-dimensional systems, including the practically important case of D = 3, we find that the phase transition from incoherence to partially coherent states occurs via a discontinuous, nonhysteretic transition as the coupling strength K increases through 0 (Sec. III A, also see Fig. 6). In contrast, even-dimensional systems, like D = 2, numerically appear to undergo continuous transitions of the coherence at a critical coupling strength K c > 0 ( Fig. 4 (a)). We also note that, unlike the two-dimensional Kuramoto model, the state of the system is not always completely classified by the magnitude of the order parameter. In particular, for the twodimensional Kuramoto model there is a single stable incoherent steady-state distribution in the infinite size limit (f = 1/(2π)), whereas the three-dimensional Kuramoto model has an infinite number of such distributions (for example, Eq. (24)) each with different stability properties (see Eq. (32)). By considering a setup of extendedbody agents, in Sec. IV we further motivated our choice of model Eq. (2) in the context of swarms of drones or flocks of birds. In particular, we demonstrated that our qualitative results relating to the difference between oddand even-dimensional systems continue to hold for models that use a different choice of the dispersal term. This study of extended-body agents in Sec. IV also explains why the choice of the dispersal term Wσ in the context of the qualitative phase transitions observed for D = 3 is justified .
While other authors [26,27] have also studied the Kuramoto model generalized to higher dimensions, their consideration has been limited to the case of identical natural rotations. Our setup of the problem (i.e., with heterogeneous natural rotations) by setting G(W) = δ(W − W 0 ) reproduces the results in Refs. [26,27] for the case of globally coupled systems (here we interpret the Dirac delta function acting on the antisymmetric matrix W as the product of Dirac delta functions acting on each of the upper-triangular elements of the matrix individually). This heterogeneous setup of the problem now describes the interplay of two opposing tendencies, i.e., the tendency for agent states to align due to the interagent coupling, and the tendency for agents to disperse themselves in opposition to such alignment. This leads to the possibility of new and interesting phenomena such as the difference between the odd and even dimensionality described in this paper.
In addition to the variant described in Sec. IV, the setup of the generalized Kuramoto model given by Eqs. (6) and (5) can be modified and generalized in various ways. An interesting question for possible future study is whether a striking difference between odd and even dimensions (as we have found for the generalized Kuramoto model and its variant in Sec. IV) manifests in these modifications. For example, beyond the globally coupled systems we have considered, one might consider network-based coupling, wherein agent j influences agent i with a strength A ij . This is equivalent replacing ρ in Eq. (6) with ρ i , where In the context of swarms of drones, a further natural generalization would be to have the network-based coupling A ij depend on the spatial distance and relative orientation between the i th and j th swarm agent. As discussed earlier, for positive K the dynamics of each σ i are attracted towards the average state of the system, ρ. This could be interpreted as a target direction for each σ i , and can be generalized by replacing Eq. (5) by other definitions of ρ. For example, in the context of swarms of drones, it could be desirable for the orientation of the drones to be biased towards the plane of the horizon, or to be biased toward the direction of a given target destination. To achieve this, the 'target direction', ρ in Eq. (6) could be modified from the average state of the system to the average state biased towards a given target. Studying the dependence of the dynamics of such swarms of agents on modifications to ρ (via either the presence of network dependent interaction, or other bias targets) would be an interesting line of future research.
In a future paper [38] we will present a mathematical formulation for studying the D-dimensional Kuramoto model in the infinite size limit via a generalization of the Ott-Antonsen ansatz [30,39], wherein we will also address the issue of generalization of ρ.
Keepingω fixed, we can independently choose K, and hence µ. Thus the term in Eq. (A6) in the square brackets must be independently zero.
Taking motivation from the shape of the domain Γ shown in Fig. (10), we express Γ as the disjoint union of Γ 0 , Γ 1 , . . . , Γ Λ , where Γ 0 is the component of Γ within the dashed circle of radius L in Fig. 10, and for j ≥ 1, Γ j is the region for which |µ j | 1 and |µ k |≥ L for all k = j.
Note that the left-hand side of Eq. (41) is ρ 2 , hence we can ignore terms on the right-hand side of order smaller than O(ρ 2 ). We now show that the contribution from Γ 0 is of a smaller order than this. By construction, in the subdomain Γ 0 , |µ i |≤ L, and hence µ i Kρ i ∼ O( √ ρ) → 0 as ρ → 0. Further, γ = ρ·σ F ≤ ρ. Thus the contribution I Γ0 to the integral in Eq. (C3) from the subdomain Γ 0 will be Since L 1, the volume of Γ 0 will scale as O(L Λ−1 ) ∼ O(ρ −(Λ−1)/2 ). Thus I Γ0 ∼ O(ρ (Λ+3)/2 ), which for D > 2 is negligible compared to ρ 2 and the contributions to the integrals in Eq. (C3) from the subdomains Γ j . Since K c for D = 2 is already known (e.g. Ref. [29]), we focus on the cases D ≥ 4, and hence will ignore the contribution from the subdomain Γ 0 . By symmetry, each Γ i will give the same contribution. Hence, without loss of generality, we will look at the contribution from the subdomain Γ 1 , and will append a factor of Λ. We will also only look at µ i > 0 and will hence append a factor of 2 Λ .  In the subdomain Γ 1 , µ i µ 1 for all i ≥ 2, and thus we can use the above equation in the small ρ 1 approximation, Also, µ 1 < 1 implies µ 1 Kρ 1 ∼ O(ρ) → 0 as ρ → 0. Thus, We then change variables back to ω i for each of the µ i integrals for i = 2 . . . Λ, and explicitly evaluate the integral over µ 1 . The lower limits of the integrals change from L to LKρ i which goes to zero in the limit of small ρ, since L ∼ O(ρ −1/2 ). Thus Since we are in the limit of small but nonzero ρ, we can cancel ρ 2 from both sides to obtain the desired result in Eq. (42)