# Body-attitude alignment: first order phase transition, link with rodlike polymers through quaternions, and stability

Amic Frouvelle

► To cite this version:

Amic Frouvelle. Body-attitude alignment: first order phase transition, link with rodlike polymers through quaternions, and stability. Recent Advances in Kinetic Equations and Applications, 2021, 978-3-030-82945-2. 10.1007/978-3-030-82946-9_7. hal-03027574

HAL Id: hal-03027574

https://hal.science/hal-03027574

Submitted on 27 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Body-attitude alignment: first order phase
transition, link with rodlike polymers through
quaternions, and stability

Amic Frouvelle*†

Abstract

We present a simple model of alignment of a large number of rigid bod-
ies (modeled by rotation matrices) subject to internal rotational noise. The
numerical simulations exhibit a phenomenon of first order phase transition
with respect the alignment intensity, with abrupt transition at two thresh-
olds. Below the first threshold, the system is disordered in large time: the
rotation matrices are uniformly distributed. Above the second threshold, the
long time behaviour of the system is to concentrate around a given rotation
matrix. When the intensity is between the two thresholds, both situations
may occur.

We then study the mean-field limit of this model, as the number of particles
tends to infinity, which takes the form of a nonlinear Fokker-Planck equation.
We describe the complete classification of the steady states of this equation,
which fits with numerical experiments. This classification was obtained in
a previous work by Degond, Diez, Merino-Aceituno and the author, thanks
to the link between this model and a four-dimensional generalization of the
Doi-Onsager equation for suspensions of rodlike polymers interacting through
Maier-Saupe potential.

This previous study concerned a similar equation of BGK type for which
the steady-states were the same. We take advantage of the stability results
obtained in this framework, and are able to prove the exponential stability of
two families of steady-states: the disordered uniform distribution when the
intensity of alignment is less than the second threshold, and a family of non-
isotropic steady states (one for each possible rotation matrix, concentrated
around it), when the intensity is greater than the first threshold. We also
show that the other families of steady-states are unstable, in agreement with
the numerical observations.

<footnote>*CEREMADE, CNRS, Université Paris-Dauphine, Université PSL, 75016 Paris, France. E-mail: frouvelle@ceremade.dauphine.fr</footnote>

<footnote>†CNRS, Université de Poitiers, UMR 7348 — Laboratoire de Mathématiques et Applications (LMA), 86000 Poitiers, France. E-mail: amic.frouvelle@math.univ-poitiers.fr</footnote>
# Introduction

The mathematical study of active matter, such as aligning self-propelled particles, is now a well established field of research, inspired for instance by phase transition phenomena that appear in the Vicsek model [VCBJ+95, CGGR08]. Following the kinetic approach introduced in [DM08], a simple model of alignment of unit vectors subject to internal rotational noise gives rise to a continuous phase transition at the kinetic level [FL12]. When the alignment intensity (that we call $\rho$, since it is related to the local density $\rho$ of particles in the inhomogeneous version [DFL13], where the unit vectors represent the velocities of self-propelled particles) is below a threshold $\rho_c$, the only stable steady-state is the uniform distribution on the unit sphere. On the other hand, when $\rho > \rho_c$, this isotropic equilibria becomes unstable and a family of stable equilibria arises: von Mises distributions with concentration parameter depending on $\rho$, around a given unit vector. When setting the intensity of alignment as a nonlinear function of the order parameter of the system [DFL15], this continuous phase transition may become a discontinuous one (or first order), with hysteresis phenomenon: a second threshold $\rho^* < \rho_c$ appears, the uniform equilibrium distribution being stable for $\rho < \rho_c$ and the concentrated distributions being stable for $\rho > \rho^*$. Around those thresholds, the order parameter cannot vary continuously from a family of equilibria to the other.

Recently, in a work with Degond and Merino-Aceituno [DFMA17] we extended the model of self-propelled particles of Degond and Motsch [DM08] to the case where the orientation of particles are not only given by their velocity (a unit vector) but by their whole body attitude (an orthonormal frame, given by a rotation matrix). Then, still with Degond and Merino-Aceituno, together with Trescases [DFMAT18] we proposed a similar model based on quaternion representation for rotation matrices, and the models appeared to be equivalent. In these models, the interaction was normalized and no phenomenon of phase transition could occur, but we remarked that the non-normalized version may lead to such a phenomena. Finally, with Degond, Diez and Merino-Aceituno [DDFMA20] we managed to treat this phenomenon of phase transition in a homogeneous Bhatnagar-Gross-Krook (BGK) model, thanks to this link with unit quaternions and an analogy with a four-dimensional generalization of the Doi-Onsager equation for suspensions of rodlike polymers interacting through Maier-Saupe potential. Indeed, the compatibility equation we need to solve to determine the possible steady-states can be reformulated in this quaternionic formulation, and leads to a compatibility equation for the Maier-Saupe potential in dimension 4, which was solved in [WH08]. We obtain a discontinuous phase transition with two thresholds $\rho^* < \rho_c$, still with the same two types of stable equilibria: the uniform distribution for $\rho < \rho_c$, and a family of generalized von Mises distributions, concentrated around a given rotation matrix when $\rho > \rho^*$.

The aim of this paper is twofold. We first want to introduce the model of alignment of rigid bodies through numerical simulations of the particle system, in order to present the first order phase transition that we observe numerically. And then we want to provide a rigorous mathematical description of this phase transition phenomenon, at the kinetic level: the mean-field limit of the particle system when the number of particles is large is given by a nonlinear Fokker-Planck equation, for
which the steady states are the same as those characterized in [DDFMA20] for the BGK equation. The main result of this article is that we have a fine description of the long-time behaviour of the solution to the Fokker-Planck equation: we classify all the families of equilibria regarding their stability, and prove the exponential stability of the uniform equilibrium when $\rho < \rho_c$ and of the concentrated von Mises distributions when $\rho > \rho^*$.

In Section 1, we present the framework of our model: a system of coupled stochastic differential equations for $N$ matrices in $SO_3(\mathbb{R})$. We present a time discretization scheme of Euler-Naruyama type, and provide numerical simulations which illustrate the phenomenon of first order phase transition. In Section 2, we describe the mean-field limit of this system, which takes the form of a nonlinear Fokker-Planck equation. We give general results on the behaviour of the solution of this evolution equation, and we show that the determination of its steady states amounts to solve a matrix compatibility equation. Thanks to the free energy associated to the Fokker-Planck equation, the uniform equilibria is shown to be unstable for $\rho > \rho_c = 6$, proving that in that case there are others solutions than 0 for the compatibility equation. Section 3 is a summary of the results of [DDFMA20] to solve this compatibility equation: we present the link between rotation matrices and unit quaternions, and the fact that the compatibility equation can be transformed to a compatibility equation for $Q$-tensors which was solved in [WH08]. We therefore get a precise description of all the steady-states of the equation, and a way to obtain the second threshold $\rho^*$ (as the minimum of a one-dimensional function) such that for $\rho \ge \rho^*$ there exists non-trivial steady-states. In Section 4, we summarize the results of [DDFMA20] regarding the stability of these equilibria in the framework a BGK equation (which shares the same steady-states), and we are able to use these results to obtain the classification of the steady-states, as critical points of the free energy. In particular we show that three families of equilibria are unstable, and the remaining two other types are local minimizers of the free energy: the uniform distribution when $\rho < \rho_c$ and the concentrated von Mises distributions for $\rho > \rho^*$. Finally, Section 5 is devoted to the main new result of this paper: the exponential stability of these two types of steady-states. In Theorem 3, we prove that if a function $f_0$ is sufficiently close to the set of equilibria (in relative entropy), then there exist such an equilibrium $f_\infty$ such that the solution of the Fokker-Planck equation converges exponentially fast towards $f_\infty$ (still in relative entropy). We finish this last section by some comments and perspectives.

# 1 Numerical evidence of a first-order phase transition in a system of interacting particles

## 1.1 A simple SDE on $SO_3(\mathbb{R})$ and its time-discretization

First of all let us recall some basic facts about $SO_3(\mathbb{R})$.

**Definition 1.1.** For any $\mathbf{u} = \begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix} \in \mathbb{R}^3$ we denote by $[\mathbf{u}]_\times = \begin{pmatrix} 0 & -u_3 & u_2 \\ u_3 & 0 & -u_1 \\ -u_2 & u_1 & 0 \end{pmatrix}$
the (antisymmetric) matrix associated to the linear map $\mathbf{v} \in \mathbb{R}^3 \mapsto \mathbf{u} \times \mathbf{v}$ in the canonical basis.

**Proposition 1.1** (Rodrigues' formula). Any special orthogonal matrix $A \in SO_3(\mathbb{R})$ can be written as a rotation around an axis in $\mathbb{R}^3$. More precisely, there exists a unique angle $\theta \in [0, \pi]$ and a unit vector $\mathbf{n} \in S_2$ such that $A$ is the rotation $R(\theta, \mathbf{n})$ of angle $\theta$ around the axis directed by $\mathbf{n}$, given by the following formula:

$$R(\theta, \mathbf{n}) = \exp(\theta[\mathbf{n}]_{\times}) = \cos \theta I_3 + \sin \theta[\mathbf{n}]_{\times} + (1 - \cos \theta)\mathbf{n}\mathbf{n}^{\top}. \quad (1)$$

where $I_3$ is the identity matrix. When $\theta \in (0, \pi)$, the unit vector $\mathbf{n}$ is unique. When $\theta = \pi$ there are two such vectors $\mathbf{n}$, opposite one to the other. And when $\theta = 0$, any unit vector $\mathbf{n}$ can be used.

To introduce the model and some important notations, we first start with a simple stochastic differential equation (SDE) modeling a rotation matrix $A(t) \in SO_3(\mathbb{R})$ trying to align with another fixed rotation matrix $A_0 \in SO_3(\mathbb{R})$, with strength of alignment $\nu > 0$, and subject to angular noise of intensity $\tau > 0$:

$$dA = -\nu \nabla_A (\frac{1}{2} \|A - A_0\|^2) dt + 2\sqrt{\tau} P_{T_A} \circ dB_t. \quad (2)$$

To give a meaning to the previous equation, let us describe the terms one by one, from left to right. We need to define a metric on $SO_3(\mathbb{R})$ in order to define the gradient $\nabla_A$. As it is usually the case in $SO_3(\mathbb{R})$, we will take the metric induced by the scalar product in $M_3(\mathbb{R})$ given by

$$A \cdot B = \frac{1}{2} \mathrm{Tr}(AB^{\top}). \quad (3)$$

One of the reasons to take this metric is that the geodesic distance between a matrix $A \in SO_3(\mathbb{R})$ and its composition by a rotation matrix of angle $\theta \in [0, \pi]$ is exactly $\theta$. Said differently, if $\mathbf{n} \in S_2$, then the curve $\theta \in \mathbb{R} \mapsto R(\theta, \mathbf{n})A$ given by the formula (1) is a geodesic travelled at unit speed. The other reason is that the map $\mathbf{u} \in \mathbb{R}^3 \mapsto [\mathbf{u}]_{\times}$ given by Definition 1.1 is an isometry from $\mathbb{R}^3$ to the antisymmetric matrices (which is the Lie algebra of $SO_3(\mathbb{R})$). The norm $\|A - A_0\|$ in the SDE (2) is the one associated to this scalar product. The operator $P_{T_A}$ is the orthogonal projection on the tangent space of $SO_3(\mathbb{R})$ at $A$, given by $P_{T_A}H = \frac{1}{2}(H - AH^{\top}A)$. The notation $\circ$ in the SDE (2) means that it must be understood in the Stratonovich sense, and the Brownian motion $B_t$ is a $3 \times 3$ matrix whose entries are independent real standard Brownian motions<fnref n="1" />. This ensures that the matrix $A$ stays on $SO_3(\mathbb{R})$ for all time, and this is the usual way of defining SDEs on manifolds (we refer to [Hsu02] for a reference on this topic). Therefore the first term in the right-hand

<footnote n="1">¹Note that this does not give a standard Brownian motion on the euclidean space $M_3(\mathbb{R})$, equipped with this scalar product, but $\tilde{B}_t = \sqrt{2}B_t$ is such a standard Brownian motion. The SDE for a standard Brownian motion on the manifold, with generator $\frac{1}{2}\Delta_A$, would be $dA = P_{T_A} \circ d\tilde{B}_t$, which explain the choice of $2\sqrt{\tau}$ instead of the usual $\sqrt{2\tau}$ in the SDE (2) so that the Fokker-Planck equation (4) has the simplest coefficients.</footnote>
side of (2) may be written $\nu\nabla_A(A \cdot A_0)$ since $\|A\|^2 = \frac{3}{2}$ whenever $A \in SO_3(\mathbb{R})$. Finally, the law $t \mapsto \mu(t, \cdot)$ (with values in $\mathcal{P}(SO_3(\mathbb{R}))$, the set of probability measures on $SO_3(\mathbb{R})$) of such a process satisfies the following Fokker-Planck equation:

$$ \partial_t \mu + \nu \nabla_A \cdot (\nabla_A (A \cdot A_0) \mu) = \tau \Delta_A \mu, \quad (4) $$

where $\nabla_A \cdot$ and $\Delta_A$ are the divergence and Laplace-Beltrami operators on $SO_3(\mathbb{R})$. Up to a time rescaling, we see that the important parameter is $\kappa = \frac{\nu}{\tau}$, and we can then without loss of generality study the following PDE, obtained by replacing $\tau$ by 1 and $\nu$ by $\kappa$ in (4):

$$ \partial_t \mu = -\kappa \nabla_A \cdot (\nabla_A (A \cdot A_0) \mu) + \Delta_A \mu = \nabla_A \cdot \left[ \exp(\kappa A \cdot A_0) \nabla_A \left( \frac{\mu}{\exp(\kappa A \cdot A_0)} \right) \right]. \quad (5) $$

In view of the above formulation, we now define the generalized von Mises distribution (a probability measure) on $SO_3(\mathbb{R})$ of parameter $J \in M_3(\mathbb{R})$ by

$$ M_J(A) = \frac{1}{\mathcal{Z}(J)} \exp(J \cdot A), \text{ where } \mathcal{Z}(J) = \int_{SO_3(\mathbb{R})} \exp(J \cdot A) dA, \quad (6) $$

the normalized volume form on $SO_3(\mathbb{R})$ being its Haar probability measure (this comes from invariance of the metric with respect to left or right multiplication by a given rotation matrix). Therefore it is for instance easy to see that $\mathcal{Z}(\kappa A_0)$ only depends on $\kappa$ when $A_0 \in SO_3(\mathbb{R})$. With this notation, we can multiply the PDE (5) by $\frac{\mu}{M_{\kappa A_0}}$, integrate by parts and take advantage of the fact that the integral of $\mu$ on $SO_3(\mathbb{R})$ remains constant in time, to obtain

$$ \frac{1}{2} \frac{d}{dt} \int_{SO_3(\mathbb{R})} \left| \frac{\mu}{M_{\kappa A_0}} - 1 \right|^2 M_{\kappa A_0} dA = - \int_{SO_3(\mathbb{R})} \left\| \nabla_A \left( \frac{\mu}{M_{\kappa A_0}} - 1 \right) \right\|^2 M_{\kappa A_0} dA. \quad (7) $$

Together with a weighted Poincaré inequality on $SO_3(\mathbb{R})$, this shows that the solution to the PDE (5) converges exponentially fast to the von Mises distribution $M_{\kappa A_0}$. Let us remark that when $\kappa$ is small (strong noise, or weak alignment), this distribution tends to be uniform on $SO_3(\mathbb{R})$, and when $\kappa$ is large (strong alignment or low level of noise), it is concentrated around the maximizer of $A \mapsto A \cdot A_0$, which is exactly $A_0$, as expected.

Let us finish this subsection by describing a numerical discretization of the SDE (2). By using the fact that $\nabla_A(A \cdot A_0) = P_{T_A} A_0$, and denoting by $\Pi$ the orthogonal projection on $SO_3(\mathbb{R})$ (well-defined in a neighborhood of the manifold), a naive projected Euler-Naruyama scheme would read as follows:

$$ A(t + \Delta t) \approx \Pi(A(t) + \nu \Delta t P_{T_{A(t)}} A_0 + \sqrt{\Delta t} 2\sqrt{\tau} P_{T_{A(t)}} N_9), \quad (8) $$

where $N_9$ is a three by three matrix whose 9 entries are independent samples of standard Gaussian distribution. One could even remove the projections on the tangent plane and use this model, easy to describe as a starting point: “Start from $A \in SO_3(\mathbb{R})$, move with step $\nu \Delta t$ in the direction of the target $A_0$, add some noise of intensity $2\sqrt{\tau \Delta t}$ and project the result back on $SO_3(\mathbb{R})$”. However, there is a way to avoid sampling 9 entries per step and to take advantage of the Lie group
structure of $SO_3(\mathbb{R})$ instead of computing the projection on $SO_3(\mathbb{R})$ (which is the polar decomposition of matrices and may have some cost). Indeed, the right-hand side of the scheme (8) can be written

$$ \Pi(I_3 + \frac{1}{2}\nu\Delta t [A_0A(t)^{\top} - A(t)A_0^{\top}] + \sqrt{\tau}\Delta t [\mathcal{N}_9A(t)^{\top} - A(t)\mathcal{N}_9^{\top}])A(t). $$

Since a rotation of a standard Gaussian vector is still a standard Gaussian vector, one can see that the matrix $\mathcal{N}_9A(t)^{\top}$ is also a matrix whose 9 entries are independent samples of standard Gaussian distribution. Therefore $\mathcal{N}_9A(t)^{\top} - A(t)\mathcal{N}_9^{\top}$ is an antisymmetric matrix whose independent entries are samples of centered Gaussian distribution of variance 2. It is then a matrix of the form $\sqrt{2}[\eta]_{\times}$ (see Definition 1.1), where $\eta$ is a standard Gaussian vector in $\mathbb{R}^3$. When $H$ is a small antisymmetric matrix, a consistent approximation to $\Pi(I_3+H)$ is given by $\exp(H)$ and can be computed thanks to Rodrigues' formula (1). Therefore a numerical scheme consistent with the naive scheme (8) is given by

$$ A(t + \Delta t) \approx \exp(\frac{1}{2}\nu\Delta t [A_0A(t)^{\top} - A(t)A_0^{\top}] + \sqrt{2\tau}\Delta t [\eta]_{\times})A(t), \quad (9) $$

where $\eta$ is a standard Gaussian vector in $\mathbb{R}^3$.

## 1.2 A system of SDEs and its numerical simulations

We are now ready to introduce our model. In the article [DFMA17], we considered $N$ individuals located at positions $X_i \in \mathbb{R}^3$ for $1 \le i \le N$ and with body orientations $A_i \in SO_3(\mathbb{R})$, moving at unit speed in the direction of their first vector $A_i e_1$ and aligning their orientations with their neighbours, as in the simple SDE (2). This could take the following form<fnref n="2" />:

$$ \left\{ \begin{aligned} \mathrm{d}X_k &= A_k e_1 \mathrm{d}t \\ \mathrm{d}A_k &= -\sum_{j=1}^{N} \nu_{j,k} \nabla_{A_k} (\frac{1}{2}\|A_k - A_j\|^2) \mathrm{d}t + 2\sqrt{\tau} P_{T_{A_k}} \circ \mathrm{d}B_{t,k}, \end{aligned} \right. \quad (10) $$

where $\nu_{j,k}$ is the intensity at which particle $k$ aligns with particle $j$, and which may depend for instance on the distance $\|X_j - X_k\|$ between the particles. We consider here a much simpler model, homogeneous in space, so we only look at $N$ rotation matrices $(A_i)_{1 \le i \le n} \in SO_3(\mathbb{R})$, with the same intensity $\frac{\rho}{N}$ of alignment between any pair of particles. We are therefore interested in the following system of SDEs, using the fact that $\nabla_A(\frac{1}{2}\|A - A_0\|^2) = -\nabla_A(A \cdot A_0) = -P_{T_A}A_0$:

$$ \forall k \in 1 \dots N, \quad \mathrm{d}A_k = \frac{\rho}{N} \sum_{j=1}^{N} P_{T_{A_k}} A_j \mathrm{d}t + 2\sqrt{\tau} P_{T_{A_k}} \circ \mathrm{d}B_{t,k}. $$

In this model, when all the rotation matrices are close to a given one $A_0$, the behaviour of the system can be expected to be similar to the one of the simple

<footnote n="2">²Actually, the model studied in [DFMA17] (which does not present the phenomenon of phase transition we are studying here) is a little bit more involved: each particle first chose an average target and aligns with it, instead of averaging the “forces of alignment” as it is the case in the system of SDEs (10).</footnote>
SDE (2), and we may expect the matrices to concentrate if the alignment intensity $\rho$ is high (or $\tau$ is low). Conversely, if they are not concentrated around some target, the average of the alignment forces is small and the noise level may prevent the matrices to align if $\rho$ is low (or $\tau$ is high). From now on, up to rescaling time (and dividing $\rho$ by $\tau$), we consider the case $\tau = 1$ and we denote by $J$ the average “flux”, so our system has the following form:

$$ \left\{ \begin{aligned} dA_k &= P_{T_{A_k}} J dt + 2P_{T_{A_k}} \circ dB_{t,k}, \quad (1 \le k \le N) \\ J(t) &= \frac{\rho}{N} \sum_{j=1}^{N} A_j(t). \end{aligned} \right. \qquad (11) $$

We are then interested in the different behaviours of the system (11) for different values of $\rho$. One way to measure how much matrices are concentrated is to compute the variance $\langle \|A - \langle A \rangle\|^2 \rangle$ (where we denote $\langle h(A) \rangle = \frac{1}{N} \sum_{j=1}^N h(A_j)$ for any function $h$). This nonnegative quantity is equal to $\langle \|A\|^2 \rangle - \|\langle A \rangle\|^2 = \frac{3}{2} - \|\frac{J}{\rho}\|^2$, which implies that if we define the order parameter $c(t)$ by

$$ c(t) = \frac{\sqrt{2}}{\sqrt{3\rho}} \|J(t)\|, \qquad (12) $$

we obtain a quantity between 0 (when the variance is maximal) and 1 (the variance is 0, all matrices are the same). To give a numerical illustration of the phenomenon we are interested in, we use a scheme similar to the scheme (9) of the previous subsection: we take $N$ matrices $A_k \in SO_3(\mathbb{R})$ for $1 \le k \le N$), a time step $\Delta t$, and at each time iteration, we compute $J = \frac{\rho}{N} \sum_{j=1}^N A_j$ and we update each $A_k$ for $1 \le k \le N$ with the matrix

$$ \exp(\frac{1}{2}\Delta t [JA_k^\top - A_k J^\top] + \sqrt{2\Delta t}[\eta_k]_\times)A_k, $$

where $(\eta_k)_{1 \le k \le N}$ are independent samples of a standard Gaussian vector in $\mathbb{R}^3$.

Figure 1 depicts the time evolution of the order parameter $c(t)$ given by the formula (12) for two realisations of this numerical scheme. In both cases the number of particles is $N = 500$, the time step is $\Delta t = 0.04$ and we run the simulation for 100 time iterations. In the top-left part of Figure 1 where we took $\rho = 1$, even if we started with all the particles in the same position (order parameter equal to 1), as time evolves, the order parameter becomes very small. In the top-right part, with $\rho = 10$, even if the particles were uniformly sampled on $SO_3(\mathbb{R})$ (order parameter close to 0), as time evolves, the order parameter stabilizes around a quite high value, indicating that the matrices are concentrated around a given rotation matrix. This indicates that a phase transition phenomenon is occurring with respect to the parameter $\rho$. However, for some intermediate values of $\rho$, as in the bottom part of Figure 1 where $\rho = 5$, two different behaviours may happen: starting with concentrated particles lead to an order parameter stabilizing around a non-zero value, while the configuration starting with particles uniformly sampled on $SO_3(\mathbb{R})$ stays with an order parameter close to 0 as time evolves.

In order to obtain a more precise illustration of this phenomenon, we ran 500 such
<figure><img src="image_2.png" /><figcaption>Figure 1: Time evolution of the order parameter in four situations.</figcaption></figure>

simulations with various values of the parameter $\rho$ and different initial conditions<fnref n="3" />, still with $N = 500$ and $\Delta t = 0.04$, for 500 time iterations. Figure 2 depicts the initial order parameters $c$ and strengths $\rho$, and their value after 500 iterations ($t = 20$). We clearly see two thresholds for $\rho$. The first threshold that we will denote $\rho^*$, is such that for all simulations with $\rho < \rho^*$, the order parameter seems to be close to 0 for large times. The second threshold, that we will denote $\rho_c$ (with $\rho^* < \rho_c$), is such that for all simulations with $\rho > \rho_c$, the order parameter does not stay close to 0 for large times, and stabilizes around a quite high value. In the intermediate regime $\rho^* < \rho < \rho_c$, both behaviours occur. This is what is called first-order (or discontinuous) phase transition: the order parameter does not vary continuously when going from one behaviour to the other.

The aim of the next sections is to present a rigorous mathematical description of this phenomenon in the framework of a kinetic equation corresponding to the limiting behaviour of the system of SDEs (11) when $N \to \infty$.

<footnote n="3">³For a better illustration, the parameter $\rho$ and the initial order parameter $c$ are not uniformly sampled, in order to see more points in the region of interest.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Numerical illustration of a first-order phase transition.</figcaption></figure>

## 2 Mean-field limit and compatibility equation

Let us first consider the first part of the system (11), as if $t \in \mathbb{R} \mapsto J(t) \in M_3(\mathbb{R})$ was a prescribed regular function:

$$dA = P_{T_A} J dt + 2P_{T_A} \circ dB_t. \qquad (13)$$

As before for the simple SDE (2), the law $t \mapsto \mu(t, \cdot)$ of such a stochastic process would satisfy the following (linear) Fokker-Planck equation:

$$\partial_t \mu = -\nabla_A \cdot (\mu P_{T_A} J) + \Delta_A \mu = \nabla_A \cdot \left[ M_J(A) \nabla_A \left( \frac{\mu}{M_J(A)} \right) \right], \quad (14)$$

where the definition of the generalized von Mises distribution $M_J$ is given by the formula (6). Let us now suppose that several such processes $A_k$ satisfying the SDE (13) were independently drawn, with different independent Brownian motions $B_{t,k}$, and independent initial conditions following a probability measure $\mu_0$ on $SO_3$. Their law at time $t$ would be given by $\mu(t, \cdot)$, solution of the Fokker-Planck equation (14) with initial condition $\mu_0$ by the law of large numbers the average $\frac{1}{N} \sum_{k=1}^N A_k(t)$ would converge to the expectation of one of this process, that we call $\mathcal{J}[\mu(t, \cdot)]$. More generally, we define $\mathcal{J}[f]$ for any finite measure $f$ on $SO_3(\mathbb{R})$ (not necessarily a probability measure, it may also be a signed measure):

$$\mathcal{J}[f] = \int_{SO_3(\mathbb{R})} A f(A) dA. \qquad (15)$$

To deal with the system (11), where $J(t) = \frac{\rho}{N} \sum_{k=1}^N A_k(t)$ is not prescribed but depends on all the particles, we cannot expect the particles $A_k$ to behave independently. However one can show that in the limit $N \to \infty$, their behaviour is close to independent particles. This is called the propagation of chaos property, and we refer to [Szn91] for an introduction on this subject. One of the typical results in this theory is that the empirical measure of the particle system converges to a solution to the (now nonlinear) Fokker-Planck equation corresponding to (14) with $J(t) = \rho \mathcal{J}[\mu(t, \cdot)]$:
**Proposition 2.1.** If $A_{k,0}$ are independent random rotation matrices distributed according to the probability measure $\mu_0$, then the empirical measure $\mu^N(t) = \frac{1}{N} \sum_{k=1}^N \delta_{A_k(t)}$ associated to the solution of the system of SDEs (11) converges (in Wasserstein distance) to the solution $\mu$ of the following nonlinear Fokker-Planck equation, with initial condition $\mu_0$:

$$ \partial_t \mu = -\rho \nabla_A \cdot (\mu P_{T_A} \mathcal{J}[\mu]) + \Delta_A \mu. \quad (16) $$

The convergence is uniform on $[0, T]$ for all $T > 0$.

*Proof.* We will not provide the proof in detail here, as it follows the classical theory of propagation of chaos for coupled drift-diffusion processes, but we will recall some important steps. It has to be adapted to the framework of SDEs on a manifold, but this is not a real problem in this compact case (see for instance [BCC12] in the case of the Vicsek model on the sphere). Let us recall the coupling argument such as the one in [Szn91]. We start by proving the well-posedness of this following SDE (the coupling process):

$$ \begin{cases} dA = \rho P_{T_A} \mathcal{J}[\mu]dt + 2P_{T_A} \circ dB_t, \\ \mu(t, \cdot) \text{ is the law of } A(t). \end{cases} \quad (17) $$

The proof of this well-posedness, seen as a fixed point problem (either for the function $J(t) = \rho \mathcal{J}[\mu]$ or directly on the law $\mu$) is done thanks to a Picard iteration which leads to a contraction in the appropriate Wasserstein metric.

We then construct independent solutions to this coupling process $\bar{A}_k$ with independent Brownian motions $B_{t,k}$ and initial conditions $A_{k,0}$: the same as the Brownian motions and initial conditions used for the original system of SDEs (11). All these processes $\bar{A}_k$ have the same law, which is the solution $\mu$ of the Fokker-Planck equation (16) starting with $\mu_0$. By the law of large numbers, the empirical distribution $\bar{\mu}^N$ of the coupling processes converges to $\mu$, and therefore it is enough to estimate the distance between $\bar{\mu}^N$ and $\mu^N$. This can be done by obtaining estimates of the form

$$ \mathbb{E}[\|A_k - \bar{A}_k\|^2] \le \frac{\exp(CT)}{N}, \quad (18) $$

for all $1 \le k \le N$, which gives control on the 2-Wasserstein distance between $\bar{\mu}^N$ and $\mu^N$ on the time interval $[0, T]$. $\square$

We know want to study the long time behaviour of the nonlinear Fokker-Planck equation (16), that we will rewrite in function of $f = \rho \mu$ (in that case, $\rho$ represents the total “mass” of $f$). Since $\rho \mathcal{J}[\mu] = \mathcal{J}[f]$, it therefore has the following form, without any parameter on the equation:

$$ \partial_t f = -\nabla_A \cdot (f P_{T_A} \mathcal{J}[f]) + \Delta_A f. \quad (19) $$

This is an equation of the form $\partial_t f = C[f]$ where $C[f]$ can also be written, using the definition (6) of the von Mises distribution $M_J$, under the following factorized form:

$$ C[f] = \nabla_A \cdot \left[ M_{\mathcal{J}[f]}(A) \nabla_A \left( \frac{f}{M_{\mathcal{J}[f]}(A)} \right) \right]. $$

In order to understand the long time behaviour of the solution, let us first look at stationary solutions.
**Proposition 2.2.** A measure $f$ on $SO_3(\mathbb{R})$ is a stationary solution of the Fokker-Planck equation (19) if and only if it is of the form $f = \rho M_J$, where $J$ satisfies the following compatibility equation

$$J = \rho \mathcal{J}[M_J]. \qquad (20)$$

*Proof.* Since we have, by integration by parts,

$$\int_{SO_3(\mathbb{R})} \frac{f}{M_{\mathcal{J}[f]}(A)} C[f] dA = - \int_{SO_3(\mathbb{R})} \left\| \nabla_A \left( \frac{f}{M_{\mathcal{J}[f]}(A)} \right) \right\|^2 M_{\mathcal{J}[f]}(A) dA,$$

we immediately get that if $C[f] = 0$ then $f$ has to be proportional to $M_{\mathcal{J}[f]}$, and the total mass of $f$, denoted by $\rho$, gives the coefficient of proportionality. Then, taking the average on $SO_3(\mathbb{R})$ against $A$, thanks to the definition (15) of $\mathcal{J}$, we obtain, denoting $J = \mathcal{J}[f]$:

$$J = \mathcal{J}[f] = \mathcal{J}[\rho M_{\mathcal{J}[f]}] = \rho \mathcal{J}[M_J],$$

which is the compatibility equation for $J$. Conversely, if $J$ is a fixed point of this map $J \mapsto \rho \mathcal{J}[M_J]$, then setting $f = \rho M_J$, we get $\mathcal{J}[f] = J$, and then $C[f] = 0$. $\square$

Before obtaining a simple characterization of the solutions of the compatibility equation (20), which is the object of the next section, let us give some more results on the solutions to the Fokker-Planck equation (19).

**Proposition 2.3.** For all nonnegative measure $f_0$ on $SO_3(\mathbb{R})$, with total mass $\rho > 0$, there exists a unique weak solution $f$ to the nonlinear Fokker-Planck equation (19) such that $f(t, \cdot)$ converges to $f_0$ (in Wasserstein distance) as $t \to 0$. This solution belongs to $C^\infty((0, +\infty), SO_3(\mathbb{R}))$ and is positive for any positive time. Furthermore, we have the following uniform estimates in time: for all $t_0 > 0$, and $s \in \mathbb{R}$, the solution $f$ is uniformly bounded on $[t_0, +\infty)$ in the Sobolev space $H^s(SO_3(\mathbb{R}))$.

The proof of this proposition can be obtained through simple energy estimates in $H^s(SO_3(\mathbb{R}))$, using Poincaré inequalities for high modes and the fact that the low modes are uniformly bounded in time. Indeed, the nonlinearity in the Fokker-Planck equation (19) is only through $\mathcal{J}[f]$, which is uniformly bounded thanks to its definition (15) and the fact that $SO_3(\mathbb{R})$ is compact, together with the fact that the total mass $\rho$ is preserved. The positivity comes from the maximum principle. We refer to [FL12] to a detailed proof of such results on the unit sphere instead of $SO_3(\mathbb{R})$, for which all the arguments may be used similarly.

Let us now describe the free energy associated to this Fokker-Planck equation, which may be rewritten

$$\partial_t f = \nabla_A \cdot (f \nabla_A (\ln f - A \cdot \mathcal{J}[f]))$$

Multiplying by $\ln f - A \cdot \mathcal{J}[f]$ and integrating over $SO_3(\mathbb{R})$, the left-hand side of the equality can be seen as a time derivative, and the right-hand side can be integrated by parts, to obtain the following dissipation relation:

$$\frac{d}{dt} \mathcal{F}[f] + \mathcal{D}[f] = 0, \qquad (21)$$
where

$$ \mathcal{F}[f] = \int_{SO_3(\mathbb{R})} f(A) \ln f(A) dA - \frac{1}{2} \|\mathcal{J}[f]\|^2, \quad (22) $$

$$ \mathcal{D}[f] = \int_{SO_3(\mathbb{R})} f(A) \|\nabla_A(\ln f - A \cdot \mathcal{J}[f])\|^2 dA. \quad (23) $$

We can then prove, as in [FL12] that being a stationary state of the Fokker-Planck equation (see Proposition 2.2) is equivalent to be a critical point of $\mathcal{F}$ under the constraint of mass $\rho$, and that is also equivalent to be a function with no dissipation ($\mathcal{D}[f] = 0$).

We then have a decreasing free energy $\mathcal{F}[f]$, and thanks to a kind of LaSalle's principle, we obtain that the solution converges to a set of equilibria:

**Proposition 2.4.** Let $f_0$ be a nonnegative measure on $SO_3(\mathbb{R})$ with mass $\rho > 0$. We denote by $\mathcal{F}_\infty$ the limit of $\mathcal{F}[f(t, \cdot)]$ as $t \to +\infty$, where $f$ is the solution to the Fokker-Planck equation (19) with initial condition $f_0$. Then the set of equilibria $\mathcal{E}_\infty$, given by

$$ \mathcal{E}_\infty = \{ \rho M_J \text{ such that } J = \rho \mathcal{J}[M_J] \text{ and } \mathcal{F}[\rho M_J] = \mathcal{F}_\infty \}, $$

is not empty. Furthermore, the solution $f$ converges in any Sobolev space $H^s$ to this set of equilibria in the following sense:

$$ \liminf_{t \to \infty} \inf_{g \in \mathcal{E}_\infty} \|f(t, \cdot) - g\|_{H^s} = 0. $$

Once more, the proof of this proposition follows exactly the one given in [FL12]. The important point of this proposition is that once the structure of the solutions of the compatibility equation (20) is known (which is the aim of the next section), it gives a lot of information on the large time behaviour of the solutions to the Fokker-Planck equation.

Before giving a precise description of these solutions, let us remark that $J = 0$ is always a solution to the compatibility equation, since $\mathcal{J}[\rho] = 0$, therefore the uniform distribution with mass $\rho$ is a steady-state. We want to expand the free energy $\mathcal{F}$ around this steady-state. We will need the following lemma (Lemma 3.3 of [DDFMA20]):

**Lemma 1.** For all $J \in M_3(\mathbb{R})$,

$$ \int_{SO_3(\mathbb{R})} (J \cdot A) A \, dA = \frac{1}{6} J. \quad (24) $$

Consequently, if $f$ is a finite measure, the orthogonal projection of $f$ on the space of functions of the form $A \mapsto J \cdot A$ for $J \in M_3(\mathbb{R})$ is given by $A \mapsto 6\mathcal{J}[f] \cdot A$. Now, let us take a nonnegative measure $f$ with mass $\rho$, we write $J = \mathcal{J}[f]$ and $g(A) = 6J \cdot A$. We suppose that $\|J\|$ is sufficiently small, so that $\rho + g > 0$ on $SO_3(\mathbb{R})$. We write $h = f - \rho - g$, so $h$ is a finite measure with zero average and $\mathcal{J}[h] = 0$. Then
we obtain, by convexity of $x \mapsto x \ln x$ on $\mathbb{R}_+$:

$$ 
\begin{aligned}
\mathcal{F}[f] &\geq \int_{SO_3(\mathbb{R})} [(\rho + g(A)) \ln(\rho + g(A)) + h(A)(\ln(\rho + g(A)) + 1)] dA - \frac{1}{2} \|J\|^2 \\
&\geq \mathcal{F}[\rho + g] + \int_{SO_3(\mathbb{R})} h(A) \left(1 + \ln \rho + \frac{g(A)}{\rho}\right) dA - O(\|g\|_{\infty}^2) \int_{SO_3(\mathbb{R})} |h(A)| dA. \\
&\geq \mathcal{F}[\rho + g] - O(\|J\|^2) \left( \int_{SO_3(\mathbb{R})} |f(A) - \rho| dA + O(\|J\|) \right). \quad (25)
\end{aligned}
$$

Next we compute

$$ 
\begin{aligned}
\mathcal{F}[\rho + g] &= \rho \ln \rho + \frac{1}{2\rho} \int_{SO_3(\mathbb{R})} (6A \cdot J)^2 dA + O(\|g\|_{\infty}^3) - \frac{1}{2} \|J\|^2 \\
&= \mathcal{F}[\rho] + \frac{6 - \rho}{2\rho} \|J\|^2 + O(\|J\|^3), \quad (26)
\end{aligned}
$$

thanks to Lemma 24. We therefore see that the sign of $6 - \rho$ plays a role to study the nature, as a critical point of $\mathcal{F}$, of the uniform distribution of mass $\rho$:

**Proposition 2.5.** We set $\rho_c = 6$.

* If $\rho < \rho_c$, then the uniform distribution with mass $\rho$ is a local strict minimizer of the free energy $\mathcal{F}$ under the constraint of total mass $\rho$.

* If $\rho > \rho_c$, the uniform distribution with mass $\rho$ is not a local minimizer of the free energy $\mathcal{F}$ under the constraint of total mass $\rho$.

*Proof.* When $\rho > 6$, it is clear thanks to (25) and (26) that if $\|J\|$ and $\int_{SO_3(\mathbb{R})} |f - \rho|$ are sufficiently small and $J \neq 0$, then $\mathcal{F}[f] > \mathcal{F}[\rho]$. If $J = 0$ but $f \neq \rho$, then by strict convexity of $x \mapsto x \ln x$ on $\mathbb{R}_+$, we get

$$ \mathcal{F}[f] = \int_{SO_3(\mathbb{R})} f(A) \ln f(A) dA > \int_{SO_3(\mathbb{R})} (\rho \ln \rho + (f(A) - \rho) \ln \rho) dA = \mathcal{F}[\rho]. $$

The second point follows directly from (26).

<figure><img src="image_4.png" /></figure>

This last proposition gives an insight on the stability of the uniform steady-state (we will indeed see later that this uniform steady-state is isolated). In summary, we have shown that there is a phenomenon of phase transition at the threshold $\rho = \rho_c$, and we know thanks to Proposition 2.4 that there must exist other types of steady-states, at least when $\rho > \rho_c$. We are now ready to give a precise description of those non-isotropic equilibria.

## 3 Link with higher dimensional polymers, solutions to the compatibility equation

This section is the summary of the results we obtained in [DDFMA20] to solve the compatibility equation (20) (in a slightly different context, see Section 4), therefore we will omit the proofs.
Let us first recall some definitions. We denote by $\mathbb{H}$ the set of quaternions: objects of the form $q = a + bi + c\mathbf{j} + d\mathbf{k}$, where $(a, b, c, d) \in \mathbb{R}^4$ and the imaginary quaternions satisfy $\mathbf{i}^2 = \mathbf{j}^2 = \mathbf{k}^2 = \mathbf{ijk} = -1$. For such a quaternion $q$, we denote by $q^* = a - bi - c\mathbf{j} - d\mathbf{k}$ its conjugate. It satisfies $qq^* = q^*q = a^2 + b^2 + c^2 + d^2 = |q|^2$, if we identify the Euclidean space $\mathbb{R}^4$ with $\mathbb{H}$. We denote then by $\mathbb{H}_1$ the set of units quaternions: those for which $|q|^2 = 1$.

We say that a quaternion $q$ of the previous form is purely imaginary if its real part $a$ is zero. It allows now to identify $\mathbb{R}^3$ with the set of purely imaginary quaternions. We will use boldface letters when using this identification.

The first proposition is a link between $SO_3(\mathbb{R})$ and $\mathbb{H}_1/\{\pm 1\}$.

**Proposition 3.1.** For any $q \in \mathbb{H}_1$, the linear map $\mathbf{u} \mapsto q\mathbf{u}q^*$ sends purely imaginary quaternions on purely imaginary quaternions of the same norm. It is therefore identified as a rotation of $\mathbb{R}^3$, and the corresponding rotation matrix is denoted $\Phi(q)$. Conversely for any rotation matrix $A \in SO_3(\mathbb{R})$, there exists a unit quaternion $q$ such that $A = \Phi(q)$ (this quaternion is not unique, the only other possibility being $-q$). The map $\Phi$ can then be seen as a group isomorphism between $SO_3(\mathbb{R})$ and $\mathbb{H}_1$ (this is actually a local isometry between the manifolds). In practice, the matrix $R(\theta, \mathbf{n})$ given by Rodrigues' formula (1) corresponds to the quaternion $q = \cos(\theta/2) + \sin(\theta/2)\mathbf{n}$ (remember that vectors in $\mathbb{R}^3$ are seen as purely imaginary quaternions, and remark that if we replace $\theta$ by $\theta + 2\pi$, we get the same rotation matrix, but the opposite quaternion).

This allows to represent a rotation matrix by a unit quaternion up to multiplication by $\pm 1$. This is reminiscent of describing rodlike polymers as unit vectors up to multiplication by $\pm 1$, but generalized in dimension 4. This analogy was the starting point of our work [DFMAT18], where we used those unit quaternions for the modeling of alignment of rigid bodies. In the present case, we will see that this analogy will actually be very helpful, by transforming the compatibility equation (20) into another one which has already been solved in [WH08], in the context of suspensions of diluted polymers.

We denote by $S_4^0(\mathbb{R})$ the space of symmetric and trace-free matrices of dimension 4, which are called $Q$-tensors. To a unit quaternion $q$, we can associate the $Q$-tensor given by $q \otimes q - \frac{1}{4}I_4$. Remark that two unit quaternions $q$ and $\tilde{q}$ are associated to the same $Q$-tensor if and only if $q = \pm \tilde{q}$ (this is a unit vector in the eigenspace of this $Q$-tensor associated to the eigenvalue $\frac{3}{4}$, which is one-dimensional). So we have another way to represent unit quaternions up to multiplication by $\pm 1$ in this space. The important fact to notice is that those two embeddings are actually the same, up to a linear isomorphism between the spaces $M_3(\mathbb{R})$ and $S_4^0(\mathbb{R})$, which has nice properties.

**Proposition 3.2.** There exists a linear isomorphism $\phi$ between the spaces $M_3(\mathbb{R})$ and $S_4^0(\mathbb{R})$ (both of dimension 9) with the following properties:

$$ \forall q \in \mathbb{H}_1, \quad \phi(\Phi(q)) = q \otimes q - \frac{1}{4}I_4, \qquad (27) $$

$$ \forall J \in M_3(\mathbb{R}), \forall q \in \mathbb{H}_1, \quad \frac{1}{2}J \cdot \Phi(q) = q \cdot \phi(J)q, \qquad (28) $$
where the map $\Phi$ is given by Proposition 3.1. The dot product in the left-hand side of (28) is the metric in the space $M_3(\mathbb{R})$ given in (3), while the one in the right-hand side is the canonical scalar product of $\mathbb{R}^4$. Furthermore, the isomorphism $\phi$ preserves the diagonal structure: $J \in M_3(\mathbb{R})$ is diagonal if and only if $\phi(J)$ is diagonal in $S_4^0(\mathbb{R})$.

The proof of this proposition is done in [DDFMA20]. The expression (28) is actually the definition of $\phi$: the left-hand side is a quadratic form in $q$ (seen as an element of $\mathbb{R}^4$), defined for any unit quaternion, which defines a symmetric bilinear form on all quaternions, the matrix of which is $\phi(J)$. The expression of $\phi(J)$ is given in the appendix of [DDFMA20], which gives the fact that it is bijective and with values in trace-free matrices, and the provides the property (27). With this isomorphism, we can rewrite the compatibility equation in the framework of $Q$-tensors. For a finite measure $f$ on $SO_3(\mathbb{R})$

$$ \mathcal{Q}[f] = \int_{\mathbb{H}_1} f(q)(q \otimes q - \frac{1}{4}I_4)dq. $$

Therefore, thanks to the definition (15) of $\mathcal{J}$ and the fact that $\Phi$ is a local isometry, we obtain, for a finite measure $f$ on $SO_3(\mathbb{R})$

$$ \phi(\mathcal{J}[f]) = \int_{SO_3(\mathbb{R})} \phi(A)f(A)dA = \int_{\mathbb{H}_1} \phi(\Phi(q))f(\Phi(q))dq = \mathcal{Q}[f \circ \Phi]. $$

Finally, we also define the generalized von Mises associated to $Q \in S_4^0(\mathbb{R})$ by

$$ M_Q(q) = \frac{1}{\mathcal{Z}(Q)} \exp(q \cdot Qq), \text{ where } \mathcal{Z}(Q) = \int_{\mathbb{H}_1} \exp(q \cdot Qq)dq, $$

where we use the same notation as in (6) for the generalized von Mises on $SO_3(\mathbb{R})$, but it will always be clear following the context which definition is concerned. Using the property (28), it is then clear that $M_J(\Phi(q)) = M_{2\phi(J)}(q)$. Therefore, the compatibility equation (20) becomes, writing $Q = 2\phi(J)$:

$$ Q = 2\phi(J) = 2\rho \phi(\mathcal{J}[M_J]) = 2\rho \mathcal{Q}[M_Q]. $$

It happens that this equation is exactly the compatibility equation that we obtain when we try to obtain the steady states of the following Fokker-Planck equation, for a probability measure $\mu$ on $\mathbb{H}_1$:

$$ \partial_t \mu = -2\rho \nabla_q \cdot (\mu \nabla_q (q \cdot \mathcal{Q}[\mu]q)) - \Delta_q \mu. $$

This corresponds to the Smoluchowski (or Doi-Onsager) equation for suspensions of dilute rodlike polymers with Maier-Saupe potential of strength $2\rho$, and is nothing else than our Fokker-Planck equation (14), up to a change of variable thanks to the map $\Phi$. It happens that this compatibility equation has been studied a lot in dimension 3 (instead of 4 here), with the independent works [CKT04, FS05, LZZ05]. And in the work [WH08], a unified approach has been proposed, which allows to treat the case of higher dimensional space. The main result is that a solution $Q \in S_n^0(\mathbb{R})$ of
the compatibility equation $Q = \alpha Q[M_Q]$ can have at most two different eigenvalues. In dimension 4, it means that if $Q$ is different from zero, there are only two cases: either one eigenvalue is simple and the other one is triple, or both are double. In the first case, if we take $q$ a unit quaternion in the eigenspace of dimension one, we get that $Q$ is proportional to $q \otimes q - \frac{1}{4}I_4$, which means that $J = \phi^{-1}(Q)$ is proportional to the rotation matrix $\Phi(q)$. And indeed it is possible to see that if $A_0$ is a rotation matrix and $\alpha \in \mathbb{R}$, then $J[M_{\alpha A_0}]$ is proportional to $A_0$, with a coefficient $c_1(\alpha)$ (that can be expressed using an appropriate volume form on $SO_3(\mathbb{R})$ and will be given later on). Therefore the compatibility equation (20) becomes the one-dimensional equation $\alpha = \rho c_1(\alpha)$. For the second case, it is a little bit more subtle, but it still leads to a one-dimensional equation of the form $\alpha = \rho c_2(\alpha)$. The results are summarized in the following proposition (corresponding to Theorem 5 of [DDFMA20]):

**Proposition 3.3.** *The solutions to the compatibility equation (20) are:*

* The matrix $J = 0$,

* the matrices of the form $J = \alpha A_0$ with $A_0 \in SO_3(\mathbb{R})$ and where $\alpha \in \mathbb{R} \setminus \{0\}$ satisfies the scalar compatibility equation

$$ \alpha = \rho c_1(\alpha), \qquad (29) $$

* the matrices of the form $J = \alpha \sqrt{3} \mathbf{a}_0 \otimes \mathbf{b}_0$ where $\mathbf{a}_0$ and $\mathbf{b}_0$ are two unit vectors of $\mathbb{R}^3$ and $\alpha > 0$ satisfies the scalar compatibility equation

$$ \alpha = \rho c_2(\alpha), \qquad (30) $$

with the functions $c_1$ and $c_2$ given by

$$ c_1(\alpha) = \frac{\int_0^\pi \frac{1}{3}(2 \cos \theta + 1) \sin^2(\frac{\theta}{2}) \exp(\alpha \cos \theta) d\theta}{\int_0^\pi \sin^2(\frac{\theta}{2}) \exp(\alpha \cos \theta) d\theta}, \\ c_2(\alpha) = \frac{1}{\sqrt{3}} \frac{\int_0^\pi \cos \varphi \sin \varphi \exp(\frac{\sqrt{3}}{2} \alpha \cos \varphi) d\varphi}{\int_0^\pi \sin \varphi \exp(\frac{\sqrt{3}}{2} \alpha \cos \varphi) d\varphi}. $$

Compared to the convention taken in [DDFMA20], we chose to add the constant $\sqrt{3}$ in the last type of solutions (changing accordingly the expression of $c_2(\alpha)$). The reason is that if $J$ is a solution to the compatibility equation (20), where $\alpha$ satisfies (29) or (30), then $\|J\|^2 = \frac{3}{2}\alpha^2$. The order parameter $c$ associated to the steady state $\rho M_J$ by the formula (12) is then equal to $\frac{|\alpha|}{\rho}$ which is $|c_1(\alpha)|$ or $|c_2(\alpha)|$. These functions $c_1$ and $c_2$ then provide the values of the order parameter of the considered steady-state. The study of these functions (and more precisely the behaviour of $\frac{\alpha}{c_1(\alpha)}$ and $\frac{\alpha}{c_2(\alpha)}$) is the key to provide a complete description of the possible steady-states. Once more, the following proposition is taken from [DDFMA20].

**Proposition 3.4.** *The functions $c_1$ and $c_2$ are both strictly increasing on $\mathbb{R}$ having value 0 at 0. Therefore 0 is always a solution to the scalar compatibility equations (29) and (30). If we set $\rho_c = 6$, then when $\alpha \to 0$, the functions $\rho_1 : \alpha \mapsto \frac{\alpha}{c_1(\alpha)}$ and $\rho_2 : \alpha \mapsto \frac{\alpha}{c_2(\alpha)}$ both have a limit equal to $\rho_c$. Furthermore:*
* There exists $\alpha^* > 0$ such that $\rho_1$ is decreasing on $(-\infty, \alpha^*]$ and increasing on $[\alpha^*, +\infty)$, converging to $+\infty$ at $\pm\infty$. We set $\rho^* = \rho_1(\alpha^*)$ (which is less than $\rho_c$). For all $\rho \ge \rho^*$, we define $\alpha_1^>(\rho)$ (resp. $\alpha_1^-(\rho)$) to be the unique value of $\alpha \ge \alpha^*$ (resp $\alpha \le \alpha^*$) such that $\rho_1(\alpha) = \rho$. Finally, we define $\tilde{c}_1^>(\rho) = c_1(\alpha_1^>(\rho))$ and $\tilde{c}_1^-(\rho) = c_1(\alpha_1^-(\rho))$. Setting $c^* = c_1(\alpha^*)$, the function $\tilde{c}_1^-$ (resp. $\tilde{c}_1^+$) is increasing (resp. decreasing) on $[\rho^*, +\infty)$, with value $c^*$ at $\rho^*$, and converging to 1 (resp. $-\frac{1}{3}$) at $+\infty$.

Numerically, we obtain $\alpha^* \approx 1.9395$, $\rho^* \approx 4.5832$, and $c^* \approx 0.4232$.

* The function $\rho_2$ is (even and) increasing on $[0, +\infty)$, converging to $+\infty$ at $+\infty$. For all $\rho \ge \rho_c$, we define $\alpha_2(\rho)$ to be the unique value of $\alpha \ge \alpha^*$ such that $\rho_2(\alpha) = \rho$. Finally, we define $\tilde{c}_2(\rho) = c_2(\alpha_2(\rho))$. The function $\tilde{c}_2$ is increasing on $[\rho_c, +\infty)$, with value 0 at $\rho_c$ and converging to $\frac{1}{\sqrt{3}}$ at $+\infty$.

Figure (3) depicts a plot of these functions $\tilde{c}_1^>$ (solid), $\tilde{c}_1^-$ (dashed), and $\tilde{c}_2$ (dashed-dot line), in log-scale for $\rho \in [2, 40]$. They represent the order parameters (up to sign) of the different families of steady-states. We also drew a solid line at level 0 for $\rho < \rho_c$ and a dotted line at level 0 for $\rho > \rho_c$, corresponding to the order parameter of the uniform steady-state (and illustrating the result of Proposition 2.5 regarding its stability).

<figure><img src="image_5.png" /><figcaption>Figure 3: Behaviors of the functions $\tilde{c}_1^>$ (solid line), $\tilde{c}_1^-$ (dashed line) and $\tilde{c}_2$ (dashed-dot line).</figcaption></figure>
We can therefore describe more precisely the long time behaviour of the solution to the Fokker-Planck equation according to the value of $\rho$, thanks to Proposition 2.4.

**Theorem 1.** Let $f_0$ be a nonnegative measure with mass $\rho > 0$, and $f$ the solution to the Fokker-Planck equation (19) with initial condition $f_0$. For the following statements, the notion of convergence is with respect to any $H^s$ norm on $SO_3(\mathbb{R})$.

* If $\rho < \rho^*$, the only steady-state is the uniform distribution on $SO_3(\mathbb{R})$, and the solution $f(t, \cdot)$ converges to this steady state as $t \to +\infty$.

* If $\rho^* \le \rho \le \rho_c$, there are three families of steady-states (two of which are equal when $\rho = \rho^*$ or $\rho = \rho_c$), and $f(t, \cdot)$ converges to one of these families:

- either there exists $A_0(t) \in SO_3(\mathbb{R})$ such that $f(t, \cdot) - \rho M_{\alpha_1^\uparrow(\rho)A_0(t)}$ converges to zero,

- either $f(t, \cdot)$ converges to the uniform distribution on $SO_3(\mathbb{R})$,

- or there exists $A_0(t) \in SO_3(\mathbb{R})$ such that $f(t, \cdot) - \rho M_{\alpha_1^\downarrow(\rho)A_0(t)}$ converges to zero, as $t \to +\infty$.

* If $\rho > \rho_c$, there is an additional family of steady-states, and $f(t, \cdot)$ converges to one of these four families:

- either there exists $A_0(t) \in SO_3(\mathbb{R})$ such that $f(t, \cdot) - \rho M_{\alpha_1^\uparrow(\rho)A_0(t)}$ converges to zero,

- either $f(t, \cdot)$ converges to the uniform distribution on $SO_3(\mathbb{R})$,

- either there exists $A_0(t) \in SO_3(\mathbb{R})$ such that $f(t, \cdot) - \rho M_{\alpha_1^\downarrow(\rho)A_0(t)}$ converges to zero,

- or there exist unit vectors $\mathbf{a}_0(t), \mathbf{b}_0(t)$ such that $f(t, \cdot) - \rho M_{\alpha_2^\uparrow(\rho)\sqrt{3}\mathbf{a}_0(t)\otimes\mathbf{b}_0(t)}$ converges to zero, as $t \to +\infty$.

*Proof.* This result is a summary of the possible steady-states according to Proposition 2.2 and 3.3. The convergence of $f$ to one of these families comes from Proposition 2.4 and from the fact that, even if the limit set $\mathcal{E}_\infty$ of equilibria may consist of several distinct such families, they would belong to different connected components of $\mathcal{E}_\infty$. $\square$

Let us now try to understand the stability of each of these families of equilibria. Figure 4 is a zoom on the region $\rho \in [3, 8]$ of the plots of the functions $\tilde{c}_1^\uparrow, |\tilde{c}_1^\downarrow|$ and $\tilde{c}_2$ (remember that these functions are the order parameters of the corresponding steady-states), on top of the final values of the order parameters of the numerical simulations which were given in the right part of Figure 2. It suggests the only stable equilibria, apart from the uniform one when $\rho < \rho_c$, are those corresponding to the curve $\tilde{c}_1^\uparrow$. This is indeed what we will show in the next section.
<figure><img src="image_6.png" /><figcaption>Figure 4: Behaviors of the functions $\tilde{c}_1^\dagger$, $|\tilde{c}_1^\dagger|$ and $\tilde{c}_2$ and final order parameters of the numerical simulations.</figcaption></figure>

# 4 Stability results thanks to a BGK model

Instead of the Fokker-Planck equation (19), let us consider the following BGK equation:

$$ \partial_t f = \rho M_J \mathcal{J}[f] - f. \qquad (31) $$

This is still an equation where the total mass is preserved and for which the steady states satisfy the same compatibility equation: if $f$ is a steady-state, it has to be of the form $\rho M_J$ where $J = \mathcal{J}[f] = \rho \mathcal{J}[M_J]$. Therefore these two evolution equations share the same steady-states, which were determined in [DDFMA20] and summarized in the previous section. Let us now give a summary of the results of stability of these equilibria which were obtained in [DDFMA20]. It happens that these two evolution equations (BGK and Fokker-Planck) also share the same property of dissipation of the free energy $\mathcal{F}$: if $f$ is a positive solution to (31), then by multiplying both sides by $\ln f(A) - A \cdot \mathcal{J}[f]$ and integrating on $SO_3(\mathbb{R})$, we obtain

$$ \frac{d}{dt} \mathcal{F}[f] + \tilde{\mathcal{D}}[f] = 0, $$
where $\mathcal{F}[f]$ is given by (22) and

$$ \tilde{\mathcal{D}}[f] = \int_{SO_3(\mathbb{R})} (f - \rho M_{\mathcal{J}[f]})(\ln f - \ln(\rho M_{\mathcal{J}[f]}))dA \ge 0. $$

Then, by writing $J(t) = \mathcal{J}[f(t, \cdot)]$ where $f$ is a solution of the BGK equation (31), we obtain that $J$ satisfies an ordinary differential equation:

$$ \frac{d}{dt} J = \rho \mathcal{J}[M_J] - J. \qquad (32) $$

The long-time behaviour of the solution of the BGK equation is much simpler to study, since it can be reduced to the study of a finite dimensional ODE.

A further reduction can be done through the special singular value decomposition, for which we state a result which will be useful in the following.

**Proposition 4.1.** If $J \in M_3(\mathbb{R})$, we call Special Singular Value Decomposition (SSVD) of $J$ a decomposition of the form $J = PDQ$ where $D = \text{diag}(d_1, d_2, d_3)$ is a diagonal matrix satisfying $d_1 \ge d_2 \ge |d_3|$ and $P, Q \in SO_3(\mathbb{R})$.

Such a SSVD always exists, and the matrix $D$ is unique (the rotations $P$ and $Q$ may not be unique). Furthermore, we have

$$ \min_{A \in SO_3(\mathbb{R})} \|J - A\| = \|J - PQ\| = \|D - I_3\|. \qquad (33) $$

*Proof.* The existence and uniqueness can be obtained through the singular value decomposition, and modifying the orthogonal matrices if necessary to change the sign of the last entry of the diagonal part and get special orthogonal matrices, see [DDFMA20]. We now compute

$$ \|J - A\|^2 = \|D - P^\top A Q^\top\|^2 = \|D\|^2 - 2 B \cdot D + \frac{3}{2}, $$

where $B = P^\top A Q^\top$. Therefore minimizing $\|J - A\|$ for $A \in SO_3(\mathbb{R})$ amounts to maximizing $B \cdot D$, for $B \in SO_3(\mathbb{R})$. The set of diagonal parts of rotation matrices (seen as vectors of $\mathbb{R}^3$) is given by Horn's tetrahedron [Hor54]: this is the convex hull $\mathcal{T}$ of the points $(\pm 1, \pm 1, \pm 1)$ with an even number of minus signs. Therefore we want to maximize $x \cdot d$ for $x \in \mathcal{T}$ and $d = (d_1, d_2, d_3)$. This convex function reaches its maximum on extremal points of $\mathcal{T}$, that is to say on one of the vertices of $\mathcal{T}$. Since we have

$$ d_1 + d_2 + d_3 \ge d_1 - d_2 - d_3 \ge -d_1 + d_2 - d_3 \ge -d_1 - d_2 - d_3, $$

we see that the maximum is reached for $x = (1, 1, 1)$. Therefore the maximum of $B \cdot D$ for $B \in SO_3(\mathbb{R})$ is reached for $B = I_3$, which ends the proof$^4$. $\square$

<footnote>$^4$Let us remark that if $d_2 > -d_3$, the maximum of $x \cdot d$ is unique on $\mathcal{T}$ and since the only rotation matrix for which the diagonal part is $(1, 1, 1)$ is the identity matrix $I_3$, we get that the minimizer $PQ$ of (33) is unique. So even if $P$ and $Q$ may not be unique, in that case the matrix $PQ$ is unique, and could be seen as a Special Polar Decomposition of $J$ (with the analogy with the fact that if $\det J > 0$, then $J = PDQ$ is the singular value decomposition of $J$ and $PQ$ is the polar decomposition of $J$ [DFMA17]).</footnote>
With this definition of the SSVD, the reduction that can be done is that the flow of the ODE (32) preserves the SSVD: if a SSVD of the initial condition is given by $J(0) = PD_0Q$, then for all time $t$, we have the following SSVD: $J(t) = PD(t)Q$, with the same rotation matrices $P$ and $Q$, and where $D(t) = (d_1(t), d_2(t), d_3(t))$ is a diagonal matrix satisfying the same ODE (32) as $J$, with initial condition $D_0$ (the fact that the matrix is diagonal and the inequalities $d_1(t) \ge d_2(t) \ge |d_3(t)|$ are preserved by the flow of this ODE). We therefore only have to study a three-dimensional ODE. Finally, the last observation we can do is that the flow of the ODE (32) is actually a gradient flow of a potential: if we write

$$V(J) = \frac{1}{2} \|J\|^2 - \rho \ln \mathcal{Z}(J), \text{ where } \mathcal{Z}(J) = \int_{SO_3(\mathbb{R})} \exp(J \cdot A) dA, \quad (34)$$

as in the definition (6) of the generalized von Mises distribution, we obtain

$$\nabla V(J) = J - \rho \mathcal{J}[M_J], \quad (35)$$

where the gradient is taken with respect to the inner product of $M_3(\mathbb{R})$ given by (3).

Therefore the ODE (32) is simply $\frac{d}{dt}J = -\nabla V(J)$, and one can prove that any solution will converge to a critical point of $V$, which corresponds to a solution of the compatibility equation (20). We then obtain the same type of convergence as in Theorem 1. The main difference is that we have convergence to a unique steady-state (and not to a set of steady-states), that can be determined by knowing a special singular value decomposition of $\mathcal{J}[f_0]$. The other difference is that the convergence does not takes place in any Sobolev space $H^s$: the BGK equation is not regularizing in time. The following proposition is a summary of results in [DDFMA20]:

**Proposition 4.2.** Let $f_0$ be a finite nonnegative measure with mass $\rho > 0$, and $f$ the solution to the BGK equation (31) with initial condition $f_0$. We write the decomposition $\mathcal{J}[f_0] = P_0D_0Q_0$, where $P_0, Q_0 \in SO_3(\mathbb{R})$ and $D_0 = \text{diag}(d_{1,0}, d_{2,0}, d_{3,0})$, with $d_{1,0} \ge d_{2,0} \ge |d_{3,0}|$ (special singular value decomposition). Then for all $t \in \mathbb{R}$, we have $\mathcal{J}[f(t, \cdot)] = P_0D(t)Q_0$, where $D(t) = \text{diag}(d_1(t), d_2(t), d_3(t))$ is the solution to the ODE (32) with initial condition $D_0$, satisfying $d_1(t) \ge d_2(t) \ge |d_3(t)|$. In the following statements, the notion of convergence of $f(t, \cdot)$ is in the space of measures (or any normed space for which $f_0$ is an element and for which the map $f \mapsto \mathcal{J}[f]$ is continuous).

* If $\rho < \rho^*$, then $D(t) \to 0$ and $f(t, \cdot)$ converges to the uniform distribution as $t \to +\infty$.

* If $\rho^* \le \rho \le \rho_c$, there are three families of steady-states (two of which are equal when $\rho = \rho^*$ or $\rho = \rho_c$), and $f(t, \cdot)$ converges to one of these steady-states, as $t \to +\infty$:

- either $D(t) \to 0$, and $f(t, \cdot)$ converges to the uniform distribution,

- either $D(t) \to \alpha_1^\uparrow(\rho)I_3$, and $f(t, \cdot) \to \rho M_{\alpha_1^\uparrow(\rho)A_0}$ where $A_0 = P_0Q_0$,

- or $D(t) \to \alpha_1^\downarrow(\rho)I_3$, and $f(t, \cdot) \to \rho M_{\alpha_1^\downarrow(\rho)A_0}$ where $A_0 = P_0Q_0$.
* If $\rho > \rho_c$, there is an additional family of steady-state, and $f(t, \cdot)$ converges to one of these steady-states, as $t \to +\infty$:

- either $D(t) \to 0$, and $f(t, \cdot)$ converges to the uniform distribution,

- either $D(t) \to \alpha_1^\uparrow(\rho)I_3$, and $f(t, \cdot) \to \rho M_{\alpha_1^\uparrow(\rho)A_0}$ where $A_0 = P_0Q_0$,

- either $D(t) \to \alpha_1^\downarrow(\rho) \text{diag}(-1, -1, 1)$, and $f(t, \cdot)$ converges to $\rho M_{\alpha_1^\downarrow(\rho)A_0}$, where $A_0 = P_0 \text{diag}(-1, -1, 1)Q_0$

- or $D(t) \to \alpha_2(\rho) \text{diag}(\sqrt{3}, 0, 0)$, and $f(t, \cdot)$ converges to $\rho M_{\alpha_2(\rho)\sqrt{3}\mathbf{a}_0 \otimes \mathbf{b}_0}$, with $\mathbf{a}_0 = P_0\mathbf{e}_1$ and $\mathbf{b}_0 = Q_0^\top \mathbf{e}_1$ (where $\mathbf{e}_1$ is the first element of the canonical basis of $\mathbb{R}^3$).

We now turn to stability results. For convenience, we will denote $\bar{V}$ the restriction of $V$ to the space of diagonal matrices. Its Hessian Hess $\bar{V}$ is then a symmetric bilinear form on a space of dimension 3. Thanks to the study of the signature of this Hessian, we obtained in [DDFMA20] the characterization of the stability of all steady-states. The next proposition is a summary of these results (without details on the domains of convergence):

**Proposition 4.3.** The uniform steady-state for the BGK equation (31) corresponds to the critical point 0 of the potential $\bar{V}$ (and $V$).

* If $0 < \rho < \rho_c$, the Hessian Hess $\bar{V}(0)$ has signature (++++) (and so 0 is a local minimizer of $V$). Therefore the uniform steady-state is locally asymptotically stable (with exponential rate of convergence).

* If $\rho > \rho_c$, the signature is (---) (therefore 0 is not a local minimizer of $V$), and the uniform steady-state is unstable.

When $\rho \ge \rho^*$, the steady-states of the form $\rho M_{\alpha_1^\uparrow(\rho)A_0}$ (resp. $\rho M_{\alpha_1^\downarrow(\rho)A_0}$) with $A_0$ in $SO_3(\mathbb{R})$ (see Theorem 1) correspond to the critical points of the form $\alpha_1^\uparrow(\rho)A_0$ (resp. $\alpha_1^\downarrow(\rho)A_0$) of $V$. Their nature can be reduced to the study of the critical point $D_\infty^\uparrow = \alpha_1^\uparrow(\rho)I_3$ (resp. $D_\infty^\downarrow = \alpha_1^\downarrow(\rho)I_3$) of $\bar{V}$.

* If $\rho > \rho^*$, the Hessian Hess $\bar{V}(D_\infty^\uparrow)$ has signature (++++) (and so $\alpha_1^\uparrow(\rho)A_0$ is a local minimizer of $V$). Therefore the steady-states of the form $\rho M_{\alpha_1^\uparrow(\rho)A_0}$ are locally asymptotically stable (with exponential rate of convergence).

* If $\rho^* < \rho < \rho_c$ (resp. $\rho > \rho_c$), the Hessian Hess $\bar{V}(D_\infty^\downarrow)$ has signature (---++) (resp. (+---)) (therefore $\alpha_1^\downarrow(\rho)A_0$ is not a local minimizer of $V$), and the steady-states of the form $\rho M_{\alpha_1^\downarrow(\rho)A_0}$ are unstable.

When $\rho > \rho_c$, the steady-states of the form $\rho M_{\alpha_2(\rho)\sqrt{3}\mathbf{a}_0 \otimes \mathbf{b}_0}$ with $\mathbf{a}_0, \mathbf{b}_0 \in \mathbb{S}^2$ (see Theorem 1) correspond to the critical points of the form $\alpha_2(\rho)\sqrt{3}\mathbf{a}_0 \otimes \mathbf{b}_0$ of $V$, which reduces to the study of the critical point $D_\infty = \alpha_2(\rho)\text{diag}(1, 0, 0)$ of $\bar{V}$.

* The Hessian Hess $\bar{V}(D_\infty)$ has signature (++-) (therefore $\alpha_2(\rho)\sqrt{3}\mathbf{a}_0 \otimes \mathbf{b}_0$ is not a local minimizer of $V$), and the steady-states of the form $\rho M_{\alpha_2(\rho)\sqrt{3}\mathbf{a}_0 \otimes \mathbf{b}_0}$ are unstable.
Furthermore, the critical cases are unstable: the uniform steady-state is unstable for $\rho = \rho_c$, and the steady-states of the form $\rho M_{\alpha^*A_0}$ are unstable when $\rho = \rho^*$ (the corresponding matrices $J = 0$ or $J = \alpha^*A_0$ are not local minimizers of $V$).

The main object of this section is to show, as it was claimed in Remark 5.5 of [DDFMA20], that we can directly use these results of (in)stability for the BGK equation (and more precisely for the potential $V$) to obtain (in)stability results for the Fokker-Planck equation, in order to complete the results around the uniform distribution given by Proposition 2.5. We provide a proposition and a theorem which give details on this statement.

The first proposition allows to compare the behaviours of $V$ and of $J \mapsto \mathcal{F}[\rho M_J]$.

**Proposition 4.4.** Let us define for $J \in M_3(\mathbb{R})$

$$W(J) = \mathcal{F}[\rho M_J],$$

Then, we have that $\nabla W(J) = 0$ if and only if $\nabla V(J) = 0$, that is to say $J$ is a solution to the compatibility equation (20). Furthermore, if $J$ is such a critical point, the Hessian Hess $W$ has the same signature as Hess $V$ (and more precisely, if $\overline{W}$ is the restriction of $W$ to the diagonal matrices, then Hess $\overline{W}$ and Hess $\overline{V}$ have the same signature).

*Proof.* We first compute

$$ 
\begin{aligned}
W(J) &= \int \rho(\ln \rho + A \cdot J - \ln \mathcal{Z}(J))M_J(A)dA - \frac{\rho^2}{2} \|\mathcal{J}[M_J]\|^2 \\
&= \rho \ln \rho - \ln \mathcal{Z}(J) + \frac{1}{2}\|J\|^2 - \frac{1}{2}\|J - \rho \mathcal{J}[M_J]\|^2 \\
&= V(J) - \frac{1}{2}\|\nabla V(J)\|^2 + \rho \ln \rho,
\end{aligned}
$$

thanks to (35). Therefore we obtain

$$\nabla W(J) = \nabla V(J) - \text{Hess } V(J)(\nabla V(J)). \quad (36)$$

We want to compute the Hessian of $V$, seen as a linear mapping from $M_3(\mathbb{R})$ to $M_3(\mathbb{R})$, symmetric with respect to the inner product of $M_3(\mathbb{R})$. Let us take $H$ small in $M_3(\mathbb{R})$. We first have that $\mathcal{Z}(J+H) = (1+\mathcal{J}[M_J] \cdot H)\mathcal{Z}(J) + O(\|H\|^2)$. Thus we get $M_{J+H}(A) = (1+A \cdot H - \mathcal{J}[M_J] \cdot H)M_J(A) + O(\|H\|^2)$. Finally we obtain

$$\mathcal{J}[M_{J+H}] = \mathcal{J}[M_J] - (\mathcal{J}[M_J] \cdot H)\mathcal{J}[M_J] + \int_{SO_3(\mathbb{R})} A(A \cdot H)M_J(A)dA + O(\|H\|^2).$$

Now, using the expression (35) of $\nabla V$, we get

$$\text{Hess } V(J)(H) = H - \rho \left[ (\mathcal{J}[M_J] \cdot H)\mathcal{J}[M_J] - \int_{SO_3(\mathbb{R})} A(A \cdot H)M_J(A)dA \right].$$
Said differently, seeing now Hess V as a symmetric bilinear form on $M_3(\mathbb{R})$:

$$ \begin{aligned} \text{Hess } V(J)(H, H) &= \|H\|^2 - \rho \left[ (\mathcal{J}[M_J] \cdot H)^2 - \int_{SO_3(\mathbb{R})} (A \cdot H)^2 M_J(A) dA \right] \\ &= \|H\|^2 - \rho \int \left[ (A - \mathcal{J}[M_J]) \cdot H \right]^2 M_J(A) dA, \end{aligned} \quad (37) $$

and we see that all the eigenvalues of Hess V are strictly less than 1. Therefore the (symmetric) linear mapping $Id - \text{Hess } V$ from $M_3(\mathbb{R})$ to $M_3(\mathbb{R})$ has only strictly positive eigenvalues, and is therefore an isomorphism. The expression (36) of $\nabla W$ then provides the equivalence between critical points for V and for W.

Finally, at a point J for which $\nabla V(J) = 0$, we obtain

$$ \text{Hess } W(J) = \text{Hess } V(J) - [\text{Hess } V(J)]^2. $$

Therefore, the eigenvalues of Hess $W(J)$ are given by $\lambda(1 - \lambda)$, where $\lambda$ are the eigenvalues of Hess $W(J)$, which all satisfy $\lambda < 1$. Therefore their signs are the same. And this is also true when restricted to the space of diagonal matrices. $\square$

We can now state the final theorem of this section.

**Theorem 2.** The nature of all the critical points of the free energy $\mathcal{F}$ is given by the following statements.

* For $\rho < \rho_c$, the uniform equilibrium of mass $\rho$ is a local strict minimizer of the free energy $\mathcal{F}$.

* For $\rho > \rho^*$, the set $\mathcal{E} = \{\rho M_{\alpha_1(\rho)A_0}, A_0 \in SO_3(\mathbb{R})\}$ is a local strict minimizer of the free energy $\mathcal{F}$, in the sense that there exists a neighborhood $\mathcal{V}$ of $\mathcal{E}$ (in the space of nonnegative measures of mass $\rho$) such that if $f \in \mathcal{V} \setminus \mathcal{E}$, then $\mathcal{F}[f] > \mathcal{F}_\infty$, where $\mathcal{F}_\infty$ is the common value of $\mathcal{F}$ on $\mathcal{E}$.

* For $\rho \ge \rho_c$, the uniform equilibrium of mass $\rho$ is not a local minimizer of the free energy $\mathcal{F}$.

* For $\rho \ge \rho^*$ (and $\rho \ne \rho_c$), any steady-state of the form $\rho M_{\alpha_1(\rho)A_0}$ for $A_0 \in SO_3(\mathbb{R})$ is not a local minimizer of the free energy $\mathcal{F}$.

* For $\rho > \rho_c$, any steady-state of the form $\rho M_{\alpha_2(\rho)\sqrt{3}a_0\otimes b_0}$ for $A_0 \in SO_3(\mathbb{R})$ is not a local minimizer of the free energy $\mathcal{F}$.

Therefore, the last three families of steady-states are unstable for the Fokker-Planck equation (19): there exist initial conditions arbitrarily close to these families (in any $H^s$ norm), such that the solution to the Fokker-Planck equation converges in long time towards another family of equilibria (see Theorem 1).

*Proof.* The first point has been proven in Proposition 2.5. For the second one, if it was not true, there would exist $f_0$ as close as we want from $\mathcal{E}$ such that $\mathcal{F}(f_0) \le \mathcal{F}_\infty$, and $f_0 \notin \mathcal{E}$. Since the different families of steady-states are isolated, $f_0$ cannot be
a steady-state. By letting $f$ be the solution of the BGK equation with initial condition $f_0$, we would have $\tilde{Q}[f_0] > 0$ and therefore $\mathcal{F}[f(t, \cdot)] < \mathcal{F}_\infty$ for all $t > 0$. Combined with the fact that $\mathcal{F}[f(t, \cdot)]$ is nonincreasing in time, this would be in contradiction with the fact that $f(t, \cdot)$ converges towards the set $\mathcal{E}$, thanks to the asymptotic stability of those steady-states for the BGK equation given by Proposition 4.3. Let us remark that the first point of the theorem could be proven in the same way, without having to expand the free energy, but only using the known results for the BGK equation and the fact that $\mathcal{F}$ is nonincreasing.

To prove the last three points, let us take such a steady state, of the form $\rho M_{J_0}$. We want to prove that $J_0$ is not a local minimizer of $W$, therefore $\rho M_{J_0}$ is not a local minimizer of $\mathcal{F}$. We write a SSVD of the form $J_0 = PD_0Q$ where $D_0$ is a diagonal matrix and $P, Q \in SO_3(\mathbb{R})$. If $J = PDQ$ where $D$ is a diagonal matrix close to $D_0$, then $W(J) = \bar{W}(D)$. Therefore we only need to prove that $D_0$ is not a local minimizer of $\bar{W}$. In the case where $\rho \neq \rho_c$ and $\rho \neq \rho^*$, since the signature of Hess $\bar{W}(D_0)$ has negative components (thanks to Propositions 4.3 and 4.4), we directly get the results. In the critical cases we will use a mountain-pass lemma argument. In the case where $\rho = \rho_c$, suppose that $0$ is a local minimizer of $\bar{W}$. Then it is a local strict minimizer, since this critical point is isolated. Therefore by looking at the other local strict minimizer $\alpha_1^\uparrow(\rho_c)I_3$ of $\bar{W}$ (for which the signature of the Hessian is $(+++)$, thanks again to Propositions 4.3 and 4.4), we would obtain, by the mountain-pass lemma, a third critical point $D$ of $\bar{W}$, which would satisfy $\bar{W}(D) > \max(\bar{W}(0), \bar{W}(\alpha_1^\uparrow(\rho_c)I_3))$. This is in contradiction with the fact that we only have two families of equilibria for this value of $\rho$. The same argument can be used to show that when $\rho = \rho^*$, the point $\alpha^*I_3$ is not a local minimizer of $\bar{W}$, using as other local strict minimizer the point $0$.

The conclusion of the statement of the theorem comes from the fact that we actually proved that the critical points were not local minimizers of $W$, which is the evaluation of $\mathcal{F}$ on smooth functions of the form $\rho M_J$, so the $H^s$ norm of $\rho M_J - \rho M_{J_0}$ is small when $J$ is close to $J_0$. $\square$

For the first two points of Theorem 2, we did not provide the corresponding stability results. Indeed, in the next section, a more detailed study will show that they are exponentially stable.

# 5 Exponential convergence for the stable steady-states

We will now show that the two families of steady-states that correspond to what we observe in the numerical simulations are locally exponentially attracting. In particular, when $f$ is a solution to the Fokker-Planck equation in the neighborhood of those steady-states, we will show that $\mathcal{J}[f(t, \cdot)]$ will converge to a solution $J_\infty$ of the compatibility equation (20). However, since this $J_\infty$ (if it is non-zero) is not known from the initial condition (contrary to the case of the BGK equation), it is not easy to control directly the distance between $f$ and $\rho M_{J_\infty}$, but we will see that controlling the distance from $f$ and $\rho M_{\mathcal{J}[f]}$, even if this last one is not a steady-
state, will be the key to our analysis. A convenient framework is to use the relative entropy, for which we will need the following results.

**Proposition 5.1.** Let $\rho > 0$. If $f, g$ are two measurable nonnegative functions on $SO_3(\mathbb{R})$ with total mass $\rho$ and with $g > 0$, we define the relative entropy and Fisher information by

$$ \mathcal{H}(f|g) = \int_{SO_3(\mathbb{R})} f(A) \ln \left( \frac{f(A)}{g(A)} \right) dA, \quad \mathcal{I}(f|g) = \int_{SO_3(\mathbb{R})} f(A) \left\| \nabla \ln \left( \frac{f(A)}{g(A)} \right) \right\|^2 dA. $$

Then, for two such functions, we have the Csiszár-Kullback-Pinsker inequality:

$$ \int_{SO_3(\mathbb{R})} |f(A) - g(A)| dA \le \sqrt{2\rho \mathcal{H}(f|g)}. \quad (38) $$

Finally, we have the following families of (weighted) logarithmic Sobolev inequalities: there exists a constant $\lambda > 0$ such that for all $J \in M_3(\mathbb{R})$ with $\|J\| \le \frac{\sqrt{3}}{\sqrt{2}}\rho$, and all measurable nonnegative function $f$ with total mass $\rho$, we have

$$ \mathcal{H}(f|\rho M_J) \le \frac{1}{2\lambda} \mathcal{I}(f|\rho M_J). \quad (39) $$

*Proof.* The Csiszár-Kullback-Pinsker inequality is well-known [Csi67,Pin64], we just notice the factor $\rho$ since we do not work with probability measures here. The logarithmic Sobolev inequality (39) in the case $J = 0$ (uniform measure on $SO_3(\mathbb{R})$) come for instance from the Bakry-Émery criterion [BE85] since $SO_3(\mathbb{R})$ has positive Ricci curvature (this is the same as the curvature of $\mathbb{S}^3$, thanks to the local isometry $\Phi$ given in Proposition 3.1)$^5$. Then, we use the fact that the logarithmic Sobolev inequality is stable by bounded perturbation [HS87,Vil03]. Since $\|J\| \le \frac{\sqrt{3}}{\sqrt{2}}\rho$, then $M_J$ is bounded above and below, uniformly in $J$, which ends the proof. $\square$

Let us now compute the relative entropy of $f$ with respect to $\rho M_J$ for $J \in M_3(\mathbb{R})$. Using the definition (6), we obtain

$$ \begin{aligned} \mathcal{H}(f|\rho M_J) &= \int_{SO_3(\mathbb{R})} (f(A) \ln f(A) - f(A) A \cdot J) dA + \rho \ln \mathcal{Z}(J) - \rho \ln \rho \\ &= \mathcal{F}[f] + \frac{1}{2} \|J - \mathcal{J}[f]\|^2 - V(J) - \rho \ln \rho, \end{aligned} \quad (40) $$

thanks to the definitions (22) and (34) of $\mathcal{F}$ and $V$. Therefore, if $J_{eq}$ is a solution to the compatibility equation and $f_{eq} = \rho M_{J_{eq}}$, we apply (40) with $f = f_{eq}$ and $J = J_{eq}$ to obtain $\rho \ln \rho = \mathcal{F}[f_{eq}] - V(J_{eq})$. Now applying (40) with $J = J_{eq}$ or with $J = \mathcal{J}[f]$, we obtain

$$ \mathcal{F}[f] - \mathcal{F}[f_{eq}] = \mathcal{H}(f|\rho M_{\mathcal{J}[f]}) + V(\mathcal{J}[f]) - V(J_{eq}), \quad (41) $$

$$ \mathcal{H}(f|f_{eq}) = \mathcal{F}[f] - \mathcal{F}[f_{eq}] + \frac{1}{2}\|J_{eq} - \mathcal{J}[f]\|^2. \quad (42) $$

<footnote>$^5$Actually, as already stated by Bakry and Émery [BE85], this criterion does not give the optimal constant in $\mathbb{S}^3$, which was given by Mueller and Weissler in [MW82], but here even the optimal constant in $\mathbb{S}^3$ would not be necessarily optimal in $SO_3(\mathbb{R})$, since we only want the logarithmic Sobolev inequality for even functions on $\mathbb{S}^3$.</footnote>
Furthermore, it is straightforward to see, thanks to the definition (23) of $\mathcal{D}[f]$, that

$$ \mathcal{D}[f] = \mathcal{I}(f|\rho M_{\mathcal{J}[f]}). \qquad (43) $$

These links between the free energy, its dissipation, the relative entropy, the Fisher information, and the potential $V$ associated to the BGK equation are the key points to prove the stability of the steady-states associated to solutions of the compatibility equation corresponding to local minimizers of $V$.

**Theorem 3.** Let $\rho > \rho^*$ (resp. $\rho < \rho_c$).

We define the set of equilibria $\mathcal{E}_\infty = \{\rho M_{\alpha_1^\dagger(\rho)A_0}, A_0 \in SO_3(\mathbb{R})\}$ (resp. $\mathcal{E}_\infty$ reduced to the uniform distribution on $SO_3(\mathbb{R})$ of mass $\rho$).

Then there exists $\delta > 0$, $\tilde{\lambda} > 0$ and $C > 0$ such that for all nonnegative measurable function $f_0$ with mass $\rho$, if there exists $f_{\text{eq},0} \in \mathcal{E}_\infty$ such that $\mathcal{H}(f_0|f_{\text{eq},0}) < \delta$, then there exists $f_\infty \in \mathcal{E}_\infty$ such that for all time $t \ge 0$, we have

$$ \mathcal{H}(f(t, \cdot)|f_\infty) \le C e^{-2\tilde{\lambda}t} \mathcal{H}(f_0|f_{\text{eq},0}). $$

*Proof.* For convenience, we write $\alpha = \alpha_1^\dagger(\rho)$ (resp. $\alpha = 0$ for the study of stability of the uniform equilibrium) and $V_\infty = V(\alpha I_3)$. We also denote by $E_\infty$ the set of matrices $J_\text{eq}$ solutions to the compatibility equation (20) corresponding to the family of equilibria we are interested in, that is to say $E_\infty = \{\alpha A_0, A_0 \in SO_3(\mathbb{R})\}$.

Since the signature of Hess $\overline{V}(\alpha I_3)$ is $(+++)$ (thanks to Proposition 4.3), by continuity of Hess $\overline{V}$ (and of its smallest eigenvalue), there exists $\delta_0 > 0$ and $\eta > 0$ such that for all diagonal matrix $D$ with $\|D - \alpha I_3\| < \delta_0$, Hess $\overline{V}(D)$ is positive definite with lowest eigenvalue being greater than or equal to $\eta$ (we recall that thanks to (37), its highest eigenvalue is always less than 1). By the following Taylor formulas, for all such $D$, we have

$$ \|\nabla \overline{V}(D)\|^2 = (D - \alpha I_3) \cdot \left( \int_0^1 \text{Hess} \overline{V}(\alpha I_3 + t(D - \alpha I_3)) dt \right)^2 (D - \alpha I_3), \\ V(D) - V_\infty = \int_0^1 (1-t)(D - \alpha I_3) \cdot \text{Hess} \overline{V}(\alpha I_3 + t(D - \alpha I_3))(D - \alpha I_3) dt $$

and therefore

$$ \|\nabla \overline{V}(D)\|^2 \ge \eta \|D - \alpha I_3\|^2, \\ \frac{\eta}{2} \|D - \alpha I_3\|^2 \le V(D) - V_\infty \le \frac{1}{2} \|D - \alpha I_3\|^2 \le \frac{1}{2\eta} \|\nabla \overline{V}(D)\|^2. $$

Therefore, we write $U = \{J \in M_3(\mathbb{R}), \min_{J_{\text{eq}} \in E_\infty} \|J - J_{\text{eq}}\| < \delta_0\}$, which is a neighborhood of $E_\infty$. If $J \in U$ and we write the SSVD $J = PDQ$, we obtain by Proposition 4.1 that $\min_{J_{\text{eq}} \in E_\infty} \|J - J_{\text{eq}}\| = \|D - \alpha I_3\| \le \delta_0$ (when $\alpha > 0$, and the result is still true if $\alpha = 0$ since $E_\infty = \{0\}$ in that case). Therefore, since $V(J) = V(D)$ we obtain that there exists $J_{\text{eq}} \in E_\infty$ (which is equal to $\alpha PQ$) such that

$$ \frac{\eta}{2} \|J - J_{\text{eq}}\|^2 \le V(J) - V_{\infty} \le \frac{1}{2\eta} \|\nabla V(J)\|^2 = \frac{1}{2\eta} \|J - \rho \mathcal{J}[M_J]\|^2. \quad (44) $$
By the Csiszár-Kullback-Pinsker inequality (38), we have that if $g$ is a nonnegative measure with mass $\rho$:

$$ \|\mathcal{J}[f] - \mathcal{J}[g]\| \le \int_{SO_3(\mathbb{R})} \|A\| |f(A) - g(A)| dA \le \frac{\sqrt{3}}{\sqrt{2}} \sqrt{2\rho \mathcal{H}(f|g)}, \quad (45) $$

and therefore for $g = \rho M_{\mathcal{J}[f]}$, we obtain

$$ \|\mathcal{J}[f] - \mathcal{J}[\rho M_{\mathcal{J}[f]}]\| \le \sqrt{3\rho \mathcal{H}(f|\rho M_{\mathcal{J}[f]})}. $$

Combining this with (41) and (44) with $J = \mathcal{J}[f]$, we get that if $\mathcal{J}[f] \in U$, then

$$ \mathcal{F}[f] - \mathcal{F}_{\infty} \le (1 + \frac{3\rho}{2\eta}) \mathcal{H}(f|\rho M_{\mathcal{J}[f]}). $$

Therefore, as soon as $\mathcal{J}[f] \in U$, we have by (43) and the logarithmic Sobolev inequality (39) (we recall that $\|\mathcal{J}[f]\| \le \frac{\sqrt{3}}{\sqrt{2}}\rho$ if the total mass of $f$ is $\rho$):

$$ \mathcal{D}[f] \ge \frac{2\lambda}{1 + \frac{3\rho}{2\eta}} (\mathcal{F}[f] - \mathcal{F}_{\infty}). $$

By the dissipation of the free energy (21), writing $\tilde{\lambda} = \frac{\lambda}{1 + \frac{3\rho}{2\eta}}$ we obtain that as long as $\mathcal{J}[f] \in U$,

$$ 0 \le \mathcal{F}[f] - \mathcal{F}_{\infty} \le e^{-2\tilde{\lambda}t} (\mathcal{F}[f_0] - \mathcal{F}_{\infty}) \le e^{-2\tilde{\lambda}t} \mathcal{H}(f_0|f_{\text{eq},0}), \quad (46) $$

the first inequality coming from (41) and the fact that $V(\mathcal{J}[f]) - V_{\infty} \ge 0$ thanks to (44), and the last inequality coming from (42). Finally, thanks to (44), (41) and (46), we obtain that still as long as $\mathcal{J}[f(t, \cdot)] \in U$, there exists $J_{\text{eq}}(t)$ such that

$$ \|\mathcal{J}[f(t, \cdot)] - J_{\text{eq}}(t)\| \le \sqrt{\frac{2}{\eta}(V(\mathcal{J}[f(t)]) - V_{\infty})} \le \frac{\sqrt{2}}{\sqrt{\eta}} e^{-\tilde{\lambda}t} \sqrt{\mathcal{H}(f_0|f_{\text{eq},0})}. \quad (47) $$

Therefore, by taking $\delta = \min(\frac{\eta}{2}\delta_0^2, \frac{1}{3\rho}\delta_0^2)$, and using (45) with $g = f_{\text{eq},0}$, we obtain that if $\mathcal{H}(f_0|f_{\text{eq},0}) < \delta$, then $\|\mathcal{J}[f_0] - \mathcal{J}[f_{\text{eq},0}]\| < \delta_0$, so $\mathcal{J}[f_0] \in U$, and for all positive time $\|\mathcal{J}[f(t, \cdot)] - J_{\text{eq}}(t)\| < \delta_0$ (and therefore $\mathcal{J}[f(t, \cdot)]$ stays in $U$) thanks to (47). Indeed, if it was not the case, for the first exit time $t_0 > 0$ of $U$, we would have $\|\mathcal{J}[f(t_0, \cdot)] - J_{\text{eq}}(t_0)\| \le \frac{\sqrt{2}}{\sqrt{\eta}} \sqrt{\mathcal{H}(f_0|f_{\text{eq},0})} < \delta_0$ which is a contradiction. From now on we suppose that $\mathcal{H}(f_0|f_{\text{eq},0}) < \delta$, so that (47) and (46) are valid for all time $t \ge 0$.

Let us now find a way to control the displacement of $\mathcal{J}[f]$. For $J \in M_3(\mathbb{R})$, using the Fokker-Planck equation (19) and integrating by parts, we have

$$ \frac{d}{dt} J \cdot \mathcal{J}[f] = \int_{SO_3(\mathbb{R})} [\nabla_A(A \cdot J) \cdot \nabla_A(A \cdot \mathcal{J}[f]) - \Delta_A(A \cdot J)] f(A) dA, $$

which can be written

$$ \frac{d}{dt} \mathcal{J}[f] = \mathcal{M}[f](\mathcal{J}[f]) - \mathcal{L}[f], \quad (48) $$
where, when $g$ is an integrable function on $SO_3(\mathbb{R})$, we define $\mathcal{M}[g]$ as the linear operator from $M_3(\mathbb{R})$ to $M_3(\mathbb{R})$ given by the fact that for any $J, J' \in M_3(\mathbb{R})$,

$$J \cdot \mathcal{M}[g](J') = \int_{SO_3(\mathbb{R})} \nabla_A(A \cdot J) \cdot \nabla_A(A \cdot J')g(A) \, dA,$$

and $\mathcal{L}[g]$ as the matrix<fnref n="6" /> such that for all $J \in M_3(\mathbb{R})$,

$$J \cdot \mathcal{L}[g] = \int_{SO_3(\mathbb{R})} \Delta_A(A \cdot J)g(A) \, dA.$$

We therefore see that since the functions under the integral are smooth and bounded, there exists $C_0 > 0$ such that for all $J \in M_3(\mathbb{R})$ and for any integrable function $g$ on $SO_3(\mathbb{R})$,

$$\|\mathcal{L}[g]\| \le C_0 \int_{SO_3(\mathbb{R})} |g(A)| \, dA \quad \text{and} \quad \|\mathcal{M}[g](J)\| \le C_0 \|J\| \int_{SO_3(\mathbb{R})} |g(A)| \, dA. \quad (49)$$

Therefore, defining $f_{\text{eq}}(t, \cdot) = \rho M_{J_{\text{eq}}(t)}$, and using the fact that it is a stationary solution, thus giving by (48) that $\mathcal{M}[f_{\text{eq}}](\mathcal{J}[f_{\text{eq}}]) - \mathcal{L}[f_{\text{eq}}] = 0$, we obtain

$$\begin{aligned} \left\| \frac{d}{dt} \mathcal{J}[f] \right\| &= \| \mathcal{M}[f](\mathcal{J}[f]) - \mathcal{M}[f_{\text{eq}}](\mathcal{J}[f_{\text{eq}}]) - \mathcal{L}[f] + \mathcal{L}[f_{\text{eq}}] \| \\ &\le \| \mathcal{M}[f](\mathcal{J}[f - f_{\text{eq}}]) \| + \| \mathcal{M}[f - f_{\text{eq}}](\mathcal{J}[f_{\text{eq}}]) \| + \| \mathcal{L}[f - f_{\text{eq}}] \|. \end{aligned}$$

Therefore, by using (49) and the Csiszár-Kullback-Pinsker inequalities (38) and (45), we get that there exists a constant $C_1 > 0$ (only depending on $\rho$) such that

$$\left\| \frac{d}{dt} \mathcal{J}[f] \right\| \le \sqrt{C_1 \mathcal{H}(f|f_{\text{eq}})}.$$

Combining this with (42), (47), and (46), we then get that there exists a constant $C_2$ (not depending on $f_0$) such that for all $t \ge 0$

$$\left\| \frac{d}{dt} \mathcal{J}[f] \right\| \le e^{-\tilde{\lambda}t} \sqrt{C_2 \mathcal{H}(f_0|f_{\text{eq},0})}.$$

Finally, this gives that $\mathcal{J}[f]$ converges exponentially fast with rate $\tilde{\lambda}$ towards a given matrix $J_{\infty} \in M_3(\mathbb{R})$ and since the distance between $\mathcal{J}[f]$ and $E_{\infty}$ converges to 0 thanks to (47), we obtain that $J_{\infty} \in E_{\infty}$. More precisely, we have

$$\| \mathcal{J}[f(t, \cdot)] - J_{\infty} \| \le \int_{t}^{+\infty} \left\| \frac{d}{ds} \mathcal{J}[f(s, \cdot)] \right\| dt \le \frac{e^{-\tilde{\lambda}t}}{\tilde{\lambda}} \sqrt{C_2 \mathcal{H}(f_0|f_{\text{eq},0})}. \quad (50)$$

Defining $f_{\infty} = \rho M_{J_{\infty}}$ and using (42) with $f_{\text{eq}} = f_{\infty}$, (46), and (50), we then get that there exists a constant $C_3 > 0$ (not depending on $f_0$) such that

$$\mathcal{H}(f(t, \cdot)|f_{\infty}) \le C_3 e^{-2\tilde{\lambda}t} \mathcal{H}(f_0|f_{\text{eq},0}),$$

which ends the proof. Let us remark that this proof covers the case $\alpha = 0$, but if we only want to do this case, it can be simplified a lot since $E_{\infty} = \{0\}$. □

<footnote n="6">⁶We can actually show (but we do not need it here) that $\mathcal{L}[f]$ is proportional to $\mathcal{J}[f]$. Indeed, since $q \mapsto q \cdot Qq$ is an eigenfunction of the Laplacian on the unit sphere of $\mathbb{R}^4$ (more precisely a spherical harmonic of degree 2) when $Q$ is a symmetric trace-free matrix, we get, thanks to the local isometry $\Phi$ and Proposition (3.2), that $A \mapsto A \cdot J$ is also an eigenfunction of the Laplacian on $SO_3(\mathbb{R})$.</footnote>
Let us finish this section by some comments. The proof of Theorem 3 has been
done here in relative entropy. It may look similar in some points to [FKM18], but
the main idea is above all based on the fact that we measure the relative entropy
with respect to a target measure $\rho M_{J[f]}$ which is not itself a steady-state. The
fine control of the potential $V$ around the solutions of the compatibility equation
is the key to link all these different quantities. The proof would have worked the
same in $L^2$, by using the regularizing effect of the equation (and $L^\infty$ bounds), as
was done in [FL12] for the Vicsek model, but the main difference is again that
we would compare $D[f]$ and $F[f] - F_\infty$ with $\|f - \rho M_{J[f]}\|_2^2$. This proof seems to
be adaptable to a lot of different models of Fokker-Planck type, such as the Doi-
Onsager theory for suspensions of rodlike polymers, for which, as far as we know,
no proof of exponential convergence is available (but the analog to the potential $V$
has been studied, therefore the nature of the critical points is well-known). This is
left for future work.

Finally, now that we have a good understanding of the long time behaviour of
the Fokker-Planck equation (19), we could try to further understand the limit of the
particle system as $N \to \infty$. Since the mean-field limit is essentially a law of large
numbers, we expect fluctuations of order $\frac{1}{\sqrt{N}}$, which explains why the order parame-
ters of the numerical simulations in Figure 4 are not so close to 0 for what is expected
to be the uniform distribution. More precisely, as indicated by the estimate (18),
the distance between the empirical measure and the solution to the Fokker-Planck
equation can be bounded by $\frac{e^{CT}}{\sqrt{N}}$, for all $t$ in $[0, T]$. Therefore if we want such an
estimate for a large time $T$, we cannot do better than $T$ of order $\ln N$. However,
since the equilibria are exponentially stable, the fluctuations that would push the
empirical distribution away from the family of stable equilibria, are compensated
by the deterministic dynamics of the Fokker-Planck equation. Therefore the only
remaining fluctuations would cause the solution to fluctuate mainly in the tangential
component of the family of equilibria. This approach has been made rigorous in the
case of identical Kuramoto oscillators in [BGP14] (which corresponds to the Vicsek
model studied in [FL12] in dimension two), where it is proved that the solution
stays close to the set of equilibria up to times of order $N$, but with the center of
synchronization of the distribution performing a Brownian motion on the circle at
these time scales. In analogy with this result, we could expect in our case that, close
to the family of von Mises distributions $\rho M_{\alpha A}$ with $\alpha > 0$ and $A \in SO_3(\mathbb{R})$, the long
time behaviour at time $t = sN$ of the empirical measure of particle system would
be close to $\rho M_{\alpha A(s)}$, where $A(s)$ performs a Brownian motion on $SO_3(\mathbb{R})$. This is
also left for future work.

Acknowledgements

The author wants to thank his collaborators Pierre Degond, Sara Merino-Aceituno,
Ariane Trescases and Antoine Diez for all the work done together on body-attitude
coordination models [DFMA17,DFMAT18,DFMAT19,DDFMA20,Die20], which in-
spired the talk given at $[iN\delta A\mathbb{M}]$ in November 2019, and finally led to the present
paper. This work has been supported by the Project EFI ANR-17-CE40-0030 of the
References

[BCC12] F. Bolley, J. A. Cañizo, and J. A. Carrillo. Mean-field limit for the stochastic Vicsek model. *Appl. Math. Lett.*, 3(25):339–343, 2012.

[BE85] D. Bakry and Michel Émery. Diffusions hypercontractives. In *Séminaire de probabilités, XIX, 1983/84*, volume 1123 of *Lecture Notes in Math.*, pages 177–206. Springer, Berlin, 1985.

[BGP14] Lorenzo Bertini, Giambattista Giacomin, and Christophe Poquet. Synchronization and random long time dynamics for mean-field plane rotators. *Probab. Theory Related Fields*, 160(3-4):593–653, 2014.

[CGGR08] H. Chaté, F. Ginelli, G. Grégoire, and F. Raynaud. Collective motion of self-propelled particles interacting without cohesion. *Phys. Rev. E*, 77(4):046113, 2008.

[CKT04] P. Constantin, I. G. Kevrekidis, and E. S. Titi. Asymptotic states of a Smoluchowski equation. *Arch. Ration. Mech. Anal.*, 174(3):365–384, 2004.

[Csi67] I. Csiszár. Information-type measures of difference of probability distributions and indirect observations. *Studia Sci. Math. Hungar.*, 2:299–318, 1967.

[DDFMA20] P. Degond, A. Diez, A. Frouvelle, and S. Merino-Aceituno. Phase transitions and macroscopic limits in a BGK model of body attitude coordination. *J. Nonlinear Sci.*, 2020.

[DFL13] Pierre Degond, Amic Frouvelle, and Jian-Guo Liu. Macroscopic limits and phase transition in a system of self-propelled particles. *J. Nonlinear Sci.*, 23(3):427–456, 2013.

[DFL15] P. Degond, A. Frouvelle, and J.-G. Liu. Phase transitions, hysteresis, and hyperbolicity for self-organized alignment dynamics. *Arch. Ration. Mech. Anal.*, 216(1):63–115, 2015.

[DFMA17] P. Degond, A. Frouvelle, and S. Merino-Aceituno. A new flocking model through body attitude coordination. *Math. Models Methods Appl. Sci.*, 27(06):1005–1049, 2017.

[DFMAT18] P. Degond, A. Frouvelle, S. Merino-Aceituno, and A. Trescases. Quaternions in collective dynamics. *Multiscale Mod. Simul.*, 16(1):28–77, 2018.
[DFMAT19] P. Degond, A. Frouvelle, S. Merino-Aceituno, and A. Trescases. Alignment of self-propelled rigid bodies: from particle systems to macroscopic equations. In G. Giacomin, S. Olla, E. Saada, H. Spohn, and Stoltz G., editors, *Stochastic Dynamics Out of Equilibrium. IHPStochDyn 2017*, volume 282 of Springer proceedings in mathematics and statistics, pages 28–66. Springer, 2019.

[Die20] Antoine Diez. Propagation of chaos and moderate interaction for a piecewise deterministic system of geometrically enriched particles. *Electron. J. Probab.*, 25:Paper No. 90, 38, 2020.

[DM08] P. Degond and S. Motsch. Continuum limit of self-driven particles with orientation interaction. *Math. Models Methods Appl. Sci.*, 18:1193–1215, 2008.

[FKM18] Alessio Figalli, Moon-Jin Kang, and Javier Morales. Global well-posedness of the spatially homogeneous Kolmogorov-Vicsek model as a gradient flow. *Arch. Ration. Mech. Anal.*, 227(3):869–896, 2018.

[FL12] A. Frouvelle and J.-G. Liu. Dynamics in a kinetic model of oriented particles with phase transition. *SIAM J. Math. Anal.*, 44(2):791–826, 2012.

[FS05] I. Fatkullin and V. Slastikov. Critical points of the Onsager functional on a sphere. *Nonlinearity*, 18:2565–2580, 2005.

[Hor54] Alfred Horn. Doubly stochastic matrices and the diagonal of a rotation matrix. *Amer. J. Math.*, 76:620–630, 1954.

[HS87] Richard Holley and Daniel Stroock. Logarithmic Sobolev inequalities and stochastic Ising models. *J. Statist. Phys.*, 46(5-6):1159–1194, 1987.

[Hsu02] E. P. Hsu. *Stochastic Analysis on Manifolds*, volume 38 of *Graduate Series in Mathematics*. American Mathematical Society, Providence, Rhode Island, 2002.

[LZZ05] H. Liu, H. Zhang, and P. Zhang. Axial symmetry and classification of stationary solutions of Doi-Onsager equation on the sphere with Maier-Saupe potential. *Commun. Math. Sci.*, 3(2):201–218, 2005.

[MW82] Carl E. Mueller and Fred B. Weissler. Hypercontractivity for the heat semigroup for ultraspherical polynomials and on the $n$-sphere. *J. Functional Analysis*, 48(2):252–283, 1982.

[Pin64] M. S. Pinsker. *Information and information stability of random variables and processes*. Translated and edited by Amiel Feinstein. Holden-Day, Inc., San Francisco, Calif.-London-Amsterdam, 1964.

[Szn91] A.-S. Sznitman. Topics in propagation of chaos. In *École d'Été de Probabilités de Saint-Flour XIX — 1989*, volume 1464 of Lecture Notes in Mathematics, pages 165–251. Springer, Berlin, 1991.
[VCBJ+95] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet. Novel type of phase transition in a system of self-driven particles. *Phys. Rev. Lett.*, 75(6):1226–1229, 1995.

[Vil03] Cédric Villani. *Topics in optimal transportation*, volume 58 of *Graduate Studies in Mathematics*. American Mathematical Society, Providence, RI, 2003.

[WH08] H. Wang and P. J. Hoffman. A unified view on the rotational symmetry of equilibria of nematic polymers, dipolar nematic polymers and polymers in higher dimensional space. *Commun. Math. Sci.*, 6(4):949–974, 2008.