Pierre Degond
email: pierre.degond@math.univ-toulouse.fr

Antoine Diez
email: toine.diez18@imperial.ac.uk

Mingye Na
email: mingye.na18@imperial.ac.uk

Bulk topological states in a new collective dynamics model *

Keywords: individual-based model, macroscopic model, self-organization, topological phase transition, winding number, order parameter AMS subject classifications. 22E70, 35Q70, 37B25, 60J76, 65C35, 70F10

In this paper, we demonstrate the existence of topological states in a new collective dynamics model. This individual-based model (IBM) describes self-propelled rigid bodies moving with constant speed and adjusting their rigid-body attitude to that of their neighbors. In previous works, a macroscopic model has been derived from this IBM in a suitable scaling limit. In the present work, we exhibit explicit solutions of the macroscopic model characterized by a non-trivial topology. We show that these solutions are well approximated by the IBM during a certain time but then the IBM transitions towards topologically trivial states. Using a set of appropriately defined topological indicators, we reveal that the breakage of the non-trivial topology requires the system to go through a phase of maximal disorder. We also show that similar but topologically trivial initial conditions result in markedly different dynamics, suggesting that topology plays a key role in the dynamics of this system.

1. Introduction. Systems of particles (or agents) which exhibit self-organized collective behavior are ubiquitous in the living world at all scales, from bird flocks [START_REF] Lukeman | Inferring individual rules from collective behavior[END_REF] to sperm [START_REF] Creppy | Symmetrybreaking phase transitions in highly concentrated semen[END_REF] or bacterial colonies [START_REF] Czirók | Formation of complex bacterial colonies via self-generated vortices[END_REF]. Examples are also found in social sciences [START_REF] Castellano | Statistical physics of social dynamics[END_REF][START_REF] Degond | Continuum dynamics of the intention field under weakly cohesive social interaction[END_REF] or for inert matter [START_REF] Bricard | Emergent vortices in populations of colloidal rollers[END_REF]. In such systems, the agents interact locally with a limited number of neighbors through rather simple rules such as attraction, repulsion or alignment [3,[START_REF] Couzin | Collective memory and spatial sorting in animal groups[END_REF][START_REF] Gautrais | Deciphering interactions in moving animal groups[END_REF] without any leader or centralized control. When the number of agents becomes large, vast structures encompassing many agents appear, such as clusters [START_REF] Martín-Gómez | Collective motion of active brownian particles with polar alignment[END_REF][START_REF] Vicsek | Collective motion[END_REF], traveling bands [START_REF] Chaté | Collective motion of self-propelled particles interacting without cohesion[END_REF], vortices [START_REF] Costanzo | Spontaneous emergence of milling (vortex state) in a Vicsek-like model[END_REF][START_REF] Czirók | Formation of complex bacterial colonies via self-generated vortices[END_REF], lanes [START_REF] Couzin | Self-organized lane formation and optimized traffic flow in army ants[END_REF], etc. As there is no direct or apparent relation between these structures and the nature of the agents interactions, such a phenomenon is named "emergence". Its study has stimulated a vast literature (see e.g. [START_REF] Vicsek | Collective motion[END_REF] for a review).

There are mainly two levels of description of particle systems: the most detailed one consists of individual based models (IBM) where the agents dynamics are described by coupled ordinary or stochastic differential equations. When the number of agents becomes large, a macroscopic description in terms of average quantities such as the agents mean density or velocity is preferred. The rigorous link between these two levels of description involves two successive limits by which the number of agents is first sent to infinity (mean-field limit) and then, the system size relative to the typical interaction distance between the agents is also sent to infinity (hydrodynamic limit), see e.g. [START_REF] Cercignani | The Mathematical Theory of Dilute Gases[END_REF][START_REF] Degond | Macroscopic limits of the Boltzmann equation: a review[END_REF]. In collective dynamics, particles are capable of self-propulsion by transforming an internal source of chemical energy into motion [START_REF] Vicsek | Collective motion[END_REF]. There are two main classes of IBM of self-propelled particles. The first class is based on the Cucker-Smale model [4,[START_REF] Cucker | Emergent behavior in flocks[END_REF][START_REF] Ha | A simple proof of the Cucker-Smale flocking dynamics and mean-field limit[END_REF][START_REF] Ha | From particle to kinetic and hydrodynamic descriptions of flocking[END_REF] where self-propulsion is treated as an external force. The second class is based on the Vicsek model [2,[START_REF] Caussin | Emergent spatial structures in flocking models: a dynamical system insight[END_REF][START_REF] Chaté | Collective motion of self-propelled particles interacting without cohesion[END_REF][START_REF] Czirók | Formation of complex bacterial colonies via self-generated vortices[END_REF][START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF][START_REF] Martín-Gómez | Collective motion of active brownian particles with polar alignment[END_REF][START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] where self-propulsion is modeled by imposing the norm of the particle velocity to be a constant. At the mean-field or hydrodynamic levels, the two frameworks give rise to corresponding models (see e.g. [1,5] for Cucker-Smale type models and [10,[START_REF] Degond | Phase transitions, hysteresis, and hyperbolicity for selforganized alignment dynamics[END_REF][START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF][START_REF] Peruani | A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms[END_REF][START_REF] Toner | Flocks, herds, and schools: A quantitative theory of flocking[END_REF] for Vicsek type models). The two categories are linked by an asymptotic limit [12,[START_REF] Bostan | Reduced fluid models for self-propelled particles interacting through alignment[END_REF]. Of course, there are many variants of these models and we refer to [8,9,[START_REF] Carrillo | Double milling in self-propelled swarms from kinetic theory[END_REF][START_REF] Cavagna | Flocking and turning: a new model for self-organized collective motion[END_REF][START_REF] Degond | A macroscopic model for a system of swarming agents using curvature control[END_REF][START_REF] D'orsogna | Self-propelled particles with soft-core interactions: patterns, stability, and collapse[END_REF][START_REF] Motsch | A new model for self-organized dynamics and its flocking behavior[END_REF] for a non-exhaustive set of examples.

Recently, a series of studies has investigated the existence of topological states in collective dynamics. Topological states have appeared with the quantum Hall effect [START_REF] Klitzing | New Method for High-Accuracy Determination of the Fine-Structure Constant Based on Quantized Hall Resistance[END_REF][START_REF] Laughlin | Quantized Hall conductivity in two dimensions[END_REF][START_REF] Foundation | The Nobel Prize in Physics[END_REF][START_REF] Thouless | Quantization of particle transport[END_REF] which relies on so-called conducting chiral edge states: when a sample of a 2-dimensional insulator is placed in a magnetic field, its bulk conductance is nil but a current can flow around its edges in only one direction (hence the 'chiral' terminology). Then, materials that exhibit chiral edge states without a magnetic field have been discovered, the so-called "topological insulators" [START_REF] Hasan | Colloquium: topological insulators[END_REF][START_REF] Foundation | The Nobel Prize in Physics[END_REF][START_REF] Qi | Topological insulators and superconductors[END_REF]. Chiral edge states are robust against perturbations because of their non trivial topology which can be characterized by a integer, the winding number. Any destruction of the chiral edge state would require a finite jump of this integer, which consumes a finite amount of energy. Hence lower energy perturbations will fail to destroy the chiral edge state. This property is of strategic interest for various applications such as quantum computers. Recently a series of works have explored the occurrence of topological states in collective dynamics (see e.g. [START_REF] Shankar | Topological sound and flocking on curved surfaces[END_REF][START_REF] Sone | Anomalous topological active matter[END_REF][START_REF] Souslov | Topological sound in active-liquid metamaterials[END_REF]). They are based on numerical simulations of the Toner and Tu model [START_REF] Toner | Flocks, herds, and schools: A quantitative theory of flocking[END_REF], which is a continuum analog of the Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF]. Investigating appropriate geometrical configurations (a sphere in [START_REF] Shankar | Topological sound and flocking on curved surfaces[END_REF], a network of rings in [START_REF] Sone | Anomalous topological active matter[END_REF][START_REF] Souslov | Topological sound in active-liquid metamaterials[END_REF]), they show that linearized perturbations of the stationary state (i.e. sound waves) generate chiral edge states which propagate uni-directionally, revealing an underpinning non-trivial topology. However, the question of whether this effect could be realized with a finite (even large) number of discrete particles and whether the topological states would survive the noise induced by this finite particle number long enough is not investigated.

In this paper, we demonstrate the existence of non-trivial bulk topological states in a new collective dynamics model. Bulk states propagate in the whole domain, by opposition to edge states which are localized at the boundary. The collective dynamics model studied here has first been proposed in [START_REF] Degond | A new flocking model through body attitude coordination[END_REF] and later analyzed and expanded in [START_REF] Degond | Phase transitions and macroscopic limits in a BGK model of body-attitude coordination[END_REF][START_REF] Degond | Quaternions in collective dynamics[END_REF][START_REF] Degond | Alignment of self-propelled rigid bodies: from particle systems to macroscopic equations[END_REF]. Referred to below as the "Body-Alignment Individual-Based Model" (BA-IBM or IBM for short), it describes self-propelled rigid bodies moving with constant speed and trying to adjust their rigid body attitude to that of their neighbors. In [START_REF] Degond | Quaternions in collective dynamics[END_REF][START_REF] Degond | A new flocking model through body attitude coordination[END_REF] the BA-IBM was based on Stochastic Differential Equations (SDE) and a macroscopic model named the "Self-Organized Hydrodynamics for Body-orientation (SOHB)" was derived. In [START_REF] Degond | Alignment of self-propelled rigid bodies: from particle systems to macroscopic equations[END_REF][START_REF] Degond | Phase transitions and macroscopic limits in a BGK model of body-attitude coordination[END_REF], SDE were replaced by Piecewise Deterministic Markov Processes (PDMP) in the IBM but the macroscopic model remained the SOHB model (with possibly different coefficients). In [START_REF] Degond | Phase transitions and macroscopic limits in a BGK model of body-attitude coordination[END_REF], a variant of the BA-IBM was shown This manuscript is for review purposes only.

to exhibit phase transitions which were rigorously studied. In the present work, we derive explicit solutions of the SOHB model which exhibit striking non-trivial topologies revealed by non-zero winding numbers. We explore how these non-trivial topologies are maintained at the level of the IBM by solving the PDMP of [START_REF] Degond | Alignment of self-propelled rigid bodies: from particle systems to macroscopic equations[END_REF]. In particular, we observe that, due to noise induced by the finite particle number, topological phase transitions from states with nontrivial topology to states with trivial one may occur and we study these phase transitions in detail. Using a set of appropriately defined topological indicators, we reveal that the breakage of the non-trivial topology requires the system to go through a phase of maximal disorder.

We also show that similar but topologically trivial initial conditions result in markedly different dynamics, suggesting that topology plays a key role in the dynamics of this system.

We are led to question the possible existence of topological protection against perturbations as mentioned above for topological insulators. Compared to previous works on topological states in collective dynamics, we deal with bulk states instead of edge states and we explore them at the level of the IBM and not just at the continuum level, which is closer to realistic particle systems. The present work adds a new item to the list of collective dynamics models exhibiting topological states. The topological protection concept could bring new perspectives to poorly understood questions such as the robustness of morphogenesis or the emergence of symmetries in growing organisms.

The present model belongs to the category of Vicsek-like models in the sense that it introduces a geometrical constraint within the degrees of freedom of the particles. In the Vicsek model, the particle velocities were constrained to belong to the unit sphere (after convenient normalization). In the present IBM, the particles carry an orthonormal frame, or equivalently, a rotation matrix, that describes their body attitude. Thus their degrees of freedom are constrained to belong to the manifold SO 3 (R) of 3 × 3 rotation matrices. Fig. 1 highlights the difference between the Vicsek and body orientation models. The left picture shows alignment of two agents in the Vicsek sense, while the right picture shows alignment in the body-alignment sense. We mention that models involving full body attitudes have already been considered in [START_REF] Cavagna | Flocking and turning: a new model for self-organized collective motion[END_REF][START_REF] Hemelrijk | Schools of fish and flocks of birds: their shape and internal structure by self-organization[END_REF][START_REF] Hemelrijk | Emergence of oblong school shape: models and empirical data of fish[END_REF][START_REF] Hildenbrandt | Self-organized aerial displays of thousands of starlings: a model[END_REF] in the context of flocking, but the alignment rules were different and essentially based on a velocity orientation (and not full body attitude) alignment.

not aligned aligned Polar alignment

Body-orientation not aligned aligned

Figure 1: Vicsek model versus body-alignment model. Left: polar alignment of velocity orientations (red vectors) of two agents. Right: alignment of body-orientations: in addition to its velocity orientation (red), each agent has two other axes (green and blue), the three vectors forming a direct orthogonal frame.

We complete this introduction by a review of the mathematical literature on the Vicsek model and the BA-IBM. The mean-field limit of the IBM has been proven in [10] for the Vicsek model and in [START_REF] Diez | Propagation of chaos and moderate interaction for a piecewise deterministic system of geometrically enriched particles[END_REF] for the body orientation model. Existence theory for the mean-field Vicsek model is available in [START_REF] Briant | Cauchy theory and mean-field limit for general Vicsek models in collective dynamics[END_REF][START_REF] Figalli | Global well-posedness of the spatially homogeneous Kolmogorov-Vicsek model as a gradient flow[END_REF][START_REF] Gamba | Global weak solutions for Kolmogorov-Vicsek type equations with orientational interactions[END_REF] but the corresponding theory for the mean-field body orientation model is still open. The mean-field kinetic models exhibit phase transitions which have been studied in [START_REF] Degond | Macroscopic limits and phase transition in a system of self-propelled particles[END_REF][START_REF] Degond | Phase transitions, hysteresis, and hyperbolicity for selforganized alignment dynamics[END_REF][START_REF] Frouvelle | Dynamics in a kinetic model of oriented particles with phase transition[END_REF] and [START_REF] Degond | Phase transitions and macroscopic limits in a BGK model of body-attitude coordination[END_REF] for the Vicsek and body orientation models respectively.

The numerical approximation of the mean-field kinetic model has been undertaken for the Vicsek model only in [START_REF] Gamba | Spectral method for a kinetic swarming model[END_REF][START_REF] Griette | Kinetic equations and self-organized band formations[END_REF]. The derivation of macroscopic equations from the mean-field

Vicsek kinetic equations has first been formally achieved in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] and later rigorously proved in [START_REF] Jiang | Hydrodynamic limits of the kinetic self-organized models[END_REF]. Corresponding works for the body alignment model are only formal [START_REF] Degond | A new flocking model through body attitude coordination[END_REF][START_REF] Degond | Quaternions in collective dynamics[END_REF][START_REF] Degond | Alignment of self-propelled rigid bodies: from particle systems to macroscopic equations[END_REF]. Existence theory for the hydrodynamic models derived from the Vicsek model can be found in [START_REF] Degond | Hydrodynamic models of self-organized dynamics: derivation and existence theory[END_REF][START_REF] Zhang | A local existence of viscous self-organized hydrodynamic model[END_REF] and numerical methods in [START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF][START_REF] Gamba | Spectral method for a kinetic swarming model[END_REF][START_REF] Motsch | Numerical simulations of a nonconservative hyperbolic system with geometric constraints describing swarming behavior[END_REF]. Both questions are still open for the body orientation model.

The organization of this paper is as follows. Section 2 is devoted to the exposition of the IBM and macroscopic models. Then explicit solutions of the macroscopic model are derived in Section 3 and are shown to exhibit non-trivial topology. They also serve as benchmarks to

show that the macroscopic model is an accurate approximation of the IBM. But after a some time, the IBM departs from the special solutions of the macroscopic model and undergoes a topological phase transition. The study of these phase transitions require appropriate topological indicators which are developed in Section 4. Then, the topological phase transitions are analyzed in Section 5. A discussion and some open questions raised by these observations can be found in Section 6. The supplementary material (SM) collects additional information: a list of supplementary videos (Section SM1), a summary of the quaternion framework (Section SM2), a description of the numerical methods (Section SM3), a summary of the derivation of the macroscopic models (Section SM4) and finally a derivation of the explicit solutions presented in Section 3 (Section SM6).

Models.

2.1. The Individual-Based body-alignment Model.

2.1.1. Description of the model. In this section, we present the Individual-Based bodyalignment Model (IBM). This model was first proposed in [START_REF] Degond | Alignment of self-propelled rigid bodies: from particle systems to macroscopic equations[END_REF]. We consider N particles (or individuals, or agents) indexed by k ∈ {1, . . . , N } whose spatial locations are denoted by

X k (t) ∈ R 3 where t ∈ [0, ∞) is the time. A direct orthonormal frame {Ω k (t), u k (t), v k (t)} is attached to each particle (i.e. Ω k , u k , v k ∈ S 2 , Ω k • u k = 0 and v k = Ω k × u k ). Likewise, if (e 1 , e 2 , e 3 )
is a fixed direct orthonormal reference frame, we define A k (t) to be the unique element of the special orthonormal group SO 3 (R) which maps (e 1 , e 2 , e 3 ) onto (Ω k (t), u k (t), v k (t)). We will choose (e 1 , e 2 , e 3 ) once for all and write

A k (t) = [Ω k (t), u k (t), v k (t)
]. This will be referred to as the local particle frame or as the particle's body orientation. Ω k (t) is the self-propulsion direction: Particle k moves in straight line in the direction of Ω k with unchanged local frame

A k except at exponentially distributed times at which the local frame jumps and adjusts itself to the average neighbors' local frame up to some noise. The motion of the particles is thus

described by the functions [0, ∞) ∋ t → (X k (t), A k (t)) ∈ R 3 × SO 3 (R) for k ∈ {1, . . . , N }.
We first describe how the average neighbors' local frame is defined. We introduce a fixed

This manuscript is for review purposes only.

observation (or sensing) kernel K: R 3 ∋ x → K(x) ∈ [0, ∞). We assume that K is a radial function (i.e. there exists K

: [0, ∞) ∋ r → K(r) ∈ [0, ∞) such that K(x) = K(|x|),
where |x| is the euclidean norm of x). For a collection of N particles {(X k , A k )} k∈{1,...,N } ∈ (R 3 × SO 3 (R)) N , we define the local flux as the following 3 × 3 matrix:

J k = 1 N N j=1 K(X k -X j ) A j .
Typically, we can think of K(x) as the indicator function of the ball centered at zero with radius R. In this case, J k is just the sum of the matrices A j of all particles j located within a distance R to Particle k, divided by the total number of particles N . However, more sophisticated sensing functions can be used to account for the fact that e.g. distant particles will contribute to J k less than neighboring particles. In general, J k is not a rotation matrix.

To recover a rotation matrix, we need to map J k back onto the manifold SO 3 (R). To do so, the space M 3 (R) of 3 × 3 matrices, is equipped with the inner product:

(2.1)

A • B := 1 2 Tr(A T B),
where Tr denotes the trace operator and A T is the transpose of the matrix A. Now, we define the average neighbors' local frame A k of Particle k as follows:

(2.2)

A k := arg max A∈SO 3 (R) A • J k .

This expression stands for the element

A k ∈ SO 3 (R) that maximizes the function SO 3 (R) ∋ A → A • J k ∈ R. The maximization procedure (2.
2) has a unique solution as soon as J k is not singular, i.e. det J k ̸ = 0 where det stands for the determinant. Since the singular matrices form a zero-measure set in M 3 (R) it is legitimate to assume that, except for a zero-measure set of initial data, this situation will not occur. Furthermore, when det J k > 0, A k is nothing but the unique rotation matrix involved in the polar decomposition of J k .

We let the particles evolve according to the following Piecewise Deterministic Markov Process (PDMP).

• To each agent k ∈ {1, . . . , N } is attached an increasing sequence of random times (jump times) T 1 k , T 2 k , . . . such that the intervals between two successive times are independent and follow an exponential law with constant parameter ν > 0 (Poisson process). At each jump time T n k , the function X k is continuous and the function A k has a discontinuity between its left and right states respectively denoted by

A k (T n k -0) and A k (T n k + 0). • Between two jump times (T n k , T n+1 k
), the evolution is deterministic: the orientation of Agent k does not change and it moves in straight line at speed c 0 > 0 in the direction

A k (T n k + 0) e 1 , i.e. for all t ∈ [T n k , T n+1 k ), we have (2.3) X k (t) = X k (T n k ) + c 0 (t -T n k ) A k (t) e 1 , A k (t) = A k (T n k + 0).
This manuscript is for review purposes only.

• To compute A k (T n k + 0) from A k (T n k -0), we compute the local flux defined at time

T n k -0 given by:

(2.4) (2.5)

J n- k := 1 N N j=1 K X k (T n k ) -X j (T n k ) A j (T n k -0), having in mind that A j (T n k -0) = A j (T n k ) for j ̸ = k. From J n- k ,
A k (T n k + 0) ∼ M A n k .
The von Mises distribution on SO 3 (R) with parameter A ∈ SO 3 (R) is defined to be the probability density function:

(2.6) M A (A) := e κA•A SO 3 (R) e κA•A ′ dA ′ ,
where κ > 0 is a supposed given parameter named concentration parameter, or inverse of the noise intensity. The von Mises distribution, also known in the literature as the matrix Fisher distribution [START_REF] Kent | A new unified approach for the simulation of a wide class of directional distributions[END_REF][START_REF] Lee | Bayesian attitude estimation with the matrix Fisher distribution on SO(3)[END_REF], is an analog (in the case of SO 3 (R)) of the Gaussian distribution in a flat space. The new orientation of Agent k at time T n can therefore be interpreted as a small random perturbation of the average local orientation given by A n k , where the perturbation size is measured by 1/ √ κ.

In Formula (2.6) and in the remainder of this paper, the manifold SO 3 (R) is endowed with its unique normalized Haar measure defined for any test function φ by:

(2.7)

SO 3 (R) φ(A) dA := 2 π π 0 S 2 φ(A(θ, n)) sin 2 (θ/2) dθ dn,
where dn is the uniform probability measure on the sphere S 2 . Here, a rotation matrix A ≡ A(θ, n) is parametrized by its rotation angle θ ∈ [0, π] and its axis n ∈ S 2 through Rodrigues' formula:

(2.8)

A(θ, n) := I 3 + sin θ [n] × + (1 -cos θ) [n] 2 × = exp(θ[n] × )
with n = (n 1 , n 2 , n 3 ) T and I 3 is the 3 × 3 identity matrix. For any vector w = (w 1 , w 2 , w 3 ) T ∈ R, [w] × is the antisymmetric matrix of the linear map R 3 ∋ u → w × u (where × denotes the cross product) which has the following expression:

(2.9) [w] × :=   0 -w 3 w 2 w 3 0 -w 1 -w 2 w 1 0   .
Additional details on the structure of SO 3 (R) can be found for instance in [START_REF] Huynh | Metrics for 3D rotations: Comparison and analysis[END_REF]. The IBM (2.3), (2.5) is schematically represented in Fig. 2.

This manuscript is for review purposes only.

A n k region time T n-1 k time T n k A k (T n k + 0) A k (T n k -0) Observation time T n-2 k 1
Figure 2: Schematic representation of the PDMP described in the text: the motion of Particle k is represented in physical space as the black broken dotted line. The body frame A k is represented with Ω k in red, u k in green and v k in blue. Each angular point of the trajectory corresponds to one of the jump times T n k . Between two jump times, the trajectory is the straight line spanned by Ω k and the body frame stays constant. The jump dynamics is depicted at time T n k . At this time, the observation region is colored in yellow and body frames of the other particles present in this region are depicted in light blue. The averaged body frame A n k is depicted with thick lightly colored arrows. The body frame before the jump A k (T n k -0) is drawn in broken lines whereas that after the jump A k (T n k + 0) is drawn in plain lines. A k (T n k + 0) is close, but not equal to A n k because of the noise intensity proportional to 1/κ. For clarity, the frames involved in the description of the jump are magnified.

2.1.2. Numerical simulations of the IBM. Unless otherwise specified, throughout this paper, a square box of side length L with periodic boundary conditions is used. As sensing kernel K, we use the indicator function of the ball centered at 0 and of radius R. Thus, an agent interacts with all its neighbors at a distance less than R (radius of interaction). Table 1 summarizes the model parameters.

For the numerical simulations presented in this paper, we have used the convenient framework offered by quaternions. Indeed, there is a group isomorphism between SO 3 (R) and H/{±1} where H is the group of unit quaternions. We can express the IBM (2.3), (2.5) using this representation (see [START_REF] Degond | Alignment of self-propelled rigid bodies: from particle systems to macroscopic equations[END_REF] and Section SM2). Roughly speaking, body-alignment as described here is equivalent to nematic alignment of the corresponding quaternions (nematic This manuscript is for review purposes only. alignment of a unit quaternions q to the mean direction Q is unchanged if q is replaced by -q, as opposed to polar alignment where the result depends on the sign of q). This is because a given rotation can be represented by two opposite quaternions and thus, the outcome of the alignment process should not depend of the choice of this representative. The numerical algorithm is described in Section SM3. Additionally, the quaternion framework also suggests to use order parameters derived from nematic alignment dynamics (such as in liquid crystal polymers). We shall use this analogy to define appropriate order parameters in Section 4.1.

All the simulations were written in Python using the SiSyPHE library [START_REF] Diez | SiSyPHE: A Python package for the Simulation of Systems of interacting mean-field Particles with High Efficiency[END_REF] specifically developed for the simulation of large-scale mean-field particle systems by the second author.

The implementation is based on the PyTorch [START_REF] Paszke | PyTorch: An imperative style, high-performance deep learning library[END_REF] library and more specifically on the GPU routines introduced by the KeOps [START_REF] Charlier | Kernel operations on the GPU, with autodiff, without memory overflows[END_REF] library. The computational details as well as the source code are freely available on the documentation website https://sisyphe.readthedocs.

io/. The outcomes of the simulations were analyzed and plotted using the NumPy [START_REF] Harris | Array programming with NumPy[END_REF] and Matplotlib [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF] libraries. The 3D particle plots were produced using VPython [START_REF] Scherer | VPython: 3D interactive scientific graphics for students[END_REF]. All the particle simulations have been run on a GPU cluster at Imperial College London using an Nvidia GTX 2080 Ti GPU chip.

A typical outcome of the IBM is shown in Figure 3 (see also Section SM1, Video 1) for a moderate number of particles (N = 3000). Throughout this paper, in the plots, we will represent each agent graphically by an elongated tetrahedron pointing in the direction of motion. The three large faces around the height will be painted in blue, green and magenta and the base will be in gold, as described in Fig. 3a. We notice that, starting from a uniformly random initial state (Fig. 3b), the system self-organizes in small clusters (Fig. 3c) and finally reaches a flocking equilibrium where all the agents have roughly the same body-orientation (Fig. 3d). We will see below that flocking is not necessarily the ultimate fate of the system, because it may be trapped in a so-called topologically protected state. To better understand these aspects, we first need to develop the continuum (or macroscopic) description of the system. This is done in the next section.

2.1.3. Relation with other collective dynamics models. We finally make a comparison with previous models. First, there is a version of the IBM where particles follow a stochastic differential equation (SDE) instead of a jump process [START_REF] Degond | A new flocking model through body attitude coordination[END_REF][START_REF] Degond | Quaternions in collective dynamics[END_REF]. Both the current and previous models have the same hydrodynamic model as macroscopic limit (see forthcoming section).

There are two reasons for us to prefer the jump process. First, its simulation is slightly easier

This manuscript is for review purposes only. and second, the coefficients of the macroscopic model are explicit, which is not so in the SDE case where they require the resolution of an auxiliary elliptic problem [START_REF] Degond | A new flocking model through body attitude coordination[END_REF][START_REF] Degond | Quaternions in collective dynamics[END_REF].

Beyond the present body-orientation model, numerous models of self-propelled particles have been proposed in the literature (see the review [START_REF] Vicsek | Collective motion[END_REF]). The most closely related one is the celebrated Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF]. There are several versions of this model: time-discrete ones [START_REF] Chaté | Collective motion of self-propelled particles interacting without cohesion[END_REF][START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF], time-continuous ones relying on an SDE description of the particle trajectories [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] and time-continuous ones using a jump process instead [START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF]. The latter version is the most closely related to the present work. In [START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF], the difference is that particles carry a single direction vector Ω k instead of a whole body frame. This vector gives the direction of self-propulsion. The particles follow a similar PDMP, namely

• The random jump times are defined in the same way: they follow an exponential law with constant parameter ν > 0. At jump times, the position is continuous and the direction vector Ω k is discontinuous with left and right states respectively denoted by Ω k (T n k -0) and Ω k (T n k + 0).

• Between two jump times T n k , T n+1 k , the direction vector Ω k does not change and the particle moves in straight line at speed c 0 > 0 in the direction given by Ω k (T n k + 0). 3 and, assuming that it is non-zero, the mean direction Ωn

• To pass from Ω k (T n k -0) to Ω k (T n k + 0), we compute the local flux given by J n- k = 1 N N j=1 K X k (T n k ) -X j (T n k ) Ω j (T n k -0) ∈ R
k = J n- k /|J n- k | ∈ S 2 at time T n k -0. Then, Ω k (T n k + 0) is drawn from a von Mises distribution on S 2 : Ω k (T n k + 0) ∼ MΩn k , with MΩ(Ω) = e κ( Ω•Ω) / S 2 e κ( Ω•Ω) dΩ, for Ω and Ω in S 2 .
So, the current model is an elaboration of [START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF] replacing self-propulsion directions by whole body frames and polar alignment of unit vectors (as expressed by the von Mises distribution on the sphere) by alignment of rotations matrices. Outcomes of numerical simulations of the This manuscript is for review purposes only.

Vicsek model do not show striking differences whether one uses any of the above mentioned versions (time-discrete, time-continuous with SDE or time-continuous with jump process).

Results given in [START_REF] Chaté | Collective motion of self-propelled particles interacting without cohesion[END_REF][START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] for the time-discrete version display the emergence of a global alignment together with the formation of clusters when the noise intensity 1/κ is not too big. The outcome strongly resembles what is shown in Fig. 3 for the body-orientation model, but for the depiction of the body orientation itself which is not provided by the Vicsek model. So, it is legitimate to wonder whether the inclusion of the full body orientation instead of the mere self-propulsion direction makes any change in the dynamics of the particle positions and direction vectors. In particular, do the particle positions and directions follow the same dynamics in the Vicsek and body orientation model? We will see below that this is not the case and that in certain circumstances, striking differences between the two models are obtained.

To show this, the use of the macroscopic limit of the IBM, as developed in the forthcoming section, will be of crucial importance.

2.2. The macroscopic body-alignment model.

Description of the model.

As soon as N is not very small, the IBM (2.3), (2.5) involves a large number of unknowns which makes its mathematical analysis virtually impossible. A reduced description, more amenable to mathematical analysis, is obtained through the macroscopic limit of the IBM, and consists of a system of partial differential equations.

This reduced description gives a valid approximation of the IBM in an appropriate range of parameters, namely

(2.10) N ≫ 1, R L ∼ c 0 ν L ≪ 1.
Throughout the remainder of this paper, we will focus on this regime. The macroscopic limit of the IBM (2.3), (2.5) has first been proposed in [START_REF] Degond | Alignment of self-propelled rigid bodies: from particle systems to macroscopic equations[END_REF] and leads to a model called "Self-Organized Hydrodynamics for Body orientation (SOHB)". The derivation relies on earlier work [START_REF] Degond | A new flocking model through body attitude coordination[END_REF][START_REF] Degond | Quaternions in collective dynamics[END_REF]. This derivation is "formally rigorous" in the sense that, if appropriate smoothness assumptions are made on the involved mathematical objects, the limit model can be identified rigorously as being the SOHB. For the reader's convenience, we summarize the main steps of this mathematical result in Section SM4.

The unknowns in the SOHB are the particle density ρ(t, x) and mean body-orientation A(t, x) ∈ SO 3 (R) at time t and position x = (x, y, z) ∈ R 3 . They satisfy the following set of equations:

∂ t ρ + c 1 ∇ x • (ρ Ae 1 ) = 0, (2.11a) ∂ t + c 2 (Ae 1 ) • ∇ x A + (Ae 1 ) × (c 3 ∇ x log ρ + c 4 r) + c 4 δ Ae 1 × A = 0. (2.11b)
The quantities r and δ have intrinsic expressions in terms of A [START_REF] Degond | A new flocking model through body attitude coordination[END_REF]. However, it is more convenient to write the rotation field A in terms of the basis vectors

Ω = Ae 1 , u = Ae 2 , v = Ae 3 .
With these notations, the vector r(t, x) ∈ R 3 and scalar δ(t, x) ∈ R fields are defined by

r := (∇ x • Ω) Ω + (∇ x • u) u + (∇ x • v) v, (2.12) δ := [(Ω • ∇ x ) u] • v + [(u • ∇ x )v] • Ω + [(v • ∇ x )Ω] • u. (2.13)
This manuscript is for review purposes only.

Here, for a vector field B(x) ∈ R 3 and a scalar field λ(x) ∈ R we denote by ∇ x • B, and ∇ x × B the divergence and curl of B respectively, by ∇ x λ, the gradient of λ and we set (B • ∇ x )λ = B • ∇ x λ with • the inner product of vectors in R 3 . We remind that × denotes the cross product and we refer to formula (2.9) for the definition of [w] × when w is a vector in R 3 . Alternate expressions of δ can be found in Section SM5 of the Supplementary Material.

The quantities c 1 , c 2 , c 3 , c 4 are functions of κ and c 0 given as follows:

c 1 c 0 = 2 3 1 2 + cos θ exp(κ( 1 2 +cos θ)) sin 2 ( θ 2 ) , (2.14)
c 2 c 0 = 1 5 ⟨2 + 3 cos θ⟩ exp(κ( 1 2 +cos θ)) sin 4 ( θ 2 ) cos 2 ( θ 2 ) , (2.15) c 3 c 0 = 1 κ , (2.16) c 4 c 0 = 1 5 ⟨1 -cos θ⟩ exp(κ( 1 2 +cos θ)) sin 4 ( θ 2 ) cos 2 ( θ 2 ) , (2.17)
where, for two functions f and g: [0, π] → R, we write

⟨f ⟩ g = π 0 f (θ) g(θ) dθ π 0 g(θ) dθ
∂ t ρ + c 1 ∇ x • (ρ Ω) = 0, (2.18a) D t A + [w] × A = 0, (2.18b)
where the convective derivative D t and the vector w are given by:

D t = ∂ t + c 2 Ω • ∇ x , (2.19) w = -Ω × F + c 4 δ Ω, with F = -c 3 ∇ x log ρ -c 4 r, (2.20)
This manuscript is for review purposes only.

Eq. (2.18a) is the mass conservation equation of the fluid. The vector Ω gives the direction of the fluid motion. The fluid velocity deduced from (2.18a) is c 1 Ω. Since c 1 /c 0 ∈ [0, 1] as can be seen from Fig. 4 (see also [START_REF] Degond | A new flocking model through body attitude coordination[END_REF] for a rigorous proof), the fluid motion is oriented positively along Ω and its magnitude is smaller than the particles self-propulsion velocity c 0 . This is because the average of vectors of identical norms has smaller norm. The quantity c 1 /c 0 can be seen as an order parameter [START_REF] Degond | Phase transitions and macroscopic limits in a BGK model of body-attitude coordination[END_REF] but we will not dwell on this issue here.

Eq. (2.18b) provides the rate of change of A with time along the integral curves of the vector field c 2 Ω as expressed by the convective derivative D t . Note that this vector field is not the fluid velocity c 1 Ω since c 2 ̸ = c 1 . It can be interpreted as the propagation velocity of A when w is zero. Since D t A is the derivative of an element of SO 3 (R), it must lie in the tangent space to SO 3 (R) at A which consists of all matrices of the form W A with W antisymmetric. This structure is indeed satisfied by Eq. (2.18b) since, from the definition (2.9), the matrix [w] × is antisymmetric. It can be shown that the SOHB system is hyperbolic [START_REF] Degond | Some properties of Self-Organized Hydrodynamics for body orientation[END_REF].

In fact, Eq. (2.18b) shows that the vector w is the instantaneous rotation vector of the frame A(t, X(t)), where t → X(t) is any solution of dX dt = c 2 Ω(t, X(t)). Indeed, Eq. (2.18b) can be equivalently written as a system of equations for (Ω, u, v) of the form Fig. 5 gives a graphical representation of the actions of these two infinitesimal rotations.

D t Z = w × Z, with Z = Ω, u, v

Relation with other models.

To better understand how the SOHB model (2.11) relates to other models, we re-write the equation for Ω as follows:

(2.21)

D t Ω = P Ω ⊥ F,
where P Ω ⊥ is the 3 × 3 projection matrix on the orthogonal plane to the vector Ω and is written P Ω ⊥ = I 3 -Ω ⊗ Ω with ⊗ standing for the tensor (or outer) product. Eq. (2.21) bears similarities and differences with the momentum equation of isothermal compressible fluids.

The latter is exactly recovered if the following three modifications are made:

1. the projection matrix P Ω ⊥ is removed from (2.21) (i.e. it is replaced by I 3 );

2. c 2 = c 1 in the convective derivative D t (see (2.19)); 3. c 4 = 0 in the expression of F (see (2.20)).
Indeed, under these three modifications, we get the following system for (ρ, U) where U = c 1 Ω is the fluid velocity:

∂ t ρ + ∇ x • (ρU) = 0, (∂ t + U • ∇ x )U = -Θ ∇ x log ρ.
This is the isothermal compressible Euler equations with the fluid temperature Θ = c 1 c 3 .

We now investigate what consequences follow from undoing the above three modifications, one by one.

This manuscript is for review purposes only. (Ω, u, v) denotes the position of the frame at time t while (Ω ′ , u ′ , v ′ ) is its position at time t + dt with dt ≪ 1. The frame at time t is denoted in plain colors (red for Ω, green for u and blue for v) while that at time t + dt is in light colors. The motion of the vectors is indicated by a segment of circle in black color. (a) Action of Ω × F: the vectors F and Ω × F are in plain and light black respectively. The vector F is shown with unit norm for the ease of the representation but could be of any norm in reality. The passage from (Ω, with Ω (and assuming that all functions are smooth) leads to D t |Ω| 2 = 0, which guarantees that |Ω| is constant along the integral curves of the vector field c 2 Ω. Thus, if |Ω| = 1 at time t = 0, it will stay so at any time.

F u Ω v ′ v Ω × F u ′ Ω ′ 1 (a) Action of Ω × F u ′ -c 4 δΩ Ω v ′ Ω ′ u v 1 (b) Action of -c4δΩ
u, v) to (Ω ′ , u ′ , v ′ ) is via
2. Having c 2 ̸ = c 1 is a signature of a loss of Galilean invariance. This is consistent with the fact that the microscopic system is not Galilean invariant as well, Indeed, there is a distinguished reference frame where the particle speed is c 0 . Of course, this speed does not remain equal to c 0 in frames that translate at constant speed with respect to this frame.

So far, with the introduction of P Ω ⊥ and different constants c 2 ̸ = c 1 but still with c 4 = 0, the system for (ρ, Ω) is decoupled from the equations for u and v and is written (see Eqs. (2.18a), (2.21) with F given by (2.20) in which c 4 = 0):

∂ t ρ + c 1 ∇ x • (ρ Ω) = 0, (2.22a) D t Ω = -c 3 P Ω ⊥ ∇ x log ρ. (2.22b)
This is nothing but the hydrodynamic limit of the Vicsek particle model (known as "Self-Organized Hydrodynamics (SOH)") as established in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF]. This system has This manuscript is for review purposes only.

been shown to be hyperbolic [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] and to have local-in-time smooth solutions [START_REF] Degond | Hydrodynamic models of self-organized dynamics: derivation and existence theory[END_REF].

3. When c 4 ̸ = 0, in addition to the pressure gradient, a second component of the force F appears. This component depends on the full rotation matrix A through Ω, u, v and their gradients (see Eq. 2.12). It is thus truly specific of the body orientation model.

We are now going to compare the IBM and the SOHB models on a set of explicit stationary solutions of the SOHB model described in the next section.

3. Special solutions of the macroscopic model.

3.1. Three classes of explicit solutions. In this section, we exhibit three different classes of global-in-time solutions of the SOHB model (2.18). They are special classes of a larger family of solutions which will also be introduced. All these solutions are characterized by uniform (i.e. independent of the spatial coordinate) fields ρ, r and δ. From now on we fix a wave-number (inverse of the length) ξ ∈ R \ {0} and define

(3.1) ω = ξ c 4 , λ = c 2 + c 4 .
We denote by x = (x, y, z) T the coordinates of x in the basis (e 1 , e 2 , e 3 ).

3.1.1. Flocking state. The flocking state (FS) is a trivial but important special solution of the SOHB model (2.18) where both the density and rotation fields are constant (i.e. independent of time) and uniform:

ρ(t, x) ≡ ρ 0 = constant, A(t, x) ≡ A 0 = constant, ∀(t, x) ∈ [0, ∞) × R 3 .
3.1.2. Milling orbits. We have the following Lemma 3.1. The pair (ρ, A) consisting of a constant and uniform density ρ(t, x) = ρ 0 = constant and the following rotation field:

A(t, x) = Ãmill (t, z) =   cos(ωt)
sin(ωt) cos(ξz)sin(ωt) sin(ξz) sin(ωt) cos(ωt) cos(ξz)cos(ωt) sin(ξz) 0 sin(ξz) cos(ξz)

  (3.2) = A(-ωt, e 3 ) A(ξz, e 1 ), (3.3)
is a solution of the SOHB system (2.18), where ω and ξ are given by (3.1). We recall that A(θ, n) is the rotation of axis n ∈ S 2 and angle θ ∈ R defined by (2.8). This solution will be referred to as a milling orbit (MO).

The proof of this lemma is deferred to Section SM6. The MO is independent of x and y.

Its initial condition is (3.4) A mill (0, z) = A(ξz, e 1 ) =   1 0 0 0 cos(ξz) -sin(ξz) 0 sin(ξz) cos(ξz)   .
The initial direction of motion (the first column of A mill (0, z)) is independent of z and aligned along the x-direction, i.e. Ω(0, z) ≡ e 1 . As z varies, the body-orientation rotates uniformly This manuscript is for review purposes only.

about the x-direction with spatial angular frequency ξ. As the rotation vector is perpendicular to the direction of variation, (3.4) is called a "perpendicular twist". As time evolves, the rotation field is obtained by multiplying on the left the initial perpendicular twist by the rotation A(-ωt, e 3 ). This means that the whole body frame undergoes a uniform rotation about the z-axis with angular velocity -ω. As a consequence, the direction of motion is again independent of z. It belongs to the plane orthogonal to z and undergoes a uniform rotation about the z-axis. Consequently, the fluid streamlines, which are the integral curves of c 1 Ω, are circles contained in planes orthogonal to z of radius

c 1 ω = c 1 c 4 1
ξ traversed in the negative direction if ξ > 0. These closed circular streamlines motivate the "milling" terminology. It can be checked that the MO satisfies:

r = ξ (sin(ωt), cos(ωt), 0) T , δ = 0.
As announced, r and δ are uniform but r depends on time. Actually, Ω × r = ξe 3 is independent of time. The MO is depicted in Fig. 6 and its dynamics is visualized in Video 2 (see Many examples of milling (also known as vortex) solutions have been observed in the collective dynamics literature as well as in biological systems [START_REF] Calovi | Swarming, schooling, milling: phase diagram of a data-driven fish school model[END_REF][START_REF] Couzin | Self-organized lane formation and optimized traffic flow in army ants[END_REF][START_REF] Vicsek | Collective motion[END_REF]. On the modelling side, milling states have not been observed so far in alignment models without the inclusion of an additional process such as an attraction-repulsion force between the agents [START_REF] Carrillo | Double milling in self-propelled swarms from kinetic theory[END_REF], a bounded cone of vision [START_REF] Costanzo | Spontaneous emergence of milling (vortex state) in a Vicsek-like model[END_REF] or an anticipation mechanism [START_REF] Gerlee | Impact of anticipation in dynamical systems[END_REF]. The body-orientation framework is, to the best of our knowledge, a new situation in which milling can be observed just with

Section SM1). 1 (a) t = 0 1 (b) t > 0
This manuscript is for review purposes only.

alignment assumptions. Milling states can also be found in physical systems. A typical and important example is the motion of a charged particle in a uniform magnetic field, resulting in the formation of so-called cyclotron orbits. Once again, in the body-orientation framework, an external field is not needed and self-induced cyclotron orbits emerge only from the variations of the internal body-orientation. Here, the analog of the magnetic field would be Ω × r and the cyclotron frequency would be ω. Note that ω is under the control of coefficient c 4 which depends on the noise intensity 1/κ.

A(t, x) = Ãhtw (t, x) =   1 0 0 0 cos (ξ(x -λt)) -sin (ξ(x -λt)) 0 sin (ξ(x -λt)) cos (ξ(x -λt))   (3.5) = A(ξ(x -λt), e 1 ), (3.6)
is a solution of the SOHB system (2.18) where ξ and λ are defined by (3.1). This solution will be referred to as a helical traveling wave (HW).

The proof of this lemma is given in Section SM6.2. The HW is independent of y and z.

Its initial condition is (3.7) A htw (0, x) = A(ξx, e 1 ) =   1 0 0 0 cos(ξx) -sin(ξx) 0 sin(ξx) cos(ξx)   .
Here the self-propulsion direction is still independent of x and equal to e 1 . Also, the body orientation still rotates uniformly about e 1 with spatial angular frequency ξ but when x is varied instead of z. This means that the body orientation is now twisted when varied along the propagation direction. So, this initial condition is called a "parallel twist". In the HW, the self propulsion direction Ω remains constant in time and uniform in space. The initial twist is propagated in time in this direction at speed λ and gives rise to a traveling wave Ãhtw (t, x) = Ãhtw (0, xλt).

Note that the traveling wave speed λ depends on the noise intensity 1/κ and is different from the fluid speed c 1 . So, the frame carried by a given fluid element followed in its motion is not fixed but rotates in time. Since Ω does not change, the fluid streamlines are now straight lines parallel to e 1 . So, as a fluid element moves, the ends of the frame vectors u and v follow a helical trajectory with axis e 1 , hence the terminology "helical traveling waves" for these solutions. It can be checked that r = 0, δ = ξ, and again, r and δ are spatially uniform as announced. The HW is depicted graphically in Fig. 7. Its dynamics is visualized in Video 3 (see Section SM1). The HW belongs to a larger class of solutions described in Section SM6.2.

This manuscript is for review purposes only. The pair (ρ, A) consisting of a constant and uniform density ρ(t, x) = ρ 0 = constant and the following rotation field:

(3.8) A(t, x) = A ξ,θ (t, z) := A(-ωt, e 3 ) A θ - π 2 , e 2 A(ξ(z -λt), e 1 ),
is a solution of the SOHB system (2.18). We recall that A(θ, n) is the rotation of axis n ∈ S 2 and angle θ ∈ R. This solution will be referred to as a Generalized topological Solution (GS).

The proof of this lemma is deferred to the Supplementary Material SM6.3. Each of the three previous classes of solutions can be obtained for specific values of the parameters ξ and θ.

• When ξ = 0, the solution A 0,θ is constant for any θ, which corresponds to a FS.

• When θ = π 2 and ξ ∈ R, then λ = 0 and the rotation with respect to the y-axis is equal to the identity: the solution A ξ,π/2 is therefore equal to the MO (3.3).

• When θ = 0 and ξ ∈ R then λ = c 2 and the solution A ξ,0 is equal to

A ξ,0 =   0 -sin(ξ(z -λt)) -cos(ξ(z -λt)) 0 cos(ξ(z -λt)) -sin(ξ(z -λt)) 1 0 0   , λ = c 2 + c 4 ,
which is an HW along the z-axis. The situation is analogous when θ = π.

This manuscript is for review purposes only.

All these solutions have a non-zero gradient in the body-orientation variable which is always along the z-axis. This gradient is controlled by the parameter ξ. However, in the GS, the direction of motion Ω (or fluid velocity) is not necessarily parallel nor perpendicular to this gradient. Specifically, Ω has a constant polar angle equal to the parameter θ. The behavior of the solution is then a combination of the two previously introduced phenomena:

milling around the z-axis and a travelling wave of the body-orientation variable along the same axis. The applet accessible at https://www.glowscript.org/#/user/AntoineDiez/folder/ MyPrograms/program/BOfield provides a graphical representation of the GS for arbitrary polar angles using VPython [START_REF] Scherer | VPython: 3D interactive scientific graphics for students[END_REF] and with the same conventions as in Fig. 6.

In the following, we will focus on each of these two elementary behaviors, i.e. the standard milling and helical travelling wave solutions, and in particular on their topological properties.

The study of the full continuum of generalized solutions is left for future work. However, we will encounter GS obtained from a perturbed milling solution in Section 5.4.

3.2. Some properties of these special solutions. Clearly, in the definitions of the MO and HW, the choice of reference frame is unimportant. So, in the whole space R 3 , such solutions exist in association with any reference frame. In a square domain of side-length L with periodic boundary conditions, periodicity imposes some constraints on the direction of the reference frame. For simplicity, we will only consider the case where the reference frame has parallel axes to the sides of the square and ξ is linked to L by an integrality condition

L ξ = 2π n, with n ∈ Z \ {0}.
The study of the stability of the MO and the HW is left for future work. By contrast, the FS is linearly stable as the SOHB system is hyperbolic [START_REF] Degond | Some properties of Self-Organized Hydrodynamics for body orientation[END_REF]. However, there is no guarantee that the FS at the level of the IBM is stable. Indeed, there are strong indications that the FS is not stable for the Vicsek model [START_REF] Chaté | Collective motion of self-propelled particles interacting without cohesion[END_REF] for some parameter ranges and a similar trend is likely to occur here.

We can now answer the question posed at the end of Section 2. This means that information (which triggers adjustments in Ω) propagates downstream the fluid by contrast to the Vicsek case where it propagates upstream. While the reason for this difference is unclear at this stage, we expect that it may induce large qualitative differences in the behavior of the system in some cases. This point will be investigated in future work.

Numerical simulation of the SOHB will be subject to future work. Here, we will restrict ourselves to the MO and HW for which we have analytical formulas. In the next section, using these two special solutions, we verify that the SOHB model and the IBM are close in an appropriate parameter range.

3.3. Agreement between the models. In this section we use the MO and HW to demonstrate the quantitative agreement between the SOHB model (2.11) and the IBM (2.3), (2.5) in the scaling (2.10). In the simulations below, we consider a periodic cube of side-length L and choose

(3.9) R = 0.025, ν = 40, c 0 = 1, L = 1, ξ = 2 π, so that R L = c 0 ν L = 0.025 ≪ 1,
ensuring that the scaling (2.10) is satisfied. Furthermore, we see that the choice of ξ is such that the twists in the MO or HW have exactly one period over the domain size. We sample N particles according to the initial condition (3.4) of the MO and simulate the IBM (2.3), (2.5). We recall that the average direction Ω(t) of the exact MO (3.2) is spatially uniform at any time and undergoes a uniform rotation motion about the z-axis. So, we will compare Ω(t) with the average direction Ω(t) of all the particles of the IBM, where

Ω(t) = (Ω 1 , Ω 2 , Ω 3 ) T is defined by: Ω = N k=1 Ω k (t) | N k=1 Ω k (t)| ,
(provided the denominator is not zero, and where we recall that Ω k (t) = A k (t) e 1 ). To ease the comparison, we compute the azimuthal and polar angles of Ω respectively defined by:

(3.10) φ := arg(Ω 1 + iΩ 2 ) ∈ [0, 2π), θ = arccos(Ω 3 ) ∈ [0, π],
where arg(x + iy) stands for the argument of the complex number x + iy. We note that the corresponding angles φ and θ of Ω(t) are given by

(3.11) φ(t) = -ω t = -2π c 4 (κ) t, θ = π/2,
This manuscript is for review purposes only.

where we have used (3.1) and (3.9) to compute the value of ω.

Fig. 8a shows the azimuthal angle φ as a function of time over 5 units of time, for increasing particle numbers: N = 5 10 4 (green curve), N = 1.5 10 5 (orange curve) and N = 1.5 10 6 (blue curve). Note that for very small values of N , the macroscopic model loses its relevance: below a few thousand particles we only observe a noisy behavior, not shown in the figure.

For the considered range of particle numbers, we notice that the angle φ decreases linearly with time, which shows that the behavior of the IBM is consistent with the exact solution (3.11). However, quantitatively, we see that |d φ/dt| depends on the particle number and decreases with increasing particle number. We investigate this behavior in more detail in This manuscript is for review purposes only.

the two classes of special solutions, the MO and the HW. To compare the models, we compute the following macroscopic quantities:

• For the MO: starting from a sampling of the initial condition (3.4), we measure the angular velocity |d φ/dt| in a similar way as in the previous section. Given the parameter choice (3.9), the theoretical value of |dφ/dt| predicted by (3.2) is |ω| = 2πc 4 (κ)

where the function c 4 is given by (2.17).

• For the HW, starting from a sampling of the initial condition (3.7), we measure the wave speed. To this aim, using (2.2), we compute the mean body-orientation A of the agents in a slice of size 10 IBM and finite particle number effects as explained in Section 3.3.1 may eventually destabilize the IBM. Then, in most cases, its solution is seen to transition towards an approximation of the FS after some time. This transition implies a change of the topology of the solution which, from initially non-trivial, becomes trivial, since the winding number of the FS is zero. One may wonder whether the evolution towards a FS is slower if the initial state has non-trivial topology and exhibits some kind of "topological protection" against noise-induced perturbations. To test this hypothesis quantitatively, we first need to develop appropriate indicators. This is done in the next section.

4. Order parameters and topological indicators. We will use two types of indicators.

The first one is the global order parameter which will discriminate between the various types of organization of the system (disorder, MO or HW and FS). The second type of indicators are based on analyzing the roll angle. They will enable a finer characterization of topological phase transitions.

4.1. Global order parameter. We first introduce the following scalar binary order parameter which measures the degree of alignment between two agents with body-orientations A,

à ∈ SO 3 (R) : (4.1) ψ(A, Ã) := 1 2 A • à + 1 4 .
In the quaternion framework (see Section 2.1.2 and SM2 for details), we have

(4.2) ψ(A, Ã) = (q • q) 2 ,
where q and q are two unit quaternions respectively associated to A and Ã, and q • q indicates the inner product of two quaternions. This expression makes it clear that ψ(A, Ã)

∈ [0, 1]. The
This manuscript is for review purposes only.

square exponent in (4.2) indicates that ψ(A, Ã) measures the nematic alignment of the two associated unit quaternions, as it should because two opposite quaternions represent the same rotation. We note that ψ(A, Ã) = 1 if and only if à = A. On the other hand, ψ(A, Ã) = 0 if and only if A • à = -1/2, which corresponds to the two rotation axes being orthogonal and one rotation being an inversion about its axis.

The Global Order Parameter (GOP) of a system of N agents at time t > 0 is the average of all binary order parameters over all pairs of particles:

(4.3) GOP N (t) = 1 N (N -1) k̸ =ℓ ψ A k (t), A ℓ (t) .
From (4.3) we have GOP N (t) ∈ [0, 1]. A small GOP N indicates large disorder and a large one, strong alignment. This is a global measure of alignment, by contrast to a local one where ψ would be averaged over its neighbors only (and the result, averaged over all the particles).

This global measure of alignment allows us to separate the MO and HW from the FS as shown below, which would not be possible with a local one.

The GOP (4.3) can also be defined at the continuum level. As shown in Section SM4, in the macroscopic limit, the particles become independent and identically distributed over R 3 ×SO 3 (R), with common distribution ρ M A where (ρ, A) satisfies the SOHB system (2.11) and M A is the von Mises distribution (2.6). Therefore, the GOP of a solution of the SOHB system (ρ, A) is obtained as (4.3) where the sum is replaced by an integral, A k (t) is replaced by A distributed according to the measure (ρ M A )(t, x, A) dx dA and A ℓ (t) is replaced by à distributed according to the same measure, but independently to A. Therefore,

GOP(ρ, A) := (R 3 ×SO 3 (R)) 2 ψ(A, Ã) ρ(x) ρ(x) M A(x) (A) M A(x) ( Ã) dx dx dA d Ã.
Using (2.7) and (2.8) one can prove that for any A ∈SO 3 (R), we have (4.4)

SO 3 (R) A M A (A) dA = c 1 (κ) c 0 A,
with c 1 (κ) defined by (2.14) and c 0 being the particle speed. Using (4.1), we obtain:

(4.5) GOP(ρ, A) = 1 2 c 1 (κ) c 0 2 R 3 ×R 3 A(x) • A(x) ρ(x) ρ(x) dx dx + 1 4 .
From now on, we let ρ be the uniform distribution on a square box of side-length L. We can compute the GOP corresponding to each of the three solutions defined in Section 3.1.

For the MO (3.2), HW (3.5) and GS (3.8), for all time t > 0, in all cases, the GOP remains equal to:

(4.6) GOP 1 = 1 4 c 1 (κ) c 0 2 + 1 4 .
For the FS, A(x) ≡ A = constant and the GOP is equal to

(4.7) GOP 2 = 3 4 c 1 (κ) c 0 2 + 1 4 .
This manuscript is for review purposes only.

Note that the GOP:

GOP 0 = 1 4 ,
corresponds to a disordered state of the IBM where the body-orientations of the particles are chosen independently and randomly uniformly (or equivalently to the SOHB case κ → 0 in (4.6) and (4.7)). For the typical value κ = 10 used in our simulations, one can compute that:

(4.8) GOP 1 ≃ 0.45, GOP 2 ≃ 0.85.
The GOP values between GOP 1 and GOP 2 can be reached by generalized HW as shown in Section SM6.4.

4.2. Roll angle.

Definition. Let

A = [Ω, u, v] ∈ SO 3 (R) be a body-orientation. Let θ ∈ [0, π], φ ∈ [0, 2π
) be the spherical coordinates of Ω defined by (3.10) (omitting the bars). We let {Ω, e θ , e φ } be the local orthonormal frame associated with the spherical coordinates (θ, φ) and we define p(Ω) = e φ and q(Ω) = -e θ . Then we define the rotation matrix

R(Ω) := [Ω, p(Ω), q(Ω)] =   sin θ cos φ -sin φ -cos θ cos φ sin θ sin φ cos φ -cos θ sin φ cos θ 0 sin θ   .
Since u and v belong to the plane spanned by p(Ω) and q(Ω), we let ζ ∈ [0, 2π) be the angle between p(Ω) and u. Then, it is an easy matter to show that A = R(Ω) A(ζ, e 1 ). In aircraft navigation, θ, φ and ζ are respectively called the pitch, yaw and roll angles: the pitch and yaw control the aircraft direction with respect to the vertical and in the horizontal plane respectively, while the roll controls the plane attitude (see Fig. 10a). These angles are related to the Euler angles. The construction of the roll angle ζ is summarized in Figure 10b.

Pursuing the analogy with aircraft navigation, we see from Fig. 5 that F controls variations of pitch and yaw while δ controls variations of roll.

As an example, we examine the pitch, yaw and roll of the three solutions of the SOHB model (2.11) described in Section 3.1.

FS:

A is constant and uniform. Then, the pitch, yaw and roll are also constant and uniform.

MO:

A is given by (3.2) (see Figs.

HW:

A is given by (3.5) (see Fig. 7). Then, we have R(Ω) = I 3 and ζ = ξ (xλ t).

The pitch and yaw are constant and uniform while the roll is uniform on planes of constant x. It depends on x and time through the traveling phase xλ t. Here, the non-trivial topology results from the roll making a complete turn when x increases by the quantity 2π/ξ.

The goal of the next section is to see how we can recover the roll field from the simulation of a large particle system.

This manuscript is for review purposes only. compute an average roll on planes of constant z. To cope with the discreteness of the particles, we will rather consider slices comprised between two planes of constant z. If the distance ∆z between the planes is chosen appropriately, we can access to both the average and the variance of the roll. They will be collected into one single vector, the Roll Polarization in planes of constant z or RPZ. A similar quantity characterizes the HW, the Roll Polarization in planes of constant x or RPX. Below, we detail the construction of the RPZ. Obviously the procedure is the same (changing z into x) for the RPX.

We assume that the domain is a rectangular box of the form D := [0, L x ] × [0, L y ] × [0, L z ], and L z = n (2π/ξ) with n ∈ Z \ {0}. The domain D is partitioned into M slices of fixed size across z, where M is a fixed integer. For m ∈ {1, . . . , M }, the slice S m is defined by:

S m := [0, L x ] × [0, L y ] × m -1 M L z , m M L z .
Let us consider a system of N agents with positions and body-orientations (X k , A This manuscript is for review purposes only.

1. Disordered state: the particles are drawn independently uniformly randomly in the product space D× SO 3 (R). For each m, the RPZ (4.9) is an average of uniformly distributed vectors on the circle and its norm is therefore close to 0. The RPZ-curve is thus reduced to the origin, as shown in Figure 12a; 2. FS: the positions of the particles are drawn independently uniformly in D and their body-orientations independently according to a von Mises distribution M A 0 with a fixed mean body orientation A 0 ∈ SO 3 (R). In this case, for all slices, the corresponding RPZ (4.9) is an average of identically distributed vectors on the circle whose distribution is peaked around the same point of the unit circle, and the peak is narrower as κ is larger. Therefore, the RPZ vectors (4.9) concentrate on a point near the unit circle (Figure 12b). The RPZ-curve reduces to a single point different from the origin;

3. MO: the positions of the particles are drawn independently uniformly in D. Then for a particle at position x, its body-orientation is drawn independently according to a von Mises distribution M A mill (0,z) with A mill (0, z) defined by (3.4) (with ξ = 2π/L z ). This time, the von Mises distribution is peaked around a point which depends on z. For each slice, the position of the RPZ (4.9) depends on m. Since A mill (0, z) is L z -periodic, the RPZ-curve is a discrete closed circle (Figure 12c). Note that the RPX-curve of a HW is similar.

From Figure 12, we realize that three quantities of interest can be extracted from the RPZ-curve:

1. the distance of its center of mass to the origin d z :

(4.10)

d z = 1 M M m=1 ūm ,
2. its mean distance to the origin rz :

w z = 1 2π M m=1 δβ m+1/2 ,
(see e.g. [62, p. 176]).

The subscript z indicates that the slicing has been made across z. Similar quantities with an index 'x' will correspond to the slicing made across x. We have a similar conclusion with (d x , rx , w x ) for a disordered state or an FS. For an HW, we have (d x , rx , w x ) ≈ (0, 1, w) with w ̸ = 0. Thus, monitoring either or both triples (according to the situation) will give us an indication of the state of the system in the course of time.

In particular, non-trivial topological states are associated with non-zero winding numbers w x or w z . In practice, we will use the nonzero-rule algorithm to compute the winding numbers

This manuscript is for review purposes only.

numerically [62, p. 176].

5. Topological phase transitions: are the MO and HW topologically protected?. As pointed out in Section 3.4, for the IBM, the MO and HW are only metastable: they typically persist for a finite time before degenerating into a FS. This is in stark contrast with the macroscopic model for which they persist for ever. The transition of a MO or HW to a FS implies a topological change. To analyze whether the MO or HW are more robust due to their non-trivial topological structure (i.e. are topologically protected), we will compare them with similar but topologically trivial initial conditions (Sections 5.1, 5.2 and 5.3). We also test their robustness against perturbed initial conditions and show that, in this case, MO may transition to GS (Section 5.4). In the Supplementary Material SM8, we investigate rarer events, where an MO does not transition directly to an FS but through a HW.

5.1. Initial conditions. In Section 5.2, we will compare the solutions of the IBM with different initial conditions using the perpendicular or parallel twists as building blocks. Some will have a non-trivial topology and the others, a trivial one. Specifically we define the following initial conditions.

5.1.1. Milling orbit. Let D = [0, L]×[0, L]×[0, 2L
] be a rectangular domain with periodic boundary conditions and let ξ = 2π/L. We consider the following two initial conditions:

• Double mill initial condition MO1:

(5.1)

A m,1 (0, z) = A(ξ z, e 1 ), z ∈ [0, 2L],
where we recall again that A(θ, n) is the rotation of axis n ∈ S 2 and angle θ ∈ R defined by (2.8). This initial condition has non-trivial topology: the curve generated by the end of the vector u in the (y, z)-plane as z ranges in [0, 2L] makes two complete turns around the origin in the same direction. Thus, this initial condition has winding number equal to 2.

• Opposite mills initial condition MO2:

(5.2)

A m,2 (0, z) = A(ξ z, e 1 ), z ∈ [0, L], A(-ξ z, e 1 ), z ∈ [L, 2L].
This initial condition has trivial topology: starting from z = 0, the curve generated by the end of the vector u makes one complete turn around the origin in the counterclockwise direction until it reaches z = L but then reverses its direction and makes a complete turn in the clockwise direction until it reaches z = 2L. Thus, this initial condition has winding number equal to 0 and has trivial topology.

• Perturbed double mill initial condition MO3:

(5.3) A m,3 (0, z) = A(ξ z + √ σB z , e 1 ), z ∈ [0, 2L],
where (B z ) z is a given one-dimensional standard Brownian motion in the z variable and σ > 0 is a variance parameter which sets the size of the perturbation. The Brownian motion is subject to B 0 = B 2L = 0 (i.e. it is a Brownian bridge). Similarly to the initial condition MO1 (5.1), this initial condition has a nontrivial topology, in this case a winding number equal to 2.

This manuscript is for review purposes only.

Helical traveling wave. Let now

D = [0, 2L] × [0, L] × [0, L].
Compared to the previous case, the domain has size 2L in the x-direction instead of the z-direction. Let again ξ = 2π/L. We consider now the following two initial conditions:

• Double helix initial condition HW1:

(5.4)

A h,1 (0, x) = A(ξ x, e 1 ), x ∈ [0, 2L],
This initial condition has non-trivial topology and has winding number equal to 2 by the same consideration as for initial condition MO1.

• Opposite helices initial condition HW2:

(5.5)

A h,2 (0, x) = A(ξ x, e 1 ), x ∈ [0, L], A(-ξ x, e 1 ), x ∈ [L, 2L].
Again, by the same considerations as for MO2, this initial condition has trivial topology, i.e. winding number equal to 0.

Observation of topological phase transitions.

We initialize the IBM by drawing N positions independently uniformly randomly in the spatial domain and N body-orientations independently from the von Mises distribution M A(0,x) where A(0, x) is one of the initial conditions MO1 or MO2. Then, we run the IBM and record the various indicators introduced in Section 4 as functions of time. The results are plotted in Fig. 13, as plain blue lines for the solution issued from MO1 (the topologically non-trivial initial condition), and as broken orange lines for that issued from MO2 (the topologically trivial one). We proceed similarly for the two initial conditions HW1 and HW2 and display the results in Fig. 14. See also Videos 4 to 7 in Section SM1 supplementing Fig. 13 and Videos 8 to 11 supplementing Fig. 14.

Figs. 13a and 14a display the GOP. We observe that, for all initial conditions, the GOP has initial value GOP 1 , which is consistent with the fact that the initial conditions are either MO or HW. Then, again, for all initial conditions, at large times, the GOP has final value GOP 2 which indicates that the final state is a FS. This is confirmed by the inspection of the second line of figures in Figs. [START_REF] Bostan | Reduced fluid models for self-propelled particles interacting through alignment[END_REF] and HW2 initial conditions, no topological phase transition is needed to reach the FS. The differences in the initial topology of the solutions induce strong differences in the trajectories followed by the system.

For the topologically non-trivial initial conditions MO1 or HW1, the system remains in the MO or HW state for some time; hence it follows the macroscopic solution during this phase. Indeed, the GOP displays an initial plateau at the value GOP 1 , while the triplet of topological indicators stays at the value (0, 1, 2), which characterize the MO or HW state. For

This manuscript is for review purposes only.

MO1, this is also confirmed by the yaw φ (Fig. 13c, blue curve), which varies linearly in time and by the pitch θ (Fig. 13b blue curve) which is constant in time, consistently with the MO solution of the macroscopic model (Section 3.1.2) (see also Fig. 8a for the linear variation of the yaw). The duration of this initial phase, also referred to as the persistence time, is significantly longer for HW1 than for MO1. In our experiments, the former can reach several hundred units of time and sometimes be infinite (up to our computational capacity). By contrast, the latter is random and of the order of ten units of time. After this initial plateau, the GOP decreases until it reaches a minimum at a time highlighted in Figs. 13, 14 and subsequent figures by a gray shaded zone, showing that the system passes through a state of maximal disorder. Around that time, r has a sharp drop which is another confirmation of an increased disorder. The topological transition precisely occurs at this time with a transition of the winding number from 2 to 0 through a short sequence of oscillations. However, r has not reached 0 and d has already started to increase, which suggests that disorder is not complete.

At this time also, the linear variation of φ suddenly stops and φ remains constant afterward, while θ shows a small oscillation and jump. For HW1, θ and φ are initially plateauing with small oscillations. At the time when the system leaves the HW state (around t ≃ 178), we observe a sudden drop of φ from 2π to π which indicates that the system suddenly reverses its average direction of motion. The GOP starts to decrease significantly before this time so we can infer that during the time period between t ≃ 125 and t ≃ 178, even though the mean direction of motion Ω remains constant, groups of particles of almost similar proportions are moving in opposite directions, which preserves the average direction of motion (and may explain the oscillations during the initial persistence phase). This is confirmed by Video 8

(see description in Section SM1). Then, once this minimum is reached, the GOP increases quickly to finally reach the value GOP 2 of the FS. Likewise, r and d quickly reach the value 1 while the winding number stays at the value 0.

By contrast to the previous case, the system immediately leaves the topologically trivial initial conditions MO2 or HW2 as shown by the GOP immediately leaving the value GOP 1 .

For HW2 the GOP increases right after initialization and smoothly reaches the value GOP 2 , at a much earlier time than HW1. The trend is different for MO2. In this case, the GOP first decreases. Then, after a minimum value, it increases again and smoothly reaches the value GOP 2 at a time similar to MO1. The initial decay of the GOP for the MO2 solution can be explained by the fact that the macroscopic direction Ω turns in opposite directions for the two opposite mills, thus decreasing the global order. For HW2, the macroscopic direction stays constant and uniform. So, it is the same for the two opposite helices, giving rise to a larger GOP. The mean radii rz and rx stay constant it time, showing that the evolutions of MO2 and HW2 do not involve phases of larger disorder. The quantity d x increases monotonically towards the value 1 while d z is subject to some oscillations close to convergence. This is due to the fact that the RPZ or RPX curves stay arcs of circles with decreasing arc length for the RPX and with some arc length oscillations for the RPZ as displayed in Videos 7 and 11.

Of course, the winding number stays constant equal to 0 as it should for topologically trivial solutions. In both the MO2 and HW2 cases, θ and φ remain constant throughout the entire simulation. In the MO2 case, this is the consequence of the two counter-rotating mills which preserve the direction of motion on average. In the HW2 case, this is due to the fact that there is no variation of the direction of motion for HW solutions in general (see also Video 6 and This manuscript is for review purposes only.

Video 10). Again, we observe that the convergence towards the FS takes more time for HW2 than for MO2. This points towards a greater stability of the HW-type solutions compared to the MO ones. GOP which reaches the value GOP 2 of the FS. A similar reproducibility of the results has been observed for the other initial conditions (MO2, HW1, HW2) (not shown). this final steady behavior is linked to a winding number equal to 1 after a transition around t ≃ 12. Consequently, a milling behavior is observed in Fig. 16g for the mean azimuthal angle. This angle evolves linearly but with a slower speed, approximately divided by 2, after the transition, as expected since the winding number has dropped from 2 to 1. However, the final mean polar angle θ shown in Fig. 16f is not equal to π/2. Since the gradient in bodyorientation is along the z-axis, this indicates that the final state corresponds to a GS rather than a standard MO. This demonstrates that the family of generalized topological solutions enjoys some greater stability. The transition between MO and GS has not been observed when starting from a non-perturbed initial state. However, starting with perturbed initial conditions, the MO and GS with winding number 1 seem stable during several tens of units of time.

(a) (b) (c) (d) (e) (f)
The transition between MO and GS with different winding numbers happens when the perturbation size is large enough and seems to be the typical behavior: out of 6 independent simulations for values of σ evenly spread between 258 and 876, 5 simulations led to a MO or a GS with winding number 1 stable during more than 50 units of time. The other one

This manuscript is for review purposes only.

led to a FS. We can think that the perturbation brings the system to a state closer to the MO with winding number 1, in particular due to the stochastic spatial inhomogeneities of the perturbation. On the particle simulations, we observe that the density of agents does not remain uniform, which creates different milling zones with possibly different milling speeds depending on the local gradient of body-orientations. The denser region then seems to attract the other particles before expanding into the full domain. The global direction of motion is not necessarily preserved during this process. In comparison, starting from an unperturbed MO with winding number 2, the density remains uniform and the system is globally subject to numerical errors which homogeneously degrade the topology up to the point that the system becomes closer to a FS. The situation is analogous when the size of the perturbation is too large as shown in Figs. 16i, 16k, 16l for σ = 1000 : the MO is preserved during less than 5 units of time and after an immediate drop of the GOP, the system quickly reaches a FS.

5.5. Critique. The existence of a persistence time for the MO1 and HW1 solutions suggests that they enjoy some kind of topological protection against the noisy perturbations induced by the IBM and that MO2 and HW2 do not have such protection. However, since explicit solutions of the SOHB model for the initial conditions MO2 and HW2 are not available, it is not possible to assess the role of noise in the observed evolutions of the MO2 and HW2 solutions. So, further investigations are needed to confirm that non-trivial topology actually provides increased robustness against perturbations. Moreover, the MO1 is robust against perturbed initial conditions. The MO and GS with winding number 1 seem to be much more more stable than with winding number 2. This manuscript is for review purposes only. 4. Beyond the mean-field limit N → ∞, it would be interesting to quantify the fluctuation about the mean-field, for instance through a large deviation approach (see e.g. [6,7,11,[START_REF] Dawson | Large deviations from the McKean-Vlasov limit for weakly interacting diffusions[END_REF][START_REF] Fernandez | A Hilbertian approach for fluctuations on the McKean-Vlasov model[END_REF][START_REF] Lancellotti | On the fluctuations about the Vlasov limit for N-particle systems with mean-field interactions[END_REF]).

5. Direct numerical simulations of the macroscopic model need to be developed to answer some of the questions raised by the study of topological protection (see Section 5).

This manuscript is for review purposes only.

6. It is desirable to develop more sophisticated topological indicators to gain better insight into the topological structure of the solutions.

7. The multiscale approach developed here could be extended to other geometrically structured systems involving e.g. a wider class of manifolds which would enlarge the applicability of the models.

(see Fig. 3a for details on the representation of the body orientation). We notice the counterrotation of the particle directions about the z axis in the bottom and top halves of the domain, corresponding to the opposite mills. These two counter-rotations gradually dissolve while the solution approaches the FS. about the z axis until an instability initiated at time t ≈ 10 gradually disrupts this organization.

However, the disruption does not drive the system to an FS, but rather to a HW as shown by the alternations of blue, green and magenta colors propagating along the particle orientations. Initially, the RPX-curve is reduced to the origin, showing total disorder across the x direction.

Then, after some chaotic transient, a closed curve enclosing the origin is formed. This curve initially stays close to the origin, still showing strong disorder. But gradually, the radius of the curve increases and approaches the unit circle. Thus, across x, the topology is initially undefined, but when it builds up, it shows its non-trivial character, the emerging RPX-curve having non-zero winding number about the origin. SM2. Quaternion framework. Despite its formal simplicity, the SO 3 (R)-framework used in the definition of the Individual Based Model is not well suited to numerical simulations due to the high computational cost required to store and manipulate rotation matrices. A more efficient representation of rotations in R 3 is the quaternion representation based on the group isomorphism Φ :

H/ ± 1 -→ SO 3 (R) q -→ Φ(q) : w ∈ R 3 → {q[w]q * } ∈ R 3 ,
where the 3-dimensional vector w = (w 1 , w 2 , w 3 ) T ∈ R 3 is identified with the pure imaginary quaternion denoted by [w] = iw 1 + jw 2 + kw 3 and q * denotes the conjugate quaternion

This manuscript is for review purposes only.

to q. Conversely, the pure imaginary quaternion q = iq 1 + jq 2 + kq 3 is identified with the 3-dimensional vector denoted by {q} := (q 1 , q 2 , q 3 ) T . Note that for any quaternion q and any vector w ∈ R 3 , the quaternion q[w]q * is a pure imaginary quaternion. The group of unit quaternions is denoted by H and is homeomorphic to the sphere S 3 ⊂ R 4 .

We refer the reader to [SM5, Section 2] and [SM4, Appendix A] where details about the equivalence between the two representations can be found. Note that [SM4] studies a model in a full quaternion framework. Table SM1 below summarizes how the different objects can be computed in either of the two representations.

Matrix Quaternion

Orientation A ∈ SO 3 (R) q ∈ H/ ± 1 such that Φ(q) = A Flux J k = j K(X k -X j )A j Q k = j K(X k -X j ) (q j ⊗ q j -1/4I 4 ) Mean orientation A = arg max{A → A • J} q ∈ H eigenvector associated to the largest eigenvalue of Q Von Mises distribu- tion M A (A) = exp(κA • A) Z M q (q) = exp(2κ(q • q) 2 ) Z
Table SM1: Matrix vs quaternion formulation SM3. Numerical methods. The IBM (2.3), (2.5) has been discretized within the quaternion framework using the time-discrete algorithm described in Table SM2 below. This table

shows one iteration of the algorithm during which the positions X n k ∈ R 3 and orientations q n k ∈ H for k ∈ {1, . . . , N } are updated into X n+1 k and q n+1 k respectively.

Algorithm: Iteration n → n + 1 of the time-discrete algorithm 1. Update the positions:

for k ∈ {1 . . . , N }, set X n+1 k = X n k + c 0 {q n k [e 1 ](q n k ) * } ∆t 2
Q n k = 1 N N j=1 K(X n k -X n j ) (q n j ⊗ q n j -1 4 I 4 ).
4. Update the orientations: for k ∈ I compute one unit eigenvector q n k of Q n k of maximal eigenvalue and draw q n+1 k ∼ M q n k .

Table SM2: One iteration of the time-discrete algorithm

This manuscript is for review purposes only.

At step 2, the Poisson process is discretized with a time step ∆t during which the indices of the jumping agents are recorded. In the simulations ∆t has to be chosen small enough so that the event that an agent jumps twice or more during a time interval of size ∆t is negligible.

In all the simulations, we take ∆t such that ν ∆t = 10 -2 .

At step 3, a random quaternion q sampled from a von Mises distribution with prescribed mean orientation q can be obtained as q = qr where r ∈ H is sampled from a von Mises distribution with mean orientation 1 (see [SM5, Proposition 9]). An efficient rejection algorithm to sample von Mises distributions can be found in [SM9].

All the simulations in this paper take place in a periodic box of size L = (L x , L y , L z ). The observation kernel K is the indicator of the ball centered at 0 and of radius R > 0. The six parameters of the simulations are summarized in Table 1.

Finally, we would like to stress that the quaternion formulation is not only a convenient numerical trick. The equivalence it provides between body-orientation models and models of nematic alignment of polymers in dimension four has been exploited in [SM2] to study phase transitions in the body alignment model. the agents, seen as random variables in R 3 × SO 3 (R) become independent and identically distributed. Their law is given by the kinetic particle distribution f which satisfies the following PDE:

∂ t f + c 0 Ae 1 • ∇ x f = ν (ρ f M A K * f -f ),
where ρ f ≡ ρ f (t, x) is the local spatial density:

ρ f (t, x) = SO 3 (R) f (t, x, A) dA,
and A K * f ≡ A K * f (t, x) is the local average body-attitude defined by

A K * f (t, x) := arg max A∈SO 3 (R) A • J K * f (t, x),
computed from the local flux:

J K * f ≡ J K * f (t, x) := R 3 ×SO 3 (R) K(x -y) A f (t, y, A) dy dA.
From a mathematical point of view, the probability distribution f ≡ f (t, x, A) is obtained as the limit in law of the empirical measure of the N -particle system. We refer to [SM7] where a rigorous proof of this result is presented for a similar model, and to [SM1] for a related work on the Vicsek model.

In the macroscopic regime the agent interactions become strong, which is expressed by the following hydrodynamic scaling:

ε ∼ c 0 ν L ∼ R L ≪ 1,
This manuscript is for review purposes only.

where L is a typical macroscopic length-scale of the system (such as the typical size of the flock). We define c0 = ενL = O(1) and c ′ 0 = c 0 /c 0 . Then, defining dimensionless time and space variables t ′ and x ′ such that x = Lx ′ and t = (L/c 0 )t ′ , we obtain (dropping the primes for simplicity):

(SM4.1)

∂ t f ε + c 0 Ae 1 • ∇ x f ε = 1 ε (ρ f ε M A f ε -f ε ) + O(ε), where A f ε ≡ A f ε (t, x) := arg max A∈SO 3 (R) A • J f ε (t, x), and
J f ε ≡ J f ε (t, x) := SO 3 (R) A f ε (t, x, A) dA.
This last expression is obtained by Taylor expanding

J K * f ε = J f ε + O(ε 2
) and means that the interactions between the agents become spatially localized in the macroscopic regime.

The macroscopic model is obtained by formally taking the limit ε → 0 in (SM4.1). If such a limit exists, it is necessarily of the form

(SM4.2) f ε -→ ε→0 ρ M A
where ρ ≡ ρ(t, x) and A ≡ A(t, x) depend on t and x. Thus, the limiting distribution is fully described by the spatial density of agents and their average orientation. To obtain a system of equations for (ρ, A), we first use the local conservation of mass: integrating (SM4.1) over SO 3 (R) and noting the right-hand side vanishes, it holds that,

∂ t SO 3 (R) f ε dA + c 0 SO 3 (R) A e 1 • ∇ x f ε dA = O(ε).
When ε → 0, assuming (SM4.2) and using (4.4), we obtain (2.11a).

To obtain an equation for A, it could be tempting to pursue this approach and multiply (SM4.1) by A before integrating it over SO 3 (R). However, the term resulting from the righthand side of (SM4.1) does not vanish but equals (using (4.4) again):

1 ε SO 3 (R) A (ρ f ε M A f ε -f ε ) dA = 1 ε c 1 c 0 ρ f ε A f ε -J f ε ̸ = 0.
Due to the factor ε -1 , its limit as ε → 0 is unknown. An easy fix can be found if, instead of multiplying Eq. (SM4.1) by A before integrating it over SO 3 (R), we multiply it by the quantity

ψ A f ε (A) := A T f ε A -A T A f ε .
The rationale for using this quantity is because we aim to find an equation for the time-derivative of A. Such a derivative must lie in the tangent space to SO 3 (R) at A, denoted by T A . This suggests to multiply (SM4.1) by an element of T A . Given an arbitrary matrix A, a natural way to obtain an element of T A is to take its orthogonal projection on T A , which is given by 1 2 (A -AA T A). We could therefore choose to multiply (SM4.1) by this quantity. But a further simplification is possible by noting that this quantity is equal to 1 2 A ψ A (A) and that 1 2 A does not depend on A and so can be factored

This manuscript is for review purposes only.

SM6.1. Proof of Lemma 3.1. Starting from the initial condition (3.4), we are looking for solutions of (2.11b) of the form

A(t, x) =   cos(ωt) u 1 (t, z) v 1 (t, z) -sin(ωt) u 2 (t, z) v 2 (t, z) 0 u 3 (t, z) v 3 (t, z)   ,
where ω ∈ R is an angular velocity which will be related to the parameters of the problem later and where the basis vectors u = (u 1 , u 2 , u 3 ) T and v = (v 1 , v 2 , v 3 ) T depend only on the z variable and time. In this situation, Equation (2.11a) is trivially satisfied which means that the system stays homogeneous in space. Solutions of this form have to satisfy three geometrical constraints which ensure that A ∈ SO 3 (R). The first two ones are Ω × u = v and v × Ω = u, which lead to (SM6.1)

A(t, x) =   cos(ωt) sin(ωt)v 3 (t, z) -sin(ωt)u 3 (t, z) -sin(ωt) cos(ωt)v 3 (t, z) -cos(ωt)u 3 (t, z) 0 u 3 (t, z) v 3 (t, z)   .
The third one is a normalization constraint:

(SM6.2) ∀t > 0, ∀z ∈ R, u 3 (t, z) 2 + v 3 (t, z) 2 = 1.
Using (SM6.2), we define a function α ≡ α(t, z) such that

u 3 (t, z) = sin(α(t, z)), v 3 (t, z) = cos(α(t, z)).
A direct computation shows that for A of the form (SM6.1), we have

r = (∂ z u 3 ) u + (∂ z v 3 ) v, δ = 0.
Therefore, Eq. (2.11b) can be rewritten more concisely into:

(SM6.3) ∂ t A + c 4 [Ω × r] × A = 0,
where we recall Eq. (2.9) for the definition of [ ] × . A direct computation shows that

(SM6.4) Ω × r = (v 3 ∂ z u 3 -u 3 ∂ z v 3 ) e 3 = (∂ z α) e 3 .
Inserting this in (SM6.3) implies that u 3 (t, z) ≡ u 3 (z) and v 3 (t, z) ≡ v 3 (z) are independent of time. We then observe that:

(SM6.5) A(t, x) = A(-ωt, e 3 ) A(α(z), e 1 ), where we recall Eq. (2.8) for the meaning of A. Therefore, using (SM6.3) and (SM6.4), we obtain:

-

ω [e 3 ] × A + c 4 (∂ z α) [e 3 ] × A = 0,
from which we deduce that A satisfies (2.11b) if and only if α and ω satisfy:

c 4 ∂ z α = ω, which implies (SM6.6) α(z) = ω c 4 z + ᾱ,
where ᾱ is a constant, which can be interpreted as the phase at the origin z = 0. To recover Eq. (3.2), we just need to take ᾱ = 0 and define ξ = ω/c 4 . Eq. (3.3) follows from (SM6.5).

This manuscript is for review purposes only.

SM6.2. Generalized HW and proof of Lemma 3.2. Starting from the initial condition (3.7), we are looking for solutions of (2.11b) of the form

A(t, x) =   1 0 0 0 cos(α(t, x)) -sin(α(t, x)) 0 sin(α(t, x)) cos(α(t, x))   ,
for a real-valued function α of the t and x variables only. In this case, Ω is a constant vector and Equation (2.18a) is trivially satisfied. Moreover a direct computation shows that:

r = 0, δ = (∂ x α)(t, x).
As a consequence, Eq. (2.21) is trivially satisfied and straightforward computations show that Eq. (2.11b) reduces to

∂ t α + (c 2 + c 4 ) ∂ x α = 0.
This last equation is a linear transport equation with velocity c 2 + c 4 , the solutions of which are given by (SM6.7) α(t, x) = α 0 (x -(c 2 + c 4 )t)

for any initial condition α 0 ∈ L 1 loc (R). In the case of (3.7), α 0 (x) = ξ x. However, we see that there are as many different solutions as functions in L 1 loc (R). Such general solutions are called "generalized HW". SM6.4. GOP of the MO and generalized HW. The GOP (given by Eq. (4.5)) of the MO and HW do not depend on time and only depend on the function α defined respectively by (SM6.6) and (SM6.7). Using Eq. (4.5), we can compute that the GOP is equal to:

GOP = 1 2 c 1 (κ) c 0 2 1 + 2 |⟨u⟩| 2 + 1 4 ,
This manuscript is for review purposes only.

where ⟨u⟩ denotes the spatial average of the vector u with respect to ρ (here the with respect to the uniform measure on the domain since ρ is constant and uniform). With the previous notations, we obtain |⟨u⟩| 2 = ⟨cos α⟩ 2 + ⟨sin α⟩ 2 , For the generalized HW, depending on the choice of α, the GOP can take any value between where A(0, x) is given by (3.4) with ξ = 2π/L (winding number equal to 1).

SM8.1. From milling orbit to helical wave. Here, we report on the occurrence of transitions from a MO to a HW. Among twenty independent simulations, this transition occurred only once (the other cases being a transition from a MO to a FS). We run the IBM and record the time-evolution of a set of indicators as shown in Fig. SM2 (see also supplementing videos 12 to 14 in Section SM1).

As shown in Fig. SM2a, the GOP does not converge towards GOP 2 characterizing the FS, but towards an intermediate value between GOP 1 (which characterizes MO or HW) and GOP 2 . As explained in Section SM6.4, such values of the GOP can be attained by a generalized helical wave solution (as can be observed in Video 12). The pitch θ (Fig. SM2b) and yaw φ (Fig. SM2c) behave like in the milling-to-flocking transition (see Figs. 13b and13c) except for small-amplitude, slow-frequency oscillations appearing after the topological transition time.

This may be due to some competition between two attractors, the FS and the HW, which being

This manuscript is for review purposes only. also display large oscillations. The quantity rz drops, and at the same time, d z remains small, while the winding number w z has strong oscillations, indicating a state of large disorder across z, which is consistent with the fact that the temporary HW order is organized in a different direction. However, we see that w z has a calmer period between two series of oscillations.

This calmer period corresponds to the interval of time during which the temporary HW order prevails. Eventually the triplet converges to the value (1, 1, 0) characterizing the FS.

A k e 1 A k e 1 1 (

 11

Figure 3 :

 3 Figure 3: (a) Graphical representation of particles and their body orientations as elongated tetrahedra pointing towards the self-propulsion direction with blue, magenta and green large faces and gold bases. (b,c,d) Snapshots of a typical output of the simulation at three different times (b) Time=0, (c) Time=4 and (d) Time=40. Parameters: N = 3000, L = 1, R = 0.075, κ = 20, ν = 5, c 0 = 0.2. see also Section SM1, Video 1.

.

Fig. 4

 4 Fig.4provides a graphical representation of these functions.

Figure 4 :

 4 Figure 4: Dimensionless coefficients c i /c 0 as functions of the inverse of concentration parameter 1/κ. Blue curve c 1 /c 0 , orange curve c 2 /c 0 , green curve c 3 /2c 0 and red curve c 4 /c 0 . At the crossover value κ * ≃ 2.58, the sign of c 2c 1 changes (see Section 3.2).

2. 2 . 2 .

 22 Interpretation of the model. To better understand what the SOHB system (2.11) does, we re-write it as follows:

.

  This describes a rigid body rotation of the frame {Ω, u, v} with angular velocity w. The rotation vector w has two components. The first one is Ω × F and tends to relax Ω towards F. Due to its expression (2.20), the force F includes two contributions: that of the pressure gradient -c 3 ∇ x log ρ and that of gradients of the body orientation through the vector -c 4 r. The second component of the rotation vector is -c 4 δΩ and corresponds to a rotation of the body frame about the self propulsion direction Ω driven by gradients of the body orientation through the scalar -c 4 δ. The contributions of gradients of body orientation in the two components of the rotation vector are under the control of the single coefficient c 4 .

Figure 5 :

 5 Figure 5: Graphical representations of the two components of the infinitesimal rotation.(Ω, u, v) denotes the position of the frame at time t while (Ω ′ , u ′ , v ′ ) is its position at time t + dt with dt ≪ 1. The frame at time t is denoted in plain colors (red for Ω, green for u and blue for v) while that at time t + dt is in light colors. The motion of the vectors is indicated by a segment of circle in black color. (a) Action of Ω × F: the vectors F and Ω × F are in plain and light black respectively. The vector F is shown with unit norm for the ease of the representation but could be of any norm in reality. The passage from (Ω, u, v) to (Ω ′ , u ′ , v ′ ) is via an infinitesimal rotation of axis Ω × F. (b) Action of δ: the vector -c 4 δΩ is shown in black. The vectors Ω and Ω ′ are identical and collinear to -c 4 δΩ. The passage from (Ω, u, v) to (Ω ′ , u ′ , v ′ ) is via an infinitesimal rotation of axis Ω.

1 .

 1 Figure 5: Graphical representations of the two components of the infinitesimal rotation.(Ω, u, v) denotes the position of the frame at time t while (Ω ′ , u ′ , v ′ ) is its position at time t + dt with dt ≪ 1. The frame at time t is denoted in plain colors (red for Ω, green for u and blue for v) while that at time t + dt is in light colors. The motion of the vectors is indicated by a segment of circle in black color. (a) Action of Ω × F: the vectors F and Ω × F are in plain and light black respectively. The vector F is shown with unit norm for the ease of the representation but could be of any norm in reality. The passage from (Ω, u, v) to (Ω ′ , u ′ , v ′ ) is via an infinitesimal rotation of axis Ω × F. (b) Action of δ: the vector -c 4 δΩ is shown in black. The vectors Ω and Ω ′ are identical and collinear to -c 4 δΩ. The passage from (Ω, u, v) to (Ω ′ , u ′ , v ′ ) is via an infinitesimal rotation of axis Ω.

Figure 6 :

 6 Figure 6: Graphical representation of the milling orbit (MO) at (a): initial time, and (b): time t > 0. The frame vectors Ω, u and v are represented at a certain number of points of the (O, x, y) and (O, y, z) planes. In (b), the rotation motion of the frame vectors is depicted by dotted circles of the color of the corresponding frame vector. The red dotted circle can be seen as a depiction of the fluid streamlines. See also Section SM1, Video 2.

3. 1 . 3 .

 13 Helical traveling wave. We have the following Lemma 3.2. The pair (ρ, A) consisting of a constant and uniform density ρ(t, x) = ρ 0 = constant and the following rotation field:

1 (Figure 7 :

 17 Figure 7: Graphical representation of the helical traveling wave (HW) at (a): initial time, and (b): time t > 0. See Fig. 6 for captions. See also Section SM1, Video 3.

  1.3 namely whether the inclusion of the full body orientation makes any change in the dynamics of the particle positions and directions compared to the Vicsek model. To this end, we consider the corresponding macroscopic models, i.e. the SOH model (2.22) for the Vicsek model and the SOHB model (2.11) for the body-orientation dynamics. If we initialize the SOH model with uniform initial density ρ and mean direction Ω, inspection of (2.22) shows that the solution remains constant in time and thus corresponds to a flocking state of the Vicsek model. In the SOHB model, the three classes of solutions described in the previous sections (the FS, MO and HW) also have uniform initial density ρ and mean direction Ω. If the dynamics of the particle positions and directions in the body orientation model was the same as in the Vicsek model, these three classes of solutions should have a constant mean direction Ω. However, it is not the case for the MO, where Ω changes with time and is subject to a planar rotation. This means that gradients of body attitude do have a non-trivial influence on the direction of motion of the particles and that the body orientation model does not reduce to a Vicsek model for the particle positions and directions.There is another, more subtle, difference between the two models concerning the dynamics of Ω. It does not concern the MO and HW but we discuss it here in relation with the previous paragraph. Indeed, Fig.4reveals that the velocities c 1 and c 2 for the SOHB model crossover This manuscript is for review purposes only. at a certain value κ * of the concentration parameter. The coefficients c 1 and c 2 for the SOH model can be found in[START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF], Fig.A1(b) and appear to satisfy c 1 > c 2 for the whole range of values of κ, i.e. do not exhibit any crossover. In particular, at large noise, the propagation velocity c 2 of Ω in the SOHB model is larger than the mass transport velocity c 1 .

3. 3 . 1 .

 31 The IBM converges to the macroscopic model as N → ∞. In this section, we numerically demonstrate that the solutions of the IBM converge to those of the macroscopic model in the limit N → ∞ and investigate the behavior of the IBM at moderately high values of N .

Fig. 8b where

  Fig.8bwhere the difference between the measured angular velocity |d φ/dt| and the theoretical prediction 2πc 4 (κ) is plotted as a function of N . Each data point (blue dot) is an average of 10 independent simulations. This figure confirms that, as N increases, |d φ/dt| decreases and converges towards 2πc 4 (κ). The inset in Fig.8bshows the same data points in a loglog-scale with the associated regression line (orange solid line). We observe that the error between the measured and theoretical angular velocities behaves like N -α with a measured exponent α ≃ 1.01 which is close to the theoretical value α = 1 derived in Section SM7 of the Supplementary Material.

Figure 8 : 3 . 3 . 2 .

 8332 Figure 8: (a) Time evolution of the angle φ for three values of N : N = 0.05 10 6 (green curve), N = 0.15 10 6 (orange curve) and N = 1.5 10 6 (blue curve). (b) Difference between the measured angular velocity |d φ/dt| and the theoretical value 2πc 4 (κ). Each data point (blue dot) is an average of 10 independent simulations with the error bar showing one standard deviation. Solid black horizontal line at 0 for convenience. Inset: same data in log-log scale and regression line (solid orange line). Parameters: L = 1, ξ = 2π, R = 0.025, ν = 40, c 0 = 1, κ = 10.

  -3 along the x-axis (which is the global direction of motion) as a function of time. As predicted by(3.5) the coefficient A 22 of the mean orientation is a periodic signal. The inverse of the period of this signal (obtained through a discrete Fourier transform) gives the traveling wave speed of the HW. The theoretical value predicted by (3.5) is given by λ = c 2 (κ) + c 4 (κ) where the function c 2 is given by(2.15).The output of these simulations is shown in Figs. 9a for the MO and 9b for the HW.They respectively display the angular velocity and traveling wave speed obtained by running the IBM for a discrete set of values of κ (big blue dots). By comparison, the black dotted curves show the theoretical values as functions of κ. For the parameters of Fig.9, the order of magnitude of the standard deviation of 10 independent simulations is 10 -3 . The relative error between the average measured value and its theoretical prediction varies between 2% and 5% on the whole range of concentration parameters considered.These figures show an excellent agreement between the prediction of the macroscopic SOHB model and the results obtained by running the IBM when the number of particles is large. This confirms that the SOHB model provides an excellent approximation of the IBM, at least during a certain period of time which is a function of the particle number. We will see below that fluctuations induced by the finite number of particles may eventually destabilize the MO and lead to a HW or a FS. As these solutions are associated with different topological structure, these transitions will be analyzed as topological phase transitions in the forthcoming sections.3.4. Topology. Both the MO and HW have non-trivial topology: inspecting the perpendicular twist (3.4) (see also Fig.6a), we observe that the two-dimensional curve generated by the end of the vector u in the (y, z)-plane as one moves along the z-axis is a closed circle. A similar observation can be made on the parallel twist (3.7) (see Fig.7a) as one moves along the x-axis. Both curves have therefore non-zero winding numbers about the origin. When the domain is R 3 , these winding numbers are ±∞ (where the sign corresponds to that of ξ) as these curves make an infinite number of turns. If the domain has finite extension L along the z-axis (in the MO case) or the x-axis (in the HW case) and, due to the periodic boundary conditions, L is related to ξ by L = n 2π/ξ with n ∈ Z \ {0}, then the winding numbers are equal to n. As observed on Formulas (3.2) and (3.5) (or on Figs6b and 7b), this initial non-trivial topological structure is propagated in time.When we initialize particles by sampling the initial conditions (3.4) or (3.7), we expect that the solution of the IBM remains an approximation of the MO (3.2) or HW (3.5) respectively as evidenced in Section 3.3.2. However, noise induced by both the inherent stochasticity of the This manuscript is for review purposes only.

Figure 9 :

 9 Figure 9: (a) MO: angular velocity |dφ/dt| as a function of 1/κ. (b) HW: traveling wave speed λ as a function of 1/κ. Measured values from the IBM at discrete values of κ (big blue dots) and theoretical prediction from the SOHB model (dotted black curve). Parameters: N = 1.5 10 6 , L = 1, ξ = 2π, R = 0.025, ν = 40, c 0 = 1.

  6). Using Eq. (3.3), we have R(Ω) = A(-ω t, e 3 ) and the roll is given by ζ = ξz. The pitch and yaw are constant and uniform. The roll is constant in time and is also uniform on planes of constant z. The non-trivial topology of the MO results from the roll making a complete turn when z increases by the quantity 2π/ξ.

Figure 10 : 4 . 2 . 2 .

 10422 Figure 10: (a) Pitch, yaw and roll angles of an aircraft with body orientation [Ω, u, v] (original picture released under the Creative Commons CC0 license by https://pixabay.com). (b) Construction of the roll angle of A = [Ω, u, v], where the vectors Ω, u and v are respectively in red, green and blue. The local frame is (Ω, p(Ω), q(Ω)) where p(Ω) and q(Ω)) and the plane generated by them are in purple. u and v belong to this plane. ζ is the angle between p(Ω) and u.

  k ), indexed by k ∈ {1, . . . , N }. Each body orientation A k has roll ζ k ∈ [0, 2π). We define the discrete This manuscript is for review purposes only. RPZ for Slice m, ūm , by (4.9) ūm := 1 N m k∈Im (cos ζ k , sin ζ k ) T ∈ R 2 , where I m = {k ∈ {1, . . . , N }, X k ∈ S m } and N m is the cardinal of I m . Note that the RPZ ūm has norm smaller than one. The unit vector ūm /|ū m | or equivalently, its angle with the vector (1, 0) T gives the average roll in S m . The euclidean norm |ū m | is a measure of the variance of the set of roll angles {ζ k } k∈Im . If this variance is small, then |ū m | ∼ 1, while if the variance is large, |ū m | ≪ 1. When plotted in the plane R 2 , the set of RPZ {ū m } m=1,...,M forms a discrete curve referred to as the RPZ-curve. It will be used to characterize the topological state of the particle system. A summary of this procedure is shown in Figure 11.

Figure 11 : 4 . 2 . 3 .

 11423 Figure 11: Construction of the RPZ and graphical representation. The spatial domain D is partitioned into M slices represented in different colors (top left). In each slice S m , we have I m particles with roll ζ k each of them plotted in the particle's local plane spanned by p(Ω k ),q(Ω k ) (top right: we plot 3 particles in the slice S 1 ). Note that the local planes of different particles of the same slice may not coincide when imbedded in R 3 . For this given slice, the RPZ ūm is computed and plotted in R 2 (bottom right). The RPZ has norm smaller than 1 and belongs to the unit disk, whose boundary, the unit circle, is plotted for clarity. The RPZ of each slice is then plotted on a single figure in the same color as the slice it corresponds to (bottom left). This collection of points forms a discrete curve (here a fragment of a circle): the RPZ-curve.

3 .

 3 its winding number about the origin w z : for m ∈ {1, . . . , M }, letβ m = arg (ū m ) 1 + i(ū m ) 2 ∈ [0, 2π) (with ūm = ((ū m ) 1 , (ū m ) 2 ) T ) and δβ m+1/2 ∈ [-π, π) be such that δβ m+1/2 ≡ β m+1β m modulo 2π, where we let β M +1 = β 1 . Then:

  Fig. 12d provides a graphical illustration of the triple (d z , rz , w z ). For the examples given above, this triple has the following values: Disordered state: (d z , rz , w z ) = (0, 0, ND), where ND stands for "undefined", (4.12) FS: (d z , rz , w z ) ≈ (1, 1, 0), (4.13) MO: (d z , rz , w z ) ≈ (0, 1, w), with w ̸ = 0. (4.14)This manuscript is for review purposes only.

Figure 12 :

 12 Figure 12: Examples of RPZ-curves: in each figure, the roll Polarization RPZ vectors corresponding to M = 1000 slices are plotted. The color bar to the right of each figure assigns a unique color to each slice. The same color is used to plot the corresponding RPZ. In each figure the unit circle and its center are represented in blue. (a) Disordered state: all RPZ concentrate near the origin. (b) FS: all RPZ concentrate on a point close to the unit circle. (c) MO (3.4): the RPZ-curve is a discrete circle centered at the origin and of radius close to unity. The total number of particles is N = 1.5 • 10 6 . Note that in Figs. (a) and (b), all RPZ are superimposed and only the last one (in magenta color) is visible. (d) Quantifiers of RPZ curve morphology: point G (in red) is the center-of-mass of the RPZ curve and d z is its distance to the origin O (shown in blue). The mean radius rz of the RPZ curve is illustrated by the circle in black broken line which has same radius. The winding number, which is the number of turns one makes following the spectrum of colors in the same order as in the color bar from bottom to top (the green arrow indicates the direction of progression along the RPZ curve) is w z = -1 in this example.

Figure 13 :

 13 Figure 13: Examples of solutions of the IBM for initial conditions sampled from the double mill MO1 (plain blue curves) and the opposite mills MO2 (boken orange curves). The following indicators are plotted as functions of time: (a) Global Order Parameter (GOP) (see Eq. (4.3)). Horizontal lines at GOP values 0.25, 0.45 and 0.85 materialize the special values GOP 0 , GOP 1 and GOP 2 respectively corresponding to totally disordered states, MO or HW, and FS (see Eqs. (4.6)-(4.8)). (b) Pitch angle θ of the global particle average direction Ω (see (3.10)). (c) Yaw φ of Ω. (d) Distance of center of mass of RPZ curve to the origin d z (see (4.10)). (e) Mean distance of RPZ curve to the origin rz (see (4.11)). (f) Winding number of RPZ curve w z (see (4.10)). Gray shaded zones highlight a small region around the time of minimal GOP for the MO1 solution. Parameters: N = 3 10 6 , R = 0.025, κ = 10, ν = 40, c 0 = 1, L = 1, ξ = 2π. See also Videos 4 to 7 in Section SM1.

Figure 14 :

 14 Figure 14: Examples of solutions of the IBM for initial conditions sampled from the double helix HW1 (plain blue curves) and the opposite helices HW2 (broken orange curves). The following indicators are plotted as functions of time: (a) Global Order Parameter (GOP). (b) Pitch angle θ of Ω. (c) Yaw φ of Ω. (d) Distance of center of mass of RPX curve to the origin d x . (e) Mean distance of RPX curve to the origin rx . (f) Winding number of RPX curve w x . Gray shaded zones highlight a small region around the time of minimal GOP for the HW1 solution. The HW2 and HW1 solutions are computed during 200 and 250 units of time respectively. The two simulations have reached equilibrium by their final time. Parameters: N = 3 10 6 , R = 0.025, κ = 10, ν = 40, c 0 = 1, L = 1, ξ = 2π. See caption of Fig. 13 for further indications. See also Videos 8 to 11 in Section SM1.

5. 4 .

 4 Robustness against perturbations of the initial conditions. In this section, we study the robustness of the MO when the initial condition is randomly perturbed as described by the initial condition MO3 (5.3). Three typical outcomes for three different values of the perturbation size σ are shown in Fig. 16. For each value of σ, the temporal evolution of the four main indicators are shown: the GOP (Figs. 16a, 16e, 16i), the mean polar angle or pitch (Figs. 16b, 16f, 16j), the mean azimuthal angle or yaw (Figs. 16c, 16g, 16k) and the winding number along the z-axis (Figs. 16d, 16h, 16l). For small to moderate values (approximately σ < 100), the outcomes of the simulation are the same as in Fig. 13 and areThis manuscript is for review purposes only.

Figure 15 :

 15 Figure 15: GOP as a function of time for 20 independent simulations of the transition from a MO to a FS starting from MO1. The parameters are the same as the ones on Figure 13.

6 . 2 . 3 .

 623 Discussion and conclusion. An Individual Based Model describing the alignment of body-orientations in 3D and its macroscopic limit have been presented. The model involves new kinds of internal degrees of freedom involving geometrical constraints, here due to the manifold structure of SO 3 (R), leading to new types of self-organized phenomena. In particular, the macroscopic model has been shown to host special solutions with non-trivial topological structures. Corresponding solutions of the Individual Based Model have been computed and their non-trivial topological structure, shown to persist for a certain time before being destroyed by noise-induced fluctuations. Quantitative estimates of the agreement between the Individual Based Model and the Macroscopic model have been given. This study provides one more evidence of the role of geometry and topology in the emergence of self-organized behavior in active particle systems. The model presented in this article opens many new research directions. Some of them are listed below. 1. The stability of the MO (3.2), HW (3.5) and GS (3.8) solutions as well as those of the generalized HW solutions described in Section SM6 is an open problem. It would enable us to investigate the potential link between topological structure and stability. Numerical simulations have been carried out in a periodic setting. Real systems though are confined by solid walls. To model the influence of confinement, it is necessary to explore wider classes of boundary conditions. Most topological states in physical systems consist of linear perturbations of bulk states that propagate on the edges of the system (edge states). It would be interesting to determine whether linear perturbations of the MO or HW solutions could host such edge states.

Figure 16 :

 16 Figure 16: Different outcomes of the simulation of the IBM starting from perturbed initial MO. Only the four main indicators are shown: from left to right, the GOP, the mean polar angle (or pitch) θ, the mean azimuthal angle (or yaw) φ and the winding number w z . (a)-(d) For σ = 134, the system stays a MO for a long time (t ≃ 20) but eventually converges to a FS; (e)-(h) for σ = 753, the system converges towards a generalized solution with a polar angle not equal to π/2 and a winding number equal to 1 along the z-axis; (i)-(l) for σ = 1000, the MO is quickly disrupted (at t ≃ 5) and converges almost immediately towards a FS. Parameters: N = 3 10 6 , R = 0.025, κ = 10, ν = 40, c 0 = 1, L = 1, ξ = 2π.

Video 7 .Video 8 .Video 9 .Video 10 .Video 11 .Video 12 .

 789101112 It supplements Fig. 13 in Section 5.2. It provides the time-evolution of the RPZ curve for the initial condition MO2 (5.2). The circle formed by the initial RPZ curve immediately opens. The opening width constantly increases, until the arc is reduced to a single point opposite to the opening point at time t ≈ 10. Then there is a bounce and the arc forms again and increases in size until it reaches a maximum and decreases again. Several bounces are observed with decreasing amplitudes. These bounces result in the non-monotonous behavior of the quantity d z displayed on Fig. 13d. It supplements Fig. 14 in Section 5.2. It shows the time-evolution of the particles for the initial condition HW1 (5.4) (see Fig. 3a for details on the representation of the body orientation). For clarity, only a sample of 5000 particles are shown. Before time t ≃ 125, we observe a steady HW state. Then, after time t ≈ 125, the particles show an undulating wave-like behavior, with slowly increasing frequency and amplitude, which causes the decrease of the GOP. Around time t ≈ 178, the particles are divided into two groups with pitch angles θ ≃ 0 and θ ≃ π, which suddenly reverses the global direction of motion. After time t ≈ 178, the particles quickly adopt the same body-orientation. Shortly after time t = 178, the particles still have an undulating behavior but it quickly fades away until a FS is reached. It supplements Fig. 14 in Section 5.2. It shows the time-evolution of the RPXcurve for the initial condition HW1. Unlike in the MO case, the RPX curve does not shrinks to the center of the circle before migrating to its limiting point. In this case, the limiting point near the unit circle towards which the RPX curve is converging attracts the RPX. During this transition, the circular shape of the RPX curve is preserved until it becomes a point. It supplements Fig. 14 in Section 5.2. It shows the time-evolution of the particles for the initial condition HW2 (5.5). For clarity, only a sample of 5000 particles are shown (see Fig.3afor details on the representation of the body orientation). At the beginning, we see two opposite alternations of the three side colors of the tetrahedra (green-blue-magenta followed by green-magenta-blue), which signals a double parallel twist. Then, gradually, the green color is eaten up by the blue and magenta ones and only one alternation of the blue and magenta colors remains. Then the color alternation shades away and gives room to a homogeneous color showing that the body orientations have stopped rolling and a FS is attained. It supplements Fig.14in Section 5.2. It provides the time-evolution of the RPX curve for the initial condition HW2 (5.5). The circle formed by the initial RPX curve immediately opens. The opening width constantly increases, although at a slower pace than for MO2 (see Video 7). Here, also contrasting with the MO2 case, the monotonous opening of the arc results in a monotonously increasing quantity d x as shown in Fig. 14d. It supplements Fig. SM2 in Section SM8.1. It shows the time-evolution of the particles for a MO initial condition (5.1) in a rare case where it evolves into a HW. For clarity, only a sample of 5000 particles are shown (see Fig. 3a for details on the representation of the body orientation). It starts like Video 4 with the ensemble rotation of the particle directions This manuscript is for review purposes only.

Video 13 .Video 14 .

 1314 It supplements Fig.SM2in Section SM8.1. It provides the time-evolution of the RPZ curve for a MO initial condition (5.1) in a rare case where it evolves into a HW. The behavior is essentially the same as in Video 5 except that the RPZ-curve shrinks to a single point far away from the unit circle. This shows that the end state of the RPZ-curve is closer to disorder than for a milling to flocking transition. Before that, the non-trivial topology across z is lost following a similar scenario as for the milling-to-flocking transition. It supplements Fig.SM2in Section SM8.1. It provides the time-evolution of the RPX curve for a MO initial condition (5.1) in a rare case where it evolves into a HW.

Video 15 .Video 16 .

 1516 It supplements Fig. SM3 in Section SM8.2. It shows the time-evolution of the particles for a MO initial condition (5.1) in a rare case where it evolves into a FS through a transient HW. For clarity, only a sample of 5000 particles are shown (see Fig. 3a for details on the representation of the body orientation). The point of view is changed from Video 12 to better visualize the transient HW moving along the diagonal, appearing around time t ≈ 16. At the beginning we witness the ensemble rotation of the particles and its disruption by an instability. After some chaotic behavior, the transient HW establishes as shown by the alternations of blue, green and magenta colors propagating along the diagonal.But after some time, the HW structure is disrupted again and the system eventually establishes a FS. It supplements Fig.SM3in Section SM8.2. It provides the time-evolution of the RPZ curve for a MO initial condition (5.1) in a rare case where it evolves into a FS through a transient HW. The behavior is essentially the same as in Video 5 except that the RPZ-curve undergoes a longer-lasting chaotic dynamics before shrinking to a point which migrates towards the unit circle.

3 .

 3 . Draw a subset I ⊂ {1, . . . , N } of jumping agents: for each agent k ∈ {1 . . . , N }, draw a random number r k uniformly in [0, 1]. If r k > exp(-ν ∆t), then k ∈ I. Compute the local flux: for k ∈ I, compute

SM4.

  Derivation of the macroscopic model. The derivation of the continuum theory presented in Section 2.2 has been achieved in [SM5] (see also [SM2]) following earlier works [SM3, SM4]. It consists of two steps. The first step is the derivation of a mean-field kinetic model in the limit N → ∞ showing that the system satisfies the propagation of chaos property:

SM6. 3 .

 3 Proof of Lemma 3.3. The three rotation matrices are given by A(-ωt, e 3 ) =   cos(ωt) sin(ωt) 0 sin(ωt) cos(ωt) ξ(z -λt))sin(ξ(z -λt)) 0 sin(ξ(z -λt)) cos(ξ(z -λt))   , and a direct computation shows that the three column vectors Ω, u and v of the matrix A ξ,θ θ sin(ξ(z -λt)) cos(ωt) + cos(ξ(z -λt)) sin(ωt) cos θ sin(ξ(z -λt)) sin(ωt) + cos(ξ(z -λt)) cos(ωt) sin θ sin(ξ(z -λt)) θ cos(ξ(z -λt)) cos(ωt)sin(ξ(z -λt)) sin(ωt) cos θ cos(ξ(z -λt)) sin(ωt)sin(ξ(z -λt)) cos(ωt) sin θ cos(ξ(z -λt))r = ξ sin θ cos(ξ(z -λt))uξ sin θ sin(ξ(z -λt))u = ξ sin θ(sin(ωt), cos(ωt), 0) T , δ = cos θ∂ z u • v + u 3 δ z v • Ω = ξ cos θ,where we have used that ∂ z u = ξv and ∂ z v = -ξu. It remains to check that Eq. (2.11b) holds true. We split this equation into three equations, one for each vector Ω, u and v. The first equation on Ω reads(∂ t + c 2 (Ω • ∇ x ))Ω + c 4 P Ω ⊥ r = 0.This equation holds true because∂ t Ω = -ω   sin θ sin(ωt) sin θ cos(ωt) 0   , (Ω • ∇ x )Ω = 0, P Ω ⊥ r = r -(r • Ω)Ω = ξ sin θ ω = c 4 ξ. The second equation on u reads (∂ t + c 2 (Ω • ∇ x ))uc 4 (u • r)Ω + c 4 δv = 0.Because λ = c 2 cos θ, we have∂ t + c 2 Ω • ∇ x = ∂ t + c 2 cos θ∂ z = ∂ t +λ∂ z and ∂ t + λ∂ z (z -λt) = 0. Thus (∂ t + c 2 (Ω • ∇ x ))u = ω   cos θ sin(ξ(z -λt)) sin(ωt) + cos(ξ(z -λt)) cos(ωt) cos θ sin(ξ(z -λt)) cos(ωt)cos(ξ(z -λt)) sin(ωt) 0   , and using ω = c 4 ξ, it can be checked that (∂ t + c 2 (Ω • ∇ x ))uc 4 (u • r)Ω = -c 4 ξ cos θv = -c 4 δv, which yields the result. The equation on v is analogous.

GOP 1

 1 and GOP 2 , these two extreme values being attained respectively when |⟨u⟩| = 0 and |⟨u⟩| = 1.SM7. Convergence rate of |d φ/dt| as N → ∞. The fact that the convergence rate of |d φ/dt| is close to N -1 agrees with previously documented observations in spherical statistics.Indeed, it has been shown in [SM11, Theorem 3(e)] that the estimation of the concentration parameter of a (spherical) von Mises distribution obtained from a crude averaging procedure from N independent samples produces a biased estimator with a (nonnegative) bias of order N -1 (see also [SM10,Section 10.3]). In the present case, a similar reasoning can be applied, which we now briefly develop. The key observation is that all the measured quantities are functions of empirical averages of the form (2.4). Under the chaos assumption (see Section SM4), when N is large, the body-orientations of the particles behave as N independent samples with common law M A , where A solves the SOHB model (2.11) and M A is defined by(2.6). In [SM3, Theorem 4.1], it has been shown that c 4 (κ) can actually be expressed as a function of a certain number p of averaged quantitiesc 4 (κ) = F (⟨g 1 ⟩ M A , . . . , ⟨g p ⟩ M A ),whereg i : SO 3 (R) → M 3 (R) and F : M 3 (R) p → Rare smooth functions. The IBM simulation thus defines an estimator κ of the concentration parameter such that c 4 (κ) = F (ĝ 1 , . . . , ĝp ), where ĝi is the average of g i obtained by replacing M A by the empirical measure of the N body-orientations of the particles. We can then measure the bias by taking the expectation of the Taylor expansion of the previous expression around the point (⟨g 1 ⟩ M A , . . . , ⟨g p ⟩ M A ) : c 4 (κ) = c 4 (κ) + δĝ • ∇F + (δĝ) T (Hess F )δĝ + R, where δĝ = (ĝ 1 , . . . , ĝp ) T -(⟨g 1 ⟩ M A , . . . , ⟨g p ⟩ M A ) T and R is a remainder. The gradient ∇ and Hessian Hess are defined within the Euclidean framework given by (2.1). By the chaos hypothesis E[δĝ] = 0 and by the central limit theorem, the term of order two behaves as N -1 . Since SO 3 (R) is compact, higher order moments of δĝ can be controlled by a classical argument based on Hoeffding's inequality [SM12, Lemma 5.5 and Theorem 5.29]. This ensuresthat E[R] is O(N -2 ).We therefore obtain a biased estimator:E[c 4 (κ)] = c 4 (κ) + a N + O(N -2 ),where a ∈ R depends on the derivatives of the considered functions and on the variance of the estimator(2.4) where the particles are replaced by independent identically distributed This manuscript is for review purposes only.samples with law M A . The fact that a > 0 can be empirically verified on Fig.8bbut has not been proved yet. For each N , the fluctuations around the average (biased) value can be monitored by computing the standard deviation of the 10 independent simulations. Fig.SM1shows this standard deviation as a function of N in a log-log-scale (blue dots). Although fluctuations remain significant with only 10 simulations per data point, by a standard linear regression (solid orange line) we obtain that the size of the standard deviation behaves as N -β with β ≃ 0.54. which is close to the value β = 1/2 which we expect from an application of the central limit theorem.

Figure SM1 :

 SM1 Figure SM1: Standard deviation of the 10 independent simulations as a function of N (blue dots) and regression line (solid orange line) in log-log scale. Parameters: L = 1, ξ = 2π, R = 0.025, ν = 40, c 0 = 1, κ = 10.

  and θ. The topological indicators (d z , rz , w z ) shown in the second line of figures of Fig.SM3

Figure SM2 :

 SM2 Figure SM2: Transition from a MO to a HW: example of a solution of the IBM for an initial condition sampled from (5.1) in the rare case where it leads to a HW. The following indicators are plotted as functions of time: (a) GOP (b) Pitch θ of Ω. (c) Yaw φ of Ω. (d) Distance of center of mass of RPZ curve to the origin d z . (e) Mean distance of RPZ curve to the origin rz . (f) Winding number of RPZ curve w z . (g) Distance of center of mass of RPX curve to the origin d x . (h) Mean distance of RPX curve to the origin rx . (i) Winding number of RPX curve w x . Gray shaded zones highlight a small region around the time of minimal GOP. Parameters: N = 1.5 • 10 6 , R = 0.025, L = 1, D = 0.1, ν = 40, c 0 = 1. See caption of Fig. 13 for further indications. See also Videos 12 to 14 in Section SM1.

Figure SM3 :

 SM3 Figure SM3: Transition from a MO to a FS via an unstable HW: example of a solution of the IBM for an initial condition sampled from (5.1) in the rare case where it leads to a FS through a transient HW. The following indicators are plotted as functions of time: (a) GOP (b) Pitch θ of Ω. (c) Yaw φ of Ω. (d) Distance of center of mass of RPZ curve to the origin d z . (e) Mean distance of RPZ curve to the origin rz . (f) Winding number of RPZ curve w z . Gray shaded zones highlight a small region around the time of minimal GOP. Parameters: N = 1.5 • 10 6 , R = 0.025, L = 1, D = 0.1, ν = 40, c 0 = 1. See caption of Fig. 13 for further indications. See also Videos 15 and 16 in Section SM1.

Table 1 :

 1 Parameters of the IBM (2.3), (2.5).

	Parameter	Symbol
	Number of particles	N
	Computational box side length	L
	Interaction radius	R
	Particle speed	c 0
	Concentration parameter	κ
	Alignment frequency	ν

This manuscript is for review purposes only.

Acknowledgments. Part of this research was done when PD was affiliated to Department of Mathematics, Imperial College London, London, SW7 2AZ, United Kingdom. Data statement:. No new data were collected in the course of this research.

PD acknowledges support by the Engineering and Physical Sciences Research Council (EPSRC) under grants no. EP/M006883/1 and EP/P013651/1, by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048. The work of AD is supported by an EPSRC-Roth scholarship cofounded by the Engineering and Physical Sciences Research Council and the Department of Mathematics at Imperial College London.

We refer to Fig. 3a for details on the representation of the body orientation using four-colored tetrahedra. We notice the ensemble rotation of the particle directions about the z axis until an instability disrupts the body orientation twist along the z axis (around time t ≈ 13) and eventually drives the system to a FS. where its radius shrinks down. Then, the RPZ-curve shows a fairly chaotic dynamics during which the topology is lost. This happens around time t ≈ 13 which is the first time when the RPZ-curve passes through the origin; at this time, the winding number is not defined. Then, the RPZ-curve slowly migrates towards the unit circle while shrinking to a single point which signals a FS. From time t ≈ 15 on, it remains a single immobile point. † Institut de Mathématiques de Toulouse; UMR5219; Université de Toulouse; CNRS; UPS; F-31062 Toulouse Cedex 9, France (pierre.degond@math.univ-toulouse.fr, https://sites.google.com/site/degond/Home).

‡ Department of Mathematics, Imperial College London, London, SW7 2AZ, United Kingdom (antoine.diez18@imperial.ac.uk, mingye.na18@imperial.ac.uk).

SM1

This manuscript is for review purposes only.

out of the integral with respect to A. These considerations naturally lead to the choice of the antisymmetric matrix ψ A f ε (A) as a multiplier. Because A f ε is obtained as the polar decomposition of J f ε , there exists a symmetric matrix S such that J f ε = A f ε S. Using this remark and (4.4), we easily find that

Then, multiplying (SM4.1) by ψ f ε , taking the limit ε → 0 and assuming (SM4.2) leads to:

Eq. (2.11b) of the SOHB model follows from this equation through tedious but straightforward computations detailed in [SM3, SM4].

Note that the simple form of the multiplier ψ A f is due to a particular simple expression of the collision operator. In more general cases, the obtention of the multiplier (referred to as the generalized collision invariant in [SM6]) is more involved (see e.g. [SM3, SM4, SM5]).

A rigorous convergence result for the limit ε → 0 is not available to date. In the case of the Vicsek model, such a rigorous result has been proved in [SM8].

SM5. Alternate expressions of δ. The following lemma provides alternate expressions for δ:

Lemma SM5.1. We have

Proof. Eq. (SM5.1) follows from inserting the formula

and similar formulas after circular permutation of {Ω, u, v}into (2.13). Eq. (SM5.2) follows from taking the half sum of (2.13) and (SM5.1) and applying the formula

and similar formulas after circular permutation of {Ω, u, v}. SM6. MO, HW, GS and generalized HW solutions. In this section, we provide proofs of Lemmas 3.1, 3.2 and 3.3. The prototypical helical traveling wave (HW) presented in Lemma 3.2 belongs to a more general class of solutions called generalized HW solutions described in Section SM6.2 below.

This manuscript is for review purposes only.

alternately stronger and weaker, generate this oscillatory behavior. Note that a transition to a HW cannot occur when the global direction of motion at the transition time is not one of the principal axes of the square domain since a HW along another direction is not compatible with the periodic boundary conditions (see Section SM8.2). This is confirmed by the final values of φ and θ (both equal to π/2) which correspond to a global direction of motion oriented along the y-axis (in what follows, in reference to (5.4) and to avoid confusion, we will still call that direction, the x direction).

The second and third lines of figures in Fig. SM2 show the triplets of topological indicators (d z , rz , w z ) and (d x , rx , w x ) which materialize the MO and HW structures respectively. The mean distance of the RPZ-curve to the origin rz (Figs. SM2e) decreases, revealing an increase of the disorder. Simultaneously, the distance of its center of mass to the origin d z increases (Figs. SM2d) showing a transition trend to a FS. The winding number w z (Fig. SM2f) jumps from 1 to 0 at the time of maximal disorder. However, d z and rz do not reach zero, showing that complete disorder across z is not reached. Since the final state of the system is a generalized helical wave state (see Section SM6.4), we do not necessarily expect that complete disorder will be reached along the z-direction. In the mean time, rx starts from 0 (complete disorder) and increases up to a value close to unity, showing the build-up of a HW.

The quantity d x increases during some time but eventually decreases to 0 (not shown in the figure) as it should for a HW. Finally, the winding number w x is undefined in the initial stage, as it should for complete disorder, but builds up to 1 at the time where the winding number w z drops to 0. There is a transfer of non-trivial topology from an MO structure to a HW structure. SM8.2. From milling to flocking via a helical wave state. In some rare cases an intermediate unstable HW can be observed. Note that due to the periodic setting, an HW cannot be stable for most of the the global directions of motion. Although stable or unstable HW typically appear in one over twenty of our simulations, it should be kept in mind that the occurrence frequency also depends on the geometry of the domain and that this phenomena may be more frequent for other simulation settings. The procedure is the same as in the previous section. Fig. SM3 shows the results (see also supplementing videos 15 and 16 in Section SM1).

The transition stage between the MO and FS is significantly longer than in the previous situations. During that phase, the GOP (Fig. SM3a) oscillates between the value Ψ 1 characterizing the MO and lower values, i.e. lower order. Likewise, there are significant variations of the pitch θ (Fig. SM3b) and yaw φ (Fig. SM3c). As in the previous section, this could be explained by antagonist effects of different attractors (the MO and HW) and subsequent oscillations of the system between them. Video 15 reveals large scale band structures similar to a HW except that the global direction of motion is not one of the principal axes of the square domain. As, in most cases, this cannot be compatible with the periodic boundary conditions, such state cannot persist in time. The relatively long-time persistence of this stage could be explained in the present case by the fact that the global direction of motion seems to oscillate around the direction given by e 1 + e 2 (i.e. φ = π/4 and θ = π/2) which is theoretically compatible with the periodic boundary conditions, provided the wave length ξ is changed from 2π/L to √ 2π/L. This state does not seem to be stable as shown by the large oscillations of φ

This manuscript is for review purposes only.