A new flocking model through body attitude coordination

We present a new model for multi-agent dynamics where each agent is described by its position and body attitude: agents travel at a constant speed in a given direction and their body can rotate around it adopting different configurations. In this manner, the body attitude is described by three orthonormal axes giving an element in $SO(3)$ (rotation matrix). Agents try to coordinate their body attitudes with the ones of their neighbours. In the present paper, we give the Individual Based Model (particle model) for this dynamics and derive its corresponding kinetic and macroscopic equations.


Introduction
In this paper we model collective motion where individuals or agents are described by their position and body attitude. The body attitude is given by three orthonormal axis; one of the axes describes the direction in which the agent moves at a constant speed; the other two axis indicate the relative position of the body with respect to this direction. Agents try to coordinate their body attitude with those of near neighbours (see Figure 1). Here we present an Individual Based Model (particle model) for body attitude coordination and derive the corresponding macroscopic equations from the associated mean-field equation, which we refer to as the Self-Organized Hydrodynamics for body attitude coordination (SOHB), by reference to the Self-Organized Hydrodynamics (SOH) derived from the Vicsek dynamics (see Ref. 24 and discussion below).
There exists already a variety of models for collective behaviour depending on the type of interaction between agents. However, to the best of our knowledge, this is the first model that takes into account interactions based on body attitude coordination. This has applications in the study of collective motion of animals such as birds and fish and it is a stepping stone to model more complex agents formed by articulated bodies (corpora). 13,14 In this section we present related results in the literature and the structure of the document.
There exists a vast literature on collective behaviour. In particular, here we deal with the case of self-propelled particles which is ubiquitous in nature. It includes, among others, fish schools, flocks of birds, herds 8,9,41 ; bacteria 5,47 ; human walking behaviour. 36 The interest in studying these systems is to gain understanding on the emergent properties: local interactions give rise to large scale structures in the form of patterns and phase transitions (see the review in Ref. 46). These large scale structures take place when the number of individuals or agents is very large. In this case 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 results presented here are inspired from those in Ref. 24. There the authors consider the Vicsek model which is a particular type of model for self-propelled particles. 1,15,34,45 The Vicsek model describes collective motion where agents travel at a constant speed in a given direction. At each time step the direction of movement is updated to the averaged one of the neighbouring agents, with some noise. The position is updated considering the distance travelled during that time step.
Our results here are inspired by the Self-Organized Hydrodynamics (SOH) model (the continuum version of the Vicsek model), where we have substituted velocity alignment by body attitude coordination. 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 the following ones as examples: in Ref. 10 an individual-based model is proposed to better describe collective motion of turning birds; in Ref. 25 agents are considered to have the shape of discs and volume exclusion is included in the dynamics; a description of nematic alignment in rods is done in Ref. 23. In Ref. 24 the authors investigate the mean-field limit and macroscopic limit of the Vicsek model. The mean-field limit gives a kinetic equation that takes the form of a Fokker-Planck equation referred to as the mean-field limit Vicsek model. (c) Dolphins moving in the same direction but with different body attitude. In this example one can see that the body attitude coordination model gives more information than the Vicsek model (which only describes the direction of movement). To obtain the macroscopic equations (the SOH model), the authors in Ref. 24 use the well-known tools of hydrodynamic limits, first developed in the framework of the Boltzmann equation for rarefied gases. 11,16,42 Since its first appearance, hydrodynamics limits have been used in other different contexts, including traffic flow modeling 3,35 and supply chain research. 2,26 However, in Ref. 24 a new concept is introduced: the Generalized Collision Invariant (GCI). Typically to obtain the macroscopic equations we require as many conserved quantities in the kinetic equation as the dimension of the equilibria (see again Ref. 46). In the mean-field limit 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 Refs. 17,18,21,22,27,29. After this introduction and the following discussion about the main result, the next part of the paper is dedicated to the modelling. In Section 3.1 we explain the derivation of the Individual Based Model for body coordination dynamics: given N agents labelled by k = 1, . . . , N the positions and body attitudes (X k , A k ) ∈ R 3 × SO(3) over time are given by the Stochastic Differential Equations (3.10)- (3.11). In Section 3.2 we give the corresponding (formal) mean-field limit (Prop. 3.2) for the evolution of the empirical measure when the number of agents N → ∞.
The last part concerns the derivation of the macroscopic equations (Theorem 4.1) for the total density of the particles ρ = ρ(t, x) and the matrix of the mean body attitude Λ = Λ(t, x). To obtain these equations we first study the rescaled mean-field equation (Eq. (4.1) in Section 4.1), which is, at leading order, a Fokker-Planck equation. We determine its equilibria, which are given by a von Mises distribution on SO(3) (Eq. (4.4), Section 4.3). Finally in Section 4.4 we obtain the Generalized Collision Invariants (Prop. 4.6), which are the main tool to be able to derive the macroscopic equations in Section 4.5.
2. Discussion of the main result: the Self-Organized Hydrodynamics for body attitude coordination (SOHB) The main result of this paper is Theorem 4.1 which gives the following macroscopic equations for the density of agents ρ = ρ(t, x) and the matrix of the mean body attitude Λ = Λ(t, x) ∈ SO(3) (i.e., the Self-Organized Hydrodynamics for body attitude coordination (SOHB)): In the equations above c 1 , c 2 , c 3 and c 4 are explicit constants which depend on the parameters of the model (namely the rate of coordination and the level of noise). The expressions of the constants c 2 , c 3 and c 4 depend on the Generalized Collision Invariant mentioned in the introduction (and computed thanks to Prop. 4.6). The constant c 1 is obtained as a "consistency" relation (Lemma 4.4). In (2.2), the operation [·] × transforms a vector v in an antisymmetric matrix such that [v] × u = v × u for any vector u (see (3.2) for the exact definition). The scalar δ x (Λ) and the vector r x (Λ) are first order differential operators intrinsic to the where ∇ x × is the curl operator. These operators are well-defined as long as Λ is smooth: as we will see in the next section, we can always express a rotation matrix as exp ([b] × ) for some vector b ∈ R 3 , and this function b → exp ([b] × ) is a local diffeomorphism between a neighborhood of 0 ∈ R 3 and the identity of SO(3). This gives a unique smooth representation of b in the neighborhood of 0 when x is in the neighborhood of x 0 since then Λ(x)Λ(x 0 ) −1 is in the neighborhood of Id.
We express Eq. (2.2) in terms of the basis vectors {Ω = Λe 1 , u = Λe 2 , v = Λe 3 } and we write Λ = Λ(Ω, u, v). System (2.1)-(2.2) can be expressed as an evolution system for ρ and the basis {Ω, u, v} as follows: The operator P Ω ⊥ denotes the projection on the orthogonal of Ω. We easily see here that these equations preserve the constraints |Ω| = |u| = |v| = 1 and Ω · u = Ω · v = u · v = 0. The expressions of δ and r are: Eq. (2.1) (or equivalently Eq. (2.3)) is the continuity equation for ρ and ensures mass conservation. The convection velocity is given by c 1 Λe 1 = c 1 Ω and Ω gives the direction of motion. Eq. (2.2) (or equivalently Eqs. (2.4)-(2.6)) gives the evolution of Λ. We remark that every term in Eq. (2.2) belongs to the tangent space at Λ in SO(3); this is true for the first term since (∂ t + c 2 (Λe 1 ) · ∇ x ) is a differential operator and it also holds for the second term because it is the product of an antisymmetric matrix with Λ (see Prop. Appendix A.2). Alternately, this means that (Ω(t), u(t), v(t)) is a direct orthonormal basis as soon as (Ω(0), u(0), v(0)).
The term corresponding to c 3 in (2.2) gives the influence of ∇ x ρ (pressure gradient) on the body attitude Λ. It has the effect of rotating the body around the vector directed by (Λe 1 ) × ∇ x ρ at an angular speed given by c3 ρ (Λe 1 ) × ∇ x ρ , so as to align Ω with −∇ x ρ. Indeed the solution of the differential equation dΛ dt +γ[n] × Λ = 0, when n is a constant unit vector and γ a constant scalar, is given by Λ(t) = exp(−γ t[n] × )Λ 0 , and exp(−γ t[n] × ) is the rotation of axis n and angle −γ t (see (3.4), it is called Rodrigues' formula). Since we will see that c 3 is positive the influence of this term consists of relaxing the direction of displacement Λe 1 towards ∇ x ρ. Alternately, we can see from (2.4) that Ω turns in the opposite direction to ∇ x ρ, showing that the ∇ x ρ term has the same effect as a pressure gradient in classical hydrodynamics. We note that the pressure gradient has also the effect of rotating the whole body frame (see influence of ∇ x ρ on u and v) just to keep the frame orthonormal. Similarly to what happens with the ∇ x ρ term in Eq. (2.2), the term containing c 4 ρ r in Eq. (2.4) has the effect of relaxing the direction of displacement Ω towards −r (we will indeed see that c 4 is positive). Finally, the last terms of Eqs. (2.5)-(2.6) have the effect of rotating the vectors u and v around Ω along the flow driven by D t at angular speed c 4 δ.
If we forget the term proportional to r in (2.4), System (2.3)-(2.4) is decoupled from (2.5)-(2.6), and is an autonomous system for ρ and Ω, which coincides with the Self-Organized Hydrodynamic (SOH) model. The SOH model provides the fluid description of a particle system obeying the Vicsek dynamics. 24 As already discussed in Ref. 24, the SOH model bears analogies with the compressible Euler equations, where (2.3) is obviously the mass conservation equation and (2.4) is akin to the momentum conservation equation, where momentum transport ρD t Ω is balanced by a pressure force −P Ω ⊥ ∇ x ρ. There are however major differences. The first one is the presence of the projection operation P Ω ⊥ which is there to preserve the constraint |Ω| = 1. Indeed, while the velocity in the Euler equations is an arbitrary vector, the quantity Ω in the SOH model is a velocity orientation and is normalized to 1. The second one is that the convection speed c 2 in the convection operator D t is a priori different from the mass convection speed c 1 appearing in the continuity equation. This difference is a signature of the lack of Galilean invariance of the system, which is a common feature of all dry active matter models.
The major novelty of the present model, which can be referred to as the Self-Organized Hydrodynamic model with Body coordination (or SOHB) is that the transport of the direction of motion Ω involves the influence of another quantity specific to the body orientation dynamics, namely the vector r. The overall dynamics tends to align the velocity orientation Ω, not opposite to the density gradient ∇ x ρ but opposite to a composite vector (c 3 ∇ x ρ + c 4 ρ r). The vector r is the rotational of a vector b locally attached to the frame (namely the unit vector of the local rotation axis multiplied by the local angle of rotation around this axis). This vector gives rise to an effective pressure force which adds up to the usual pressure gradient. It would be interesting to design specific solutions where this effective pressure force has a demonstrable effect on the velocity orientation dynamics.
In addition to this effective force, spatial inhomogeneities of the body attitude also have the effect of inducing a proper rotation of the frame about the direction of motion. This proper rotation is also driven by spatial inhomogeneities of the vector b introduced above, but are now proportional to its divergence.

Modelling: the Individual Based Model and its mean-field limit
The body attitude is given by a rotation matrix. Therefore, we work on the Riemannian manifold SO(3) (Special Orthogonal Group), which is formed by the subset of matrices A such that AA T = Id and det(A) = 1, where Id stands for the identity matrix.

ws-m3as
A new flocking model through body attitude coordination 7 In this document M indicates the set of square matrices of dimension 3; A is the set of antisymmetric matrices of dimension 3; S is the set of symmetric matrices of dimension 3. Typically we will denote by A, Λ matrices in SO(3) and by P matrices in A. Bold symbols n, v, e 1 indicate vectors.
We will often use the so-called Euler axis-angle parameters to represent an element in SO(3): to A ∈ SO(3) there is associated an angle θ ∈ [0, π] and a vector n ∈ S 2 so that A = A(θ, n) corresponds to the anticlockwise rotation of angle θ around the vector n. It is easy to see that tr(A) = 1 + 2 cos θ (3.1) (for instance expressing A in an orthonormal basis with n), so the angle θ is uniquely defined as arccos( 1 2 (tr(A)−1)). Notice that n is uniquely defined whenever θ ∈ (0, π) (if θ = 0 then n can be any vector in S 2 and if θ = π then the direction of n is uniquely defined but not its orientation). For a given vector u, we introduce the antisymmetric matrix [u] × , where [·] × is the linear operator from R 3 to A given by so that for any vectors u, v ∈ R 3 , we have [u] × v = u × v. In this framework, we have the following representation for A ∈ SO(3) (called Rodrigues' formula): We also have n × (n × v) = (n · v) n − (n · n)v, therefore when n is a unit vector, we have : [n] 2 × = n ⊗ n − Id, (3.5) where the tensor product a ⊗ b is the matrix defined by (a ⊗ b)u = (u · b)a for any u ∈ R 3 . Finally, SO(3) has a natural Riemannian metric (see Ref. 38) induced by the following inner product in the set of square matrices of dimension 3: This normalization gives us that for any vectors u, v ∈ R 3 , we have that Moreover, the geodesic distance on SO(3) between Id and a rotation of angle θ ∈ [0, π] is exactly given by θ (the geodesic between Id and A is exactly t ∈ [0, θ] → exp(t[n] × )). See Appendix Appendix A for some properties of SO(3) used throughout this work.
Seeing SO(3) as a Riemannian manifold, we will use the following notations: T A is the tangent space in SO(3) at A ∈ SO(3); P T A denotes the orthogonal projection onto T A ; the operators ∇ A , ∇ A · are the gradient and divergence in SO(3), respectively. These operators are computed in Section 4.2 in the Euler axis-angle coordinates.

The Individual Based Model
Consider N agents labelled by k = 1, . . . , N with positions X k (t) ∈ R 3 and associated matrices (body attitudes) A k (t) ∈ SO(3). For each k, the three unit vectors representing the frame correspond to the vectors of the matrix A k (t) when written as a matrix in the canonical basis (e 1 , e 2 , e 3 ) of R 3 . In particular, the direction of displacement of the agent is given by its first vector A k (t)e 1 .
Evolution of the positions. Agents move in the direction of the first axis with constant speed v 0 Evolution of the body attitude matrix. Agents try to coordinate their body attitude with those of their neighbours. So we are facing two different problems from a modelling viewpoint, namely to define the target body attitude, and to define the way agents relax their own attitude towards this "average" attitude.
As for the Vicsek model, 24 we consider a kernel of influence K = K(x) ≥ 0 and define the matrix (3.8) This matrix corresponds to the averaged body attitude of the agents inside the zone of influence corresponding to agent k. Now M k (t) / ∈ SO(3), so we need to orthogonalize and remove the dilations, in order to construct a target attitude in SO(3). We will see that the polar decomposition of M k (t) is a good choice in the sense that it minimizes a weighted sum of the squared distances to the attitudes of the neighbours. We also refer to Ref. 40 for some complements on averaging in SO(3).
We give next the definition of polar decomposition:  (i) The matrix A minimizes the quantity 1 Proof. We get the equivalence between the first two assertions by expanding: since A and A i (t) are both orthogonal matrices. So minimizing the weighted sum of the squares distances amounts to maximizing inner product of A and the weighted sum M k of the matrices A i given by (3.8). Therefore if det M k > 0, and A is the polar decomposition of M k , we immediately get that det A > 0, hence A ∈ SO(3). We know that S can be diagonalized in an orthogonal basis : S = P T DP with P T P = Id and D is a diagonal matrix with positive diagonal elements λ 1 , So the matrixB = P B T AP T maximizes tr(BD) = λ 1b11 + λ 2b22 + λ 3b33 among the elements of SO(3) (the map B → P B T AP T is a one-to-one correspondence between SO(3) and itself). But sinceB is an orthogonal matrix, all its column vectors are unit vectors, and so b ii 1, with equality for i = 1, 2 and 3 if and only ifB = Id, that is to say P B T AP T = Id, which is exactly B = A.
We denote by PD(M k ) ∈ O(3) the corresponding orthogonal matrix coming from the Polar Decomposition of M k .
We now have two choices for the evolution of A k . We can use the second point of Proposition 3.1 and follow the gradient of the function to maximize : (see (A.2) for the last computation, P T A k is the projection on the tangent space, this way the solution of the equation stays in SO (3)).
Or we can directly relax to the polar decomposition PD(M k ), in the same manner: We can actually see that the trajectory of this last equation, when PD(M k ) belongs to SO(3) and does not depend on t, is exactly following a geodesic (see Prop. Appendix A.4). Therefore in this paper we will focus on this type of coordination. The positive coefficient ν gives the intensity of coordination, in the following we will assume that it is a function of the distance between A k and PD(M k ) (the angle of the rotation A T k PD(M k )), which is equivalent to say that ν depends on A k · PD(M k ). We assume that the motion of the agents is smooth enough so that the average M k stays 'close' to SO(3) and that, in particular, det(M k ) > 0. A simple example is when we only average two different matrices 3) and to (3.5). Since the matrix S = cos θ 2 Id + (1 − cos θ 2 )n ⊗ n is a positive-definite symmetric matrix as soon as θ ∈ [0, π), we have that det(M) > 0. The polar decomposition of M is then A, which is the midpoint of the geodesic joining A 1 to A 2 (which corresponds to the curve t ∈ [0, θ] → A 1 exp(t[n] × )). As soon as we average more than two matrices, there exist cases for which det(M) < 0: for instance if we take Noise term. Agents make errors when trying to coordinate their body attitude with that of their neighbours. This is represented in the equation of A k by a noise term: From all these considerations, we obtain the IBM where the Stochastic Differential Equation is in Stratonovich sense (see Ref. 32). The projection P T A k and the fact that we consider the SDE in Stratonovich sense ensures that the solution A k (t) stays in SO (3). The normalization constant 2 √ D ensures that the diffusion coefficient is exactly D : the law p of the underlying process given by which is encountered when considering diffusion process on manifolds isometrically embedded in the euclidean space R n , because we are here considering SO(3) embedded in M (isomorphic to R 9 ), but with the metric (3.6), which corresponds to the canonical metric of R 9 divided by a factor 2. We refer to the book 37 for more insight on such stochastic processes on manifolds.

Mean-field limit
We assume that the kernel of influence K is Lipschitz, bounded, with the following properties: In Ref. 7 the mean-field limit is proven for the Vicsek model. Using the techniques there it is straightforward to see that for of the empirical measure associated to the Stratonovich Stochastic Differential Equation (SDE): converges weakly f N → f as N → ∞. The limit satisfies the kinetic equation: The equations we are dealing with (3.10)-(3.11), since we consider the Polar Decomposition of the averaged body attitude M k , are slightly different from (3.13)-(3.14), which would correspond to the modelling point of view of Eq. (3.9). As a consequence, the corresponding coefficient of the SDE is not Lipschitz any more and the known results for existence of solutions and mean-field limit (see Th. 1.4 in Ref. 43) fail. More precisely, the problem arises when dealing with matrices with determinant zero; the orthogonal matrix of the Polar Decomposition is not uniquely defined for matrices with determinant zero and, otherwise, . A complete proof of the previous results in the case of Eq. (3.10)-(3.11) would involve proving that solutions to the equations stay away from the singular case det(M k ) = 0. This is an assumption that we make on the Individual Based Model (see the second point of Remark 3.1). This kind of analysis has been done for the Vicsek model (explained in the introduction) in Ref. 28 where the authors prove global well-posedness for the kinetic equation in the spatially homogeneous case.
In our case one expects the following to hold: . When the number of agents in (3.10)-(3.11) N → ∞, its corresponding empirical distribution converges ,

Hydrodynamic limit
The goal of this section will be to derive the macroscopic equations (Theorem 4.1). From now on, we consider the kinetic equation given in (3.15).

Scaling and expansion
We express the kinetic Eq. (3.15) in dimensionless variables. Let ν 0 be the typical interaction frequency scale so that ν( We introduce also the typical time and space scales t 0 , x 0 such that t 0 = ν −1 0 and x 0 = v 0 t 0 ; the associated variables will be t = t/t 0 and x = x/x 0 . Consider the dimensionless diffusion coefficient d = D/ν 0 and the rescaled influence kernel K (|x |) = K(x 0 |x |). Skipping the primes we get , Here d, ν and K are assumed to be of order 1.
Remark 4.1. Notice in particular that before and after scaling the ratio Now, to carry out the macroscopic limit we rescale the space and time variables by settingt = εt,x = εx to obtain (skipping the tildes): Assuming that f is sufficiently smooth (with bounded derivatives), we have the expansionM where Proof. This is obtained by performing the change of variable x = x + εξ in the definition of M ε [f ] and using a Taylor expansion of f (x+εξ, A , t) with respect to ε. We use that K is isotropic and with bounded second moment by assumption (see Eq. (3.12)).
From the lemma, we rewrite , .
, Q(f ) and F 0 [f ] are non-linear operators of f , which only acts on the attitude variable A.

Preliminaries: differential calculus in SO(3)
In the sequel we will use the volume form, the gradient and divergence in SO(3) expressed in the Euler axis-angle coordinates (θ, n) (explained at the beginning of Section 3). In this section we give their explicit forms; the proofs are in appendix Appendix A.  (3)). Let f : is the expression of f in the Euler axisangle coordinates by Rodrigues' formula (3.3), we have then where A = A(θ, n) and ∇ n is the gradient on the sphere S 2 .
The volume form in SO (3) is left invariant (it is the Haar measure), due to the fact that the inner product in M is also left invariant: (3). We give its expression in the Euler axis-angle coordinates (θ, n) : where dn is the Lebesgue measure on the sphere S 2 , normalized to be a probability measure, and We have seen in Prop. 4.1 that the gradient is decomposed in the basis (3)), and suppose that for some smooth function b and smooth vector function v such that v(θ, n) ⊥ n. Then Now we can compute the Laplacian in SO (3): (3) can be expressed as where ∆ n is the Laplacian on the sphere from here we just need to apply Prop. 4.2 knowing that (n × ∇ nf ) × n = ∇ nf since ∇ nf is orthogonal to n.

Equilibrium solutions and Fokker-Planck formulation
We define a generalization of the von-Mises distributions on SO (3) by and we also obtain that M Λ (A) is actually M Id (Λ T A).
We are now ready to describe the properties of Q in terms of these generalized von-Mises distributions. i) The operator Q can be written as and we have ii) The equilibria, i.e., the functions f = f (x, A, t) such that Q(f ) = 0 form a 4dimensional manifold E given by where ρ is the total mass while Λ is mean body attitude of ρM Λ (A), i.e., Furthermore, H(f ) = 0 iff f = ρM Λ for arbitrary ρ ∈ R + and Λ ∈ SO (3).
To prove Lemma 4.3 we require the following one, which is of independent interest and for which we introduce the following notation: for any scalar function g : (0, π) → R and a given integrable scalar function h : (0, π) → R which remains positive (or negative) on (0, π), we define for m(θ) = exp(d −1 σ( 1 2 + cos θ)). (4.8) Proof. Using the fact that the measure on SO (3) is left invariant, we obtain We now write B = Id + sin θ [n] × + (1 − cos θ)[n] 2 × thanks to Rodrigues' formula (3.3). Therefore, using Lemma 4.2, we get Next, we see that since the function n → [n] × is odd, we have S 2 [n] × dn = 0. We also have (see (3.5)) that [n] 2 × = n ⊗ n − Id. Since we know that S 2 n ⊗ ndn = 1 3 Id (by invariance by rotation), it is easy to see that the integral in S 2 has to be proportional to Id, the coefficient is given by computing the trace), we get that which gives the formula (4.7) for c 1 .
One can actually check that the probability measure M Λ on SO(3) converges in distribution to the uniform measure when d → ∞ (by Taylor expansion) and it converges to a Dirac delta at matrix Λ when d → 0 (this can be seen for M Id thanks to the decomposition of the volume form and the Laplace method, since the maximum of σ( 1 2 + cos θ) is reached only at θ = 0 which corresponds to the identity matrix, and we then get the result for any Λ since M Λ (A) = M Id (Λ T A)). So for small diffusion, at equilibrium, agents tend to adopt the same body attitude close to Λ.
With these asymptotic considerations, we have in particular the behaviour of c 1 :

Generalized Collision Invariants
To obtain the macroscopic equation, we start by looking for the conserved quantities of the kinetic equation: we want to find the functions ψ = ψ(A) such that By Lemma 4.3, this can be rewritten as This happens if ∇ A ψ ∈ T ⊥ A which holds true only if ∇ A ψ = 0, implying that ψ is constant.
Consequently, our model has only one conserved quantity: the total mass. However the equilibria is 4-dimensional (by Lemma 4.3). To obtain the macroscopic equations for Λ, a priori we would need 3 more conserved quantities. This problem is sorted out by using Generalized Collision Invariants (GCI) a concept first introduced in Ref. 24.

Definition and existence of GCI
Define the operator Using this operator we define: In particular, the result that we will use is : which ends the proof since the expression of the adjoint is L * (ψ) = 1 We prove the existence and uniqueness of the solution ψ satisfying Eq. (4.10) in the following: Our goal is to prove that there exists a unique ψ ∈ H 1 (SO(3)) such that a(ψ, ϕ) = b(ϕ) for all ϕ ∈ H 1 (SO (3)).
To begin with we apply the Lax-Milgram theorem on the space In this space the H 1 -norm and the H 1 semi-norm are equivalent thanks to the Poincaré inequality, i.e., there exists C > 0 such that |ϕ| 2 dA for some C > 0, for all ϕ ∈ H 1 0 (SO (3)).
Notice that since the mapping P → ψ Λ0 P is linear and injective from A (of dimension 3) to H 1 0 (SO(3)), the vector space GCI(Λ 0 ) is of dimension 4.

The non-constant GCIs
From now on, we omit the subscript on Λ 0 , and we are interested in a simpler expression for ψ Λ P . Rewriting expression (4.10) using (4.12), for any given P ∈ A we want to find ψ such that (4.13) Proposition 4.5. Let P ∈ A and ψ be the solution of (4.13) belonging to H 1 0 (SO (3)). If we denoteψ(B) := ψ(ΛB), thenψ is the unique solution in H 1 0 (SO(3)) of: (4.14)

Proof. Let ψ(A) =ψ(Λ T A). Consider A(ε) a differentiable curve in SO(3) with
Then, by definition and therefore we have that We conclude that implying (since this is true for any δ A ∈ T A ) that Now to deal with the divergence term, we consider the variational formulation. Consider ϕ ∈ H 1 (SO(3)), then our equation is equivalent to (3)). The left hand side can be written as: and the right hand side is equal to where we define analogouslyφ(B) = ϕ(ΛB). This concludes the proof.
Going back to the GCI ψ(A), we can write it as  , n), the symmetric part of B(θ, n) gives no contribution when computing P · B and we get ψ(B) = P · Bψ 0 ( 1 2 tr(B)) = sin θ ψ 0 (θ)P · [n] × = sin θ ψ 0 (θ)(p · n), where the vector p is such that P = [p] × and this leads to Using that the Laplacian in the sphere has the property ∆ n (p · n) = −2(p · n) (p · n corresponds to the first spherical harmonic), we conclude that expression (4.16) is satisfied. In the computation we used the same procedure as for the proof of the expression of the Laplacian in SO(3) (Corollary 4.1), but (using the same notations) we have taken b(θ, n) = m(θ)∂ θ (sin θ ψ 0 (θ))(p · n).
To conclude the proof we just need to check that ψ 0 exists and corresponds to a functionψ in H 1 0 (SO(3)). Using the expression of the volume form, since S 2 p · n dn = 0, we get that if ψ 0 is smooth, we have SO(3)ψ (A)dA = 0, and using the expression of the gradient, we get that This Hilbert space is equipped with the corresponding norm: |∂ θ (sin θψ(θ))| 2 sin 2 (θ/2) dθ.

The macroscopic limit
In this section we investigate the hydrodynamic limit. To state the theorem we first give the definitions of the first order operators δ x and r x . For a smooth function Λ from R 3 to SO(3), and for x ∈ R 3 , we define the following matrix D x (Λ) such that for any w ∈ R 3 , we have Notice that this first-order differential equation D x is well-defined as a matrix; for a given vector w, the matrix (w · ∇ x )Λ is in T Λ and thanks to Prop. Appendix A.3, it is of the form P Λ, with P an antisymmetric matrix. Therefore there ex- We now define the first order operators δ x (scalar) and r x (vector), by We first give an invariance property which allows for a simple expression for these operators. Consequently, in the neighborhood of x 0 ∈ R 3 , we can write and therefore Proof. For any w ∈ R 3 , we have, since A is constant: This proves that D x (ΛA) = D x (Λ), and by (4.20), the same is obviously true for δ x and r x . We now write, in the neighborhood of We then get, since exp([b(x 0 )]) = Id, that We are now ready to state the main theorem of our paper (see Section 2 for a discussion on this result).
Proof. Suppose that f ε → f as ε → 0, then using (4.1) we get Q(f ε ) = O(ε), which formally yields Q(f ) = 0 and by Lemma 4.3 we have that Using the conservation of mass (integrating (4.1) on SO(3)), we have that Ae 1 f ε dA, and in the limit (formally) thanks to Lemma 4.4. This gives us the continuity equation (4.21) for ρ. Now, we want to obtain the equation for Λ. We write Λ ε = Λ[f ε ], and we take P ∈ A a given antisymmetric matrix. We consider the non-constant GCI associated to Λ ε and corresponding to P in (4.18): ψ ε (A) = P · ((Λ ε ) T A)ψ 0 (Λ ε · A). Since we have ψ ε ∈ GCI(Λ[f ε ]), we obtain, thanks to the main property (4.9) of the GCI, that Multiplying (4.1) by ψ ε , integrating w.r.t. A on SO(3) and using the expression of ψ ε as stated above, we obtain Assuming the convergence f ε → f is sufficiently strong, we get in the limit (4.24) Since (4.24) is true for any P ∈ A, the matrix is orthogonal to all antisymmetric matrices. Therefore, it must be a symmetric matrix, meaning that we have We have with the definition of M Λ in (4.4) that Inserting the two previous expressions into (4.25), we compute separately each component of X defined by: For the first term we have (changing variables B = Λ T A): For the term X 2 we make the change of variables B = Λ T A and compute where we have used the expression of the Haar measure dB = 2 π sin 2 (θ/2)dθdn (see Lemma 4.2) and that writing B = B(θ, n) = Id + sin θ[n] × + (1 − cos θ)[n] 2 × thanks to Rodrigues' formula (3.3), we have B − B T = 2 sin θ[n] × . Removing odd terms with respect to the change n → −n, we obtain for some vector λ t λ t λ t . Therefore So using the definition (4.23) of m(θ), we get because the mapping w → [w] × is linear, and S 2 n ⊗ n dn = 1 3 Id. Denote by then we conclude that Now, for the term X 3 we compute the following, starting again by the change of variables B = Λ T A: where we used similar considerations as for X 2 , as well as that

Denote by
We now compute X 4 in the same way, with the change of variables B = Λ T A: We now use the definition of D x (Λ) given in (4.19) to get Using the fact that To simplify the notations, we denote L = Λ T D x (Λ)Λ. Since the symmetric part of B does not contribute to the scalar product B · [LBe 1 ] × , we get Therefore we obtain, in the same manner as before, n · L (cos θe 1 + (1 − cos θ)(n · e 1 ) n) n dn where the term involving [n] × vanishes since its integrand is odd with respect to n → −n.
To compute the second term of this expression we will make use of the following lemma proved at the end of this section: π 0 m(θ) sin 2 (θ/2)(1 + 4 cos θ) dθ, Finally putting all the terms together we have that In particular ΛX = 0 and from the fact that Λ[w] × = [Λw] × Λ for any w ∈ R 3 we get Since we have taken L = Λ T D x (Λ)Λ, we get that tr(L) = tr D x (Λ) = δ x (Λ) and, thanks to (4.20): thanks to the definition of D x given in (4.19). Finally, inserting these expressions into (4.26) and dividing by C 2 , we get the equation which ends the proof.
from which we conclude the lemma. In the computations we used that Finally, we consider the orthonormal basis given by where {e 1 , e 2 , e 3 } is the canonical basis of R 3 . We can have an expression of the operators δ x and r x in terms of these unit vectors {Ω, u, v}, which allows to rewrite the evolution equation of Λ as three evolution equations for these vectors.
Proposition 4.8. We have Consequently, we have the following evolution equations for Ω, u, and v, corresponding to the evolution equation of Λ given in (4.22): where D t := ∂ t + c 2 (Ω · ∇ x ), and where δ x (Ω, u, v) is the expression of δ x (Λ) given by (4.27).
Proof. We first prove (4.27). We have thanks to the definition of D x given in (4.19). Now we use the fact that for two matrices A, B, we have A · B = 1 2 tr(A T B) = 1 2 i Ae i · Be i (half the sum of the scalar products of the corresponding columns of the matrices A and B), to get For this last equality we used the fact that since u ⊥ v and analogously for the other components.
We proceed next to proving the expression of r x (Λ) given by (4.28). We first prove that r x (Λ) · Ω = ∇ x · Ω. We have (recall that [r and that for all w in Since (Ω · ∇ x )Ω is orthogonal to Ω, we therefore get since ΛΛ T = Id (the first line is actually the expression of the divergence of Ω in the basis {Ω, u, v}). For the other two components of r x (Λ), we perform exactly the same computations with a circular permutation of the roles of Ω, u, v to get r Therefore we obtain (4.28). Finally we rewrite the equation for Λ as the evolution of the basis {Ω, u, v}. To obtain the evolution of Λe k for k = 1, 2, 3, we multiply the Eq. (4.22) by e k and compute to obtain: where D t = ∂ t +c 2 (Ω·∇ x ). To perform the computations we have used, for w = ∇ x ρ or w = r that [w × Ω] × Ω = −P Ω ⊥ (w) and (w × Ω) × u = (u · w)Ω since Ω ⊥ u (analogously for v). From here, using (4.28) we obtain straightforwardly Eqs. (4.29), (4.30), and (4.31) for Ω, u and v respectively.

Conclusions and open questions
In the present work we have presented a new flocking model through body attitude coordination. We have proposed an Individual Based Model where agents are described by their position and a rotation matrix (corresponding to the body attitude). From the Individual Based Model we have derived the macroscopic equations via the mean-field equations. We observe that the macroscopic equation gives rise to a new class of models, the Self-Organized Hydrodynamics for body attitude coordination (SOHB). This model does not reduce to the more classical Self-Organized Hydrodynamics (SOH), which is the continuum version of the Vicsek model. The dynamics of the SOHB system are more complex than those of the SOH ones of the Vicsek model. In a future work, we will carry out simulations of the Individual Based Model and the SOHB model and study the patterns that arise to compare them with the ones of the Vicsek and SOH model.
Also, there exist yet many open questions on the modelling side. For instance, one could consider that agents have a limited angle of vision, thus the so-called influence kernel K (see Section 3.1) is not isotropic any more, see Ref. 29 for the case of the Vicsek and SOH models. One could also consider a different interaction range for the influence kernel K that may give rise to a diffusive term in the macroscopic equations, see Ref. 19. Moreover, in the case of the SOH model, when the coordination frequency and noise intensity (quantities ν and D in the Individual Based Model (3.10)-(3.11)) are functions of the flux of the agents, then phase transitions occur at the macroscopic level, 19 (see also Refs. 4,6,20,44). An analogous feature is expected to happen in the present case. Finally, one could think of elaborating on the model by adding repulsive effects at short range and attraction effects at large range.
On the analytical side, this model opens also many questions like making Prop. 3.2 rigorous, which means dealing with Stochastic Differential Equations with non-Lipschitz coefficients. In the context of the Vicsek model, the global wellposedness has been proven for the homogeneous mean-field Vicsek equation and also its convergence to the von Mises equilibria in Ref. 28, see also Ref. 31; an analogous result for our model will be desirable. The convergence of the Vicsek model to the model which was formally done in Ref. 24 has been recently achieved rigorously in Ref. 39. Again, one could also think of generalizing these results to our case.

Appendix A. Special Orthogonal Group SO(3)
Throughout the text, we used repeatedly the following properties: So if M ∈ T A , we must have (A T M +M T A) = 0, that is to say that P = A T M ∈ A.
Conversely if M = AP with P ∈ A, the solution of the linear differential equation Λ (t) = Λ(t)P with Λ(0) = A is given by Λ(t) = Ae tP so it is a curve in SO (3). Notice that then Proof. It suffices to verify that the expression given for To compute the gradient in SO(3) of a function ψ : SO(3) → R we will consider A(ε) a differentiable curve in SO(3) such that In particular, one can check that We now show that the differential equation corresponding to following this gradient has trajectories supported on geodesics. thanks to (A.1) and we have A(0) = A 0 . Since θ ∈ [0, θ 0 ] → exp(θ[n] × ) is a geodesic between Id and B T A 0 , then θ → B exp(θ[n] × ) is a geodesic between B and A, and the solution A(t) is supported on this geodesic. It is also easy to see that, except in the case θ 0 = π or θ 0 = 0, for which the solution is constant, the function t → θ(t) (solution of the one-dimensional differential equation θ = −ν( 1 2 + cos θ) sin θ) is positive, decreasing, and converge exponentially fast to 0, with an asymptotic exponential rate ν( 3 2 ). Therefore, as time goes to infinity, the trajectory covers the whole geodesic from A 0 to B (excluded).

Proposition Appendix
We now turn to the proofs of the expressions of the gradient, the volume form and the divergence in SO(3) in the so-called Euler axis-angle coordinates, that were presented in section 4.2. With these notations, for a function f = f (A(θ, n)) it holds: On the other hand, it holds true that where L [n]× is given in (A.4).