Coupled Self-Organized Hydrodynamics and Stokes models for suspensions of active particles

We derive macroscopic dynamics for self-propelled particles in a fluid. The starting point is a coupled Vicsek-Stokes system. The Vicsek model describes self-propelled agents interacting through alignment. It provides a phenomenological description of hydrodynamic interactions between agents at high density. Stokes equations describe a low Reynolds number fluid. These two dynamics are coupled by the interaction between the agents and the fluid. The fluid contributes to rotating the particles through Jeffery's equation. Particle self-propulsion induces a force dipole on the fluid. After coarse-graining we obtain a coupled Self-Organised Hydrodynamics (SOH)-Stokes system. We perform a linear stability analysis for this system which shows that both pullers and pushers have unstable modes. We conclude by providing extensions of the Vicsek-Stokes model including short-distance repulsion, finite particle inertia and finite Reynolds number fluid regime.


Introduction
Self-organised motion is ubiquitous in nature. It corresponds to the formation of large-scale coherent structures that emerge from the many-interactions between individuals without leader. Well-known examples are bird flocks, fish schools or insect swarms. However, self-organisation also takes place at the microscopic level, for example in bacterial suspensions and sperm dynamics (see e.g. Refs. [7,40] and the reviews [19,29,32]). In these cases, the environment, typically a viscous fluid, plays a key role in the dynamics.
In this paper we investigate self-organised motion of self-propelled particles (which we will refer to as 'swimmers') in a viscous fluid. The main difficulty in studying these systems comes from the complex mechanical interplay between the swimmers and the fluid. Particularly, highly non-linear interactions occur between neighbouring swimmers through the perturbations that their motions create in the surrounding fluid. While these interactions may be treated through far-field expansions in dilute suspensions [23], they require a much more complex treatment when the density of swimmers is high. Here we assume that, as a result of these swimmer-swimmer interactions, the swimmers align their direction of motion. In view of this, we adopt the Vicsek model for self-propelled particles undergoing local alignment to account for these swimmer-swimmer interactions in a phenomenological way. We then couple this model with the Stokes equation for the surrounding viscous fluid by taking into account the interactions between the swimmers and the fluid. The main goal is the derivation of macroscopic equations for this coupled system in terms of the time-evolution of the velocity of the fluid, on the one hand, and the swimmers' density and mean direction of motion, on the other hand.
The coupling terms considered here coincide with the ones in the kinetic Doi-Saintillan-Shelley model, which models active and passive rod-like dilute suspensions, [6,34,35]. This kinetic equation extends the Doi model [17,18] for liquid crystals (corresponding to passive rod-like or ellipsoidal particle suspensions) 0123456789().: V,-vol (space-independent) swimmer density and mean orientation fields as well as uniform fluid velocity and pressure. We linearise the SOH-Stokes system around these stationary solutions, meaning that we consider small perturbation of a uniform state. Note that the SOH model describes aligned states as the swimmer distribution is given by a von Mises distribution. So this analysis gives access to the stability of suspensions in their aligned state only. The investigation of the stability of the isotropic state is deferred to future work.
Here, we consider the two main types of swimmers present in nature: pullers and pushers. Pullers swim by using 'arms' to pull on the fluid (such as e.g., the green alga Chlamydomonas), see Fig. 1. Pushers typically have a tail with which they push the fluid behind them (such as e.g., spermatozoa), see Fig. 3. The stability analysis reveals that both pusher and puller types of swimmers have unstable modes. This is consistent with previous studies [34] which showed that both pushers and pullers are unstable to perturbations of an aligned state. Note that in [34] nematic interactions were considered, while we deal with polar interactions. Polar alignment corresponds to bringing the two swimming velocities parallel with the same orientation while nematic alignment brings them parallel but with either same or opposite orientation according to what angular change is the smallest. Steric interaction corresponds to the interaction of particles occupying finite volumes through volume exclusion. See an illustration of these various interactions in Fig. 2. Additionally, the aligned state in [34] is a Dirac delta in the orientation while ours is a von Mises distribution; we show that instability happens for all values of the angular dispersion around the alignment direction. We also notice that pullers are stable if they are slender rod particles. For both pushers and pullers, the instability only prevails for small |k| modes (or large wavelength). The largest growth rate takes place in the limit when the mode k → 0, which means that patterns induced by the instability will have roughly the same size as the system. Alignment interaction is not sufficient to prevent the appearance of large-concentration clusters in general [8]. So, in cases where such clusters are not observed, it is likely that short-range repulsion effects take place in addition to alignment. Following Ref. [10], we will investigate how both the micro and macroscopic models can be extended through the introduction of a short-range repulsion force. Besides, when the particle mass and size are larger, for example for fish, it is not legitimate to neglect the particle inertia and the fluid Reynolds number any longer. Therefore, we will show how to extend the micro and macroscopic models to include such finite size effects.
The document is structured as follows. In the next section we present the individual based model for the Vicsek-Stokes coupling and discuss the main result corresponding to its hydrodynamic limit. In Sect. 3 we present the mean-field limit, the scaling considered, and the Generalised Collision Invariant concept. In Sect. 4 we prove the main result. Sect. 5 shows the stability analysis. In Sect. 6 we extend the model to account for short-range repulsion and finite inertia and Reynolds number. Finally, we conclude in Sect. 7 discussing some perspectives on this problem.

The Vicsek-Stokes Coupled Dynamics
The dynamics of the viscous fluid follow Stokes equations. We couple these two models by incorporating the interaction mechanisms between the agents and the fluid. The dynamics of N agents are given by the evolution of (X i (t), ω i (t)) i∈{1,...,N } as a function of time t ≥ 0, where X i (t) ∈ R 3 is the position of the i-th agent and ω i (t) ∈ S 2 (the 2-dimensional sphere) is a unitary vector giving its direction of motion. We denote by v = v(x, t) ∈ R 3 the fluid velocity at position x ∈ R 3 at time t and p(x, t) ∈ R its pressure. Here we assume that the fluid density remains constant. In Sect. 6.2 (see Remark 6.2) we derive the following Vicsek-Stokes coupled dynamics, where all the quantities are dimensionless and where the stochastic differential equation (2.1b) must be understood in the Stratonovich sense, where the unknowns In this system a, ν, D, λ, R and b are constants. The symbol '⊗' denotes the tensorial product and 'Id' the 3×3 identity matrix. The symbol P ω ⊥ i = Id − ω i ⊗ ω i gives the orthonormal projection operator onto the sphere S 2 at ω i ; the '•' symbol following it indicates that the Stochastic Differential Equation (2.1b) has to be understood in the Stratonovich sense. The terms (B i t ) t≥0 , i = 1, . . . , N are independent Brownian motions in R 3 . The terms S, A are matrices that will be defined later. The operators Δ x , ∇ x , ∇ x · indicate the Laplacian, the gradient and the divergence in R 3 , respectively. The symbol δ X is the delta distribution in R 3 at X ∈ R 3 . Finally, K = K(r) ≥ 0, r ≥ 0, is a given sensing function.
We explain first the meaning of the equations without the coupling terms. Equations (2.1a)-(2.1c) without the terms involving the velocity of the fluid v correspond to the Vicsek model: each agent i moves at a constant speed a > 0 in the direction ω i while trying to adopt the average direction of motion of its neighbours. This averaged direction is given byω i in Eq. (2.1c). The positive kernel K weights the influence of the neighbouring agents and the constant R > 0 gives the typical interaction range between agents. The intensity of alignment is given by ν > 0. While trying to align, agents make errors. This is modelled via a noise term √ 2DdB i t , where D > 0 is the standard deviation of this random motion per unit of time. The presence of projection operator P ω ⊥ i ensures that |ω i (t)| = 1 for all times (since the stochastic differential equations (2.1a)-(2.1b) are interpreted in the Stratonovich sense, see Ref. [25]).
Equations (2.1d) and (2.1e) are the Stokes equation for the velocity of the fluid v with a time-dependent force term at the right-hand side of Eq. (2.1d) whose meaning will be explained in the following section. The hydrostatic pressure of the fluid p = p(x, t) is the Lagrange multiplier of the incompressibility constraint (2.1e).
The coupling terms.

Effect of the fluid on the agents.
(i) Effect on the agents' velocity In the limit of zero particle inertia there is an instantaneous relaxation of the passive part of the particles' velocity (i.e., the particle velocity minus the self-propulsion velocity) to the velocity of the fluid (see Appendix 6.2). As a consequence, the term u i in Eq. (2.1a), giving the total velocity of agent i, is the sum of the fluid velocity v and the agent's self-propelled velocity aω i . (ii) Effect on the agents' orientation This is expressed by the term (λS where the matrices A and S are the antisymmetric and symmetric parts of the linear flow ∇ , respectively: where the exponent 'T ' indicates the transpose of the matrix. This term encompasses Jeffery's equation, which describes the effect of a viscous fluid on the orientation of a spheroidal passive particle. In a spatially homogeneous flow where ∇ x v is constant, these equations give the motion of the principal axis of spheroidal particles, as follows: where ∇ ω is the gradient on the sphere S 2 ; ∇ x × denotes the curl; and the symbols '·', '×' denote the inner product and the cross product in R 3 , respectively. The matrix S describes straining forces in the fluid which forces a passive particle to orient in a preferred direction called 'local extensional axis', given by the eigenvector of maximal eigenvalue of S. The matrix A describes shear effects in the fluid that have the effect of rotating the suspended particle around an axis parallel to the vorticity ∇ x ×v. The parameter λ ∈ [−1, 1] is a shape parameter: for a spheroidal particle with aspect ratio χ, we have λ = (χ 2 − 1)/(χ 2 + 1). The limit λ → 1 corresponds to a slender rod-like particle, the limit λ → −1 corresponds to a thin disk and the case λ = 0 corresponds to a sphere. For an explanation of Jeffery's equation see e.g Refs. [6,26]. 2. Effect of the agents on the fluid. We consider two forces produced by agents that act on the fluid: (i) Drag force exerted on the fluid by the motion of the agents The motion of the agents creates a drag force on the fluid. However, in the limit of zero particle inertia this force vanishes, as detailed in Sect. 6.2. So we do not take it into account here. There is a symmetric effect of the fluid on the particles which in this limit produces instantaneous relaxation of the passive part of the velocity (the particle velocity minus the self-propulsion velocity) to the fluid velocity hence justifying Eq. (2.1a), see above. (ii) The self-propulsion force The source term that appears on the right-hand side of the equation for v (2.1d) describes the influence of the self-propulsion force of the agents on the fluid. The term ∇ x δ Xi(t) denotes the gradient of the Dirac delta δ Xi(t) and it is defined in weak form for any vector test function ϕ by: where T, ϕ denotes the duality bracket between a distribution T and a test function ϕ.
Suppose that an agent swims by pushing with its tail with a force F in the opposite direction of motion. Then the head also exerts a force − F on the fluid. This type of swimmer is called a 'pusher'. If the centre of the swimmer is in location X i and it has a length , then the pushing force is applied at location X i − 2 ω i while the force of the head is at X i + 2 ω i , see Fig. 3. Since F and − F are applied at different points they do not cancel, this is referred as a 'force dipole' and in this case From a Taylor expansion for small , this leads to the right-hand side term in Eq. (2.1d), with b > 0 (after dividing by the fluid viscosity). We do not give here the details of this derivation but refer the reader to Refs. [6,34,35] and references therein. Another common way of swimming is by using the arms. Then the swimmer is called 'puller' and in this case b < 0 in Eq. (2.1d). For a drawing of a puller see Ref.
[6, Fig. 3]. 3. Other interactions and effects. The model presented here is a simplification of the actual dynamics.
Therefore, it could be enriched by taking into consideration other mechanical effects, like noise in the spatial variable (see the Doi-Saintillan-Shelley model [6]) and extra forces acting on the fluid due to the inextensibility of the particles (resistance to stretching and compression) when the particle sizes are not supposed infinitesimally small, see Ref. [20]. We will consider refined versions of this model in Sect. 6.1 where we add volume exclusion between agents, or in Sect. 6.2 where we include both fluid and particle inertia.
Remark 2.1. Jeffery's equations have been derived by Jeffery under the assumption of spheroidal particles (i.e. ellipsoids with a circular cross-section). In [28], these equations are generalized to a body of arbitrary shape. In this reference, it is also shown that Jeffery's equations are the leading order of an asymptotic expansion of the Navier-Stokes equation in powers of the body size and the next order of this expansion is derived. However, these terms do not have an analytical formula and their use in the forthcoming theory would represent a significant challenge which is left aside here.

Macroscopic Coupled Dynamics: The SOH-Stokes Model
In Sect. 4 (Theorem 4.1), from the Vicsek-Stokes dynamics (2.1a)-(2.1e), we derive the following macroscopic system , that we refer to as the 'Self-Organised Hydrodynamics-Stokes model' (SOH-Stokes). It gives the time-evolution of the spatial density of agents ρ = ρ(x, t), the mean direction of motion Ω = Ω(x, t), the velocity of the fluid v = v(x, t) and the pressure p = p(x, t):  [11,15,27] (the diffusive term is derived in Ref. [11]). The SOH system resembles a fluid dynamics equation, particularly, a compressible Navier-Stokes equation, with the differences that c 1 = c 2 and the 'velocity' Ω is constrained to be of norm one with the presence of the projection operator P Ω ⊥ in the pressure ∇ x ρ and the diffusion Δ x (ρΩ). This projection operator precisely ensures that |Ω| = 1 at all times (provided that |Ω| t=0 = 1), but, as a consequence, the equation is not conservative, meaning that the terms involving spatial derivatives cannot be written as the spatial divergence of a flux function. Eqs. (2.3c)-(2.3d) without the right-hand side in (2.3c) correspond to the Stokes equation.
The coupling terms. The terms involving the velocity of the fluid v in Eqs. (2.3a)-(2.3b) express the effect of the fluid on the particles. As expected, the continuity equation (2.3a) for the density ρ has a velocity U which is the sum of the local average self-propulsion velocity (cc 1 Ω) and the velocity of the fluid v. Also the convective velocity in Eq. (2.3b) for Ω is a weighted sum of Ω and the velocity of the fluid. The last term in Eq. (2.3b) reflects how Jeffery's equation on individual swimmers (2.2) translate to the population level. The terms resulting from Jeffery's equation describe the propensity of Ω to align in the direction of the so-called local extensional axis, given by the eigenvector of maximal eigenvalue of S, as well as to rotate about the vorticity axis, parallel to ∇ x × v. Notice, though, that the shape parameterλ = λ. Numerically, (see Fig. 4a) we observe that λ 0 ∈ [0, 1]. This implies thatλ ∈ [−1, 1];λ has the same sign as λ; and |λ| ≤ |λ|. Therefore, Jeffery's equation for the agents' individual orientations is coarse-grained into another Jeffery's equation for the local mean orientation, but the 'mean particle shape' associated withλ is different from the individuals' shapes associated with λ. Particularly, when κ → 0 (large noise regime), we getλ = 0, which corresponds to the shape of a sphere, and when κ → ∞ (low noise regime), we getλ = λ, and we recover the original shape parameter.
Finally, the right-hand side in Eq. (2.3c) gives the influence of the agents on the fluid. It involves the divergence of the deviatoric stress tensor Q(Ω), i.e., the contribution of the swimmers to the extra-stress, which provides its non-Newtonian character to the fluid. This term results from the coarse-graining of the right hand side of Eq. (2.1d). Numerically we observe (Fig. 4b) that c 4 > 0, which shows that the coarse-graining of a population of pushers preserves the 'pusher' behaviour, as it should.
The reader is referred to Sect. 6 for some extensions of this model.

The Mean-Field Limit Equation
As explained in the introduction, the derivation of the macroscopic equations is carried out with an intermediate step: the kinetic or mean-field equations. The mean-field limit of System (2.1a)-(2.1e) provides the time-evolution of the distribution function f = f (x, ω, t) in space and orientation of a typical agent. From the equation on f , we will derive the macroscopic equations in Sect. 4. For the case of the Vicsek model alone, a rigorous proof of the mean-field limit has been obtained in Ref. [3] when there is no normalisation ofω i in Eq. (2.1c), i.e., whenω i = J i . Following the proof in [3] formally we have the:

Proposition 3.1 ((Formal) Mean-field limit). Consider the empirical distribution associated to the dynamics of the agents in Eqs
where δ xi(t) (x) and δ ωi(t) (ω) denote the Dirac delta at x i (t) and ω i (t) on R 3 and S 2 , respectively. Assume that f N converges weakly to f = f (x, ω, t) as the number of agents N → ∞. Then, the limit f satisfies the following system: where ∇ ω · and Δ ω stand for the divergence and the laplacian in S 2 , respectively; and where

Scaling and Expansion
We scale the alignment intensity and the variance of the noise by setting ν =ν/ε, D =D/ε, whereν, D are given fixed quantities. Considering the classical Vicsek model (without the coupling terms), this rescaling corresponds to i.e., it corresponds to a time-rescaling t = t/ε which, as ε → 0, gives the long-time dynamics for ω i . In order words, with this rescaling we express the fact that the self-propulsion velocity of the agents ω is a fast-varying variable while the velocity of the fluid v is a slow-varying variable. Notice, however, the invariance of the quotient that we denote by κ. We also scale the radius of influence R in Eq. (3.4) by setting R = √ εR, which is the rescaling considered in Ref. [11]. This rescaling expresses that the interactions between agents become localized in space as ε → 0. After rescaling the kinetic equation (3.1) in this way, we obtain (after skipping the tildes): We simplify this system by considering the following expansion: and Proof. The result is a direct consequence of the Taylor expansion for after performing the change of variables z = (x − y)/( √ εR), which gives: , Thanks to the previous Lemma 3.2, we can rewrite the rescaled system (3.6) as follows: with

Equilibria and Generalised Collision Invariants
In Ref. [15] the authors studied the operator Q given in Eq. (3.14). They proved that it can be recast into a Fokker-Planck form: where the density on the sphere is the so-called von Mises distribution (Z is a normalizing constant). The equilibria of Q as a function of ω are given by the set of functions Moreover, in Ref. [15] it is proven that for showing the consistency relationship Details can be found in Eq. (4.1). Collision invariants are fundamental in the derivation of macroscopic equations. They are defined as the scalar functions ψ such that In the present case, ψ = constant clearly satisfies this relation. This is a consequence of the conservation of mass during the interactions between agents. It can be shown that there are no other conserved quantities. This implies, particularly, that the dimension of the space of collision invariants is smaller than the dimension of the kernel Q in (3.15), which is 3-dimensional. Classical methods require the dimension of the two spaces to be the same in order to derive a full system of macroscopic equations. The collision invariant corresponding to the constants will allow us to derive the equation for the spatial density ρ = fdω (as we will see in the next section), but it will not be enough to determine the equation for the mean orientation Ω. To sort out this problem, the authors in Ref. [15] introduce the concept of Generalised Collision Invariant (GCI) defined as follows: where .
Notice that with this definition It has been proven that the GCI has the following properties: For a given function f : S 2 → R, we consider the associated Ω f given by and consider Then

Macroscopic Limit: The SOH-Stokes System
In this section we investigate the hydrodynamic limit as ε → 0 for the system (3.10)-(3.13). We will use the following change of variables: for Ω ∈ S 2 fixed, we decompose any given vector ω ∈ S 2 uniquely as We take the convention S 2 dω = S dw = 1. One can check that (see Ref. [21, Ap. A2]) for any measurable function a(ω) =ā(θ, w): and S w dw = 0, and We will also use the notations: where h is the function appearing in Eq. (3.19).
Theorem 4.1 ((Formal) macroscopic limit). Consider the rescaled system (3.10)- (3.13). When ε → 0, it holds (formally) that where ρ = ρ(x, t) ≥ 0 and Ω = Ω(x, t) ∈ S 2 are the limits of the local density ρ ε = S 2 f ε dω and the local mean orientation Ω f ε in Eq. (3.9), respectively. Moreover, if the convergence is strong enough and Ω, ρ, v and p are smooth enough, they satisfy the coupled system (2.3a)-(2.3d) with explicit constants where we used the following notation: for any functions g, The constants a, b, κ = D/ν correspond to the ones in the individual based model (2.1a)-(2.1e) and the value of k 0 is given in Eq. (3.7).
Proof. Suppose that f ε converges to f as ε → 0. Then, from Eq. (3.10), f belongs to the kernel of Q, i.e., Q(f ) = 0. Therefore, f = ρM Ω by Eq. (3.15), with ρ = ρ(x, t) ≥ 0 and Ω = Ω(x, t) ∈ S 2 . We start by computing the equations for these two macroscopic quantities. We obtain the continuity equation (2.3a) for ρ by integrating the kinetic equation (3.10) with respect to ω; dividing by ε; taking the limit ε → 0; and using the consistency relationship in Eq. (3.16). Notice that the integral of the right hand side of the kinetic equation (3.10) vanishes since ψ = 1 is a collision invariant in Eq. (3.17) , i.e., We compute next Eq. (2.3b) for the mean direction of the agents Ω. We multiply the kinetic equation (3.10) by 1 is the Generalised Collision Invariant given by Proposition 3.4, and integrate with respect to ω: Notice that the term involving Q vanishes thanks to Eq. (3.20). Taking the limit ε → 0 on the previous expression, we obtain: or, equivalently, To compute this last expression we decompose X into X = X 1 + X 2 + X 3 + X 4 for Notice that in X 3 we used the incompressibility condition The term P Ω ⊥ X 1 has been computed in Ref. [15]: with the notations in Eq. (4.4). The term P Ω ⊥ X 2 has been studied in Ref. [11]. One observes that in the limit j ρMΩ = c 1 ρΩ using Eq. (3.16) and, therefore, for any vectors A, B ∈ R 3 (see Ref. [21]). With this expression we have that where in the second equality we used the change of variable (4.1), as well as, (this formula is a consequence of Eqs. (4.2)-(4.3)); in the third equality, the odd integrands in w vanish; and in the last equality we used Eq. (4.3). Now, the terms X 3 and X 4 correspond to the coupling terms. Firstly, for X 3 we have that where in the second equality the term (v · ∇ x )ρ vanishes since and in the last equality we used that where we used that ∇ ω M Ω = κP ω ⊥ ΩM Ω and that ∇ ω · (P ω ⊥ Bω) = B : (Id − 3ω ⊗ ω) for any matrix B independent of ω. The notation B : C indicates the contractions of the two matrices B = (B) ij , . In this way we can decompose X 4 into X 4 = X 41 + X 42 + X 43 with To compute the term X 41 , notice that, if C is an antisymmetric matrix, then C : (Id − 3ω ⊗ ω) = 0 (since the second matrix is symmetric), therefore In the following computation, in the second equality we use the change of variables (4.1); in the third equality the odd terms in w vanish from the integral; in the fourth equality we use that S : (w ⊗ Ω + Ω ⊗ w) = 2w · SΩ (since S is symmetric); and the last equality is consequence of Eq. (4.3): For the term X 42 it is immediate to obtain proceeding analogously as in previous computations (remember Eq. (4.9)). The term X 43 is computed similarly as for X 41 : where in the last equality we substituted (B + B T )/2 = λS; and where Grouping terms we conclude: Finally, putting all the terms together, we obtain Dividing the previous expression by κC 0 we obtain Eq. (2.3b) for Ω with To conclude the theorem, we are left with computing the limit for Stokes equation (3.12). For this, we just need to compute the limit of the right hand side term, which in the limit ε → 0 corresponds to We compute next the value of the integral: Id + c 5 Id , P. Degond et al. JMFM where in the third equality we have disregarded the odd terms in w; and where A computation shows that c 5 = 0. This implies that Id .

Linearised Stability Analysis of the SOH-Stokes System
In this section, we investigate the linearised stability of the SOH-Stokes system (2.3). We linearize the SOH-Stokes system about constant (space-independent) functions ρ, Ω, v, p and study the stability of the resulting linear system. The main result of this section is that the SOH-Stokes model exhibits unstable modes for both pushers (b > 0) and pullers (b < 0). Since the SOH model describes aligned states (as the particle distribution function is non-isotropic, given by a von Mises distribution with non-zero parameter κ), this corresponds to analyzing the stability of the suspension near an aligned state. A previous analysis performed in [34] in the case of nematic interactions, Fig. 2 (see also [24]) has shown that both pushers and pullers are unstable to perturbations of an aligned state. We show that this instability still prevails for both pushers and pullers interacting though polar alignment. However, we show that pullers can be stable if they are slender rods (λ = 1). We will also see that the unstable modes for pushers and pullers are not the same. In the case of pullers, these are transverse modes (the perturbation to Ω is normal to the wave-vector) propagating along the unperturbed orientation vector Ω. For pushers, these are longitudinal modes propagating transversely to the unperturbed orientation vector Ω. The former have vanishing density perturbation while the latter have non-trivial density perturbation. For both pushers and pullers, the instability only develops at small values of |k| (i.e. for large wavelengths) and has maximal growth rate at k = 0. Therefore, we can expect that the typical spatial extension of the instability patterns will be set up by the system size.
Here we assume that a = 1 to simplify the analysis. Let be a uniform steady state for the SOH-Stokes system with |Ω 0 | = 1. We expand it with a small perturbation parameter τ : Dropping the higher order terms O(τ 2 ) and using (ρ, Ω, v, p) to represent the first order perturbation (rather than (ρ 1 , Ω 1 , v 1 , p 1 )) we obtain the linearised system: The first equation is consequence of |Ω 0 | = 1. To deduce the first term in the right hand side of Eq. (5.1c) we used that P Ω ⊥ 0 Ω 0 Δ x ρ = 0. Finally, to obtain the last term in Eq. (5.1d) we used that The main result of this section is the following: In this caseρ = 0,Ω is an arbitrary unit vector orthogonal to Ifb < 0 (puller case), the modes are unstable in the range 3)

4)
and corresponds to the limit of Im(α) when k → 0.

5)
with η = ±1 and (α, k) are linked by the following dispersion relation, In the particular case where k 0 = 0, the dispersion relation simplifies tõ The corresponding modes are stable (Im(α) ≤ 0) if and the perturbation is given bȳ Ifb > 0 (pusher case), the modes are unstable in the range

The supremum of Im(α) in this range is
and corresponds to the limit k → 0.
Remark 5.1. Notice that Case (B) when k 0 = 0 corresponds tov ⊥Ω, while Case (A) (b) corresponds tov Ω . This is the signature that these two cases correspond to different modes.

Case (A)(a):
This case corresponds to the simple propagation of a density perturbation along Ω 0 at speed with no perturbation of the orientation sinceΩ = 0.

Case (A)(b):
Notice thatλ − 1 ∈ [−2, 0], sinceλ ∈ [−1, 1]. Therefore, ifb > 0 (pushers), the mode is stable and ifb < 0 (pullers) the mode is unstable for small values of |k|. In this last case, the coupling with Stokes equation destabilizes the model given that in the SOH model alone all modes are stable, see Ref. [15]. Notice that the diffusion term helps to stabilize the modes by damping them, but since it involves a second order derivative, the damping is proportional to |k| 2 and is very small for small values of |k| but dominates for large values of |k|. Therefore, for large values of |k|, the diffusion damping is enough to compensate the instability due to the coupling with the Stokes equation, which is independent of |k|. This is why the model is stable for large values of |k|. However, for small values of |k|, the diffusion damping is not strong enough and the instability of the Stokes coupling is predominant. Consequently, small |k|-modes (large wavelength) are unstable. Moreover, the supremum of Im(α) corresponds to the limit k → 0, which means that the typical spatial extension of the fastest growing unstable mode will be of the size of the system.

Case (B):
We analyse the particular case k 0 = 0. We observe analogous phenomena as in Case (A)(b) but reversing the roles of pullers and pushers sinceλ + 1 ≥ 0: ifb < 0 (pullers), the constant solution is always stable but in the pusher case (b > 0), the coupling with Stokes equation destabilizes the mode. The supremum value of Im(α) also corresponds to the limit k → 0.
Remark 5.3. Figures 5, 6 and 7 provide a schematical explanation of the instability mechanisms. Figure 5 depicts the perturbation velocity field generated by pushers and pullers. Figures 6 and 7 provide a description of the instability mechanisms for pullers and pushers respectively.
Remark 5.4. The instability shown here is a classical case of long wave-length instability where the unstable modes are found only for a bounded range of k. The mode with largest growth rate (most negative imaginary value of α) will grow and shortly dominate but then the solution will leave the domain of validity of the linear approximation. Nonlinear terms will take control. A likely scenario is that the energy from the unstable long wave-length modes will be cascaded to the stable short wave-length ones by the nonlinearity and will be ultimately dissipated, leading to the stabilization of the system. However, in the absence of a rigorous analysis of the full nonlinear model, this scenario remains a conjecture. Note however that this situation is classical in fluid mechanics and gives rise to patterns such as the von Karman vortex streets. The situation is likely to be similar here and should give rise to the appearance of patterns that could be explored numerically. In order to take advantage of the stabilization by the short wave-length modes, numerical simulations should use a fine enough mesh so that the shortest wave-length resolved by the mesh (of the order of the mesh size) lies in the stability region. This indeed may lead to severe numerical constraints when the diffusion coefficient γ is small. For larger mesh sizes, numerical diffusion may provide a stabilization mechanism, but at the expense of a loss of accuracy. Finally, note that the stability analysis of the macroscopic equations is analytically solvable. It would not be the case for the particle model for which stability can only be explored numerically. This is an advantage of the macroscopic model over the particle model. P. Degond et al. JMFM Fig. 6. Puller unstable mode: a geometric configuration. The puller unstable mode is a transverse mode (Ω ⊥ k) propagating parallel to the unperturbed mean orientation Ω 0 . b Schematics of the instability mechanism. There is no density perturbation involved. The instability is due to a reinforcement of the mis-alignment between the swimmer mean orientation Ω = Ω 0 + τΩ, (τ 1) and the unperturbed orientation Ω 0 . This reinforcement results from the torque applied to a given swimmer by the velocity perturbation generated by the neighbouring swimmers ahead and behind it. This torque is materialised in the picture by the double arrows Fig. 7. Pusher unstable mode: a geometric configuration. The pusher unstable mode is a longitudinal mode. The perturbationΩ is parallel to the propagation direction (Ω k) and both are perpendicular to the unperturbed mean orientation Ω 0 . b Schematics of the instability mechanism. Due to the configuration of the perturbation velocity field that they generate, swimmers are attracted by regions of higher swimmer density, thereby amplifying density perturbations Proof of Theorem 5.1. Substituting the plane-wave solution (5.2) in the linearised system (5.1) we obtain (5.10e) where in (5.10b) we used that k ·v = 0 thanks to (5.10e); in (5.10c) we used that P Ω ⊥ 0Ω =Ω thanks to (5.10a), as well as, that ∇ x v = ik ⊗v and therefore and so Now, we look for the existence of a non-trivial solution of system (5.10). From the last two equations we deduce thatp We note that we can divide by |k| 2 since k = 0. Otherwise, if k = 0, then in the case α = 0 (which is the case of a non-trivial perturbation we are interested in), this implies thatρ = 0,Ω = 0, i.e., the perturbation is null. Next, we obtain an expression forv by decomposing it intov = P Ω ⊥ 0v + P Ω0v , since this will be useful in the sequel. Doing the inner product of Eq. (5.10d) with Ω 0 and using Eq. (5.10a), we obtain:v Projecting now Eq. (5.10d) on the orthogonal to Ω 0 we obtain We insert these expressions in (5.10c) to obtain: The condition for stability is Im(ω) ≤ 0, i.e., In this case one can check thatρ = 0,Ω is arbitrary withΩ, Ω 0 = 0,p = 0 and Moreover, forb < 0, it is straightforward to see that the range for which |k| is unstable is given by (5.3) and the supremum of Im(α) is attained at (5.4) in the limit k → 0.

Case (B)
Suppose that k ⊥ = 0. The coefficient on the right-hand side of Eq. (5.12) is written (thanks to (5.11)) as First we check that X = 0. Suppose thatρ = 0. From Eq. (5.10a), this implies thatk = 0. So X = 0 and we conclude thatΩ = 0. So the perturbation is null. Therefore it cannot be thatρ = 0, which implies Im(X) = 0, so that X = 0. Therefore,Ω = 0 and from (5.12), it should be given by Eq. (5.5). Then, from (5.12) again and the fact that k ⊥ = (|k| 2 − k 2 0 ) 1/2 , the dispersion relation is given by: and from Eq. (5.10b) we have the relation We check that − α + U 0 · k = 0 by contradiction. Suppose that − α + U 0 · k = 0, then, from the previous equation, we deduce thatk =Ω · k =Ω · (k ⊥ + k 0 Ω 0 ) = 0. From Eq. (5.5), we get thatΩ k ⊥ , which implies that |k ⊥ | = 0. This contradicts our assumption that k ⊥ = 0 and therefore, we conclude that, effectively, − α + U 0 · k = 0. So, multiplying Eq. (5.13) by −α + U 0 · k = 0 and using Eq. (5.14), we get the dispersion relation in Eq. (5.6). Now, to simplify the analysis we will restrict ourselves to the case where k ⊥ = k, i.e. k 0 = k · Ω 0 = 0. This implies, in particular, that U 0 · k = V 0 · k = v 0 · k, as well as, With these considerations one can simplify the dispersion relation (5.6) intõ Using the variable X = α − v 0 · k we can recastD(α, k) = 0 into: after multiplying by i. Now, changing variables X = iY , we have that Y solves Stability in this case means Im(X) < 0, i.e. Re(Y ) < 0. The polynomial P has real coefficients. There are two possibilities: • If P has positive discriminant, its two roots are real. In this case, to have stability we require them both to be negative, i.e. their product π has to be positive and their sum σ negative. The product is given by and their sum is So in this case the stability criteria corresponds to σ ≤ 0, which leads to Eq. (5.8). • If the polynomial P has negative discriminant, the two roots are complex conjugate. Their real part is half their sum σ. So again the stability criterion reduces to asking that σ is negative, and we are left with the same stability criterion (5.8) as before.
We suppose now thatb > 0, and we want to determine the supremum on the instability range |k| ∈ ⎡ ⎣ 0, ρb(λ + 1) 2γ and the corresponding value of |k|. This supremum corresponds to k max = argmax Re(Y ) = argmax Im(α). We have seen that in the case where the roots are real, they have the same sign. Therefore, any root is less than the sum σ(|k|) given by (5.15). The maximum value of σ(|k|) is at k = 0, i.e., One can easily check that σ(0) is a root of P for |k| = 0 (the other root being 0) and therefore one of the roots is maximal at |k| = 0.

Extensions of the Model
In the next two sections we consider two ways of extending the Vicsek-Stokes model: (i) adding shortrange repulsion to account for volume exclusion interactions and, (ii) considering the Vicsek-Navier-Stokes coupling to account for regimes with high Reynolds number and finite particle inertia. Further to this, one could consider also a combination of these two extensions. Stability analysis is not performed here for these two extensions. Indeed, while the methodology is a straightforward extension of that used for the Vicsek-Stokes model, the analytical resolution of the dispersion relation is far from obvious and requires further developments. It would be risky to make conjectures based on the results for the Vicsek-Stokes model as there are examples where a model and its singular limits do not enjoy the same stability properties.

Adding Short-Range Repulsion
The Vicsek-Stokes coupling (2.1a)-(2.1e) presented here can be extended towards different directions. Particularly, in regions where agents become highly packed, a repulsion force can be enforced between neighbouring particles to better account for volume exclusion (or steric interactions, see Fig. 2). This can be easily done following Ref. [10] where repulsion is introduced in the Vicsek model and coarse-grained into the Self-Organised Hydrodynamic model with Repulsion (SOHR). Particularly the individual based model corresponds to: with the same notations as for the system (2.1a)-(2.1e), where μ, ξ > 0 and the repulsive potential Φ is defined as where φ = φ(|x|) is a binary repulsion potential that only depends on the distance, and where r > 0 is the typical repulsion range. We assume that x → φ(|x|) is smooth, as well as, Degond et al. JMFM which implies, in particular, that φ(|x|) → 0 as |x| → ∞. The only differences between System (6.1)-(6.5) with the original Vicsek-Stokes system (2.1a)-(2.1e) are the addition of two new terms: the last term to the evolution of X i (t) in Eq. (6.1), which expresses the repulsion force, and the second term in the evolution of ω i (t) in Eq. (6.2), which is a relaxation term of ω i towards the force ∇ x Φ(X i (t), t). This terms models the fact that particles tend to actively align their directions of motion with the force.
The presence of these new terms modifies the coarse-grained equations. To begin with, the mean-field equations correspond to (following Sect. 3 and Ref. [10]): following the notations of Proposition 3.1 and where The difference with respect to the mean-field system in (3.1) is the extra term ξ∇ x Φ f (x, t) in the equation for f and the term −μ∇ x Φ f in the equation for the velocity u (f,v) .
To perform the macroscopic limit, we rescale the mean-field equations (6.6) analogously as in Sect. 3 adding the rescaling of r = εr forr > 0. Remember that the alignment interaction range R is rescaled as R = √ εR, therefore the alignment interaction range is larger than the repulsive range. Skipping the tildes, we obtain the rescaled system: , From here we obtain the macroscopic equations as ε → 0: Theorem 6.1 (Macroscopic equations with volume exclusion). Consider the rescaled system (6.7). When ε → 0, it holds (formally) that where ρ = ρ(x, t) ≥ 0 and Ω = Ω(x, t) ∈ S 2 are the limits of the local density ρ ε and the local mean orientation Ω f ε in Eqs. (3.2), (3.9), respectively. Moreover, if the convergence is strong enough and Ω, ρ, v and p are smooth enough, they satisfy the coupled system and where the constants c 1 , . . . , c 4 , k 0 , γ and λ are given by Eqs. The proof of this theorem is direct from the Proof of Theorem 4.1 and the results in Ref. [10].
Remark 6.1 (Discussion of the result.). The repulsive force intensity is given by the parameter μΦ 0 . Observe that when μΦ 0 = 0 we recover the SOH-Stokes system (2.3a)-(2.3d). Notice that the presence of the repulsion modifies the velocity U of ρ in (6.8) and the convective velocity V of Ω in (6.9) by adding a term −μΦ 0 ∇ x ρ. This term in (6.8) gives rise to a diffusion-type term for ρ of the form −μΦ 0 ∇ x · (ρ∇ x ρ), which resembles a porous-medium equation and that prevents the formation of high particle concentrations. In the case of the convective velocity of Ω, this term indicates the tendency of particles to change their orientation towards regions of lower concentration. The other important difference is the presence of a non-linear term in the pressure ∇ x p(ρ) for Ω in (6.9) which increases the pressure effects, due to the repulsion forces, when the concentrations become high.

The Individual Based Model.
In a finite Reynolds number regime, fluid dynamics is described by the Navier-Stokes equations rather than by the Stokes equation. In this section, we propose a Vicsek-Navier-Stokes coupling also assuming finite particle inertia and derive the coarse-grained equations. We will also see how the Vicsek-Stokes coupling in Eqs. (2.1a)-(2.1e) is obtained from this Vicsek-Navier-Stokes coupling by assuming a low Reynolds number regime and negligible particle inertia. We consider the following coupled system: Most of the terms have previously been explained for Eqs. (2.1a)-(2.1e) in Sect. 2. The term ρ 0 is the density of the fluid and σ > 0 its viscosity; η is a friction coefficient; m i is the mass of agent i; and F i is the force generating its acceleration. Notice that in the present case the individuals' velocity u i relaxes towards v(X i , t) + aω i (t), while in the Vicsek-Stokes coupling we considered directly the relaxed system. The influence of the force of particle i on the fluid is given by the term F i δ Xi(t) in Eq. (6.12f) (this is an application of Newton's third law of action and reaction). A sanity check of our model consists of ensuring that the momentum and the angular momentum are conserved by the dynamics, as expected in a closed system with no dissipation at the boundaries: Proposition 6.2. Suppose that in the system (6.12) the domain has no boundaries and the solution vanishes at large distances, then the total momentum and angular momentum are conserved.
The proof can be found in the Appendix.

Dimensional Analysis and Simplifications.
Next we check the orders of magnitude of the coefficients in Eqs. (6.12) by a dimensional analysis. We assume that each agent has the same mass m = m i . We consider dimensionless variables x = x/x 0 , t = t/t 0 such that x 0 /t 0 = u 0 is the typical speed of an agent. With this, we define the dimensionless quantities Now, we assume that the range of interaction of K is given by R, so we can write We introduce the dimensionless variable R = R/x 0 so that K =K(|x − y |/R ). Changing variables and expressing the system (6.12) in the prime variables we obtain, after skipping the primes, the following system: where all the variables and parameters are now dimensionless and Notice that the constant λ remains unchanged with respect to the original equation; it is already a dimensionless quantity. The parameter c is a measure of the particle inertia, whereas Re is a measure of the fluid inertia.
Remark 6.2 (Reduction to Vicsek-Stokes coupling). The Vicsek-Stokes coupling (2.1a)-(2.1e) is obtained from the previous system in the regime where Re 1 as well as c 1, η 1. This corresponds to physical systems where the size (and mass) of the agents is very small (microscopic sizes). Therefore, as a simplification we can consider directly Re = 0, c = 0, 1/η = 0 thus removing the inertial terms in the Navier-Stokes equation and the force term −c N i=1 F i δ xi(t) ; as well as relaxing the velocity of the particles to u i = v(X i , t) + aω i . Typically the coefficient b will not be small and should not be simplified.

The Mean-Field Limit.
From now on, we will consider the large friction limit regime defined as follows: Definition 6.3 (Large friction limit regime). The large friction limit regime corresponds to the friction coefficient η → ∞ in the system (6.13) (but leaving c and Re to be O(1)). Then Eq. (6.13b) is replaced by and the rest of equations in (6.13a)-(6.13g) remain unchanged.
This section is devoted to proving the following: Proposition 6.4 (Mean-field limit at finite Reynolds number and finite particle inertia). Given N particles, consider the following scaling of the constant c in Eq. (6.13f): Consider also the empirical distribution associated to the dynamics of the agents in (6.13) in the regime of large friction coefficient (Definition 6.3) with the previous scaling for c, i.e.: where δ xi(t) (x) and δ ωi(t) (ω) denote the Dirac delta at x i (t) and ω i (t) on R 3 and S 2 , respectively. Assume that f N converges weakly to f = f (x, ω, t) as the number of agents N → ∞. Then, the limit f satisfies the following system: where the density ρ f , the flux j f and the Q-tensor Q f are given by Id f dω, (6.16) andω f is given in Eq. (3.3).
Remark 6.3. We must assume that c = O(1/N ) as the number of particles N → ∞. This is because in a mean-field limit interacting terms scale like 1/N so that their sum acting on a given particle remains finite.
Proposition 6.4 is consequence of the following two lemmas: P. Degond et al. JMFM Lemma 6.5. Consider the large friction limit regime in Definition 6.3. The density ρ f and the flux j f given in Eq. (6.16) satisfy the following equations: Lemma 6.6. Consider the large friction limit regime in Definition 6.3. Consider also the scaling for the constant c given in Eq. (6.14). Then, the mean-field limit of the force term in Eq.
The proof of these two Lemmas is given at the end of this section. We prove first Proposition 6.4: Proof of Proposition 6.4. The mean-field limit equation for the density f is computed analogously as in Proposition 3.1. We just need to compute the mean-field limit equation for the fluid velocity v in Eq. (6.13f) and this is done in Lem. 6.6.
Proof of Lemma 6.5. As in Proposition 3.1, we have that the density f satisfies Eq. (6.15a). To obtain Eq. (6.17) for ρ f we integrate this equation with respect to ω. To obtain Eq. (6.18) for the flux j f we multiply the kinetic equation (6.15a) by ω and integrate over ω: where we used that Next, we recast the last two terms of this equation. Firstly, it holds that using integration by parts and the fact that the laplacian in the sphere satisfies Δ ω (ω · u) = −2(ω · u) for any vector u ∈ R 3 (this is the spherical harmonic of degree 1 in S 2 , see [21]). Secondly, it holds that )ω]f dω.
A proof of the last equality can be found in Proposition A.1. Substituting this last expression and Eq. (6.22) into Eq. (6.21) we conclude Eq. (6.18) for j f .
Proof of Lemma 6.6. We consider the following decomposition: For the limit of T N 1 as N → ∞, we have, using (6.13f) and ignoring the Dirac deltas (in Newtonian mechanics, self-forces are ignored to keep the expressions finite) that To compute the limit of T N 2 we first recast the stochastic differential equation (6.13c) for ω i , which is expressed in Stratonovich sense, in its equivalent Itô's form (see [33,Th. (30.14) p. 185], also [3]): With this, we consider the decomposition of T N 2 into To compute the limit of T N 21 we define With these notations we rewrite: and whereω f is given in Eq. (3.3). This leads to For the term T N 22 we have that For any test function ϕ = ϕ(x, ω) we have that where ·, · denotes the duality brackets. Now it holds that The term dB t = B t+dt − B t denotes Brownian motion increments, by the properties of Brownian motion, we have that dB t is normally distributed with mean 0 and variance dt, i.e., dB t ∼ N (0, dt). For fixed t, the following term is a gaussian random variable since it is the sum of independent gaussian random variables (notice that for fixed t, ω i (t) takes a particular fixed value and it is not random). Particularly, its expectation E is zero: ϕ(X i (t), ω i (t)) = 0, (since E(dB 1 t ) = 0) and, moreover, since the Brownian motions are independent (and hence E(dB i t dB j t ) = 0 if i = j), it holds that the variance is zero too in the limit N → ∞: where we used that E(dB 1 t ) 2 = dt and the fact that One can show analogously that the term (ω i · dB i t )ω i in Eq. (6.23) satisfies the same properties since each component of B i t is also a Brownian motion (in 1-dimension). From this we conclude, that Finally, one can see following analogous computations to the previous ones that Putting all the terms together we conclude the proof of statement (6.19). We prove next Eq. (6.20). Using Eq. (6.19), the mean-field limit for the fluid velocity v (6.13f) corresponds to: Now, using that ∇ x · v = 0, as well as ∇ x · (v ⊗ v) = (v · ∇ x )v and the equation for the density ρ f in Eq. (6.17), the previous expression is recast into Finally, from this expression we obtain Eq. (6.20) using Eq. (6.18) for the flux j f and the fact that The second term in (6.25c) corresponds to the momentum flux and it is divided in two contributions. Firstly, corresponds to the momentum flux generated by the fluid and by the passive transport of the particles by the fluid. Secondly, the term corresponding to gives the momentum flux through the exchange between fluid velocity v and particles velocity c 1 aΩ. Notice that the momentum flux is given by a symmetric matrix. The term (a 2c + b)Q gives an extrastress tensor coming from the active nature of the particles and splits into a contribution coming from the dipolar force exerted by the particles (corresponding to the contribution given by the constant b), on the one hand, and from their net force (corresponding to the contribution given by the product a 2c ), on the other hand.

Conclusions
In this paper we have presented the macroscopic derivation of a coupled Vicsek-Stokes system. This coupling describes collective motion in a fluid in a low Reynolds number regime. The fluid is described by Stokes system and the collective motion by the Vicsek model, which represents phenomenologically the interactions between neighbouring agents mediated by the fluid. The coupling is obtained by taking into account the interactions between the agents and the fluid. This involves, particularly, Jeffery's equation that expresses the influence of a viscous fluid on spheroidal particles, on the one hand, and the force exerted by the agents on the fluid due to the dipolar force created by their self-propulsion motion, on the other hand.
The coarse-grained model corresponds to a Self-Organised Hydrodynamics and Stokes coupling. Interestingly, we have shown that Jeffery's equation is coarse-grained into Jeffery's equation but with a different value for the shape parameter. The linear stability analysis shows that both pullers and pushers have unstable modes, but the instability of pullers disappears in the case of rod-like particles.
At the end, we have extended the Vicsek-Stokes coupling into two directions: firstly, we take into account volume exclusion to avoid concentration effects in the dynamics; secondly, we consider a finite Reynolds number and finite particle inertia regime to model systems where the particles' mass and size is large such as fish.
Finally, these results open many exciting paths to be explored, for example, one could consider the coupling of the Vicsek model with other types of fluid dynamics (given e.g. by Darcy's law, Brinkmann law, non-Newtonian fluids). Also, it would be interesting to perform numerical simulations of the dynamics to confirm the stability analysis and to apply these models to the investigation of real-life systems like sperm and bacterial suspensions.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Data Statement
No new data was generated in the course of this research.

A Some Proofs and Properties
Proof of Proposition 6.2. The total momentum is given by ρ 0 v(x, t)dx + N i=1 m i u i . It is a direct computation to check that its derivative is zero. The total angular momentum for the system is given by: We have that Id ∇ x δ Xi dx =: I 1 + I 2 + I 3 + I 4 + I 5 .
One can check directly with the help of the Lévy-Civita symbol to compute the vector products (and integration by parts in some cases) that I 1 = I 2 = I 3 = I 5 = 0, Notice, that I 5 = 0 thanks to (ω i ⊗ ω i − Id/3) being a symmetric matrix. Therefore, the only term that does not vanish is I 4 and it is compensated by the angular momentum of the agents.
Proposition A.1. For any vector u ∈ R 3 , it holds Proof. This can be proven as follows: for any vector q ∈ R 3 q · given that ∇ ω (ω · q) = P ω ⊥ q and P ω ⊥ q · P ω ⊥ u = q · P ω ⊥ u for any pair of vectors q, u ∈ R 3 . From this we conclude the result.