SELF-ORGANIZED HYDRODYNAMICS WITH DENSITY-DEPENDENT VELOCITY

. Motivated by recent experimental and computational results that show a motility-induced clustering transition in self-propelled particle systems, we study an individual model and its corresponding Self-Organized Hydrodynamic model for collective behaviour that incorporates a density-dependent velocity, as well as inter-particle alignment. The modal analysis of the hydrodynamic model elucidates the relationship between the stability of the equilibria and the changing velocity, and the formation of clusters. We ﬁnd, in agreement with earlier results for non-aligning particles, that the key criterion for stability is ( ρv ( ρ )) (cid:48) ≥ 0, i.e. a nondecreasing mass ﬂux ρv ( ρ ) with respect to the density. Numerical simulation for both the individual and hydrodynamic models with a velocity function inspired by experiment demonstrates the validity of the theoretical results.


Dedicated to Peter Markowich to celebrate 30 years of friendship.
Abstract. Motivated by recent experimental and computational results that show a motility-induced clustering transition in self-propelled particle systems, we study an individual model and its corresponding Self-Organized Hydrodynamic model for collective behaviour that incorporates a density-dependent velocity, as well as inter-particle alignment. The modal analysis of the hydrodynamic model elucidates the relationship between the stability of the equilibria and the changing velocity, and the formation of clusters. We find, in agreement with earlier results for non-aligning particles, that the key criterion for stability is (ρv(ρ)) ≥ 0, i.e. a nondecreasing mass flux ρv(ρ) with respect to the density. Numerical simulation for both the individual and hydrodynamic models with a velocity function inspired by experiment demonstrates the validity of the theoretical results.
1. Introduction. The study of flocking is inspired by the natural behaviour of animal groups, such as flocks of birds and schools of fishes. Natural flocks exhibit a range of states, including moving swarms, compact flocks, correlated turning and enhanced density fluctuations. To capture flock properties, numerical flocking models such as the Vicsek model [29] have been developed, and subsequently been studied in great detail; see e.g. [5,6,8]. The Vicsek model exhibits a complex first order transition [5] between an aligned flocking state at low noise levels and a disordered state at high noise levels, via a band state that depends sensitively on the detailed implementation [22]. In parallel, hydrodynamic models of the aligned state have been proposed [28,12], which show that the aligned state exhibits critical scaling fluctuations. In particular, enhanced transverse diffusion is responsible for stabilizing true long-ranged order in the Vicsek model, in contrast to the quasi-longranged order found in the XY model and cemented in the Mermin-Wagner theorem [4].
Missing from the Vicsek model are the effects of excluded volume, and repulsion or attraction between individual animals or agents. A flurry of recent numerical, analytical and experimental work in the physics community has begun to investigate the effects of non point-like agents. Using soft self-propelled particles, that is a particle model where in addition to short-range repulsive forces, self-propulsion is introduced as a force into fully overdamped Langevin dynamics, several groups showed [17,23] that the mix of self-propulsion and volume exclusion has a profound effect on the system properties. These results were first obtained for non-aligning active particles; in this paper we will investigate the effect of additional alignment. At intermediate densities, the chief effect of volume exclusion is a slowdown of the effective hydrodynamic velocity v(ρ) where ρ is the density and v (ρ) < 0. Experiments with artificially created active colloids [21] also show a pronounced clustering tendency, likely due to the same slowdown of velocity with density. Pedestrian dynamics is one practical area of application with the combination of alignment and repulsion between actors that we focus on here. Encouragingly, the empirical flow-density relation in a range of crowd situations shows the same decrease with density that was observed numerically [25,24] in what is a field where experiments are very difficult. The flow-density relationship in pedestrian dynamics, (aka the fundamental diagram), displays nonmonotonic behaviour of the flux ρv that first increases with respect to ρ in the free-flow branch and then decreases due to the formation of jams.
Crowding effects in self-propelled aligning particles have previously been taken into account by various ways. A first one is by enhancing the pressure as a result of crowding; see e.g. the derivation of the corresponding Self-Organized Hydrodynamic (SOH) model in [9] and the application to crowd dynamics in [11]. A second one consists of supposing that the actual velocity of the particle has a component resulting from repulsion by other particles that adds up to the self-propulsion velocity [26,19]. This second approach leads to an SOH model with other transport terms [1,10,14]. The two approaches have been compared in [10]. Here, the present approach aims at mimicking the effect observed in [1,10] with a simpler model allowing for analytical investigation of its stability. Another way of taking into account crowding in randomly diffusive particle systems is by reducing the diffusion coefficient in the crowded areas. Such a model has been derived in [1] and has been used for modelling cross-diffusion in [2]. The investigation of a density-dependent diffusivity on self-propelled aligning particles has been investigated in [18] in the framework of the SOH model. It has been shown that a decreasing angular diffusivity leads to a loss of hyperbolicity of the corresponding SOH model, which is a signature of clustering.
Fits to simulations of self-propelled hard and soft particles and collision-based models suggest a universal form at low and intermediate densities where c depends only on the Péclet number P e = v 0 /aν r , where a is the particle radius and ν r is the rotational diffusion constant, or more generally the ratio of persistence rate to diffusion rate. This density-dependent velocity leads to a density instability, and finally to a clustering transition where the system phase-separates into a (single) cluster and a low density gas phase. This motility induced phase separation (MIPS) transition appears to be of a spinodal decomposition type, in a direct analogy to the liquidgas transition. The transition line is determined by the Péclet number, and it is hypothesized that it terminates in a critical point around P e ≈ 10 [23]. At high density, a second transition branch separates the cluster phase from a dense liquid phase, and ultimately a high density, low driving active glassy phase [16].
Analytically, MIPS was first proposed in a one-dimensional model of interacting run-and-tumble particles [27]. By mapping the Fokker-Planck equation onto an equivalent equilibrium equation, Tailleur and Cates were able to define an effective free energy with a spinodal transition analogous to the liquid-gas transition. This theoretical approach was later extended to fully brownian particles and tested numerically [3].
The effect of alignment on MIPS was first studied by Farrell et al. [15] using a combination of hydrodynamic equations derived from a microscopic particle model, and direct numerical simulation. The phase diagram contains both homogeneous and MIPS phases, but also travelling bands and finite clusters. However, a full understanding of the phase separation mechanism in the presence of alignment is still lacking.
Instead of the usually constant speed, we introduce a density-dependent velocity v(ρ) to the Vicsek model with alignment between individuals, and study the dynamics of the high density system both through direct simulation and the Self-Organized Hydrodynamic (SOH) formalism. In this paper, we focus on the deeply aligned phase, and study in detail the location of the instability line and its angular dependence. In the unstable phase, we determine the unstable eigenmode as a function of wave vector and orientation, and determine its growth rate with perturbations. These results are then compared to a numerical solution of the full SOH equations, and a direct solution of the particle model. This paper is organized as follows: Section 2 introduces the particle model with the density dependent velocity. Section 3 presents the derived hydrodynamic model and studies the stability of the inviscid and viscous cases. Finally, Section 4 presents numerical results from both the particle and hydrodynamic models and ends with a discussion of the growth rate of the instability. The appendices detail the derivation of the hydrodynamic model (Sections A.1 and A.2), and the numerical scheme used to integrate the SOH equations (Section B). Finally we detail in Section C the discrete Fourier transform used for the numerical evaluation of the growth rate of the instability.
To derive the SOH equations, we employ the mathematical theory developed in [12] and provide the proof in Appendix A.1 and A.2. One of the key components of the proof is the concept of "Generalized Collision Invariance" (GCI) which allows passing to the hydrodynamic limit in spite of the lack of momentum conservation at the particle level. Further generalization and elaboration of this method can be found in [18,9].
Our macroscopic model describes two quantities, the density ρ(x, t) ≥ 0 and the mean orientation Ω(x, t) at position x ∈ R n and time t ≥ 0, and is referred to as the Self-Organized Hydrodynamic (SOH) model: where c 1 , c 2 and d are dimensionless mobility parameters coarse-grained from the ones at the particle level. γ has the dimensions of a diffusion constant and expresses that polarization diffuses as a result of microscopic alignment interactions. v(ρ) is the function which specifies the relation between the speed and the density, and therefore non-negative. Note that the projection operator P Ω ⊥ = Id − Ω ⊗ Ω will preserve the geometric constraint |Ω| = 1. Using a two-dimensional system, we show below that the stability of the equilibrium is related to the behaviour of the mass flux ρv(ρ). More precisely, the inviscid model (γ = 0) is hyperbolic and hence stable against perturbations if (ρv(ρ)) , the derivative of ρv(ρ) with respect to ρ, is non-negative or certain constraints on the equilibria are satisfied. The viscous model (γ = 0) is stable if and only if (ρv(ρ)) ≥ 0.
This result shows that a density instability associated to clustering emerges below a threshold mass flux derivative and also that we obtain steady states with constant density for mass fluxes above this threshold. Our result agrees with the conclusions of [27,17,23,3,16] for the soft particle and hydrodynamics models where self-propulsion decreases with local density. In common with Farrell et al. [15], we find that the alignment and clustering transitions are largely independent of each other; however we find correlations between orientation of a perturbation and stability which have not previously been explored. This analysis on the SOH model is supported by numerical simulations with forms of v(ρ) inspired by experiment for both the individual-based model and the SOH model. In particular, in the high concentration limit, we can decompose the instability as a sum of unstable eigenmodes, the growth rate of which we then study numerically. Given a steady state for the density, we observe a positive correlation between the growth rate and the stationary polarisation orientation for both the particle and hydrodynamic models.

2.
Particle model with density-dependent velocity. In this section, we introduce the particle model which is the starting point for the derivation of the SOH model (1). Consider a system of N self-propelled particles in R n . Let t be the time, X i (t) the position of the i-th particle and ω i (t) its velocity orientation. Then the time evolution for the i-th particle is given by Here P ω ⊥ i represents the projection on the plane which is perpendicular to ω i , which allows the geometrical constraint on ω i (2c) to be preserved. B i t is a Brownian motion with noise strength D andω i is the mean velocity oriention in the neighborhood where K 1 (X) is the kernel for the alignment with a compact support in B R1 , and the constant parameter ν is the alignment rate.
The new ingredient of this individual-based model compared to the Vicsek model is the dependency of the velocity on the mass in the neighborhood, described by the function v(m i ) where m i is the density of the neighbourhood of X i . Let R 2 be the interaction range and |B R2 | the volume of the ball Then m i is defined as where K 2 (X) is a kernel defined on a compact support in B R2 . For instance, K 2 can be chosen as χ Xi (|X i − X|), the characteristic function on B R2 (X i ). An intuitive motivation for the density-dependent velocity is as follows: Particles that are not point-like will experience collisions. While in thermal systems this just randomises directions, for active, persistent motion the combination of collision and persistence just slows the particle down, much like a crowd slows a pedestrian who passes through. A simple scaling argument shows that the collision rate is proportional to density, leading to v(ρ) ≈ v 0 (1−cρ). For aligning systems, like here, the situation is a bit more complex since colliding particles will also align, eventually removing many of the collisions if the system has strong overall polarisation. Here, we keep a generic form of v(ρ), and investigate the onset of instability as a function of wave vector and angle with the polarisation direction.
3. Linearized stability analysis of the SOH model with density-dependent velocity. The derivation of the SOH model (1) is analogous with the work in [12] and [10]; for the details please see Appendix A.
We derive the linearized stability criteria for the inviscid and viscous cases in a two-dimensional space in the following subsections. Consider the SOH model for x = (x, y) ∈ R 2 . Since |Ω| = 1, we define Ω(x, y, t) = (cos θ(x, y, t), sin θ(x, y, t)) through the angle function θ(x, y, t) of the vector Ω. Then the projection operator P Ω ⊥ becomes a matrix operator sin 2 θ − sin θ cos θ − sin θ cos θ cos 2 θ . Without loss of generalization, we scale the system (1) such that c 1 = 1, and arrive at a system of ρ and θ: where the two matrices A x (ρ, θ) and A y (ρ, θ) are given by andṽ(ρ) = ρv(ρ) andṽ (ρ) is the derivative ofṽ with respect to ρ.
To simplify, suppose that (ρ, θ) depends only on x, i.e., we are interested in the propagation of waves with arbitrary orientation Ω in the horizontal direction. Let (ρ s (x), θ s (x)) denote the equilibrium solutions and they must satisfy the following system: Let us then expand around them with a small perturbation parameter σ: Dropping the higher order terms O(σ 2 ) and using (ρ σ , θ σ ) to represent the first order perturbation, we arrive at the system below 3.1. The inviscid model, i.e., γ = 0. The condition for the linearised system (6) with γ = 0 to be hyperbolic, i.e. with perturbations that decay back to stability, is Proof. We examine the hyperbolicity of the above system by looking at the eigenvalues λ of the matrix A x (ρ s , θ s ). Neglecting the subscript s of (ρ s , θ s ), the equation where the discriminant Ifṽ (ρ) ≥ 0, in other words, if ρv(ρ) is an nondecreasing function, the system (6) with γ = 0 possesses two real eigenvalues and is hence always hyperbolic. Otherwise we obtain two real eigenvalues if Remark 1. An analogous relation can be derived for a perturbation in only the ydirection, where the right-hand side of the inequality (9) is replaced by cot 2 θ. The conclusion here is that the result above is generic, i.e. independent of the direction of perturbation, if we define θ as the angle between the direction of perturbation and the direction of Ω.
Proof. We apply a Fourier transform from position variable x to wave number ξ to (ρ σ , θ σ ): The system of (ρ σ ,θ σ ) in a matrix form is where The stability of the system (6) is equivalent to requiring that Imλ, the imaginary part of the eigenvalues of A ξ , and ξ have the opposite signs, so that the decay constant of the perturbation is negative. Note that the eigenvalues λ can be written as where the discriminant ∆ is and √ ∆ denotes the square root of the complex number ∆. Next we examine the sign of Imλ by comparing the values |Im √ ∆| and |γξ|.
where Sign(·) is the sign function for any real number. We introduce several notations: and √ ∆ = α + iβ. Using the expressions for ∆ and √ ∆, we obtain the equalities for the real and imaginary parts of ∆: There exist three cases: (i) β = 0. Then √ ∆ = α. Hence the system possesses two real eigenvalues and is alway stable.
To summarize, ifṽ (ρ) ≥ 0, ±|β| − |b| ≤ 0 for any ξ and Imλ has a different sign from ξ. Ifṽ (ρ) < 0, there is always one eigenvalue λ that Imλ has the same sign as ξ and sustains the unstable eigenmode, and hence the overall system will be unstable.
Remark 2. Usually, diffusive models are stable for large values of ξ because diffusion becomes stronger. Note that this is not the case here. Indeed, when the system is unstable, i.e.,ṽ (ρ s ) < 0, the unstable mode associated with the eigenvalue λ (where ξImλ > 0) grows with an exponential rate ξImλ. Moreover, In order to prove the above limit as ξ → ∞, we study the Taylor expansion of √ ∆ around 1 ξ . We employ the notations introduced in the proof of Theorem 3.2 and write ∆ in the polar representation: where θ ∆ = arctan 2aγξ a 2 + e − γ 2 ξ 2 + π.

Re
We are only interested in the imaginary part as ξ → ∞. The Taylor expansion at 1 ξ gives For the unstable mode, we have In addition, Eq. (12) tells us that in the limit we already took, ξImλ is increasing with respect to θ s for ρ s fixed. As ξ → 0, note that Im √ ∆ is approaching zero as well, which indicates a slow growth rate of the unstable modes for the long wavelength. 4. Numerical results of the microscopic and macroscopic models. The particle model is solved using the circle method that updates the angle of ω i at each discrete time step; the reader can refer to [20] for details. The challenge in the numerical resolution of the SOH model is caused by the geometric constraint |Ω| = 1 and hence the non-conservative nature of the system (1). The splitting scheme proposed in [20] solves the relaxation problem of (1) in the sense that the norm of Ω is initially not restricted to be 1, but then takes the zero limit of an expansion parameter η which realizes the geometric constraint on Ω. We will extend this idea to the SOH model, starting with the following proposition of the relaxation model. Proposition 1. Let η be a scalar parameter and (ρ η , Ω η ) the solutions to the relaxation model Then (ρ η , Ω η ) converges to the solutions of the SOH model (1) as η goes to zero.
Proof. The main idea is based on the fact that the right-hand side of (13b) is parallel to Ω η . Then the proof is analogous to the one in [20].
The numerical method, the so-called splitting scheme, is based on the work in [20] and details are provided in Appendix B.
We test the particle model and the SOH model on a rectangular domain [0, L x ] × [0, L y ] with L x = L y = 10 and impose periodic boundary conditions in both directions. We choose a time step ∆t = 0.001 to discretize time. For the function v(ρ), we consider a monotonically decaying power law with exponent α inspired by experimental and numerical results, and with an lower density threshold ρ * : below ρ * , the velocity is essentially constant, and then it rapidly decreases with power law α, v(ρ) = β ρ ρ * + 1 −α with three positive parameters ρ * , α, β.
Note thatṽ The parameters (ρ * , α, β) will be carefully chosen to place the system in the different regimes of stability derived in Section 3, in order to validate the corresponding stability results on the SOH model.
and d = 0.1, which are computed using the formulas derived in Appendix A.2 with the parameters of the particle model: N = 10 5 , ν = 100, D = 10, R 1 = R 2 = 0.1, so that we can perform the comparison between the two types of models. Indeed, this choice of parameters would correspond, after the scaling of Appendix A.1, to a value of ε = 0.01. The three parameters for the velocity function are chosen as (ρ * , α, β) = (0.12, 10, 1) which results in a stable model. This choice is only made for demonstration purposes.
The domain is uniformly partitioned using two integers N x and N y which represent the mesh size for the numerical integration in two dimensions. We apply the splitting method to the SOH model starting with the same initial condition and iteratively calculate the numerical errors with gradually decreasing mesh spacing (∆x = Lx Nx , ∆y = Ly Ny ). Fig. 1 is the profile of the errors as a function of mesh spacing ∆x in log scale and indicates that as expected the method is accurate to the first order.  Figure 1. The accuracy test at t = 1 shows that the splitting scheme is of first order. The initial data is given by ρ 0 = ρ s (1 + 0.1 sin(πx)), θ 0 = θ s (1 + 0.1 sin(πx)) with (ρ s , θ s ) = (0.01, π 4 ) and the mesh sizes are iteratively N x = N y = 32, 64, 128, 256.
We show that the SOH model agrees with the particle model using the following example. We construct an initial configuration with ρ 0 (x, y) ≡ 0.01 and a velocity field given by the Taylor-Green type function Ω 0 (x, y) = (sin π 5 x cos π 5 y , − cos π 5 x sin π 5 y ) T . Fig. 2(a) shows the contours of the density ρ and the velocity field Ω at t = 0.5, which are the average of 40 simulations using the particle model with N = 10 5 , ν = 100, D = 10, R 1 = R 2 = 0.1. This ensemble average is necessary due to the stochastic nature of both the initial particle positions and the angular dynamics. Fig. 2(b) is the numerical solution produced by the SOH model with (14) and d = 0.1. As explained at the the beginning of Section 4.1, the parameters for both the Vicsek and SOH models are consistent. The agreement is fairly good.

4.2.
Numerical results for the stability of the models. We vary the three parameters in the function v(ρ) to demonstrate the stability of the models around uniform steady states. The model parameters are given by (14). We choose an initial condition with a sinusoidal perturbation along x in both density and angle, in phase with each other, ρ 0 = ρ s (1 + σ sin(πx)), θ 0 = θ s (1 + σ sin(πx)) with (ρ s , θ s ) = (0.01, π 4 ). The model stability will be measured by the Root Mean Square Fluctuation (RMSF) of (ρ, θ), i.e. the L 2 norm of ρ − ρ s and θ − θ s : (a) t = 0.5: the particle model.   (9) does not hold. Therefore the resulting model is unstable and RMSF of (ρ, θ) must grow in time, which is demonstrated in Fig. 3 with RMSF in log-scale with respect to time t. An exponential growth of the perturbation, like we expect, would translate to a straight line in these graphs.     Let ρ * = 0.005 where we haveṽ (ρ) < 0 at the initial time. Hence the numerical solutions are expected to evolve away from (ρ s , θ s ), indicating the instablility of the model. Fig. 5 shows the RMSF of the numerical solutions (ρ, θ) for the time interval t ∈ [0, 5]. One can observe that they grow significantly; compare with the numerical solutions in Fig. 4 whose RMSF decreases to zero.

4.3.
Growth rate of the instability. The modal analysis of the hydrodynamic system (10) shows that the magnitude ofρ σ possesses an exponential growth rate ξImλ. In this section, we will examine the numerical solutions provided by the particle and SOH models and compare the numerical growth rate of the perturbation to the modal analysis result of the viscous system (10). For the sake of convenience, we perform a Discrete Fourier Transform (details can be found in Appendix C).
The initial configuration is taken as: where σ = 0.01, and ρ s = 0.01 is fixed. ρ σ takes the form and θ σ is given in a similar way. Here a 1 (ξ), a 2 (ξ) in ρ σ and θ σ are different random numbers generated from a uniform distribution on the interval [0, 1]. The parameters of the SOH model are fixed as c 1 = 0.975, c 2 = 0.925, d = 0.05 and γ = 0.12188. They correspond to the particle parameters where N = 10 5 , ν = 100, D = 5, and R 1 = R 2 = 0.1. The three parameters for v(ρ) are chosen as (ρ * , α, β) = (0.005, 2, 5) in which case the viscous model is unstable. The steady state for the density is chosen as ρ s = 0.01. We vary the steady state orientation θ s in the interval [0, π 2 ] and plot the growth rate ξImλ as a function ξ ∈ [0, 6] for both the linearized and SOH models in Fig. 6. Fig. 6(a) is computed using the fomula (11). To obtain Fig. 6(b), we proceed for each θ s in the following way. We compute the numerical solutions of the SOH model up to time t = 1, and perform the Discrete Fourier transform on the perturbed part ρ σ = ρ − ρ s to getρ σ (ξ, t). With different initial data, i.e. different a 1 and a 2 , we collect N sam samples of ρ σ (ξ, t) and apply a simple linear regression to the averaged quantities 1 N sam Nsamρ σ (ξ, t) ρ σ (ξ, 0) with respect to time t. This will generate the growth rate for each ξ; Fig. 6(b) shows the results for ξ = 0, 1, . . . , 6 with N sam = 100. The motivation here is as follows: by choosing random coefficients for the different modes we generate a set of initial conditions that have statistically equal weight for each mode. The growth rate is interpreted as the slope of the function t →ρ σ (ξ, t) in log scales. The similarity of these two sets of contours is the increase of the growth rate with respect to both ξ and θ s . The difference is that the fully nonlinear SOH model rapidly develops much larger growth rates compared to the linearised solution. There are also fluctuations at certain θ s and ξ likely due to finite size effects. Fig. 7 shows the growth rate of the perturbation ρ σ given by the particle model. The number of particles is N = 10 5 . The parameters are given as ν = 100, D = 5, R 1 = R 2 = 0.1. And they match the parameters for the SOH model in Fig.  6. The three parameters for v(ρ) are chosen as (ρ * , α, β) = (0.005, 2, 5) to match the SOH model. For each θ s , the growth rate is computed using the average of 10 simulations, in order to reduce the effects of noise. Although the contours do not exhibit the monotonic behaviour of the growth rate of the SOH model with respect to the steady state angle θ s and the eigenmode ξ due to the nonlinearity and stochastic effects, one can observe the stronger instability for larger θ s and ξ.    . Growth rate of the perturbation ρ σ given by the particle model with N = 10 5 . The parameters are ν = 100, D = 5, R 1 = R 2 = 0.1. The three parameters for v(ρ) are chosen as (ρ * , α, β) = (0.005, 2, 5). For each θ s , the growth rate is computed using the average of 10 simulations, in order to reduce the effects of noise.

Conclusion.
We have studied a Vicsek model where the velocity depends on the local density and then derived the corresponding SOH model. At the hydrodynamic scale, we analyse the stability of the two-dimensional inviscid and viscous models around their steady states. In summary, the stability of the SOH models is determined by the behaviour of the mass flux ρv(ρ). The theoretical results are illustrated by the numerical simulations with different choices of the velocity function. In general, we find good agreement between the theoretical prediction of the onset of instability and the numerical results. In the unstable regime, while our numerical results are qualitatively compatible with the predictions, the strong non-linearities present in the SOH model and especially the particle model quickly dominate the response. The SOH model we have developed here is a useful model to describe the semen flow in the experiments designed in [7]. These experiments record the correlation between the averaged velocity and the density of the sperm cells, and the velocity as a function of the density was fitted using the SOH model. Further study of vortices observed in the collective behaviour of the semen flow is under way.
Appendix A. Derivation of the SOH model with nonconstant velocity. The SOH model can be explicitly coarse-grained from the particle model introduced in Eq. (2). For completeness, we provide the main steps here and readers are referred to [12] for details.
A.1. The mean field model. We consider the limit of the system when N → ∞.
Sending N → ∞ and scaling out to obtain dimensionless parameters, the formal mean-field system for the probability distribution function f (x, ω, t) on R n ×S n−1 × (0, ∞) is given by where ∆ ω denotes the Laplace-Beltrami operator on the sphere and Here we assume that K 1 and K 2 only depend on the distance between particles, characterized by the dimensionless parametersȒ 1 andȒ 2 . We also assume that both of the ranges of the interaction kernels K 1 and K 2 are small. Let ε be a small positive number. We will perform the explicit hydrodynamic limit by introducing rescaled parameters that individually tend to zero in the limit → 0. ThenȒ 1 = √ εR 1 ,Ȓ 2 = √ εR 2 . We also assume that the alignment strength (a.k.a the social forces), and the diffusion coefficient are large, but of similar magnitude i.e.,ν For simplicity, we drop the hats and have the following result: Lemma A. 1. The density f ε (x, ω, t) satisfies the following equation: where , We drop the higher order term of ε, O(ε 2 ) and define the collisional operator Q(f ε ) by The rescaled system can be written as where A.2. The hydrodynamics model. This section derives the hydrodynamic model for the local density ρ f and the local mean orientation Ω f which will be valid at the macroscopic scale by taking the limit of the system (19)-(22) as ε → 0. We first introduce the von Mises-Fisher (VMF) probability distribution M Ω (ω) on S n−1 associated to a given Ω ∈ S n−1 : The main result in this section is the following theorem: Theorem A.2. Let f ε be the solution of (19)- (22). Assume that there exists f such that pointwise and the limit holds for its derivatives. Then there exist ρ(x, t) and Ω(x, t) such that and (ρ, Ω) are the solutions of where Here · M Ω denotes the integration with the weight function M Ω with respect to θ on the domain [0, π].
Proof. The proof is divided into three steps: (i) the determination of the equilibria; (ii) the Generalized Collision Invariants; (iii) the hydrodynamic limit. The three subsections below give a sketch of the proof.
Step 1. The equilibrium states, i.e., the null space of Q.
We will prove that the set E consists of the VMF distribution. More precisely, we have the following result: Lemma A.4. Under certain regularity assumptions, the set of the equilibria E is given by E = {f (ω) = ρM Ω (ω) for arbitrary ρ ≥ 0}.
Proof. Please refer to [12] for the proof of Lemma A.4.
Definition A.5. A collision invariant (CI) is a function ψ(ω) such that for any function f (ω) ≥ 0 with sufficient regularity we have We denote by C the set of CI, the set C is a vector space.
Due to the lack of physical conservation laws except for the total mass, the set of CI is not large enough to allow us to derive the evolution of the macroscopic quantities ρ and Ω. To overcome this difficulty, a weaker concept of collision invariant, the so-called "Generalized Collision Invariant" (GCI) has been introduced in [12]. We define the collision operator Q(Ω, f ) such that for a given vector Ω ∈ S n−1 , we have Notice that Q(f ) = Q(Ω f , f ).
where the numerical flux F i+ 1 2 ,j is defined as G i,j+ 1 2 is defined in the similar manner. Here P 2 ( ∂F ∂Q ) is a second degree polynomial of a matrix at the intermediate state of (Q i,j , ∂ x Q i,j ) and (Q i+1,j , ∂ x Q i+1,j ); see [13] for more details.
Appendix C. The discrete Fourier transform for the viscous system (10). Let N x × N y denote the mesh size over the domain [0, L x ] × [0, L y ] and (∆x, ∆y) the uniform mesh spacing. Each nonoverlapping computational cell is centered at (x j , y k ) = ( j − 1 2 ∆x, k − 1 2 ∆y) for 1 ≤ j ≤ N x , 1 ≤ k ≤ N y . We study the spatial variable x only and apply the Discrete Fourier Transform on (ρ σ (x j , t), θ σ (x j , t)): It follows that where the matrix where √ ∆ denote the square root of the discriminant, the complex number ∆: ∆ = c 1ṽ (ρ s ) − c 2ṽ (ρ s ) ρ s cos θ s + i 2πξγ L x 2 + 4c 1 dṽ (ρ s )ṽ (ρ s ) ρ s sin 2 θ s .