Self-organization on Riemannian manifolds

We consider an aggregation model that consists of an active transport equation for the macroscopic population density, where the velocity has a nonlocal functional dependence on the density, modelled via an interaction potential. We set up the model on general Riemannian manifolds and provide a framework for constructing interaction potentials which lead to equilibria that are constant on their supports. We consider such potentials for two specific cases (the two-dimensional sphere and the two-dimensional hyperbolic space) and investigate analytically and numerically the long-time behaviour and equilibrium solutions of the aggregation model on these manifolds. Equilibria obtained numerically with other interaction potentials are also presented.

1. Introduction. The literature on self-collective behaviour of autonomous agents (e.g., biological organisms, robots, nanoparticles, etc) has been growing very fast recently. One of the main interests of such research is to understand how swarming and flocking behaviours emerge in groups with no leader or external coordination. Such behaviours occur for instance in natural swarms, e.g., flocks of birds or schools of fish [14,25]. Also, swarming and flocking of artificial mobile agents (e.g., robots) in the absence of a centralized coordination mechanism is of major interest in engineering [37,38]. Consequently, there exists a variety of models for swarming or flocking, ranging from difference equations (discrete in both time and space) [27,25,26] to ordinary/partial differential equations (continuous in time and discrete/continuous in space) [23,45,28].
In this paper we consider an aggregation model that consists in an integrodifferential equation for the evolution of a population density ρ(x, t) on a Riemannian manifold M : into the model, which in this class of models amounts to incorporating it into the interaction potential. Similar to the setup in Euclidean spaces, we assume that there are no limitations on the mutual sensing of individuals, that is, all individuals sense each other. For the example above (coordination of mobile agents), the assumption is valid provided the agents are set to communicate globally with each other via a central unit. For coordination of biological agents, one can assume that individuals possess a sensing mechanism (such as smell for instance) which enables them to communicate with each other regardless of the geometry of the space they live in. We note here that from a mathematical point of view, including local sensing and/or asymmetry/anisotropy in the model is expected to bring up major challenges (see for instance [29] for a study of anisotropic interactions in model (1) posed on R 2 ).
Model (1) is in the form of a continuity equation for the density ρ. Note that this is an active transport equation, as the velocity field defined in (1b) depends on ρ. The interpretation of (1) as an aggregation model is in fact encoded in (1b): by interacting with a point mass at location y, the point mass at x moves either towards or away from y. The velocity v at x computed by (1b) takes into account all contributions from interactions with point masses y ∈ M . Also, in geometric terms, the continuity equation (1a) represents the transport of the volume form ρµ along the flow on M generated by the tangent vector field v [1, Chapter 8]. Equivalently, the mass in each subregion of the manifold remains constant through time, as the subregion evolves by the flow. In particular, the total mass m = M ρ dµ is conserved. This geometric interpretation of the aggregation equation will play a major role in the paper.
In the present research we investigate solutions for model (1) posed on two simple manifolds: the two-dimensional sphere and the two-dimensional hyperbolic space. We design interaction potentials for which the equilibrium densities have a simple structure that is amenable to analytic investigations (specifically, equilibria are constant on their supports). The strategy for designing such potentials is inspired by [32,31], where similar goals were pursued for the aggregation model in Euclidean spaces. Our work is the first to demonstrate that model (1) set up on manifolds leads to swarming behaviour that can be studied mathematically. We also perform numerical investigations for other interaction potentials and find a diverse set of equilibrium solutions. Given the outstanding interest shown recently in the aggregation model on R n , we hope to have set up the stage for further studies on swarming and self-collective behaviour on general manifolds, opening new perspectives and motivating applied mathematicians to expand the research on this class of models to novel applications.
We finally note that the continuum model (1) has an immediate discrete/ODE analogue, in which one considers the time evolution of a fixed number of individuals/particles on a manifold. While the discrete model is used in this paper only for numerical purposes (Sections 3.3 and 4.3), it has an interest in its own (e.g., for applications in biology or robotics). Specifically, equilibria of constant densities, as achieved in this paper, translate in the discrete setup to agents covering uniformly (with respect to the metric) a certain space modelled as a Riemannian manifold. This relates to the coverage problem in robotics, where the goal is to have a group of robots well-distributed over a region/area, so that it achieves an optimal coverage needed for surveillance/tracking. Such configurations are referred to as balanced or anti-consensus in the control literature [49]; see also [24] for a recent survey on coordinated control of robot systems. The ODE model has been extensively employed for the aggregation model in Euclidean space, either as a numerical tool [41,53,5,32], or in rigorous analysis studies that establish the passage from discrete to continuum by mean-field limits [21]. Consequently, we also expect a rich use and interest for the discrete aggregation model on Riemannian manifolds, as set up in the present paper.
The summary of the paper is as follows. In Section 2 we present some preliminaries and provide motivation to the present research. In Section 3 we set up the aggregation model on the two-dimensional sphere, with a certain choice of the interaction potential, and investigate the long-time behaviour of solutions. A similar study is done in Section 4 for the model on the two-dimensional hyperbolic space. In Section 5 we showcase some numerical simulations with other interaction potentials and discuss a new application, namely the aggregation model on the rotation group SO(3).

Preliminaries and motivation.
In this section we present briefly some standard concepts from differential geometry that will be used in the sequel, as well as bring motivation to the present work.
Local coordinates. Consider a generic n-dimensional Riemannian manifold (M, g) with local coordinates (x 1 , x 2 , . . . , x n ) and local metric coefficients g ij .
Expressed in the local basis ∂ ∂x 1 , ∂ ∂x 2 , . . . , ∂ ∂x n of the tangent space, the gradient (with respect to the metric g) of a scalar function f on M is the tangent vector ∇ M f given by where g ij are the entries of the metric's inverse, and we used the Einstein convention on index summation [42]. Given a tangent vector field F = F i ∂ ∂x i on M , its divergence ∇ M · F is given in local coordinates by where |g| denotes the determinant of the metric. By combining (4) and (5) one then finds the Laplace-Beltrami operator ∆ M in local coordinates: Aside from the operators above, we will also use the representation in local coordinates of the canonical volume form µ on (M, g) [1]: Continuity equation on Riemannian manifolds. As discussed in [1,Chapter 8.2], the continuity equation (1) posed on Riemannian manifolds has the following interpretation.
Consider the flow Φ t : M → M on M generated by the vector field v, that is: For a fixed time t, denote by ρ t (x) = ρ(x, t), and also recall that µ represents the canonical volume form on M . Then, the continuity equation (1) is equivalent to the transport of the volume form ρ t µ along the flow Φ t , that is: Here, Φ * t denotes the pull-back by Φ t . Moreover, this is equivalent to: where J(Φ t ) denotes the Jacobian of Φ t (with respect to the canonical volume form µ).
We also note that for a smooth, invertible flow map Φ t , the Jacobian J(Φ t ) satisfies (see [1,Proposition 7.1.10]): Other properties of the model. Energy. The aggregation model (1) in Euclidean spaces can be formulated as a gradient flow on the space of probability measures with finite second moments, equipped with the 2-Wasserstein metric [4]. Such interpretation exists as well for the model on Riemannian manifolds with Euclidean pairwise interactions from [54]. The general model (1) also has an energy associated to it, which decays with time. Indeed, define: By formally computing the evolution of the energy in time, we find: For the first equal sign in the derivation above we used the symmetry of the potential and equations (1) and (2). For the second equal sign we used the formula [1, Proposition 6.5.17]: where ∇ ρv (K * ρ) denotes the covariant derivative of the scalar function K * ρ along the vector field ρv. Integrating the equation above over M , by divergence theorem (either assuming that M has no boundary or that M is non-compact, but density vanishes at infinity) the first term in the right-hand-side yields zero. Then, (12) follows from ∇ ρv (K * ρ) = ρ∇ v (K * ρ) (linearity of the covariant derivative) and the definition of gradient on M by which: Though we do not make any further energy considerations in this work, we believe that there is a rich potential for applications of the theory on gradient flows as developed in [4,48] to the model investigated here. In addition, the study of equilibria of model (1) as minimizers of energy E seems a very interesting direction to pursue as well. This approach has been proven very successful for the model in R n [22,50,15,5].

RAZVAN C. FETECAU AND BERIL ZHANG
Centre of mass. An important property of model (1) in Euclidean space R n is the conservation of the centre of mass. This can be derived easily as follows: multiply equation (1) by x and integrate over R n , then use integration by parts and the symmetry of K to conclude that R n ρ(x)dx remains constant in time. In the context of Riemannian manifolds, as a generalization of the usual centre of mass in R n , we consider the L 2 Riemannian centre of mass (also known as the Karcher mean) [39,2]. The L 2 Riemannian centre of mass (simply referred throughout as the Riemannian centre of mass) of a subset A ⊂ M is a minimizer in M of the function One can check indeed that for M = R n , f has a unique minimizer which coincides with the usual centre of mass of set A. For general manifolds, existence and uniqueness of the Riemannian centre of mass, along with numerical methods for finding it, are delicate issues [2,3]. We make the important observation here that model (1) does not necessarily conserve the Riemannian centre of mass. Explaining this fact in detail would be an unnecessary detour for the purpose of this paper. To give some intuition on why such result is not expected to hold in general however, we point our that in R n , by symmetry of K (see (3), which in R n it amounts to K(x, y) = K(|x − y|)), ∇K is antisymmetric in x, y, i.e., ∇ x K(x, y) = −∇ y K(x, y). This property is used in an essential way to show conservation of centre of mass in R n . On the other hand, such a property would not even make immediate sense on general Riemannian manifolds, as the manifold gradients with respect to x and y lie in tangent spaces at different points (∇ M,x K(x, y) ∈ T x M and ∇ M,y K(x, y) ∈ T y M ). The lack of conservation of centre of mass brings an additional challenge for analytical investigations of model (1). We return to this point at the end of Section 4.
The aggregation model in Euclidean space. To motivate and put the present work in context, we present first some key ideas and calculations for the aggregation model in Euclidean space, that is, model (1) set up on M = R n endowed with the standard Euclidean metric. We refer here to the research in [32,31], where one of the main goals was to design interaction potentials for the aggregation model in R n that yield equilibrium states which are biologically relevant and at the same time, are simple enough to be investigated analytically.
Of particular importance in [32] is an attractive-repulsive interaction potential that yields equilibria of constant densities and compact support. Specifically, this potential consists of Newtonian repulsion and quadratic attraction, which in two dimensions (n = 2) it amounts to: where G(x; x ) denotes the Green's function for the negative Laplacian in R 2 : Note that and hence, from (1b), where m is the (constant) total mass. The key conclusion from (16) is that with this choice of interaction potential, ∇ · v is a local quantity, despite the fact that v itself is nonlocal, as given by (1b) through a convolution.
and write the continuity equation (1a) as The left-hand-side of (18) is the material derivative of the density ρ along the flow generated by v. Hence, by (18) and (16), along the particle path Φ t (α) that originates from location α (see equation (8) Note that the right-hand-side of (19) depends only on ρ evaluated along the carrying characteristic Φ t (α), and not on values of ρ along any other other characteristic. Therefore, this ODE for ρ(Φ t (α), t) can be investigated individually.
From (19) one infers immediately that the solution ρ(Φ t (α), t) approaches the value 2m as t → ∞. Consequently, equilibria of model (1) in R 2 , with potential (13), have constant densities on their supports. In fact, it was shown in [32,10] that the constant density supported on a disk, i.e., is a global attractor for solutions of model (1) with potential (13) in R 2 ; see Figure  1 for a numerical simulation using a particle method. The considerations above extend to the aggregation model set in Euclidean R n of arbitrary dimension. Indeed, one can construct an interaction potential that leads to similar long-time behaviour of solutions (for G(x, x ) one would have to use the fundamental solution of the negative Laplacian in R n instead).
Motivation for the present work. The ideas above can be extended to model (1) posed on arbitrary Riemannian manifolds (M, g). Indeed, similar to (15), consider an interaction potential that satisfies where C > 0 a constant. Then, as in (16), we can calculate from (1b), (21), and the conservation of mass: On a Riemannian manifold (M, g), analogous to (17), we have where ∇ v ρ denotes the covariant derivative of ρ along the tangent vector field v.
Then, one can proceed as in the Euclidean case (see (18) and (19)), and get from (1a), (22) and (23) an ODE for the evolution of ρ(Φ t (α), t): Here, the left-hand-side denotes the material derivative of the density along the flow Φ t defined in (8): By (24), densities along particle paths approach the constant value Cm. Therefore, as in the Euclidean case, equilibrium states have constant densities on their supports. The primary purpose of the current paper is to apply the considerations above for two specific geometries: the two-dimensional sphere in R 3 and the two-dimensional hyperbolic plane. In particular, we find interaction potentials that satisfy (21) and then investigate the dynamics and equilibria of the aggregation model (1) in these setups. 3. Aggregation model on the sphere. In this section we set up the aggregation model (1) on the 2-dimensional sphere in R 3 . We then construct a certain interaction potential and investigate analytically and numerically the long time behaviour of the solutions.
3.1. Model setup. Let x = xe 1 + ye 2 + ze 3 denote the position of a particle in R 3 , where {e 1 , e 2 , e 3 } is the standard orthonormal Cartesian basis.
The tangent vectors at (θ, φ) on S are x θ = cos θ cos φ e 1 + cos θ sin φ e 2 − sin θ e 3 , x φ = − sin θ sin φ e 1 + sin θ cos φ e 2 , and the first fundamental form (the metric) is given by: The metric matrix has determinant |g| = sin 2 θ and its inverse has entries: Gradient and Laplace-Beltrami operator on sphere. For a scalar function f on the sphere, its surface gradient is given by (see (4)): Also, by (6), we have the expression of the Laplace-Beltrami operator in spherical coordinates: Choice of interaction potential. We look for an interaction potential on the sphere that satisfies (see (21)): where ∆ S is the Laplace-Beltrami operator given by (29), and C > 0 a constant. We remark from the start that since the sphere is a closed manifold, the constant C cannot be arbitrary; by the solvability condition for (30), C = 1 4π . Based on the observation above, we choose K to be the Green's function in the generalized sense of ∆ S [40], i.e., we set where Here, d(x, x ) denotes the spherical distance between the points x and x . It is a simple exercise to check indeed that G S satisfies: In local coordinates (θ, φ), (θ , φ ) (corresponding to x and x , respectively), the spherical distance is given by the law of cosines on sphere: Also, using local coordinates in (2) we have: where by an abuse of notation, we used K(θ, φ, θ , φ ) for K(x, x ). The interaction potential K is purely repulsive. Indeed, since K(x, x ) only depends on the geodesic distance between the points, it is enough to take x as the North pole of the sphere ((0, 0) in local coordinates). Then, d(x, x ) = θ and by (1b), (28) and the choice of K in (31)- (32), upon interacting with the North pole, the point x moves in the direction Hence, all points on S sense a repelling force (except the South pole for which the expression above vanishes) -see Figure 2(a) for an illustration.

3.2.
Asymptotic convergence to equilibrium. For simplicity, set the total mass m = 1. From (1b), (31), (33) and the conservation of mass, we find (see also (22)): Hence, along the flow Φ t generated by the vector field v on S, the density ρ(Φ t (α), t) satisfies (see (24)): As noted above, along particle paths, densities ρ(Φ t (α), t) approach a constant value (here, the constant density at equilibrium is 1 4π ). The ODE (38) can be solved exactly in fact, with solution: where ρ 0 is the initial density.
Also, from (9) and (39), one can find an exact explicit expression for the Jacobian of the flow map: where for notational convenience we have denoted Note that J(α, t) > 0 for all t, guaranteeing that the particle map is invertible for as long as it exists.
Theorem 3.1 (Global attractor for sphere). Consider model (1) set up on the sphere S, with the interaction potential K given by (31)- (32). Assume that the model has a global in time C 1 solution ρ(x, t) and that the flow map Φ t : S → S is C 1 and invertible for all t > 0. Then, solutions ρ(x, t) approach asymptotically, as t → ∞, a constant equilibrium density supported over the entire sphere.
Proof. The proof is essentially provided by the considerations above. From (38) and (39), along particle paths that originate from the support of ρ 0 (i.e., ρ 0 (α) = 0), one has Consequently, equilibria must have constant density 1 4π on their support. And since the total mass m = 1 is conserved through the time evolution, then necessarily the support of the equilibrium state must be the entire sphere (up to a zero measure set). This shows thatρ is a global attractor.
Symmetric initial density. When the initial density is symmetric with respect to rotations about a North-South axis, we can find explicit expressions for the evolution of the particle paths, and hence, show directly the result in Theorem 3.1. While the calculation is interesting in itself, it is key to investigating attractors on the hyperbolic plane (Section 4.2).
For simplicity, assume the initial density is symmetric with respect to rotations about the z-axis. In spherical coordinates, this amounts to considering an initial density ρ 0 that depends only on θ, but not on φ; by an abuse of notation we will denote ρ 0 (α) by ρ 0 (θ). Note that a symmetric initial density results in a symmetric solution for all times. Indeed, for the initial density ρ 0 (θ), evaluate the velocity at a generic point (θ, φ) using (1b) and (35). One finds that the velocity only depends on θ and also, that it points in the direction x θ , i.e., the point moves along the meridian φ = const. The same argument applies to a symmetric density at an arbitrary time instance, hence the symmetry is preserved through time evolution.

RAZVAN C. FETECAU AND BERIL ZHANG
Then, by (10) and (43) (see also notation (41)), we find in local coordinates that: Note that, as expected, the Jacobian depends only on θ and not on φ, as the latter coordinate remains constant along the flow. Now, we get from (40) and (44) the following differential equation for λ(θ, t): sin θ.
To find the asymptotic behaviour of the particle trajectories, one can send t → ∞ in equation (45). We find: where Note that the (conserved) unit mass can be written in spherical coordinates as: Now consider a symmetric domain θ 0 ≤ θ ≤ θ 1 that contains the support of ρ 0 . Then, for any particle trajectory that originates from θ < θ 0 (outside the support of ρ 0 ), we have by (46) that Λ θ = 0. In other words, all trajectories starting from θ < θ 0 approach the North pole as t → ∞. Similarly, for any particle trajectory that originates from θ > θ 1 (also outside the initial support), one finds by (46) and (47) that Λ θ = π, so the trajectory approaches the South pole as t → ∞. On the other hand, by (46), Λ θ is monotonic and continuous in θ, so trajectories starting from inside the initial support will spread over the entire sphere, as expected by the result in Theorem 3.1 -see Figure 3 for a numerical illustration using a particle method. Remark 1. The considerations above hold for initial densities that are symmetric with respect to rotations about any North-South axis of the sphere. The same exact expressions of the particle trajectories hold, upon a rotation of the coordinate axes.
3.3. Numerical results. We will use the discrete particle system associated to the macroscopic model (1). Set the total mass m = 1 and consider N particles of equal mass (so each particle has mass 1 N ). Let x i (t) represent the location in R 3 of the i-th particle. Equation (1b) can be written in discrete form to express the velocity of particle x i in terms of the locations x j (j = i) of the other particles. Hence, one arrives at the discrete particle system: The configuration remains symmetric for all times and evolves into a uniform particle distribution supported over the entire sphere -see (42). The solid lines represent the trajectories of the particles indicated by stars in figure  (a).
In Euclidean settings, the rigorous mean-field limit of the particle system (48) was established in [21]. Specifically, it was shown that the empirical distribution associated to (48), i.e., µ N (t) = 1 N N i=1 δ xi(t) , converges weak- * as measures to a weak solution ρ(t) of the macroscopic model (1). The result holds for a general class of potentials, including repulsive-attractive potentials that have a (strictly better than Newtonian) singularity at origin.
We write both sides of (48) in the local spherical basis at x i . By chain rule, the left-hand-side can be expanded as: while for the right-hand-side we use (28) to get: By matching the coefficients on each side, we obtain the following ODE system for the spherical coordinates: The results presented below are obtained by solving numerically the particle system (49) with the classical (4th order) Runge-Kutta method. We note that the particle method is very suitable here, as it complements the Lagrangian approach from Section 3.2. Simulating the discrete system is in fact the main numerical tool used for the model in R n .
Numerical simulations on sphere. We use the particle method to confirm the theoretical findings from Section 3.2. First, we consider an initial configuration which is symmetric about the z-axis (Figure 3(a)). To generate this initial density, we first placed particles randomly along a given meridian (φ = const., π 8 < θ < 3π 8 ), and then rotated the configuration about the z-axis. We find indeed that particles evolve along meridians (see the trajectories of the particles indicated by stars) into a uniform particle distribution supported over the entire sphere.
The equilibrium (42) is a global attractor (Theorem 3.1). Figure 4 shows a numerical validation of this result: a random initial particle distribution (Figure 4(a)) evolves into the expected equilibrium state (Figure 4(b)). The θ and φ coordinates of the initial configuration were drawn randomly from the intervals π 8 , 3π 8 and 0, π 2 , respectively. The solid lines in Figure 4(b) represent several individual trajectories, corresponding to the particles indicated by stars.  4. The aggregation model on the hyperbolic space. In this section we set up the aggregation model on the 2-dimensional hyperbolic space and investigate the dynamics of its solutions for interaction potentials that lead to equilibria of constant densities.

Model setup.
We use the hyperboloid model of the two dimensional hyperbolic space [16]. Specifically, we consider the upper sheet of the two-sheeted hyperboloid: H = {(x, y, z) ∈ R 3 | x 2 + y 2 − z 2 = −1 and z > 0}, embedded in R 3 endowed with the Minkowski inner product x, x = xx + yy − zz .
The tangent vectors at (θ, φ) on H are x θ = cosh θ cos φ e 1 + cosh θ sin φ e 2 + sinh θ e 3 , x φ = − sinh θ sin φ e 1 + sinh θ cos φ e 2 , and the metric coefficients are given by: The determinant of the metric is |g| = sinh 2 θ and its inverse given by: The hyperbolic distance d(x, x ) between two points x and x on H, of local coordinates (θ, φ) and (θ , φ ), respectively, can be found from the hyperbolic law of cosines: Gradient and Laplace-Beltrami operators on H. From (4) and (6) we find: and (55) Choice of interaction potential. The Green's function of the negative Laplacian on H is given by (see [40]): Indeed, using (55) and (53), one can check that −∆ H G(x; x ) = δ x holds in the sense of distributions. Motivated by the considerations made in Section 2, we look now for a function with positive constant Laplacian. By elementary methods (using (55)), one finds that satisfies ∆ H A(θ) = 1. Working in more generality, and using the hyperbolic distance between points, it also holds that: ∆ H A(d(x, x )) = 1. Consequently, we choose the interaction potential as

RAZVAN C. FETECAU AND BERIL ZHANG
which by the considerations above satisfies Remark 2. We note that one can multiply the function A by any positive constant and reach a model with similar properties. Indeed, the interaction potential Consequently, any C > 0 serves the purpose of the present study, which is to have a model that evolves into an equilibrium of constant density (the constant density in this case would be Cm, see (24)).
To determine the nature of the potential (58) we proceed as in the sphere case. With no loss of generality, we take x to be the vertex of the hyperboloid (of local coordinates (0, 0)) and x = (θ, φ) a generic point on H, so that d(x, x ) = θ. From (1b), (54) and (58) (see also (56) and (57)), by interacting with the vertex, the point x moves in the direction Specifically, the first term in the right-hand-side of (36), resulting from the Green's function component of K, is repulsive, while the second term, corresponding to term A(θ) of K, is attractive -see Figure 2(b) for an illustration. Repulsion and attraction balance each other at θ = cosh −1 1 + 1 2π , while the net effect of the two interactions is repulsive or attractive at distances θ smaller or larger than this value, respectively (short range repulsion and long range attraction).
Remark 3. The attractive-repulsive character of the interactive potential (58) mimics very well that of the interaction potential (13) from the Euclidean case. Namely, it contains a repulsive component due to the Green's function of the negative Laplacian (on H), and an attractive component that is smooth and unbounded (growing at infinity). Remarkably, as shown next, this similar structure results in similar long-time behaviours of solutions of the two models, that is, to approach asymptotically equilibria of constant densities supported on geodesic disks.

4.2.
Asymptotic convergence to equilibrium. Set the (conserved) total mass m = 1. Similar to the Euclidean and the sphere cases, from (1b), (58), (59) and the conservation of mass, we have (see (22)): Hence, along the flow Φ t (α), the density ρ(Φ t (α), t) evolves according to (see (24)): We conclude from (62) that densities ρ(Φ t (α), t) approach the constant value 1 along particle trajectories, and consequently, the equilibrium is constant (= 1) on its support. As for the sphere, (62) can be solved exactly, for an initial density ρ 0 : Also, from (9) and (63), we find the Jacobian of the flow map: where we have used the notation (41). Note again that J(α, t) > 0 for all times, so the particle map is invertible for as long as it exists.
Regarding the equilibrium configuration, by conservation of mass, we can infer immediately that the surface area of the equilibrium density is 1. However, (62) does not provide any information about the shape of the equilibrium's support. As opposed to the sphere case, where we showed the asymptotic convergence for arbitrary initial densities (Theorem 3.1), for the hyperboloid we will only prove asymptotic convergence for symmetric initial data. The result is the following theorem.
Theorem 4.1 (Attractor for hyperbolic plane: symmetric initial density). Consider model (1) on H, with the interaction potential (58). Take a compactly supported initial density ρ 0 on H that is symmetric 1 about the vertex of the hyperboloid, and assume that model (1) with this initial data has a global in time C 1 solution ρ(x, t), with a flow map Φ t : H → H that is C 1 and invertible for all t > 0. Then, ρ(x, t) approaches asymptotically, as t → ∞, a constant equilibrium density supported over a geodesic disk centred at the vertex.
Proof. We follow very similar considerations as in Section 3.2. The flow map Φ t : H → H (assumed to be smooth and invertible) has a Jacobian J(Φ t ) that satisfies (10) [1, Proposition 7.1.10]. Here, µ represents the canonical volume form on H which in local coordinates (50) is given by: Using an argument similar to that for the sphere, one infers that an initial density that is symmetric about the vertex results in a symmetric solution for all times (as points move along meridians with velocities which only depend on their θ coordinate). Hence, the flow map in local coordinates takes a point α ∈ H of coordinates (θ, φ) into Φ t (α) ∈ H of coordinates (λ(θ, t), φ), for some function λ. Following exactly the same steps as in Section 3.2, we find from (10) and (65): and further, by using (64), we get a differential equation for λ: Note that in (67) we abused notation and wrote (due to symmetry of initial data) ρ 0 (θ) for ρ 0 (α). The left-hand-side in (67) is a total derivative. Upon integration, we find: where we used the symmetry of the flow, by which the vertex of the hyperboloid (θ = 0) remains fixed (i.e., λ(0, t) = 0 for all t). Equation (68) provides an exact explicit expression for the flow map (θ, φ) → (λ(θ, t), φ).
For the asymptotic behaviour we pass t → ∞ in (68). We find: where Recall that we have set the total mass to be 1, which in local coordinates yields: The initial density was assumed of compact support, so consider θ 0 > 0 such that the support of ρ 0 is contained in 0 ≤ θ ≤ θ 0 . Then, for any particle trajectory that originates from θ > θ 0 , by (69) and (70), we infer that R θ = R, where By elementary geometric considerations [16], R represents the radius of the geodesic disk of unit area centred at the vertex of the hyperboloid. Hence, all trajectories starting from θ > θ 0 approach the geodesic circle of radius R as t → ∞. Finally, by (69), R θ is monotonic and continuous in θ, so trajectories starting from inside θ < θ 0 end up inside the geodesic circle of radius R. Combining these facts with the conservation of mass and the fact that ρ → 1 along all particle paths as t → ∞, we infer that symmetric solutions of model (1) approach asymptotically an equilibrium of constant density (= 1) supported on the geodesic disk of radius R centred at the vertex. That is, initial densities that are symmetric about the vertex, are globally attracted to: where x = (θ, φ) and R is given by (71) -see Figure 5 for a numerical illustration using a particle method.
Remark 4. Theorem 4.1 holds more generally, for initial densities that are symmetric with respect to an arbitrary point on H. Since the interaction potential (and consequently, the interaction forces) only depend on the geodesic distance between points, by symmetry, particles flow along geodesic rays through the symmetry point. The exact expressions for the particle paths, as well as the considerations on the asymptotic behaviour could be then adapted to this more general context.

4.3.
Numerical results. We use the particle method detailed in Section 3.3. Specifically, we consider N particles x i on H and evolve them according to the discrete model

SELF-ORGANIZATION ON RIEMANNIAN MANIFOLDS 415
In local coordinates (see the derivation of system (49) for the sphere), this amounts to solving: (74) Again, the particle method is very appropriate for numerical simulations, given the method of characteristics used in the theoretical results (Section 4.2). The results presented below correspond to numerical solutions of (74) using the classical Runge-Kutta method.
Symmetric initial data. We first validate numerically Theorem 4.1 and consider an initial particle configuration that is symmetric about the vertex (Figure 5(a)). This initial configuration was generated by first placing particles randomly along a given meridian (φ = const., 0.2 < θ < 1.25), and then rotating the configuration about the z-axis. We find that particles evolve along meridians (see the trajectories of the particles indicated by stars) into the symmetric configuration shown in Figure  5(b); see also the zoom-in plot in Figure 5(c). The distance from the vertex of the hyperboloid to the particles on the equilibrium's boundary is 0.5742 which is within 3% of the geodesic radius computed in (71) (R ≈ 0.5570).
We have also performed numerical experiments with larger numbers of particles, to confirm that the particle method gives better numerical approximations of the continuum equilibrium as the number of particles increases. Indeed, with N = 400 and N = 900 particles we found symmetric equilibria supported in geodesic disks of radii 0.5692 and 0.5661, respectively, which are within 2.1% and 1.6% of the continuum geodesic radius in (71).
Arbitrary initial data. Our numerical simulations have indicated that all (not just the symmetric) solutions to model (1) on H, with interaction potential given by (58), approach asymptotically constant equilibrium densities supported on geodesic disks of radius R given by (71). For the model in Euclidean space the analogous result was illustrated and proved in [32,10]. To have such a result hold for the model on the hyperbolic plane suggests that the construction of the interaction potential in (58) can be potentially extended, with similar outcomes, to other non-compact manifolds. Figure 6 corresponds to a numerical simulation with N = 100 particles initiated at a randomly-generated configuration on H, with the θ and φ coordinates of the particles drawn randomly in the intervals (0.3, 2.3) and (0, π/2), respectively. Figure  6(a) shows the initial configuration and Figure 6(b) shows the equilibrium state, along with some particle trajectories (see also the zoom-in plot in Figure 6(c)). The equilibrium density is constant within its support (see (62) and (63)), which for a particle simulation translates into a uniform distribution (with respect to the metric on H). We found very similar results with larger number of particles and different initialization procedures. Equilibria obtained from simulations with N = 400 and N = 900 will be discussed below.
Next we provide evidence that the support of the equilibrium is a geodesic disk. We illustrate the procedure for the equilibrium in Figure 6  equilibrium's shape, we calculate numerically the Riemannian centre of mass of the equilibrium configuration. The Riemannian centre of mass is guaranteed to exist in such case, given that H has negative curvature everywhere [2]. To locate the centre of mass we use the intrinsic gradient descent algorithm investigated in [3]; recall that the Riemannian centre of mass of a set of points on a manifold minimizes the sum of squares of the geodesic distances to the data points, so using a gradient decent method is very natural here. Figure 7(a) shows the centre of mass C (red diamond) of the equilibrium configuration in Figure 6(c) computed with this method.
After we locate the particles on the equilibrium's boundary (see filled blue circles in Figure 7(a)) we compute the geodesic distances between the Riemannian centre of mass and the boundary points. We denote these distances by R i (note that the particles on the boundary have been relabelled so that they have consecutive indices starting from 1) -see Figure 7(a) for an illustration. Figure 7(b) shows these distances for three simulations. For the simulation with N = 100 discussed above, there are 31 particles on the boundary. Their distances to the Riemannian centre of mass of the equilibrium are shown in the figure as magenta circles connected by dotted lines, where the thick dotted line represents their mean value. The mean value is 0.5109 with a relative standard deviation of 0.55%. The blue circles connected by dash-dotted lines correspond to a simulation with N = 400 particles (67 particles on the equilibrium's boundary). The mean value (thick dash-dotted line) is 0.5330 with a relative standard deviation of 0.47%. Finally, the red circles connected by dashed lines represent the distances from the centre to the boundary points for a simulation with N = 900 particles (98 of which on the boundary). Their mean value (thick dashed line) is 0.5415 with a relative standard deviation of 0.20%.
The results presented in Figure 7 strongly suggest that the continuum equilibrium is supported on a geodesic disk of radius R given by (71) (this value has been indicated as a thick black solid line in the figure; recall R ≈ 0.5570). Indeed, as the number of particles increases, not only that the mean value of the distances R i approach R, but also their relative standard deviation decreases. In other words, the larger the number of particles, the closer the boundary of the particle equilibrium is to a geodesic circle of radius R. Similar results were obtained with a variety of initial configurations and different number of particles.
We conclude this section by posing the following conjecture: Geodesic disks of constant density are global attractors for model (1) on H with the interaction potential (58). Recall that the analogous statement holds for the model in Euclidean space with potential (13) [32,10]. One major obstacle for proving such a global convergence result is that the aggregation model on manifolds (and in particular on H) does not necessarily conserve the Riemannian centre of mass (see Section 2). In Euclidean spaces the centre of mass is conserved and the general proof in [10] relies fundamentally on this fact. Consequently, while the centre of the attracting geodesic disk is a priori known in Euclidean spaces (it is the centre of mass of the initial density), it is not known for the aggregation model on H. For this reason the proof for global attractors in the plane does not immediately extend to H.

5.
Other potentials and a new application. In this section we present numerical results for the model on the hyperbolic plane using other interaction potentials. We also present a completely new application, namely the aggregation model set up on the rotation group SO(3).

5.1.
Other interaction potentials. Potentials in power-law form have been frequently considered in the literature on the aggregation model (1) in Euclidean spaces [6,31,32,41,53]. In our context, by using the geodesic distance between points, a power-law potential reads: where the exponents p and q (with p < q) correspond to repulsive and attractive interactions, respectively. As shown in various works [31,41,53], the delicate balance between attraction and repulsion often leads to complex equilibrium configurations, supported on sets of  Figure 7) suggest that the equilibrium configuration consists of a uniform (with respect to the metric on H) particle distribution supported over a geodesic disk of radius R (see (71)).
various dimensions. The existence and characterization of equilibria as minimizers of the interaction energy has been an active research topic lately (existence was investigated in [22,50,15], the dimensionality of local minimizers was studied in [5]). Below, we consider the aggregation model (1) on the hyperbolic plane H and demonstrate numerically that a variety of equilibria can also be obtained with a power-law potential on such manifold. similar equilibria have been found for the aggregation model in Euclidean plane with power-law interaction potentials [31,41,22,5]. We also point out that to establish whether certain boundaries/shapes consist of geodesic circles, we followed a very similar procedure to that discussed in the previous section (see Figure 7). Namely, we located the Riemannian centre of mass of the equilibrium, along with the particles on the boundaries, and computed the mean distance (radius) from the centre of mass to these particles. In all cases, including the simulations presented below, the radii of these various circles had a relative standard deviation within 1% (in most cases much smaller in fact). Another class of interaction potentials which has been widely used in the literature on model (1) in R n consists of (generalized) Morse-type potentials [18]: where V (r) = −e − r s s , with s > 0.
(77) Here, C and l are positive constants, which control the relative size and range of the repulsive interactions. In one dimension with s = 1, V (r) is a multiple of the Green's function of the differential operator ∂ 2 r − Id. This property enables explicit calculations of the equilibrium solutions by converting integral equations for equilibria into differential equations [7]. Potentials of form (76)-(77) have been also used in other models for swarming and flocking [23].
The aggregation model (1) in the Euclidean plane with Morse-type interaction potentials exhibits a wide variety of possible equilibria [18]. In Figure 8(d)-(f) we showcase some equilibria obtained numerically for the model in the hyperbolic plane with potential (76)-(77). We picked values of the parameters C, l and s that have been used for the Euclidean model [18]. In each plot we observe mixed dimensionality of the equilibrium's support: (d) geodesic circle (ring) with a delta accumulation at the centre, (e) concentration on a ring with a continuous density supported on a concentric geodesic disk, (f) concentration on a ring with a continuous density inside. A rigorous investigation of equilibria shown in Figure 8 is challenging. One possible direction could be to study the stability of some of these equilibria (e.g., the ring, the annular region), extending similar analyses in Euclidean spaces [41,53,6]. Also, the results in Figure 8 suggest that the dimensionality of the equilibria on H relates to the strength of repulsion (value of the exponent p for potential (75)), as for the model in the Euclidean plane (see [5]). Investigating this connection further is also an interesting future direction that can be pursued.

5.2.
Aggregation model on the rotation group. Rotation group as a Riemannian manifold. We present briefly the Riemannian structure on the 3D rotation group SO(3), defined by The tangent space to SO(3) at a rotation R ∈ SO(3) is given by where so(3) is the space of 3 × 3 skew symmetric matrices. The Riemannian metric on the tangent space T R SO(3) is given by: (3), where ·, · F denotes the Frobenius inner product. Any rotation R ∈ SO(3) can be identified via the exponential map with a pair (θ, v) ∈ [0, π] × S 2 , where S 2 denotes the unit sphere in R 3 . The pair (θ, v) is referred to as the angle-axis representation of the rotation, where the unit vector v indicates the axis of rotation and θ represents the angle of rotation (by the righthand rule) about the axis. The representation of R in terms of (θ, v) is given by Rodrigues's formula. To list it we need the following common notation: Then, the angle-axis representation of a rotation R is: withv given by (78). The inverse of the representation (79), θv = log R, is given explicitly by: Here, exp and log represent the matrix exponential and logarithm, respectively. We also list here some standard facts about geodesics and the exponential map on the rotation group. Given two rotation matrices R 1 , R 2 ∈ SO(3), the shortest path between R 1 and R 2 is along the geodesic curve Q : [0, 1] → SO(3) given by: Note that Q (t) = Q(t) log(R T 1 R 2 ). From here, one infers easily that the Riemannian distance on SO(3) between R 1 and R 2 is where θ 12 v 12 = log(R T 1 R 2 ). From (81) one finds explicitly the exponential map at R 1 : and its inverse:

RAZVAN C. FETECAU AND BERIL ZHANG
We also note that, in particular, the formulas (79) and (80) can be expressed using the exponential map at the identity I. Indeed, R = exp I (θv), θv = exp −1 I R. Aggregation model on SO(3). We consider now the aggregation model on the Riemannian manifold SO (3). For the purpose of numerical simulations, similar to what we did in Sections 3.3 and 4.3 for the sphere and the hyperbolic plane, we consider the evolution of N rotation matrices R i (t) on SO(3), using a discrete form of (1b). The discrete particle system on M = SO(3) (see also (48) and (73)) reads: Note that in these considerations we may refer to a rotation matrix as "particle", and also use M to denote SO(3) when notations would become cumbersome otherwise. Given that the interaction potential depends only on d(R i , R j ), the distance in SO(3) between R i and R j , one has: where we used the explicit form of the gradient of the distance function in terms of the exponential map [47]. By (82) and (83), one finds from (84) and (85) the following particle system: where θ ij v ij = log(R T i R j ). We solve numerically the ODE system (86) for the angle-axis representation θ ivi = log R i , i = 1, . . . , N , with the 4th order Runge-Kutta method. In Figure  9 we show results for several simulations using N = 100 particles and power-law potentials. The rotation matrices R i are plotted in angle-axis representation: the angles θ i are shown on the left and the unit vectors v i are shown on the right plots, respectively.
For plots in Figure 9(a)-(b) the initial values for θ i were drawn randomly in the interval (0, π/2), while the unit vectors v i were generated in spherical coordinates, with the polar angles drawn randomly in the interval (π/6, π/2) and the azimuthal angles equispaced in (0, π/2). For this simulation we used an attractive-repulsive power-law potential with p = 5 and q = 10 (see (75)). In the figure we show the states at initial time and at time t = 40. At time t = 40 the configuration has evolved into an aggregation of four points on SO(3) (see the red stars indicating four values of θ and four locations on the unit sphere). Note that we sorted the angles of the final state for a better visualization.
The simulation in Figure 9(c)-(d) is for a purely attractive interaction potential in power-law form with q = 2 (see (75), but with no repulsion term). The initial state of N rotation matrices was randomly generated across the entire space SO(3), in the following sense: the angles θ i were drawn randomly in the interval (0, π), and the polar and the azimuthal angles of v i were drawn randomly in the intervals (0, π) and (0, 2π), respectively. In other words, all variables were drawn randomly from within their full ranges. The final state is an aggregation in one point on SO(3) (see the red stars indicating one value of θ and one unit vector), corresponding to a full synchronization of the rotation matrices. Achieving consensus in a network of rigidbodies (modelled as a multi-agent system on SO (3)) is a very important problem in the robotic control community [49,51]. In future work we plan further investigations of consensus on SO(3), hoping to provide new insights and perspectives (from a dynamical systems approach) to this problem. In closing, we believe that the aggregation model on general manifolds that has been proposed in this paper, the general construction of an interaction potential that leads to constant density equilibria on the sphere and the hyperbolic plane (along with the analytical considerations that can be made in such cases), as well as the various numerical illustrations that demonstrated swarming with other interaction potentials, will have set up a framework for and motivate further research and developments on self-organization models and their applications.