Phase transitions and macroscopic limits in a BGK model of body-attitude coordination

In this article we investigate the phase transition phenomena that occur in a model of self-organisation through body-attitude coordination. Here, the body-attitude of an agent is modelled by a rotation matrix in $\mathbb{R}^3$ as in [Degond, Frouvelle, Merino-Aceituno, 2017]. The starting point of this study is a BGK equation modelling the evolution of the distribution function of the system at a kinetic level. The main novelty of this work is to show that in the spatially homogeneous case, self-organisation may appear or not depending on the local density of agents involved. We first exhibit a connection between body-orientation models and models of nematic alignment of polymers in higher dimensional space from which we deduce the complete description of the possible equilibria Then, thanks to a gradient-flow structure specific to this BGK model, we are able to prove the stability and the convergence towards the equilibria in the different regimes. We then derive the macroscopic models associated to the stable equilibria in the spirit of [Degond, Frouvelle, Merino-Aceituno, 2017] and [Degond, Frouvelle, Liu, 2015].


Introduction
The model studied in the present work is a new elaboration of the work initiated in [14] to model collective behaviour of agents described by their position and body-attitude. New results about emergence of phenomena of body-attitude coordination are presented in the context of a Bhatnagar-Gross-Krook (BGK) model. Such models can be applied to many biological systems such as flocking birds [37], fish school [35,36] or sperm motion [20]. These systems are constituted by a large number of self-propelled agents which move at a constant speed and try to imitate their neighbours by moving in the same direction and trying to coordinate their body attitude. The agents are modelled by a moving frame in dimension 3, i.e. three orthogonal axes, one of which gives the direction of the motion and the two others the body orientation. In this work, as in [14], the body attitude is modelled by a rotation matrix in dimension 3, i.e. an element of the special orthogonal group SO 3 (R). In [17] agents are modelled by quaternions.
Collective behaviour in many-agent systems has been a thoroughly studied subject in the mathematical literature, from both theoretical and applied points of view. Among the models which have received the most attention, one can cite the Cucker-Smale model [10,31,47], attractive-repulsive models [7] or the Vicsek model for self-propelled particles [54]. The present work belongs to the class of Vicsek-inspired models. Such models have two main 1 Introduction distinctive features, first the assumption that the particles are self-propelled and secondly a geometrical constraint: in the original work of Vicsek, the velocities of the particles have constant norm and the dynamics therefore takes place on the sphere S n−1 in dimension n (n = 2 in [54]). Here the dynamics takes place on the Riemannian manifold SO 3 (R).
The tools used to study models of collective behaviour are generally borrowed from the mathematical kinetic theory of gases which gives a mathematical framework to study manyparticle systems. At a microscopic scale, the motion of each particle is detailed (Individual Based Model, IBM) through Ordinary Differential Equations (ODE) coming from Newton's laws or through stochastic processes. When the number of particles is large, the whole system is described at a mesoscopic scale by a kinetic partial differential equation such as the Boltzmann, Fokker-Planck or BGK equation. Finally, large-scale dynamics is described by macroscopic equations (Euler, Navier-Stokes. . . ). A review of the main results of kinetic theory of gases can be found in [11]. In particular the BGK equation (for Bhatnagar-Gross-Krook) was introduced in [3] as a substitute for the Boltzmann equation in the context of gas dynamics. The BGK operator is a relaxation operator towards a Maxwellian having the same moments as the distribution function of the system. Its mathematical properties and relevance in the mathematical kinetic theory of gases have been studied in particular in [49] and [51]. The BGK operator has been used in a model of collective dynamics of self-propelled particles in [22]. However, together with [16], it is the first time that it is rigorously studied in a body-attitude coordination model.
The main mathematical challenge in classical kinetic theory is the rigorous derivation of the kinetic equations from the IBM and of the macroscopic models from the kinetic equations. These questions are at the core of Hilbert's sixth problem and have received much attention in the last decades. Many different techniques have been developed to derive kinetic equations from hard-sphere gases (Boltzmann-Grad limit [25,44]), from systems of interacting particles (mean-field limit and propagation of chaos [34,40,53]) or from stochastic processes (and in particular jump processes [43,46]). Some of these techniques have been adapted to problems arising in the study of collective behaviour [4,5,9]. The passage from kinetic equations to macroscopic models generally depends on physical constraints and in particular on conservation laws (hydrodynamic limits, Hilbert and Chapman-Enskog methods, see [8,11] for a review) and is still an active research field [6,23,29,30]. In the context of self-propelled particles, due to the lack of conservation laws which normally hold in the classical kinetic theory of gases, specific tools are needed. In [18], a methodological breakthrough has been achieved by introducing the so called Generalised Collisional Invariants (GCI) to rigorously link kinetic and macroscopic equations in the context of collective behaviour of self-propelled particles. This technique is now rigorously justified [42] and has already been successfully applied to a wide range of problems [41,56]. It will be the key here to derive the macroscopic model in Section 6. This will lead to a system of partial differential equations on the mean density and body attitude, referred as the Self-Organised Hydrodynamics for Body-attitude coordination (SOHB) in [14].
The aim of this work is to show the emergence of collective behaviour and self-organisation which give rise to macroscopic scale patterns such as clusters, travelling bands etc. These patterns emerge from the collective interactions and are not directly encoded in the behaviour of the individual particles as described by the IBM. The continuum version of the Vicsek model [18] named the Self-Organised Hydrodynamics (SOH) model is an exemple of a model able to describe such emergence of self-organised dynamics. The Vicsek model describes a system where agents try to imitate their neighbours by adapting their direction of motion to the average direction of their neighbours. It has been shown that, in a certain scaling and when the equilibrium of the system is reached, the directions of motion of the agents are not uniformly distributed but follow a von Mises distribution. For κ ∈ R + and Ω ∈ S n−1 the 1 Introduction von Mises distribution of parameters κ and Ω is the Probability Density Function (PDF) on S n−1 defined by: where the dot product is the usual dot product in R n . This model [18] has been the starting point of many other models of self-organised dynamics, including [14] for the body-attitude coordination. In this context, we define the von Mises distribution of parameter J ∈ M 3 (R) (a 3 × 3 real matrix) as the following PDF on SO 3 (R): where the dot product and the measure on SO 3 (R) come from the Riemannian structure of SO 3 (R) detailed in Section 3.
In the present work, we focus on phase transition phenomena between non-organised and organised dynamics (collective motion). We will prove that the spatial density of agents is the key parameter which encodes the main features of phase transitions: in low density regions, no self-organised dynamics appears but when the density crosses a critical value, self-organised dynamics, given by a von Mises distribution for the body-attitude, becomes a stable equilibria of the system. This phase transition in the dynamics is purely an emergent phenomena, in the sense that at the macroscopic scale, different equations are required to describe the dynamics for different values of the density of agents, whereas for the IBM and at a mesoscopic level, the dynamics is described by one unique (system of) equation(s).
The starting point of this study is the BGK equation where f (t, x, A) is a probability measure which gives the distribution of agents at position x ∈ R 3 with body-orientation A ∈ SO 3 (R) at time t ∈ R + and where: The left-hand side of the equation models the transport phenomenon: an agent with body orientation A ∈ SO 3 (R) moves in the direction Ae 1 where e 1 is the first vector of the canonical basis of R 3 . The right-hand side of the equation is the BGK operator which models the interactions between the agents: here we assume that f relaxes towards a "moving equilibrium" which takes the form of a von Mises distribution. In particular, the von Mises distribution appears as the analog of the Maxwellian distribution of the classical gas dynamics. The flux J f plays the same role as the momentum density for gas dynamics or the average flux for the Vicsek model. The term ρ f M J f can therefore be seen as the analog of the "Maxwellian distribution with same moments as f " in the context of the BGK equation for gas dynamics.
The main results of this work are (informally) summarised in the two following theorems. Theorem 1. Let us consider the spatially homogeneous BGK equation: where ρ ∈ R + is a given density of agents.
1. The equilibria f eq of the spatially homogeneous BGK model are either the uniform equilibrium f eq = ρ or of the form f eq = ρM αΛ or f eq = ρM α p⊗q where Λ ∈ SO 3 (R) and p, q ∈ S 2 and where α ∈ R and ρ are linked by a compatibility equation to be defined later (see Section 4 and equations (20) and (21)).
2. Depending on the density of agents ρ ∈ R + , the only stable equilibria are either the uniform equilibrium f eq = ρ or the equilibria of the form f eq = ρM αΛ where Λ ∈ SO 3 (R) and where α ∈ R + is linked to ρ by a compatibility equation to be defined later.
The first point of this theorem is detailed in Section 4 (see in particular Theorem 5 and Corollary 4.2). The second point is detailed in Section 5 (see in particular Theorem 7). We will then prove the following result.
Theorem 2 (Formal). Let us consider the rescaled spatially inhomogeneous problem 1. We assume that in a disordered region, f ε converges as ε → 0 towards a density ρ = ρ(t, x) uniform in the body-attitude variable. Then the density ρ ε ≡ ρ f ε satisfies at first order the following diffusion equation: 2. We assume that in an ordered region, f ε converges as ε → 0 towards an equilibrium of the form ρM αΛ with ρ ∈ R + , α ∈ R + and Λ ∈ SO 3 (R) defined above. Then the density ρ = ρ(t, x) and mean body attitude Λ = Λ(t, x) ∈ SO 3 (R) satisfy the SOHB model given by the following system of partial differential equations: where α = ρc 1 (α) andc 2 ,c 3 , c 4 are functions of ρ to be defined later and δ and r are the "divergence" and "rotational" operators defined in [14] (see Section 6) This theorem is detailed in Section 6 (see in particular Proposition 6.1 and Theorem 9).
The phase transition problem has been completely treated in the space-homogeneous case for the Vicsek model in [13] but the geometrical structure inherent to body-orientation models requires specific tools and techniques. In particular, the rotation group SO 3 (R) is a compact Lie group, endowed with a Haar measure. The links between this topological structure and the Riemannian structure (detailed in Section 3 and Appendix B) will be the key to reduce the problem to a form that shares structural properties with the models of nematic alignment of polymers, studied in a completely different context to model liquid crystals [1,2,32,55,57]. These two worlds will be formally linked through the isomorphism between SO 3 (R) and the group of unit quaternions detailed in Section 4.2 and Appendix A. It will lead to the first point of Theorem 1 (the complete description of the equilibria, Section 4). As in [55] we will see that there exist a class of equilibria which cannot be interpreted as equilibria around a mean-body orientation. These equilibria were not studied in [14,17]. A key point of the proof will be the reduction to a problem for diagonal matrices which will be a consequence of the left and right invariance of the Haar measure together with an adapted version of the Singular Value Decomposition of a matrix (Definition 3.2).
The stability of the different equilibria, are studied in Section 5.2. We will show that our model has an underlying gradient-flow structure which will allow us to determine the asymptotic behaviour of the system after a reduction to an ODE in R 3 . This is a specificity of the BGK model which doesn't hold for the other models of body-attitude coordination [14,17] and allows us to use different and simpler techniques. In particular, we will prove that the equilibria which cannot be interpreted as equilibria around a mean body-orientation are always unstable, which tends to justify the analysis carried out in [14] for a model where only equilibria around a mean body-orientation were considered.
Finally, the SOHB model (Section 6) will be obtained as in [14] by using the GCI. However, compared to [14], additional terms appear which require a specific treatment and in particular the coefficientc 3 that appears in Theorem 2 is different from the one that appears in [14]. The SOHB model (42) raises many questions, most of which are still open, and its mathematical and numerical analyses are still in progress. In particular, the hyperbolicity of the model is currently under study [15] and has been shown whenc 3 is constant.
The organisation of the work is the following: in Section 2 we will give a review of the existing models at a microscopic and mesoscopic scales and motivate the study of the BGK equation among them. In Section 3, we gather the main technical results we will constantly use throughout this work. In Section 4, we will describe, depending on the density, all the possible equilibria of the system. We will use the tools developed to mathematically study the alignment of polymers [55,57]. In Section 5 we will describe the asymptotic behaviour of the system and in particular which equilibria are attained, leading to a self-organised dynamics or not. This will be based on a specific underlying gradient-flow structure of the BGK equation. Finally in Section 6 we will write the macroscopic models for the stable equilibria.
Notations. For the convenience of the reader, we collect here the main notations we will use in the following.
• M n (R) is the set of n × n real matrices.
• D n (R) ⊂ M n (R) is the subspace of n × n diagonal real matrices.
• Tr(M ) denotes the trace of the matrix M ∈ M n (R) and M T its transpose.
• I n denotes the identity matrix in dimension n.
• diag : R n → D n (R) is the vector space isomorphism such that for (d 1 , . . . , d n ) ∈ R n , • S n (R) and A n (R) denote respectively the sets of symmetric and skew-symmetric matrices of dimension n.
• SO n (R) is the special orthogonal group in dimension n, i.e. the group of matrices P ∈ M n (R) such that P P T = I n and det P > 0.
• S n ⊂ R n+1 is the sphere of dimension n.
• H is the group of unitary quaternions.
• · g denotes the mean for the probability density g on SO n (R). We will simply write · when g is the uniform probability (g ≡ 1).
• A will generically be a rotation matrix in SO 3 (R) and a ij its (i, j) coefficient. In this section we give a review of the different existing models of collective dynamics at both microscopic and mesoscopic levels and emphasise the singularity of the BGK model among them.

A review of the different IBM
The rigorous proofs of the two following theorems (Theorems 3 and 4) can be found in [21] in a more general framework.
At a microscopic level, we fix a reference frame given by the canonical basis (e 1 , e 2 , e 3 ) of R 3 . The agents are described by their position X ∈ R 3 and their body-attitude A ∈ SO 3 (R) which can be seen as a moving frame. We assume that an agent with body attitude A ∈ SO 3 (R) moves at a constant speed in the direction of the first vector of A : the instantaneous velocity of the agent is Ae 1 .
In the following we consider an increasing sequence of jump times (T n ) n such that the increments between two jumps are independent and follow an exponential law of parameter N ∈ N * (their expectation is 1/N ). The N agents are described at time t ∈ R + by their positions and body-attitudes The interactions between the agents can be modelled by the following Piecewise Deterministic Markov Process (PDMP) which has already been described heuristically in [16,22]: 1. Between two jump times (T n , T n+1 ), the systems evolves in a deterministic way: 2. At time T n+1 , a particle i ∈ {1, . . . , N } is chosen uniformly among the N particles. At time T + n+1 , the new body-orientation of particle i is sampled from the PDF M where for Z N ∈ (R 3 × SO 3 (R)) N we define the flux: and where K is a smooth observation kernel.
The following theorem describes the limiting behaviour of the laws of the particles when N → +∞ under the assumption that the empirical measure of the N processes converges (weakly) towards a smooth function f (propagation of chaos property). The equation on f can be derived as in [16,Section 4.2].
Theorem 3. Let f 0 be a probability measure on the space R 3 × SO 3 (R) and let Z N 0 ∈ (R 3 × SO 3 (R)) N be an initial state given by N independent random variables, identically distributed with law f 0 . Then for any t ∈ R + , the law f N t of any of one of the processes (X i,N t , A i,N t ) t at time t converges weakly towards the solution f t of the following BGK equation with initial condition f 0 : where J K * f is a matrix-valued function of the space variable x ∈ R 3 defined by In the previous works, the IBM were typically given as in [14] by a system of stochastic differential equations such as the following: where P T A k denotes the projection on the tangent space of SO 3 (R) at A k ∈ SO 3 (R) (see Section 3). In this case, the resulting equation when N → +∞ is a non-linear Fokker-Planck equation (see [5] for a rigorous proof in the Vicsek case).
In the spatially homogeneous case, we can take the observation kernel K to be constantly equal to 1 to prove the mean-field limit. The agents are described at time t ∈ R + only by their body-attitudes A i,N i∈{1,...,N } ∈ SO 3 (R) N and they follow the following jump process: at each jump time T n , compute the flux Theorem 4. Let A i,N 0 i∈{1,...,N } ∈ SO 3 (R) N be an initial state given by N independent random variables, identically distributed according to a law f 0 on SO 3 (R). Then for any t ∈ R + , the law f N t of any of one of the processes (A i,N t ) t at time t converges weakly towards the solution f t of the following spatially homogeneous BGK equation with initial condition f 0 :

A review of the different kinetic equations
The model studied in the present article belongs to a class of models, the study of which has been initiated in [18] as a continuum version of the Vicsek model [54]. These models can be classified in two types. First, in the Vicsek-type models, the agents are described by their orientation defined as a unit vector in S n−1 . In the second type of models, we take into account their body-orientation, defined as a rotation matrix in SO 3 (R). Our study enters into this second framework.
The kinetic version of the Vicsek-type or Body-Orientation-type models is given either by a Fokker-Planck equation or by a BGK equation. In this work we will focus on the BGK equation where The Fokker-Planck version of our model corresponds to: where ∇ A and ∇ A · are respectively the gradient and the divergence in SO 3 (R) for the Riemannian structure detailed in Section 3. Apart from the fact that the underlying interaction process [16,22] which leads to the BGK model is different from the one that leads to the Fokker-Planck model, the BGK model is structurally different and can be treated independently by using specific and simpler mathematical techniques presented in the next sections. Nevertheless, the BGK and Fokker-Planck models share important properties. For instance, the following functional is a free-energy for both the spatially homogeneous BGK equation and the spatially homogeneous Fokker-Planck equation (though with a different dissipation term): It satisfies in both cases: where D[f ] is the dissipation term which is equal for the BGK model to: In the context of the Vicsek model, this free energy was the key to study the phase transition phenomena [13] and we believe that the same kind of study can be made in the body-attitude coordination dynamics modelled by a Fokker-Planck equation (5). Moreover, in the Fokker-Planck case this dissipation inequality implies a gradient flow structure in the Wasserstein-2 distance which has been studied (in the Vicsek case) in [24]. However, the BGK model has another underlying gradient-flow dynamics (studied in Section 5) on which the present study will be based, and we will therefore not use this free-energy in the present work.
Both models (BGK and Fokker-Planck) have a normalised and a non-normalised version. The model (4) will be referred as the non-normalised BGK model. A normalised model is a model where the flux J f is replaced by the orthogonal part of its polar decomposition Λ f := P D(J f ) as defined in the introduction and under the assumption that det J f > 0. The normalised Fokker-Planck model is the model studied in [14]: This terminology comes from the continuum version of the Vicsek model [18] where either the total flux or its normalisation is considered. A mathematical analysis of the normalised Vicsek model can be found in [24,26]. The importance of this distinction in the context of phase transitions has been shown in [13] and [12]: phase transitions appear only in non-normalised models.
The following chart ( Figure 1) shows the different models and gives references where they are studied (when such references exist).
Model [18,24,26] [5, 12,13] [ 22] In progress [14,16] In progress [16] Present work Vicsek Fokker-Planck  Finally, in Sections 4 and 5, we will focus on the spatially homogeneous version of the BGK model (4) given by: where the probability distribution f (t, A) only depends on the body-orientation variable and time. In the spatially homogeneous case, the local density of agents previously denoted by ρ f does not depend on f in the sense that an initial density ρ f0 ∈ R + associated to the initial distribution f 0 is preserved by the dynamics: as it can be seen by integrating the equation over SO 3 (R). We therefore take ρ ∈ R + as a fixed parameter of the problem. Note also that the well-posedness of (7) directly follows from Duhamel's formula: since J f is given as the solution of the following differential equation on M 3 (R): as it can be seen by multiplying (7) by A and integrating over SO 3 (R). Note that it contrasts with the Fokker-Planck case where even the well-posedness of the spatially-homogeneous equation would require further investigations. This will be part of future work.
This paragraph collects the main properties of the Riemannian manifold SO n (R) and other technical results. In this paragraph n ≥ 3 denotes the dimension, we will mainly consider the case n = 3 in the next sections.
3.1 Structure and Haar measure on SO n (R) Lemma 3.1. The following is an inner product on M n (R) : and the following properties hold: • Endowed with this metric, SO n (R) is a topological group and a Riemannian manifold.
• The sets S n (R) and A n (R) of symmetric and skew-symmetric matrices are orthogonal and M n (R) = S n (R) ⊕ A n (R).
• For A ∈ SO n (R), the tangent space to SO n (R) at A is denoted by T A and M ∈ T A if and only if there exists P ∈ A n (R) such that M = AP.
The norm on M n (R) associated to the inner product (8) will be denoted by · .
The general theory of locally compact topological groups ensures the existence of a Haar measure µ on SO n (R) which satisfies for all P ∈ SO n (R) and all Borel set E of the Borel σ-algebra of SO n (R): where P E = {P A, A ∈ E} and EP = {AP, A ∈ E}. We will assume that µ is the unique Haar measure which is a probability measure and simply write As a consequence if P ∈ SO n (R), A → P A and A → AP are two changes of variable with unit Jacobian. We will constantly use the following changes of variable : Definition 3.1 (Useful changes of variable). Let us define the following matrices: • For i = j ∈ {1, . . . , n}, D ij ∈ SO n (R) is the diagonal matrix such that all its coefficients are equal to 1 except at positions i and j where they are equal to −1.
The other coefficients are equal to 0. Then we define the following changes of variable with unit Jacobian: • A = D ij A multiplies the rows i and j by −1. Everything else remains unchanged.
• A = AD ij multiplies the columns i and j by −1. Everything else remains unchanged.
Everything else remains unchanged • A = P ij A multiplies row i by −1 and permutes the rows i and j.
• A = P ij A(P ij ) T exchanges the diagonal coefficients (i, i) and (j, j) (and involves other changes).
The two following lemmas are important applications of these results.
Proof. Let k = and m = k, . The change of variable A → D km AD km gives: Lemma 3.3. For any n ≥ 3 and any J ∈ M n (R), for given a, b, c ∈ R depending on g and on the dimension, the expressions of which can be found in the proof.
The proof of these lemmas and other technical results about SO 3 (R) and SO n (R) are postponed to Appendix B.

Volume forms in SO 3 (R)
When an explicit calculation will be needed, we will use one of the two following parametrisations of SO 3 (R) which give two explicit expressions of the normalised Haar measure in dimension 3.
• To a matrix A ∈ SO 3 (R) there is an associated angle θ ∈ [0, π] and a vector n ∈ S 2 such that A is the rotation of angle θ around the axis n. Rodrigues' formula gives a representation of A knowing θ and n = (n 1 , n 2 , n 3 ) : where and we have: With the usual parametrisation of the sphere S 2 we can take n = (n 1 , n 2 , n 3 ) T with    n 1 = sin ψ cos ϕ, n 2 = sin ψ sin ϕ, where ψ ∈ [0, π] and ϕ ∈ [0, 2π]. The volume form for the sphere is given by: • We have the following one to one map : where and φ 2 ∈ [0, π], we define: The matrix A a performs an arbitrary rotation of the first 2 coordinates and the matrix M (p) ∈ SO 3 (R) maps the vector e 3 to p ∈ S 2 . A matrix A ∈ SO 3 (R) can thus be written as the product: With this parametrisation: and the volume form on the sphere is given by: This parametrisation can be extended in any dimension and comes from the Lie groups quotient:

Singular Value Decomposition (SVD)
We recall the following classical result proved in [50, Section 1.9]. where P, Q ∈ O n (R) and D diagonal with nonnegative coefficients listed in decreasing order.
In order to use the properties of the Haar measure, we will need the matrices P and Q to belong to SO 3 (R) (not only O 3 (R)) and we define therefore another decomposition, called the Special Singular Value Decomposition (SSVD) in the following.
The existence of a SSVD follows from Proposition 3.1. Let us start from a SVD M = P D Q .
• If det M > 0, either P , Q ∈ SO 3 (R) and the SVD is a SSVD or P , Q have both negative determinant and in this case we can take . Assume without loss of generality that Q ∈ SO 3 (R). Then we can take: Remark 3.1. As for the polar decomposition and the standard SVD, the matrix D is always unique. However the matrices P and Q may not be unique.
The subset D ⊂ M 3 (R) of the diagonal matrices which are the diagonal part of a SSVD is the cone delimited by the image by the isomorphism diag of the three planes {d 1 = d 2 }, {d 2 = d 3 } and {d 2 = −d 3 } in R 3 and depicted in Figure 3 :

Equilibria of the BGK operator
In this section we determine the equilibria for the BGK operator: that is to say the distributions f such that Q BGK (f ) = 0. In Section 4.1 we characterise these equilibria (Theorem 5) and show that for them to exist, compatibility equations must be fulfilled. These compatibility equations depend on the density ρ. Therefore, for different values of the density ρ, there exists different equilibria. These will be determined in Section 4.2 by studying the compatibility equations. A full description of the equilibria of the BGK operator is finally given in Corollary 4.2.

Characterisation of the equilibria and compatibility equations
The main result of this section is Theorem 5 which gives all the equilibria of the BGK operator (13). Before stating and proving it we will need the following lemma which is the analog of lemma 4.4 in [14]. The proof of this lemma is an application of the results presented in Section 3. (i) There exists a function c 1 = c 1 (α) defined for all α ∈ R such that for all Λ ∈ SO 3 (R), The function c 1 can be explicitly written c 1 (α) = 1 3 (2 cos θ + 1) α where {·} α denotes the mean with respect to the probability density There exists a function c 2 = c 2 (α) defined for all α ∈ R such that for all B ∈ B, The function c 2 can be explicitly written: Remark 4.1. The relevance of the set B will become apparent in Proposition 4.2.
Proof. (i) Using the left invariance of the Haar measure, it is enough to prove the result for Λ = I 3 , since When Λ = I 3 , Lemma 3.2 first ensures that A M αI 3 is diagonal, then the change of variable A = P 12 A(P 12 ) T (see Definition 3.1) shows that: Proceeding analogously with the other coefficients we have that A M αI 3 is proportional to The parametrisation of SO 3 (R) using Rodrigues' formula (9) then gives the explicit expression of c 1 by taking the trace in Equation (18) and using that for A = A(θ, n), Tr(A) = 2 cos θ + 1.
(ii) As before, using the left and right invariance of the Haar measure it is enough to prove the result for  (16) holds. The parametrisation of SO 3 (R) coming from the isomorphism (10) then gives the explicit expression of c 2 by taking B = diag(1, 0, 0) in Equation (16). First, using the change of variable A → P 13 A(P 13 ) T it holds that, Then, using the parametrisation (11), it follows that: We could alternatively use one of the two parametrisations of SO 3 (R) given in Section 3.2 or the quaternion formulation to prove that A αI3 and A αB are proportional to I 3 and B. However, the proof that we have just presented here holds in any dimension (the value of the constants c 1 (α) and c 2 (α) depends on the dimension but not the form of the matrices) whereas the volume forms and the quaternion formulation strongly depend on the dimension n = 3.
We can now state the main result of this section: Theorem 5 (Equilibria for the homogeneous Body-Orientation BGK equation). Let ρ ∈ R + be a given density. The equilibria of the spatially homogeneous BGK equation (7) are the distributions of the form f = ρM J where J ∈ M 3 (R) is a solution of the matrix compatibility equation: The solutions of the compatibility equation (19) are: 1. the matrix J = 0, 2. the matrices of the form J = αΛ with Λ ∈ SO 3 (R) and where α ∈ R satisfies the scalar compatibility equation 3. the matrices of the form J = αB where B ∈ B and where α ∈ R satisfies the scalar compatibility equation where the set B and the functions c 1 and c 2 are defined in Lemma 4.1.

Remark 4.3.
Notice that the existence of a non-zero solution for the scalar compatibility equations (20) and (21) is not guaranteed for all values of ρ > 0 . The existence of nonzero solutions for these equations will be explored in Section 4.2. They will determine the existence of equilibria for Equation (7) for a given value of ρ (Corollary 4.2).
Remark 4.4. The fact that these matrices are solutions of the matrix compatibility equation (19) follows directly from the consistency relations (14) and (16) as it will be shown in the proof of Theorem 5. The main difficulty of the proof is therefore the necessary condition: we will prove that a solution of the matrix compatibility equation (19) is necessarily of one of the forms listed in Theorem 5.
The proof of this theorem will use the two following propositions. The first one and its corollary (Proposition 4.1 and Corollary 4.1) show that the compatibility equation (19) can be reduced to a compatibility equation on diagonal matrices (equation (22)). The second one (Proposition 4.2) provides a necessary condition for a diagonal matrix to be a solution of (22). The proof of Proposition 4.2 is deferred to the next section.
Proposition 4.1 (Orbital reduction). The following equivalence holds: J ∈ M 3 (R) is a solution of the matrix compatibility equation (19) if and only if for all J ∈ Orb(J), J is a solution of the matrix compatibility equation (19).
Proof. This is a consequence of the left and right invariance of the Haar measure which ensures that for any J ∈ M 3 (R) and any P, Q ∈ SO 3 (R) : Since the diagonal part of the SSVD of a matrix J is in the orbit of J, we obtain the following corollary: The following equivalence holds: J is a solution of (19) if and only if D is a solution of (19).
We will therefore consider only the following problem in dimension 3: find all the diagonal where the set D is the subset of diagonal matrices which are the diagonal part of a SSVD and is defined by (12). Notice that Equation (22) is just Equation (19) restricted to the set D.
Remark 4.5. The diagonal part D ∈ M 3 (R) of a SSVD of a matrix J ∈ M 3 (R) is unique so the problems (19) and (22) are equivalent. Notice that there might be other diagonal matrices in Orb(J) (take for example J diagonal which does not satisfy the conditions (12)). However the diagonal part of any SSVD of these matrices is D : the diagonal part of the SSVD characterises the orbit of a matrix. In the following, we will find all the diagonal solutions of (19) (i.e. the solutions of (22) without the restriction D ∈ D) and then only consider the ones which belong to D. For instance we will see that there are solutions of (19) of the form diag(0, −α, 0) where α > 0. The diagonal part of their SSVD is diag(α, 0, 0) and is a solution of (22).
Remark 4.6. A diagonal solution D of the matrix compatibility equation (19) verifies that D/ρ belongs to the set: where P(SO 3 (R)) is the set of probability measures on SO 3 (R). The set diag −1 (Ω) ⊂ R 3 is exactly the tetrahedron T defined as the convex hull of the points (±1, ±1, ±1) with an even number of minuses (which we will call Horn's tetrahedron). It is a consequence of Horn's theorem [39,Theorem 8] which states that T is exactly the set of vectors which are the diagonal of an element of SO 3 (R). It ensures that if f is a probability measure, we have by convexity of T : and therefore diag −1 (Ω) ⊂ T . Conversely, taking the Dirac deltas δ I3 and similarly for the other vertices of T , we see that the four vertices of Horn's tetrahedron belong to diag −1 (Ω).
The diagonal solutions of the matrix compatibility equation (19) satisfy the following necessary condition.
Proposition 4.2. The diagonal solutions of the compatibility equation (19) are necessarily of one of the following the types : with an even number of minus signs and where α ∈ R \ {0}.
If α ∈ (0, +∞), the diagonal part of the SSVD of these diagonal matrices is equal to D = αI 3 . If α ∈ (−∞, 0), the diagonal part of the SSVD of these diagonal matrices is equal to (c) D = α diag(±1, 0, 0) and the matrices obtained by permutation of the diagonal coefficients and where α ∈ R \ {0}. The diagonal part of the SSVD of these diagonal matrices is equal to D = diag(|α|, 0, 0). Section 4.2 will be devoted to the proof of this proposition. We are now ready to prove Theorem 5.

Proof (of Theorem 5). An equilibria of the BGK equation is of the form
It is straightforward to check that J = 0 is a solution of (19). Now, let D a matrix of one the form described in Proposition 4.2 with a parameter α ∈ R. For instance, for a matrix of type (c) like D = α diag(0, −1, 0), thanks to Lemma 4.1 we have: Similarly for the other diagonal matrices of type (c), we prove that they are solution of the matrix compatibility equation (19) if and only if their parameter α ∈ R is solution of the scalar compatibility equation (21). Analogously one can check that the diagonal matrices of type (b) are solutions of the matrix compatibility equation (19) if and only if their parameters α ∈ R are solutions of the scalar compatibility equation (20). This yields all the diagonal solutions of (19). Now, the solutions of (19) are exactly the matrices J ∈ Orb(D) where D is a diagonal solution of (19)

Proof of Proposition 4.2
The proof of Proposition 4.2 is based on two results. The first one has been proved in [55,Section 4] to study the nematic alignment of polymers in higher dimensional spaces: Theorem 6 ([55]). Let n ≥ 3, b ∈ R + and s = (s 1 , s 2 , . . . , s n ) ∈ R n a solution of the nonlinear system s j = m 2 j g s,b , j = 1, . . . , n, where the average is taken with respect to the PDF on the sphere S n−1 : where Z is the normalisation constant which ensures that g s,b is a PDF on the sphere S n−1 .
The second tool that we will use to prove Proposition 4.2 is an isomorphism between SO 3 (R) and the space of unitary quaternions which transforms the compatibility equation (22) into the compatibility equation (23) studied in Theorem 6. 1. There is an isomorphism between the group SO 3 (R) and the quotient group H/ ± 1, where H is the group of unit quaternions. Since H is homeomorphic to S 3 , there is an isomorphism Φ : Moreover Φ is an isometry in the sense that it maps the volume form of S 3 /±1 (defined as the image measure of the usual measure on S 3 by the projection on the quotient space) to the volume form on SO 3 (R): for all measurable function f on SO 3 (R), 2. There is a linear isomorphism between the vector space M 3 (R) and the vector space S 0 4 (R) of trace free symmetric matrices of dimension 4: such that for all J ∈ M 3 (R), and q ∈ H/ ± 1, The first dot product is defined by Equation (8) and the second one is the usual dot product in R 4 .

For all
The proof of this proposition can be found in appendix A. We are now ready to prove Proposition 4.2.

Proof (of Proposition 4.2). Using the first and second points of Proposition 4.3, it holds that
The compatibility equation (22) then becomes: Applying the isomorphism φ defined in Proposition 4.3 to this last equation, we obtain thanks to the third point of Proposition 4.3 : Using the fourth point of Proposition 4.3, we then obtain the following equivalent problem: find all the trace-free diagonal matrices Q = diag(s 1 , s 2 , s 3 , s 4 ) of dimension 4 such that where Z is a normalisation constant: Equivalently, defining for i ∈ {1, 2, 3, 4} : we want to solve the system of compatibility equations: where s = (s 1 , s 2 , s 3 , s 4 ) and g s ,2ρ is given by (24). Thanks to Theorem 6, we conclude that if s is a solution of (25), then the coefficients s 1 , s 2 , s 3 , s 4 can take at most two distinct values. So, the same result holds for the coefficients s 1 , s 2 , s 3 , s 4 . Now thanks to the fourth point of Proposition 4.3, we only have the following possibilities: • if s 1 = 3α/4 and s 2 = s 3 = s 4 = −α/4 for α ∈ R, then • if s 2 = 3α/4 and s 1 = s 3 = s 4 = −α/4 for α ∈ R, then

Determination of the equilibria for each density ρ
In Theorem 5 we saw that the BGK operator can have three types of equilibria. The uniform equilibria f = ρ (corresponding to J = 0) is always an equilibrium. However, the existence of the other two types of equilibria depends on Equations (20) and (21) having a solution for a given ρ. Therefore the existence of these types of equilibria will depend on the value of ρ. In this section we will determine the existing equilibria for each value of ρ. In particular, we will draw the phase diagram for ρ and α, that is to say the parametrised curves defined by Equations (20) and (21) in the plane (ρ, α) (see Figure 2). We first prove the following proposition. (i) The function α → α/c 1 (α) is well-defined on R, its value at zero is ρ c . Moreover, there exists α * > 0 such that this function is decreasing on (−∞, α * ] and increasing on [α * , +∞). Defining ρ * := α * /c 1 (α * ), it holds that ρ * < ρ c .
(iii) We have the following asymptotic behaviours: Proof. The idea of the proof is taken from [55].
an integration by parts shows that: It proves that α and c 1 (α) have the same sign for all α ∈ R. Then we define function m : α → {sin 2 θ} α = 3c 1 (α)/α which satisfies the property: where Var α is the variance for the probability density (15). This property implies that α/c 1 (α) has only one critical point which is a global minimum. This minimum is attained at a point α * > 0 as a simple computation shows that m (0) > 0 and consequently ρ * < ρ c . A simple computation gives m(0) = 1 2 so ρ c = 6. (ii) We have similarly: from which we can easily see that α → α/c 2 (α) is even and has only one minimum attained at α = 0. A simple computation shows that its value at 0 is ρ c = 6.
(iii) The behaviour at infinity is obtained by Laplace's method: with the change of variable s = 1 − cos θ on [0, π], we get With the same method we have: Thanks to Proposition 4.4 and Theorem 5 we can now fully describe the equilibria of the BGK operator. A graphical representation of this result is given by the phase diagram depicted in Figure 2 : Corollary 4.2 (Equilibria of the BGK operator, depending on the density ρ). The set of equilibria of the BGK operator (13) depends on the value of ρ. In particular we need to distinguish three regions ρ ∈ (0, ρ * ), ρ ∈ (ρ * , ρ c ) and ρ > ρ c where ρ * and ρ c are defined in Proposition 4.4. We have the following equilibria in each region: • For 0 < ρ < ρ * , α = 0 is the unique solution of Equations (20) and (21) and therefore the only equilibrium is the uniform equilibrium f eq = ρ.
When an equilibrium is of the form f eq = ρM αΛ with parameters α > 0 and Λ ∈ SO 3 (R) then these parameters can respectively be interpreted as a concentration parameter and a mean body-orientation. They are analogous to the equilibria found in [13] in the Vicsek case. However, in SO 3 (R), there exist other equilibria which are not of this form. We will see in Section 5 that these latter equilibria are always unstable.  (13). Depending on the density, there are one, two, three or four branches of equilibria (α 2 and −α 2 give the same orbit). The uniform equilibrium f eq = ρ is always an equilibrium (corresponding to α = 0, depicted in green). The equilibria of the form f eq = ρM αΛ , Λ ∈ SO 3 (R) exist for ρ > ρ * and correspond to the two branches of the red curve α = ρc 1 (α). Finally the equilibria of the form f eq = ρM αB , B ∈ B exist for ρ > ρ c and correspond to the two branches of the blue curve α = ρc 2 (α). The dotted and dashed lines correspond to unstable equilibria (as shown in Section 5). The signs are the signature of the Hessian matrix Hess V (D) defined in Section 5 taken at an equilibrium point. The elements α * , ρ * and ρ c are defined in Proposition 4.4; the elements α + , α − , α 1 , α 2 and α 3 are given in Corollary 4.2.
Finally the following picture ( Figure 3) is a representation in the space R 3 of the diagonal parts of the SSVDs of the solution of the matrix compatibility equation (19) when ρ > ρ c . They all belong to the domain D defined by (12) and depicted in orange in Figure 3.

Convergence to equilibria
Now that we know all the equilibria of the spatially homogeneous BGK equation (7) we proceed to investigate the asymptotic behaviour of f (t, A) as t → +∞. This problem can be reduced to looking at the asymptotic behaviour of J f since, if J f → J ∞ ∈ M 3 (R), then f (t) will converge as t → +∞ towards ρM J∞ as it can be seen by writing Duhamel's formula for equation (7) : The asymptotic behaviour of J f is much simpler than the one of f since J f is the solution of the following ODE as it can be seen by multiplying (7) by A ∈ SO 3 (R) and integrating over SO 3 (R). Since J ∈ M 3 (R) → M J ∈ L ∞ (SO 3 (R)) is locally Lipschitz, the flow of Equation (28) is defined globally in time as a bounded Lipschitz perturbation of the linear system d dt J = −J.
Notice that the solutions of the compatibility equation (19) are exactly the equilibria of the dynamical system (28). We therefore obtain the following proposition: Proposition 5.1 (Equilibria of the BGK operator, equilibria of the ODE (28)). A distribution f eq = ρM J is an equilibrium of the BGK operator (13) if and only if J ∈ M 3 (R) is an equilibrium of the dynamical system (28).
We will call stable/unstable an equilibrium of the BGK operator (13) such that the associated matrix J ∈ M 3 (R) is a stable/unstable equilibrium of the ODE (28). This section is devoted to the proof of the following theorem: Theorem 7 (Convergence towards equilibria). Let ρ ∈ R + be such that ρ = ρ * and ρ = ρ c (as defined in Proposition 4.4). Let f 0 be an initial condition for (7) and let J f0 = P D 0 Q be a SSVD. Let f (t) be the solution at time t ∈ R + of the spatially homogeneous BGK equation (7) with initial condition f 0 . Let D(t) be the solution at time t ∈ R + of the ODE (28) with initial condition D 0 ∈ D. It holds that: is a SSVD and there exists a subset N ρ ⊂ R 3 of zero Lebesgue measure such that: converges as t → ∞ towards an equilibrium f eq of the BGK operator (13) of the form f eq = ρM J eq , where J eq ∈ M 3 (R) is of one of the forms described in Theorem 5. The convergence is locally exponentially fast in the sense that there exists constants δ, K, µ > 0 such that if J f0 − J eq ≤ δ then for all t > 0, 2. If D 0 / ∈ N ρ , we have the following asymptotic behaviours depending on the density ρ : (i) if 0 < ρ < ρ * , then D(t) → 0 as t → +∞ and, consequently, f eq = ρ, (ii) if ρ * < ρ < ρ c , then D(t) → 0 or D(t) → α + I 3 as t → +∞ and, consequently, f eq = ρ or f eq = ρM α+Λ respectively, where α + > 0 is defined in corollary 4.2 and Λ := P Q ∈ SO 3 (R), (iii) if ρ > ρ c , then D(t) → α 1 I 3 as t → +∞ and, consequently, f eq = ρM α1Λ where α 1 > 0 is defined in corollary 4.2 and Λ := P Q ∈ SO 3 (R).
Remark 5.1. The subset N ρ ⊂ R 3 depends on the density ρ. This subset will be made explicit in the three cases ρ < ρ * , ρ * < ρ < ρ c and ρ > ρ c in Section 5.3. If D 0 ∈ N ρ , then there is convergence towards an unstable equilibrium at a rate which may not be exponential.
Remark 5.2 (Phase transitions). Theorem 7 demonstrates a phase transition phenomenon triggered by the density of agents ρ : when ρ < ρ * , the system is disordered (asymptotically in time) in the sense that J f → 0 and we therefore cannot define a mean body-attitude. When the density increases and exceeds the critical value ρ c , the system is self-organised (asymptotically in time and for almost every initial data), in the sense that J f → αΛ where α ∈ R + and Λ ∈ SO 3 (R) can be respectively interpreted as a concentration parameter and a mean body attitude. When ρ * < ρ < ρ c the self-organised and disordered states are both asymptotically stable and the convergence towards one or the other state depends on the initial data. Such "transition region" also appears in the Vicsek model, as studied in [13], and gives rise to an hysteresis phenomenon.
The proof of this theorem can be found at the end of Section 5.2. It is based on a gradient flow structure for the flux J f studied in Section 5.1. This structure ensures the convergence of J f towards a matrix J eq ∈ M 3 (R) as t → +∞ and consequently the convergence of f (t) as t → +∞ towards an equilibrium. The stability of the equilibria determines which equilibrium can be attained. This question is addressed in Section 5.2. Additional details about the subset N ρ as well as the study of the critical cases ρ = ρ * and ρ = ρ c are provided in Section 5.3.

A gradient-flow structure in R 3
In this section we show that the ODE (28) can be reduced to a gradient-flow ODE in R 3 . We first show (Section 5.1.1) how (28) can be reduced to an ODE in R 3 , the equilibria of which are linked to the equilibria of (28) (and therefore of (7)). Then we show that this ODE in R 3 has a gradient-flow structure which will allow us to conclude on the asymptotic behaviour of the solution of (7) (Section 5.1.2).

Reduction to a nonlinear ODE in R 3 and equilibria
The ODE (28) is a matrix-valued nonlinear ODE (in dimension 9) but, as in the previous section (Proposition 4.1 and Corollary 4.1), we will use the left and right invariance of the Haar measure and the SSVD to reduce the problem to a vector-valued nonlinear ODE in dimension 3, as explained in Proposition 5.2 and Corollary 5.1. Proof. (i) Using the left and right invariance of the Haar measure, we see that if the matrix J(t) ∈ M 3 (R) is a solution of (28), then for any P, Q ∈ SO 3 (R), P J(t)Q is also a solution (and conversely).
For any matrix J 0 , such a diagonal matrix D 0 always exists: we can take the diagonal part of the SSVD of J 0 . We therefore only have to study the following ODE for diagonal matrices D 3 (R). Since this vector space is isomorphic to R 3 through the isomorphism diag defined in the introduction, we obtain two equivalent ODEs: where D 0 = diag −1 (D 0 ) and the following equivalence holds : D(t) = diag −1 (D(t)) is the solution of (29b) if and only if D(t) = diag( D(t)) is the solution of (29a).
Note that it is not clear that if J 0 = P D 0 Q is a SSVD, then J(t) = P D(t)Q is a SSVD for all t > 0. The following proposition and corollary ensure that the SSVD is preserved by the dynamical system which will allow us to restrict the domain on which the ODE (29a) is posed.

Proposition 5.3 (Invariant manifolds).
The following subsets of R 3 are invariant manifolds of the dynamical system (29b) : • the intersections of two of these planes and in particular the lines Proof. For i = 2 and j = 3, the result has already been proved in the second point of Lemma 4.1. The other cases are similar.
Corollary 5.1. Let J 0 ∈ M 3 (R) and J 0 = P D 0 Q be a SSVD with P, Q ∈ SO 3 (R) and D 0 ∈ D 3 (R) diagonal. Let J(t) be the solution of the ODE (28) with initial condition J(t = 0) = J 0 . Let D(t) the solution of the ODE (29a) with initial condition D(t = 0) = D 0 . Then the decomposition J(t) = P D(t)Q is a SSVD for J(t).
Proof. The fact that J(t) = P D(t)Q is a consequence of the first point of Proposition 5.2. The fact that it is a SSVD is a consequence of Proposition 5.3 which ensures that the conditions (12) remain true for all t > 0 : D is stable in the sense that if D 0 ∈ D, then D(t) ∈ D for all time t > 0. This follows from the fact that the image by the isomorphism diag of the invariant manifolds of (29b) described in Proposition 5.3 are invariant manifolds of the dynamical system (29a). These manifolds in D 3 (R) form the boundary of the subset D.
In conclusion, the study of the asymptotic behaviour of f (t) as t → +∞ can be reduced to the study of the asymptotic behaviour of the solutions of the ODE (29a) posed on the domain D (see Figure 3).
The following Proposition describes the equilibria of the dynamical system (32a) and is a consequence of the results of Section 4.
Proposition 5.4 (Equilibria of the dynamical system (29a)). The equilibria of the dynamical system (29a) are, depending on the density ρ : • the matrix D = 0 for any density ρ, • the diagonal matrices of type (b) with parameters α + and α − when ρ * < ρ < ρ c , Proof. The equilibria of (29a) are the diagonal matrices D such that which is exactly Equation (19) for the diagonal matrices. This equation has been solved in Theorem 5, Remark 4.7 and Corollary 4.2.
Remark 5.3. As in the previous section, we have found all the diagonal equilibria of (29a). However, thanks to Corollary 5.1, only the ones which belong to D are needed.
The following proposition is a straightforward consequence of the previous results. (1) The distribution f eq = ρM J is an equilibrium of the BGK operator (13).
(2) The matrix D is an equilibrium of the dynamical system (29a) on the domain D.

A gradient-flow structure
Lemma 5.1 (Gradient-flow structure). We define the partition function of a matrix J ∈ M 3 (R) : and the potentials V (J) and V ( D) respectively on M 3 (R) and on R 3 : where · and | · | are the Euclidean norms respectively on M 3 (R) and on R 3 . Then we can rewrite equations (29a) and (29b) into a gradient flow structure as follows: where ∇ is the gradient operator in M 3 (R) endowed with the Riemaniann structure (8) or the gradient operator in R 3 endowed with the usual Euclidean structure.
Proof. The partition function satisfies that for all J ∈ M 3 (R), since ∇(e J·A ) = Ae J·A . The result follows in M 3 (R). The result in R 3 follows from the fact that for any w 1 , w 2 ∈ R 3 , it holds that: where · denotes the dot product on R 3 and on M 3 (R) as defined in (8).
Remark 5.4. This gradient-flow structure on J f is specific to the BGK equation and does not hold for the Fokker-Planck operator (as shown in [13], the differential equation satisfied by J f in the Vicsek case involves the spherical harmonics of degree 2 and higher of f ; here the equation for J f is closed).
In particular the gradient-flow structure (32a) implies that the system (32a) will converge towards an equilibrium. When all the equilibria of the dynamical system (29a) are hyperbolic the convergence towards the stable equilibria is exponentially fast (see [38,Section 9.3]). The goal of the next section is to find which equilibria among the ones found in Section 4 are stable and to completely describe the asymptotic behaviour of the system depending on the initial condition and the local density ρ. We will see that phase transitions appear between ordered and disordered dynamics when the density ρ increases.

Stability of the equilibria and conclusion
Since the flow of (29b) is the image by the isomorphism diag −1 of the flow of (29a), an equilibria D eq of (29a) is stable (resp. unstable) if and only if diag −1 (D eq ) is a stable (resp. unstable) equilibria of (29b). The stability properties of the equilibria of (29b) are much simpler to study than the stability properties of the equilibria of (29a) since they are given by the signature of the Hessian matrix Hess V ( D eq ) ∈ S 3 (R) of the potential V given in equation (31b) (the linearisation of ODE (29b) around an equilibrium D eq is indeed d dt H = − Hess V ( D eq ) H). In particular, an equilibrium D eq is stable if and only if the signature of Hess V ( D eq ) is (+ + +).
Note that in the matrix framework (32a), the Hessian of the potential (31a) in the Euclidean space M 3 (R) would be a rank 4 tensor. Here we are reduced to the computation of the signature of 3×3 symmetric matrices. For a diagonal matrix D ∈ D 3 (R) and D = diag −1 (D) we will write with a slight abuse of notations: When D ∈ D 3 (R) is an equilibrium of (29a) we call signature of the Hessian matrix Hess V (D) the signature of Hess V ( D) where D = diag −1 (D). A simple computation shows that the Hessian matrix Hess V (D) is given by : The following theorem is an extension of Corollary 4.2 and gives a full description of the equilibria of the BGK operator with their stability.
Theorem 8 (Stability of the equilibria of the ODE (29a)).
• For ρ * < ρ < ρ c , the equilibrium D = 0 and the equilibria of type (b) with parameter α + are stable. The equilibria of type (b) with parameter α − are unstable and the signature of the Hessian matrix is (− + +). The proof of Theorem 8 will be based on the following lemma which states an important orbital invariance principle for the signature of the Hessian matrix. Proof. Let D 1 and D 2 be two equilibria of (29b) such that diag( D 1 ) and diag( D 2 ) are of the same type. Then there exist P, Q ∈ SO 3 (R) and R ∈ O 3 (R) such that Proof (of Theorem 8). Thanks to lemma 5.2, we therefore do not have to compute the signatures of the Hessian matrix taken at all equilibria, it is enough to choose one matrix in each orbit: for the equilibria of type (b) (see Proposition 4.2) we will compute the signature of Hess(αI 3 ) where α = ρc 1 (α) and for the equilibria of type (c) we will compute the signature of Hess(α diag (1, 0, 0)) where α = ρc 2 (α). We first start with the case of the uniform equilibrium.
Case 1. Uniform equilibrium D = 0 For the uniform equilibrium D = 0, a ij = 0 for all (i, j). Moreover, by the change of variable A = D ik A where k = i, j, we have when i = j : It proves that Γ is diagonal. Then with the changes of variables A = P ij A or A = AP ij it can be seen that all the 3 2 quantities a 2 ij are equal. Since their sum is equal to n = 3 we get that where ρ c = 6. In conclusion, the signature of Hess V (0) is ( • 1 − 1 2 ρν − ργ of order 1 with eigenvector (1, 1, 1) T . But taking the derivative with respect to α of (14) with Λ = I 3 we obtain the relation : and using ρ = α/c 1 (α) we can rewrite : Its sign is then given by Proposition 4.4 and the fact that c 1 (α) has the same sign as α : c 1 (α) id c1 (α) > 0 when α < 0, c 1 (α) id c1 (α) < 0 when 0 ≤ α < α * and (α) > 0 when α > α * .
We get that : We have : Linearizing sin 2 (θ/2) and expanding everything gives : so that : Now for x ≥ 0, let θ 0 = arccos(−2/3). We cut the integral at θ 0 and we get The Hessian matrix is therefore equal to : with eigenvector (0, 1, 1) T . We have as before : so the first eigenvalue can be rewritten where we have used Proposition 4.4 to determine the sign. The two other eigenvalues are equal to: Using the change of variable A → P 13 A(P 13 ) T and the parametrisation (11), one can see that : The conclusion of the proof follows from the study of the roots of Equations (20) and (21) provided by Corollary 4.2 and depicted in Figure 2. In particular, the cases α = α * and α = 0 correspond to the critical cases ρ = ρ * and ρ = ρ c which are studied in the next section.
Remark 5.5. The above calculations are similar to the computation of the equilibria and their stability of where F is the free-energy (6). In particular, it can be shown that the equilibria of Ψ and their stability are the same as the ones described in Theorem 8. In particular, since (6) is also a free energy in the Fokker-Planck case, this analysis shows the instability of equilibria of the Fokker-Planck of the form ρM D when D is of one of the unstable equilibria of (29a) described in Theorem 8. This technique is similar to the one which was used in [57] in the case of 3D polymers. However, it does not provide global or local convergence of the solution of the Fokker-Planck equation towards an equilibrium. This requires further investigations which will be left for future work. In the Vicsek case [13], it was based on a LaSalle's principle and on estimates for the dissipation term.
We can finally prove Theorem 7 : Proof (of Theorem 7). Thanks to the Duhamel's formula (27), the asymptotic behaviour of f (t) as t → +∞ is given by the asymptotic behaviour of J f (t) . Thanks to Proposition 5.2 and Corollary 5.1, we only have to study the asymptotic behaviour of the solution of the ODE (29a) where the initial condition is the diagonal part of a SSVD of J f0 . Since this equation has a gradient-flow structure (32a), we know that the solution D(t) converges as t → +∞ towards an equilibrium D eq and consequently J f (t) → P D eq Q := J eq . Moreover the equilibrium D eq is a stable equilibrium provided that D 0 does not belong to the stable manifold of an unstable equilibrium. Since these manifolds have dimension at most 2, the union of these manifolds, called N ρ is of zero measure and there is convergence towards a stable equilibrium for all D 0 / ∈ N ρ . In this case, the convergence is locally exponentially fast in the sense that there exist constants δ, λ, C > 0 such that if J f0 − J eq ≤ δ then for all Let f eq = ρM J eq . It follows from Duhamel's formula (27) that for all A ∈ SO 3 (R) : Since SO 3 (R) is compact and J f (t) is bounded uniformly in t, there exists a constant L > 0 such that for all t > 0, the following Lipschitz bound holds: Reporting (35) into (34) and using (33), the first point of Theorem 7 follows with constants K = CL/|1 − λ| and µ = min(1, λ) when λ = 1 and K = CL and µ = 1 − ε for any ε > 0 when λ = 1.
The stability of the equilibria of the dynamical system (29a) is given in Theorem 8. Finally the conclusions of points 2.(ii) and 2.(iii) follow from the fact that the diagonal parts of the SSVD of the equilibria of type (b) with parameters α + and α 1 are respectively α + I 3 and α 1 I 3 .

Final remark: characterisation of the stable manifold of the unstable equilibria and critical cases
This section is devoted to a more precise description of the subset N ρ of zero measure of initial conditions which do not necessarily lead to one of the behaviours detailed in the previous section (see Theorem 7).
When ρ = ρ * and ρ = ρ c , all the equilibria of (29b) are hyperbolic and we can apply the stable manifold theorem ([48] Section 2.7) which states that for a given hyperbolic equilibrium D eq of the nonlinear equation (29b), the stable set of D eq is a smooth manifold. Its dimension is given by the number of minuses in the signature of − Hess V ( D eq ) and its tangent space at D eq is the stable subspace of the linearized system around D eq . We are now ready to describe the behaviour of the system (29b) for a given density ρ and an initial condition D 0 ∈ R 3 . Note that we can restrict ourselves to the case D 0 ∈ diag −1 (D) and that our goal is to describe diag −1 (N ρ ) (in the R 3 -framework). We write D(t) = diag D(t) the solution of (29a) with initial condition D 0 = diag D 0 . Case 1. When ρ < ρ * When ρ < ρ * , there is only the uniform equilibrium : D(t) −→ t→+∞ 0 and N ρ = ∅.
• The uniform equilibrium is unstable and − Hess( V )(0) has signature (+ + +) therefore, there is no stable direction and D(t) cannot converge towards 0 unless D 0 = 0 (and then D(t) = 0 for all t ∈ R + ).
• The equilibrium α 3 diag(−1, −1, 1) is unstable and its stable manifold M 1 has dimension 1: it is the half line R * • The equilibrium α 2 (1, 0, 0) T is unstable and its stable manifold has dimension 2 with a tangent plane generated by the two vectors (1, 0, 0) T and (0, 1, −1) T : it is the plane {d 2 + d 3 = 0}. We note that this plane is an invariant manifold, so if D 0 belongs to this plane, the limit when t → +∞ also belongs to this plane. But since the half-lines , we can restrict ourselves to the eighth of plan M 2 delimited by the half lines R * + (1, 1, −1) T excluded and R * + (1, 0, 0) T included. • In every other cases, D(t) converges towards α 1 I 3 which is the only stable equilibrium.

The critical cases
We end this section by an informal description of the expected behaviour in the critical cases.
• When ρ = ρ * the uniform equilibrium is stable and there is an other equilibrium of the form α * I 3 with α * > 0. This equilibrium is non hyperbolic: the kernel of − Hess V (α * (1, 1, 1) T ) is one dimensional, spanned by the vector (1, 1, 1) T . The two other eigenvalues are negative. Therefore, for the system (29b), there exists a center manifold of dimension 1 and a stable manifold of dimension 2, which tangent plane being the orthogonal of (1, 1, 1) T .
If D 0 belongs to the stable manifold, D(t) converges exponentially fast towards α * I 3 ([33] Theorem 3.22). The stable manifold of dimension 2 delimits two domains and one of them is included in the subset {x ∈ R 3 , x · (1, 1, 1) T ≥ α * }. If D 0 belongs to this domain, then it belongs to a center manifold and D(t) is attracted exponentially fast towards the the line R(1, 1, 1) T and converges at rate 1/t towards α * (1, 1, 1). In every other cases, D(t) converges exponentially fast towards 0.
• When ρ = ρ c there is one stable equilibrium of the form αI 3 with α > 0. The uniform equilibrium is non hyperbolic and − Hess V (0) = 0.
In the case when D(t) converges towards 0, with a formal computation, the rate of convergence is expected to be either 1/t or 1/ √ t depending on D 0 . The convergence towards the stable anisotropic equilibrium is exponentially fast.

Macroscopic limit for the stable equilibria
We now go back to the spatially inhomogeneous model. We want to investigate (at least formally) the hydrodynamic models derived from the BGK equation (4). To do so, we introduce the scaling t = εt and x = εx for ε > 0 and we define f ε (t , x , A) := f (t, x, A). After this change of variables in (4) and dropping the primes, we see that f ε satisfies the following equation: We want to investigate the macroscopic limit ε → 0 with the assumption that f ε converges towards a stable equilibrium. Thanks to the results of the last section, we will assume that f ε → ρM αΛ where ρ = ρ(t, x), α = ρc 1 (α), α ∈ R + and Λ = Λ(t, x) ∈ M 3 (R) (with a notion of convergence as strong as needed). Since the equilibrium is assumed to be stable there are two cases: either Λ ∈ SO 3 (R) (and therefore ρ > ρ * ) or Λ = 0 that is to say f ε is uniform in the body-orientation variable and converges towards ρ = ρ(t, x) (and therefore ρ < ρ c ). For a given time t ∈ R + , we will say that x ∈ R 3 belongs to a disordered region when Λ(t, x) = 0. Otherwise, when Λ(t, x) ∈ SO 3 (R), we will say that x ∈ R 3 belongs to an ordered region.
The purpose of the two next sections is to write at least formally the hydrodynamic equations satisfied by ρ = ρ(t, x) and Λ = Λ(t, x). First notice that integrating (36) over SO 3 (R) leads to the conservation law: where ρ ε ≡ ρ f ε and The macroscopic model then depends on the region considered.
1. In a disordered region, j[f ε ] → 0 and assuming that the convergence is sufficiently strong, we get that ∂ t ρ = 0. To obtain more information we will look at the next order in the Chapman-Enskog expansion (Section 6.1).
2. In an ordered region, j[f ε ] → αΛe 1 where α = ρc 1 (α) and therefore, assuming that the convergence is strong enough: However due to the lack of conserved quantities, we will need specific tools to write an equation satisfied by Λ = Λ(t, x) in order to obtain a closed system of equations on (ρ, Λ). This is the purpose of Section 6.2.

Diffusion model in a disordered region
We consider a region where f ε converges as ε → 0 to a density ρ(t, x) uniform in the bodyattitude variable. The following proposition gives the diffusion model obtained by looking at the next order in the Chapman-Enskog expansion.
Proposition 6.1 (Formal). In a disordered region, the density ρ ε satisfies formally at first order the following diffusion equation: Proof. We follow the same calculations as in [13] : we write f ε = ρ ε + εf ε 1 (where f ε 1 is defined by this relation) and notice that: . Inserting this in (36), multiplying by A and integrating over SO 3 (R) leads to: Using Lemma 3.3, it holds that To compute the second term, we note that Ae 1 · ∇ x ρ ε = 2A · R ρ ε where R ρ ε is the matrix, the first column of which is equal to ∇ x ρ ε and the others are equal to zero. Using Lemma 3.3 we obtain SO3(R) By multiplying (39) by e 1 , it follows from (40) and (41) that : which gives the result by inserting this in (37).
Remark 6.1. This analysis does not depend on the dimension. In SO n (R) the same formal result holds: , ρ c = 2n.

Self-organised hydrodynamics in an ordered region
In the following, for a given density ρ ∈ R + , α(ρ) denotes the maximal nonnegative root of α = ρc 1 (α). We are going to prove the following theorem.
The first equation (42a) is the conservation law (37) and the goal is to obtain the equation (42b) for Λ = Λ(t, x). However, here and contrary to the classical gas dynamics, the total momentum is not conserved: d dt SO3(R) f A dA = 0 and we therefore cannot deduce easily a closed system of equations. This lack of conserved quantities is specific to self-propelled particle models such as the Vicsek model. The main tool to tackle the problem will be the Generalised Collision Invariants method introduced in [18] for the study of the hydrodynamic limit of the continuum Vicsek model. Its precise setting in the context of our body-attitude model is detailed Section 6.2.1. The formal proof of Theorem 9 can be found in Section 6.2.2.

Generalised collision invariants
To obtain an equation on Λ, the main tool are the Generalised Collisional Invariants (GCI) first introduced in [18]. For a given J ∈ M 3 (R), we define first the linear collision operator: with det J > 0 and let Λ := P D(J) ∈ SO 3 (R) be the orthogonal part of its polar decomposition. The set of GCI associated to J is defined as: The condition ψ ∈ C J is equivalent to: Therefore, following the ideas of the proof of [14, Proposition 4.3], we have: that is to say: or equivalently since B ∈ T Λ means that there exists P ∈ A 3 (R) such that B = ΛP : where ψ Λ P (A) = −P · Λ T A. Now for any P ∈ A 3 (R), denoting Λ f ε ≡ Λ ε = P D(J f ε ), we get by multiplying the equation (36) by ψ Λ ε P : The right-hand side vanishes by Definition (43) of the GCI. Note that if f ε → ρM J with det(J) > 0, we have also det(J f ε ) > 0 for ε sufficiently small and Λ ε ∈ SO 3 (R).

Hydrodynamic limit (formal proof of Theorem 9)
Taking formally the limit ε → 0 in (44) and since it is true for all P ∈ A 3 (R), we obtain: All the terms that involve the time derivative of ρ are equal to zero since ∂ t ρ does not depend on A and with the change of variable A = ΛA T Λ which has unit jacobian, it holds that A · Λ = A · Λ and therefore We thus have: where and With the change of variable A → Λ T A these terms become and they can be computed using the same techniques as in [14] or lemma 3.4 in the appendix. More precisely, we can write where R 1,ρ is the matrix, the first column of which is equal to 2Λ T ∇ x ρ and the others are all equal to zero. It satisfies: By the change of variable B → B T , we have This integral is of the form (46) where is invariant by transposition and conjugation and is a skew-symmetric matrix. From (48) and (49) we get The first term Y 1 can therefore be written: and using Rodrigues' formula we obtain: where {·} α has been defined in Proposition 4.1. Similarly the second term Y 2 can be written: where the coefficient C 3 is the same as in [14] : Finally, we obtain: Putting all the terms together, we can conclude as in [14]. First we notice that: Therefore we obtain by multiplying (45) by Λ and dividing by αC 2 : where the coefficients {sin 2 θ(2 + 3 cos θ)} α {sin 2 θ} α and c 4 := C 5 C 2 = 1 5 are respectively equal to the coefficients c 2 and c 4 in [14] and the coefficient c 3 in [14] (which is equal to 1) becomes: (1 + 2 cos θ) sin 2 θ α {sin 2 θ} α   .

Conclusion
In this work, we have presented a new BGK model of body-attitude coordination where agents are described by a rotation matrix. Starting from the kinetic level (a space homogeneous BGK equation) we have drawn a parallel between our Vicsek-type model and the models of nematic alignment of polymers. We then have deduced the equilibria of the system and have shown a phase transition phenomenon triggered by the density of agents. Thanks to a gradient-flow structure specific to the BGK equation we have been able to describe the asymptotic behaviour of the system. Finally, we have derived the macroscopic models (SOHB) in the spatially inhomogeneous case.
On the modelling side, a rigorous mean-field limit which leads to the BGK equation is currently under study and will be the object of future work. However, many other questions remain open. At the kinetic level, our study relies on the dimension 3 and it would be interesting to extend the ideas developed here in SO n (R), n ≥ 3, for example by drawing new parallels with higher dimensional polymers models or other similar models [27,28]. In addition, the mathematical and numerical analyses of the macroscopic SOHB model are still in progress.
The BGK model studied here is a step towards the full description of the models of collective behaviour depicted in Figure 1. The tools and ideas that we have presented here may help to analyse other models of body-attitude coordination such as the non-normalised Fokker-Planck model in SO 3 (R). Other models which take into account curvature control in addition to body-orientation, in the spirit of [19], could also be considered.

A Quaternions and rotations
These appendix is devoted to the proof of Proposition 4.3. We also give additional results about quaternions. The following lemma gives a link between quaternions and the theory of Q-tensors.
Lemma A.1. Let S 0 4 (R) be the space of symmetric 4 × 4 trace free matrices. If Q ∈ S 0 4 (R) has two eigenvalues with eigenspaces of dimensions 1 and 3, then Q can be written for a given unit quaternion q seen as a vector of R 4 . A matrix of this form is called a uniaxial Q-tensor. When α = 1 we will say that Q is a normalised uniaxial Q-tensor.
Proof. Let Q ∈ S 0 4 (R) such that Q has two eigenvalues with eigenspaces of dimensions 1 and 3. By the spectral theorem, there exists P ∈ O 3 (R) such that for a given α > 0 : Proof of Proposition 4.3.
1. The group isomorphism Φ is explicitly computed in [52]. In particular, let q ∈ S 3 / ± 1. The matrix A = Φ(q) is defined for all purely imaginary quaternion u ∈ H by A[u] = [quq * ] where (u 1 , u 2 , u 3 ) T =: [u] ∈ R 3 is the vector associated to u = u 1 i + u 2 j + u 3 k. More explicitly, if q = x + iy + zj + tk, then A =   x 2 + y 2 − z 2 − t 2 2(yz − xt) 2(xz + yt) 2(xt + yz) x 2 − y 2 + z 2 − t 2 2(zt − xy) 2(yt − xz) 2(xy + zt) Note that we have identified q ∈ H and its equivalence class in S 3 / ± 1 and that Φ is well defined since only quadratic expressions are involved. The fact that this group isomorphism is an isometry follows from [ is a quadratic form for q. We take Q the matrix associated to this quadratic form. For J = (J ij ) i,j , using the explicit form of A = Φ(q) with q = x + yi + zj + tk we obtain: ∈ Span(SO n (R)).
The SSVD (Definition 3.2) gives the result for any matrix.
Corollary B.1. For n ≥ 3, a matrix that commutes with any matrix of SO n (R) is of the form λI n .
We can now prove Lemma 3.3 and its generalisation Lemma 3.4. Therefore, Φ(P ) = λP for any P ∈ SO n (R) and the same is true for any matrix J ∈ M n (R) by lemma B.1. To compute λ, notice that for the matrix e i ⊗ e j : 1 2 SOn(R) (e i · Ae j ) 2 dA = λ.
Summing this equality for all i, j gives: 1. We first note that ψ is self-adjoint for the dot product A · B = Tr(A T B): for any K ∈ M n (R), ψ(J) · K = SOn(R) (J · A)(K · A) g(A) dA = J · ψ(K).
2. We prove that Span(I n ) is a stable supspace for ψ : for any P ∈ SO n (R) we have: P ψ(I n )P T = SOn(R) (I n · A)P AP T g(A) dA = ψ(I n ), and we conclude with corollary B.1 that ψ(I n ) = αI n with: α = 2 n I n · ψ(I n ) = 1 2n SOn(R) Tr(A) 2 g(A) dA.
3. Since ψ is a self adjoint operator, the orthogonal subspace Span(I n ) ⊥ is also a stable subspace. Moreover, using the change of variable A → A T , we see that ψ(J T ) = ψ(J) T and we have the decomposition: where S 0 n (R) and A n (R) are respectively the subspace of trace free symmetric matrices and the subspace of skew-symmetric matrices. They are both stable subspaces. 4. We prove now that ψ : S 0 n (R) → S 0 n (R) is a uniform scaling. By the spectral theorem, every matrix J ∈ S 0 n (R) can be written where P ∈ SO n (R) and D is diagonal. Since ψ(P DP T ) = P ψ(D)P T , it is enough to prove that there exists λ ∈ R such that for all diagonal matrices D, ψ(D) = λD.

5.
We prove similarly that ψ : A n (R) → A n (R) is a uniform scaling. Every J ∈ A n (R) can be written And there are of course many other ways to write the coefficients a, b and c.
Remark B.1. In dimension 4, the result still holds for symmetric matrices. For general matrices, the result can be proved in particular cases, for instance when g is a function of the trace, by using an explicit parametrisation of SO 4 (R) such as the 4-dimensional version of (10).