Quaternions in collective dynamics

We introduce a model of multi-agent dynamics for self-organised motion; individuals travel at a constant speed while trying to adopt the averaged body attitude of their neighbours. The body attitudes are represented through unitary quaternions. We prove the correspondance with the model presented in a previous work by the three first authors where the body attitudes are represented by rotation matrices. Differently from this previous work, the individual based model (IBM) introduced here is based on nematic (rather than polar) alignment. From the IBM, the kinetic and macroscopic equations are derived. The benefit of this approach is twofold: firstly, it allows for a better understanding of the macroscopic equations obtained and, secondly, these equations are prone to numerical studies, which is key for applications.

. The reference frame given by the orthonormal basis \{ e \bfone , e \bftwo , e \bfthree \} is rotated according to the unitary quaternion q and gives a new frame \{ e \bfone (q), e \bftwo (q), e \bfthree (q)\} describing the body attitude of the bird (where e \bfi (q) denotes the rotation of the vector e \bfi by q \in \BbbH 1 ). The direction of motion is given by the vector e \bfone (q), while the pair (e \bftwo (q), e \bfthree (q)) gives the relative position of the body around this direction.
Agents move at a constant speed while attempting to coordinate their body attitude with those of near neighbors. Here we present an individual-based model (particle model) for body attitude coordination. We derive the corresponding macroscopic equations from the associated mean-field equation, which we refer to as the Self-Organized Hydrodynamics based on Quaternions (SOHQ), by reference to the Self-Organized Hydrodynamics (SOH) derived from the Vicsek dynamics (see [21] and the discussion below). Our model is inspired by the one in [17] where the body attitude is represented by elements of the rotation group SO (3). The macroscopic equations obtained in [17] present various drawbacks. First, some of the terms in the equations do not have a clear interpretation. Second, and most importantly, the macroscopic equations are impractical for numerical simulations due to their complexity, especially since some terms are defined implicitly (see (16)).
In the future, we aim to investigate and compare numerically the microscopic and macroscopic dynamics. The individual-based dynamics require great computational power given the large number of agents in the system. A quaternion formulation has better computational performance than a matrix formulation in terms of storage and computational efficiency. In terms of memory, quaternions only require 4 units of memory, while matrices require 9 units; i.e., for a given amount of memory, we can store more than double the number of agents with quaternions than with matrices. In terms of computational efficiency, at each time step, due to the accumulation of numerical error, we will need to orthogonalize the matrices to make sure that they belong to SO (3), which is computationally expensive [13]. The equivalent operation for quaternions corresponds to just normalizing a 4-dimensional vector. Even though the rotation of a vector by a unitary quaternion is more expensive than by a matrix (17 additions and 24 multiplications versus 6 additions and 9 multiplications [26]), the savings in memory and in computational power required for the other operations largely outweigh this advantage of the matrix representation. Also, more specifically to our models, the individual-based model for matrices requires a polar decomposition to be performed at each time step for each agent (as will be seen in (101)). This is in general computed using the singular value decomposition of the matrix, which makes this process computationally expensive [36]. A new method for computing the polar decomposition of a matrix is proposed in [33], which claims to be more efficient. This method is based on computing the corresponding Q matrix in quaternion form (see (20) below) and computing the leading eigenvector of this matrix. With the quaternion formulation, we do this directly without adding the extra computational costs required to convert back and forth into the matrix formulation. Analogously, the simulation of the macroscopic equations also benefits from the computational efficiency of the quaternions: each grid point will require less than half the memory to store the information of the system than is required when using matrices; we will avoid orthogonalizing matrices by just normalizing 4-dimensional vectors; and, importantly, all the terms in the quaternion formulation are explicit, while in the matrix formulation they are not and require the inversion of an operator (see (16)). In summary, by considering a quaternionic representation, we render the results in [17], which uses the matrix representation, amenable to a numerical study.
In contrast with the use of rotation matrices in [17], the use of the quaternion representation makes the modeling more difficult on the individual-based level; first, it is not clear how to define a mean body attitude based on quaternions and, second, we need to consider nematic alignment rather than polar alignment. However, the macroscopic equations obtained are easier to interpret than in [17] and provide the right framework to carry out numerical simulations. The main contributions of the present paper include (i) deriving the macroscopic equations, (ii) finding the right modeling at the individual-based level, and (iii) proving the equivalence of the models and results obtained here with those in [17] for the rotation-matrix representation.
There exist already a variety of models for collective behavior depending on the type of interaction between agents. In the case of body attitude coordination, apart from [17], other models have been proposed; see [40] and references therein. This has applications in the study of collective motion of biological systems such as sperm dynamics, or of animals such as birds and fish, and it is a stepping stone to modeling more complex agents formed by articulated bodies (corpora) [10,11]. In the rest of the section we present related results in the literature and the structure of the document.
The literature on collective behavior is extensive. Such systems are ubiquitous in nature: fish schools, flocks of birds, herds [6,7,37], bacteria [4,45], and human walking behavior [32] are some examples. The main benefit to studying collective motion and self-organization is to gain an understanding of their emergent properties: local interactions between a large number of agents give rise to large scale structures (see the review in [44]). Given the large number of agents, a statistical description of the system is more pertinent than an agent-based one. With this in mind, mean-field limits are devised when the number of agents tend to infinity. From them, macroscopic equations can be obtained using hydrodynamic limit techniques (as we explain below).
The body attitude coordination model presented here and the one in [17] are inspired by the Vicsek model. The Vicsek model is a particular type of model for self-propelled particles [1,12,30,43] where agents travel at a constant speed while attempting to align their direction of motion with their neighbors. Other refinements and adaptations of the Vicsek model (at the particle level) or the SOH model (at the continuum level) have been proposed in the literature; we just mention a couple as examples: in [8] an individual-based model is proposed to better describe collective motion of turning birds; in [22] agents are considered to have the shape of discs, and volume exclusion is included in the dynamics.
One key difference in the modeling with respect to [21] is that we consider nematic alignment rather than polar alignment: given q \in \BbbH 1 , q and - q represent the same rotation. Collective dynamics based on nematic alignment is not, however, new; see, for example, [20,22] and references therein. Nematic alignment also appears extensively in the literature of liquid crystals and colloids, like suspensions of polymers; see [20] and the reference book [25]. Our results are inspired by the SOH model (the continuum version of the Vicsek model) presented in [21], where we have substituted velocity alignment by body attitude coordination. The macroscopic equations are obtained from the mean-field limit equation, which takes the form of a Fokker--Planck equation.
To obtain the macroscopic equations, the authors in [21] use the well-known tools of hydrodynamic limits, first developed in the framework of the Boltzmann equation for rarefied gases [9,14,41]. Since its first appearance, hydrodynamic limits have been used in other contexts, including traffic flow modeling [3,31] and supply chain research [2,23]. However, in [21] a methodological breakthrough is introduced: the Generalized Collision Invariant (GCI), which will be essential to the present study (section 4.3). Typically, to obtain the macroscopic equations we would require as many conserved quantities in the kinetic equation as the dimension of the equilibria (see again [44]). In the mean-field limit of the Vicsek model this requirement is not fulfilled, and the GCI is used to sort out this problem. For other cases where the GCI concept has been used, see [15,16,17,18,19,24,27].
After this introduction, we discuss the main results in section 2. In section 3 we explain the derivation of the individual-based model for body coordination dynamics and show its equivalence to the model in [17] in section 5.2. Then in section 3.2 we give the corresponding (formal) mean-field limit for the evolution of the empirical measure when the number of agents goes to infinity.
The following part concerns the derivation of the macroscopic equations (Theorem 4.1) for the macroscopic density of the particles \rho = \rho (t, x) and the quaternion of the mean body attitude \= q = \= q(t, x). To obtain these equations we first study the rescaled mean-field equation ((27) in section 4.1), which is, at leading order, a Fokker--Planck equation. We determine its equilibria (see (31)). In section 4.3 we obtain the GCIs (Proposition 4.12), which are the main tools used to derive the macroscopic equations in section 4.4. Finally, in section 5 we prove the equivalence of our equations and results with the ones obtained in [17].
2. Discussion of the main results.
The quaternions q and - q represent the same rotation. The product of unitary quaternions corresponds to the composition of rotations. The interested reader can find further information on the theory of quaternions in [34,38] Remark 2.1 (identification between purely imaginary quaternions and vectors in \BbbR 3 ). Notice that in (3) we abuse notation: the product qvq \ast must be understood in the quaternion sense (therefore we consider v = (v 1 , v 2 , v 3 ) as a quaternion which is purely imaginary, i.e., v = v 1 \vec{} \imath + v 2 \vec{} \jmath + v 3 \vec{} k). Conversely, \= v is understood as a vector in \BbbR 3 rather than a purely imaginary quaternion. This abuse of notation where we identify vectors in \BbbR 3 with purely imaginary quaternions (and the converse) will be used throughout the text. The sense should be clear from the context. We will also use in general q to denote a unitary quaternion and p to denote an arbitrary quaternion.
2.2. Self-organized hydrodynamics based on quaternions (SOHQ). In section 3 we introduce an individual-based model for collective dynamics where individuals are described by their location in space and the position of their body (body attitude). Individuals move at a constant speed while trying to adopt the same body attitude, up to some noise; see (24)-- (25). The body attitude is given by three orthonormal vectors, where one of the vectors indicates the direction of motion and the other two represent the relative position of the body around this direction. In this manner, the body frame of a given individual is characterized by the rotation of a fixed reference frame. This rotation will be represented here by elements of the group of unitary quaternions, denoted by \BbbH 1 (see Figure 1). The main result of this paper is Theorem 4.1, which gives the macroscopic equations for these dynamics, i.e., the time-evolution equations for the macroscopic mass of agents \rho = \rho (t, x) and the mean quaternion \= q = \= q(t, x), which corresponds, as explained, to the mean body attitude. Here t \geq 0 is the time and x \in \BbbR 3 denotes a point of the physical space. We will refer to this system as the Self-Organized Hydrodynamics based on Quaternions (SOHQ).
With this notation the SOHQ corresponds to where e \bfone is a vector in \BbbR 3 and e \bfone (\= q) denotes the rotation of e \bfone by the quaternion \= q, that is, and where we used the (right) relative space differential operators where ((\partial xi \= q)\= q \ast ) i indicates the ith component of (\partial xi \= q)\= q \ast . In (7) and in the last three terms of (6) we use the abuse of notation explained in Remark 2.1. The matrix product in the fourth term of (6) has to be understood as a matrix product, giving rise to a scalar product in \BbbH : In (5)--(6), c 1 , c 2 , c 3 , and c 4 are explicit constants (given in Theorem 4.1) that depend on the parameters of the model, namely, the rate of coordination and the level of the noise. The constants c 2 , c 3 , and c 4 depend on the Generalized Collision Invariant (GCI; see the introduction and section 4.3). Interestingly, c 1 had a special meaning as a``(polar) order parameter"" in [17,21] (see Remark 4.14). Here it has the same meaning, but as a``nematic"" order parameter.
Equation (5) gives the continuity equation for the mass \rho and ensures mass conservation. The convection velocity is given by c 1 e \bfone (\= q), where the direction is given by e \bfone (\= q), a unitary vector (since e \bfone is unitary), and the speed is c 1 . Notice that the convection term is quadratic in \= q. This is a new structure with respect to [17,21]. We consider next (6) for \= q. Observe first that all the terms in the equation belong to the tangent space at \= q in \BbbH 1 , i.e., to \= q \bot . This is true for the first term since (\partial t + c 2 (e \bfone (\= q) \cdot \nabla x ) is a differential operator (giving the transport of \= q), and it also holds for the rest of the terms since they are of the form u\= q with u purely imaginary (see Proposition A.2 in the appendices).
The term corresponding to c 3 gives the influence of \nabla x \rho (pressure gradient) on the body attitude \= q. It has the effect of rotating the body around the vector directed by e \bfone (\= q) \times \nabla x \rho at an angular speed given by c3 \rho \| e \bfone (\= q) \times \nabla x \rho \| , so as to align e \bfone (\= q) with - \nabla x \rho . Indeed the solution to the differential equation when u is a constant purely imaginary unitary quaternion and \gamma a constant scalar, is given by q(t) = exp( - \gamma ut)q(0), and exp( - \gamma ut) is the rotation of axis u and angle - \gamma t (see (2)). Since c 3 is positive, the influence of this term consists of relaxing the direction of movement e \bfone (\= q) toward - \nabla x \rho , i.e., making the agents move from places of high concentration to low concentration. In this manner, the \nabla x \rho term has the same effect as a pressure gradient in classical hydrodynamics. In the present case the pressure gradient provokes a change in the full body attitude \= q. Finally, notice that in regions where \rho > 0 we can divide (6) by \rho , and this gives us the influence of each term depending on the local density. After division by \rho , we observe that the only term depending on the density \rho in (6) is the third term in the form Therefore, for small densities, this term may take large values and become dominant, while for large densities it becomes small and the other terms in the equation prevail for reasonably large \nabla x \rho . The fact that agents tend to relax their direction of motion toward regions of lower concentration creates dispersion; this term is a consequence of the noise at the microscopic level. However, the relaxation becomes weaker once agents are in regions of high density or areas with small variations of density. The last two terms in (6) are unique to the body attitude coordination model and are the main difference with respect to the SOH equations for the Vicsek model. Analogously to the discussions in [17,21] for the body attitude model based on rotation matrices and for the Vicsek model, the SOHQ model bears some similarities to the compressible Euler equations, where (5) is the mass conservation equation and (6) is akin to the momentum conservation equation, where momentum transport is balanced by a pressure force. There are, however, major differences. First, the pressure term belongs to \= q \bot in order to ensure that \= q \in \BbbH 1 for all times; in the Euler equations the velocity is an arbitrary vector, not necessarily normalized. Second, the convection speed c 2 is a priori different from the mass conservation speed c 1 . This difference signals the lack of Galilean invariance of the system, which is a common feature of all dry active matter models (models for collective motion not taking place in a fluid); see [42]. Finally, the last two terms of (7) do not have a clear analogue to the compressible Euler equations: they seem quite specific to our model. and the interpretation in terms of a relative variation is standard.
Remark 2.2. When \partial = \partial t is the time derivative, for a function q = q(t, x) with values in \BbbH 1 , the vector \partial t,rel q = \partial t q q - 1 is half of the angular velocity of a solid of orientation represented by q. By analogy, the vector \partial xi,rel q = \partial xi q q - 1 for i = 1, 2, 3 is half of the angular variation in space of a solid of orientation represented by q.
Multiplying from the right the evolution equation (6) for \= q by \= q \ast , we obtain the following equivalent equation: In this equation we notice that all the differential operators naturally appear under their (right) relative form. Notice also that all other nonlinearities in \= q are expressed in terms of e \bfone (\= q). Therefore, the previous system can be interpreted as the evolution of the relative changes of \= q. In terms of r, the previous system can be recast into \Bigl[ \rho \partial t r + \rho c 2 (e \bfone (q) \cdot \nabla x )r + c 3 e \bfone (q) \times \nabla x \rho + c 4 \rho [\nabla x r e \bfone (q) + (\nabla x \cdot r)e \bfone (q)] \Bigr] \bigm| \bigm| \bigm| (t0,x0) = 0.
For an interpretation (again, local ) in terms of b := log r we refer the reader to section 5.3.

2.3.
Equivalence with the previous body attitude model. In [17] a model for body attitude coordination is presented where the body attitude is represented by a rotation matrix (element in SO(3), orthonormal group) rather than by a quaternion. In section 5.2 we will prove the equivalence between the individual-based model presented in [17] and the one here, in the sense that the two stochastic processes are the same in law (Theorem 5.6).
In [17] also the macroscopic equations are obtained for the mean body attitude \Lambda = \Lambda (t, x) \in SO(3) and spatial density of agents \rho = \rho (t, x) \geq 0, called self-organized hydrodynamics for body attitude coordination (SOHB): \rho \Bigl( \partial t \Lambda + \c 2 \bigl( (\Lambda e \bfone ) \cdot \nabla x \bigr) \Lambda \Bigr) \Lambda t + \bigl[ (\Lambda e \bfone ) \times \bigl( \c 3 \nabla x \rho + \c 4 \rho r x (\Lambda ) \bigr) + \c 4 \rho \delta x (\Lambda )\Lambda e \bfone \bigr] \times = 0, with explicit constants \c i , i = 1, . . . , 4, where \Lambda t indicates the transpose matrix of \Lambda , and where for a vector u = (u 1 , u 2 , u 3 ) the antisymmetric matrix [u] \times is defined by (15) [u] \times := The scalar \delta x (\Lambda ) and the vector r x (\Lambda ) are first order differential operators intrinsic to the dynamics. We define them next. For a smooth function \Lambda = \Lambda (x) from \BbbR 3 to SO(3), we define the matrix D x (\Lambda ) by the equality (16) for all w \in \BbbR 3 , (w \cdot \nabla x )\Lambda = [D x (\Lambda )w] \times \Lambda (the matrix D x (\Lambda ) is well defined; see [17]). The first order operators \delta x (\Lambda ) and r x (\Lambda ) are then defined by Since the individual-based model formulated in terms of quaternions is equivalent (in law) to the one formulated with rotation matrices, we expect their respective macroscopic limits to also be equivalent. This is the case, as expressed in Theorem 5.13; i.e., if at time t = 0 \Lambda (0) and \= q(0) represent the same rotation, then \Lambda (t) and \= q(t) represent the same rotation for all t where the solutions are well defined. There are, however, important differences between the SOHB and SOHQ macroscopic equations. On one hand, notice that the operators \delta x and r x cannot be expressed under a simple explicit form, which makes the meaning of these operators less clear. In the quaternion case, all the elements in (5)--(6) are explicit. Moreover, quaternions give the right framework for numerical simulations (in terms of memory and operation efficiency), as explained in the introduction. On the other hand, when using rotation matrices, we obtain clear equations for the evolution of each one of the orthonormal vectors that define the body frame (see [17]). However, the expressions for these vectors in the quaternion formulation is complicated and not very revealing, due to the quadratic structure of the rotation (see (7)).

Modeling:
The individual-based model and its mean-field limit.
3.1. The individual-based model. Consider a reference frame in \BbbR 3 given by the orthonormal basis \{ e \bfone , e \bftwo , e \bfthree \} . Consider also N agents labeled by k = 1, . . . , N with position X k (t) \in \BbbR 3 and body attitude given by the unitary quaternion q k (t) \in \BbbH 1 . As explained in the introduction, the body frame of agent k corresponds to \{ e \bfone (q k ), e \bftwo (q k ), e \bfthree (q k )\} , where e \bfi (q k ) denotes the rotation of e \bfi by q k for i = 1, 2, 3. The vector e \bfone (q k ) gives the direction of motion of agent k, and the other two vectors give the position of the body relative to this direction (see Figure 1). Our goal is to model collective dynamics where agents move at a constant speed while adopting the body attitude of their neighbors, up to some noise.
Evolution of the positions, (X k ) k=1,...,N . The fact that agent k moves in direction e \bfone (q k ) at constant speed v 0 > 0 simply corresponds to the equation dX k dt = v 0 e \bfone (q k ), e \bfone (q k ) := Im(q k e \bfone q \ast k ) (recall (3) and Remark 2.1). Notice that the speed for all agents is constant and equal to v 0 > 0.
Evolution for the body attitudes, (q k ) k=1,...,N . Agents try to coordinate their body attitudes with those of their neighbors. To model this phenomenon, we need to first define an``average"" body attitude around a given agent k and, second, express the relaxation of the body attitude of agent k towards this average. Remark 3.1 (nematic alignment, sign invariance). The body attitude is uniquely defined by quaternions up to a sign since q k and - q k represent the same rotation. This implies, first, that the time evolution equation for q k = q k (t) must be sign invariant and, second, that the average must take this sign invariance into account. This is called``nematic alignment"" (as opposed to``polar"" alignment), and it appears in other collective models [20,22] and in liquid crystals [25]. Therefore, we cannot define the average analogously as in [21] since the alignment is polar in this case. For example, if one considers the normalized averaged quaternion defined in the same way as in the Vicsek model, (18) we obtain a unitary quaternion that can be interpreted as a rotation. However, the meaning of this rotation is unclear, and it is not invariant under changes of sign of any of the vectors q i . Nor can we use the nematic average used in [20] since it is valid only in \BbbR 2 .
We define (up to a sign) the average around q k by \= q k := unitary eigenvector of the maximal eigenvalue of Q k = arg max\{ q \cdot Q k q, q \in \BbbH 1 \} , with where the nonnegative-valued function K is a kernel of influence. It is in the definition that \= q k \in \BbbH 1 ; one can check that if \= q k is an average, so is - \= q k (so it is sign invariant); \= q k remains invariant under the change of sign of any of the arguments q 1 , . . . , q N ; and \= q k maximizes over \BbbH 1 : where \widehat \langle q i , q\rangle denotes the angle between q i and q (seen as elements of the hypersphere \BbbS 3 ). Now to express the relaxation of q k toward this average we define first and write (23) dq where P \bfq \bot \bfk indicates the projection on the orthogonal space to q k (which corresponds to the tangent space of q k in \BbbH 1 ); \nabla \bfq k indicates the gradient in \BbbH 1 (the second equality in (23) is proven in Proposition A.1). Equation (23) relaxes q k toward the maximizer of q \mapsto \rightar q \cdot \bigl( \= q k \otimes \= q k -1 4 Id \bigr) q, which corresponds precisely to \= q k or - \= q k . Finally, putting everything together, we obtain the evolution equations: where \nu > 0 indicates the intensity of the relaxation. The evolution for the body attitudes results from two competing phenomena: body attitude coordination (the F k -term) and noise due to errors that the agents make when trying to coordinate. The noise term is given by (B k t ) k=1,...,N independent 4-dimensional Brownian motions. The constant D > 0 gives the intensity of the noise term (by modifying the variance of the Brownian motion). The projection P \bfq \bot k in (25) ensures that q k (t) \in \BbbH 1 for all times t where the solution is defined. The stochastic differential equation (24)--(25) must be understood in the Stratonovich sense; see [35].
Remark 3.2. Some comments: (i) Typically the noise term would be scaled by \surd 2D, because then the generator of the process is the Laplacian with coefficient D. However, we chose the scaling \sqrt{} D/2 to make the model equivalent with the one based on rotation matrices in [17]. This will be discussed in section 5.2.
(ii) The operator Q k in (20) corresponds to the de Gennes Q-tensor that appears in the theory of liquid crystal [25] and which is also related to the so-called nematic order coefficient. Notice that in the definition of Q k in (20) the 1/4-factor could be ignored and the definition of the average \= q k would remain unchanged. Also in F k in (22) the 1/4-factor can be ignored since that term disappears in the projection in (23). We keep it here for the parallelisms that it bears with the theory of liquid crystal and because it will appear when we define the equilibrium distribution in (31). (iii) Notice that to define the average \= q k in (21), we assume that the maximal eigenvalue is simple. At the formal level, this assumption is reasonable since for general symmetric matrices the event of a multiple maximal eigenvalue is negligible. Of course, for a rigorous analysis, we would need to ensure carefully that this event can actually be neglected, or we would need to add an extra rule to determine uniquely the average. (iv) Notice that we could have defined the relaxation F k by considering directly F k = Q k q k instead of (22), since, in this case, (23) relaxes also to \= q k . However, for this case the relaxation is weaker. This is a modeling choice. We will prove in section 5.2 that our choice here is the one that corresponds to the model presented in [17], where the body attitude is described with rotation matrices.
(v) One can check that the particle system (24)--(25) is frame invariant in the sense that if \X k = R frame (X k ) for k = 1, . . . , N , with R frame the frame change associated to the quaternion q frame , and \q k = q frame q k , then the pair ( \X k , \q k ) satisfies the same system (with the appropriate initial conditions).

Mean-field limit.
We now obtain formally the mean-field limit for (24)--(25) as the number of particles N \rightar \infty . The rigorous mean-field limit has been proven for the Vicsek model in [5]. A key difference between the Vicsek model and the system (24)--(25) is the way we compute the average in (21). Consider the empirical distribution in (x, q) \in \BbbR 3 \times \BbbH 1 over time ..,N satisfy (24)-- (25). Assume that f N converges weakly to f = f (t, x, q) as N \rightar \infty . It is standard to show (formally) that f satisfies where \nabla \bfq and \Delta \bfq denote the gradient and the Laplacian in \BbbH 1 , respectively, and dq is the Lebesgue measure on \BbbH 1 .

Hydrodynamic limit.
The goal of this section is the derivation of the macroscopic equations for (26). After a dimensional analysis and a time and space scaling, described next in section 4.1, we recast (26) into with \= q f defined by where dq being the Lebesgue measure on \BbbH 1 . (Note that after the dimensional analysis and the rescaling the values of the parameters D and \nu as well as the variables t and x have changed; see section 4.1 for details.) We then analyze in section 4.2 the collisional operator \Gamma in (28); particularly, we determine its (von Mises--like) equilibria, given by (for \= q \in \BbbH 1 ) where d is a parameter given by and where Z is a normalizing constant (such that \int \BbbH 1 M \= \bfq dq = 1). We then describe the structure of the Generalized Collision Invariants (GCIs) for \Gamma in section 4.3.
With this information we are ready to prove our main result.
Theorem 4.1 ((formal) macroscopic limit). When \varepsi \rightar 0 in the kinetic equation (27) it holds (formally) that Moreover, if the convergence is strong enough and the functions \= q and \rho are regular enough, then they satisfy the system (5)--(6) that we recall here: where the (right) relative differential operator \nabla x,rel is defined in section 2.2, where and where c i , i = 1, . . . , 4, are explicit constants. To define them we use the following notation: For two real functions g, w consider (36) \langle g\rangle w := Then the constants are given by where d is given in (32), where and where h is the solution of the differential equation (64).
We recall that in section 2 we provided a discussion of this main result. The proof is given in section 4.4. We conclude this study with section 5.3, where we compare the macroscopic limit obtained here with the corresponding one for the body attitude model with rotation matrices from [17].

Scaling and expansion.
We assume that the kernel of influence K is Lipschitz, bounded, and such that We express the kinetic equation (26) in dimensionless variables. Let \nu 0 be the typical interaction frequency, i.e., \nu = \nu 0 \nu \prime with \nu \prime = \scrO (1). We consider also the typical time and space scales t 0 , x 0 with t 0 = \nu - 1 Skipping the primes, we get the same equation as (26) except that v 0 = 1, all the quantities are dimensionless, and D, \nu , and K are assumed to be of order 1. Notice, in particular, that \nu \prime D \prime is the same before and after the dimensional analysis.
To perform the macroscopic limit we rescale space and time by x \prime = \varepsi x and t \prime = \varepsi t. After skipping the primes we obtain Proof. The result is obtained by a Taylor expansion in \varepsi and by using that (recall (42)) In particular, we have Proof. Let \lambda \varepsi max , respectively, \lambda max , be the maximal eigenvalue of Q \varepsi f , respectively, Q f (we assume them to be uniquely defined). From Lemma 4.2, we have and, multiplying by \= q \varepsi f on both sides, By maximality of \lambda max , we have that \= By symmetry, we also have \lambda max \leq \lambda \varepsi max + \scrO (\varepsi 2 ), and therefore On the other hand, we have that By our crucial assumption that \lambda max is a single eigenvalue with eigenvector \= q f , we can invert the matrix (Q f - \lambda max Id) on the 3-dimensional space \= q \bot f . By a small abuse of notation we write its inverse which proves (45). Taking the scalar product with \= q \varepsi f and using the fact that \= q \varepsi f is unitary, we have that 1 -(\= q f \cdot \= q \varepsi f ) 2 = \scrO (\varepsi 2 ), so using (45) we can finally show that Using Proposition 4.3 we recast the rescaled kinetic equation (43) as (27).

Equilibrium solutions and
Fokker--Planck formulation. Define d = D/\nu and consider the generalization of the von Mises distribution in \BbbH 1 : where Z is a normalizing constant. Observe that Z < \infty is independent of \= q since the volume element in \BbbH 1 is left-invariant, i.e., Note that we can recast \int \pi 0 m(\theta \prime ) sin 2 (\theta \prime /2)d\theta \prime , with \= q \cdot q = cos(\theta /2), with m(\theta ) given by (41). Indeed, and by Proposition A.3, Proposition 4.4 (properties of \Gamma ). The following hold: (i) The operator \Gamma in (28) can be written as and we have where \rho is the macroscopic mass, i.e., and \= q is the eigenvector corresponding to the maximum eigenvalue of Furthermore, H(f ) = 0 if and only if f = \rho M \= \bfq for some \rho \geq 0 and \= q \in \BbbH 1 .
Remark 4.5 (comparison with the equilibria considered in [17]). Thanks to (95), one can check that the equilibria M \= \bfq represent the same equilibria as for the kinetic model corresponding to the body attitude model with rotation matrices in [17], which is given by (where Z \prime is a normalizing constant); i.e., as long as \Lambda = \Phi (\= q) and A = \Phi (q) (\Phi is defined in (90)), we have Z \prime M \Lambda (A) = ZM \= \bfq (q). Note that the normalizing constants Z and Z \prime are not equal since the measures chosen on SO(3) and \BbbH 1 are identical only up to a multiplicative constant.
Proof of point (ii). From inequality (51) We are left with proving that \= q is the eigenvector corresponding to the maximum eigenvalue of \int (this will not change if multiplied by \rho since it is positive). For any quaternion p 0 \in \BbbH , the left multiplication by p 0 , that is, p \in \BbbH \mapsto \rightar p 0 p \in \BbbH , is an endomorphism on \BbbH . We write E l (p 0 ) the associated matrix, so that for all p \in \BbbH the (quaternion) product p 0 p is equal to the (matrix) product E l (p 0 )p. Using the change of variable q \prime = \= q \ast q, we compute To compute the value of the integral in the term above, first note that M 1 (q) depends only on Req. We use a change of variable that switch q i and q j (for i \not = j) to check that the off-diagonal terms (i, j) and (j, i) are zero. Then we compute the diagonal terms: the zeroth diagonal term is clearly given by \int \BbbH 1 (Re(q)) 2 M 1 (q) dq, while with the same changes of variable that switch q i and q j (for i \not = j) we check that the first to third diagonal terms are identical and equal to 1 3 \int \BbbH 1 Im 2 (q)M 1 (q) dq. Using the fact that Re 2 (q) + Im 2 (q) = 1, we obtain where we defined Note that for any p \in \BbbH , we have Therefore, equality (52) is a diagonalization of the matrix \int \BbbH 1 q \otimes qM \= \bfq (q) dq in an orthonormal basis. It is direct to check that \= q is an eigenvector corresponding to the first eigenvalue I 2 . It is the maximum eigenvalue if and only if that is, if and only if (54) We compute (using Proposition A.3) and, writing w(r) = (1 - by Jensen's inequality.
Therefore, we conclude that By the dominated convergence theorem, and using an integration by parts, we have lim d - \rightar\infty so that (54) holds true, which completes the proof. which particularly holds for \nabla \bfq \psi = 0, i.e., when \psi is a constant. Consequently, we only know one conserved quantity for our model corresponding to the macroscopic mass \rho . To obtain the macroscopic equation for \= q, a priori we would need three more conserved quantities. To sort out this problem, we use the method of the Generalized Collision Invariants (GCIs) introduced in [21].
Definition of the GCI. Define the operator for a function f and \= q \in \BbbH 1 . Notice that We write GCI(\= q) to denote the set of GCIs associated with \= q.
Characterization of the GCI. The main result of this section is the following description of the set of GCIs. where, for \beta \in \= q \bot , the function \psi \beta is defined by with h = h(r) the unique solution of the following differential equation on ( - 1, 1): Furthermore, the function h is odd, h( - r) = - h(r), and it satisfies for all r \geq 0 h(r) \leq 0.
This proposition will be crucial to compute the hydrodynamical limit in section 4.4. The proof is done in the two subsections below.
We define the Hilbert space (1 -r 2 ) 3/2 h 2 (r) dr < \infty and By a Lax--Milgram argument, we obtain the existence and uniqueness of the solution h in H ( - 1,1) . By uniqueness of h, we see that h is an odd function of its argument. By a maximum principle, we furthermore obtain that h(r) \leq 0 for r \geq 0.
To conclude, it only remains to show that the solution h corresponds to a function \psi \in H 1 0 (\BbbH 1 ). Since we know by Proposition 4.9 that the GCI exists and is unique in H 1 0 (\BbbH 1 ), this proves that \psi given by (63) is the GCI. For that, we first compute the L 2 (\BbbH 1 ) norm of gradient of \psi with We see directly that a sufficient condition for the first term of the last expression above to be finite is that Since Re(\beta \= q \ast ) = 0 (remember that by definition \beta \in \= q \bot ), the second term is finite as soon as the 3 \times 3 submatrix corresponding to the imaginary coordinate of the integral \int \BbbH 1 (q \otimes q)(1 -Re 2 q)(h \prime (Req)) 2 dq is finite, that is, as soon as \int \BbbH 1 (Imq \otimes Imq)(1 -Re 2 q)(h \prime (Req)) 2 dq < \infty , that is, when the diagonal terms are finite (the off-diagonal elements being null by the changes of variables that change the sign of the coordinate q i for i = 1, 2, 3), i.e., (1 -Re 2 q) 2 (h \prime (Req)) 2 dq < \infty .
Since h is in H ( - 1,1) , the second condition is true. Let us check the first condition.
By a computation similar to that for \nabla \bfq \psi , we see that which is null since \beta \in \= q \bot , so that \psi has mean zero on \BbbH 1 .
We are now ready to prove Proposition 4.7.
Proof of Proposition 4.7. The statement is a direct consequence of Proposition 4.9, Corollary 4.10, and Proposition 4.12.

The macroscopic limit.
This section is devoted to the proof of Theorem 4.1. We will use the following. where the positive constant c 1 is given in (37).
Remark 4.14 (comments on the constant c 1 ).
(i) The value for the constant c 1 obtained here is the same one as in the body attitude coordination model based on rotation matrices in [17]. This will allow us to prove the equivalence between the respective continuity equations (see section 5.3). (ii) In the case of the Vicsek model in [21] and in the body attitude coordination model based on rotation matrices in [17], the constant c 1 played the role of order parameter."" Particularly, it holds that c 1 \in [0, 1], and the larger its value, the more organized (coordinated/aligned) the dynamics are (and, conversely, the smaller c 1 , the more disordered the dynamics are). The extreme cases take place, for example, when D \rightar \infty , and then c 1 = 0, and when D \rightar 0, giving c 1 = 1. Here we have the same properties and interpretations for c 1 .
Using the changes of variable that switch the coordinates q i and q j for i \not = j, we see that the diagonal elements corresponding to q 2 i for i = 1, 2, 3 give rise to the same value for the integral, and therefore \biggr] e \bfone (\= q) M 1 (Re(q))dq, so that equality (74) holds for where we used Proposition (A.3) on the volume element.
We are now ready to prove Theorem 4.1.

Comparison with the results in
After an introductory presentation of the links between SO(3) and \BbbH 1 (section 5.1), we present the two main results of this section: the equivalence between the individual-based models (Theorem 5.6 in section 5.2) and the equivalence between the macroscopic models (Theorem 5.13 in section 5.3).

Relation between unitary quaternions and rotation matrices.
We first introduce some notation. Rotations in \BbbR 3 can be described mathematically in different ways. In this section we consider three particular descriptions, namely, the group of orthonormal matrices corresponding to the rotation group SO(3); the description via unitary quaternions q \in \BbbH 1 ; and, finally, rotations described by the pair (\theta , n) \in [0, \pi ] \times \BbbS 2 , where n indicates the axis of rotation and \theta the angle of rotation counterclockwise around n. For A \in SO(3), q \in \BbbH 1 , and (\theta , n) \in [0, \pi ] \times \BbbS 2 corresponding to the same rotation, we have the identities (Rodrigues's formula), q = q(\theta , n) = cos(\theta /2) + sin(\theta /2)n, where \= v \in \BbbH with Re(\= v) = 0 and Im(\= v) = v, and where we abuse notation in (88) and understand n as written in the Hamiltonian basis n = n 1 \vec{} \imath + n 2 \vec{} \jmath + n 3 \vec{} k rather than in the canonical basis n = (n 1 , n 2 , n 3 ). Notice that when \theta = 0, the vector n is not defined, but this does not pose a problem in the sense that there is an unambiguous correspondence with A = Id and q = 1.
Identities (ii) and (iv) are consequences of (89); identity (iii) is a consequence of (iv), noticing that q \ast = q - 1 and A t = A - 1 .
First, we show the relation between the inner products in SO(3) and \BbbH 1 .
Lemma 5.1. Let A, B \in SO(3) and q, r \in \BbbH 1 to be such that A = \Phi (q) and B = \Phi (r). Then where A \cdot B = Tr(AB t )/2, with Tr denoting the trace.
Next, we establish the correspondence between integrals in SO(3) and \BbbH 1 . Proof. We apply Proposition A.3 to the (even) function f (q) := g(\Phi (q)), to get  Finally, one can check that \Phi is continuously differentiable on \BbbH 1 , given that it is a quadratic function on \BbbH 1 . The following holds.

Equivalence between individual-based models.
In this section we check that the flocking dynamic considered in [17] corresponds with that of (24)-- (25).
In [17] the authors describe an individual-based model for body attitude coordination given by the evolution over time of (X k , A k ) k=1,...,N of N agents, where X k \in \BbbR 3 is the position of agent k and A k \in SO(3) is a rotation matrix giving its body attitude. The evolution of the system is given by the following equations: where the stochastic differential equation is in the Stratonovich sense (see [28]); W k t is the Brownian motion in the space of squared matrices; M k is defined as where K is a positive interaction kernel; \nu , v 0 , and D are positive constants; e \bfone is a vector; and P T A is the projection in SO(3) to the tangent space to A. The term P D(M ) denotes the orthogonal matrix obtained from the polar decomposition of M, which is defined as follows. The vector A k e \bfone in (100) gives the direction of movement of agent k and is obtained as the rotation of the vector e \bfone by A k . Equivalently, we can express it as A k e \bfone = e \bfone (q k ) (in the notation of (24)) as long as A k and q k represent the same rotation. Therefore, (24) and (100) represent the same dynamics, and we are left to check that q k = q k (t) and A k = A k (t) in (25) and (101) represent the same rotation for each time t where the solutions are defined.
The goal of this section will be to prove that the solution to the stochastic differential equation (24)-- (25) and the solution of the stochastic differential equation (100)--(101) are the same in law (in a precise way that will be given later).
The main result of this section is the following. The proof is done at the end of this section. First, we remark that in the absence of randomness (Brownian motion) the equations for the evolution of the body attitude are equivalent.
Proposition 5.7. Let A 0 \in SO(3) and q 0 \in \BbbH 1 represent the same rotation. Consider the matrix M k given in (102), the matrix Q k given in (20), and \= q k \in \BbbH 1 given in (19). Then, if det(M k ) > 0, the following two Cauchy problems are equivalent (in the sense that A k = A k (t) and q k = q k (t) represent the same rotation for all t where the solution is uniquely defined): Note that these two Cauchy problems can also be written, respectively, as where \nabla A and \nabla \bfq are the gradients in SO(3) and \BbbH 1 , respectively.
To prove this proposition we first check that the average orientation of the neighbors is the same in the two models, in the sense described below.
Proof. Assume for simplicity that K \equiv 1 (the general case can be proven equally). First, notice that for A = \Phi (q) (q \in \BbbH 1 ) it holds that for Q k given in (20) and where we used Lemma 5.1 to compute the inner product. Therefore, for any q \in \BbbH 1 , 2q \cdot Q k q = M k \cdot \Phi (q).
Now, the definition of \= q k implies that it maximizes q \mapsto \rightar q\cdot Q k q in \BbbH 1 . Since q\cdot Q k q = 1 2 M k \cdot A, this implies that \Phi (\= q k ) maximizes A \mapsto \rightar M k \cdot A in SO (3), which is a property that characterizes the matrix P D(M k ) (see [17,Prop. 3.1]). Therefore, it holds that \Phi (\= q k ) = P D(M k ), as long as det(M k ) > 0, and in this case, using again Lemma 5.1, we have We are now ready to prove Proposition 5.7.
Proof of Proposition 5.7. The fact that we can rewrite the first pair of Cauchy problems as the second pair comes from the equalities We conclude the equivalence thanks to Lemma 5.8 and Proposition B.1.
To prove Theorem 5.6 we need the following result.
Let \\sigma > 0, and let \H be a time-dependent tangent vector field on SO(3): Suppose that the following relations hold: Let \p t be the law over time of a stochastic process in SO(3) defined by dA = \H (A, t)dt + \\sigma P T A \circ d \B t for \B t a 9-dimensional Brownian motion. Then, if \p t is an absolutely continuous measure, the absolutely continuous measure p t defined by (108) \p t (\Phi (q)) = 2\pi 2 p t (q) for all q \in \BbbH 1 , t \geq 0, is the law over time of a stochastic process in \BbbH 1 defined by (109) dq = H(q, t)dt + \sigma P \bfq \bot \circ dB t for B t a 4-dimensional Brownian motion.
Proof. First, notice that for any Borel set B \subset \BbbH 1 it holds that thanks to Lemma 5.2. Note that this is the reason why we introduce the factor 2\pi 2 in (108), which allows us to have this equivalence of integrals. We start from the equation for \p t : Notice that the fact that we obtain a factor \\sigma 2 /4 is a consequence of considering the inner product A \cdot B = trace(A t B)/2 (see [17]). By Proposition B.2 we have that and by Proposition B.3 we have that \Delta A \p t (\Phi (q)) = 1 4 \Delta \bfq (2\pi 2 p t (q)).
Consequently, p t is the law of the process dq = H(q, t) dt + \\sigma which is exactly (109) thanks to (107).
Finally we are ready to prove Theorem 5.6.
Proof of Theorem 5.6. Using Proposition B.1, we deduce from (103) that which, thanks to (105), can be rewritten as The condition (106) in Proposition 5.9 is then satisfied with \H (A) = P T A (P D(M k )) and H(q) = P \bfq \bot F (q), so that, proceeding similarly as in Proposition 5.9 for the N -particles system, we conclude the result. (13)-- (14). In this section we show the equivalence between the macroscopic system (33)--(34) (or, equivalently, (12) for the last expression), expressed in terms of unitary quaternions, and the system (13)-- (14), expressed in term of rotation matrices from [17].

Comparison of the macroscopic model with
Recall \Phi , the natural map between unitary quaternions and rotation matrices defined in (90). We first notice that if q and \Lambda represent the same rotation (that is, if \Phi (q) = \Lambda ), then \Lambda e \bfone = e \bfone (\= q).
Therefore, the continuity equations (5) and (13) represent the same dynamics (it is direct from their definitions in (37) and in [17] that the constants c 1 and \c 1 are identical). We are left with comparing the various differential operators in \= q and \Lambda in (12) and (14).

Proof. Equation
where we have used successively the definition of \partial rel , the fact that in components w = (w 1 , w 2 , w 3 ), and (113). Recall that since \partial xi \= q is orthogonal to \= q, the product \partial rel,xi = (\partial xi \= q)\= q \ast is purely imaginary and can be identified with a vector in \BbbR 3 , so all the above terms make sense. Recall now the definition of D x (\Lambda ) in (16); we have just proved (114).
To check that (122)--(123) are equivalent, it suffices to show the correspondence between the constants since the equivalence of the terms is already given by (121). Therefore, we are left to check that \c 2 -\c 4 = c 2 and \c 4 = c 4 .

Conclusion.
In the present work we have introduced a flocking model for body attitude coordination where the body attitude is described through rotations represented by unitary quaternions. The deliberate choice of representing rotations by unitary quaternions is based on their numerical efficiency in terms of memory usage and operation complexity. This will be key for future applications of this model. At the modeling level, we introduce an individual-based model where agents try to coordinate their bodies' attitudes with those of their neighbors. To express this we needed to define an appropriate``averaged"" quaternion based on nematic alignment. This average is related to the Gennes Q-tensor that appears in liquid crystal theory. From the individual-based model we have derived the macroscopic equations (SOHQ) via the mean-field equations. We also show the equivalence between the SOHQ and the macroscopic equations (SOHB) of [17] where the body attitude is expressed through rotation matrices. However, we observe that the SOHQ is simpler to interpret than the equivalent SOHB. In particular, all the terms in the SOHQ are explicit. We have also seen that the dynamics of the SOHQ system are more complex than those of the SOH system (macroscopic equations corresponding to the Vicsek model). The body attitude coordination model presented here opens many questions and perspectives. We refer the reader to [17, Conclusions and open questions] for an exposition.
One may wonder why we did not consider translating the results in [17] for rotation matrices directly into quaternions. The answer is that, first, for the individualbased model, it is not possible to obtain a direct translation, in the sense that we need to consider some particular modeling choices (like the average in (21) and the relaxation in (23)) and check a posteriori the equivalence with the model in [17]. Second, the relation at the macroscopic level is not easy to obtain a priori. It is the macroscopic limit that gives us the necessary information and intuition to establish the link between both results.
In a future work, we will carry out simulations of the individual-based model and the SOHQ model, study the patterns that arise, and compare them with those of the Vicsek and SOH models. where \nabla \bfq is the gradient in \BbbH 1 .
Proof. To make the proof clearer we will use the notation \langle \cdot , \cdot \rangle rather than the symbol``\cdot "" to indicate the inner product (in the sense of matrices as well as in the sense of vectors and quaternions). We first check that (129) and (130) are equivalent: Indeed, since D \bfq \Phi (\nabla \bfq f (q)) belongs to T \Phi (\bfq ) = \{ [u] \times \Phi (q), u \in \BbbR 3 \} , (129) can be rewritten as, for all u \in \BbbR 3 , so that we recover (130).
We now prove (130): Fix some q \in \BbbH 1 , u \in \BbbR 3 and let \q = \q(s) \in \BbbH 1 be a differentiable path in This concludes the proof.
Proposition B.2 (comparison of the divergence operator). Let G be a vector field tangent to SO(3) and H a vector field tangent to \BbbH 1 such that (139) G(\Phi (q)) = D \bfq \Phi (H(q)) for all q \in \BbbH 1 .
Then we can compute for all f such that f (q) = f ( - q). This implies that (140) holds.