Anton Peshkov

Eric Bertin

Francesco Ginelli

Hugues Chaté

Boltzmann-Ginzburg-Landau approach for continuous descriptions of generic Vicsek-like models

We describe a generic theoretical framework, denoted as the Boltzmann-Ginzburg-Landau approach, to derive continuous equations for the polar and/or nematic order parameters describing the large scale behavior of assemblies of point-like active particles interacting through polar or nematic alignment rules. Our study encompasses three main classes of dry active systems, namely polar particles with 'ferromagnetic' alignment (like the original Vicsek model), nematic particles with nematic alignment ("active nematics"), and polar particles with nematic alignment ("self-propelled rods"). The Boltzmann-Ginzburg-Landau approach combines a low-density description in the form of a Boltzmann equation, with a Ginzburg-Landau-type expansion close to the instability threshold of the disordered state. We provide the generic form of the continuous equations obtained for each class, and comment on the relationships and differences with other approaches.

Introduction

Active matter is currently a very... active field. Under this vocable one usually has in mind collections of "active units" able to extract and dissipate energy from their surrounding to produce systematic and often persistent motion [1]. The far-fromequilibrium collective dynamics of systems as diverse as vertebrate groups [2], insects swarms [3], colonies of bacteria [START_REF] Zhang | Proc. Natl. Acad. Sci. USA[END_REF], motility assays [5,[START_REF] Sumino | [END_REF], as well as driven granular matter [7,8] are thus nowadays routinely modeled by collections of locally interacting active particles.

Many such models have now been proposed and studied with various degrees of sophistication and intended realism, but the Vicsek model stands out for both its simplicity and its historical role in marking the irruption of physicists in collective motion arXiv:1404.3275v1 [cond-mat.stat-mech] 12 Apr 2014 studies. In 1995, Vicsek et al. considered constant-speed point particles aligning with neighbors in the presence of noise. This driven-overdamped, first-order dynamical rule, in spite of its minimality, bears some relevance in all situations where the fluid surrounding the particles can be safely neglected, e.g. when the particles are moving on a substrate acting as a momentum sink. Many of the examples listed above fall into this category of "dry active matter".

Even within the seemingly restrictive setting of competition between (effective) alignment and noise, a wealth of important results has been gathered over the recent years, the most generic ones being the presence of long-range correlations and anomalously-strong "giant" number fluctuations in homogeneous, orientationally ordered, active phases, and the spontaneous emergence of large-scale high-density/high order structures in the region bordering the onset of order. Furthermore, a picture of three possible classes of models has been suggested, defined by the symmetries of particles and alignment rules. The "polar class" is that of the original Vicsek model: it deals with polar particles aligning ferromagnetically. In "active nematics", particle carrying a uniaxial director move back and forth along it and align nematically. Finally, polar particles aligning nematically can be thought of as self-propelled "rods".

Continuous theories have been proposed for these three classes, aiming to capture the long-wavelength behavior of the slow modes of the dynamics. This program was first carried over by Toner and Tu [9,10] who formulated a phenomenological hydrodynamic theory for the polar class based on conservation laws and symmetry considerations. Extension to active nematics was performed in Ref. [11] along the same lines. While successful at describing fluctuations of homogeneous ordered phases, this phenomenological approach is of course unable to build accurate connections to particle-based models, such as the values of the transport coefficients and, crucially, their dependence on hydrodynamic fields. As a result, they typically miss the density-segregated regime and its nonlinear structures.

Other approaches, rooted in a direct coarse-graining of the microscopic dynamics [12,13,14,15,40,16,41], on the other hand, are able to determine most transport coefficients and their functional dependence on the fields. The "Boltzmann-Ginzburg-Landau" (BGL) framework, which we detail below, stems from the early Boltzmann approach of [13] to the polar case, and combines it with the traditional Ginzburg-Landau weakly nonlinear analysis. We argue below that it offers better overall control, something needed since some confusion remains: the equations obtained by different methods often differ not only in details but also in structure. For instance, the equations for rods derived by Baskaran and Marchetti from a Smoluchowski equation contain more and different terms than those derived in the BGL framework. The Chapman-Enskog formalism put forward in [16] for the polar case yields many more terms than those retained in other approaches, and their effect on the dynamics, especially at the nonlinear level, remains unclear.

Recent work has shown, we believe, that the Boltzmann-Ginzburg-Landau framework is near-ideal for deriving minimal, well-behaved, nonlinear hydrodynamic descriptions from Vicsek-like models, without aiming at a quantitative agreement with what remains, after all, rather unrealistic starting points. The equations derived for the three basic classes however have, we believe, universal value and have been shown to account remarkably well for most features observed, including highly nonlinear solutions and chaotic regimes.

In the following, we present a unified formulation of the BGL approach which is applicable to generic dry active matter systems in their dilute limit and encompasses both ballistic propagation and positional (anisotropic) diffusion. This allows one to treat both propagative (polar case, rods) and diffusive (active nematics) active systems, as well as systems where non-equilibrium activity and "thermal" (random) agitation compete.

In Section 2, we define the generalized Vicsek model we start from and derive the associated Boltzmann kinetic equation. In Section 3, we derive hydrodynamic equations for the three basic classes. In Section 4, we discuss the BGL framework, the obtained results, and compare them to other works.

Vicsek-like models and their Boltzmann description

In dry systems, motion takes place over a substrate or in a fluid which acts as a momentum sink, so that overall (active particles plus fluid) momentum conservation is not a concern and long-ranged hydrodynamic interactions between particles can be neglected (hence the term dry). Over the elementary dissipative timescale ∆t of the dynamics, particle motion is coupled to particle orientation, leading either to self-propulsion or to active diffusion. From a modeling viewpoint, active particles are represented here by point-like objects interacting locally through alignment 'forces' in the presence of noise. The alignment and the noise are a simple means to model the (possibly) complicated interactions between active particles in real systems. Due to the strongly dissipative nature of the substrate/surrounding fluid, our point particles follow an effective first order dynamics, both positionally and rotationally.

Definition of a generalized Vicsek model

We consider active particles moving in a two dimensional continuous space. Particles are described by their position r, and a (unit) heading vector defined by an angle θ ∈ [-π, π]. In a generalization of the classical VM streaming rule, particles positions are updated according to r(t

+ ∆t) = r(t) + v ∆t (1)
where v is a random vector drawn from a displacement distribution Φ(v, θ v -θ), with v = |v| and θ v being the angle defining the orientation of v. In the absence of chirality one has Φ(v, θ v -θ) = Φ(v, θ -θ v ). Purely ballistic motion with speed v 0 (standard VM dynamics), for instance, is recovered considering the delta distributions

Φ(v, θ v -θ) = δ(v -v 0 )δ(θ v -θ).
The angular part of the distribution Φ(v, θ v -θ) is obviously 2π-periodic, but higher order symmetries may reflect further symmetries in particles displacement. Active nematic particles, for instance, are characterized by π-periodicity (and of course no higher order symmetries). For compactness, we shall say that Φ is mπ-periodic, with m = 1 in the nematic case and m = 2 in the polar case.

The heading angle θ of a given particle evolves according to the simple stochastic dynamics,

θ = Ψ (p) (θ, θ i1 , θ i2 , . . . , θ ip ) + η , (2)
where η is a symmetrically-distributed (zero mean) random angle (delta correlated in time and between different particles) and θ i1 , θ i2 , . . . , θ ip are the angles of the p interacting particles which, in metric models, are typically defined as the particles closer than the metric interaction range r 0 . When no neighbours are present, particles simply experience self-diffusion events

Ψ (0) (θ) = θ + η.
Contrary to the orientational degree of freedom, the modulus of the displacement vector v has essentially no dynamics, and is assumed to be a hard mode1 with only fast fluctuations around its constant mean value.

Boltzmann equation

The above generic model can be described in the standard framework of kinetic theory, through the introduction of a (generalized) Boltzmann equation, which describes the time evolution of the one-particle distribution function f (r, θ, t), namely the probability to find a particle at position r at time t. We choose the normalization

1 V V dr π -π dθf (r, θ, t) = ρ 0 (3)
with V the system volume and ρ 0 its mean density.

The Boltzmann approach relies on the molecular chaos hypothesis and the binary collision approximation, which, to a certain extent, are justifiable in low density systems. We will comment on these approximations in section 4.

Positional part

To obtain the Boltzmann equation, we first consider the (generalized) "streaming" part of the dynamics, described by Eq. ( 1). One has

f (r, θ, t + ∆t) = dv Φ(v, θ v -θ)f (r -v∆t, θ, t). (4)
On time scales much larger than ∆t, applying Itô calculus [START_REF]Oksendal Stochastic Differential Equations[END_REF] to second order, one obtains the master equation

∂f ∂t = -v α ∂ α f (r, θ, t) + 1 2 δv α δv β ∆t ∂ α ∂ β f (r, θ, t) (5)
where summation over the repeated spatial indices α, β is understood and the brackets . . . denote an average over the distribution

Φ(v, θ v -θ) v α = ∞ 0 dv π -π dθ v vn α (θ v ) Φ(v, θ v -θ) (6)
δv α δv β = ∞ 0 dv π -π dθ v δv α (θ v )δv β (θ v )Φ(v, θ v -θ) (7)
with n(θ) = (cos θ, sin θ) T the unit vector in the direction θ, v = vn(θ) and

δv α = v α -v α .
For polar particles, the vector v points in the direction defined by θ. One can then write v = v 0 n(θ), which defines the average speed of particles v 0 = v cos δθ (δθ ≡ θ v -θ). For nematic particles v = 0, and one can thus also formally write v = v 0 n(θ), with v 0 = 0. The covariance term δv α δv β is less straightforward to evaluate. After some algebra we find

δv α δv β = 1 2 v 2 (1 -cos 2δθ) δ αβ + v 2 cos 2δθ -v 2 0 n α (θ)n β (θ) (8)
which can be recast as the sum of an isotropic (∼ δ αβ ) and an anisotropic (∼ g αβ = n α n β -δ αβ /2) term. Here δ αβ is a Kronecker delta and g αβ is the nematic tensor familiar to liquid crystal physics [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF]. We write the positional part of the Boltzmann equation as

∂f ∂t = -v 0 n α ∂ α f + D 0 ∆f + D 1 g αβ ∂ α ∂ β f ( 9
)
where ∆f is the Laplacian of f and

D 0 = 1 4 v 2 -v 2 0 ∆t (10)
D 1 = 1 2 v 2 cos 2δθ -v 2 0 ∆t (11)
are, respectively, the isotropic and anisotropic diffusion constants. Fluctuations in the displacement vector v, either in modulus or in heading, lead to a finite isotropic diffusion term D 0 > 0. Anisotropic diffusion (D 1 = 0), on the other hand, is produced by anisotropies in the fluctuations of the displacement vector v. In this regard, it is instructive to consider the simple "4 directions" case

Φ(v, δθ) = δ(v -v s ) 4 j=1 p j δ δθ -(j -1) π 2 ( 12
)
where v s is the speed of the particles, p j ≥ 0, and j p j = 1. For illustrative reasons, we can temporarily relax the non-chiral condition and allow p 2 = p 4 . One then has a drift term with mean velocity v = v 0 n(θ + ψ 0 ) where

v 0 = v s (p 1 -p 3 ) 2 + (p 2 -p 4 ) 2 (13)
and

ψ 0 = Arg [(p 1 -p 3 ) + i(p 2 -p 4 )] (14)
(i being the imaginary unit). Diffusion constants can be written as

D 0 = (D + D ⊥ )/2 and D 1 = (D -D ⊥ )
, where

D = p 1 + p 3 -(p 1 -p 3 ) 2 v 2 s 2 ∆t (15)
D ⊥ = p 2 + p 4 -(p 2 -p 4 ) 2 v 2 s 2 ∆t
highlighting the proportionality of D 1 to the difference between parallel and perpendicular fluctuations (w.r.t. local heading).

Full Boltzmann equation

In the low density limit, interactions between more than two particles are unlikely, and the heading angle dynamics is captured by self-diffusion events and binary collisionlike events. In self-diffusion events, the heading angle θ is changed by a random amount η, drawn from a zero-mean, symmetric distribution P σ (η) with standard deviation σ. As a slight simplification, we assume that such self-diffusion events occur with a constant rate λ (typically of the order of 1/∆t), independent of the density. In binary collisions, the incoming angles θ 1 and θ 2 of the two particles are changed, after collision, into

θ 1 = Ψ (θ 1 , θ 2 ) + η 1 , θ 2 = Ψ (θ 2 , θ 1 ) + η 2 (16)
where η 1 and η 2 are random variables drawn from P σ (η). In principle, the collision and self-diffusion noise distributions could be different. To simplify the presentation, we shall consider in what follows that the two distributions are the same; however calculations can be easily carried out with distinct distributions.

The collision rate is encoded in a collision kernel K(∆ = θ 2 -θ 1 ) ≥ 0, which depends only on the angle difference due to a global rotational invariance of the problem. In metric models, where the collision probability can be modeled as a scattering process [START_REF] Huang | Statistical Mechanics[END_REF], the collision kernel K(∆) inherits the same mπ-symmetries as that of the displacement distribution Φ.

Interaction symmetries also impose constraints on the interaction rule Ψ (θ 1 , θ 2 ), which can typically be either polar or nematic depending on its periodicity with respect to each angle θ 1 and θ 2 . Formally, Ψ (θ 1 , θ 2 ) is nπ-periodic, modulo nπ, with respect to both arguments θ 1 and θ 2 independently, which defines the symmetry index n of the interaction rule. In the following, we call "ferromagnetic" such an alignment rule with 2π-periodicity (n = 2) and nematic an alignment rule with π-periodicity (n = 1). The interaction rule must also be consistent with the displacement symmetry of particles (reflected in the periodicity of the displacement distribution Φ and of kernel K), which implies n ≤ m.

In addition, isotropy imposes that for an arbitrary rotation of angle φ,

Ψ (θ 1 + φ, θ 2 + φ) = Ψ (θ 1 , θ 2 ) + φ [nπ]. ( 17
)
Choosing φ = -θ 1 and θ 2 = θ 1 + ∆, one obtains

Ψ (θ 1 , θ 1 + ∆) = θ 1 + Ψ (0, ∆) [nπ]. ( 18
)
Hence the interaction rule Ψ (θ 1 , θ 2 ) is parameterized by a single variable interaction function H(∆) ≡ Ψ (0, ∆), which is nπ-periodic. We further assume that there is no chirality in the problem and that the particle exchange symmetry is respected by point particles collisions, implying that

K(-∆) = K(∆) and H(-∆) = -H(∆) [nπ].
Due to these symmetry properties, we need only define K(∆) and H(∆) over the positive half of their definition interval.

The full Boltzmann equation is obtained by taking into account the self-diffusion and binary collision rules (assumed to be valid in the low density regime considered here),2 under the molecular chaos assumption which approximates the two-particle distribution f (2) as the product of two one-particle distributions

f (2) (r, θ 1 , θ 2 , t) ≈ f (r, θ 1 , t)f (r, θ 2 , t). (19)
Physically, this assumption means that the heading angle decorrelates between consecutive collisions. One then finds for the Boltzmann equation

∂f ∂t + v 0 n α ∂ α f = D 0 ∆f + D 1 g αβ ∂ α ∂ β f + I dif [f ] + I col [f ] ( 20
)
where the self-diffusion integral I dif [f ] and the collision integral I col [f ] are defined as

I dif [f ] = -λf (r, θ, t) + λ π -π dθ ∞ -∞ dη P σ (η)δ mπ (θ + η -θ)f (r, θ , t) (21)
I col [f ] = -f (r, θ, t) π -π dθ K(θ -θ) f (r, θ , t) (22)
+ π -π dθ 1 π -π dθ 2 ∞ -∞ dη P σ (η) K(θ 2 -θ 1 )f (r, θ 1 , t)f (r, θ 2 , t) × δ mπ Ψ (θ 1 , θ 2 ) + η -θ
where δ mπ is a generalized Dirac delta imposing that the argument is equal to zero modulo mπ.

Without any loss of generality, we can set λ = 1 (this amounts to a rescaling of time). For convenience, we further consider the noise distribution P σ to be Gaussian. The physics described by the Boltzmann equation is determined by just three functions encoding the microscopic dynamics: the displacement distribution Φ, the collision kernel K and the interaction function H.

Three basic classes

A trivial case for the Boltzmann equation ( 20) is given by a displacement distribution Φ with an isotropic angular part, for which one has v 0 = 0, D 0 ∼ d 2 0 /4 and D 1 = 0. This corresponds to an isotropic diffusive motion, which can be mapped onto a passive, equilibrium system. Here displacement is strictly random and completely decoupled from orientational dynamics.

Anisotropies in Φ, on the other hand, lead to non-equilibrium activity. Although this does not cover all possible cases, it is relatively easy to identify the three simplest classes depending on the symmetries of the particles displacement (mπ-periodicity of the displacement distribution Φ) and of the interactions (nπ-periodicity of the interaction rule Ψ ). These three classes arise from the fact that the symmetry (polar or nematic) of both the particles and the interaction rule can be varied, with however the constraint that the interaction rule must obey at least the particle symmetry (n ≤ m).

Polar particles with ferromagnetic interaction

This class, labeled as m = n = 2, is exemplified by the Vicsek model for self-propelled particles. Displacement distributions breaking the nematic symmetry θ v → θ v + π result in a nonzero drift, v 0 = 0. Such active particles are typically described as self-propelled. In metric systems dominated by drift, K(∆) is 2π-periodic (polar symmetry); for an isotropic interaction range, it is given by the flux of incoming particles through the cross-section 2r 0 of a target particle [13],

K(∆) = 2v 0 r 0 |n(θ 2 ) -n(θ 1 )| = 4v 0 r 0 sin ∆ 2 . (23)
However, more complex forms of K(∆) (typically involving higher harmonics in ∆) can be used to describe anisotropic interaction ranges (as the one used to model collisions of elongated objects in a point-like framework [START_REF] Weber | [END_REF]).

The interaction function H(∆) describing ferromagnetic interactions is 2π-periodic. For the "canonical" Vicsek model studied in [22,23], the function H(∆) is given by H(∆) = ∆/2 and diffusion terms strictly vanish, with Φ(v, δθ) = δ(v -v 0 )δ 2π (δθ) implying D 0 = D 1 = 0. If the displacement distribution is characterized by a finite variance, non-zero diffusion is added to drift. Isotropic diffusion is characterized by D 0 > 0 and can be interpreted as having a thermal origin. The anisotropic part of diffusion, on the other hand, arises from fluctuation anisotropies and has (in the present framework of point-like particles) a strictly non-equilibrium origin [11]. Indeed, speed fluctuations in the presence of drift are enough to generate anisotropic diffusion even in the absence of heading fluctuations.

Polar particles with nematic interaction

This class, labeled as m = 2, n = 1, corresponds to self-propelled (v 0 = 0) "rods", namely self-propelled particles with an elongated shape making their alignment nematic. The kernel K(∆) is 2π-periodic (to zeroth order in particle anisotropy), and for metric systems is still given by Eq. ( 23), but H(∆) is now π-periodic. For the representative model of this class, studied in [24,25], one has H(∆) = ∆ 2 , and Φ(v, δθ) = δ(v -v 0 )δ 2π (δθ).

Nematic particles with nematic interaction

This case m = n = 1 correspond to the "active nematics" class. Here the displacement distribution is π-symmetric and drift is strictly zero, with the system dominated by finite isotropic and anisotropic (non-equilibrium) diffusion. Both K(∆) and H(∆) are π-periodic. In the simple case of isotropic interaction range, the polar kernel (23) can be modified to account for nematic symmetry

K(∆) = 2v 0 r 0 sin ∆ 2 + cos ∆ 2 . ( 24
)
A representative model of this class studied in the literature [26,27] is defined by

H(∆) = ∆ 2 and Φ(v, δθ) = δ(v -v 0 )δ π (θ).

Metric-free models

So far, we have mostly discussed models characterized by metric interaction ranges. Motivated by recent results on the observation of starling flocks [28], and fish schools [29], models where the interaction range is not defined by a metric distance, but rather by a typical number of neighbors, are of interest. Metric-free versions of the Vicsek model, where neighbours are not chosen inside a metric range but rather by some topological criteria (for instance, making use of Voronoi tessellations) have been introduced in [30]. From a theoretical perspective, a difficulty with such models is that a particle always interacts with its neighbors, whatever their metric distance. In order to make such a situation suitable for a kinetic theory approach, it has been proposed to introduce a low interaction rate such that at each time step, two neighboring particles (in a topological sense) have only a small probability to interact [31].

Once two particles have interacted, their velocity angles are updated according to Eq. ( 16). The difference with the metric cases appears in the collision kernel K. Here interactions cannot be described by a scattering process, and the probability of collision does not depend on the angles of the two interacting particles. In metric models, K is defined under the implicit assumption that the probability of collision is proportional to the local particles density ρ(r, t) ≡ π -π dθf (r, θ, t). On the other hand, for metric-free models, the probability of interaction is independent of the local density and to take this crucial effect into account, one uses the framework of metric models, but with an effective interaction kernel proportional to 1/ρ.

Hydrodynamic equations

In order to derive hydrodynamic equations from the Boltzmann equation, we proceed as follows. First, we expand the Boltzmann equation in angular Fourier modes, yielding an infinite hierarchy of equations for these modes. Second, we look for the first linearly unstable angular mode when decreasing the noise and increasing the density, and truncate the hierarchy of equations in the vicinity of this linear instability. The continuous description thus includes the (conserved) local density field and the relevant angular modes (typically the linearly unstable mode).

The specificity of the "Boltzmann-Ginzburg-Landau" approach is to rely on a systematic scaling ansatz to truncate and close the hierarchy to a set of minimal equations.

Angular Fourier modes expansion

Here we chose to restrict ourselves to two spatial dimensions, where the use of Fourier transforms and complex notation greatly simplifies calculations. In order to define a framework consistent with the most general 2π-symmetry, we work with an angular Fourier expansion defined over the 2π interval:

f (r, θ, t) = 1 2π ∞ k=-∞ fk (r, t) e -ikθ , (25)
where the Fourier coefficients fk are defined as

fk (r, t) = π -π dθ f (r, θ, t) e ikθ . ( 26
)
Note that f0 (r, t) is nothing but the local density ρ(r, t). Note also that for all k, f-k = f * k (the star denotes the complex conjugate). An additional π-periocity of the one particle distribution f then manifests itself by the nullity of odd Fourier coefficients, f2l+1 = 0. In the following, we drop the 'hat' over the Fourier coefficients in order to lighten notations.

It is easy to verify that the Fourier modes f k (k ≥ 1) can be interpreted as order parameter fields associated to specific spontaneous symmetry breakings. In particular, f 1 encodes the momentum field, and f 2 the nematic tensorial field:

ρP = Ref 1 Imf 1 , ρQ = 1 2 Ref 2 Imf 2 Imf 2 -Ref 2 , ( 27
)
where P is a polarity field of components P α = n α l , and Q the traceless tensorial field of components Q αβ = n α n β l -δ αβ /2 (here . . . l denotes a local average over the distribution f .) Introducing the complex derivatives ∇ = ∂ x + i∂ y and ∇ * = ∂ x -i∂ y as well as the Laplacian ∆ = ∇∇ * (note that in this notation ∇ 2 = (∂ x + i∂ y ) 2 is not the Laplacian), the angular Fourier expansion of the Boltzmann equation then yields the infinite hierarchy of coupled equations

∂f k ∂t + v 0 2 (∇f k-1 + ∇ * f k+1 ) = -(1 -P k )f k + D 0 ∆f k + D 1 4 (∇ 2 f k-2 + ∇ * 2 f k+2 ) + ∞ q=-∞ (P k I k,q -I 0,q )f q f k-q (28)
where the coefficient I k,q is defined by the integral

I k,q = 1 2π π -π d∆ K(∆) e -iq∆+ikH(∆) (29)
and

P k ≡ P k (σ) = ∞ -∞ dη P σ (η) e ikη
is the Fourier transform of the noise distribution (restricted here to integer values of k). One has 0 ≤ P k (σ) ≤ 1 and P k (0) = 1 ∀k. For a Gaussian noise distribution, the Fourier transform has the simple form

P k (σ) = e -σ 2 k 2 /2 . ( 30
)
Note also that due to the parity properties of K(∆) and H(∆), I k,q is real. From Eq. ( 28), we find in particular the continuity equation, obtained for k = 0

∂ρ ∂t + v 0 Re(∇ * f 1 ) = D 0 ∆ρ + D 1 2 Re(∇ * 2 f 2 ). ( 31
)
Note that this equation is valid without any further assumption, due to the fact that the Fourier transform of the integral terms in Eq. ( 20) vanish for k = 0.

3.2 Linear instability of the disordered state

General considerations

The homogeneous disordered isotropic state

f 0 = ρ 0 , f k = 0 (k ≥ 1
) is a trivial solution of Eq. ( 28). Assuming spatial homogeneity, we linearize Eq. ( 28) around the isotropic solution, yielding for k > 0

∂f k ∂t = [-(1 -P k ) + ω k ρ 0 ] f k (32)
where we have defined

ω k = P k (I k,k + I k,0 ) -(I 0,k + I 0,0 ). (33)
The linear stability of f k is governed by the sign of the linear coefficient,

µ k (σ, ρ 0 ) ≡ -(1 -P k ) + ω k ρ 0 . ( 34
)
Analyzing the full noise and density dependence of this expression for arbitrary kernel and interaction rule is a complicated task. One can however notice that at low enough noise P k → 1 and the sign of µ k is then given by that of ω k , since µ k ≈ ω k ρ 0 in this limit. By studying ω k instead of µ k , one thus gets rid of the ρ 0 dependence. Therefore, the sign of ω k is controled by the zero noise limit ω

(0) k = (I k,k + I k,0 ) -(I 0,k + I 0,0 ): if ω (0) k = 0,
there exists a finite range of σ, in the vicinity of σ = 0, where ω k has the same sign as ω (0) k . Hence, the study of the quantity ω (0) k , which depends only on the kernel and on the interaction rule (and not on the noise and density) should be enough to detect an instability of f k . Of course, this simplification imposes several limitations. First, by construction, the study of ω (0) k does not allow in itself for the identification of an instability line in the parameter plane (σ, ρ 0 ). Second, it provides no information on which mode is the most unstable close to the instability threshold in case several modes would be unstable in the high density and low noise limit. However, if a single mode is found to be unstable in this limit, one is sure that this mode should be taken as the order parameter of the transition.

An explicit expression of ω

(0)
k is given by

ω (0) k = 1 π π 0 d∆K(∆) cos k(H(∆) -∆) + cos kH(∆) -cos k∆ -1 . (35)
For any given specific kernel and interaction rule, one can compute (analytically or numerically) the coefficients ω (0)

k for k = 1, 2, . . . , and evaluate their sign to determine the first unstable mode. In particular, note that in metric-free models, due to the dependence of the effective kernel on density one has ω (0) k ∼ 1/ρ, resulting in a linear coefficient µ k (σ) independent of density.

First unstable mode: study of generic classes

In the following, we try to determine the sign of ω (0) 1 and ω (0) 2 as a function of the symmetries of the kernel K(∆) and the interaction rule H(∆). The fully general case is hard to analyze, but we are able to obtain fairly general results for each of the three classes introduced in Sect. 2.3, under the additional assumption that the interaction rule is indeed an alignment rule for all values of the angle difference ∆, which translates in mathematical terms as 0 ≤ H(∆) ≤ ∆ for 0 < ∆ < nπ 2 , with 0 < H(∆) < ∆ over at least a finite subinterval of ∆. The main results can be summarized as follows.

-Polar particles with ferromagnetic alignment One generically finds that ω (0)

1 > 0, while the sign of ω (0) 2 cannot be determined without further assumptions on the interaction rule. If H(∆) is 'close enough' to the standard rule

H(∆) = ∆ 2 , one recovers ω (0) 2 < 0.
It is likely that close to the transition line ω 1 = 0 in the noise-density plane, one has ω 2 < 0 so that only polar order is unstable, but checking this for arbitrary rules within this class is a difficult task. In any case, this result shows that any generic ferromagnetic aligning rule which respects the particle exchange symmetry H(-∆) = -H(∆) is enough to make the disordered solution unstable towards the growth of polar order.

-Polar particles with nematic alignment

This case somehow mirrors the previous one, exchanging the roles of ω (0) 1

and ω (0) 2 . Under the generic assumptions of this class, one finds that ω (0) 2 > 0, pointing to a nematic instability for generic nematic alignment rules. The sign of ω (0) 1 is not completely fixed within the class. However, under the fairly natural assumption that K(∆) is a non-decreasing function in [0, π 2 ], as is the case of the kernel (23), one can show that ω (0) 1 < 0. This result remains true if K(∆) decreases only slightly.

-Nematic particles with nematic alignment

Similarly to the case of polar particles with nematic interaction, one finds ω (0) 2 > 0 for generic nematic alignment rules. The situation is however simpler here, since f 1 = 0 by symmetry, so that its dynamics need not be studied.

Let us emphasize that the 'contracting' condition on H(∆), namely 0 ≤ H(∆) ≤ ∆, is only a sufficient condition allowing for easier calculations, but instabilities of the disordered state can also occur if H(∆) is not 'contracting' everywhere. As an illustrative example of a polar rule that is not 'contracting' everywhere, let us consider the function H(∆) defined as H(∆) = 2∆ for 0 < ∆ < aπ and H(∆) = ∆ 2 for aπ < ∆ < π, with 0 < a < 1. One then finds for ω (0) 1 , assuming for simplicity that

K(∆) is a constant K 0 , independent of ∆, ω (0) 1 = K 0 π 4 -π + 1 2 sin 2aπ + sin aπ -4 sin aπ 2 (36)
which reduces for small a to

ω (0) 1 ≈ K 0 π 4 -π - 3 4 a 3 π 3 . (37)
For small enough a, ω

1 is thus positive, and the polar order parameter f 1 is unstable at low noise and/or high density.

Scaling ansatz close to instability threshold

Our goal is to reduce the infinite hierarchy (28) of equations to a small set of minimal equations for the relevant hydrodynamic fields. The number of particles being conserved, the density ρ(r, t) is a relevant field, for which we already have derived the evolution equation (31). Additional relevant fields correspond to broken symmetries, and can be identified from the linear instabilities studied above. In order to obtain a finite number of equations, one needs to truncate the infinite hierarchy of equations, which requires some truncation criterion. In the spirit of Ginzburg-Landau equations, we shall use a scaling ansatz close to the instability threshold of the first unstable mode in which the order parameter fields are small quantities and, being interested in the long-wavelength long-timescale dynamics, also time and space derivatives of the fields are smaller than the field themselves.

In the following, we shall first review how to obtain this scaling ansatz in each of the three basic classes outlined above, and next we will use what we have learned to formulate a systematic scaling ansatz valid in the general case. The simplest case is that of active nematics, so we start by considering this case.

Nematic particles with nematic interaction

The only relevant order parameter here is f 2 . We assume that close to the instability threshold, f 2 is of the order of a small parameter , related to the small value of the linear coefficient µ 2 > 0. The continuity equation (31), here with f 1 = 0, imposes two constraints: (i) density fluctuations δρ = ρ -ρ 0 with respect to the average density ρ 0 are of order ; (ii) if spatial derivatives ∇ are of order α , then the time derivative ∂ t scales as 2α . Nonlinear saturation should result from the interaction between f 2 and the next non zero mode f 4 . Assuming that the coefficient µ 4 is negative, f 4 should be slaved to f 2 in the equation for

f 4 since |∂ t f 4 | |µ 4 f 4 |.
Hence the linear term in f 4 should be balanced by the term in f 2 2 in the equation for f 4 , yielding f 4 ∼ 2 . Coming back to the equation for f 2 , one has to balance (again in a Ginzburg-Landau spirit) the diffusive term ∆f 2 with the nonlinear term f * 2 f 4 , which fixes α = 1. In summary, one has:

f 2 ∼ , f 4 ∼ 2 , ρ -ρ 0 ∼ , ∇ ∼ , ∂ t ∼ 2 . ( 38
)
Balancing the linear term µ 2 f 2 with the nonlinear term f * 1 f 2 also yields µ 2 ∼ 2 , which determines the relation between and the distance to the instability threshold.

Polar particles with ferromagnetic interaction

This case is slightly more complicated, as will appear below. Now the relevant order parameter is f 1 . We have to balance the following terms in the continuity equation (31) and in the equation for f 1 [see Eq. ( 28)]

∂ t ρ ∼ Re∇ * f 1 , ∂ t f 1 ∼ ∇ρ (39)
which implies ρ -ρ 0 ∼ f 1 ∼ and the propagative scaling ∂ t ∼ ∇ ∼ α . If µ 2 < 0, f 2 is slaved to f 1 and the linear term µ 2 f 2 should be balanced by the nonlinear term f 2 1 , yielding f 2 ∼ 2 . Finally, balancing the nonlinear term f * 1 f 2 with the diffusive term ∆f 1 in the equation for f 1 leads to α = 1, resulting in the following scaling relations

f 1 ∼ , f 2 ∼ 2 , ρ -ρ 0 ∼ , ∇ ∼ , ∂ t ∼ . ( 40
)
In addition, the relation

µ 1 f 1 ∼ f * 1 f 2 yields µ 1 ∼ 2 .
This case is apparently quite similar to active nematics apart from the fact that the diffusive scaling ∂ t ∼ ∇ 2 is replaced by the propagative one ∂ t ∼ ∇, but a major difference is that the resulting scaling does not allow all "important terms" to be balanced. This is true both in the continuity equation, where for instance ∆ρ ∼ 3 while Re∇ * f 1 ∼ 2 , and most importantly in the equation for f 1 , where terms like µ 1 f 1 and f * 1 f 2 are of order 3 while ∇ρ ∼ 2 . If one neglects the terms µ 1 f 1 and f * 1 f 2 , the resulting equation does not lead to collective order, so that one has to keep terms of different order in the same equation to account for the relevant phenomenology.

Although this unbalance of terms may look surprising at first sight, the situation is similar for instance to what happens in the simple Fokker-Planck equation describing a biased random walk. The drift dynamics and the diffusive dynamics occur on different time scales, which implies that the drift and diffusive terms in the equation are generically of different order in an expansion in the lattice spacing. In the presence of a constant drift, one can reabsorb the drift term by going to a moving frame. In more complicated situations like the ones we are dealing with here, it is not easy to reabsorb the terms ∇ρ and Re∇ * f 1 by a simple transformation. A proper way to deal with this issue would probably be to perform a multiscale expansion, introducing fast and slow time variables. Such a derivation however goes beyond the scope of the present paper.

Polar particles with nematic interaction

This last case shares similarities with both previous cases. The linearly unstable mode is f 2 , which is thus the main order parameter of the problem, but one may also wish to keep f 1 as a relevant field because of the polarity of particles. In this case, two different scalings can to some extent be considered as consistent. The first one starts from the same reasoning as in the active nematics case (sect. 3.3.1), leading again to the scaling Eq. (38). Then, to determine the scaling of f 1 (which was absent from the active nematics case), we balance the term ∂ t ρ with Re∇ * f 1 , or equivalently ∂ t f 2 with ∇f 1 , leading to f 1 ∼ 2 . Further, assuming that µ 3 and µ 4 are negative so that f 3 and f 4 are slaved to the nonlinear terms, one has f 3 ∼ f 1 f 2 and f 4 ∼ f 2 2 , yielding f 3 ∼ 3 and f 4 ∼ 2 . Finally, balancing the lowest order nonlinear term coupling to higher modes, f 4 f * 2 , with the diffusive term ∆f 2 in the equation for f 2 leads to α = 1. Altogether, we have

f 1 ∼ 2 , f 2 ∼ , f 3 ∼ 3 , f 4 ∼ 2 , ρ -ρ 0 ∼ , ∇ ∼ , ∂ t ∼ 2 . ( 41
)
This scaling ansatz then leads to equations of the same form as the active nematics case. Under this scaling assumption, f 1 is purely slaved to f 2 (in the sense that ∂ t f 1 is a higher order term), and can eventually be eliminated from the equations. An alternative view is to consider that the propagative nature of the particles suggests a ballistic scaling between time and space, namely ∂ t ∼ ∇. To fulfill this scaling in the continuity equation, one then needs to assume that f 1 ∼ , making it a relevant field. We thus end up with the following scaling ansatz,

f 1 ∼ , f 2 ∼ , f 3 ∼ 2 , f 4 ∼ 2 , ρ -ρ 0 ∼ , ∇ ∼ , ∂ t ∼ (42)
leading to a richer set of closed equations, as shown in section 3.4.3.

Systematic general scaling ansatz

The previous examples suggest the definition of a systematic scaling ansatz valid in the generic case. First of all, density fluctuations are of order , as well as all Fourier modes up to the lowest k relevant order parameters, which we will denote as k = h and which is typically determined by the alignment symmetry, h = 2/n (n = 1 or 2). The next h modes then scale as 2 , and so on and so forth. Space derivatives scaling is fixed by considering the equation for the relevant order parameter(s). Balancing the diffusion term (of order 1+2α ) with the lowest order nonlinear terms coupling to higher modes (always of order 3 ), leads to α = 1 and ∇ ∼ , consistently with the typical Ginzburg-Landau scaling ansatz. Finally, the scaling of time derivatives is fixed by the symmetry of the displacement distribution Φ, or equivalently, by the dominant propagative mode in the continuity equation (31). Non zero drift imposes the propagative scaling ∂ t ∼ , while zero drift, diffusive dynamics leads to the diffusive scaling ∂ t ∼ 2 .

To summarize, we have for our systematic ansatz (with

h = 2/n) ρ -ρ 0 ∼ , f k =0 ∼ 1+ (|k|-1)/h , ∇ ∼ , and (43)
∂ t ∼ (ballistic) or ∂ t ∼ 2 (diffusive),
where . . . denotes the integer part.

Truncation and closed hydrodynamic equations

We now discuss the general structure of the closed equations emerging from a truncation of Eqs. (28) by the general scaling ansatz (43).

Independently of the truncation scheme, the dynamics of the density field fluctuations is given by the continuity equation (31), which couples δρ to the polar (drift) and the nematic (anisotropic diffusion) order parameter fields. The isotropic diffusion term ∼ ∆ρ, in the absence of anisotropic diffusion, can be interpreted as of pure thermal origin; otherwise, it carries a contribution from nematic activity. Note that in 'pure' self-propelled models, like the original VM, where displacement fluctuations are ignored, diffusion contributions (both isotropic and anisotropic) are absent.

Since additional symmetries in the system, and thus in the one-particle distribution f manifest themselves by some zero Fourier coefficients, one can check from Eq. ( 31) that global symmetries of higher than nematic order results in both f 1 = 0 and f 2 = 0, completely decoupling density from order parameter dynamics. This suggests that, at least in two spatial dimensions, no out-of-equilibrium dynamics may emerge from active particles with global symmetries higher than the nematic one, so that m = 1, 2 are the only relevant symmetries.

We now turn our attention to the order parameter equations of order k ≥ 1. They couple the mode f k to terms of order k ± 1, 2 via linear spatial derivatives, and to terms of all orders via the nonlinear coupling term

C k = ∞ q=-∞ (P k I k,q -I 0,q )f q f k-q .
(44)

Terms of order q = 0, k in the sum can be singled out to yield the density dependent linear coefficient

µ k = -(1 -P k ) + ω k ρ (45)
which controls the local stability of the disordered phase3 . In the following we assume that all linear coefficients µ k are negative except mode h = 2/n, associated to the alignment symmetry and relevant order parameter, which for high density and/or low noise may turn positive, µ h > 0.

We have seen above that the first non-linear terms appearing in the relevant order parameter equation are typically of order 3 . As a result, we expand all equations to order 3 . This results in a well-controled truncation close to the instability threshold, in such a way that equations for additional modes which are not retained as order parameters provide closure relations for the order parameter equations.

The relevant order parameter equation couples f h to higher order terms via order 3 terms of C h of the kind f h+j f * j , with j = 1, . . . , h. Higher order terms, in turn, are slaved to order modes since µ k < 0 and

|∂ t f k | |µ k f k |.
In particular, for k = 2h (the first higher non-zero mode when m = n) one has to order 2 (cubic terms may be ignored since inserting them in the nonlinear coupling term in the equation for f h would result in terms ∼ 4 )

f 2h ≈ - 1 µ 2h (P 2h I 2h,h -I 0,h ) f 2 h -δ h,1 v 0 2 ∇f 1 (46)
with the gradient term appearing only for the polar case h = 1. Substituting the slaved mode f 2h into the k = h equation one obtains the saturating cubic term -ξ h f h |f h | 2 with the coefficient

ξ h = P 2h I 2h,h -I 0,h µ 2h [P h (I h,2h + I h,-h ) -I 0,2h -I 0,-h ] (47)
being typically positive in the parameter region where µ h > 0.

If m = n the hydrodynamic equations have a relatively simple structure, consisting of the continuity equation (31) and the relevant order parameter equation of the form

∂f h ∂t = µ h (σ, ρ) -ξ h |f h | 2 f h + ∇(. . . ) + ∇∇(. . . ) (48)
where ∇ and ∇∇ indicate generic spatial derivatives (respectively convective and diffusive) of the density ρ and order parameter f h fields. ∇ terms appear only for m = 1 (convective terms are forbidden in a nematic m = 2 or higher symmetry) and, due to the substitution of the slaved higher order term (46), are both linear and nonlinear. On the other hand, ∇∇ terms are linear to order 3 . Obviously, all terms in Eq. ( 48) need to respect the rotational invariance of the system. The familiar Ginzburg-Landau term µ h (σ, ρ) -ξ h |f h | 2 f h appearing in Eq. (48) locally describes the spontaneous breaking of a continuous symmetry, with the ordered solution

|f h | = µ h ξ h . ( 49
)
On the other hand, if n > m, the order parameter fields k < h are not zero by symmetry, and by our scaling ansatz (43) their dynamics may be relevant. Their equations need to be added to the hydrodynamic description. Moreover, the modes f j with j = 1, . . . , h are now generically coupled to the enslaved modes f h+j . Substituting the latter in the former, generates a larger number of nonlinear terms, which now also include different combinations of the fields without spatial derivatives. However, the local symmetry breaking phenomenon is still described by the Ginzburg-Landau term appearing in the f h equation.

For completeness, in the following, we derive the simplest well-behaved equations describing the physics at hand for the three basic classes introduced above.

Nematic particles with nematic interactions

In the case of nematic particles, equations since v 0 = 0 and f 1 = 0. The scaling of the different fields is given by Eq. (38). First of all, the continuity equation (31) reads in this case, given all non-zero terms are of order 3 ,

∂ρ ∂t = D 0 ∆ρ + D 1 2 Re(∇ * 2 f 2 ). ( 50
)
The equation for f 2 reads, to order 3 ,

∂f 2 ∂t = µ 2 f 2 + D 0 ∆f 2 + D 1 4 ∇ 2 ρ + [P 2 (I 2,4 + I 2,-2 ) -I 0,4 -I 0,-2 ]f * 2 f 4 . (51)
Similarly, the equation for f 4 up to order 3 is given by

∂f 4 ∂t = µ 4 f 4 + D 1 4 ∇ 2 f 2 + (P 4 I 4,2 -I 0,2 )f 2 2 . (52)
We recall that the coefficients µ k are defined in Eq. (45). Assuming µ 4 < 0, which has to be checked case by case (it is indeed true for the representative member of the class with H(∆) = ∆ 2 [27]), f 4 can be slaved to f 2 . This slaving procedure is indeed consistent with scaling considerations. Since f 4 appears in Eq. (51) only in the product f * 2 f 4 , only terms up to order 2 need to be kept in the expression of f 4 , leading to

f 4 = - 1 µ 4 (P 4 I 4,2 -I 0,2 )f 2 2 . ( 53
)
Injecting this expression of f 4 into Eq. ( 51), one finds the following closed equation for

f 2 , ∂f 2 ∂t = µ 2 f 2 + ν∆f 2 + χ∇ 2 ρ -ξ|f 2 | 2 f 2 , (54)
where the coefficients ν, χ and ξ are defined as

ν = D 0 , χ = D 1 4 (55) ξ = P 4 I 4,2 -I 0,2 µ 4 [P 2 (I 2,4 +I 2,-2 )-I 0,4 -I 0,-2 ].
Note that here, all terms allowed by symmetry4 up to order 3 have been generated.

Polar particles with ferromagnetic interactions

For polar particles, the scaling of the different quantities is given in Eq. (40). Retaining terms up to order 3 , the continuity equation reads

∂ρ ∂t + v 0 Re(∇ * f 1 ) = D 0 ∆ρ. ( 56
)
Truncating the equation for f 1 to order 3 yields

∂f 1 ∂t + v 0 2 (∇ρ + ∇ * f 2 ) = µ 1 f 1 + D 0 ∆f 1 + D 1 4 ∇ 2 f * 1 (57) + [P 1 (I 1,-1 + I 1,2 ) -I 0,-1 -I 0,2 ]f * 1 f 2 .
Since f 2 appears only in the terms ∇ * f 2 and f * 1 f 2 , only terms of order up to 2 need to be retained in the expansion of f 2 , yielding

v 0 2 ∇f 1 = µ 2 f 2 + (P 2 I 2,1 -I 0,1 )f 2 1 . (58)
Note that the time derivative has been dropped for being of order 3 . As already outlined in the active nematics case, this procedure is however only consistent when the linear coefficient µ 2 is negative, so that f 2 can be slaved to f 1 . From Eq. ( 58), f 2 can be expressed as a function of f 1 and ∇f 1 . Injecting the resulting expression for f 2 in Eq. ( 57), one eventually finds

∂f 1 ∂t = - v 0 2 ∇ρ + µ 1 f 1 -ξ|f 1 | 2 f 1 + ν∆f 1 + χ∇ 2 f * 1 + κ 1 f 1 ∇ * f 1 + κ 2 f * 1 ∇f 1 (59)
where the coefficients are given by

ν = D 0 + v 2 0 4 |µ 2 | , χ = D 1 4 κ 1 = v 0 µ 2 (P 2 I 2,1 -I 0,1 ) (60)
κ 2 = v 0 2µ 2 [P 1 (I 1,-1 + I 1,2 ) -I 0,-1 -I 0,2 ] ξ = P 2 I 2,1 -I 0,1 µ 2 [P 1 (I 1,-1 +I 1,2 )-I 0,-1 -I 0,2 ]
the coefficient µ 2 being given by Eq. ( 45) with k = 2. Note that contrary to the case of particles with continuous time ballistic motion along their polarity vector, as studied in [13,31], the presence of diffusion at the microscopic level of the dynamics generates an anisotropic diffusion term ∇ 2 f * 1 . Only the term f 1 ∇f * 1 , that would be allowed by symmetry, is not present here, and can thus not be obtained from point-like particles with (purely local) binary interactions. frame axis. Under a rotation of angle θ0 of the frame axis, the distribution f (θ) is changed into f (θ -θ0), so that in turn the Fourier coefficients f k are changed into e ikθ 0 f k . More generally, if a quantity is multiplied by a factor e isθ 0 under a rotation of axis of angle θ0, we shall call 'spin' the integer s. Hence f k has a spin s = k (note that s can be negative). As expected, the polar order parameter f1 has spin 1 and maps to a vector, while the nematic order parameter has spin 2 and maps to a tensor. The derivation operator ∇ has a spin s = 1, and complex conjugation reverses the sign of s. The total spin of a product of factors is the sum of the spin of each factor. In the evolution equation for f k , the only terms allowed by symmetry are thus those with a total spin equal to k.

Polar particles with nematic interactions

For polar particles with nematic interactions ("rods"), the scaling of the different quantities is given in Eq. (42). From this scaling, one sees that all terms in the continuity equation (31) have to be kept. After truncation to order 3 , the equation for f 1 reads, using Eq. ( 28),

∂f 1 ∂t + v 0 2 (∇ρ + ∇ * f 2 ) = µ 1 f 1 + D 0 ∆f 1 + χ∇ 2 f * 1 + [P 1 (I 1,-1 + I 1,2 ) -I 0,-1 -I 0,2 ]f * 1 f 2 + [P 1 (I 1,-2 + I 1,3 ) -I 0,-2 -I 0,3 ]f * 2 f 3 ( 61
)
with χ = D1 4 as above. Similarly, one finds for f 2

∂f 2 ∂t + v 0 2 (∇f 1 + ∇ * f 3 ) = µ 2 f 2 + D 0 ∆f 2 + χ∇ 2 ρ + (P 2 I 2,1 -I 0,1 )f 2 1 + [P 2 (I 2,-1 + I 2,3 ) -I 0,-1 -I 0,3 ]f * 1 f 3 + [P 2 (I 2,-2 + I 2,4 ) -I 0,-2 -I 0,4 ]f * 2 f 4 . (62)
One also needs to consider the equations for f 3 and f 4 . Since f 3 and f 4 appear combined with f 1 , f 2 or derivative operators, we only need to expand these two modes to order 2 , leading to

f 3 = v 0 2µ 3 ∇f 2 - 1 µ 3 [P 3 (I 3,1 +I 3,2 )-I 0,1 -I 0,2 ]f 1 f 2 (63) f 4 = - 1 µ 4 (P 4 I 4,2 -I 0,2 )f 2 2 . ( 64
)
Note that, here again, enslaving f 3 and f 4 to f 1 and f 2 is only consistent if µ 3 and µ 4 are negative. Injecting these expressions of f 3 and f 4 in Eqs. ( 61) and (62) for f 1 and f 2 , we end up with

∂f 1 ∂t = - v 0 2 (∇ρ + ∇ * f 2 ) + µ 1 f 1 + β|f 2 | 2 f 1 + ζf * 1 f 2 (65) +ν 1 ∆f 1 + χ∇ 2 f * 1 + γf * 2 ∇f 2 ∂f 2 ∂t = - v 0 2 ∇f 1 + (µ 2 -ξ|f 2 | 2 )f 2 + ωf 2 1 + τ |f 1 | 2 f 2 +ν 2 ∆f 2 + χ∇ 2 ρ + κ 1 f * 1 ∇f 2 + κ 2 ∇ * (f 1 f 2 ) ( 66
)
where the coefficients are given by

γ = v 0 2µ 3 [P 1 (I 1,-2 + I 1,3 ) -I 0,-2 -I 0,3 ] β = - 2γ v 0 [P 3 (I 3,1 + I 3,2 ) -I 0,1 -I 0,2 ] ζ = P 1 (I 1,-1 + I 1,2 ) -I 0,-1 -I 0,2 ν 1 = D 0 , χ = D 1 4
ν 2 = D 0 + v 2 0 4 |µ 3 | , ω = P 2 I 2,1 -I 0,1 (67)
κ 1 = v 0 2µ 3 [P 2 (I 2,-1 + I 2,3 ) -I 0,-1 -I 0,3 ] κ 2 = v 0 2µ 3 [P 3 (I 3,1 + I 3,2 ) -I 0,1 -I 0,2 ] ξ = P 4 I 4,2 -I 0,2 µ 4 [P 2 (I 2,-2 +I 2,4 )-I 0,-2 -I 0,4 ] τ = - 2κ 1 v 0 [P 3 (I 3,1 + I 3,2 ) -I 0,1 -I 0,2 ] .
Although Eqs. ( 65) and (66) already contain many terms, not all terms allowed by symmetry are present, especially in the equation for f 1 . Note however that if one further enslaves f 1 to f 2 , new terms will be generated in the equation for f 2 .

Discussion

ρ-dependence

In Section 3, we derived a number of transport coefficients appearing in the hydrodynamic equations for the three main classes of models. Although this was not made explicit, most of the coefficients we have derived actually depend on the local density ρ, through the linear coefficients µ k . A natural question is then that of the qualitative influence of this dependence of the transport coefficients on the density field. Indeed, it would be more convenient to have constant coefficients, depending only on microscopic parameters of the model and on the mean density. As we will see below, the ρ dependence of the linear term governing the instability of the order parameter is at the origin of the linear instability of the homogeneous ordered state ultimately leading to the emergence of the inhomogeneous solutions in the transition region. It is thus essential and must be kept. The density-dependence of the other linear terms µ k f k , in contrast, is only felt indirectly in the definition of the coefficients of the other terms. It is tempting to expand them around the mean density ρ 0 , setting ρ = ρ 0 +δρ. For terms of order 3 , corrections in δρ of the coefficients lead to terms of order 4 or higher, and can thus a priori be discarded.

The question however becomes more subtle if one wishes to assess the qualitative validity of the derived hydrodynamic equations further away from the threshold. Strictly speaking, these equations are of course not expected to be valid away from the threshold, which could in a sense close the debate. Yet, the example of Ginzburg-Landau equations shows that in many cases the qualitative behavior away from the linear instability threshold remains well described by the obtained equations. Unfortunately, there is no well-defined general procedure to discriminate between different formulations of the same terms, if these formulations become equivalent close The color codes for the angle between the most unstable wavevector and the direction of order. The line σ || indicates the limit above which longitudinal (i.e., with wavevector parallel to order) perturbations are unstable (the upper limit being σt). The line σ ⊥ indicates the limit below which transverse (i.e., with wavevector perpendicular to order) perturbations are unstable.

to threshold. This is the case in particular for the density-expansion of the transport coefficients, and one has to rely on numerical integrations of the hydrodynamic equations to check the effect of the ρ dependence of the different terms.

For the coefficient ξ of the cubic term, some theoretical justifications of the necessity to keep the full ρ-dependence of ξ (as opposed to a dependence on the mean density ρ 0 ) close to threshold can also be provided. For instance, the role of the density-dependence of ξ appears in the dynamical system analysis of the shape of non-linear patterns arising for polar particles interacting ferromagnetically [33]. Note also that the dependence of ξ over the average density ρ 0 , that we keep when doing the expansion of the density field around the mean density value, is crucial for the value f (0) k /ρ = ρ -1 µ k /ξ of the polarity (in the homogeneous phase) to increase with mean density and to saturate to a finite value at large mean density. Finally, when studying the linear instability of the homogeneous ordered state, one finds that the term in dξ/dρ (obtained from a density expansion of ξ(ρ)) is not negligible (with respect to other terms that have to be kept) in the computation of the growth rates, and may lead to significant difference in the linear stability diagrams. We illustrate this in the next section.

σ 2 5 4 3 2 1 0 ρ σ t σ s σ || σ u σ ⊥ Fig. 2.
Phase diagram of the linear stability analysis of the homogeneous solutions to the hydrodynamic equations with and without the ρ-dependence of the transport coefficients of the nonlinear terms (see text). Left column: polar particles aligning nematically, without and with ρ dependence (top and bottom panel respectively). Top right: nematic particles aligning nematically (no discernible influence of ρ dependence). Bottom right: polar particle aligning ferromagnetically with ρ dependence (compare with Fig. 1). No positional diffusion (D0 = D1 = 0). Same colormap, same legends as in Fig. 1.

Linear stability of homogeneous solutions

Denoting as f h the relevant order parameter, either polar (h = 1) or nematic (h = 2), an ordered solution generically exists in the region where the disordered state becomes unstable (µ h > 0). We denote this line σ t (ρ) in the basic noise-density phase diagram in which we performed a systematic study of the linear instability of the homogeneous ordered solution, not restricting ourselves to modes longitudinal or transversal to order. Results below illustrate this analysis.

In previous studies of the cases without positional diffusion, which in the case of polar particles yield slightly simpler equations (anisotropic diffusion terms like ∇ 2 f * 1 are absent due to a continuous time microscopic dynamics), a generic qualitative scenario was found in all cases, irrespective of the symmetries of the particles and of the interactions [13,25,27]. We first illustrate it in the case of polar particles aligning ferromagnetically where only the ρ dependence of the linear term has been kept (Fig. 1).

General scenario

Close to the line σ t corresponding to the onset of order, the homogeneous ordered state is unstable with respect to long-wavelength perturbations, leading to the formation of non-linear patterns (see next section). This instability is generically related to the density-dependence of the linear coefficient µ h , the latter generically being an increasing function of the local density ρ (only in the metric-free case, where µ h becomes independent of the density, does one recover a stable homogeneous phase close to the transition line [31]).

This instability region bordering the line σ t (ρ) closes at the line σ s (ρ) below which the homogeneous ordered solution is linearly stable (Fig. 1). However, there exists a second instability region below a line σ u (ρ). Preliminary studies of the stability of the Boltzmann equation in the case of ferromagnetic interactions indicate that this second instability region does not exist at the Boltzmann level. We thus conclude that this second instability region is an artifact of the truncation procedure.

For the unstable region between σ t and σ s , the unstable modes have long wavelengths and are oriented mostly along order, but not perfectly so. The line marking the limit of stability of purely longitudinal modes is marked in magenta. The second, "spurious" instability region below σ u has large growth rates and unstable wavevectors mostly transversal to order, although again not perfectly so (the line marking the onset of purely transversal unstable modes is marked in green inside this region).

For the case of particles aligning nematically, the unstable modes are almost transversal in the first region. The second region does not exist for active nematics, while its most unstable mode is at k = 0 for the self-propelled rods case (Fig. 2, top panels).

Influence of the ρ-dependence of transport coefficients

We now illustrate the influence of keeping the local density dependence in the transport coefficients, and notably in the cubic nonlinear term. For polar particles aligning ferromagnetically, this widens the region of stability of the homogeneous ordered state, and reduces the extent of the second instability region (compare Fig. 1 and Fig. 2, bottom right).

For the self-propelled rods case, keeping the ρ dependence has a dramatic effect: the second region now "invades" the first one near the origin, becoming dominant (compare top left and bottom left panels in Fig. 2).

For active nematics, the ρ dependence has a very weak effect, not discernable on the scale of Fig. 2 (see top right panel).

Influence of positional diffusion

We finally illustrate the influence of positional diffusion (coefficients D 0 and D 1 ) on the linear stability diagram. As expected, the general effect of diffusion is to stabilize the homogeneous ordered solution. The extent of the first instability region (between lines σ t and σ s ) varies only weakly, whereas the second, spurious, region can be made very small under the influence of strong diffusion. We now illustrate these findings for the case of polar particles aligning ferromagnetically.

A positive D 0 shrinks a bit the first instability region while making its unstable modes more clearly longitudinal (see Fig. 3, top right, where D 1 = 0). For any fixed positive D 0 , the extent of the second region varies greatly with D 1 : it can be made smaller for positive D 1 values (Fig. 3, bottom left), and larger for negative D 1 values (Fig. 3, top left). We find that the size of the second region, as measured by the location of the point ρ zero u where σ u intersects the horizontal axis, actually scales exactly like 1/D 0 (Fig. 3, bottom right).

Inhomogeneous solutions and dynamics

The linear stability analysis recalled above revealed, for metric models, the existence of a region bordering the linear onset of order inside which the homogeneous ordered solution is unstable. Our hydrodynamic equations then exhibit nonlinear, inhomogeneous, phase-separated solutions, which are faithful to those observed in microscopic simulations. In our opinion, this makes a strong case in favor of our approach, which appears to yield an overall very good qualitative agreement with relevant active particle models.

The inhomogeneous solutions take different form depending on the case considered: in the polar case, the hydrodynamic equations typically show trains of propagating "bands", as explained recently in [13,33] and observed in microscopic models. The stability and selection mechanisms at work in these smectic arrangements is still under study. In the active nematics and self-propelled rods cases, a dominating solution was observed and found analytically, which takes the form of a stationary ordered domain elongated along nematic order, and occupying a well-defined fraction of space [27,25]. For active nematics, it was recently shown that this band solution is actually itself always linearly unstable, and leads to a large-scale spatiotemporally chaotic regime [32].

Even though the phase-separated solutions take different forms depending on the case considered, we find that they exist, and are observed, in a parameter region wider than the linear instability domain of the homogeneous state. This important observation immediately implies that there exist two coexistence regions (in the sense that two different solutions may exist for the same values of the parameters) at the level of the deterministic hydrodynamic equations. In particular the nonlinear phaseseparated solutions coexist with the disordered phase, which explains the discontinuous character of the transition, and points to the "irrelevance" of linear stability thresholds in the presence of fluctuations. Indeed, they are much like spinodal lines in usual phase separation phenomena, i.e. they are not the transition points observed in systems with fluctuations.

On the Boltzmann-Ginzburg-Landau approximation

Boltzmann equation

The Boltzmann approach is based on the binary collision approximation and the molecular chaos hypothesis. Under these conditions, the Boltzmann equation [START_REF] Huang | Statistical Mechanics[END_REF] faithfully describes our microscopic dynamics on timescales τ B much larger than the microscopic dissipative scale ∆t and the collision timescale τ coll , and on lengthscales B much larger than the different microscopic scales, namely the elementary displacement d 0 = v s ∆t (where v s is the characteristic speed of particles), the interaction range r 0 , and the typical distance between particles ∼ 1/ √ ρ 0 . Metric-free models are briefly discussed in Section 2.4, so that here we will only discuss metric models.

In metric models, binary collisions dominate the dynamics when (i) the interaction range r 0 is much smaller then the mean particle interdistance ∼ 1/ √ ρ 0 (in two spatial dimensions) and (ii) the intercollision time τ free is much larger than the collision timescale τ coll .

In self-propelled systems, dominated by drift, one has

τ (b) free ≈ 1 v 0 √ ρ 0 (68)
while in purely diffusive systems we get

τ (d) free ≈ 1 Dρ 0 (69)
where v 0 is the drift speed and D a positional diffusion constant (either isotropic or anisotropic). The collision time, on the other hand, is given by the time needed by two interacting particles to separate; since interacting particles are aligned up to some noise, it is always diffusion dominated

τ coll ≈ r 2 0 D . (70)
In general the condition τ coll τ free (as well as r 0 1/ √ ρ 0 ) can be fulfilled at low densities in both ballistic and diffusive systems, thus justifying the binary collision approximation in drift-diffusion dry active systems. It should be noted, however, that in a ballistic system like the original VM, positional diffusion is strictly zero, and interacting particles may be only separated by angular noise, with an effective diffusion constant

D eff ∼ d 2 0 η 2 ∆t . (71)
Since in metric systems the critical noise amplitude η c corresponding to the onset of order decreases with density, this casts some doubts on the applicability range of the binary collision approximation for the original VM in the absence of positional diffusion. In particular, using a mean-field estimate [23] of the critical noise η c ∼ d 0 √ ρ 0 , Eq. (71) becomes

τ (V M ) coll ≈ ∆t r 2 0 d 4 0 ρ 0 (72)
so that, with v 0 = d 0 /∆t, the condition τ coll τ free would imply

r 2 0 d 3 0 √ ρ 0 1 (73)
to be satisfied simultaneously with condition (i)

r 0 1 √ ρ 0 . (74)
Eqs. ( 73)-(74) can only be satisfied simultaneously if (r 0 /d 0 ) 3 1, so that one expects that for the classic VM, as well as for other ballistic models, higher-order collisions will become relevant even at low densities for moderate to low speeds.

The molecular chaos hypothesis, which allows one to decompose many-particle distributions into products of single-particle ones, on the other hand, is a rather delicate issue. It was investigated in Refs. [34], using a binary collision VM, through the chaos propagation property: if the latter holds, initially decorrelated particles will stay decorrelated at future times, and thus the molecular chaos assumption will be justified. The analysis of [34] is however inconclusive. While at the kinetic timescale chaos propagation holds in the thermodynamic limit, longer time behavior remains out of reach of exact methods (while preliminary numerical simulations seem to indicate a violation of the molecular chaos hypothesis). See also [38] for related discussions using a numerical analysis.

Convergence issues in the hydrodynamic expansion

The second approximation involved in the BGL approach is the small expansion used to close the hierarchy (28) in order to derive hydrodynamic equations. While this expansion is 'exact' at the onset of order, one can wonder whether the obtained hydrodynamic equations remain valid away from the transition line, in particular for what regards linear stability results.

For polar particles, in the explicit examples where the stability analysis of the ordered phase has been performed in detail, we have found an additional instability deeper in the ordered phase, further away from the linear instability threshold of the disordered phase (see sect. 4.2.1) [31,25]. Such an additional instability was however not found for active nematics [27].

This instability is considered as spurious in the sense that it is not observed in the corresponding microscopic models. It typically occurs with very large growth rates, raising doubts as to its physical meaning. It is likely to be an artifact of the truncation procedure. Though we have not performed a detailed analysis of all three possible classes, results obtained for polar particles with ferromagnetic interactions indicate that this spurious instability is not present at the Boltzmann level [35].

The convergence of the -expansion of the Boltzmann equation at a finite distance from the instability threshold is a complicated issue. A first guess is that going to higher orders in the -expansion could improve the description further away from threshold, e.g. by removing, or at least shrinking, the region of the additional, "spurious" instability.

We have derived the hydrodynamic equations up to order 7 in for the canonical case of polar particles aligning ferromagnetically (Vicsek model) [START_REF] Peshkov | Boltzmann-Ginzburg-Landau approach for simple models of active matter[END_REF]. It is relatively easy to see that only odd orders in need to be considered, and that even orders yield intrinsically unbalanced equations. For example, at order 4, one keeps the equation for f 2 , while f 3 is slaved and expressed as a function of f 1 and f 2 , yielding equations of the form:

∂ t f 1 = a 1 f 1 + a 21 f 2 f * 1 (75) ∂ t f 2 = b 2 f 2 + b 11 f 2 1 + b 211 f 2 |f 1 | 2 (76)
(for clarity, terms involving space derivatives are not shown). The last term of the second equation is the only new, fourth-order one. It can be shown to render the homogeneous ordered solution linearly unstable everywhere. At order 5, we obtain 3 equations of the form:

∂ t f 1 = a 1 f 1 + a 21 f 2 f * 1 + a 23 f * 2 f 3 (77) ∂ t f 2 = b 2 f 2 + b 11 f 2 1 + b 211 f 2 |f 1 | 2 (78) ∂ t f 3 = c 3 f 3 + c 12 f 1 f 2 + c 113 |f 1 | 2 f 3 + c 122 f * 1 f 2 2 (79)
in which several order 5 terms appear. The homogeneous ordered solution now has to be determined numerically. We find that it exists only for a limited range of noise values, whereas no such limitation occurs at order 3. Thus order 5 is "worse" than order 3. At order 7, we now have 5 equations of the form:

∂ t f 1 = a 1 f 1 + a 21 f 2 f * 1 + a 23 f * 2 f 3 + a 34 f * 3 f 4 (80) ∂ t f 2 = b 2 f 2 + b 11 f 2 1 + b 211 f 2 |f 1 | 2 + b 24 f * 2 f 4 (81) ∂ t f 3 = c 3 f 3 + c 12 f 1 f 2 + c 14 f * 1 f 4 + c 25 f * 2 f 5 (82) ∂ t f 4 = d 4 f 4 + d 22 f 2 2 + d 13 f 1 f 3 + d 15 f * 1 f 5 (83) ∂ t f 5 = e 5 f 5 + e 14 f 1 f 4 + e 23 f 2 f 3 + e 115 |f 1 | 2 f 5 + e 133 f * 1 f 2 3 + e 124 f * 1 f 2 f 4 . (84)
Now the homogeneous solution exists over the whole noise range again. The maximal value of polarization (|f 1 |/ρ) reached at zero noise is higher than at order 3, which is an improvement. However, the linear stability study of the homogeneous ordered solution shows that it is linearly unstable everywhere: even at order 7, the situation is less satisfactory than at order 3.

Even though these results would need to be completed by a similar study performed on the altogether slightly better-behaved case of active nematics, they indicate already that the BGL approach only makes sense at the usual, third-order, Ginzburg-Landau level. Exploring higher orders seems to only reveal that an angular Fourier expansion is not suited far from the onset of order, calling for more sophisticated expansion schemes, possibly along the lines of the "reductive perturbation theory" approach introduced by Chen et al. [START_REF] Chen | [END_REF].

Extended Particles

Our framework here has been restricted to interactions that respect the particle exchange symmetry. However, it has been pointed out that 'more realistic' extended particles may show a more complicated collision dynamics [39]. In particular, collisions may locally violate the particle exchange symmetry which, in a Boltzmann approach, can only be restored in a statistical sense. In particular, Ref. [39] computes the effective Boltzmann level collision dynamics for two different extended particles models: self-propelled hard rods and a bead-spring model. In both cases, it is claimed that binary collisions do not lead to a positive linear coefficient µ 1 for the polar order parameter. Contrary to the conclusions of Ref. [39], we do not find these results significative to the point of constituting "a critical assessment" of the BGL approach, as there is no indication that the two microscopic models studied give rise to polar order on the mesoscopic scale. On the other hand, our discussion in Sec. 3.2.2 clearly shows that at sufficiently low noise, any interaction rule being indeed an alignment rule results, at the Boltzmann level, in the linear instability of the disordered solution towards the order parameter associated to the interaction symmetry.

Comparison with other approaches

Other statistical physics approaches have also been proposed to derive hydrodynamic equations for active particles, mainly the mean-field Fokker-Planck equation approach [14,40,41,42] and the Enskog-like approach [16].

Mean-field Fokker-Planck equation

In the mean-field Fokker-Planck equation case, one starts from interacting active particles described by a set of Langevin equations, from which a N -body Fokker-Planck equation is derived. This N -body equation is then treated in a mean-field way, meaning that the force experienced by a given particle is replaced by a local statistical average of the force. The N -body linear Fokker-Planck equation is then turned into a one-body nonlinear Fokker-Planck equation, under an assumption of factorization of the N -body distribution as a product of one-body distributions [14,40,41,42]. If the active particles are described by underdamped Langevin equations, some approximation scheme is further needed to eliminate the velocity degrees of freedom, then reducing the mean-field Fokker-Planck to a Smoluchowski equation [14].

Although the starting points are significantly different, the resulting Smoluchowski (or Fokker-Planck in the overdamped case) equation bears some formal similarities with the Boltzmann equation. It also describes the evolution of the one-body distribution f (r, θ, t), and takes the form of a bilinear integral equation. As a result, the angular Fourier expansion takes a form close to that obtained from the Boltzmann equation, and the same approximation scheme can be used close to the instability threshold of the disordered solution (though most papers following this approach use truncations at some lower order than we do in the BGL approach).

For point-like self-propelled polar particles with overdamped Langevin dynamics and ferromagnetic interactions, the resulting hydrodynamic equations are very similar to the ones obtained from the Boltzmann equation [42].

For rods, the hydrodynamic equations derived from a Smoluchowski equation in [14] are different from the ones we obtain here for polar particles with nematic interactions, for two main reasons. First, in [14], only f 1 and f 2 have been retained in the truncation, while we have kept modes up to f 4 , and used f 3 and f 4 , once enslaved to f 1 and f 2 , to close the hierarchy of equations. We thus have more terms (and in particular more non-linear terms, including the cubic term saturating the instability of the state f 2 = 0) in our equations. Second, the spatial extension of the rods has been explicitly taken into account in [14], leading to the emergence of new terms coupling the polar order parameter f 1 to its gradient, in the equation for f 1 , while we have considered point-like particles in our Boltzmann equation. Additional terms are thus obtained in [14] with respect to our results. More generally, the Fokker-Planck (or Smoluchowski) equations constitute a well-suited framework to include positionand/or orientation-dependent forces, as done in [14] for extended rods, or in [43] for forces depending on the relative orientation of particles. Although we have not done it explicitly up to now, such extensions should also be possible in the context of the Boltzmann equation, and are expected to yield the same type of additional terms into the resulting hydrodynamic equations.

On the other hand, one of the drawbacks of the Fokker-Planck approach is its sensitivity to some details of the dynamics. Indeed, in the case of point-like selfpropelled particles with overdamped Langevin dynamics considered in [42], it can be seen rather easily that normalizing or not the local force (exerted by neighboring particles) by the number of neighbors results in significantly different macroscopic dynamics at the level of the hydrodynamic equations, since the linear term in the equation for the polar order parameter becomes independent of the density if the force is normalized -recall that, in this case, the basic instability of the homogeneous ordered state is suppressed (as in metric-free models).

Enskog equation approach

An alternative method has been proposed, based on an Enskog-like equation, in order to deal with multiple collisions in the semi-dilute regime [16]. Its main advantage, with respect to the Boltzmann approach, is that it is not limited to binary interactions; all multiple interactions are taken into account. However, it also relies, like the other methods, on a factorization assumption of the N -body distribution fonction, or in other words, on a 'molecular chaos' assumption. From the Enskog equation, an angular Fourier expansion can also be performed. In principle, an -expansion close to the instability threshold could be performed, in the same way as for the Boltzmann equation. The derivation made in [16] is presented in a slightly different way. In particular the truncation and closure performed leads to hydrodynamic equations with many more (higher-order) terms than those retained in the BGL approach. At the level of the linear stability of the homogeneous solutions, these equations yield more quantitatively accurate thresholds than BGL, as compared to numerical simulations of the VM. This is satisfactory in principle, but in view of the generic presence of the long-wavelength instability of the ordered state leading to bands, and of the existence and stability domains of these nonlinear solutions extending beyond the linear thresholds, the transition lines observed in the presence of fluctuations (and firstly in the microscopic VM) cannot be the linear thresholds found in the (deterministic) hydrodynamic equations. Moreover, at the nonlinear level, the many terms retained in [16] seem to lead to unbounded solutions: the obtained hydrodynamic equations are ill-behaved.

Conclusion

In this paper, we have presented a generic Boltzmann-Ginzburg-Landau framework allowing us to derive, in a controlled way, continuous equations describing the largescale behavior of assemblies of point-like active particles interacting through polar or nematic alignment rules. Also varying the symmetry properties (polar or nematic) of the interactions between particles, we have identified three main classes of dry active systems, and derived for each class the corresponding continuous equations for the relevant order parameters. One of the main advantages of our approach is that the truncation and closure of the infinite hierarchy is well-controlled, using standard scaling arguments close to the linear instability threshold of the disordered state. As often the case for Ginzburg-Landau-type equations, the resulting minimal equations actually also describe, in a qualitative manner, the behavior of the system further away from the linear instability threshold, when nonlinear structures appear and possibly have a complicated spatiotemporal dynamics. We thus believe that on the one hand, our truncation method to derive continuous hydrodynamic equations is quite general, and can be applied to other starting points like Fokker-Planck or Smoluchowski equations [42,43], and on the other hand that the resulting equations are, for each class, the minimal well-behaved equations accounting for the relevant phenomenology.

Fig. 1 .

 1 Fig.1. Phase diagram of the homogeneous solutions of the hydrodynamic equations for polar particles with ferromagnetic interactions in the noise-density plane, in the absence of microscopic positional diffusion (D0 = D1 = 0), and without the ρ dependence of the transport coefficients of nonlinear terms. The full black line σt indicates the stability limit of the homogeneous disordered state (upper white area). The lower white area corresponds to the stability domain of the homogeneous ordered state, delimited from above by the line σs and from below by the line σu. In colored regions, the homogeneous ordered state is unstable. The color codes for the angle between the most unstable wavevector and the direction of order. The line σ || indicates the limit above which longitudinal (i.e., with wavevector parallel to order) perturbations are unstable (the upper limit being σt). The line σ ⊥ indicates the limit below which transverse (i.e., with wavevector perpendicular to order) perturbations are unstable.

Fig. 3 .

 3 Fig. 3. Phase diagram of the linear stability analysis of the homogeneous solutions to the hydrodynamic equations for the case of polar particles aligning ferromagnetically in the presence of positional diffusion: D1 = -0.5 (top left), D1 = 0 (top right), D1 = 0.5 (bottom left); D0 = 0.25 in all three cases, same colormap as in Fig. 1. Bottom right: intersection ρ zero u of the line σu with the horizontal axis, as a function of D0 for different cases (see legends), showing that ρ zero u ∝ 1/D0.

However, recent experimental results[17] suggest this may not be the case in starling flocks.

Note that strictly speaking, the Boltzmann equation describes a version of the microscopic model in which interactions takes place only once during a collision, when the interparticle distance reaches the interaction range. No further interaction occurs during the phase when particles are closer than the interaction range. A new interaction takes place when particles collide again.

Note that, crucially, density dependence cancels out for metric-free systems.

It is useful to specify how the "symmetry" of terms is characterized. In the same way as for vectors and tensors, it describes how a given quantity transforms under a rotation of the

Acknowledgements

We thank the Max Planck Institute for the Physics of Complex Systems, Dresden, for providing the framework of the Advanced Study Group "Statistical Physics of Collective Motion" within which part of this work was conducted. We also thank the Kavli Institute for Theoretical Physics in Santa Barbara where this work was completed, within the framework of the program "Active Matter: Cytoskeleton, Cells, Tissues and Flocks". We are also grateful to A. Baskaran and M.C. Marchetti for useful discussions. FG acknowledges support from grants EPSRC First Grant EP/K018450/1 and MC Career Integration Grant PCIG13-GA-2013-618399.