Complexity reduction ansatz for systems of interacting orientable agents: Beyond the Kuramoto model

Models of systems of many coupled dynamical agents are useful tools for studying a very wide variety of phenomena.¹ Examples include flashing fireflies,^2,3 circadian rhythms of mammals,^4,5 oscillating neutrinos,⁶ arrays of Josephson junctions,⁷ oscillation of footbridges,⁸ biochemical oscillators,^9,10 power-grids,^11,12 collections of neurons,^13–16 flocking dynamics,^17–20 and others. In many cases, the states of the individual agents can be described by a single angle-like variable, $θ$⁠. This class of model systems includes situations for which the dynamical agents are oscillators,¹ neurons,^13–16 or robots moving on a two-dimensional plane,¹⁸ among others. Many such models, possibly involving network-based interactions,^16,21,22 such as the Kuramoto model,²³ the Kuramoto–Sakaguchi model,^22,24 and models of theta neurons,¹³ among others,²⁵ reduce to the form

where $θi$ represents the state of the $ith$ agent, $ηi$ is a (possibly vector) constant parameter that is associated with the $ith$ agent, $ω(ηi,{θ},t)$ is its “natural frequency,” $N$ is the total number of agents, and $H(ηi,{θ},t)$ is a common field that acts on each agent, dependent on the agent’s parameter $ηi$⁠, and ${θ}$ indicates a dependence on the set of states ${θ1,…,θN}$ in the form of an average over $i$ of a function of the angle $θi$⁠.

For example, the well-studied Kuramoto model^23,26 can be expressed in the form of Eq. (1) by choosing $H(ηi,{θ},t)=N−1∑jexp⁡(ιθj)$⁠, independent of $ηi$⁠, and choosing $ω(ηi,{θ},t)$ independent of ${θ}$ and $t$⁠, allowing $ω(ηi,{θ},t)$ to be replaced by $ωi$⁠. Reference 27 introduced an ansatz to analytically achieve substantial reductions in the complexity of problems of the type exemplified by Eq. (1) in the limit of a large number of agents (⁠$N→∞$⁠). Subsequently, this reduction has been applied in studies of a wide variety of systems (e.g., Refs. 5,7,8,13–16,28, and 29).

Several flocking models employ the Kuramoto model (e.g., Refs. 17–19) to describe orientational alignment of the velocities of individuals in a flock. Since the standard Kuramoto model [in common with other models conforming to the general form of Eq. (1)] describes the dynamics of scalar angles, these models are restricted to describing flock dynamics in a two-dimensional plane. Other work has shown that the Kuramoto model can be generalized to flocks moving in three and higher-dimensional space.^30–35 In this case, each agent’s state is assumed to be specified by a unit vector $σi(t)$ in the $D$-dimensional space. Alternately, we may think of $σi$ as specifying a point on the unit sphere in $D$-dimensional space. Reference 30 notes that the vector $σi$ can be thought of as representing the opinion of an individual in a group or the orientation of the velocity of a member of a flock. (For the case of flocking of birds, fish, or flying drones, the generalization to $D=3$ is of most interest.) For $D=2$⁠, the unit vector $σi$ is determined by its scalar orientation angle $θi$ specifying a point on the unit circle, thus recovering the previous model, Eq. (1) (see Sec. II). References 36–39 have also studied the Kuramoto model and its generalizations to higher dimensions in the contexts of continuous-time consensus protocols, multiagent rendezvous, distributed control, and coalition formation. In this paper, we present a new technique that enables analytic treatment of the dynamics of a large class of systems with higher-dimensional agents, including the aforementioned systems. In particular, in this paper, we focus on the continuum limit of infinitely many higher-dimensional agents, allowing us to use ideas similar to those developed previously in the context of Eq. (1).

The remainder of the paper is organized as follows: In Sec. II, we construct a generalization of Eq. (1) to arbitrary dimensions and describe the infinite system size limit in such systems. Then, in Sec. III, we extend the ansatz of Ref. 27, resulting in a simplified analytic description of this generalized class of systems. In Sec. IV, we demonstrate the utility of our results to example systems, with particular focus on the Kuramoto model generalized to higher dimensions. Finally, in Sec. V, we conclude with a discussion and summary of our results.

In two recent papers,^40,41 we constructed a generalization of the Kuramoto model to $D$ dimensions. Here, we consider an even more general setup, where we consider a generalization to Eq. (1) to a system in $D$ dimensions,

where for each $i$⁠, $σi(t)$ is a real $D$-dimensional unit vector, $|σi(0)|=1$⁠, $ρ(ηi,{σ},t)$ is an arbitrary real $D$-dimensional vector, which can be thought of as a common field that affects each agent in an $ηi$ dependent fashion, $W(ηi,{σ},t)$ is a real $D×D$ antisymmetric matrix, $ηi$ is a (possibly vector) constant parameter associated with each agent, and, as earlier, ${σ}$ indicates a dependence on the set of all states ${σ1,…,σN}$ in the form of the average over $i$ of a function of the unit vectors $σi$ (we further quantify this dependence on ${σ}$ later). For example, in the context of flocking agents in $D$ dimensions, $σi$ represents the orientation of the $ith$ agent, $ρ(ηi,{σ},t)$ represents a “goal” orientation to which the $ith$ agent attempts to align itself, and $W(ηi,{σ},t)$ represents a fixed bias, or a systematic error to the agent dynamics causing the agent to head in a direction that deviates from the direction of $ρ$⁠.⁴⁰ Note from the form of Eq. (2) that the dot product of the right-hand side of Eq. (2) with $σi$ is identically zero so that $d|σ|/dt=0$⁠, as required by our identification of $σ$ as a unit vector. Thus, the dynamics of each $σi$ is restricted to the $(D−1)$-dimensional surface, $S$⁠, of the unit sphere, $|σ|=1$⁠. For $D=2$⁠, choosing $σi=(cos⁡θi,sin⁡θi)T$⁠, $ρ(ηi,{σ},t)=Re[H(ηi,{θ},t)],Im[H(ηi,{θ},t)]T$⁠, and

reduces Eq. (2) to Eq. (1), thus justifying Eq. (2) as a $D$-dimensional generalization of Eq. (1).

We now consider the limit of a large number of agents and denote by $F(σ,η,t)$ the distribution of agents on $S$⁠, such that $F(σ,η,t)dD−1σdη$ is the fraction of agents that lie in the $(D−1)$-dimensional differential element $dD−1σ$ on the surface $S$ centered at $σ$ at time $t$⁠, and have an associated parameter $η$ within the differential element $dη$ centered at $η$⁠. Since the associated parameter $η$ for each agent is time independent, we define

and

Noting that Eq. (2) specifies the vector field of the flow controlling the dynamics of the distribution $f$⁠, we write a continuity equation for $f$⁠,

where the velocity field $v(σ,η,t)$ is given by $v(σ,η,t)=(ρ(η,t)−(σ⋅ρ(η,t))σ)+W(η,t)σ$⁠, and $∇S⋅A$ represents the divergence of a vector field $A$⁠, along the surface $S$⁠. This can be done if the dependence of $ρ$ and $W$ on ${σ}$ can be specified as a functional of $F(σ,η,t)$ that is not explicitly dependent on $σ$⁠. [A simple example of such a dependence on ${σ}$ would be the average value of $p(σi)$ for some given function $p$⁠, which can be written as $∫∫p(σ)F(σ,η,t)dσdη$⁠.] Following Appendix B of Ref. 40, Eq. (3) can be rewritten as

where $∇SΦ$ is the gradient of a scalar field $Φ$ projected on the surface $S$⁠.

For $D=2$⁠, Refs. 25 and 27 demonstrated that the ansatz that $f(θ,t)$ is in the form

where $α(η,t)$ is a complex scalar function of $η$ and $t$⁠, $|α(η,0)|<1$⁠, reduces Eq. (1) to the following $θ$-independent form:

The form Eq. (5) represents an invariant manifold in the space of possible distributions $f$⁠, that satisfy the continuity equation Eq. (3) for $D=2$⁠. Furthermore, previous work^25,42 has shown that initial conditions for $f$ are attracted to the invariant manifold Eq. (5) for a large class of possible models of the form Eq. (1). Thus, Eq. (5) can be used to greatly simplify the study of the long-term dynamics of these systems.

Here, we present an ansatz demonstrating the existence of a similar invariant manifold for Eq. (3) in any dimension $D$⁠. Noting that $eιθ$ can be interpreted as a unit vector in the complex plane and that the complex quantity $α$ can similarly be interpreted as a two-dimensional vector of its real and imaginary parts, based on Eq. (5), we posit the following guess for the form of $f(σ,η,t)$ for arbitrary dimension $D$⁠:

where $α$ is a real $D$-dimensional vector such that $|α(η,0)|<1$⁠, $βD$ is a yet-to-be-determined constant, and $ND(α)$ is a scalar normalization chosen to ensure that

Inserting Eq. (7) into the continuity equation in Eq. (4), we obtain after some algebra,

For our ansatz Eq. (7) to apply, the above equation must hold for all $σ$⁠. Focusing on the term in Eq. (9) that is quadratic in $σ$⁠, i.e., $ND(α)[2(D−1)−βD](α⋅σ)(ρ⋅σ)$⁠, since, in general, $α$ and $ρ$ will not be zero for all $t$⁠, we require that

With $βD$ in Eq. (7) determined, we now obtain the normalization constant $ND(α)$⁠. To perform the integral in Eq. (8), without loss of generality we take the vector $α$ to be along the $z^$ axis. For an arbitrary point $σ$ on $S$⁠, we denote the angle between $σ$ and $z^$ by $θ$⁠. In particular, we note that the distance of the point $σ$ from the $z^$ axis is $sin⁡θ$⁠. For a coordinate system on the surface $S$⁠, we use $θ$ as one of the coordinates, denoting position with respect to $z^$ on the sphere. From the symmetry of $f$ in Eq. (7) about the direction $α$⁠, we see that the integrals over these remaining coordinates give the surface area $SD−1sinD−2⁡θ$ of the $(D−2)$ dimensional surface of a sphere with radius $sin⁡θ$ embedded in $(D−1)$ dimensions, where $SD−1=(2π)(D−1)/2/Γ((D−1)/2)$ is the area of the sphere of unit radius in $D−1$ dimensional space. Thus, Eq. (8) becomes

which can be evaluated to give

where $KD$ is a constant dependent only on $D$⁠. This results in

giving the form of the ansatz for arbitrary dimensions as

which, for $D=2$⁠, agrees with Eq. (5).

To determine whether the ansatz Eq. (14) is consistent with Eq. (9), we insert it into Eq. (9). We find that the ansatz with $βD$ given by Eq. (10) indeed is a solution of Eq. (9) and that Eq. (9) reduces to the following equation for $α$ (see the  Appendix for details):

The key point is that Eq. (15) does not involve $σ$ (and remarkably, also does not involve any dependence on $D$⁠). Thus, analogously to Eq. (6), we have a $σ$-independent description of the dynamics of $α$⁠. This is our main result.

We note that for initial conditions with $|α|<1$⁠, $|α|$ will remain less than 1 at all finite times since from Eq. (15) $∂t|α|=0$ at $|α|=1$⁠, thus verifying that $f$ given by Eq. (14) does not diverge for $t<∞$⁠.

We now consider a few examples illustrating the utility of the generalized ansatz, Eq. (14), to systems of the form given in Eq. (2). We detail the particular example of the Kuramoto model generalized to $D$ dimensions⁴⁰ as representative of the utility of our main result Eq. (15), and thereafter briefly mention applications of this result to a variety of other systems.

A generalization of the Kuramoto model with homogenous oscillators to arbitrary dimension was introduced by Olfati-Saber in 2006³⁰ in the context of flocking dynamics, consensus protocols, and opinion dynamics. This was later generalized to heterogeneous systems by Chandra et al.⁴⁰ For generalization to $D$ dimensions, a system order parameter, $z$⁠, can be defined as

The magnitude of $z(t)$ is a measure of the coherence of the set of agents ${σ}$⁠. The common field $ρ$ is then defined as the $ηi$-independent function,

where $K$ is a coupling constant. By interpreting the vector parameters $ηi$ in $W(ηi,{σ},t)$ as the $D(D−1)/2$ independent elements of a $D$-dimensional antisymmetric matrix $Wi$⁠, we can replace $g(η)dη$ in integrals by $G(W)dW$⁠, where $G(W)$ is a distribution of antisymmetric matrices. In cases such as these where $W(ηi,{σ},t)$ is independent of ${σ}$ and $t$⁠, we interpret $W(ηi)=Wi$ as the “natural rotation” of $σi$⁠.

In the limit of infinite system size, with a distribution of agents given according to Eq. (14),

For $D=2$ (i.e, the original Kuramoto model) Eq. (18) evaluates to give $ρ(t)=Kz(t)=K∫dωg(ω)α(ω,t)$⁠. Equation (15) is then equivalent to Eq. (6) from Ref. 27. For $D=3$⁠, the integral in Eq. (18) gives

This now allows us to use Eq. (15) with some given $G(W)$ to numerically integrate for the dynamics of $α$ and the dynamics of the order parameter $ρ$⁠.

Using this simplification, we can efficiently simulate the dynamics of the full system of agents governed by Eq. (2). We first focus on the case of homogenous agents, i.e., identical natural rotations for each agent, $G(W)=δ(W−W0)$⁠, where $δ(⋅)$ is the Dirac-delta function. We can then change to a rotating basis in which the natural rotation term of each agent is zero, $W0→0$⁠. This makes the $W$-integral in Eq. (18) trivial, allowing a direct representation of $ρ$ in terms of $α$⁠. Further, $α$ is only dependent on time (rather than $W$ and $t$⁠). This represents a very large simplification in the complexity of the dynamics of the system of agents, since Eq. (15) is now a single $D$-dimensional ordinary differential equation which represents the collective dynamics of the $N→∞$⁠, $D$-dimensional system of coupled differential equations in Eq. (2). The utility of this result is demonstrated for $D=3$ in Fig. 1(a), where we show (plotted in black) the time-series for $|ρ(t)|$ as generated from a system of $N=5000$ agents (approximating the $N→∞$ limit), compared with the time-series generated from the theory derived in Eq. (15) (orange dashed curve). The initial condition for the full system was chosen such that the agents were uniformly randomly distributed on the sphere. For the theory derived in Eq. (15), i.e., the reduced equations, the initial value of $α$ was chosen to have magnitude $0.01$ in an arbitrary direction. Note the remarkably close agreement between the black and the orange dashed curve, demonstrating that the dynamics on the reduced manifold of Eq. (14) indeed gives the large-$N$ dynamics of the full system of interacting agents.

For the case of heterogeneous agents, $α$ in Eq. (18) depends on $W$⁠, and we perform the integral in a Monte-Carlo fashion. We randomly choose $N{W}$ values of $W$ from the given distribution $G(W)$ and simulate the dynamics of the corresponding $α(W)$s. These randomly chosen $α(W)$s are then used as the Monte-Carlo samples to evaluate $z$ according to Eq. (18), simulating the dynamics of the system in the $N→∞$ limit by only simulating the dynamics of $N{W}$ variables. Here, we choose an isotropic distribution $G(W)$ constructed by choosing each upper triangular element from identical independent normal distributions with zero mean and unit variance, and choosing the remaining elements to make $W$ antisymmetric. Results are shown in Fig. 1(b) for $D=3$⁠, where $N{W}=500$ Monte-Carlo samples were chosen to evaluate the $|ρ(t)|$ curve via the theory in Eq. (15), and are compared with the curve obtained for simulating the dynamics of the full system of equations in the $N→∞$ limit, approximated by a simulation of $N=5000$ agents. Note how simulating the dynamics of $N{W}≪N$ Monte-Carlo samples yields a smooth curve approximating the noisy curve generated by simulating the individual dynamics of $N=5000$ agents. Initial conditions for the full system were chosen as a bimodal distribution of $σi$s, independent of the corresponding $Wi$⁠, with the two peaks being anti-podal to each other, hence representing a distribution explicitly not on the manifold dictated by Eq. (14). The initial condition for the reduced equations, i.e., Eq. (15) were chosen to be uniform on a sphere of radius $0.01$⁠, corresponding to an approximately uniform distribution of $f(σ,η,t)$ in $σ$⁠. Despite not lying on the invariant manifold described by Eq. (14), we observe that the dynamics of the full system rapidly approach the dynamics as predicted by Eq. (15) for the $N→∞$ limit for dynamics on the invariant manifold. This indicates that for the case of heterogeneous agents the invariant manifold Eq. (14) is attracting, as has been proven for the case of $D=2$⁠.²⁵ Full system simulations with initial conditions described by a uniform distribution in $σ$ [and hence lying on the invariant manifold Eq. (14) for $|α|=0$] yielded a curve that is not discernibly different from the curve presented in Fig. 1(b).

As demonstrated in Fig. 1, numerical integration of the dynamics on the invariant manifold, via Eq. (15) closely reproduces the time evolution of the order parameter of the Kuramoto model generalized to higher dimensions [Eq. (2) for $ρ$ according to Eq. (19)]. As we demonstrate in Fig. 2, this close similarity between a simulation of the full $N$-agent dynamics and the simulation of the reduced equation Eq. (15) holds at all values of the coupling constant $K$⁠. This allows us to recreate the discontinuous phase transition of the Kuramoto model generalized to 3 dimensions reported in Ref. 40. [The continuous phase transition observed through reduced equations of the form Eq. (15) for the standard Kuramoto model in two dimensions has been demonstrated and discussed in Ref. 27.]

Reference 41 demonstrated that the Kuramoto model generalized to even dimensions $D≥4$ exhibits the unusual behavior of Instability-Mediated Resetting: If the coupling strength $K$ is increased abruptly while remaining below the critical coupling strength for the onset to coherence, the Kuramoto model displays a short burst of coherence (see Ref. 41 for further details). This is illustrated by the results shown in Fig. 1(c). In Fig. 1(c), the initial conditions for the $σi$ in the full system simulation were chosen independently to be uniformly random over the sphere $|σ|=1$⁠. For the reduced equations, initial conditions for $α$ were chosen similar to the earlier discussion on heterogeneous agents, i.e., uniform on a sphere of radius $0.01$⁠, corresponding to an approximately uniform distribution of $f(σ,η,t)$ in $σ$⁠. Figure 1(c) demonstrates that the dynamics of Eq. (15), representing the dynamics on the invariant manifold described by Eq. (14) (shown in the dashed orange curve), accurately captures these short bursts of coherence, further demonstrating the capability of Eq. (15) in capturing the transient dynamics of the Kuramoto model generalized to higher dimensions.

To demonstrate the applicability of Eq. (15) in improving theoretical understanding of such systems, we present an example of a stability analysis of the Kuramoto model generalized to three dimensions via a study of the dynamics on the reduced manifold Eq. (14). In particular, we study the stability of the completely incoherent state for a system of heterogeneous agents, wherein the initial condition of each agent $σi$ is to be distributed independently and uniformly over the sphere $|σ|=1$⁠. In the $N→∞$ limit, this is the distribution $F(σ,η,t)=g(η)/(4π)$ corresponding to setting $|α|=0$ in Eq. (14). [In the context of the stability analysis presented in Sec. III B of Ref. 40, we are considering the case of $p=0$⁠.] Thus, we are interested in the stability analysis of Eq. (15) about $|α|=0$⁠. Performing a first order expansion of Eq. (15) for $|α|≪1$⁠,

Assuming $α(W,t)=α(W)est$⁠, we obtain

Note that in the limit of $|α|≪1$⁠, Eq. (19) can be written as

Multiplying Eq. (21) by $(4/3)G(W)$ and integrating, we obtain

Note that in three dimensions the linear transformation $Wσ$ can be represented as the cross product $ω×σ$⁠. Without loss of generality we may choose a basis that block-diagonalizes $W$⁠, corresponding to the choice of $ω=ωz^$ and

Thus,

This can now be inserted into Eq. (23) and written in a basis independent format as

We choose the distribution $G(ω)$ to be an “isotropic” distribution (i.e., a distribution that is invariant to orthogonal transformations) which can hence be written as $G(ω)dω=g(ω)U(ω^)dωdω^$⁠, where $U(ω^)=1/(4π)$ represents the isotropic distribution of rotation directions, and $g(ω)$ is the distribution of the magnitudes of rotation (see Ref. 40 for further discussion on the choice of this distribution and its implications). Integrating over the rotation directions $ω^$ in Eq. (24) gives us

which is identical to the result obtained in Ref. 40. As discussed in Ref. 40, the above equation implies that in the limit of small $K$⁠,

indicating that this completely incoherent state loses stability at $K=0$⁠. Thus, using Eq. (15) allows us to perform the stability analysis of a state easily without having to solve the partial differential equation of the dynamics of the distribution of agents in both $σ$ and $t$ as was necessary in Ref. 40.

Extensions appropriate to various contexts may be studied using Eq. (14).

For example, each of the agents in the model described above could have a bias toward a particular subspace, such as birds in a flock that have a preference to align parallel to the surface of the Earth. In this case, the common field of such a system is then defined similar to Eq. (17) as

where $Π$ is the operator that projects onto the preferred subspace (e.g., if $x^$⁠, $y^$⁠, and $z^$ are unit vectors in rectangular coordinates with $z^$ being vertical, then $Π=x^x^T+y^y^T$ would represent the preference to align to a horizontal surface), and $0≤c≤1$ models the strength of the preference. Writing $z$ using Eq. (18), along with Eq. (15) then represents the reduced equations for this problem.

Another extension to the Kuramoto model that is often studied is the Kuramoto–Sakaguchi model.²⁴ In this model, the $sin⁡(θj−θi)$ coupling term of the Kuramoto model is replaced with $sin⁡(θj−θi+δ)$⁠. A possible generalization of this to higher dimensions is represented by defining $ρ$ as $ρ=KRz$⁠, where $R$ is a given rotation matrix (for $D=2$⁠, $R$ is the two-dimensional rotation matrix that rotates vectors by an angle of $δ$⁠).

Another $D$-dimensional generalization whose analysis can be facilitated by the ansatz Eq. (14) is the consideration of time delay, $ρ(t)=Kz(t−τ)$⁠, as studied for $D=2$ in Ref. 43.

Also, we note that interactions between multiple communities of Kuramoto-like agents has received attention due to a variety of applications (e.g., Refs. 22,44, and 45), as well as the presence of interesting dynamics, such as chimera states.²² For example, for the case of homogenous natural rotations of $Wξ$ within each community $ξ$⁠,

where the subscript $ξ$ denotes quantities applying to community $ξ$⁠. For a case of generalizing the Kuramoto model, we define the order parameter $zξ$ for community $ξ$ as the average orientation of that community, and take $ρξ$ to be

with $Kξ,ξ′$ representing the coupling between community $ξ$ and $ξ′$⁠. The order parameters $zξ$ can be written in terms of $αξ$ using Eq. (18) by writing the distribution of rotations for the community $ξ$ as $δ(W−Wξ)$⁠.

The Kuramoto model with the order parameter defined as Eq. (16) is the globally coupled Kuramoto model, wherein each agent is coupled to every other agent. In two dimensions, network-based interaction of agents in Kuramoto-like models have been solved for by an application of the ansatz Eq. (5) for a wide range of network topologies, via a mean-field approach.^16,21 An analogous analysis will apply for our generalized ansatz, Eq. (14), for network-based interactions of $D$-dimensional Kuramoto-like units.

There are some strong differences between the case $D=2$ and the case of $D>2$ that must be considered in general. In the case of $D=2$⁠, making the additional assumption that $g(η)$ is a suitable analytic distribution of the scalar parameter $η$ (e.g., a Lorentzian distribution is often employed), allows the integral in Eq. (18) to be performed via a contour integral, and hence requiring the dynamics of $α(η)$ according to Eq. (15) to be calculated for only one or a few particular complex values of $η$⁠.²⁷ In $D=2$⁠, this implies that many problems of the form Eq. (1) with heterogeneous $ηi$ reduce to a system of a small number of ordinary differential equations in the $N→∞$ limit. For our generalization to higher dimension (where $η$ is now a vector parameter with at least two components), we are unable to straightforwardly employ contour integration. Thus, while Eq. (14) represents a strong reduction in the dimensionality of the dynamics as compared to the full system in the $N→∞$ limit, i.e., Eq. (3), it is still not a “low-dimensional system” in the sense of Ref. 27, since we must still calculate the dynamics of $α(η,t)$ as a function of the vector parameter $η$ [as opposed to integrating $η$ away via, e.g., a Lorentzian assumption for $g(η)$].

For the case of homogenous systems, i.e., where $g(η)$ is the Dirac-delta function, the dynamics of the full system Eq. (2) reduces to the single $D$ dimensional differential equation Eq. (15). For the particular case of the Kuramoto model generalized to higher dimensions in the manner given in Sec. IV A, this exactly reproduces recent results by Lohe³⁴ for the case of a finite number of homogenous agents derived in the context of a generalization of the Watanabe–Strogatz (WS) transform.^46,47 In 2014, Tanaka³⁵ considered a generalization of the Kuramoto model with higher dimensional “complex” vectors. In this setup (which is different from ours), he derived an extension of the Ott–Antonsen method in the context of a generalized WS transform. In 2008, Pikovsky and Rosenblum⁴⁸ demonstrated that there is a relationship between the WS transform and the Ott–Antonsen ansatz in the case of the original (⁠$D=2$ in our notation) Kuramoto model. It is possible that similar relationships may exist between our generalization of the Ott–Antonsen ansatz Eqs. (14) and (15) and the generalization of the WS transform described in Refs. 34 and 35—we leave the study of this relationship to future work.

In conclusion, we have developed a technique to tackle the generalization of several Kuramoto-like systems into higher dimensions. While our analysis has only demonstrated the existence of an invariant manifold to the dynamics of Eq. (3), from numerical experiments, we observe for all examined examples of systems given Eq. (2) with a continuous distribution $g(η)$ that this manifold is attracting. That is, initial conditions set up not satisfying Eq. (14) appear to be rapidly attracted toward this invariant manifold. While, in the case of $D=2$⁠, it has been shown analytically that for a broad class of models of the form given by Eq. (1), this manifold is a global attractor of the dynamics,²⁵ proof of attraction for $D>2$ remains an open problem. Given the wide applicability of Eq. (1) and its rich variety of dynamical phenomena, we expect that the generalization to higher dimensions, Eq. (2), may be a useful model system, applicable to diverse situations of interest, while remaining amenable to analysis via the methods developed in this paper.

This work was supported by ONR Grant No. N000141512134 and by AFOSR Grant No. FA9550-15-1-0171.

Inserting the form of $f(σ,η,t)$ from Eq. (14) into Eq. (9), we obtain

Remarkably, the explicit $D$ dependence of the differential equation cancels out, and a differential equation involving only terms that are linear and constant in $σ$ remains. For this equation to be identically zero for each direction $σ$⁠, the linear and constant terms must independently be zero. From the constant term, we obtain

which simplifies to

or alternately

From the $σ$ dependent portion, we get

Since $σ$ is allowed to be in any direction, we can cancel out the $σ$ and obtain a vector equation that must be satisfied. To further simplify this vector expression, we write $∂tα=∂t(|α|α^)=|α|∂tα^+α^∂t|α|$⁠, where $α^$ is a unit vector in the direction of $α$⁠. We can then use Eqs. (2) and (3) to simplify the expression to obtain

Equations (3) and (4) can then be combined to obtain Eq. (15).