Previous results have shown that a large class of complex systems consisting of many interacting heterogeneous phase oscillators exhibit an attracting invariant manifold. This result has enabled reduced analytic system descriptions from which all the long term dynamics of these systems can be calculated. Although very useful, these previous results are limited by the restriction that the individual interacting system components have one-dimensional dynamics, with states described by a single, scalar, angle-like variable (e.g., the Kuramoto model). In this paper, we consider a generalization to an appropriate class of coupled agents with higher-dimensional dynamics. For this generalized class of model systems, we demonstrate that the dynamics again contain an invariant manifold, hence enabling previously inaccessible analysis and improved numerical study, allowing a similar simplified description of these systems. We also discuss examples illustrating the potential utility of our results for a wide range of interesting situations.
The dynamics of systems with many coupled dynamical agents is a subject of increasing importance with a very broad range of applications. Much of the progress in this field has flowed from the discovery of solvable paradigmatic “toy” systems. In many such systems, the agents are assumed to have one-dimensional dynamics, and a certain class of such systems has been shown to be in some sense “solvable” via a novel analytic technique. In our work, we consider a more general class of dynamical systems of coupled agents that may have arbitrary dimension. This is a significantly broader class of systems that contains, but is not limited to the previously described class of systems with one-dimensional agents. We then demonstrate that in this broader class of dynamical systems, we can construct an analogous technique that can be used to solve such systems. Our method provides analytic techniques to allow previously inaccessible mathematical analyses of these systems. We give significant examples applying our method to large systems of interacting higher-dimensional agents, with a particular focus on the Kuramoto model generalized to higher dimensions.
Models of systems of many coupled dynamical agents are useful tools for studying a very wide variety of phenomena.1 Examples include flashing fireflies,2,3 circadian rhythms of mammals,4,5 oscillating neutrinos,6 arrays of Josephson junctions,7 oscillation of footbridges,8 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,1 neurons,13–16 or robots moving on a two-dimensional plane,18 among others. Many such models, possibly involving network-based interactions,16,21,22 such as the Kuramoto model,23 the Kuramoto–Sakaguchi model,22,24 and models of theta neurons,13 among others,25 reduce to the form
where represents the state of the agent, is a (possibly vector) constant parameter that is associated with the agent, is its “natural frequency,” is the total number of agents, and is a common field that acts on each agent, dependent on the agent’s parameter , and indicates a dependence on the set of states in the form of an average over of a function of the angle .
For example, the well-studied Kuramoto model23,26 can be expressed in the form of Eq. (1) by choosing , independent of , and choosing independent of and , allowing to be replaced by . 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 (). 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 in the -dimensional space. Alternately, we may think of as specifying a point on the unit sphere in -dimensional space. Reference 30 notes that the vector 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 is of most interest.) For , the unit vector is determined by its scalar orientation angle 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.
II. GENERALIZING KURAMOTO-LIKE AGENTS TO HIGHER DIMENSIONS
In two recent papers,40,41 we constructed a generalization of the Kuramoto model to dimensions. Here, we consider an even more general setup, where we consider a generalization to Eq. (1) to a system in dimensions,
where for each , is a real -dimensional unit vector, , is an arbitrary real -dimensional vector, which can be thought of as a common field that affects each agent in an dependent fashion, is a real antisymmetric matrix, is a (possibly vector) constant parameter associated with each agent, and, as earlier, indicates a dependence on the set of all states in the form of the average over of a function of the unit vectors (we further quantify this dependence on later). For example, in the context of flocking agents in dimensions, represents the orientation of the agent, represents a “goal” orientation to which the agent attempts to align itself, and 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 .40 Note from the form of Eq. (2) that the dot product of the right-hand side of Eq. (2) with is identically zero so that , as required by our identification of as a unit vector. Thus, the dynamics of each is restricted to the -dimensional surface, , of the unit sphere, . For , choosing , , and
We now consider the limit of a large number of agents and denote by the distribution of agents on , such that is the fraction of agents that lie in the -dimensional differential element on the surface centered at at time , and have an associated parameter within the differential element centered at . Since the associated parameter for each agent is time independent, we define
Noting that Eq. (2) specifies the vector field of the flow controlling the dynamics of the distribution , we write a continuity equation for ,
where the velocity field is given by , and represents the divergence of a vector field , along the surface . This can be done if the dependence of and on can be specified as a functional of that is not explicitly dependent on . [A simple example of such a dependence on would be the average value of for some given function , which can be written as .] Following Appendix B of Ref. 40, Eq. (3) can be rewritten as
where is the gradient of a scalar field projected on the surface .
III. ANALYTIC SOLUTION IN THE LIMIT OF LARGE SYSTEMS
where is a complex scalar function of and , , reduces Eq. (1) to the following -independent form:
The form Eq. (5) represents an invariant manifold in the space of possible distributions , that satisfy the continuity equation Eq. (3) for . Furthermore, previous work25,42 has shown that initial conditions for 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 . Noting that 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 for arbitrary dimension :
where is a real -dimensional vector such that , is a yet-to-be-determined constant, and is a scalar normalization chosen to ensure that
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., , since, in general, and will not be zero for all , we require that
With in Eq. (7) determined, we now obtain the normalization constant . To perform the integral in Eq. (8), without loss of generality we take the vector to be along the axis. For an arbitrary point on , we denote the angle between and by . In particular, we note that the distance of the point from the axis is . For a coordinate system on the surface , we use as one of the coordinates, denoting position with respect to on the sphere. From the symmetry of in Eq. (7) about the direction , we see that the integrals over these remaining coordinates give the surface area of the dimensional surface of a sphere with radius embedded in dimensions, where is the area of the sphere of unit radius in dimensional space. Thus, Eq. (8) becomes
which can be evaluated to give
where is a constant dependent only on . This results in
giving the form of the ansatz for arbitrary dimensions as
which, for , 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 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 ). Thus, analogously to Eq. (6), we have a -independent description of the dynamics of . This is our main result.
IV. EXAMPLE SYSTEMS
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 dimensions40 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. The Kuramoto model generalized to higher dimensions
A generalization of the Kuramoto model with homogenous oscillators to arbitrary dimension was introduced by Olfati-Saber in 200630 in the context of flocking dynamics, consensus protocols, and opinion dynamics. This was later generalized to heterogeneous systems by Chandra et al.40 For generalization to dimensions, a system order parameter, , can be defined as
The magnitude of is a measure of the coherence of the set of agents . The common field is then defined as the -independent function,
where is a coupling constant. By interpreting the vector parameters in as the independent elements of a -dimensional antisymmetric matrix , we can replace in integrals by , where is a distribution of antisymmetric matrices. In cases such as these where is independent of and , we interpret as the “natural rotation” of .
In the limit of infinite system size, with a distribution of agents given according to Eq. (14),
This now allows us to use Eq. (15) with some given 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, , 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, . This makes the -integral in Eq. (18) trivial, allowing a direct representation of in terms of . Further, is only dependent on time (rather than and ). This represents a very large simplification in the complexity of the dynamics of the system of agents, since Eq. (15) is now a single -dimensional ordinary differential equation which represents the collective dynamics of the , -dimensional system of coupled differential equations in Eq. (2). The utility of this result is demonstrated for in Fig. 1(a), where we show (plotted in black) the time-series for as generated from a system of agents (approximating the 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 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- dynamics of the full system of interacting agents.
For the case of heterogeneous agents, in Eq. (18) depends on , and we perform the integral in a Monte-Carlo fashion. We randomly choose values of from the given distribution and simulate the dynamics of the corresponding s. These randomly chosen s are then used as the Monte-Carlo samples to evaluate according to Eq. (18), simulating the dynamics of the system in the limit by only simulating the dynamics of variables. Here, we choose an isotropic distribution 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 antisymmetric. Results are shown in Fig. 1(b) for , where Monte-Carlo samples were chosen to evaluate the 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 limit, approximated by a simulation of agents. Note how simulating the dynamics of Monte-Carlo samples yields a smooth curve approximating the noisy curve generated by simulating the individual dynamics of agents. Initial conditions for the full system were chosen as a bimodal distribution of s, independent of the corresponding , 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 , corresponding to an approximately uniform distribution of 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 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 .25 Full system simulations with initial conditions described by a uniform distribution in [and hence lying on the invariant manifold Eq. (14) for ] yielded a curve that is not discernibly different from the curve presented in Fig. 1(b).
B. Applications of Eq. (15) to previous results on the generalized Kuramoto model
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 -agent dynamics and the simulation of the reduced equation Eq. (15) holds at all values of the coupling constant . 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 exhibits the unusual behavior of Instability-Mediated Resetting: If the coupling strength 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 in the full system simulation were chosen independently to be uniformly random over the sphere . 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 , corresponding to an approximately uniform distribution of 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 is to be distributed independently and uniformly over the sphere . In the limit, this is the distribution corresponding to setting in Eq. (14). [In the context of the stability analysis presented in Sec. III B of Ref. 40, we are considering the case of .] Thus, we are interested in the stability analysis of Eq. (15) about . Performing a first order expansion of Eq. (15) for ,
Assuming , we obtain
Note that in the limit of , Eq. (19) can be written as
Multiplying Eq. (21) by and integrating, we obtain
Note that in three dimensions the linear transformation can be represented as the cross product . Without loss of generality we may choose a basis that block-diagonalizes , corresponding to the choice of and
This can now be inserted into Eq. (23) and written in a basis independent format as
We choose the distribution to be an “isotropic” distribution (i.e., a distribution that is invariant to orthogonal transformations) which can hence be written as , where represents the isotropic distribution of rotation directions, and 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
C. Other examples
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 , , and are unit vectors in rectangular coordinates with being vertical, then would represent the preference to align to a horizontal surface), and models the strength of the preference. Writing 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.24 In this model, the coupling term of the Kuramoto model is replaced with . A possible generalization of this to higher dimensions is represented by defining as , where is a given rotation matrix (for , is the two-dimensional rotation matrix that rotates vectors by an angle of ).
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.22 For example, for the case of homogenous natural rotations of within each community ,
where the subscript denotes quantities applying to community . For a case of generalizing the Kuramoto model, we define the order parameter for community as the average orientation of that community, and take to be
with representing the coupling between community and . The order parameters can be written in terms of using Eq. (18) by writing the distribution of rotations for the community as .
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 -dimensional Kuramoto-like units.
There are some strong differences between the case and the case of that must be considered in general. In the case of , making the additional assumption that 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 .27 In , this implies that many problems of the form Eq. (1) with heterogeneous reduce to a system of a small number of ordinary differential equations in the 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 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 as a function of the vector parameter [as opposed to integrating away via, e.g., a Lorentzian assumption for ].
For the case of homogenous systems, i.e., where is the Dirac-delta function, the dynamics of the full system Eq. (2) reduces to the single 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 Lohe34 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, Tanaka35 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 Rosenblum48 demonstrated that there is a relationship between the WS transform and the Ott–Antonsen ansatz in the case of the original ( 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 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 , 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,25 proof of attraction for 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.
APPENDIX: PROOF OF EQ. (15)
Remarkably, the explicit 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
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 , where is a unit vector in the direction of . We can then use Eqs. (2) and (3) to simplify the expression to obtain