Bulk topological states in a new collective dynamics model
BBulk topological states in a new collectivedynamics model
Pierre
Degond , Antoine Diez , and Mingye Na [email protected] [email protected] [email protected] Institut de Math´ematiques de Toulouse ; UMR5219,Universit´e de Toulouse ; CNRS,UPS, F-31062 Toulouse Cedex 9, France
Department of Mathematics, Imperial College London, South Kensington Campus,London, SW7 2AZ, United Kingdom
Abstract
In this paper, we demonstrate the existence of topological states in a new collec-tive dynamics model. This individual-based model (IBM) describes self-propelledrigid bodies moving with constant speed and adjusting their rigid-body attitude tothat of their neighbors. In previous works, a macroscopic model has been derivedfrom this IBM in a suitable scaling limit. In the present work, we exhibit explicit so-lutions of the macroscopic model characterized by a non-trivial topology. We showthat these solutions are well approximated by the IBM during a certain time butthen the IBM transitions towards topologically trivial states. Using a set of appro-priately defined topological indicators, we reveal that the breakage of the non-trivialtopology requires the system to go through a phase of maximal disorder. We alsoshow that similar but topologically trivial initial conditions result in markedly dif-ferent dynamics, suggesting that topology plays a key role in the dynamics of thissystem.
Keywords: individual-based model, macroscopic model, self-organization, topologi-cal phase transition, winding number, order parameter
AMS subject classification:
Acknowledgements:
Part of this research was done when PD was affiliated to De-partment of Mathematics, Imperial College London, London, SW7 2AZ, United a r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n ingdom. PD acknowledges support by the Engineering and Physical Sciences Re-search Council (EPSRC) under grants no. EP/M006883/1 and EP/P013651/1, bythe Royal Society and the Wolfson Foundation through a Royal Society WolfsonResearch Merit Award no. WM130048. Data statement: no new data were collected in the course of this research.
Contents N → ∞ . . 213.3.2 Quantitative comparison between the models . . . . . . . . . 243.4 Topology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 List of supplementary videos 42B Quaternion framework 44C Numerical methods 45D Derivation of the macroscopic model 46E MO and generalized HW solutions 48
E.1 Proof of Lemma 3.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 48E.2 Generalized HW and proof of Lemma 3.2 . . . . . . . . . . . . . . . 50E.3 GOP of the MO and generalized HW . . . . . . . . . . . . . . . . . . 50
References 50
Systems of particles (or agents) which exhibit self-organized collective behavior areubiquitous in the living world at all scales, from bird flocks [70] to sperm [27] orbacterial colonies [29]. Examples are also found in social sciences [18, 39] or for inertmatter [15]. In such systems, the agents interact locally with a limited number ofneighbors through rather simple rules such as attraction, repulsion or alignment[3, 26, 51] without any leader or centralized control. When the number of agentsbecomes large, vast structures encompassing many agents appear, such as clusters[72, 89], traveling bands [23], vortices [24, 29], lanes [25], etc. As there is no direct orapparent 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. [89] for a review).There are mainly two levels of description of particle systems: the most detailedone consists of individual based models (IBM) where the agents dynamics are de-scribed by coupled ordinary or stochastic differential equations. When the numberof agents becomes large, a macroscopic description in terms of average quantitiessuch as the agents mean density or velocity is preferred. The rigorous link betweenthese two levels of description involves two successive limits by which the numberof agents is first sent to infinity (mean-field limit) and then, the system size rela-tive to the typical interaction distance between the agents is also sent to infinity(hydrodynamic limit), see e.g. [21, 31]. In collective dynamics, particles are ca-pable of self-propulsion by transforming an internal source of chemical energy intomotion [89]. There are two main classes of IBM of self-propelled particles. Thefirst class is based on the Cucker-Smale model [4, 28, 54, 55] where self-propulsionis treated as an external force. The second class is based on the Vicsek model[2, 19, 23, 29, 41, 44, 72, 88] where self-propulsion is modeled by imposing the normof 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-Smaletype models and [10, 34, 41, 44, 78, 86] for Vicsek type models). The two categoriesare linked by an asymptotic limit [12, 13]. Of course, there are many variants ofthese models and we refer to [8, 9, 17, 20, 42, 45, 74] for a non-exhaustive set ofexamples. ecently, a series of studies has investigated the existence of topological states incollective dynamics. Topological states have appeared with the quantum Hall effect[66, 68, 75, 85] which relies on so-called conducting chiral edge states: when a sam-ple of a 2-dimensional insulator is placed in a magnetic field, its bulk conductance isnil 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 magneticfield have been discovered, the so-called “topological insulators” [57, 76, 79]. Chiraledge states are robust against perturbations because of their non trivial topologywhich can be characterized by a integer, the winding number. Any destruction ofthe chiral edge state would require a finite jump of this integer, which consumes afinite amount of energy. Hence lower energy perturbations will fail to destroy thechiral edge state. This property is of strategic interest for various applications suchas quantum computers. Recently a series of works have explored the occurrenceof topological states in collective dynamics (see e.g. [82–84]). They are based onnumerical simulations of the Toner and Tu model [86], which is a continuum analogof the Vicsek model [88]. Investigating appropriate geometrical configurations (asphere in [82], a network of rings in [83, 84]), they show that linearized perturba-tions of the stationary state (i.e. sound waves) generate chiral edge states whichpropagate unidirectionally, revealing an underpinning non-trivial topology. How-ever, 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 thenoise induced by this finite particle number long enough is not investigated.In this paper, we demonstrate the existence of non-trivial bulk topological statesin a new collective dynamics model. Bulk states propagate in the whole domain,by opposition to edge states which are localized at the boundary. The collectivedynamics model studied here has first been proposed in [35] and later analyzed andexpanded in [32, 37, 38]. Referred to below as the “Body-Alignment Individual-Based Model” (BA-IBM or IBM for short), it describes self-propelled rigid bodiesmoving with constant speed and trying to adjust their rigid body attitude to thatof their neighbors. In [35, 37] the BA-IBM was based on Stochastic DifferentialEquations (SDE) and a macroscopic model named the “Self-Organized Hydrody-namics for Body-orientation (SOHB)” was derived. In [32, 38], SDE were replacedby Piecewise Deterministic Markov Processes (PDMP) in the IBM but the macro-scopic model remained the SOHB model (with possibly different coefficients). In[32], a variant of the BA-IBM was shown to exhibit phase transitions which wererigorously studied. In the present work, we derive explicit solutions of the SOHBmodel which exhibit striking non-trivial topologies revealed by non-zero windingnumbers. We explore how these non-trivial topologies are maintained at the levelof the IBM by solving the PDMP of [38]. In particular, we observe that, due tonoise induced by the finite particle number, topological phase transitions from stateswith non-trivial topology to states with trivial one may occur and we study thesephase transitions in detail. Using a set of appropriately defined topological indica-tors, we reveal that the breakage of the non-trivial topology requires the system togo through a phase of maximal disorder. We also show that similar but topologi-cally trivial initial conditions result in markedly different dynamics, suggesting thattopology plays a key role in the dynamics of this system. We are led to questionthe possible existence of topological protection against perturbations as mentioned bove for topological insulators. Compared to previous works on topological statesin collective dynamics, we deal with bulk states instead of edge states and we explorethem at the level of the IBM and not just at the continuum level, which is closer torealistic particle systems. The present work adds a new item to the list of collectivedynamics models exhibiting topological states. The topological protection conceptcould bring new perspectives to poorly understood questions such as the robustnessof morphogenesis or the emergence of symmetries in growing organisms.The present model belongs to the category of Vicsek-like models in the sense thatit 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 unitsphere (after convenient normalization). In the present IBM, the particles carryan orthonormal frame, or equivalently, a rotation matrix, that describes their bodyattitude. Thus their degrees of freedom are constrained to belong to the manifoldSO ( R ) of 3 × not aligned alignedPolar alignment Body−orientationnot aligned aligned Figure 1:
Vicsek model versus body-alignment model. Left: polar alignmentof velocity orientations (red vectors) of two agents. Right: alignement of body-orientations: in addition to its velocity orientation (red), each agent has two otheraxes (green and blue), the three vectors forming a direct orthogonal frame.We complete this introduction by a review of the mathematical literature onthe Vicsek model and the BA-IBM. The mean-field limit of the IBM has beenproven in [10] for the Vicsek model and in [43] for the body orientation model.Existence theory for the mean-field Vicsek model is available in [14, 47, 50] butthe 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[33, 34, 48] and [32] for the Vicsek and body orientation models respectively. Thenumerical approximation of the mean-field kinetic model has been undertaken forthe Vicsek model only in [49, 53]. The derivation of macroscopic equations fromthe mean-field Vicsek kinetic equations has first been formally achieved in [41] andlater rigorously proved in [64]. Corresponding works for the body alignment model re only formal [35, 37, 38]. Existence theory for the hydrodynamic models derivedfrom the Vicsek model can be found in [40, 90] and numerical methods in [44, 73].Both questions are still open for the body orientation model.The organization of this paper is as follows. Section 2 is devoted to the exposi-tion of the IBM and macroscopic models. Then explicit solutions of the macroscopicmodel are derived in Section 3 and are shown to exhibit non-trivial topology. Theyalso serve as benchmarks to show that the macroscopic model is an accurate ap-proximation of the IBM. But after a some time, the IBM departs from the specialsolutions of the macroscopic model and undergoes a topological phase transition.The study of these phase transitions require appropriate topological indicators whichare developed in Section 4. Then, the topological phase transitions are analyzed inSection 5. A discussion and some open questions raised by these observations canbe found in Section 6. Appendices collect various supplementary information: a listof supplementary videos (Appendix A), a summary of the quaternion framework(Appendix B), a description of the numerical methods (Appendix C), a summaryof the derivation of the macroscopic models (Appendix D) and finally a derivationof the explicit solutions presented in Section 3 (Appendix E). In this section, we present the Individual-Based body-alignment Model (IBM). Thismodel was first proposed in [38]. We consider N particles (or individuals, or agents)indexed by k ∈ { , . . . , N } whose spatial locations are denoted by X k ( t ) ∈ R where t ∈ [0 , ∞ ) is the time. A direct orthonormal frame { Ω k ( t ) , u k ( t ) , v k ( t ) } isattached to each particle (i.e. Ω k , u k , v k ∈ S , Ω k · u k = 0 and v k = Ω k × u k ).Likewise, if ( e , e , e ) is a fixed direct orthonormal reference frame, we define A k ( t )to be the unique element of the special orthonormal group SO ( R ) which maps( e , e , e ) onto (Ω k ( t ) , u k ( t ) , v k ( t )). We will choose ( e , e , e ) once for all andwrite A k ( t ) = [Ω k ( t ) , u k ( t ) , v k ( t )]. This will be referred to as the local particle frameor 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 exceptat exponentially distributed times at which the local frame jumps and adjusts itselfto the average neighbors’ local frame up to some noise. The motion of the particlesis thus described by the functions [0 , ∞ ) (cid:51) t (cid:55)→ ( X k ( t ) , A k ( t )) ∈ R × SO ( R ) for k ∈ { , . . . , N } .We first describe how the average neighbors’ local frame is defined. We introducea fixed observation (or sensing) kernel K : R (cid:51) x (cid:55)→ K ( x ) ∈ [0 , ∞ ). We assumethat K is a radial function (i.e. there exists ˜ K : [0 , ∞ ) (cid:51) r (cid:55)→ ˜ K ( r ) ∈ [0 , ∞ ) suchthat K ( x ) = ˜ K ( | x | ), where | x | is the euclidean norm of x ). For a collection of N particles { ( X k , A k ) } k ∈{ ,...,N } ∈ ( R × SO ( R )) N , we define the local flux as the ollowing 3 × J k = 1 N N (cid:88) 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 zerowith 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 factthat e.g. distant particles will contribute to J k less than neighboring particles. Ingeneral, J k is not a rotation matrix. To recover a rotation matrix, we need to map J k back onto the manifold SO ( R ). To do so, the space M ( R ) of 3 × A · B := 12 Tr( A T B ) , (1)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: A k := arg max A ∈ SO ( R ) A · J k . (2)This expression stands for the element A k ∈ SO ( R ) that maximizes the functionSO ( R ) (cid:51) A (cid:55)→ A · J k ∈ R . The maximization procedure (2) has a unique solutionas soon as J k is not singular, i.e. det J k (cid:54) = 0 where det stands for the determinant.Since the singular matrices form a zero-measure set in M ( R ) it is legitimate toassume that, except for a zero-measure set of initial data, this situation will notoccur. Furthermore, when det J k > A k is nothing but the unique rotation matrixinvolved in the polar decomposition of J k .We let the particles evolve according to the following Piecewise DeterministicMarkov Process (PDMP). • To each agent k ∈ { , . . . , N } is attached an increasing sequence of randomtimes (jump times) T k , T k , . . . such that the intervals between two successivetimes are independent and follow an exponential law with constant parameter ν > T nk , the function X k is continu-ous and the function A k has a discontinuity between its left and right statesrespectively denoted by A k ( T nk −
0) and A k ( T nk + 0). • Between two jump times ( T nk , T n +1 k ), the evolution is deterministic: the orien-tation of Agent k does not change and it moves in straight line at speed c > A k ( T nk + 0) e , i.e. for all t ∈ [ T nk , T n +1 k ), we have X k ( t ) = X k ( T nk ) + c ( t − T nk ) A k ( t ) e , A k ( t ) = A k ( T nk + 0) . (3) • To compute A k ( T nk + 0) from A k ( T nk − T nk − J n − k := 1 N N (cid:88) j =1 K (cid:0) X k ( T nk ) − X j ( T nk ) (cid:1) A j ( T nk − , (4) aving in mind that A j ( T nk −
0) = A j ( T nk ) for j (cid:54) = k . From J n − k , which weassume is a non-singular matrix, we compute A nk as the unique solution of themaximization problem (2) (with J k replaced by J n − k ). Then, A k ( T nk + 0) isdrawn from a von Mises distribution: A k ( T nk + 0) ∼ M A nk . (5)The von Mises distribution on SO ( R ) with parameter A ∈ SO ( R ) is definedto be the probability density function: M A ( A ) := e κ A · A (cid:82) SO ( R ) e κ A · A (cid:48) d A (cid:48) , (6)where κ > ( R )) of the Gaussian distribution in a flat space. The new orientation ofAgent k at time T n can therefore be interpreted as a small random perturbationof the average local orientation given by A nk , where the perturbation size ismeasured by 1 / √ κ .In Formula (6) and in the remainder of this paper, the manifold SO ( R ) isendowed with its unique normalized Haar measure defined for any test function ϕ by: (cid:90) SO ( R ) ϕ ( A ) d A := 2 π (cid:90) π (cid:90) S ϕ ( A ( θ, n )) sin ( θ/
2) d θ d n , (7)where d n is the uniform probability measure on the sphere S . Here, a rotationmatrix A ≡ A ( θ, n ) is parametrized by its rotation angle θ ∈ [0 , π ] and its axis n ∈ S through Rodrigues’ formula: A ( θ, n ) := I + sin θ [ n ] × + (1 − cos θ ) [ n ] × = exp( θ [ n ] × ) (8)with n = ( n , n , n ) T and I is the 3 × w =( w , w , w ) T ∈ R , [ w ] × is the antisymmetric matrix of the linear map R (cid:51) u (cid:55)→ w × u (where × denotes the cross product) which has the following expression:[ w ] × := − w w w − w − w w . (9)Additional details on the structure of SO ( R ) can be found for instance in [63]. TheIBM (3), (5) is schematically represented in Fig. 2. Unless otherwise specified, throughout this paper, a square box of side length L withperiodic boundary conditions is used. As sensing kernel K , we use the indicatorfunction of the ball centered at 0 and of radius R . Thus, an agent interacts with allits neighbours at a distance less than R (radius of interaction). Table 1 summarizesthe model parameters. nk regiontime T n − k time T nk A k ( T nk + 0) A k ( T nk − Observationtime T n − k Figure 2:
Schematic representation of the PDMP described in the text: the motionof Particle k is represented in physical space as the black broken dotted line. Thebody frame A k is represented with Ω k in red, u k in green and v k in blue. Eachangular point of the trajectory corresponds to one of the jump times T nk . Betweentwo jump times, the trajectory is the straight line spanned by Ω k and the bodyframe stays constant. The jump dynamics is depicted at time T nk . At this time, theobservation region is colored in yellow and body frames of the other particles presentin this region are depicted in light blue. The averaged body frame A nk is depictedwith thick lightly colored arrows. The body frame before the jump A k ( T nk −
0) isdrawn in broken lines whereas that after the jump A k ( T nk +0) is drawn in plain lines. A k ( T nk + 0) is close, but not equal to A nk because of the noise intensity proportionalto 1 /κ . For clarity, the frames involved in the description of the jump are magnified. Parameter SymbolNumber of particles N Computational box side length L Interaction radius R Particle speed c Concentration parameter κ Alignment frequency ν Table 1:
Parameters of the IBM (3), (5).For the numerical simulations presented in this paper, we have used the con- enient framework offered by quaternions. Indeed, there is a group isomorphismbetween SO ( R ) and H / {± } where H is the group of unit quaternions. We can ex-press the IBM (3), (5) using this representation (see [38] and Appendix B). Roughlyspeaking, body-alignment as described here is equivalent to nematic alignment ofthe corresponding quaternions (nematic alignment of a unit quaternions q to themean direction Q is unchanged if q is replaced by − q , as opposed to polar align-ment where the result depends on the sign of q ). This is because a given rotationcan be represented by two opposite quaternions and thus, the outcome of the align-ment process should not depend of the choice of this representative. The numericalalgorithm is described in Appendix C. Additionally, the quaternion framework alsosuggests to use order parameters derived from nematic alignment dynamics (suchas in liquid crystal polymers). We shall use this analogy to define appropriate orderparameters in Section 4.1.All the simulations were written in Python using the PyTorch [77] and KeOps[22] libraries. The computational details as well as the source code will be part ofa forthcoming work. The outcomes of the simulations were analyzed and plottedusing the NumPy [56] and Matplotlib [62] libraries. The 3D particle plots wereproduced using VPython [80].A typical outcome of the IBM is shown in Figure 3 (see also Appendix A, Video1) for a moderate number of particles ( N = 3000). Throughout this paper, in theplots, we will represent each agent graphically by an elongated tetrahedron pointingin the direction of motion. The three large faces around the height will be paintedin 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 systemself-organizes in small clusters (Fig. 3c) and finally reaches a flocking equilibriumwhere all the agents have roughly the same body-orientation (Fig. 3d). We willsee below that flocking is not necessarily the ultimate fate of the system, because itmay be trapped in a so-called topologically protected state. To better understandthese aspects, we first need to develop the continuum (or macroscopic) descriptionof the system. This is done in the next section. We finally make a comparison with previous models. First, there is a version ofthe IBM where particles follow a stochastic differential equation (SDE) instead ofa jump process [35, 37]. Both the current and previous models have the samehydrodynamic model as macroscopic limit (see forthcoming section). There are tworeasons for us to prefer the jump process. First, its simulation is slightly easier andsecond, the coefficients of the macroscopic model are explicit, which is not so in theSDE case where they require the resolution of an auxiliary elliptic problem [35, 37].Beyond the present body-orientation model, numerous models of self-propelledparticles have been proposed in the literature (see the review [89]). The mostclosely related one is the celebrated Vicsek model [88]. There are several versionsof this model: time-discrete ones [23, 88], time-continuous ones relying on an SDEdescription of the particle trajectories [41] and time-continuous ones using a jumpprocess instead [44]. The latter version is the most closely related to the presentwork. In [44], the difference is that particles carry a single direction vector Ω k k e A k e (a) Graphical represen-tation of particles (b)
Time=0 (c)
Time=4 (d)
Time=40
Figure 3: (a) Graphical representation of particles and their body orientationsas 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 outputof the simulation at three different times (b) Time=0, (c) Time=4 and (d) Time=40.Parameters: N = 3000, L = 1, R = 0 . κ = 20, ν = 5, c = 0 .
2. see alsoAppendix A, Video 1.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 expo-nential law with constant parameter ν >
0. At jump times, the position iscontinuous and the direction vector Ω k is discontinuous with left and rightstates respectively denoted by Ω k ( T nk −
0) and Ω k ( T nk + 0). • Between two jump times T nk , T n +1 k , the direction vector Ω k does not changeand the particle moves in straight line at speed c > y Ω k ( T nk + 0). • To pass from Ω k ( T nk −
0) to Ω k ( T nk + 0), we compute the local flux J n − k = N (cid:80) Nj =1 K (cid:0) X k ( T nk ) − X j ( T nk ) (cid:1) Ω j ( T nk − ∈ R and, assuming that it is non-zero, the mean direction ¯Ω nk = J n − k / | J n − k | ∈ S at time T nk −
0. Then, Ω k ( T nk +0)is drawn from a von Mises distribution on S : Ω k ( T nk + 0) ∼ ˜ M ¯Ω nk , with˜ M ¯Ω (Ω) = e κ (¯Ω · Ω) / (cid:82) S e κ (¯Ω · Ω) dΩ, for Ω and ¯Ω in S .So, the current model is an elaboration of [44] replacing self-propulsion directionsby whole body frames and polar alignment of unit vectors (as expressed by the vonMises distribution on the sphere) by alignment of rotations matrices. Outcomes ofnumerical simulations of the Vicsek model do not show striking differences whetherone uses any of the above mentioned versions (time-discrete, time-continuous withSDE or time-continuous with jump process). Results given in [23, 88] for the time-discrete version display the emergence of a global alignment together with the for-mation of clusters when the noise intensity 1 /κ is not too big. The outcome stronglyresembles what is shown in Fig. 3 for the body-orientation model, but for the de-piction 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 orientationinstead of the mere self-propulsion direction makes any change in the dynamics ofthe particle positions and direction vectors. In particular, do the particle positionsand 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, strik-ing differences between the two models are obtained. To show this, the use of themacroscopic limit of the IBM, as developed in the forthcoming section, will be ofcrucial importance. As soon as N is not very small, the IBM (3), (5) involves a large number of unknownswhich makes its mathematical analysis virtually impossible. A reduced description,more amenable to mathematical analysis, is obtained through the macroscopic limitof the IBM, and consists of a system of partial differential equations. This reduceddescription gives a valid approximation of the IBM in an appropriate range of pa-rameters, namely N (cid:29) , RL ∼ c ν L (cid:28) . (10)Throughout the remainder of this paper, we will focus on this regime. The macro-scopic limit of the IBM (3), (5) has first been proposed in [38] and leads to a modelcalled “Self-Organized Hydrodynamics for Body orientation (SOHB)”. The deriva-tion relies on earlier work [35, 37]. This derivation is “formally rigorous” in thesense that, if appropriate smoothness assumptions are made on the involved math-ematical objects, the limit model can be identified rigorously as being the SOHB.For the reader’s convenience, we summarize the main steps of this mathematicalresult in Appendix D.The unknowns in the SOHB are the particle density ρ ( t, x ) and mean body-orientation A ( t, x ) ∈ SO ( R ) at time t and position x = ( x, y, z ) ∈ R . They satisfy he following set of equations: ∂ t ρ + c ∇ x · ( ρ A e ) = 0 , (11a) (cid:0) ∂ t + c ( A e ) · ∇ x (cid:1) A + (cid:2) ( A e ) × ( c ∇ x log ρ + c r ) + c δ A e (cid:3) × A = 0 . (11b)The quantities r and δ have intrinsic expressions in terms of A [35]. However, it ismore convenient to write the rotation field A in terms of the basis vectorsΩ = A e , u = A e , v = A e . With these notations, the vector field r ( t, x ) ∈ R and the scalar field δ ( t, x ) ∈ R are defined by r := ( ∇ x · Ω) Ω + ( ∇ x · u ) u + ( ∇ x · v ) v , (12) δ := [(Ω · ∇ x ) u ] · v + [( u · ∇ x ) v ] · Ω + [( v · ∇ x )Ω] · u . (13)Here, for a vector field B ( x ) ∈ R 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 λ andwe set ( B · ∇ x ) λ = B · ∇ x λ with · the inner product of vectors in R . We remindthat × denotes the cross product and we refer to formula (9) for the definition of[ w ] × when w is a vector in R .The quantities c , c , c , c are functions of κ and c given as follows: c c = 23 (cid:10)
12 + cos θ (cid:11) exp ( κ ( +cos θ )) sin ( θ ) , (14) c c = 15 (cid:104) θ (cid:105) exp ( κ ( +cos θ )) sin ( θ ) cos ( θ ) , (15) c c = 12 κ , (16) c c = 15 (cid:104) − cos θ (cid:105) exp ( κ ( +cos θ )) sin ( θ ) cos ( θ ) , (17)where, for two functions f and g : [0 , π ] → R , we write (cid:104) f (cid:105) g = (cid:82) π f ( θ ) g ( θ ) d θ (cid:82) π g ( θ ) d θ . Fig. 4 provides a graphical representation of these functions.The following lemma provides alternate expressions for δ : Lemma 2.1.
We have δ = − (cid:8) [( u · ∇ x ) Ω] · v + [( v · ∇ x ) u ] · Ω + [(Ω · ∇ x ) v ] · u (cid:9) (18)= − (cid:8) ( ∇ x × Ω) · Ω + ( ∇ x × u ) · u + ( ∇ x × v ) · v } . (19) Proof.
Eq. (18) follows from inserting the formula0 = ∇ x (Ω · u ) = (Ω · ∇ x ) u + ( u · ∇ x )Ω + Ω × ( ∇ x × u ) + u × ( ∇ x × Ω) , and similar formulas after circular permutation of { Ω , u , v } into (13). Eq. (19)follows from taking the half sum of (13) and (18) and applying the formula ∇ x × v = ∇ x × (Ω × u ) = ( ∇ x · u ) Ω − ( ∇ x · Ω) u + ( u · ∇ x )Ω − (Ω · ∇ x ) u , and similar formulas after circular permutation of { Ω , u , v } . Dimensionless coefficients c i /c as functions of the inverse of concentrationparameter 1 /κ . Blue curve c /c , orange curve c /c , green curve c /c and red curve c /c . At the crossover value κ ∗ (cid:39) .
58, the sign of c − c changes (see Section 3.2). To better understand what the SOHB system (11) does, we re-write it as follows: ∂ t ρ + c ∇ x · ( ρ Ω) = 0 , (20a) D t A + [ w ] × A = 0 , (20b)where the convective derivative D t and the vector w are given by: D t = ∂ t + c Ω · ∇ x , (21) w = Ω × F + c δ Ω , with F = − c ∇ x log ρ − c r , (22)Eq. (20a) is the mass conservation equation of the fluid. The vector Ω gives thedirection of the fluid motion. The fluid velocity deduced from (20a) is c Ω. Since c /c ∈ [0 ,
1] as can be seen from Fig. 4 (see also [35] for a rigorous proof), the fluidmotion is oriented positively along Ω and its magnitude is smaller than the particlesself-propulsion velocity c . This is because the average of vectors of identical normshas smaller norm. The quantity c /c can be seen as an order parameter [32] butwe will not dwell on this issue here.Eq. (20b) provides the rate of change of A with time along the integral curvesof the vector field c Ω as expressed by the convective derivative D t . Note that thisvector field is not the fluid velocity c Ω since c (cid:54) = c . It can be interpreted as thepropagation velocity of A when w is zero. Since D t A is the derivative of an elementof SO ( R ), it must lie in the tangent space to SO ( R ) at A which consists of allmatrices of the form W A with W antisymmetric. This structure is indeed satisfiedby Eq. (20b) since, from the definition (9), the matrix [ w ] × is antisymmetric. Itcan be shown that the SOHB system is hyperbolic [36].In fact, Eq. (20b) shows that the vector w is the instantaneous rotation vectorof the frame A ( t, X ( t )), where t (cid:55)→ X ( t ) is any solution of d X d t = c Ω( t, X ( t )).Indeed, Eq. (20b) can be equivalently written as a system of equations for (Ω , u , v ) f the form D t Z = w × Z , with Z = Ω , u , v . This describes a rigid body rotationof the frame { Ω , u , v } with angular velocity w . The rotation vector w has twocomponents. The first one is Ω × F and tends to relax Ω towards F . Due to itsexpression (22), the force F includes two contributions: that of the pressure gradient − c ∇ x log ρ and that of gradients of the body orientation through the vector − c r .The second component of the rotation vector is − c δ Ω and corresponds to a rotationof the body frame about the self propulsion direction Ω driven by gradients of thebody orientation through the scalar − c δ . The contributions of gradients of bodyorientation in the two components of the rotation vector are under the control ofthe single coefficient c . Fig. 5 gives a graphical representation of the actions ofthese two infinitesimal rotations. Fu Ω v ′ v Ω × Fu ′ Ω ′ (a) Action of Ω × F u ′ − c δ ΩΩ v ′ Ω ′ u v (b) Action of − c δ Ω Figure 5:
Graphical representations of the two components of the infinitesimalrotation. (Ω , u , v ) denotes the position of the frame at time t while (Ω (cid:48) , u (cid:48) , v (cid:48) ) is itsposition at time t + dt with dt (cid:28)
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. Themotion of the vectors is indicated by a segment of circle in black color. (a) Actionof Ω × F : the vectors F and Ω × F are in plain and light black respectively. Thevector F is shown with unit norm for the ease of the representation but could be ofany norm in reality. The passage from (Ω , u , v ) to (Ω (cid:48) , u (cid:48) , v (cid:48) ) is via an infinitesimalrotation of axis Ω × F . (b) Action of δ : the vector − c δ Ω is shown in black. Thevectors Ω and Ω (cid:48) are identical and colinear to − c δ Ω. The passage from (Ω , u , v ) to(Ω (cid:48) , u (cid:48) , v (cid:48) ) is via an infinitesimal rotation of axis Ω. To better understand how the SOHB model (11) relates to other models, we re-writethe equation for Ω as follows: D t Ω = P Ω ⊥ F , (23)where P Ω ⊥ is the 3 × Ω ⊥ = I − Ω ⊗ Ω with ⊗ standing for the tensor (or outer)product. Eq. (23) bears similarities and differences with the momentum equation f isothermal compressible fluids. The latter is exactly recovered if the followingthree modifications are made:1. the projection matrix P Ω ⊥ is removed from (23) (i.e. it is replaced by I );2. c = c in the convective derivative D t (see (21));3. c = 0 in the expression of F (see (22)).Indeed, under these three modifications, we get the following system for ( ρ, U )where U = c Ω is the fluid velocity: ∂ t ρ + ∇ x · ( ρ U ) = 0 , ( ∂ t + U · ∇ x ) U = − Θ ∇ x log ρ. This is exactly the isothermal compressible Euler equations with the fluid temper-ature Θ = c c .We now investigate what consequences follow from undoing the above threemodifications, one by one.1. Introducing the projection P Ω ⊥ in (23) guarantees that the constraint | Ω | = 1is preserved in the course of time, if it is satisfied at time 0. Indeed, dotting Eq.(23) with Ω (and assuming that all functions are smooth) leads to D t | Ω | = 0,which guarantees that | Ω | is constant along the integral curves of the vectorfield c Ω. Thus, if | Ω | = 1 at time t = 0, it will stay so at any time.2. Having c (cid:54) = c is a signature of a loss of Galilean invariance. This is consistentwith 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 .Of course, this speed does not remain equal to c in frames that translate atconstant speed with respect to this frame.So far, with the introduction of P Ω ⊥ and different constants c (cid:54) = c but stillwith c = 0, the system for ( ρ, Ω) is decoupled from the equations for u and v and is written (see Eqs. (20a), (23) with F given by (22) in which c = 0): ∂ t ρ + c ∇ x · ( ρ Ω) = 0 , (24a) D t Ω = − c P Ω ⊥ ∇ x log ρ. (24b)This is nothing but the hydrodynamic limit of the Vicsek particle model(known as “Self-Organized Hydrodynamics (SOH)”) as established in [41, 44].This system has been shown to be hyperbolic [41] and to have local-in-timesmooth solutions [40].3. When c (cid:54) = 0, in addition to the pressure gradient, a second component ofthe force F appears. This component depends on the full rotation matrix A through Ω, u , v and their gradients (see Eq. 12). It is thus truly specific ofthe body orientation model.We are now going to compare the IBM and the SOHB models on a set of explicitstationary solutions of the SOHB model described in the next section. Special solutions of the macroscopic model
In this section, we exhibit three different classes of global-in-time solutions of theSOHB model (20). These three classes are characterized by uniform (i.e. indepen-dent of the spatial coordinate) fields ρ , r and δ . From now on we fix a wave-number(inverse of the length) ξ ∈ R \ { } and define ω = ξ c , λ = c + c . (25)We denote by x = ( x, y, z ) T the coordinates of x in the basis ( e , e , e ). The flocking state (FS) is a trivial but important special solution of the SOHBmodel (20) where both the density and rotation fields are constant (i.e. independentof time) and uniform: ρ ( t, x ) ≡ ρ = constant , A ( t, x ) ≡ A = constant , ∀ ( t, x ) ∈ [0 , ∞ ) × R . We have the following
Lemma 3.1.
The pair ( ρ, A ) consisting of a constant and uniform density ρ ( t, x ) = ρ = constant and the following rotation field: A ( t, x ) = ˜ A 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 ) (26)= A ( − ωt, e ) A ( ξz, e ) , (27) is a solution of the SOHB system (20) , where ω and ξ are given by (25) . We recallthat A ( θ, n ) is the rotation of axis n ∈ S and angle θ ∈ R defined by (8) . Thissolution will be referred to as a milling orbit (MO). The proof of this lemma is deferred to Appendix E.The MO is independent of x and y . Its initial condition is A mill (0 , z ) = A ( ξz, e ) = ξz ) − sin( ξz )0 sin( ξz ) cos( ξz ) . (28)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 . As z varies, the body-orientation rotates uniformly about the x -direction with spatial angular frequency ξ . As the rotation vector is perpendicular to the direction of variation, (28) iscalled a “perpendicular twist”. As time evolves, the rotation field is obtained by ultiplying on the left the initial perpendicular twist by the rotation A ( − ωt, e ).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 isagain independent of z . It belongs to the plane orthogonal to z and undergoesa uniform rotation about the z -axis. Consequently, the fluid streamlines, whichare the integral curves of c Ω, are circles contained in planes orthogonal to z ofradius c ω = c c ξ traversed in the negative direction if ξ >
0. These closed circularstreamlines motivate the “milling” terminology. It can be checked that the MOsatisfies: r = ξ (sin( ωt ) , cos( ωt ) , T , δ = 0 . As announced, r and δ are uniform but r depends on time. Actually, Ω × r = ξ e isindependent of time. The MO is depicted in Fig. 6 and its dynamics is visualizedin Video 2 (see Appendix A). (a) t = 0 (b) t > 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 certainnumber of points of the ( O, x, y ) and (
O, y, z ) planes. In (b), the rotation motion ofthe frame vectors is depicted by dotted circles of the color of the corresponding framevector. The red dotted circle can be seen as a depiction of the fluid streamlines. Seealso Appendix A, Video 2.Many examples of milling (also known as vortex) solutions have been observedin the collective dynamics literature as well as in biological systems [16, 25, 89]. Onthe modelling side, milling states have not been observed so far in alignment modelswithout the inclusion of an additional process such as an attraction-repulsion forcebetween the agents [17], a bounded cone of vision [24] or an anticipation mechanism[52]. The body-orientation framework is, to the best of our knowledge, a newsituation in which milling can be observed just with alignment assumptions. Millingstates can also be found in physical systems. A typical and important example is themotion of a charged particle in a uniform magnetic field, resulting in the formationof so-called cyclotron orbits. Once again, in the body-orientation framework, an xternal field is not needed and self-induced cyclotron orbits emerge only from thevariations of the internal body-orientation. Here, the analog of the magnetic fieldwould be Ω × r and the cyclotron frequency would be ω . Note that ω is under thecontrol of coefficient c which depends on the noise intensity 1 /κ . We have the following
Lemma 3.2.
The pair ( ρ, A ) consisting of a constant and uniform density ρ ( t, x ) = ρ = constant and the following rotation field: A ( t, x ) = ˜ A htw ( t, x )= ξ ( x − λt )) − sin ( ξ ( x − λt ))0 sin ( ξ ( x − λt )) cos ( ξ ( x − λt )) (29)= A ( ξ ( x − λt ) , e ) , (30) is a solution of the SOHB system (20) where ξ and λ are defined by (25) . Thissolution will be referred to as a helical traveling wave (HW). The proof of this lemma is given in Appendix E.2.The HW is independent of y and z . Its initial condition is A htw (0 , x ) = A ( ξx, e ) = ξx ) − sin( ξx )0 sin( ξx ) cos( ξx ) . (31)Here the self-propulsion direction is still independent of x and equal to e . Also,the body orientation still rotates uniformly about e with spatial angular frequency ξ but when x is varied instead of z . This means that the body orientation is nowtwisted when varied along the propagation direction. So, this initial condition iscalled a “parallel twist”. In the HW, the self propulsion direction Ω remains constantin time and uniform in space. The initial twist is propagated in time in this directionat speed λ and gives rise to a traveling wave˜ A htw ( t, x ) = ˜ A htw (0 , x − λt ) . Note that the traveling wave speed λ depends on the noise intensity 1 /κ and isdifferent from the fluid speed c . So, the frame carried by a given fluid elementfollowed in its motion is not fixed but rotates in time. Since Ω does not change, thefluid streamlines are now straight lines parallel to e . So, as a fluid element moves,the ends of the frame vectors u and v follow a helical trajectory with axis e , hencethe 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 graph-ically in Fig. 7. Its dynamics is visualized in Video 3 (see Appendix A). The HWbelongs to a larger class of solutions described in Appendix E.2. (a) t = 0 (b) t > Figure 7:
Graphical representation of the helical traveling wave (HW) at (a): initialtime, and (b): time t >
0. See Fig. 6 for captions. See also Appendix A, Video 3.
Clearly, in the definitions of the MO and HW, the choice of reference frame isunimportant. So, in the whole space R , such solutions exist in association withany reference frame. In a square domain of side-length L with periodic boundaryconditions, periodicity imposes some constraints on the direction of the referenceframe. For simplicity, we will only consider the case where the reference frame hasparallel axes to the sides of the square and ξ is linked to L by an integrality condition L ξ = 2 π n , with n ∈ Z \ { } .The study of the stability of the MO and the HW is left for future work. Bycontrast, the FS is linearly stable as the SOHB system is hyperbolic [36]. However,there is no guarantee that the FS at the level of the IBM is stable. Indeed, thereare strong indications that the FS is not stable for the Vicsek model [23] for someparameter ranges and a similar trend is likely to occur here.We can now answer the question posed at the end of Section 2.1.3 namelywhether the inclusion of the full body orientation makes any change in the dy-namics of the particle positions and directions compared to the Vicsek model. Tothis end, we consider the corresponding macroscopic models, i.e. the SOH model(24) for the Vicsek model and the SOHB model (11) for the body-orientation dy-namics. If we initialize the SOH model with uniform initial density ρ and meandirection Ω, inspection of (24) shows that the solution remains constant in timeand 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 andHW) also have uniform initial density ρ and mean direction Ω. If the dynamics ofthe particle positions and directions in the body orientation model was the same asin the Vicsek model, these three classes of solutions should have a constant meandirection Ω. However, it is not the case for the MO, where Ω changes with time andis subject to a planar rotation. This means that gradients of body attitude do have non-trivial influence on the direction of motion of the particles and that the bodyorientation model does not reduce to a Vicsek model for the particle positions anddirections.There is another, more subtle, difference between the two models concerningthe dynamics of Ω. It does not concern the MO and HW but we discuss it herein relation with the previous paragraph. Indeed, Fig. 4 reveals that the velocities c and c for the SOHB model crossover at a certain value κ ∗ of the concentrationparameter. The coefficients c and c for the SOH model can be found in [44], Fig.A1(b) and appear to satisfy c > c for the whole range of values of κ , i.e. donot exhibit any crossover. In particular, at large noise, the propagation velocity c of Ω in the SOHB model is larger than the mass transport velocity c . Thismeans that information (which triggers adjustments in Ω) propagates downstreamthe fluid by contrast to the Vicsek case where it propagates upstream. While thereason for this difference is unclear at this stage, we expect that it may induce largequalitative differences in the behavior of the system in some cases. This point willbe investigated in future work.Numerical simulation of the SOHB will be subject to future work. Here, we willrestrict ourselves to the MO and HW for which we have analytical formulas. In thenext section, using these two special solutions, we verify that the SOHB model andthe IBM are close in an appropriate parameter range. In this section we use the MO and HW to demonstrate the quantitative agreementbetween the SOHB model (11) and the IBM (3), (5) in the scaling (10). In thesimulations below, we consider a periodic cube of side-length L and choose R = 0 . , ν = 40 , c = 1 , L = 1 , ξ = 2 π, (32)so that RL = c ν L = 0 . (cid:28)
1, ensuring that the scaling (10) is satisfied. Further-more, we see that the choice of ξ is such that the twists in the MO or HW haveexactly one period over the domain size. N → ∞ In this section, we numerically demonstrate that the solutions of the IBM convergeto those of the macroscopic model in the limit N → ∞ and investigate the behaviorof the IBM at moderately high values of N .We sample N particles according to the initial condition (28) of the MO andsimulate the IBM (3), (5). We recall that the average direction Ω( t ) of the exactMO (26) is spatially uniform at any time and undergoes a uniform rotation motionabout the z -axis. So, we will compare Ω( t ) with the average direction Ω( t ) of allthe particles of the IBM, where Ω( t ) = (Ω , Ω , Ω ) T is defined by:Ω = (cid:80) Nk =1 Ω k ( t ) | (cid:80) Nk =1 Ω k ( t ) | , (provided the denominator is not zero, and where we recall that Ω k ( t ) = A k ( t ) e ).To ease the comparison, we compute the azimuthal and polar angles of Ω respectively efined by: ¯ ϕ := arg(Ω + i Ω ) ∈ [0 , π ) , ¯ θ = arccos(Ω ) ∈ [0 , π ] , (33)where arg( x + iy ) stands for the argument of the complex number x + iy . We notethat the corresponding angles ϕ and θ of Ω( t ) are given by ϕ ( t ) = − ω t = − π c ( κ ) t, θ = π/ , (34)where we have used (25) and (32) to compute the value of ω .Fig. 8a shows the azimuthal angle ¯ ϕ as a function of time over 5 units of time, forincreasing particle numbers: N = 5 10 (green curve), N = 1 . (orange curve)and N = 1 . (blue curve). Note that for very small values of N , the macroscopicmodel loses its relevance: below a few thousand particles we only observe a noisybehavior, not shown in the figure. For the considered range of particle numbers, wenotice that the angle ¯ ϕ decreases linearly with time, which shows that the behaviorof the IBM is consistent with the exact solution (34). However, quantitatively,we see that | d ¯ ϕ/ d t | depends on the particle number and decreases with increasingparticle number. We investigate this behavior in more detail in Fig. 8b wherethe difference between the measured angular velocity | d ¯ ϕ/ d t | and the theoreticalprediction 2 πc ( κ ) is plotted as a function of N . Each data point (blue dot) is anaverage of 10 independent simulations. This figure confirms that, as N increases, | d ¯ ϕ/ d t | decreases and converges towards 2 πc ( κ ). The inset in Fig. 8b shows thesame data points in a log-log-scale with the associated regression line (orange solidline). We observe that the error between the measured and theoretical angularvelocities behaves like N − α with a measured exponent α (cid:39) . | d ¯ ϕ/ d t | is close to N − agree with previ-ously documented observations in spherical statistics. Indeed, it has been shown in[81, Theorem 3(e)] that the estimation of the concentration parameter of a (spher-ical) von Mises distribution obtained from a crude averaging procedure from N independent samples produces a biased estimator with a (nonnegative) bias of or-der N − (see also [71, Section 10.3]). In the present case, a similar reasoning can beapplied, which we now briefly develop. The key observation is that all the measuredquantities are functions of empirical averages of the form (4). Under the chaos as-sumption (see Appendix D), when N is large, the body-orientations of the particlesbehave as N independent samples with common law M A , where A solves the SOHBmodel (11) and M A is defined by (6). In [35, Theorem 4.1], it has been shown that c ( κ ) can actually be expressed as a function of a certain number p of averagedquantities c ( κ ) = F ( (cid:104) g (cid:105) M A , . . . , (cid:104) g p (cid:105) M A ) , where g i : SO ( R ) → M ( R ) and F : M ( R ) p → R are smooth functions. The IBMsimulation thus defines an estimator ˆ κ of the concentration parameter such that c (ˆ κ ) = F (ˆ g , . . . , ˆ g p ) , where ˆ g i is the average of g i obtained by replacing M A by the empirical measure ofthe N body-orientations of the particles. We can then measure the bias by taking he expectation of the Taylor expansion of the previous expression around the point( (cid:104) g (cid:105) M A , . . . , (cid:104) g p (cid:105) M A ) : c (ˆ κ ) = c ( κ ) + δ ˆg · ∇ F + ( δ ˆg ) T (Hess F ) δ ˆg + R, where δ ˆg = (ˆ g , . . . , ˆ g p ) T − ( (cid:104) g (cid:105) M A , . . . , (cid:104) g p (cid:105) M A ) T and R is a remainder. The gradient ∇ and Hessian Hess are defined within the Euclidean framework given by (1). Bythe chaos hypothesis E [ δ ˆg ] = 0 and by the central limit theorem, the term of ordertwo behaves as N − . Since SO ( R ) is compact, higher order moments of δ ˆg canbe controlled by a classical argument based on Hoeffding’s inequality [87, Lemma5.5 and Theorem 5.29]. This ensures that E [ R ] is O ( N − ). We therefore obtain abiased estimator: E [ c (ˆ κ )] = c ( κ ) + aN + O ( N − ) , where a ∈ R depends on the derivatives of the considered functions and on the vari-ance of the estimator (4) where the particles are replaced by independent identicallydistributed samples with law M A . The fact that a > N , the fluctuations around theaverage (biased) value can be monitored by computing the standard deviation of the10 independent simulations. Fig. 8c shows this standard deviation as a function of N in a log-log-scale (blue dots). Although fluctuations remain significant with only10 simulations per data point, by a standard linear regression (solid orange line) weobtain that the size of the standard deviation behaves as N − β with β (cid:39) .
54. whichis close to the value β = 1 / (a) (b) (c) Figure 8: (a) Time evolution of the angle ¯ ϕ for three values of N : N = 0 .
05 10 (green curve), N = 0 .
15 10 (orange curve) and N = 1 . (blue curve). (b)Difference between the measured angular velocity | d ¯ ϕ/ d t | and the theoretical value2 πc ( κ ). Each data point (blue dot) is an average of 10 independent simulationswith the error bar showing one standard deviation. Solid black horizontal line at 0for convenience. Inset: same data in log-log scale and regression line (solid orangeline). (c) 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 . ν = 40, c = 1, κ = 10. .3.2 Quantitative comparison between the models In order to quantitatively confirm the agreement between the IBM and the macro-scopic model, we fix a large number N = 1 . of particles and we run the IBM fordifferent values of the concentration parameter κ and for the two classes of specialsolutions, the MO and the HW. To compare the models, we compute the followingmacroscopic quantities: • For the MO: starting from a sampling of the initial condition (28), we measurethe angular velocity | d ¯ ϕ/ d t | in a similar way as in the previous section. Giventhe parameter choice (32), the theoretical value of | d ϕ/ d t | predicted by (26)is | ω | = 2 πc ( κ ) where the function c is given by (17). • For the HW, starting from a sampling of the initial condition (31), we measurethe wave speed. To this aim, using (2), we compute the mean body-orientation A of the agents in a slice of size 10 − along the x -axis (which is the globaldirection of motion) as a function of time. As predicted by (29) the coefficient A of the mean orientation is a periodic signal. The inverse of the period ofthis signal (obtained through a discrete Fourier transform) gives the travelingwave speed of the HW. The theoretical value predicted by (29) is given by λ = c ( κ ) + c ( κ ) where the function c is given by (15).The output of these simulations is shown in Figs. 9a for the MO and 9b forthe HW. They respectively display the angular velocity and traveling wave speedobtained by running the IBM for a discrete set of values of κ (big blue dots). Bycomparison, 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 10independent simulations is 10 − . The relative error between the average measuredvalue and its theoretical prediction varies between 2% and 5% on the whole rangeof concentration parameters considered.These figures show an excellent agreement between the prediction of the macro-scopic SOHB model and the results obtained by running the IBM when the numberof particles is large. This confirms that the SOHB model provides an excellent ap-proximation of the IBM, at least during a certain period of time which is a functionof the particle number. We will see below that fluctuations induced by the finitenumber of particles may eventually destabilize the MO and lead to a HW or a FS. Asthese solutions are associated with different topological structure, these transitionswill be analyzed as topological phase transitions in the forthcoming sections. Both the MO and HW have non-trivial topology: inspecting the perpendicular twist(28) (see also Fig. 6a), we observe that the two-dimensional curve generated by theend of the vector u in the ( y, z )-plane as one moves along the z -axis is a closedcircle. A similar observation can be made on the parallel twist (31) (see Fig. 7a) asone moves along the x -axis. Both curves have therefore non-zero winding numbersabout the origin. When the domain is R , these winding numbers are ±∞ (wherethe 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 ξ a) (b) Figure 9: (a) MO: angular velocity | d ϕ/ d t | as a function of 1 /κ . (b) HW: travelingwave speed λ as a function of 1 /κ . Measured values from the IBM at discrete valuesof κ (big blue dots) and theoretical prediction from the SOHB model (dotted blackcurve). Parameters: N = 1 . , L = 1, ξ = 2 π , R = 0 . ν = 40, c = 1.by L = n π/ξ with n ∈ Z \ { } , then the winding numbers are equal to n . Asobserved on Formulas (26) and (29) (or on Figs 6b and 7b), this initial non-trivialtopological structure is propagated in time.When we initialize particles by sampling the initial conditions (28) or (31), weexpect that the solution of the IBM remains an approximation of the MO (26) orHW (29) respectively as evidenced in Section 3.3.2. However, noise induced byboth the inherent stochasticity of the IBM and finite particle number effects asexplained 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 sometime. This transition implies a change of the topology of the solution which, frominitially 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 statehas non-trivial topology and exhibits some kind of “topological protection” againstnoise-induced perturbations. To test this hypothesis quantitatively, we first need todevelop appropriate indicators. This is done in the next section. We will use two types of indicators. The first one is the global order parameter whichwill 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 rollangle. They will enable a finer characterization of topological phase transitions. .1 Global order parameter We first introduce the following scalar binary order parameter which measures thedegree of alignment between two agents with body-orientations A , ˜ A ∈ SO ( R ) : ψ ( A, ˜ A ) := 12 A · ˜ A + 14 . (35)In the quaternion framework (see Section 2.1.2 and Appendix B for details), wehave ψ ( A, ˜ A ) = ( q · ˜ q ) , (36)where q and ˜ q are two unit quaternions respectively associated to A and ˜ A , and q · ˜ q indicates the inner product of two quaternions. This expression makes it clear that ψ ( A, ˜ A ) ∈ [0 , ψ ( A, ˜ A ) measures thenematic alignment of the two associated unit quaternions, as it should because twoopposite quaternions represent the same rotation. We note that ψ ( A, ˜ A ) = 1 if andonly if ˜ A = A . On the other hand, ψ ( A, ˜ A ) = 0 if and only if A · ˜ A = − /
2, whichcorresponds to the two rotation axes being orthogonal and one rotation being aninversion about its axis.The Global Order Parameter (GOP) of a system of N agents at time t > N ( t ) = 1 N ( N − (cid:88) k (cid:54) = (cid:96) ψ (cid:0) A k ( t ) , A (cid:96) ( t ) (cid:1) . (37)From (37) we have GOP N ( t ) ∈ [0 , N indicates large disorder and alarge one, strong alignment. This is a global measure of alignment, by contrast to alocal one where ψ would be averaged over its neighbors only (and the result, averagedover all the particles). This global measure of alignment allows us to separate theMO and HW from the FS as shown below, which would not be possible with a localone.The GOP (37) can also be defined at the continuum level. As shown in Ap-pendix D, in the macroscopic limit, the particles become independent and identicallydistributed over R × SO ( R ), with common distribution ρ M A where ( ρ, A ) satisfiesthe SOHB system (11) and M A is the von Mises distribution (6). Therefore, theGOP of a solution of the SOHB system ( ρ, A ) is obtained as (37) where the sum isreplaced by an integral, A k ( t ) is replaced by A distributed according to the measure( ρ M A )( t, x , A ) d x d A and A (cid:96) ( t ) is replaced by ˜ A distributed according to the samemeasure, but independently to A . Therefore,GOP( ρ, A ) := (cid:90) (cid:90) ( R × SO ( R )) ψ ( A, ˜ A ) ρ ( x ) ρ (˜ x ) M A ( x ) ( A ) M A (˜ x ) ( ˜ A ) d x d˜ x d A d ˜ A. Using (7) and (8) one can prove that for any A ∈ SO ( R ), we have (cid:90) SO ( R ) A M A ( A ) d A = c ( κ ) c A , (38)with c ( κ ) defined by (14) and c being the particle speed. Using (35), we obtain:GOP( ρ, A ) = 12 (cid:18) c ( κ ) c (cid:19) (cid:90) R × R A ( x ) · A (˜ x ) ρ ( x ) ρ (˜ x ) d x d˜ x + 14 . (39) rom 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 solutionsdefined in Section 3.1. For the MO (26) and HW (29), for all time t > = 14 (cid:18) c ( κ ) c (cid:19) + 14 . (40)For the FS, A ( x ) ≡ A = constant and the GOP is equal toGOP = 34 (cid:18) c ( κ ) c (cid:19) + 14 . (41)Note that the GOP: GOP = 14 , corresponds to a disordered state of the IBM where the body-orientations of theparticles are chosen independently and randomly uniformly (or equivalently to theSOHB case κ → κ = 10 used in oursimulations, one can compute that:GOP (cid:39) . , GOP (cid:39) . . (42)The GOP values between GOP and GOP can be reached by generalized HW asshown in Appendix E.3. Let A = [Ω , u , v ] ∈ SO ( R ) be a body-orientation. Let θ ∈ [0 , π ], ϕ ∈ [0 , π ) be thespherical coordinates of Ω defined by (33) (omitting the bars). We let { Ω , e θ , e ϕ } be the local orthonormal frame associated with the spherical coordinates ( θ, ϕ ) andwe 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 θ θ . Since u and v belong to the plane spanned by p (Ω) and q (Ω), we let ζ ∈ [0 , π )be the angle between p (Ω) and u . Then, it is an easy matter to show that A = R (Ω) A ( ζ, e ). In aircraft navigation, θ , ϕ and ζ are respectively called the pitch,yaw and roll angles: the pitch and yaw control the aircraft direction with respectto the vertical and in the horizontal plane respectively, while the roll controls theplane attitude (see Fig. 10a). These angles are related to the Euler angles. Theconstruction of the roll angle ζ is summarized in Figure 10b. Pursuing the analogywith aircraft navigation, we see from Fig. 5 that F controls variations of pitch andyaw while δ controls variations of roll.As an example, we examine the pitch, yaw and roll of the three solutions of theSOHB model (11) described in Section 3.1. a) (b) Figure 10: (a) Pitch, yaw and roll angles of an aircraft with body orientation[Ω , u , v ] (original picture released under the Creative Commons CC0 licence by https://pixabay.com ). (b) Construction of the roll angle of A = [Ω , u , v ], wherethe vectors Ω, u and v are respectively in red, green and blue. The local frameis (Ω , p (Ω) , q (Ω)) where p (Ω) and q (Ω)) and the plane generated by them are inpurple. u and v belong to this plane. ζ is the angle between p (Ω) and u .1. FS: A is constant and uniform. Then, the pitch, yaw and roll are also constantand uniform.2. MO: A is given by (26) (see Figs. 6). Using Eq. (27), we have R (Ω) = A ( − ω t, e ) and the roll is given by ζ = ξz . The pitch and yaw are constantand uniform. The roll is constant in time and is also uniform on planes ofconstant z . The non-trivial topology of the MO results from the roll makinga complete turn when z increases by the quantity 2 π/ξ .3. HW: A is given by (29) (see Fig. 7). Then, we have R (Ω) = I and ζ = ξ ( x − λ t ). The pitch and yaw are constant and uniform while the roll is uniformon 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 completeturn when x increases by the quantity 2 π/ξ .The goal of the next section is to see how we can recover the roll field from thesimulation of a large particle system. As shown in the last section, the roll of the MO is uniform on planes of constant z . When simulating the MO by the IBM, we will use this property to compute anaverage 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 thedistance ∆ z between the planes is chosen appropriately, we can access to both theaverage and the variance of the roll. They will be collected into one single vector, the oll Polarization in planes of constant z or RPZ. A similar quantity characterizesthe HW, the Roll Polarization in planes of constant x or RPX. Below, we detail theconstruction 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 \ { } . The domain D is partitionedinto M slices of fixed size across z , where M is a fixed integer. For m ∈ { , . . . , M } ,the slice S m is defined by: S m := [0 , L x ] × [0 , L y ] × (cid:20) m − M L z , mM L z (cid:21) . Let us consider a system of N agents with positions and body-orientations ( X k , A k ),indexed by k ∈ { , . . . , N } . Each body orientation A k has roll ζ k ∈ [0 , π ). Wedefine the discrete RPZ for Slice m , ¯u m , by ¯u m := 1 N m (cid:88) k ∈ I m (cos ζ k , sin ζ k ) T ∈ R , (43)where I m = { k ∈ { , . . . , N } , X k ∈ S m } and N m is the cardinal of I m . Note thatthe RPZ ¯u m has norm smaller than one. The unit vector ¯u m / | ¯u m | or equivalently,its angle with the vector (1 , T gives the average roll in S m . The euclidean norm | ¯u m | is a measure of the variance of the set of roll angles { ζ k } k ∈ I m . If this varianceis small, then | ¯u m | ∼
1, while if the variance is large, | ¯u m | (cid:28)
1. When plotted inthe plane R , the set of RPZ { ¯u m } m =1 ,...,M forms a discrete curve referred to asthe RPZ-curve. It will be used to characterize the topological state of the particlesystem. A summary of this procedure is shown in Figure 11. The RPZ-curve is shown in Figure 12 (a) to (c), in the three following cases.1.
Disordered state: the particles are drawn independently uniformly randomlyin the product space D× SO ( R ). For each m , the RPZ (43) is an average ofuniformly distributed vectors on the circle and its norm is therefore close to0. 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 andtheir body-orientations independently according to a von Mises distribution M A with a fixed mean body orientation A ∈ SO ( R ). In this case, forall slices, the corresponding RPZ (43) is an average of identically distributedvectors on the circle whose distribution is peaked around the same point of theunit circle, and the peak is narrower as κ is larger. Therefore, the RPZ vectors(43) concentrate on a point near the unit circle (Figure 12b). The RPZ-curvereduces 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 independentlyaccording to a von Mises distribution M A mill(0 ,z ) with A mill (0 , z ) defined by (28)(with ξ = 2 π/L z ). This time, the von Mises distribution is peaked around apoint which depends on z . For each slice, the position of the RPZ (43) depends Construction of the RPZ and graphical representation. The spatialdomain D is partitioned into M slices represented in different colors (top left). Ineach slice S m , we have I m particles with roll ζ k each of them plotted in the particle’slocal plane spanned by p (Ω k ), q (Ω k ) (top right: we plot 3 particles in the slice S ).Note that the local planes of different particles of the same slice may not coincidewhen imbedded in R . For this given slice, the RPZ ¯u m is computed and plotted in R (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 thenplotted on a single figure in the same color as the slice it corresponds to (bottomleft). This collection of points forms a discrete curve (here a fragment of a circle):the RPZ-curve.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 extractedfrom the RPZ-curve:1. the distance of its center of mass to the origin d z : d z = (cid:12)(cid:12)(cid:12) M M (cid:88) m =1 ¯u m (cid:12)(cid:12)(cid:12) , (44)2. its mean distance to the origin ¯ r z :¯ r z = 1 M M (cid:88) m =1 | ¯u m | , (45)3. its winding number about the origin w z : for m ∈ { , . . . , M } , let β m =arg (cid:0) ( ¯u m ) + i ( ¯u m ) (cid:1) ∈ [0 , π ) (with ¯u m = (( ¯u m ) , ( ¯u m ) ) T ) and δβ m +1 / ∈ a) (b)(c) (d) Figure 12:
Examples of RPZ-curves: in each figure, the roll Polarization RPZvectors corresponding to M = 1000 slices are plotted. The colorbar to the right ofeach figure assigns a unique color to each slice. The same color is used to plot thecorresponding RPZ. In each figure the unit circle and its center are represented inblue. (a) Disordered state: all RPZ concentrate near the origin ; (b) FS: all RPZconcentrate on a point close to the unit circle ; (c) MO (28): the RPZ-curve is adiscrete circle centered at the origin and of radius close to unity. The total numberof particles is N = 1 . · . Note that in Figs. (a) and (b), all RPZ are superimposedand only the last one (in magenta color) is visible. (d) Quantifiers of RPZ curvemorphology: point G (in red) is the center-of-mass of the RPZ curve and d z is itsdistance to the origin O (shown in blue). The mean radius ¯ r z of the RPZ curve isillustrated by the circle in black broken line which has same radius. The windingnumber, which is the number of turns one makes following the spectrum of colorsin the same order as in the colorbar from bottom to top (the green arrow indicatesthe direction of progression along the RPZ curve) is w z = − − π, π ) be such that δβ m +1 / ≡ β m +1 − β m modulo 2 π , where we let β M +1 = β . Then: w z = 12 π M (cid:88) m =1 δβ m +1 / , see e.g. [61, p. 176]).The subscript z indicates that the slicing has been made across z . Similar quantitieswith an index ’ x ’ will correspond to the slicing made across x . Fig. 12d providesa graphical illustration of the triple ( d z , ¯ r z , w z ). For the examples given above, thistriple has the following values:Disordered state: ( d z , ¯ r z , w z ) = (0 , , ND) , where ND stands for “undefined” , (46)FS: ( d z , ¯ r z , w z ) ≈ (1 , , , (47)MO: ( d z , ¯ r z , w z ) ≈ (0 , , w ) , with w (cid:54) = 0 . (48)We have a similar conclusion with ( d x , ¯ r x , w x ) for a disordered state or an FS. Foran HW, we have ( d x , ¯ r x , w x ) ≈ (0 , , w ) with w (cid:54) = 0. Thus, monitoring either orboth triples (according to the situation) will give us an indication of the state ofthe system in the course of time. In particular, non-trivial topological states areassociated with non-zero winding numbers w x or w z . In practice, we will use thenonzero-rule algorithm to compute the winding numbers numerically [61, p. 176]. 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 starkcontrast with the macroscopic model for which they persist for ever. The transitionof a MO or HW to a FS implies a topological change. To analyze whether theMO or HW are more robust due to their non-trivial topological structure (i.e. aretopologically protected), we will compare them with similar but topologically trivialinitial conditions. Then, we will investigate rarer events, where an MO does nottransition directly to an FS.
In this section, we compare the solutions of the IBM with two initial conditionsusing the perpendicular or parallel twists as building blocks. But one will have anon-trivial topology and the other, a trivial one. Specifically we define the followinginitial condition:
Milling orbit
Let D = [0 , L ] × [0 , L ] × [0 , L ] be a rectangular domain withperiodic boundary conditions and let ξ = 2 π/L . We consider the following twoinitial conditions: • Double mill initial condition MO1: A m, (0 , z ) = A ( ξ z, e ) , z ∈ [0 , L ] , (49)where we recall again that A ( θ, n ) is the rotation of axis n ∈ S and angle θ ∈ R defined by (8). This initial condition has non-trivial topology: the curvegenerated by the end of the vector u in the ( y, z )-plane as z ranges in [0 , L ] akes two complete turns around the origin in the same direction. Thus, thisinitial condition has winding number equal to 2. • Opposite mills initial condition MO2: A m, (0 , z ) = (cid:26) A ( ξ z, e ) , z ∈ [0 , L ] , A ( − ξ z, e ) , z ∈ [ L, L ] . (50)This initial condition has trivial topology: starting from z = 0, the curvegenerated by the end of the vector u makes one complete turn around the originin the counterclockwise direction until it reaches z = L but then reverses itsdirection and makes a complete turn in the clockwise direction until it reaches z = 2 L . Thus, this initial condition has winding number equal to 0 and hastrivial topology. Helical traveling wave
Let now D = [0 , L ] × [0 , L ] × [0 , L ]. Compared to theprevious case, the domain has size 2 L 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: A h, (0 , x ) = A ( ξ x, e ) , x ∈ [0 , L ] , (51)This initial condition has non-trivial topology and has winding number equalto 2 by the same consideration as for initial condition MO1. • Opposite helices initial condition HW2: A h, (0 , x ) = (cid:26) A ( ξ x, e ) , x ∈ [0 , L ] , A ( − ξ x, e ) , x ∈ [ L, L ] . (52)Again, by the same considerations as for MO2, this initial condition has trivialtopology, i.e. winding number equal to 0.We initialize the IBM by drawing N positions independently uniformly randomlyin the spatial domain and N body-orientations independently from the von Misesdistribution M A (0 , x ) where A (0 , x ) is one of the initial conditions MO1 or MO2.Then, we run the IBM and record various indicators (discussed below) as functionsof time. The results are plotted in Fig. 13, as plain blue lines for the solution issuedfrom MO1 (the topologically non-trivial initial condition), and as broken orangelines for that issued from MO2 (the topologically trivial one). We proceed similarlyfor the two initial conditions HW1 and HW2 and display the results in Fig. 14.See also Videos 4 to 7 in Appendix A supplementing Fig. 13 and Videos 8 to 11supplementing Fig. 14. Figs. 13a and 14a display the GOP. We observe that, for all initial conditions,the GOP has initial value GOP , which is consistent with the fact that the initialconditions are either MO or HW. Then, again, for all initial conditions, at largetimes, the GOP has final value GOP which indicates that the final state is a FS. his is confirmed by the inspection of the second line of figures in Figs. 13 and14 which provide the triplet of topological indicators ( d z , ¯ r z , w z ) for MO solutionsand ( d x , ¯ r x , w x ) for HW solutions. Specifically, d z and d x are given in Figs. 13dand 14d respectively, ¯ r z and ¯ r x in Figs. 13e and 14e, and w z and w x in Figs. 13fand 14f. Initially both triplets have value (0 , ,
2) as they should for MO1 or HW1solutions (see (48)). Their final value is (1 , ,
0) which indicates a FS (see (47)). Thefact that the final state is a FS implies, for MO1 and HW1, first that the IBM hasdeparted from the MO and HW exact solutions of the macroscopic model describedin Sections 3.1.2 and 3.1.3, and second, that a topological phase transition has takenplace, bringing the topologically non-trivial MO1 and HW1 to a topologically trivialFS. For the topologically trivial MO2 and HW2 initial conditions, no topologicalphase transition is needed to reach the FS. The differences in the initial topology ofthe solutions induce strong differences in the trajectories followed by the system.For the topologically non-trivial initial conditions MO1 or HW1, the systemremains in the MO or HW state for some time; hence it follows the macroscopicsolution during this phase. Indeed, the GOP displays an initial plateau at the valueGOP , while the triplet of topological indicators stays at the value (0 , , ϕ (Fig. 13c, blue curve), which varies linearly in time and by the pitch ¯ θ (Fig. 13bblue curve) which is constant in time, consistently with the MO solution of themacroscopic model (Section 3.1.2) (see also Fig. 8a for the linear variation of theyaw). The duration of this initial phase, also referred to as the persistence time, issignificantly longer for HW1 than for MO1. In our experiments, the former can reachseveral hundred units of time and sometimes be infinite (up to our computationalcapacity). 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 timehighlighted in Figs. 13, 14 and subsequent figures by a grey shaded zone, showingthat 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. Thetopological transition precisely occurs at this time with a transition of the windingnumber from 2 to 0 through a short sequence of oscillations. However, ¯ r has notreached 0 and d has already started to increase, which suggests that disorder is notcomplete. At this time also, the linear variation of ¯ ϕ suddenly stops and ¯ ϕ remainsconstant afterward, while ¯ θ shows a small oscillation and jump. For HW1, ¯ θ and ¯ ϕ are initially plateauing with small oscillations. At the time when the system leavesthe HW state (around t (cid:39) ϕ from 2 π to π whichindicates that the system suddenly reverses its average direction of motion. TheGOP starts to decrease significantly before this time so we can infer that duringthe time period between t (cid:39)
125 and t (cid:39) 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 topologicallytrivial initial conditions MO2 or HW2 as shown by the GOP immediately leaving he value GOP . For HW2 the GOP increases right after initialization and smoothlyreaches the value GOP , at a much earlier time than HW1. The trend is differentfor MO2. In this case, the GOP first decreases. Then, after a minimum value, itincreases again and smoothly reaches the value GOP at a time similar to MO1.The initial decay of the GOP for the MO2 solution can be explained by the fact thatthe macroscopic direction Ω turns in opposite directions for the two opposite mills,thus decreasing the global order. For HW2, the macroscopic direction stays constantand uniform. So, it is the same for the two opposite helices, giving rise to a largerGOP. The mean radii ¯ r z and ¯ r x stay constant it time, showing that the evolutions ofMO2 and HW2 do not involve phases of larger disorder. The quantity d x increasesmonotonically towards the value 1 while d z is subject to some oscillations close toconvergence. This is due to the fact that the RPZ or RPX curves stay arcs ofcircles with decreasing arclength for the RPX and with some arclength oscillationsfor the RPZ as displayed in Videos 7 and 11. Of course, the winding number staysconstant equal to 0 as it should for topologically trivial solutions. In both the MO2and HW2 cases, ¯ θ and ¯ ϕ remain constant throughout the entire simulation. In theMO2 case, this is the consequence of the two counter-rotating mills which preservethe direction of motion on average. In the HW2 case, this is due to the fact thatthere is no variation of the direction of motion for HW solutions in general (see alsoVideo 6 and Video 10). Again, we observe that the convergence towards the FStakes more time for HW2 than for MO2. This points towards a greater stability ofthe HW-type solutions compared to the MO ones.The existence of a persistence time for the MO1 and HW1 solutions suggeststhat they enjoy some kind of topological protection against the noisy perturbationsinduced by the IBM and that MO2 and HW2 do not have such protection. However,explicit solutions of the SOHB model for the initial conditions MO2 and HW2 arenot available. The evolution of the IBM for these initial conditions could just be ashadow of the corresponding solutions of the SOHB model only weakly perturbedby noise. Even in a scenario where noise would have perceivable influence, it isnot clear if the induced disruption would arrive sooner or later than for MO1 orHW1. So, while these simulations suggest that the body orientation model exhibitstopologically protected solutions, further investigations are needed to confirm thatnon-trivial topology actually provides increased robustness against perturbations. Although the scenario described in the previous section is the most common one,the IBM sometimes leads to different, slightly more complex scenarios which aredescribed in the present section. now, the IBM is initialized by drawing N positionsindependently uniformly in the cubic domain D = [0 , L ] × [0 , L ] × [0 , L ] with periodicboundary conditions and N body-orientations independently from the von Misesdistribution M A (0 , x ) where A (0 , x ) is given by (28) with ξ = 2 π/L (winding numberequal to 1). Here, we report on the occurrence of transitions from a MO to a HW. We run theIBM and record the time-evolution of a set of indicators as shown in Fig. 15 (see a) (b) (c)(d) (e) (f) Figure 13:
Examples of solutions of the IBM for initial conditions sampled fromthe double mill MO1 (plain blue curves) and the opposite mills MO2 (boken orangecurves). The following indicators are plotted as functions of time: (a) Global OrderParameter (GOP) (see Eq. (37)). Horizontal lines at GOP values 0 .
25, 0 .
45 and 0 . , GOP and GOP respectively corresponding tototally disordered states, MO or HW, and FS (see Eqs. (40)-(42)). (b) Pitch angle¯ θ of the global particle average direction ¯Ω (see (33)). (c) Yaw ¯ ϕ of ¯Ω. (d) Distanceof center of mass of RPZ curve to the origin d z (see (44)). (e) Mean distance ofRPZ curve to the origin ¯ r z (see (45)). (f) Winding number of RPZ curve w z (see(44)). Grey shaded zones highlight a small region around the time of minimal GOPfor the MO1 solution. Parameters: N = 3 10 , R = 0 . κ = 10, ν = 40, c = 1, L = 1, ξ = 2 π . See also Videos 4 to 7 in Appendix A.also supplementing videos 12 to 14 in Appendix A).As shown in Fig. 15a, the GOP does not converge towards GOP characterizingthe FS, but towards an intermediate value between GOP (which characterizes MOor HW) and GOP . As explained in Appendix E.3, such values of the GOP canbe attained by a generalized helical wave solution (as can be observed in Video12). The pitch ¯ θ (Fig. 15b) and yaw ¯ ϕ (Fig. 15c) behave like in the milling-to-flocking transition (see Figs. 13b and 13c) except for small-amplitude, slow-frequency oscillations appearing after the topological transition time. This maybe due to some competition between two attractors, the FS and the HW, whichbeing alternately stronger and weaker, generate this oscillatory behavior. Note thata transition to a HW cannot occur when the global direction of motion at thetransition time is not one of the principal axes of the square domain since a HW a) (b) (c)(d) (e) (f) Figure 14:
Examples of solutions of the IBM for initial conditions sampled fromthe double helix HW1 (plain blue curves) and the opposite helices HW2 (brokenorange curves). The following indicators are plotted as functions of time: (a) GlobalOrder Parameter (GOP). (b) Pitch angle ¯ θ of ¯Ω. (c) Yaw ¯ ϕ of ¯Ω. (d) Distance ofcenter of mass of RPX curve to the origin d x . (e) Mean distance of RPX curve tothe origin ¯ r x . (f) Winding number of RPX curve w x . Grey shaded zones highlight asmall region around the time of minimal GOP for the HW1 solution. The HW2 andHW1 solutions are computed during 200 and 250 units of time respectively. The twosimulations have reached equilibrium by their final time. Parameters: N = 3 10 , R = 0 . κ = 10, ν = 40, c = 1, L = 1, ξ = 2 π . See caption of Fig. 13 for furtherindications. See also Videos 8 to 11 in Appendix A.along another direction is not compatible with the periodic boundary conditions(see Section 5.2.2). This is confirmed by the final values of ¯ ϕ and ¯ θ (both equalto π/
2) which correspond to a global direction of motion oriented along the y -axis(in what follows, in reference to (51) and to avoid confusion, we will still call thatdirection, the x direction).The second and third lines of figures in Fig. 15 show the triplets of topologicalindicators ( d z , ¯ r z , w z ) and ( d x , ¯ r x , w x ) which materialize the MO and HW structuresrespectively. The mean distance of the RPZ-curve to the origin ¯ r z (Figs. 15e)decreases, revealing an increase of the disorder. Simultaneously, the distance of itscenter of mass to the origin d z increases (Figs. 15d) showing a transition trendto a FS. The winding number w z (Fig. 15f) jumps from 1 to 0 at the time ofmaximal disorder. However, d z and ¯ r z do not reach zero, showing that completedisorder across z is not reached. Since the final state of the system is a generalized elical wave state (see Appendix E.3), we do not necessarily expect that completedisorder will be reached along the z -direction. In the mean time, ¯ r x starts from 0(complete disorder) and increases up to a value close to unity, showing the build-upof a HW. The quantity d x increases during some time but eventually decreases to0 (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 upto 1 at the time where the winding number w z drops to 0. There is a transfer ofnon-trivial topology from an MO structure to a HW structure. In some rare cases an intermediate unstable HW can be observed. The procedureis the same as in the previous section. Fig. 16 shows the results (see also supple-menting videos 15 and 16 in Appendix A).The transition stage between the MO and FS is significantly longer than in theprevious situations. During that phase, the GOP (Fig. 16a) oscillates betweenthe value Ψ characterizing the MO and lower values, i.e. lower order. Likewise,there are significant variations of the pitch ¯ θ (Fig. 16b) and yaw ¯ ϕ (Fig. 16c). Asin the previous section, this could be explained by antagonist effects of differentattractors (the MO and HW) and subsequent oscillations of the system betweenthem. Video 15 reveals large scale band structures similar to a HW except that theglobal 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 stagecould be explained in the present case by the fact that the global direction ofmotion seems to oscillate around the direction given by e + e (i.e. ϕ = π/ θ = π/
2) which is theoretically compatible with the periodic boundary conditions,provided the wave length ξ is changed from 2 π/L to √ π/L . This state does notseem to be stable as shown by the large oscillations of ¯ ϕ and ¯ θ . The topologicalindicators ( d z , ¯ r z , w z ) shown in the second line of figures of Fig. 16 also display largeoscillations. The quantity ¯ r z drops, and at the same time, d z remains small, whilethe winding number w z has strong oscillations, indicating a state of large disorderacross z , which is consistent with the fact that the temporary HW order is organizedin a different direction. However, we see that w z has a calmer period between twoseries of oscillations. This calmer period corresponds to the interval of time duringwhich the temporary HW order prevails. Eventually the triplet converges to thevalue (1 , ,
0) characterizing the FS.
An Individual Based Model describing the alignment of body-orientations in 3D andits macroscopic limit have been presented. The model involves new kinds of inter-nal degrees of freedom involving geometrical constraints, here due to the manifoldstructure of SO ( R ), leading to new types of self-organized phenomena. In particu-lar, the macroscopic model has been shown to host special solutions with non-trivialtopological structures. Corresponding solutions of the Individual Based Model havebeen computed and their non-trivial topological structure, shown to persist for a a) (b) (c)(d) (e) (f)(g) (h) (i) Figure 15:
Transition from a MO to a HW: example of a solution of the IBM foran initial condition sampled from (49) 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 ¯ r z . (f) Winding number of RPZcurve w z . (g) Distance of center of mass of RPX curve to the origin d x . (h) Meandistance of RPX curve to the origin ¯ r x . (i) Winding number of RPX curve w x . Greyshaded zones highlight a small region around the time of minimal GOP. Parameters: N = 1 . · , R = 0 . L = 1, D = 0 . ν = 40, c = 1. See caption of Fig. 13 forfurther indications. See also Videos 12 to 14 in Appendix A. a) (b) (c)(d) (e) (f) Figure 16:
Transition from a MO to a FS via an unstable HW: example of a solutionof the IBM for an initial condition sampled from (49) in the rare case where it leadsto a FS through a transient HW. The following indicators are plotted as functionsof time: (a) GOP (b) Pitch ¯ θ of ¯Ω. (c) Yaw ¯ ϕ of ¯Ω. (d) Distance of center of massof RPZ curve to the origin d z . (e) Mean distance of RPZ curve to the origin ¯ r z .(f) Winding number of RPZ curve w z . Grey shaded zones highlight a small regionaround the time of minimal GOP. Parameters: N = 1 . · , R = 0 . L = 1, D = 0 . ν = 40, c = 1. See caption of Fig. 13 for further indications. See alsoVideos 15 and 16 in Appendix A.certain time before being destroyed by noise-induced fluctuations. Quantitative es-timates of the agreement between the Individual Based Model and the Macroscopicmodel have been given. This study provides one more evidence of the role of ge-ometry and topology in the emergence of self-organized behavior in active particlesystems. The model presented in this article opens many new research directions.Some of them are listed below.1. The stability of the MO (26) and HW (29) solutions as well as those of thegeneralized HW solutions described in Appendix E is an open problem. Itwould enable us to investigate the potential link between topological structureand stability.2. Numerical simulations have been carried out in a periodic setting. Real sys-tems though are confined by solid walls. To model the influence of confinement, t is necessary to explore wider classes of boundary conditions.3. Most topological states in physical systems consist of linear perturbations ofbulk states that propagate on the edges of the system (edge states). It wouldbe interesting to determine whether linear perturbations of the MO or HWsolutions could host such edge states.4. Beyond the mean-field limit N → ∞ , it would be interesting to quantifythe fluctuation about the mean-field, for instance through a large deviationapproach (see e.g. [6, 7, 11, 30, 46, 67]).5. Direct numerical simulations of the macroscopic model need to be developedto answer some of the questions raised by the study of topological protection(see Section 5.1).6. It is desirable to develop more sophisticated topological indicators to gainbetter insight into the topological structure of the solutions.7. The multiscale approach developed here could be extended to other geometri-cally structured systems involving e.g. a wider class of manifolds which wouldenlarge the applicability of the models. ppendices A List of supplementary videos
This article is supplemented by several videos which can be accessed by followingthis link. They are listed and described below.
Video 1.
It supplements Fig. 3 of Section 2.1.2 and provides a visualization of thetime evolution of the system considered in this figure.
Video 2.
It supplements Fig. 6 of Section 3.1.2: it provides a visualization of thetime evolution of a MO. Several frames A = (Ω , u , v ) ∈ SO ( R ) are placed at variouslocations of space and evolve according to (26) (with arbitrary chosen parameters).The vectors Ω, u and v are displayed respectively in red, green and blue. Video 3.
It supplements Fig. 7 of Section 3.1.3: it provides a visualization ofthe time evolution of a HW. See caption of Video 2 for details on the graphicalrepresentation.
Video 4.
It supplements Fig. 13 in Section 5.1. It shows the time-evolution ofthe particles for the initial condition MO1 (49). For clarity, only a sample of 5000particles are shown. We refer to Fig. 3a for details on the representation of the bodyorientation using four-colored tetrahedra. We notice the ensemble rotation of theparticle directions about the z axis until an instability disrupts the body orientationtwist along the z axis (around time t ≈
13) and eventually drives the system to aFS.
Video 5.
It supplements Fig. 13 in Section 5.1. It provides the time-evolutionof the RPZ curve for the initial condition MO1 (49). The RPZ curve remains acircle until time t ≈ t ≈
13 which is the first time when the RPZ-curve passes through the origin; atthis time, the widing number is not defined. Then, the RPZ-curve slowly migratestowards the unit circle while shrinking to a single point which signals a FS. Fromtime t ≈
15 on, it remains a single immobile point.
Video 6.
It supplements Fig. 13 in Section 5.1. It shows the time-evolution ofthe particles for the initial condition MO2 (50). For clarity, only a sample of 5000particles are shown (see Fig. 3a for details on the representation of the body orien-tation). We notice the counter-rotation of the particle directions about the z axisin the bottom and top halves of the domain, corresponding to the opposite mills.These two counter-rotations gradually dissolve while the solution approaches theFS. Video 7.
It supplements Fig. 13 in Section 5.1. It provides the time-evolution ofthe RPZ curve for the initial condition MO2 (50). The circle formed by the initialRPZ curve immediately opens. The opening width constantly increases, until thearc is reduced to a single point opposite to the opening point at time t ≈
10. Thenthere is a bounce and the arc forms again and increases in size until it reaches maximum and decreases again. Several bounces are observed with decreasingamplitudes. These bounces result in the non-monotonous behavior of the quantity d z displayed on Fig. 13d. Video 8.
It supplements Fig. 14 in Section 5.1. It shows the time-evolution ofthe particles for the initial condition HW1 (51) (see Fig. 3a for details on therepresentation of the body orientation). For clarity, only a sample of 5000 particlesare shown. Before time t (cid:39) t ≈ t ≈ θ (cid:39) θ (cid:39) π , which suddenly reverses the global direction of motion. After time t ≈ t = 178,the particles still have an undulating behavior but it quickly fades away until a FSis reached. Video 9.
It supplements Fig. 14 in Section 5.1. It shows the time-evolution of theRPX-curve for the initial condition HW1. Unlike in the MO case, the RPX curvedoes 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 isconverging attracts the RPX. During this transition, the circular shape of the RPXcurve is preserved until it becomes a point.
Video 10.
It supplements Fig. 14 in Section 5.1. It shows the time-evolution ofthe particles for the initial condition HW2 (52). For clarity, only a sample of 5000particles are shown (see Fig. 3a for details on the representation of the body orien-tation). At the beginning, we see two opposite alternations of the three side colorsof the tetrahedra (green-blue-magenta followed by green-magenta-blue), which sig-nals a double parallel twist. Then, gradually, the green color is eaten up by theblue and magenta ones and only one alternation of the blue and magenta colorsremains. Then the color alternation shades away and gives room to a homogeneouscolor showing that the body orientations have stopped rolling and a FS is attained.
Video 11.
It supplements Fig. 14 in Section 5.1. It provides the time-evolution ofthe RPX curve for the initial condition HW2 (52). The circle formed by the initialRPX curve immediately opens. The opening width constantly increases, althoughat a slower pace than for MO2 (see Video 7). Here, also contrasting with theMO2 case, the monotonous opening of the arc results in a monotonously increasingquantity d x as shown in Fig. 14d. Video 12.
It supplements Fig. 15 in Section 5.2.1. It shows the time-evolution ofthe particles for a MO initial condition (49) in a rare case where it evolves into aHW. For clarity, only a sample of 5000 particles are shown (see Fig. 3a for details onthe representation of the body orientation). It starts like Video 4 with the ensemblerotation of the particle directions about the z axis until an instability initiated attime t ≈
10 gradually disrupts this organization. However, the disruption does notdrive the system to an FS, but rather to a HW as shown by the alternations of blue,green and magenta colours propagating along the particle orientations. ideo 13. It supplements Fig. 15 in Section 5.2.1. It provides the time-evolution ofthe RPZ curve for a MO initial condition (49) in a rare case where it evolves into aHW. The behavior is essentially the same as in Video 5 except that the RPZ-curveshrinks to a single point far away from the unit circle. This shows that the endstate 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 asfor the milling-to-flocking transition. Video 14.
It supplements Fig. 15 in Section 5.2.1. It provides the time-evolution ofthe RPX curve for a MO initial condition (49) in a rare case where it evolves into aHW. Initially, the RPX-curve is reduced to the origin, showing total disorder acrossthe x direction. Then, after some chaotic transient, a closed curve enclosing theorigin is formed. This curve initially stays close to the origin, still showing strongdisorder. But gradually, the radius of the curve increases and approaches the unitcircle. Thus, across x , the topology is initially undefined, but when it builds up, itshows its non-trivial character, the emerging RPX-curve having non-zero windingnumber about the origin. Video 15.
It supplements Fig. 16 in Section 5.2.2. It shows the time-evolution ofthe particles for a MO initial condition (49) in a rare case where it evolves into a FSthrough a transient HW. For clarity, only a sample of 5000 particles are shown (seeFig. 3a for details on the representation of the body orientation). The point of viewis changed from Video 12 to better visualize the transient HW moving along thediagonal, appearing around time t ≈
16. At the beginning we witness the ensemblerotation of the particles and its disruption by an instability. After some chaoticbehavior, the transient HW establishes as shown by the alternations of blue, greenand magenta colours propagating along the diagonal. But after some time, the HWstructure is disrupted again and the system eventually establishes a FS.
Video 16.
It supplements Fig. 16 in Section 5.2.2. It provides the time-evolutionof the RPZ curve for a MO initial condition (49) in a rare case where it evolvesinto a FS through a transient HW. The behavior is essentially the same as in Video5 except that the RPZ-curve undergoes a longer-lasting chaotic dynamics beforeshrinking to a point which migrates towards the unit circle.
B Quaternion framework
Despite its formal simplicity, the SO ( R )-framework used in the definition of theIndividual Based Model is not well suited to numerical simulations due to the highcomputational cost required to store and manipulate rotation matrices. A moreefficient representation of rotations in R is the quaternion representation based onthe group isomorphismΦ : H / ± −→ SO ( R ) q (cid:55)−→ Φ( q ) : w ∈ R (cid:55)→ { q [ w ] q ∗ } ∈ R , where the 3-dimensional vector w = ( w , w , w ) T ∈ R is identified with the pureimaginary quaternion denoted by [ w ] = iw + jw + kw and q ∗ denotes the conjugate uaternion to q . Conversely, the pure imaginary quaternion q = iq + jq + kq is identified with the 3-dimensional vector denoted by { q } := ( q , q , q ) T . Notethat for any quaternion q and any vector w ∈ R , the quaternion q [ w ] q ∗ is apure imaginary quaternion. The group of unit quaternions is denoted by H and ishomeomorphic to the sphere S ⊂ R .We refer the reader to [38, Section 2] and [37, Appendix A] where details aboutthe equivalence between the two representations can be found. Note that [37] studiesa model in a full quaternion framework. Table 2 below summarizes how the differentobjects can be computed in either of the two representations.Matrix QuaternionOrientation A ∈ SO ( R ) q ∈ H / ± q ) = A Flux J k = (cid:80) j K ( X k − X j ) A j Q k = (cid:80) j K ( X k − X j ) ( q j ⊗ q j − / )Mean orientation A = arg max { A (cid:55)→ A · J } ¯ q ∈ H eigenvector associated to thelargest eigenvalue of Q Von Mises distri-bution M A ( A ) = exp( κ A · A ) Z M q ( q ) = exp(2 κ ( q · q ) ) Z Table 2:
Matrix vs quaternion formulation C Numerical methods
The IBM (3),(5) has been discretized within the quaternion framework using thetime-discrete algorithm described in Table 3 below. This table shows one iterationof the algorithm during which the positions X nk ∈ R and orientations q nk ∈ H for k ∈ { , . . . , N } are updated into X n +1 k and q n +1 k respectively.At step 2, the Poisson process is discretized with a time step ∆ t during whichthe indices of the jumping agents are recorded. In the simulations ∆ t has to bechosen small enough so that the event that an agent jumps twice or more during atime interval of size ∆ t is negligible. In all the simulations, we take ∆ t such that ν ∆ t = 10 − .At step 3, a random quaternion q sampled from a von Mises distribution withprescribed mean orientation ¯ q can be obtained as q = ¯ qr where r ∈ H is sampledfrom a von Mises distribution with mean orientation 1 (see [38, Proposition 9]). Anefficient rejection algorithm to sample von Mises distributions can be found in [65].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 0and of radius R >
0. The six parameters of the simulations are summarized inTable 1. lgorithm: Iteration n → n + 1 of the time-discrete algorithm
1. Update the positions: for k ∈ { . . . , N } , set X n +1 k = X nk + c { q nk [ e ]( q nk ) ∗ } ∆ t
2. Draw a subset I ⊂ { , . . . , N } of jumping agents: for each agent k ∈ { . . . , N } ,draw a random number r k uniformly in [0 , r k > exp( − ν ∆ t ), then k ∈ I .
3. Compute the local flux: for k ∈ I , compute Q nk = N (cid:80) Nj =1 K ( X nk − X nj ) ( q nj ⊗ q nj − I ) .
4. Update the orientations: for k ∈ I compute one unit eigenvector q nk of Q nk ofmaximal eigenvalue and draw q n +1 k ∼ M q nk . Table 3:
One iteration of the time-discrete algorithmFinally, we would like to stress that the quaternion formulation is not only aconvenient numerical trick. The equivalence it provides between body-orientationmodels and models of nematic alignment of polymers in dimension four has beenexploited in [32] to study phase transitions in the body alignment model.
D Derivation of the macroscopic model
The derivation of the continuum theory presented in Section 2.2 has been achieved in[38] (see also [32]) following earlier works [35, 37]. It consists of two steps. The firststep is the derivation of a mean-field kinetic model in the limit N → ∞ showingthat the system satisfies the propagation of chaos property: the agents, seen asrandom variables in R × SO ( R ) become independent and identically distributed.Their law is given by the kinetic particle distribution f which satisfies the followingPDE: ∂ t f + c A e · ∇ x f = ν ( ρ f M A K ∗ f − f ) , where ρ f ≡ ρ f ( t, x ) is the local spatial density: ρ f ( t, x ) = (cid:90) SO ( R ) f ( t, x , A ) d A, 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 ( R ) A · J K ∗ f ( t, x ) , computed from the local flux: J K ∗ f ≡ J K ∗ f ( t, x ) := (cid:90) (cid:90) R × SO ( R ) K ( x − y ) A f ( t, y , A ) d y d A. rom a mathematical point of view, the probability distribution f ≡ f ( t, x , A ) isobtained as the limit in law of the empirical measure of the N -particle system. Werefer to [43] where a rigorous proof of this result is presented for a similar model,and to [10] for a related work on the Vicsek model.In the macroscopic regime the agent interactions become strong, which is ex-pressed by the following hydrodynamic scaling: ε ∼ c ν L ∼ RL (cid:28) , where L is a typical macroscopic length-scale of the system (such as the typical size ofthe flock). We define ˜ c = ενL = O (1) and c (cid:48) = c / ˜ c . Then, defining dimensionlesstime and space variables t (cid:48) and x (cid:48) such that x = L x (cid:48) and t = ( L/ ˜ c ) t (cid:48) , we obtain(dropping the primes for simplicity): ∂ t f ε + c A e · ∇ x f ε = 1 ε ( ρ f ε M A fε − f ε ) + O ( ε ) , (53)where A f ε ≡ A f ε ( t, x ) := arg max A ∈ SO ( R ) A · J f ε ( t, x ) , and J f ε ≡ J f ε ( t, x ) := (cid:90) SO ( R ) A f ε ( t, x , A ) d A. This last expression is obtained by Taylor expanding J K ∗ f ε = J f ε + O ( ε ) andmeans that the interactions between the agents become spatially localized in themacroscopic regime.The macroscopic model is obtained by formally taking the limit ε → f ε −→ ε → ρ M A (54)where ρ ≡ ρ ( t, x ) and A ≡ A ( t, x ) depend on t and x . Thus, the limiting distributionis fully described by the spatial density of agents and their average orientation. Toobtain a system of equations for ( ρ, A ), we first use the local conservation of mass:integrating (53) over SO ( R ) and noting the right-hand side vanishes, it holds that, ∂ t (cid:90) SO ( R ) f ε d A + c (cid:90) SO ( R ) A e · ∇ x f ε d A = O ( ε ) . When ε →
0, assuming (54) and using (38), we obtain (11a).To obtain an equation for A , it could be tempting to pursue this approach andmultiply (53) by A before integrating it over SO ( R ). However, the term resultingfrom the right-hand side of (53) does not vanish but equals (using (38) again):1 ε (cid:90) SO ( R ) A ( ρ f ε M A fε − f ε ) d A = 1 ε (cid:16) c c ρ f ε A f ε − J f ε (cid:17) (cid:54) = 0 . Due to the factor ε − , its limit as ε → A before integrating it over SO ( R ), we multiply t by the quantity ψ A fε ( A ) := A T f ε A − A T A f ε . The rationale for using this quantity isbecause we aim to find an equation for the time-derivative of A . Such a derivativemust lie in the tangent space to SO ( R ) at A , denoted by T A . This suggests tomultiply (53) by an element of T A . Given an arbitrary matrix A , a natural way toobtain an element of T A is to take its orthogonal projection on T A , which is givenby ( A − A A T A ). We could therefore choose to multiply (53) by this quantity. Buta further simplification is possible by noting that this quantity is equal to A ψ A ( A )and that A does not depend on A and so can be factored out of the integral withrespect to A . These considerations naturally lead to the choice of the antisymmetricmatrix ψ A fε ( A ) as a multiplier. Because A f ε is obtained as the polar decompositionof J f ε , there exists a symmetric matrix S such that J f ε = A f ε S . Using this remarkand (38), we easily find that1 ε (cid:90) SO ( R ) ψ A fε ( A ) ( ρ f ε M A fε − f ε ) d A = 0 . Then, multiplying (53) by ψ f ε , taking the limit ε → (cid:90) SO ( R ) ( ∂ t ( ρM A ) + c A e · ∇ x ( ρM A )) ψ A ( A ) d A = 0 . Eq. (11b) of the SOHB model follows from this equation through tedious butstraightforward computations detailed in [35, 37].Note that the simple form of the multiplier ψ A f is due to a particular simpleexpression of the collision operator. In more general cases, the obtention of themultiplier (referred to as the generalised collision invariant in [41]) is more involved(see e.g. [35, 37, 38]). A rigorous convergence result for the limit ε → E MO and generalized HW solutions
In this section, we provide proofs of Lemmas 3.1 and 3.2. The prototypical helicaltraveling wave (HW) presented in Lemma 3.2 belongs to a more general class ofsolutions described in Section E.2 below.
E.1 Proof of Lemma 3.1
Starting from the initial condition (28), we are looking for solutions of (11b) of theform A ( t, x ) = cos( ωt ) u ( t, z ) v ( t, z ) − sin( ωt ) u ( t, z ) v ( t, z )0 u ( t, z ) v ( t, z ) , where ω ∈ R is an angular velocity which will be related to the parameters of theproblem later and where the basis vectors u = ( u , u , u ) T and v = ( v , v , v ) T depend only on the z variable and time. In this situation, Equation 11a is triviallysatisfied which means that the system stays homogeneous in space. Solutions of this orm have to satisfy three geometrical constraints which ensure that A ∈ SO ( R ).The first two ones are Ω × u = v and v × Ω = u , which lead to A ( t, x ) = cos( ωt ) sin( ωt ) v ( t, z ) sin( ωt ) u ( t, z ) − sin( ωt ) cos( ωt ) v ( t, z ) − cos( ωt ) u ( t, z )0 u ( t, z ) v ( t, z ) . (55)The third one is a normalization constraint: ∀ t > , ∀ z ∈ R , u ( t, z ) + v ( t, z ) = 1 . (56)Using (56), we define a function α ≡ α ( t, z ) such that u ( t, z ) = sin( α ( t, z )) , v ( t, z ) = cos( α ( t, z )) . A direct computation shows that for A of the form (55), we have r = ( ∂ z u ) u + ( ∂ z v ) v , δ = 0 . Therefore, Eq. (11b) can be rewritten more concisely into: ∂ t A + c [Ω × r ] × A = 0 , (57)where we recall Eq. (9) for the definition of [ ] × . A direct computation shows thatΩ × r = ( v ∂ z u − u ∂ z v ) e = ( ∂ z α ) e . (58)Inserting this in (57) implies that u ( t, z ) ≡ u ( z ) and v ( t, z ) ≡ v ( z ) are indepen-dent of time. We then observe that: A ( t, x ) = A ( − ωt, e ) A ( α ( z ) , e ) , (59)where we recall Eq. (8) for the meaning of A . Therefore, using (57) and (58), weobtain: − ω [ e ] × A + c ( ∂ z α ) [ e ] × A = 0 , from which we deduce that A satisfies (11b) if and only if α and ω satisfy: c ∂ z α = ω, which implies α ( z ) = ωc z + ¯ α, (60)where ¯ α is a constant, which can be interpreted as the phase at the origin z = 0.To recover (26), we just need to take ¯ α = 0 and define ξ = ω/c . Eq. (27) followsfrom (59). .2 Generalized HW and proof of Lemma 3.2 Starting from the initial condition (31), we are looking for solutions of (11b) of theform A ( t, x ) = α ( 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 constantvector and Equation 20a is trivially satisfied. Moreover a direct computation showsthat: r = 0 , δ = ( ∂ x α )( t, x ) . As a consequence, Eq. (23) is trivially satisfied and straightfoward computationsshow that Eq. (11b) reduces to ∂ t α + ( c + c ) ∂ x α = 0 . This last equation is a linear transport equation with velocity c + c , the solutionsof which are given by α ( t, x ) = α ( x − ( c + c ) t ) (61)for any initial condition α ∈ L ( R ). In the case of (31), α ( x ) = ξ x . However, wesee that there are as many different solutions as functions in L ( R ). Such generalsolutions are called “generalized HW”. E.3 GOP of the MO and generalized HW
The GOP (given by Eq. (39)) of the MO and HW do not depend on time and onlydepend on the function α defined respectively by (60) and (61). Using Eq. (39), wecan compute that the GOP is equal to:GOP = 12 (cid:18) c ( κ ) c (cid:19) (cid:0) |(cid:104) u (cid:105)| (cid:1) + 14 , where (cid:104) u (cid:105) denotes the spatial average of the vector u with respect to ρ (here thewith respect to the uniform measure on the domain since ρ is constant and uniform).With the previous notations, we obtain |(cid:104) u (cid:105)| = (cid:104) cos α (cid:105) + (cid:104) sin α (cid:105) , For the generalized HW, depending on the choice of α , the GOP can take any valuebetween GOP and GOP , these two extreme values being attained respectivelywhen |(cid:104) u (cid:105)| = 0 and |(cid:104) u (cid:105)| = 1. References [1] P. Aceves-Sanchez, M. Bostan, J.-A. Carrillo, and P. Degond. Hydrodynamiclimits for kinetic flocking models of Cucker-Smale type.
Math. Biosci. Eng. ,16:7883–7910, 2019.
2] M. Aldana, H. Larralde, and B. V´azquez. On the emergence of collective orderin swarming systems: a recent debate.
Int. J. Mod. Phys. B , 23(18):3661–3685,2009.[3] I. Aoki. A simulation study on the schooling mechanism in fish.
Bull. Japan.Soc. Sci. Fish , 48:1081–1088, 1982.[4] A. Barbaro and P. Degond. Phase transition and diffusion among sociallyinteracting self-propelled agents.
Discrete Contin. Dyn. Syst. Ser. B , 19:1249–1278, 2014.[5] A. B. Barbaro, J. A. Canizo, J. A. Carrillo, and P. Degond. Phase transitionsin a kinetic flocking model of Cucker–Smale type.
Multiscale Model. Simul. ,14(3):1063–1088, 2016.[6] J. Barr´e, C. Bernardin, R. Ch´etrite, Y. Chopra, and M. Mariani. Gammaconvergence approach for the large deviations of the density in systems ofinteracting diffusion processes. arXiv preprint arXiv:1910.04026 , 2019.[7] L. Berlyand, R. Creese, P.-E. Jabin, and M. Potomkin. Continuum approxima-tions to systems of correlated interacting particles.
J. Stat. Phys. , 174(4):808–829, 2019.[8] E. Bertin, M. Droz, and G. Gr´egoire. Boltzmann and hydrodynamic descriptionfor self-propelled particles.
Phys. Rev. E , 74(2):022101, 2006.[9] E. Bertin, M. Droz, and G. Gr´egoire. Hydrodynamic equations for self-propelled particles: microscopic derivation and stability analysis.
J. Phys.A , 42(44):445001, 2009.[10] F. Bolley, J. A. Ca˜nizo, and J. A. Carrillo. Mean-field limit for the stochasticVicsek model.
Appl. Math. Lett. , 25(3):339–343, 2012.[11] L. Bortolussi and N. Gast. Mean-field limits beyond ordinary differential equa-tions. In
International School on Formal Methods for the Design of Computer,Communication and Software Systems , pages 61–82. Springer, 2016.[12] M. Bostan and J. A. Carrillo. Asymptotic fixed-speed reduced dynamics forkinetic equations in swarming.
Math. Models Methods Appl. Sci. , 23(13):2353–2393, 2013.[13] M. Bostan and J. A. Carrillo. Reduced fluid models for self-propelled particlesinteracting through alignment.
Math. Models Methods Appl. Sci. , 27(07):1255–1299, 2017.[14] M. Briant and S. Merino-Aceituno. Cauchy theory and mean-field limit forgeneral Vicsek models in collective dynamics. arXiv preprint arXiv:2004.00883 ,2020.[15] A. Bricard, J.-B. Caussin, D. Das, C. Savoie, V. Chikkadi, K. Shitara, O. Chep-izhko, F. Peruani, D. Saintillan, and D. Bartolo. Emergent vortices in popula-tions of colloidal rollers.
Nat. Commun. , 6:7470, 2015.
16] D. S. Calovi, U. Lopez, S. Ngo, C. Sire, H. Chat´e, and G. Theraulaz. Swarming,schooling, milling: phase diagram of a data-driven fish school model.
New J.Phys. , 16(1):015026, 2014.[17] J. A. Carrillo, M. R. D’Orsogna, and V. Panferov. Double milling in self-propelled swarms from kinetic theory.
Kinet. Relat. Models , 2(2):363, 2009.[18] C. Castellano, S. Fortunato, and V. Loreto. Statistical physics of social dy-namics.
Rev. Modern Phys. , 81(2):591, 2009.[19] J.-B. Caussin, A. Solon, A. Peshkov, H. Chat´e, T. Dauxois, J. Tailleur,V. Vitelli, and D. Bartolo. Emergent spatial structures in flocking models:a dynamical system insight.
Phys. Rev. Lett. , 112(14):148102, 2014.[20] A. Cavagna, L. Del Castello, I. Giardina, T. Grigera, A. Jelic, S. Melillo,T. Mora, L. Parisi, E. Silvestri, M. Viale, et al. Flocking and turning: anew model for self-organized collective motion.
J. Stat. Phys. , 158(3):601–627,2015.[21] C. Cercignani, R. Illner, and M. Pulvirenti.
The Mathematical Theory of DiluteGases , volume 106. Springer Science & Business Media, 2013.[22] B. Charlier, J. Feydy, J. A. Glaun`es, F.-D. Collin, and G. Durif. Kernel oper-ations on the GPU, with autodiff, without memory overflows. arXiv preprintarXiv:2004.11127 , 2020.[23] H. Chat´e, F. Ginelli, G. Gr´egoire, and F. Raynaud. Collective motion of self-propelled particles interacting without cohesion.
Phys. Rev. E , 77(4):046113,2008.[24] A. Costanzo and C. Hemelrijk. Spontaneous emergence of milling (vortex state)in a Vicsek-like model.
J. Phys. D: Appl. Phys. , 51(13):134004, 2018.[25] I. D. Couzin and N. R. Franks. Self-organized lane formation and optimizedtraffic flow in army ants.
Proc. Biol. Sci. , 270(1511):139–146, 2003.[26] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks. Collectivememory and spatial sorting in animal groups.
Journal of theoretical biology ,218(1):1–12, 2002.[27] A. Creppy, F. Plourabou´e, O. Praud, X. Druart, S. Cazin, H. Yu, and P. De-gond. Symmetry-breaking phase transitions in highly concentrated semen.
J.R. Soc. Interface , 13(123):20160575, 2016.[28] F. Cucker and S. Smale. Emergent behavior in flocks.
IEEE Trans. Automat.Control , 52(5):852–862, 2007.[29] A. Czir´ok, E. Ben-Jacob, I. Cohen, and T. Vicsek. Formation of complexbacterial colonies via self-generated vortices.
Phys. Rev. E , 54(2):1791, 1996.[30] D. A. Dawson and J. G¨artner. Large deviations from the McKean-Vlasov limitfor weakly interacting diffusions.
Stochastics , 20(4):247–308, 1987.
31] P. Degond. Macroscopic limits of the Boltzmann equation: a review. In P. De-gond, L. Pareschi, and G. Russo, editors,
Modeling and Computational Methodsfor Kinetic Equations , Modeling and Simulation in Science, Engineering andTechnology, pages 3–57. Birkh¨auser Basel, 2004.[32] P. Degond, A. Diez, A. Frouvelle, and S. Merino-Aceituno. Phase transitionsand macroscopic limits in a BGK model of body-attitude coordination.
J.Nonlinear Sci. , 30:2671–2736, 2020.[33] P. Degond, A. Frouvelle, and J.-G. Liu. Macroscopic limits and phase transitionin a system of self-propelled particles.
J. Nonlinear Sci. , 23(3):427–456, 2013.[34] P. Degond, A. Frouvelle, and J.-G. Liu. Phase transitions, hysteresis, and hy-perbolicity for self-organized alignment dynamics.
Arch. Ration. Mech. Anal. ,216(1):63–115, 2015.[35] P. Degond, A. Frouvelle, and S. Merino-Aceituno. A new flocking modelthrough body attitude coordination.
Math. Models Methods Appl. Sci. ,27(06):1005–1049, 2017.[36] P. Degond, A. Frouvelle, S. Merino-Aceituno, and A. Trescases. Some proper-ties of Self-Organized Hydrodynamics for body orientation.
In preparation .[37] P. Degond, A. Frouvelle, S. Merino-Aceituno, and A. Trescases. Quaternionsin collective dynamics.
Multiscale Model. Simul. , 16(1):28–77, 2018.[38] P. Degond, A. Frouvelle, S. Merino-Aceituno, and A. Trescases. Alignment ofself-propelled rigid bodies: from particle systems to macroscopic equations. InG. Giacomin, S. Olla, E. Saada, H. Spohn, and G. Stoltz, editors,
StochasticDynamics Out of Equilibrium , volume 282 of
Springer Proceedings in Mathe-matics and Statistics , pages 28–66. Institut Henri Poincar´e, Paris, France, 2017,Springer International Publishing, 2019.[39] P. Degond, J.-G. Liu, S. Merino-Aceituno, and T. Tardiveau. Continuum dy-namics of the intention field under weakly cohesive social interaction.
Math.Models Methods Appl. Sci. , 27(01):159–182, 2017.[40] P. Degond, J.-G. Liu, S. Motsch, and V. Panferov. Hydrodynamic models ofself-organized dynamics: derivation and existence theory.
Methods Appl. Anal. ,20:89–114, 2013.[41] P. Degond and S. Motsch. Continuum limit of self-driven particles with ori-entation interaction.
Math. Models Methods Appl. Sci. , 18(supp01):1193–1215,2008.[42] P. Degond and S. Motsch. A macroscopic model for a system of swarmingagents using curvature control.
J. Stat. Phys. , 143(4):685–714, 2011.[43] A. Diez. Propagation of chaos and moderate interaction for a piecewise de-terministic system of geometrically enriched particles.
Electron. J. Probab. ,25(90), 2020.
44] G. Dimarco and S. Motsch. Self-alignment driven by jump processes : Macro-scopic limit and numerical investigation.
Math. Models Methods Appl. Sci. ,26(07):1385–1410, 2016.[45] M. R. D’Orsogna, Y.-L. Chuang, A. L. Bertozzi, and L. S. Chayes. Self-propelled particles with soft-core interactions: patterns, stability, and collapse.
Phys. Rev. Lett. , 96(10):104302, 2006.[46] B. Fernandez and S. M´el´eard. A Hilbertian approach for fluctuations on theMcKean-Vlasov model.
Stochastic Process. Appl. , 71(1):33–53, 1997.[47] A. Figalli, M.-J. Kang, and J. Morales. Global well-posedness of the spa-tially homogeneous Kolmogorov–Vicsek model as a gradient flow.
Arch. Ration.Mech. Anal. , 227(3):869–896, 2018.[48] A. Frouvelle and J.-G. Liu. Dynamics in a kinetic model of oriented particleswith phase transition.
SIAM J. Math. Anal. , 44(2):791–826, 2012.[49] I. M. Gamba, J. R. Haack, and S. Motsch. Spectral method for a kineticswarming model.
J. Comput. Phys. , 297:32–46, 2015.[50] I. M. Gamba and M.-J. Kang. Global weak solutions for Kolmogorov–Vicsektype equations with orientational interactions.
Arch. Ration. Mech. Anal. ,222(1):317–342, 2016.[51] J. Gautrais, F. Ginelli, R. Fournier, S. Blanco, M. Soria, H. Chat´e, and G. Ther-aulaz. Deciphering interactions in moving animal groups.
PLoS Comput. Biol. ,2012.[52] P. Gerlee, K. Tunstrøm, T. Lundh, and B. Wennberg. Impact of anticipationin dynamical systems.
Phys. Rev. E , 96:062413, Dec 2017.[53] Q. Griette and S. Motsch. Kinetic equations and self-organized band forma-tions. In
Active Particles, Volume 2 , pages 173–199. Springer, 2019.[54] S.-Y. Ha, J.-G. Liu, et al. A simple proof of the Cucker-Smale flocking dynamicsand mean-field limit.
Commun. Math. Sci. , 7(2):297–325, 2009.[55] S.-Y. Ha and E. Tadmor. From particle to kinetic and hydrodynamic descrip-tions of flocking.
Kinet. Relat. Models , 1:415–435, 2008.[56] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virta-nen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern,M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. delR´ıo, M. Wiebe, P. Peterson, P. G´erard-Marchant, K. Sheppard, T. Reddy,W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant. Array programmingwith NumPy.
Nature , 585(7825):357–362, 2020.[57] M. Z. Hasan and C. L. Kane. Colloquium: topological insulators.
Rev. ModernPhys. , 82(4):3045, 2010.
58] C. K. Hemelrijk and H. Hildenbrandt. Schools of fish and flocks of birds: theirshape and internal structure by self-organization.
Interface Focus , 2(6):726–737, Aug 2012.[59] C. K. Hemelrijk, H. Hildenbrandt, J. Reinders, and E. J. Stamhuis. Emer-gence of oblong school shape: models and empirical data of fish.
Ethology ,116(11):1099–1112, 2010.[60] H. Hildenbrandt, C. Carere, and C. K. Hemelrijk. Self-organized aerial displaysof thousands of starlings: a model.
Behavioral Ecology , 21(6):1349–1359, 2010.[61] J. F. Hughes, A. van Dam, M. McGuire, D. F. Sklar, J. D. Foley, S. K. Feiner,and A. Kurt.
Computer Graphics: Principles and Practice . Addison-WesleyProfessional, 3rd edition edition, 2013.[62] J. D. Hunter. Matplotlib: A 2D graphics environment.
Comput Sci Eng. ,9(3):90–95, 2007.[63] D. Q. Huynh. Metrics for 3D rotations: Comparison and analysis.
J. Math.Imaging Vis. , 35(2):155–164, 2009.[64] N. Jiang, L. Xiong, and T.-F. Zhang. Hydrodynamic limits of the kineticself-organized models.
SIAM J. Math. Anal. , 48(5):3383–3411, 2016.[65] J. T. Kent, A. M. Ganeiber, and K. V. Mardia. A new unified approach forthe simulation of a wide class of directional distributions.
J. Comput. Graph.Statist. , 27(2):291–301, 2018.[66] K. v. Klitzing, G. Dorda, and M. Pepper. New Method for High-AccuracyDetermination of the Fine-Structure Constant Based on Quantized Hall Resis-tance.
Phys. Rev. Lett. , 45(6):494–497, Aug 1980.[67] C. Lancellotti. On the fluctuations about the Vlasov limit for N-particle systemswith mean-field interactions.
J. Stat. Phys. , 136(4):643–665, 2009.[68] R. B. Laughlin. Quantized Hall conductivity in two dimensions.
Phys. Rev. B ,23(10):5632–5633, May 1981.[69] T. Lee. Bayesian attitude estimation with the matrix Fisher distribution onSO(3).
IEEE Trans. Automat. Contr. , 63(10):3377–3392, 2018.[70] R. Lukeman, Y.-X. Li, and L. Edelstein-Keshet. Inferring individual rules fromcollective behavior.
Proc. Natl. Acad. Sci. USA , 107(28):12576–12580, 2010.[71] K. V. Mardia and P. E. Jupp.
Directional Statistics , volume 494. John Wiley& Sons, 2009.[72] A. Mart´ın-G´omez, D. Levis, A. D´ıaz-Guilera, and I. Pagonabarraga. Collec-tive motion of active brownian particles with polar alignment.
Soft Matter ,14(14):2610–2618, 2018.
73] S. Motsch and L. Navoret. Numerical simulations of a nonconservative hyper-bolic system with geometric constraints describing swarming behavior.
Multi-scale Model. Simul. , 9(3):1253–1275, 2011.[74] S. Motsch and E. Tadmor. A new model for self-organized dynamics and itsflocking behavior.
J. Stat. Phys. , 144(5):923, 2011.[75] Nobel Foundation. The Nobel Prize in Physics 1985. .[76] Nobel Foundation. The Nobel Prize in Physics 2016. .[77] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen,Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. De-Vito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, andS. Chintala. PyTorch: An imperative style, high-performance deep learninglibrary. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alch´e Buc, E. Fox,and R. Garnett, editors,
Advances in Neural Information Processing Systems32 , pages 8024–8035. Curran Associates, Inc., 2019.[78] F. Peruani, A. Deutsch, and M. B¨ar. A mean-field theory for self-propelledparticles interacting by velocity alignment mechanisms.
Eur. Phys. J. Spec.Top. , 157(1):111–122, 2008.[79] X.-L. Qi and S.-C. Zhang. Topological insulators and superconductors.
Rev.Modern Phys. , 83(4):1057, 2011.[80] D. Scherer, P. Dubois, and B. Sherwood. VPython: 3D interactive scientificgraphics for students.
Comput Sci Eng. , 2(5):56–62, 2000.[81] G. Schou. Estimation of the concentration parameter in von Mises–Fisherdistributions.
Biometrika , 65(2):369–377, 1978.[82] S. Shankar, M. J. Bowick, and M. C. Marchetti. Topological sound and flockingon curved surfaces.
Phys. Rev. X , 7(3):031039, 2017.[83] K. Sone and Y. Ashida. Anomalous topological active matter.
Phys. Rev. Lett. ,123:205502, Nov 2019.[84] A. Souslov, B. C. Van Zuiden, D. Bartolo, and V. Vitelli. Topological soundin active-liquid metamaterials.
Nature Phys. , 13(11):1091, 2017.[85] D. J. Thouless. Quantization of particle transport.
Phys. Rev. B , 27(10):6083–6087, May 1983.[86] J. Toner and Y. Tu. Flocks, herds, and schools: A quantitative theory offlocking.
Phys. Rev. E , 58(4):4828, 1998.[87] R. Vershynin. Introduction to the non-asymptotic analysis of random matri-ces. In Y. C. Eldar and G. Kutinyok, editors,
Compressed sensing, theory andapplications , pages 210–260. Cambridge University Press, 2012.
88] T. Vicsek, A. Czir´ok, E. Ben-Jacob, I. Cohen, and O. Shochet. Novel typeof phase transition in a system of self-driven particles.
Phys. Rev. Lett. ,75(6):1226, 1995.[89] T. Vicsek and A. Zafeiris. Collective motion.
Phys. Rep. , 517(3-4):71–140,2012.[90] T.-F. Zhang and N. Jiang. A local existence of viscous self-organized hydrody-namic model.
Nonlinear Anal. Real World Appl. , 34:495–506, 2017., 34:495–506, 2017.