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 R3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^3$$\end{document} as in Degond et al. (Math Models Methods Appl Sci 27(6):1005–1049, 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 with the stable equilibria in the spirit of Degond et al. (Arch Ration Mech Anal 216(1):63–115, 2015, Math Models Methods Appl Sci 27(6):1005–1049, 2017).


Introduction
The model studied in the present work is a new elaboration of the work initiated in Degond et al. (2017) to model collective behaviour of agents described by their position and body attitude. New results about emergence of phenomena of bodyattitude 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 Communicated by Paul Newton.
Extended author information available on the last page of the article birds , fish school Hemelrijk and Hildenbrandt 2012) or sperm motion (Degond and Navoret 2015). 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 Degond et al. (2017), 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 Degond et al. (2018b) 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 (Cucker and Smale 2007;Ha and Liu 2009;Motsch and Tadmor 2011), attractive-repulsive models (Carrillo et al. 2017) or the Vicsek model for self-propelled particles (Vicsek et al. 1995). A recent review can be found in Albi et al. (2019). The present work belongs to the class of Vicsek-inspired models. Such models have two main 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 Vicsek et al. 1995). Here the dynamics takes place on the Riemannian manifold SO 3 (R). In the statistical physics literature, such models belong to the field of (soft-)active matter physics and are widely used to model biological phenomena, from the microscopic motility of cells to animal flocking. A review can be found in Marchetti et al. (2013).
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 many-particle systems. At a microscopic scale, the motion of each particle is detailed (individual-based model, IBM) through ordinary differential equations (ODEs) 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 Degond (2004). In particular, the BGK equation (for Bhatnagar-Gross-Krook) was introduced in Bhatnagar et al. (1954) 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 Perthame (1989) and Saint-Raymond (2003). The BGK operator has been used in a model of collective dynamics of self-propelled particles in Dimarco and Motsch (2016). However, together with Degond et al. (2018a), 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 Lanford 1975;Gallagher et al. 2014), from systems of interacting particles (mean-field limit and propagation of chaos Jabin 2014; Hauray and Jabin 2007;Sznitman 1991) or from stochastic processes (and in particular jump processes Mischler and Mouhot 2013;Kac 1956). Some of these techniques have been adapted to problems arising in the study of collective behaviour (Chuang et al. 2007;Bolley et al. 2011Bolley et al. , 2012. 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 Cercignani et al. (2013), Degond (2004) for a review) and is still an active research field (Guo and Jang 2010;Esposito et al. 2018;Caflisch 1980;Golse and Saint-Raymond 2004). 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 Degond and Motsch (2008), 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 (Jiang et al. 2016) and has already been successfully applied to a wide range of problems Zhang and Jiang 2017). It will be the key here to derive the macroscopic model in Sect. 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 Degond et al. (2017).
The aim of this work is to show the emergence of collective behaviour and selforganisation which give rise to macroscopic scale patterns such as clusters and travelling bands. 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 (Degond and Motsch 2008) 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 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 (Degond and Motsch 2008) has been the starting point of many other models of self-organised dynamics, including Degond et al. (2017) 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 Sect. 3. In the present work, we focus on phase transition phenomena between nonorganised 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: are the respective local density and flux. The measure on SO 3 (R) is the normalised Haar measure, the main properties of which are summarised in Sect. 3. 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 Sect. 4 and Eqs. (19) and (20)]. 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 Sect. 4 (see in particular Theorem 5 and Corollary 4.2). The second point is detailed in Sect. 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: The phase transition problem has been completely treated in the space-homogeneous case for the Vicsek model in , 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 Sect. 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 (Han et al. 2015;Wang and Hoffman 2008;Zhou and Wang 2011;Ball and Majumdar 2010;Ball 2017). These two worlds will be formally linked through the isomorphism between SO 3 (R) and the group of unit quaternions detailed in Sect. 4.2 and Appendix A. It will lead to the first point of Theorem 1 (the complete description of the equilibria, Sect. 4). As in Wang and Hoffman (2008), 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 Degond et al. (2017), Degond et al. (2018b). 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 is studied in Sect. 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 bodyattitude coordination (Degond et al. 2017(Degond et al. , 2018b) 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 Degond et al. (2017) for a model where only equilibria around a mean body-orientation were considered.
Finally, the SOHB model (Sect. 6) will be obtained as in Degond et al. (2017) by using the GCI. However, compared to Degond et al. (2017), 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 Degond et al. (2017). The SOHB model (43) 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 (Degond et al. 2020) and has been shown whenc 3 is constant.
The organisation of the work is the following: in Sect. 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 Sect. 3, we gather the main technical results we will constantly use throughout this work. In Sect. 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 (Wang and Hoffman 2008;Zhou and Wang 2011). In Sect. 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 Sect. 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.
. . , 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 i j its (i, j) coefficient. (1) • R + := [0, +∞), R * + := (0, +∞)

The BGK Equation and Other Related Models of Self-organisation
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 Diez (2019) 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 bodyattitude 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 Degond et al. (2018a), Dimarco and Motsch (2016): 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 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 , we can formally derive the equation satisfied by f as in (Degond et al. 2018a, Section 4.2). The rigorous formulation of this so-called propagation of chaos property has been investigated in Diez (2019).
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 : In the following we will focus on a purely local version of the BGK equation obtained in the limit K → δ 0 . This limit can be made rigorous starting from the IBM and by using an appropriate rescaling of the size of the kernel (called moderate interaction) with respect to the number of particles (see Oelschläger 1985;Jourdain andMéléard 1998 for classical results andDiez (2019) for additional details regarding the present model).
In the previous works, the IBM was typically given as in Degond et al. (2017) 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 Sect. 3). In this case, the resulting equation when N → +∞ is a nonlinear Fokker-Planck equation (see Bolley et al. 2012 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 The following theorem describes analogously the limiting behaviour of the laws of the particles as N → +∞.
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 S O 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 Degond and Motsch (2008) as a continuum version of the Vicsek model (Vicsek et al. (1995)). 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 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 Sect. 3. Apart from the fact that the underlying interaction process (Dimarco and Motsch 2016;Degond et al. 2018a) 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  and we believe that the same kind of study can be made in the body-attitude coordination dynamics modelled by a Fokker-Planck equation (4). 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 Figalli et al. (2018). However, the BGK model has another underlying gradient-flow dynamics (studied in Sect. 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 (3) 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 Model [19,25,27] [6, 13,14] [23] In progress [15,17] In progress

Present work
Vicsek Fokker-Planck  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 Degond et al. (2017): This terminology comes from the continuum version of the Vicsek model (Degond and Motsch 2008) where either the total flux or its normalisation is considered. A mathematical analysis of the normalised Vicsek model can be found in Figalli et al. (2018), Gamba and Kang (2016). The importance of this distinction in the context of phase transitions has been shown in  and Degond et al. (2013): phase transitions appear only in non-normalised models.
The following chart ( Fig. 1) shows the different models and gives references where they are studied (when such references exist).
Finally, in Sects. 4 and 5, we will focus on the spatially homogeneous version of the BGK model (3) 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 ρ f 0 ∈ R + associated with 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 (6) 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 (6) by A and integrating over SO 3 (R). Note that it contrasts with the Fokker-Planck case where even the well-posedness of the spatiallyhomogeneous equation would require further investigations. This will be part of future work.

Preliminaries: Structure and Calculus in SO n (R)
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, and we will mainly consider the case n = 3 in the next sections.

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, S O 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 S O n (R) at A is denoted by T A and M ∈ T A if and only if there exists P ∈
The norm on M n (R) associated with the inner product (7) will be denoted by · .
Remark 3.1 The unusual scaling factor 1/2 in the definition of the inner product (7) is chosen so that the induced Riemannian distance between a rotation matrix and the identity matrix is exactly equal to the angle of the rotation. Note that without this scaling factor, the matrix-norm induced by (7) is the classical Frobenius norm.
The general theory of locally compact topological groups ensures the existence of a left-invariant Haar measure μ on SO n (R). Since SO n (R) is compact thus unimodular, the Haar measure is both left-and right-invariant, which means that it satisfies for all P ∈ SO n (R) and all Borel set E of the Borel σ -algebra of SO n (R): where PE = {P A, A ∈ E} and E P = {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 i j ∈ 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 i j A multiplies the rows i and j by −1. Everything else remains unchanged.
• A = AD i j multiplies the columns i and j by −1. Everything else remains unchanged.
Everything else remains unchanged • A = P i j A multiplies row i by −1 and permutes the rows i and j.
• A = P i j A(P i j ) T exchanges the diagonal coefficients (i, i) and ( j, j) (and involves other changes).
The two following lemmas are important applications of these results.

Lemma 3.2 Let D ∈ M n (R) be a diagonal matrix and M D the von Mises distribution with parameter D, then
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
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) (R) (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 , where ψ ∈ [0, π] and ϕ ∈ [0, 2π ]. The volume form for the sphere is given by: • We have the following one-to-one map: : where and for p = (sin φ 1 sin φ 2 , cos φ 1 sin φ 2 , cos φ 2 ) T in spherical coordinates φ 1 ∈ [0, 2π ] 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 (Quarteroni et al. 2010, Section 1.9). 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 without loss of generality that Q ∈ SO 3 (R). Then, we can take: Remark 3.2 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 Fig. 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 Sect. 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 Sect. 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 (12). Before stating and proving it we will need the following lemma which is the analog of Lemma 4.4 in Degond et al. (2017). The proof of this lemma is an application of the results presented in Sect. 3.

Lemma 4.1 (Consistency relations)
The following holds: 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: c 2 (α) = [cos φ] α , where [·] α denotes the mean with respect to the probability density (16)

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 I 3 , i.e. there exists c 1 = c 1 (α) ∈ R such that The parametrisation of SO 3 (R) using Rodrigues' formula (8) then gives the explicit expression of c 1 by taking the trace in Eq. (17) 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 (15) holds. The parametrisation of SO 3 (R) coming from the isomorphism (9) then gives the explicit expression of c 2 by taking B = diag(1, 0, 0) in Eq. (15). First, using the change of variable A → P 13 A(P 13 ) T it holds that Then, using the parametrisation (10), it follows that:

Remark 4.2
We could alternatively use one of the two parametrisations of SO 3 (R) given in Sect. 3.2 or the quaternion formulation to prove that A α I 3 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: The solutions of the compatibility equation (18) 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 nonzero solution for the scalar compatibility equations (19) and (20) is not guaranteed for all values of ρ > 0. The existence of nonzero solutions for these equations will be explored in Sect. 4.2. They will determine the existence of equilibria for Eq. (6) for a given value of ρ (Corollary 4.2).

Remark 4.4
The fact that these matrices are solutions of the matrix compatibility equation (18) follows directly from the consistency relations (13) and (15) 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 (18) is necessarily of one of the forms listed in Theorem 5.
Remark 4.5 (Physical interpretation of the equilibria). The first two types of equilibria can be interpreted as statistical descriptions of, respectively, a disordered state (case J = 0) or an ordered (or flocking) state where a physical average body-orientation ∈ SO 3 (R) can be identified and where α plays the role of a concentration parameter around the "mean" value . Note that α depends on ρ. In the following, it will be shown that the flocking equilibrium is stable when the density ρ is sufficiently large. Moreover, the concentration parameter α will be larger for larger values of the density (after a certain threshold of the density). The third type of equilibria is specific to the body-orientation model and can be physically understood as the statistical description of a system composed of two groups of agents moving in the same direction but where one group is oriented upside-down with respect to the other. The diagonal element diag(1, 0, 0) ∈ B can indeed be obtained as the arithmetic average of the bodyorientations diag(1, 1, 1) ∈ SO 3 (R) and diag(1, −1, −1) ∈ SO 3 (R). It corresponds to the case of two agents moving in the e 1 direction and where the body orientation of an agent can be obtained from the other by a rotation of angle π around the e 1 axis. As it can be expected, equilibria of this type will be shown to be always unstable.
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 (18) can be reduced to a compatibility equation on diagonal matrices [equation (21)]. The second one (Proposition 4.2) provides a necessary condition for a diagonal matrix to be a solution of (21). The proof of Proposition 4.2 is deferred to the next section. 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: We will therefore consider only the following problem in dimension 3: find all the diagonal matrices D ∈ M 3 (R) such that where the set D is the subset of diagonal matrices which are the diagonal part of a SSVD and is defined by (11). Notice that Eq. (21) is just Eq. (18) restricted to the set D. (18) and (21) are equivalent. Notice that there might be other diagonal matrices in Orb(J ) [take, for example, J diagonal which does not satisfy the conditions (11)]. 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 (18) (i.e. the solutions of (21) 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 (18) of the form diag(0, −α, 0) where α > 0. The diagonal part of their SSVD is diag(α, 0, 0) and is a solution of (21).

Remark 4.7
A diagonal solution D of the matrix compatibility equation (18) 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 (Horn 1954, 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 δ I 3 and similarly for the other vertices of T , we see that the four vertices of Horn's tetrahedron belong to diag −1 ( ). Since is convex, we conclude that T ⊂ diag −1 ( ).
The diagonal solutions of the matrix compatibility equation (18) satisfy the following necessary condition. If α ∈ (0, +∞), the diagonal part of the SSVD of these diagonal matrices is equal to D = α I 3 .

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 (18). 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 (18) if and only if their parameter α ∈ R is solution of the scalar compatibility equation (20). Analogously, one can check that the diagonal matrices of type (b) are solutions of the matrix compatibility equation (18) if and only if their parameters α ∈ R are solutions of the scalar compatibility equation (19). This yields all the diagonal solutions of (18). Now, the solutions of (18) are exactly the matrices J ∈ Orb(D) where D is a diagonal solution of (18) and the set Orb Remark 4.8 When applied to diagonal matrices, the last part of Theorem 5 states that the diagonal solutions of (18) are necessarily of one of the types (a), (b) or (c) defined in Proposition 4.2 and that, it holds that 1. the matrix 0 is always a solution of (18) (20).

Proof of Proposition 4.2
The proof of Proposition 4.2 is based on two results. The first one has been proved in (Wang and Hoffman 2008, Section 4) to study the nematic alignment of polymers in higher-dimensional spaces: Theorem 6 (Wang and Hoffman 2008). Let n ≥ 3, b ∈ R + and s = (s 1 , s 2 , . . . , s n ) ∈ R n a solution of the nonlinear system 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 . Then, Card{s 1 , s 2 , . . . , s n } ≤ 2.
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 (21) into the compatibility equation (22) 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 S O 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 Eq. (7), and the second one is the usual dot product in R 4 . 3. For all q ∈ H/ ± 1, it holds that φ (q) = q ⊗ q − 1 4 I 4 . 4. The isomorphism φ preserves the diagonal structure: 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 (21) then becomes: Applying the isomorphism φ defined in Proposition 4.3 to this last equation, we can exchange φ and the integral by linearity and 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 (23). Thanks to Theorem 6, we conclude that if s is a solution of (24), 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 = s 2 = s 3 = s 4 = 0, then • 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 and similarly by permuting the diagonal elements when s 3 = 3α/4 and when the other elements are equal s 1 = s 2 = s 4 = −α/4 or when s 4 = 3α/4 and s 1 = s 2 = s 3 = −α/4, • if s 1 = s 2 = α/4 and s 3 = s 4 = −α/4 for α ∈ R, then and similarly by permuting the diagonal elements when s 1 = s 3 = α/4 and when s 2 = s 4 = −α/4 or s 1 = s 4 = α/4 and s 2 = s 3 = −α/4. The computation of the SSVD for these matrices is an easy computation. This concludes the proof of Proposition 4.2.

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 Eqs. (19) and (20) 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 Eqs. (19) and (20) in the plane (ρ, α) (see Fig. 2). We first prove the following proposition. (i) The function α → α/c 1 (α) is well defined on R, and 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 . (ii) The function α → α/c 2 (α) is even. It is decreasing on (−∞, 0), increasing on (0, ∞) and its value at zero is ρ c . (iii) We have the following asymptotic behaviours: Proof The idea of the proof is taken from Wang and Hoffman (2008).
(i) Since d dθ sin 2 (θ/2) sin θ = sin 2 (θ/2)(1 + 2 cos θ), an integration by parts shows that: α c 1 (α) = 3 π 0 sin 2 (θ/2)e α cos θ dθ π 0 sin 2 (θ/2) sin 2 θ e α cos θ dθ 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 (14). 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. 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 Fig. 2: Corollary 4.2 (Equilibria of the BGK operator, depending on the density ρ) The set of equilibria of the BGK operator (12) 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 Eqs. (19) and (20) and therefore the only equilibrium is the uniform equilibrium f eq = ρ. • For ρ = ρ * , in addition to the uniform equilibrium, there is a family of anisotropic equilibria given by f eq = ρ * M α * where ∈ SO 3 (R) and α * = ρ * c 1 (α * ). • For ρ * < ρ < ρ c , the compatibility equation (19) has two solutions α + and α − with 0 < α − < α + which give, in addition to the uniform equilibrium, two families of anisotropic equilibria: Moreover, Eq. (20) has two solutions −α 2 < 0 < α 2 which give another family of equilibria: f eq = ρ M α 2 B where B ∈ B. The uniform equilibrium is always an equilibrium.
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  in the Vicsek case. However, in SO 3 (R), there exist other equilibria which are not of this form. We will see in Sect. 5 that these latter equilibria are always unstable.
Finally, the following picture (Fig. 3) is a representation in the space R 3 of the diagonal parts of the SSVDs of the solution of the matrix compatibility equation (18) when ρ > ρ c . They all belong to the domain D defined by (11) and depicted in orange in Fig. 3.

Convergence to Equilibria
Now that we know all the equilibria of the spatially homogeneous BGK equation (6); 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 Eq. (6):  (12). 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 Sect. 5). The signs are the signature of the Hessian matrix Hess V (D) defined in Sect. 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 (Color figure online) 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 (6) by A ∈ SO 3 (R) and integrating over SO 3 (R).
is locally Lipschitz, the flow of Eq. (27) is defined globally in time since the map J → ρ A M J is Lipschitz with bounded Lipschitz seminorm. Notice that the solutions of the compatibility equation (18) are exactly the equilibria of the dynamical system (27). We therefore obtain the following proposition:  We will call stable/unstable an equilibrium of the BGK operator (12) such that the associated matrix J ∈ M 3 (R) is a stable/unstable equilibrium of the ODE (27). 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 (6) and let J f 0 = P D 0 Q be a SSVD. Let f (t) be the solution at time t ∈ R + of the spatially homogeneous BGK equation (6) with initial condition f 0 . Let D(t) be the solution at time t ∈ R + of the ODE (27) 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, writing N ρ := diag(N ρ ) ⊂ M 3 (R): converges as t → ∞ towards an equilibrium f eq of the BGK operator (12) 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 f 0 − 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 Sect. 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 , and gives rise to an hysteresis phenomenon.
The proof of this theorem can be found at the end of Sect. 5.2. It is based on a gradient flow structure for the flux J f studied in Sect. 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 Sect. 5.2. Additional details about the subset N ρ as well as the study of the critical cases ρ = ρ * and ρ = ρ c are provided in Sect. 5.3.

A Gradient-Flow Structure in (R) (R) (R) 3
In this section we show that the ODE (27) can be reduced to a gradient-flow ODE in R 3 . We first show (Sect. 5.1.1) how (27) can be reduced to an ODE in R 3 , the equilibria of which are linked to the equilibria of (27) [and therefore of (6)]. 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 (6) (Sect. 5.1.2).

Reduction to a Nonlinear ODE in R 3 and Equilibria
The ODE (27) 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.
Proposition 5.2 (Reduction to a nonlinear ODE in dimension 3) Let J 0 ∈ M 3 (R) be a given matrix and let D 0 ∈ Orb(J 0 ) be diagonal. Let P, Q ∈ SO 3 (R) such that

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 (27), then for any P, Q ∈ SO 3 (R), P J (t)Q is also a solution (and conversely). (ii) Since D 0 is diagonal, the fact that D(t) is also diagonal is a consequence of Lemma 3.2 which states that A M D is diagonal when D is diagonal.
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: is the solution of (28b) if and only if D(t) = diag( D(t)) is the solution of (28a). 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 (28a) is posed.

Proposition 5.3 (Invariant manifolds).
The following subsets of R 3 are invariant manifolds of the dynamical system (28b):

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 (11) 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 (28b) described in Proposition 5.3 are invariant manifolds of the dynamical system (28a). 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 (28a) posed on the domain D (see Fig. 3).
The following Proposition describes the equilibria of the dynamical system (28a) and is a consequence of the results of Sect. 4. Proposition 5.4 (Equilibria of the dynamical system (28a)). The equilibria of the dynamical system (28a) are, depending on the density ρ: • the matrix D = 0 for any density ρ, • the diagonal matrices of type (b) with parameters α + and α − when ρ * < ρ < ρ c , • the diagonal matrices of type (b) with parameters α 1 and α 3 and the diagonal matrices of type (c) with parameters ±α 2 when ρ > ρ c , where the types (b) and (c) are defined in Proposition 4.2 and the elements α + , α − , α 1 , α 2 and α 3 are defined in Corollary 4.2.
Proof The equilibria of (28a) are the diagonal matrices D such that which is exactly Eq. (18) for the diagonal matrices. This equation has been solved in Theorem 5, Remark 4.8 and Corollary 4.2.

Remark 5.3
As in the previous section, we have found all the diagonal equilibria of (28a). 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.
Proposition 5.5 (Equilibria of the BGK operator, equilibria of the ODE (28a)). Let J ∈ M 3 (R) with a SSVD given by J = P D Q. The following assertions are equivalent.
(1) The distribution f eq = ρ M J is an equilibrium of the BGK operator (12).
(2) The matrix D is an equilibrium of the dynamical system (28a) 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, respectively, the norm on M 3 (R) induced by the inner product (7) and the Euclidean norm on R 3 . Then, we can rewrite equations (28a) and (28b) into a gradient flow structure as follows: where ∇ is the gradient operator in M 3 (R) endowed with the Riemannian structure (7) 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 (7).

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 , 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).
As an application of the gradient-flow structure, we can prove that all the trajectories of (27) and (31a) are bounded: since the potential V in (30a) is decreasing along the trajectories, we can write for all t > 0 Using the fact that SO 3 (R) is compact with volume one, by the mean-value theorem applied to (29) there existsĀ ∈ SO 3 (R) such that Z(J (t)) ≤ eĀ ·J (t) . Therefore, owing to the fact that A ≤ √ 3/2 for all A ∈ SO 3 (R), we obtain using the Cauchy-Schwarz inequality that log Z(J (t)) ≤ 3 2 J (t) . Then, by Young's inequality we get: Moreover, since the potential V is analytic, a classical result by Haraux (2012), Łojasiewicz (1984) implies that the solution to (31a) and (27) will converge towards an equilibrium. The same holds for the system (31b). When all the equilibria of the dynamical system (28a) are hyperbolic the convergence towards the stable equilibria is exponentially fast [see (Hirsch et al. 2012, Section 9.3) and the Hartman-Grobman theorem (Perko 2013, Section 2.8)]. The goal of the next section is to find which equilibria among the ones found in Sect. 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 (28b) is the image by the isomorphism diag −1 of the flow of (28a), an equilibria D eq of (28a) is stable (resp. unstable) if and only if diag −1 (D eq ) is a stable (resp. unstable) equilibria of (28b). The stability properties of the equilibria of (28b) are much simpler to study than the stability properties of the equilibria of (28a) since they are given by the signature of the Hessian matrix Hess V ( D eq ) ∈ S 3 (R) of the potential V given in Eq. (30b) (the linearisation of ODE (28b) 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 (31a), the Hessian of the potential (30a) 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 (28a), 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 (28a)) • For 0 < ρ < ρ * , the only equilibrium is D = 0. This equilibrium is stable. 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.
where R := P P ∈ O 3 (R) is the Hadamard product (or element-wise product) of P with itself (that is R is the permutation matrix obtained by removing all the minuses in P). Let H ∈ R 3 and H := diag( H ) ∈ D 3 (R). From the expression (32) of the Hessian matrix Hess V (D) ≡ Hess V ( D), it follows that: Using the left and right invariance of the Haar measure and (33), we obtain: and therefore, since diag −1 (L T H ) = L T H , it follows that Similarly, using (33), we obtain and since diag −1 (P T H P) = R T H , it follows that The conclusion follows from Sylvester's law of inertia (Lax 2007, Chapter 8 Theorem 1).
If D is an equilibrium of (28b), then diag( D) is of one type (a), (b) or (c) defined in Proposition 4.2. For the types (b) and (c), the diagonal matrix diag( D) is therefore in the orbit of, respectively, α I 3 (by the action of G D ) or α diag(1, 0, 0) (by the action of G P ) where α solves the associated compatibility equation. As a consequence of Lemma 5.2, one only needs to look at the signature of the Hessian matrix of V taken in one of these two representatives. We are now ready to prove Theorem 8.

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, and 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 V (α I 3 ) where α = ρc 1 (α) and for the equilibria of type (c) we will compute the signature of Hess V (α 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 i j = 0 for all (i, j). Moreover, by the change of variable A = D ik A where k = i, j, we have when i = j: a ii a j j = − a ii a j j = 0.
It proves that is diagonal. Then, with the changes of variables A = P i j A or A = AP i j it can be seen that all the 3 2 quantities a 2 i j are equal. Since their sum is equal to n = 3, we get that where ρ c = 6. In conclusion, the signature of Hess
And finally: We find similarly that f (x) ≤ 0 when x ≤ 0.
And therefore, we can deduce the sign of the eigenvalue: Case 3. Equilibria of type (c) D = α diag(1, 0, 0) For D = α diag(1, 0, 0) with α = ρc 2 (α), using the parametrisation (10), it holds that: where the first equality comes from the change of variable A → P 13 A(P 13 ) T . Similarly, 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 (10), one can see that: sin φ 2 (1 ± cos φ 2 ) 2 e α 2 cos φ 2 dφ 2 dφ 1 dθ so that: (25) The conclusion of the proof follows from the study of the roots of Eqs. (19) and (20) provided by Corollary 4.2 and depicted in Fig. 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 (5). 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 (5) 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 (28a) described in Theorem 8. This technique is similar to the one which was used in Zhou and Wang (2011) 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 , 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 (26), 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 (28a) where the initial condition is the diagonal part of a SSVD of J f 0 . Since this equation has a gradient-flow structure (31a), by analyticity of the potential and boundedness of the trajectories, Łojasiewicz inequality (Haraux 2012;Łojasiewicz 1984) implies 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, since all the equilibria are hyperbolic, the Hartman-Grobman theorem (Perko 2013, Section 2.8) ensures that the convergence is locally exponentially fast in the sense that there exist constants δ, λ, C > 0 such that if J f 0 − J eq ≤ δ then for all t > 0, Let f eq = ρ M J eq . It follows from Duhamel's formula (26) 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 (36) into (35) and using (34), the first point of Theorem 7 follows with constants K = C L/|1 − λ| and μ = min(1, λ) when λ = 1 and K = C L and μ = 1 − ε for any ε > 0 when λ = 1. The stability of the equilibria of the dynamical system (28a) 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 ρ = diag −1 (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 (28b) are hyperbolic and we can apply the stable manifold theorem (Perko 2013, Section 2.7) which states that for a given hyperbolic equilibrium D eq of the nonlinear equation (28b), 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 linearised system around D eq . We are now ready to describe the behaviour of the system (28b) 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 3 (R)) and that our goal is to describe N ρ = diag −1 (N ρ ) (in the R 3 -framework). We write D(t) = diag D(t) the solution of (28a) with initial condition D 0 = diag D 0 . Case 1. When ρ < ρ * When ρ < ρ * , there is only the uniform equilibrium: D(t) −→ t→+∞ 0 and N ρ = ∅.
Finally, the space R 3 is partitioned in two domains by N ρ : depending on which domain D 0 belongs to in R 3 \ N ρ , D(t) will converge towards either the uniform equilibrium or to α + (1, 1, 1) T (or towards α − (1, 1, 1) T if D 0 belongs to its stable manifold).
• 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. Finally, in the case ρ > ρ c , the subset N ρ ∩D 3 (R) is equal to diag {0}∪M 1 ∪M 2 .

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 (28b), there exists a centre 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 (Haragus and Iooss 2010, 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 centre manifold and D(t) is attracted exponentially fast towards 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 (3). 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 (3) 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 (37) over SO 3 (R) leads to the conservation law: where ρ ε ≡ ρ f ε and The macroscopic model then depends on the region considered.
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 Sect. 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 body-attitude 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 : we write f ε = ρ ε + ε f ε 1 (where f ε 1 is defined by this relation) and notice that: Inserting this in (37), 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 ρ ε = 2 A · 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 By multiplying (40) by e 1 , it follows from (41) and (42) that: which gives the result by inserting this in (38).

Remark 6.1
This analysis does not depend on the dimension. In SO n (R) the same formal result holds:

Self-organised Hydrodynamics in an Ordered Region
In the following, for a given density ρ ∈ R + , α(ρ) denotes the maximal non-negative root of α = ρc 1 (α). We are going to prove the following theorem.
Then, ρ and satisfy the following system of partial differential equations: wherec 2 ,c 3 , c 4 are functions of ρ to be defined later and δ and r are the "divergence" and "rotational" operators defined in Degond et al. (2017): where ∇ x × is the curl operator.
The first equation (43a) is the conservation law (38) and the goal is to obtain the equation (43b) for = (t, x). However, here and contrary to the classical gas dynamics, the total momentum is not conserved: 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 Degond and Motsch (2008) for the study of the hydrodynamic limit of the continuum Vicsek model. Its precise setting in the context of our bodyattitude model is detailed Sect. 6.2.1. The formal proof of Theorem 9 can be found in Sect. 6.2.2.

Generalised Collision Invariants
To obtain an equation on , the main tool are the Generalised Collisional Invariants (GCI) first introduced in Degond and Motsch (2008). 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 with J is defined as: The condition ψ ∈ C J is equivalent to: Therefore, following the ideas of the proof of (Degond et al. 2017, 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: Now for any P ∈ A 3 (R), denoting f ε ≡ ε = P D(J f ε ), we get by multiplying the Eq. (37) by ψ ε P : The right-hand side vanishes by Definition (44) 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 (45) and since it is true for all P ∈ A 3 (R), we obtain: We have: where α denotes the derivative of α(ρ) with respect to ρ and similarly for ∂ i (M α ).
With this we compute the term: where we have used that: and · (∂ t + Ae 1 · ∇ x ) = 0.
Most of the terms that appear in X are computed in Degond et al. (2017). Precisely, where X 1 , X 2 , X 3 and X 4 are computed in Degond et al. (2017): where the coefficients {sin 2 θ(1 + 4 cos θ)} α , and and the matrix are the same as in Degond et al. (2017). The matrix D x ( ) ∈ M 3 (R) is defined as the unique matrix such that for all w ∈ R 3 , and smooth functions : Degond et al. 2017, Section 4.5). Note that C 3 = C 2 since the noise and alignment parameters which were denoted by ν and d in Degond et al. (2017) have been taken equal to 1 here. Note also that these coefficients are functions of ρ (through α only). The term Y is an additional term which appears here due to the presence of the parameter α = α(ρ) which is a function of ρ. It depends also on the derivative α of α: 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: and With the change of variable A → T A, these terms become and they can be computed using the same techniques as in Degond et al. (2017) 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 (47) where is invariant by transposition and conjugation and J := R 1,ρ − R T 1,ρ is a skew-symmetric matrix. From (49) and (50) 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 Degond et al. (2017): Finally, we obtain: Putting all the terms together, we can conclude as in Degond et al. (2017). First we notice that: Therefore, we obtain by multiplying (46) 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 Degond et al. (2017) and the coefficient c 3 in Degond et al. (2017) (which is equal to 1) becomes:

Conclusion and Perspectives
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 Vicsektype 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.
This model and its analysis raise many open questions. At the kinetic level, research perspectives include the following: 1. It would be interesting to quantify the time that takes a solution to reach the basin of attraction of a stable equilibria. This question has been addressed for other types of collective-dynamics models, for instance in Morales and Poyato (2019). 2. 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 Giacomin et al. (2012b), Giacomin et al. (2012a). 3. The space-inhomogeneous case, not considered here, is much more difficult to tackle due to the transport operator which prevents the reduction to a gradient-flow (in finite dimension for the BGK operator and in infinite dimension for the Fokker-Planck operator). This situation is not specific to the body-orientation model presented here, and the question is largely open even in the Vicsek case (see the conclusion of . Numerical simulations of the individual-based model suggest that clusterisation phenomena appear: the space can be decomposed in large-density regions (the clusters) and low-density regions. Inside a cluster, flocking is observed which is consistent with the fact that it is the only stable equilibrium in large density regions. Between the clusters, a diffusive behaviour is observed, and again, it seems consistent with the analysis in the low-density (homogenous) case presented in this article. The dynamics of the clusters (which can be seen as a kind of free-boundary problem) remains an open problem.
At the macroscopic level, the mathematical and numerical analyses of the SOHB model formally derived here are still in progress. The factor ρ in front of the transport term of Eq. (43b) suggests that a singularity when ρ → 0 could prevent well-posedness in low-density regions. A detailed investigation of the domain of validity of the SOHB model (and the diffusive model) at the macroscopic scale would be desirable. In particular, it is not known whether there are cases where the dynamics creates "empty spaces" in finite time. The results at the kinetic level on the stability of the uniform and flocking equilibria could give a first indication. However, it is not clear how these quantitative results at the kinetic level can be transposed at the macroscopic level through the rescaling procedure.
The BGK model studied here is a step towards the full description of the models of collective behaviour depicted in Fig. 1. The tools and ideas that we have presented here may help to analyse other models of body-attitude coordination such as the nonnormalised Fokker-Planck model in SO 3 (R). Other models which take into account curvature control in addition to body-orientation, in the spirit of Degond and Motsch (2011), could also be considered.

A Quaternions and Rotations
This 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: Q = α 4 P diag(3, −1, −1, −1)P T = α P diag(1, 0, 0, 0)P T − α 4 I 4 , and the result follows by taking q equal to the first column of P.
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 (Degond et al. 2018b, Proposition A.3) and (Degond et al. 2017, Lemma 4.2).

The expression
is a quadratic form for q. We take Q the matrix associated with this quadratic form. For J = (J i j ) i, j , using the explicit form of A = (q) with q = x + yi + z j + tk we obtain: where Q is the normalised uniaxial Q-tensor: Finally, 1 2 R(θ, n) · R(θ , n ) = q · Qq and we obtain thanks to the previous point: φ R(θ, n) = Q, that is to say: if J ∈ SO 3 (R), then φ(J ) is a normalised uniaxial Q-tensor. 4. If D = diag(d 1 , d 2 , d 3 ), then using the explicit form of φ given in the second point: and if Q = diag(s 1 , s 2 , s 3 , s 4 ) with s 1 + s 2 + s 3 + s 4 = 0, then φ −1 (Q) = 2 ⎛ ⎝ s 1 + s 2 s 1 + s 3 s 1 + s 4 ⎞ ⎠ .
Proof First we prove that the diagonal matrices form a subset of Span(SO n (R)): it is enough to show that (the other diagonal matrices with only one nonzero coefficient can be obtained in a similar way). When n is odd: and both matrices in the sum are in SO n (R). When n ≥ 4 is even, 0 4 I n−4 = 1 2 I n + 1 2 −I 4 I n−4 ∈ Span(SO n (R)), and similarly: 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 S O 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 A.1. To compute λ, notice that for the matrix e i ⊗ e j : 1 2 SO n (R) (e i · Ae j ) 2 d A = λ.

Proof of
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 = SO n (R) (J · A)(K · A) g(A) d A = J · ψ(K ).
2. We prove that Span(I n ) is a stable subspace for ψ: for any P ∈ SO n (R) we have: and we conclude with Corollary B.1 that ψ(I n ) = α I n with: α = 2 n I n · ψ(I n ) = 1 2n SO n (R) Tr(A) 2 g(A) d A.
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 Using the fact that Tr(D) = 0 and that all the SO n (R) a ii a kk g(A) d A are equal for i = k (by using conjugation by the matrices P i j , i, j = k, see Definition 3.1), we obtain: 5. We prove similarly that ψ : A n (R) → A n (R) is a uniform scaling. Every J ∈ A n (R) can be written is a block diagonal matrix with blocks When n ≥ 3, by using conjugation by the matrices D 2 j−1,2 j for j ∈ {1, . . . , n/2 } (Definition 3.1) we can see that for each k ∈ {1, . . . , p}, the matrix is of the form (48). Moreover, when n ≥ 5, by using conjugation by matrices D 2 j−1,2 −1 for j = , j, ∈ {1, . . . , n/2 } and j, = k, we can see that all the diagonal blocks of M k are equal to zero except the one in position 2k − 1. When n = 3 there is only one block in position 1 so the result holds but when n = 4 such j and do not exist. In conclusion, when n = 4, ψ(C) is of the form (48) and each diagonal block C 2k−1 of ψ(C) is written C 2k−1 = μ 2k−1 C 2k−1 and γ = 1 2 (λ − μ) = 1 4 SO n (R) (a 2 11 − a 11 a 22 ) − (a 21 − a 12 )a 21 g(A) d A.
7. In conclusion, writing the decomposition 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 (9).