Phase transition and diffusion among socially interacting self-propelled agents

We consider a hydrodynamic model of swarming behavior derived from the kinetic description of a particle system combining a noisy Cucker-Smale consensus force and self-propulsion. In the large self-propulsion force limit, we provide evidence of a phase transition from disordered to ordered motion which manifests itself as a change of type of the limit model (from hyperbolic to diffusive) at the crossing of a critical noise intensity. In the hyperbolic regime, the resulting model, referred to as the `Self-Organized Hydrodynamics (SOH)', consists of a system of compressible Euler equations with a speed constraint. We show that the range of SOH models obtained by this limit is restricted. To waive this restriction, we compute the Navier-Stokes diffusive corrections to the hydrodynamic model. Adding these diffusive corrections, the limit of a large propulsion force yields unrestricted SOH models and offers an alternative to the derivation of the SOH using kinetic models with speed constraints.


Introduction
There is a considerable literature devoted to the observation and understanding of systems of swarming agents. Examples of such systems in nature are fish schools [6,39], bird flocks [2,40], insect swarms [14,18] or migrating cell assemblies [55] (see also the reviews [19,54]). Some simple inanimate physical systems also exhibit collective behavior [29,47]. Many of the models proposed in the literature are 'Individual-Based Models (IBM)'. They consist in following the dynamics of the agents and their interactions over the course of time. The 'three zone model' of Aoki [1,20] postulates that interactions obey longrange attraction, short-range repulsion and medium-range alignment. Vicsek et al. [53] have proposed a simplified version of this model where particles move at a constant speed and interact through alignment only. In spite of its simplicity, the Vicsek model exhibits complex features which have triggered a large literature (see [54] for a review and references). On the other hand, the Cucker-Smale model [22] is based on a large-scale velocity consensus formation and does not impose any constraint on the particle speed. The Cucker-Smale model has triggered considerable mathematical activity [15,21,35,36,37,46]. Many other kinds of IBM's of collective motion can be found and it is impossible to cite them all (see e.g. [16,17,30,42] and the review [54]). Additionally, comparisons of models with data can be found e.g. in [4,5,38].
IBM's are very successful but become computationally intensive for large systems. For this reason, macroscopic models of fluid type have been proposed in the literature. Macroscopic models of collective motion have been derived from heuristic rules and symmetry considerations in [48,49]. The rigorous derivation of continuum models usually starts from a statistical version of the IBM, the so-called kinetic model. Kinetic models of collective motion have been derived in [11,12,15] for various versions of the Cucker-Smale and Vicsek models. The convergence of the kinetic Cucker-Smale model to the kinetic Vicsek model is shown in [10]. In [7,8], a Boltzmann type kinetic model has been proposed for binary collision processes which mimics the Vicsek alignment dynamics. In the same spirit, a Boltzmann-Povzner type approach which mimics the Cucker-Smale process and its fluid limit has been developed in [32]. In [7,8], a hydrodynamic model for the binary version of the Vicsek interaction is derived from the kinetic model under an assumption of weak anisotropy of the velocity distribution function. In [45] a direct passage from the Vicsek IBM to a fluid model is attempted. The first derivation of hydrodynamic-like equations from the mean-field kinetic version of the Vicsek model has been performed in [26]. Further elaboration of the model can be found in [24,25,33,34]. Diffusive corrections to the model of [26] have been derived in [28] and bear analogies with the model proposed in [43]. Other kinds of macroscopic models can be found in [9,31,41,50,51].
The aim of the present work is twofold. The first objective is to give evidence of a phase transition from disordered to ordered motion in a hydrodynamic model of socially interacting agents with self-propulsion. Specifically, we want to emphasize the role of the self-propulsion in the emergence of the phase transition. Such evidence has been given for the first time in [48]. However, the model of [48] is based on analogy with the Vicsek IBM, not on an actual derivation from it. In [48], the techniques used to show the emergence of a phase transition are rather complex: they are based on stability analysis in the linear case and renormalization group theory in the nonlinear case. In the present work, the model we investigate is derived from a simple IBM of collective motion combining a noisy Cucker-Smale consensus force and self-propulsion. The phase transition appears simply when analyzing the behavior of the model in the limit of a large self-propulsion force. It manifests itself as a change in the type of limit model at the crossing of a critical noise intensity. Above the critical noise, the limit model is of diffusion type, while below the critical noise, it is of hyperbolic type. To the best of our knowledge, the present work is the first instance where phase transitions in particle swarms have been evidenced in this way from a hydrodynamic model. A similar approach, but at the level of the kinetic model, can be found in [24,34]. The second goal of this work is to discuss the relative merits of the Cucker-Smale and the Vicsek models for the derivation of hydrodynamic models of particle swarms. As mentioned above, the Vicsek kinetic model imposes that kinetic velocities be of constant norm, while no such constraint exists in the Cucker-Smale model. Instead, a self-propulsion force is imposed to force the particle speed to stay close to a 'comfort' velocity. In [10], it is shown that the Cucker-Smale model relaxes to the Vicsek model when the intensity of the self-propulsion force tends to infinity. The derivation of hydrodynamic models from the Vicsek kinetic model is considerably complexified by the velocity norm constraint. Indeed, momentum conservation is lost and the use of conserved quantities (or collision invariants), which is the cornerstone of the derivation of hydrodynamic models, cannot be implemented. In [26], this problem has been overcome by the introduction of a new concept of 'Generalized Collision Invariant'. But if the hydrodynamic limit could be performed equivalently on the Cucker-Smale model, these unpleasant technicalities would be proven unnecessary.
Unfortunately, performing the hydrodynamic limit on the Cucker-Smale model and then letting the self-propulsion force tend to infinity is not equivalent to performing the hydrodynamic limit on the Vicsek model. This is the second main result of the present work. Indeed, the type of the resulting model is the same, but the coefficients of the model are not the same. Specifically, the limit model has the form of a system of isothermal compressible Euler equations for the swarm density ρ and the mean velocity direction ω (also referred to as the polarization field, see e.g. [43]). The velocity direction ω is a vector of norm one. To maintain this constraint, the model includes some non-conservative terms. Additionally, the convective derivatives involved in the mass and momentum transport are not the same, a signature of a loss of Galilean invariance (see e.g. [52]). This model has been referred to in [25] as the 'Self-Organized Hydrodynamic (SOH)' model. In the present work, we show that the SOH model derived from the Cucker-Smale kinetic equation necessarily involves the same convective derivatives in the density and in the velocity equations. Therefore, with the Cucker-Smale model, we cannot access the whole range of possible hydrodynamic limits that we can access with the Vicsek model. With the Cucker-Smale model, we only get a sub-class of these models, which limits its practical applications: having different convective derivatives for ρ and ω increases the likelihood of correctly reproducing emergent phenomena in swarms, such as cluster formation, waves, etc.
This restriction, which is a significant disadvantage of the Cucker-Smale approach, can be weakened, at least partially. Indeed, in the last part of the present work, we include small diffusive corrections to the hydrodynamic limit of the Cucker-Smale model, by means of a Chapman-Enskog method (see e.g. [23] for a review). If the self-propulsion is taken to infinity in the resulting compressible Navier-Stokes system, then a more general SOH model, in particular with possibly different convective derivatives for ρ and ω, can be derived. This approach is limited by the necessity to keep the diffusive corrections small; this limits the range of the coefficients of the SOH model which can be obtained. This paper shows that most SOH models can be realized in a fairly general context as hydrodynamic limits of either Cucker-Smale or Vicsek kinetic models. For these reasons, we conclude that the two approaches are somehow equivalent in the amount of technical work: while the Vicsek model made use of generalized collision invariants, starting from the Cucker-Smale model necessitates dealing with the complex diffusion terms.
The paper is organized as follows. The problem is set up in Section 2. Then, some functional properties of the operators are recalled in Section 3. The hydrodynamic and large self-propulsion limits are derived in Section 4. The diffusive corrections are dealt with in Section 5. Finally, a conclusion and some perspectives are drawn in Section 6. Three appendices collect some of the more technical proofs.
2 Setting of the problem 2.1 Velocity consensus in self-propelled agent systems We consider a system of agents with positions x i (t) ∈ R d and velocity v i (t) ∈ R d , where d is the system dimension (in practice equal to 2 or 3), t ≥ 0 is the time and i ∈ {1, . . . , N} is the agent's label. These agents are subject to a self-propulsion force which tends to restore a comfort velocity a > 0, and to a social force which drives them to the average velocity of the neighboring agents. Addtionally, they are subject to random velocity fluctuations which account for potential misperceptions and their propensity to leave the swarm and explore a new environment. The equations of motion are given bẏ . (2. 2) The social force F i is written in (2.2) as a relaxation force towards the average velocitȳ v i in the neighborhood of particle i. The relaxation rate is σ −1 (in other words, σ is the typical time needed for agent i to align with the velocity of his neighbors). The kernel K, supposed spherically symmetric for the sake of simplicity, describes how the various partner velocities v j are combined according to the distance of j to i. For instance, if K is the indicator function of the ball of radius R, it means that the agents adopt the mean velocity of the other agents within a distance R. The second term in the expression of dv i in (2.1) is the self-propulsion force. It takes the form of a relaxation term driving |v i | towards a at rate τ −1 . In other words, it takes a time τ for the velocity |v i | to relax to the comfort speed a. Finally, the last term is the velocity fluctuation term, where B i t are independent normalized Brownian processes and D > 0 is the diffusion coefficient. The force F i has been previously proposed by Cucker and Smale [22] as a model for consensus formation in particle swarms. A noisy version of the Cucker-Smale model is proposed in [21].
In the large particle limit N → ∞, by adapting the arguments in [11], the empirical measure of the system , v i (t)), can be approximated by a continuous distribution function f (x, v, t). It solves the following Fokker-Planck equation: The left-hand side expresses particle displacement at velocity v. The right-hand side consists of three terms. The first one is the consensus force. The second term is the self-propulsion force. The last term takes into account the random velocity fluctuations.

Scaling
In order to understand the roles of the various terms, it is useful to introduce dimensionless quantities. We set x 0 and t 0 to be space and time units and deduce units of velocity v 0 = x 0 /t 0 and force F 0 = x 0 /t 2 0 . We assume that the range of the interaction kernel K is R, meaning that we can write K(|x|) =K(|x|/R) withK having second moment of order 1 (i.e. K (|x|)|x| 2 dx = O(1); we assume thatK is normalized to 1, i.e. K (|x|) dx = 1). We now introduce dimensionless variablesx = x/x 0 ,t = t/t 0 ,ṽ = v/v 0 and the change of variablesf (x,ṽ,t) = . Finally, we introduce the dimensionless parameters:R In this new system of coordinates, the system is written: where we have dropped the tildes for the sake of clarity. Now, by fixing the relations between the five dimensionless parameters (2.5), we define the regime we are interested in. We suppose that the diffusion and social forces are simultaneously large, while the range of the social force tends to zero. The parameters of the self-propulsion are kept of order 1. More specifically, we let ε ≪ 1 be a small parameter and we assume that (parameters of the social force are order unity). In order to highlight these scaling assumptions, we define constants D ♯ , σ ♯ , R ♯ , which are all O(1) and such that Then, with these new notations, and dropping all 'hats' and 'sharps', we get the following scaled system: (2.10) We will investigate the limit as ε → 0 of this system, while all other parameters (i.e. τ , a, D, σ and R) are kept fixed. Hence, we highlight the dependence of f upon ε.
We can simplify the problem by using Taylor's expansion. At leading order, we find: , which leads to: (2.12) We will drop the O(ε 2 ) terms, since, as a first step, we consider only the leading and first order terms. Then, problem (2.12) can be written as with the collision operator Q(f ): Some remarks concerning scaling (2.8) can be made. The diffusion and social forces are supposed of the same order of magnitude and much larger than all other forces. They counterbalance each other. Indeed, the social force makes the agents adopt the same velocity while diffusion tends to spread the velocities out. This balance results in a Maxwellian velocity profile (i.e. Gaussian in velocity space), as shown later on. The choice which is made here is to assume that the self-propulsion force is weaker. Another choice would have been to make the self-propulsion force as large as the social force and the diffusion. In this case, the balance would involve three different effects and would result in more complicated equilibria. This investigation is in progress [3]. The interaction range is supposed to tend to zero like the inverse of the interaction rate. It is no surprise that, at leading order, only spatially local interaction terms remain (which can be seen in (2.14) by the replacement of the non-local average velocityv ε f by the local mean velocity u f ). Again, other choices can be made. In [25], in the case of the Vicsek model (which is the limit of (2.13) when τ → 0), it is shown that the different choiceR = O( √ ε) leads to a different macroscopic limit when ε → 0. This choice takes better care of the non-local character of the interaction. It will be investigated in future work. Our plan is now to investigate the hydrodynamic limit ε → 0 in this model. To this end, we first examine the properties of the collision operator Q.

Properties of Q
When ε → 0 in (2.13), f ε formally converges to an element of the null-space of Q, i.e. a function f such that Q(f ) = 0. In this limit, the dynamics are characterized by the projection of the left-hand side of (2.13) onto the space orthogonal to the range of Q. This space is spanned by the so-called 'collision invariants'. In this section, we successively determine the null-space of Q and its collision invariants.

Null-Space
We first define the Maxwellian with mean velocity u ∈ R d and temperature T = σD > 0 as follows: Note that M u satisfies M u dv = 1 and M u v dv = u.
To proceed, we need to determine a functional setting. Let u ∈ R d and define a weighted L 2 -space H u such that with the associated norms: and inner products (·, ·) Hu , (·, ·) Vu , and ((·, ·)) Vu respectively.
Lemma 3.1. (i)The operator Q given by (2.14) can be reformulated as: In particular, we have: The proof of this lemma is postponed to Appendix A. An element ρM u of Ker(Q) is called a local thermodynamic equilibrium with density ρ and mean velocity u.

Collision invariants
for every f such that f /M u f ∈ V u f and ψ ∈ V u f . The set of CI's is denoted by C. It is a vector space.
We have the: The proof of this proposition is again postponed to Appendix A.

Hydrodynamic limit and fast relaxation
The goal of this section is to investigate the formal limit ε → 0 in (2.13) and to examine some of the properties of the limit system relative to the propulsion force. More precisely, we exhibit a phase transition when the intensity of the velocity fluctuations crosses a certain threshold dependent on the magnitude of the propulsion velocity.

Hydrodynamic limit
The goal of this section is to prove the following formal theorem: We suppose that f I is independent of ε for simplicity. We assume that solutions of (2.13) exist on any time interval [0, T ]. Assume that f ε → f 0 as ε → 0 as smoothly as needed, which means in particular that derivatives of f ε converge to the corresponding derivatives of f 0 . Then, there exist two functions ρ(x, t) > 0 and u(x, t) ∈ R d such that Furthermore, ρ and u satisfy the following system of isothermal compressible Euler equations with relaxation: associated with initial data (ρ I , u I ) such that Proof. First, we note from (2.13) that Q(f ε ) = O(ε). Therefore Q(f 0 ) = 0, which, because of (3.20), implies that f 0 is of the form (4.1). Next, using the CI's given by Proposition 3.3, we multiply (2.13) successively by 1 and v and use the fact that the right-hand side vanishes upon integration. Then, we get the following conservation relations, which are valid for any ε: with ρ ε , j ε , Σ ε , the density, flux and pressure tensor associated to f ε , given by: and right-hand side q ε given by Now, letting ε → 0, we can express j = lim j ε , Σ = lim Σ ε and q = lim q ε as functions of ρ and u: Inserting these expressions into the conservation equations leads to (4.2), (4.3). The statement about the initial conditions is obvious.
This model is nothing but the Euler system of isothermal compressible gas dynamics with a forcing term. The mass conservation equation (4.2) (also known as the continuity equation) has a standard form. The momentum balance equation (4.3) involves an isothermal pressure T ρ on the left-hand side and a self-propulsion force on the right-hand side. The temperature T = σD is proportional to the ratio of the intensities of the velocity fluctuations D and of the social force σ −1 . The self-propulsion force takes the form of a relaxation of the fluid velocity to a comfort fluid velocity a 2 − (d + 2)T . In comparison to the force acting on individual particles, the force acting on the fluid involves a term depending on the temperature. The temperature being a truly macroscopic quantity, this term cannot have any counterpart at the particle level. Here, at the fluid level, the comfort fluid velocity becomes a pure imaginary number when T is larger than a critical temperature T c = a 2 /(d + 2). When the temperature crosses T c , a phase transition occurs. This phase transition is studied in the next section when the intensity τ −1 of the self-propulsion force is taken to infinity, a limit which we refer to as the 'fast relaxation limit'.

Fast relaxation limit in the hydrodynamic model
We recall that, according to our notation, τ is the ratio of the physical relaxation time to the time unit t 0 (see Section 2.2) and is a dimensionless parameter. The goal of this section is to investigate the limit τ → 0 in the hydrodynamic model (4.2), (4.3). We denote by (ρ τ , u τ ) the solution for finite τ and assume that its limit (ρ, u) as τ → 0 exists and is as smooth as needed. We will use the non-conservative form of the model, which we recall here for the sake of convenience: Letting τ → 0 formally in (4.10) leads to |u| 2 + (d + 2)T − a 2 = 0. Therefore, there are two cases according to whether the quantity (d + 2)T − a 2 is positive or negative, i.e. according to the position of T with respect to the critical temperature T c defined by The constant s only depends on the problem data and not on the solution. In this case, equation (4.10) becomes (4.12) To examine the limit τ → 0, we need to find the equilibria of the right-hand side of (4.12), i.e. the solutions u of u(|u| 2 + s 2 ) = 0. Obviously, the only solution is u = 0. Additionally, at least in the spatially homogeneous setting, this is a stable solution. Indeed, the unique solution of satisfies u(t) → 0 as t → ∞. So, in the spatially non-homogeneous case, we formally have u τ → 0 as τ → 0. Therefore, the formal limit of (4.2), (4.3) gives In order to get a more precise description of the limit τ → 0, we need to rescale time and velocity.
With this aim, we let t ′ = τ t and u τ ( . Inserting this into (4.2), (4.3), we find (dropping the tildes): (4.14) The behavior of this system when τ → 0 is that of a diffusion. More precisely, we state: Assume that the solution (ρ τ , u τ ) of system (4.13), (4.14) is smooth and converges smoothly towards a pair (ρ, u) as τ → 0. Then, ρ satisfies the following diffusion equation: Proof. Taking τ → 0 in equation (4.14), we get u = − a 2 T s 2 ∇ x ln ρ = − T Tc T −Tc ∇ x ln ρ. The limit model therefore follows from the continuity equation (4.13) and is given by (4.15).

Case T < T c (small noise)
Let us now consider the dynamics of the system in the case of smaller noise. In this section, we aim to prove the following: Assume that the solution (ρ τ , u τ ) of system (4.9), (4.10) is smooth and converges smoothly towards a pair (ρ, u) as τ → 0. Assume additionally that (ρ τ , u τ ) is not identically equal to (ρ 0 , 0) where ρ 0 is constant in both space and time. Then, u = cω, where ω ∈ S d−1 and the pair (ρ, ω) satisfies the following system: where P = Id − ω ⊗ ω is the orthogonal projection onto the hyperplane orthogonal to ω.
Proof. Equation (4.10) is now written as: (4.18) We now look for the equilibria of the right-hand side of (4.18), namely the solutions of u c 2 − |u| 2 = 0.
There are two sets of equilibria. The first set reduces to the single point u = 0. The second set is the sphere |u| 2 = c 2 . We now show that in the spatially homogeneous setting, the first equilibrium is unstable, while the second class of equilibria is orbitally stable. Indeed, let us consider the solution u(t) of the differential equation Its solution can be analytically given by and is such that |u(t)| → c for all initial data u 0 except u 0 = 0. This shows that u 0 = 0 is an unstable equilibrium. On the other hand, if one perturbs an equilibrium |u 0 | 2 = c 2 by a small amount, the solution will relax to an element of the circle |u| 2 = c 2 (where u may be different from u 0 ). Now, we let τ → 0 in the non spatially homogeneous system (4.16), (4.18). Unless u τ is identically zero, which can only occur for a uniform density ρ 0 , u τ converges (at least formally) towards one of the stable equilibria. Therefore, we have u τ → u = cω with |ω| = 1. In order to find the equation satisfied by ω, we introduce a polar decomposition of the solution u τ . We write u τ = c τ ω τ , with c τ = |u τ | and ω τ = u τ /c τ . Let P τ the orthogonal projection of R d onto the hyperplane orthogonal to ω τ . The projection P τ can be written tensorwise as P τ = Id − ω τ ⊗ ω τ . Inserting the polar decomposition of u τ in (4.18) leads to: (4. 19) We note that, because |ω τ | = 1 and the operator ∂ t + u τ · ∇ x is a derivative, the vector (∂ t + u τ · ∇ x )ω τ is orthogonal to ω τ . Furthermore, the first term of the left-hand side and the right-hand side are parallel to ω τ . Therefore, applying P τ to the second term of the left-hand side leaves it unchanged, while it cancels the first term of the left-hand side and the right-hand side. Consequently, applying P τ to (4.19) leads to: Now, taking the limit τ → 0 formally leads to (4.17).
System (4.16), (4.17) has been referred to in the literature as the Self-Organized Hydrodynamic (SOH) system [24,25]. It has the form of a compressible gas dynamics system with isothermal equation-of-state and geometric constraint |ω| = 1. The projection operator P which multiplies the pressure term ∇ x ln ρ maintains the constraint over the course of time (see Remark 4.3). It results in a non-conservative term (since P depends on ω). This type of system has been derived for the first time in [26] and has been shown to be hyperbolic. Beyond this result, the mathematical study of such systems is in its infancy.
A local existence result is given in [25] and some special solutions are given in [44]. When T < → T c , we have c → 0 and therefore, ρ becomes constant. To find a nontrivial dynamics, we must rescale time to a slower time-scale. By the time rescaling t ′ = ct, system (4.16), (4.17) can be written (dropping the primes): The parameter (T c − T )/T plays the role of the squared Mach-number in the standard compressible Euler system. Therefore, the limit T → T c , is similar to a small Machnumber limit. When T < → T c , we formally get ρ → ρ 0 where ρ 0 is a constant in space. With appropriate boundary conditions we can assume that ρ 0 is also independent of time. Then, the limit system as T < → T c is written as follows: where π is a hydrodystatic pressure defined by This system has been already proposed in [27] for modeling gregariousness. It is a system of incompressible Euler equations subject to the geometric constraint |ω| = 1. The behavior of the (appropriately rescaled) SOH system (4.16), (4.17) when T < → T c is very different from the behavior of the (appropriately rescaled) diffusion equation (4.15) when T > → T c . The former is given by the incompressible Euler equations with geometric constraint |ω| = 1 while the latter is the standard diffusion equation (see Remark 4.2). This drastic change of type is the signature of a phase transition further elaborated on in Section 4.2.3.

Comments
The fast relaxation limit of the compressible Euler system with self-propulsion (4.2), (4.3) exhibits two different regimes for temperatures below or above the critical temperature T c . Above T c , the behavior of the system is that of a diffusion, while below T c , the system obeys a hyperbolic system, the Self-Organized Hydrodynamic (SOH) model (4.16), (4.17). Even when the temperature is close to the critical temperature, it is not possible to match the two types of models in a smooth way, as noticed in the previous section. This abrupt change in the type of the model as the temperature crosses a threshold is a manifestation of a phase transition. Since T c depends on a, the phase transition originates from the  Figure 1: Two strategies. First strategy: first take the relaxation limit τ → 0 (top horizontal arrow) and then pass to the hydrodynamic limit ε → 0 in the resulting constrained kinetic model (right vertical arrow). This is done in [26] and [10]. Second strategy: first pass to the hydrodynamic limit ε → 0 in the original kinetic model (left vertical arrow) and then pass to the fast relaxation limit τ → 0 in the resulting Euler system (bottom horizontal arrow). This is what is done here.
self-propulsion force. Indeed, when a = 0 (i.e. no self-propulsion), T c = 0 and there is no phase transition: the system is in the diffusion regime in all circumstances. Phase transitions in self-propelled particle systems have already been evidenced [53] and an abundant literature has been devoted to them (see the review [54]). Phase transition in hydrodynamic models of self-propelled particles have first been studied in [48] (see also the review [49]). However, the equations proposed in [48] are more complicated than the ones seen here. They are derived solely on heuristic principles and invariance considerations. Their analysis is based on a combination of linear and nonlinear techniques. There is no link to the underlying particle models. The link between hydrodynamic and particle models of self-propeled particle systems has been made in [7,8], but for binary interaction mechanisms instead of the mean-field interaction considered here. Here, the hydrodynamic model is derived from the underlying particle dynamics and is much simpler: it merely consists of the isothermal compressible Euler model complemented with the self-propulsion force. The phase transition manifests itself in the change of type of the PDE which describes the system under large self-propulsion.
The SOH model (4.16), (4.17) has previously been derived in [26] from a system of self-propelled particles which have constant and uniform velocity. In [10], it has been shown that the kinetic model with the velocity norm constraint of [26] is the the fast relaxation limit τ → 0 of (2.13). A natural question is then whether the imposition of the norm constraint at the particle and kinetic levels as in [26] is necessary or useful. Indeed, there are now two ways of deriving the SOH model from (2.13), which are summarized in Fig. 1. The first way is to follow the top horizontal and right vertical arrows successively. This is what is done in [26] and [10]. The second way is to follows the left vertical and bottom horizontal arrows successively. This is what is done here.
There is however one noticeable difference between the two strategies. In [26], the system is written: where the constants c 1 and c 2 are such that 0 < c 2 ≤ c 1 . In [33], it has even been shown that taking into account an anisotropic vision, (c 1 , c 2 ) can be made arbitrary. We would recover (4.16), (4.17) if we could set c 2 = c 1 = c and δ = T /c. However, it is not possible since, in [26], c 2 = c 1 . Therefore, the diagram in Fig. 1 is not commutative. This shows how important it is to enforce the norm constraint at the kinetic level, in spite of the induced difficulties. Indeed, such constrained kinetic models do not exhibit momentum conservation and lack classical collision invariants such as those derived in Section 3.2. This makes the derivation of hydrodynamic models considerably more complex. In [26], new concepts have been developed to bypass these difficulties. The fact that c 1 = c 2 is a direct consequence of these features. The question whether c 1 = c 2 or not has important consequences. Indeed, the SOH system is Galilean invariant if and only if c 1 = c 2 (which is the case here). However, in the generic case c 1 = c 2 , the SOH system is not Galilean invariant. This fact reflects a key feature of collective motion: the anisotropy of information flow. This is evidenced in car traffic, where perturbations of velocities (typically moving with speed c 2 < c 1 ) propagate upstream the flow (whose speed is c 1 ). This phenomenon is a consequence of the fact that information propagates upstream from drivers ahead to drivers behind. When c 1 = c 2 , this property is lost and information spreads in an isotropic way just like usual gas dynamics. In such a case, the model is unable to correctly reproduce the complex emerging patterns (such as congestions and waves in car traffic). We refer to [52] for a discussion of the loss of Galilean invariance in biological swarms.
The approach of [26] does not lead to phase transitions: whatever the values of the temperature T or self-propulsion velocity a, the limit system is of hydrodynamic type. The very large intensity of the self propulsion force which is needed to perform the τ → 0 limit first (top horizontal arrow in diagram 1) prevents any diffusion regime to establish. However, the emergence of a phase transition is a well-established experimental fact (see review in [54]). The inability of [26] to produce phase transitions can be seen as a major drawback. Fortunately, it has been proved in [24] that adding a dependence of the social force intensity σ −1 upon the flux j restores the phase transition. It also brings densitydependent phase transition. More precisely, the parameter controling the phase transition is the ratio T /ρ. A natural question is to explore similar features here. This question will be investigated in future work.
Here, we choose a different direction: we explore whether the inclusion of diffusion terms can cure the deficiency of the SOH models derived here (i.e. the fact that c 1 = c 2 ). This is the goal of the next section.
5 Diffusive corrections to the hydrodynamic model 5

.1 Setting of the problem
In this section, we derive O(ε) diffusive corrections to the compressible Euler model with self propulsion (4.2), (4.3). We show that these diffusive corrections lead to a system of isothermal compressible Navier-Stokes equations including self-propulsion. In a second step, we perform the fast relaxation limit τ → 0 in the resulting Navier-Stokes model and compare it to the SOH model obtained in Proposition 4.3.
We now set up the problem. We define We assume that the observation kernel K(|ξ|) is such that K(|ξ|) dξ = 1 and we denote by k > 0 the second moment of K, i.e.
For instance, if K is the indicator function of the ball of radius 1, In the cases of d = 2 and d = 3, we respectively get k = π/8 and k = 2π/15.
We first give the expansion of F ε f up to the fourth order in ε. where and k R = kR 2 .
We give the proof of this very simple lemma in Appendix B. The O(ε 2 ) correction to the mean velocity u f takes into account the non-local character of the average (2.10). This correction involves gradients of the local density and flux. They quantify how information spreads due to the fact that the agents observe their environment over a certain spatial extent. Because the observation is supposed isotropic, the correction only involves O(ε 2 ) terms and second order derivatives. In the case of non-isotropic observation kernels, O(ε) corrections involving first order gradients would be obtained. The study of this effect is postponed to future work. Now, up to terms of order ε 4 (which are dropped), equation (2.9) is written: (5.27) where again, Q is given by (2.14). Now, we look for a fluid model approximating (5.27) which includes the O(ε) terms. We will adopt Chapman-Enskog's method consisting in writing the O(ε) terms in terms of spatial derivatives only (see e.g. the review in [23]). The model is introduced and discussed in the next section.

Compressible Navier-Stokes equations with self-propulsion
The compressible Navier-Stokes system with self-propulsion is established in the following: Theorem 5.2. Let f ε be the solution of (5.27) associated to a given initial condition f I , which is supposed independent of ε for simplicity. Let (ρ f ε , ρ f ε u f ε ) be the moments of f ε defined by (5.23). Then, we can formally write , provided that (ρ ε , ρ ε u ε ) satisfy the following set of compressible Navier-Stokes-like equations (where we drop the superscripts ε upon ρ ε and u ε for the sake of clarity): , The proof of this result is fairly technical. It is given in Appendix C. The coefficient µ is the fluid viscosity, while π ε (ρ, u) is a velocity-dependent pressure.
We compare this Navier-Stokes model to the hydrodynamic model (4.2), (4.3) and provide a physical interpretation of the O(ε) correction terms. The mass conservation equation (5.28) is unchanged compared to (4.2), as it should be. Indeed, the particle density is still a conservative variable transported by the fluid velocity. By contrast, there is a wealth of new terms in the momentum balance equation (5.29). We can identify the only term which come from classical fluid viscosity: this is εµ∇ x · (ρE(u)). The present viscous term is different from that appearing in the full Navier-Stokes system, where the temperature is determined by the energy balance equation. Here, the viscous term involves the symmetrized velocity gradient tensor E(u) instead of the rate of strain tensor This difference is solely due to the isothermal character of the model and not to self-propulsion.
Then, there are two terms which come from the non-locality of the interaction and which are proportional to k R . The first one contributes to adding more viscosity and is equal to ε k R σ ρ∆ x u. The second one contributes to convecting the velocity in the direction of the gradient of ρ. It is written 2εk R σ (∇ x ρ · ∇ x )u. All other O(ε) terms originate from self-propulsion. The self-propulsion first contributes to a similar relaxation source term as in (4.3): the term − 1 τ ε ρu |u| 2 a 2 − χ ε . However, both the relaxation rate 1 τ ε and the bulk comfort velocity √ χ ε a towards which this source term is relaxing are different from those of (4.3), respectively equal to 1 τ and a 2 − (d + 2)T . We notice that the former are O(ε) corrections of the latter, as they should be. However, there is a second source term, the term which describes a force comprising three terms. The first term is just a friction proportional to the compressibility ∇ x · u. The second term is a force acting in the direction of the gradient of |u| 2 . Finally, the last term is a force acting in the direction of (u · ∇ x )u. The self-propulsion force also induces some changes in the convection terms. The first one is a change in the convection velocity of the momentum, which is not just u but λ ε u. The coefficient λ ε is close to one but not exactly equal to one, the difference being O(ε). This feature makes the model closer to the generic SOH model (4.21), (4.22), which has different convection velocities c 1 and c 2 for the density and velocity. The last contribution of the self-propulsion is a modification of the pressure. The isothermal pressure ρT is complemented by an O(ε) pressure correction επ(ρ, |u|) which depends on the norm of the velocity. This correction is positive or negative according to whether the speed is larger or smaller than the bulk comfort velocity a 2 − (d + 2)T of the non-viscous model.
We now discuss this model in view of the model proposed in [48,49], which is considered as the 'paradigmatic' model of hydrodynamic type for flocking. In this model, the mass conservation equation is the same as (5.28) but the momentum balance equation takes the following form (in the absence of extrernal force): where α, β, D L , D 1 and D 2 are positive coefficients and P = P (ρ) is the pressure, which depends nonlinearly of ρ. This model has been derived on the basis of invariance considerations. In order to compare our model with (5.35), we use the mass conservation equation (5.28) to write the momentum balance equation (5.29) as follows: There is one term in (5.35) which is missing from (5.36): the anisotropic velocity diffusion D 2 (u·∇ x ) 2 u. By contrast, there are many terms appearing in (5.36) which are not present in (5.35): the third, fourth and fifth terms of the left-hand side and the third and fifth terms of the right-hand side. Additionally, among the terms which are common to both formulas (namely the first and second terms of the left-hand side and the first, second, fourth and sixth terms of the right-hand side of (5.36)), some of them assume different forms. Indeed, the second term (the convection term (u · ∇ x )u) is multiplied by the constant 1 − ε 3λ 2 less than one in (5.36), while this coefficient is exactly one in (5.35). This difference is significant, in view of the previous discussion about mass and momentum convection velocities in swarming systems. Another difference is in the pressure term (second term at the right-hand side of (5.36)). In our model, the pressure depends on both the density and the norm of the velocity, while it depens on the density only in (5.35). By contrast, the dependence upon the density is linear in our case, while it is nonlinear in (5.35). This discussion illustrates that phenomenological models can differ significantly from first principle models when complex phenomena such as swarming behavior are concerned. The question whether these differences lead to perceivable changes in the qualitative behavior of the solution has not yet been investigated.
The mathematical properties of this system (such as e.g. the stability of this system for small perturbations of a homogeneous state) will be studied in future work. In this work, we investigate the fast relaxation limit τ → 0. This is the goal of the next section.

The fast relaxation limit in the compressible Navier-Stokes equations with self-propulsion
In this section, we examine the limit τ → 0 in the Navier-Stokes system (5.28), (5.29).
Since this system was derived under the assumption that ε is small, ε needs to tend to 0 at least as fast as τ tends to 0. Here, we decide to make ε and τ proportional. This is the borderline case, because, to be consistent, we should have linked ε and τ in this way already at the kinetic level. However, as pointed out earlier, the analysis of this scaling limit is more complex and is still under scrutiny [3]. Our conjecture is that the kind of model we get is the same in both limits. Therefore, investigating it at the level of the Navier-Stokes system is a preparation before performing the limit directly from the kinetic level. Since λ is proportional to 1/τ , we decide to relate ε and λ in such a way that ελ is a constant α, i.e. α = ελ. (5.37) We need to keep in mind that, strictly speaking, α must be ≪ 1 otherwise the derivation of the Navier-Stokes model in the previous section loses its valididty. We decide to express ε as a function of τ . We can write We also define We can introduce an α-dependent critical temperature T c (α) by: where T c (0) is the critical temperature (4.11). Then, we can write the square of the comfort velocity in the Navier-Stokes model c 1 (α) as Since α must be small, we can limit its range in such a way that ξ α > 0 and that the relaxation time τ α remains positive. Therefore, we have α ∈ [0, 2 d+8 ]. In this interval, T c (α) is an increasing function of α with values in [T c (0), 3 2 T c (0)]. We also assume that the temperature is below the critical temperature T < T c (α) in such a way that c 2 1 (α) > 0. We finally denote by (ρ τ , u τ ) the solution of the Navier-Stokes system (5.28), (5.29). We re-write the system in the new notation: Now, we formally let τ → 0 in this system, keeping all other parameters fixed, and in particular α. We get the following: ] and let T < T c (α). Assume that the solution (ρ τ , u τ ) of system (5.44), (5.45) is smooth and converges smoothly towards a pair (ρ, u) as τ → 0. Assume additionally that (ρ τ , u τ ) is not identically equal to (ρ 0 , 0) where ρ 0 is constant in both space and time. Then, u = c 1 (α)ω, where ω ∈ S d−1 and the pair (ρ, ω) satisfies the following system: where P = Id − ω ⊗ ω is the orthogonal projection onto the hyperplane orthogonal to ω and Remark 5.1. For small α, up to terms of order α 2 , we have: Therefore, the convection speeds are larger when first order corrections are included. Additionally, it is easy to prove that c 1 (α) is an increasing function of α on its interval of definition.
Proof. The proof is similar to that of Proposition 5.3 and is only sketched. First, observe that the second line of the right-hand side of (5.45), which is proportional to τ , simply vanishes in the limit. Since all the diffusion terms in the Navier-Stokes system are contained on this line, there is no diffusion in the limit system. Second, since |u τ | → c 1 (α), which is a constant, the second term of the third line of (5.45) also vanishes in the limit. For the same reason, the pressure π α (ρ τ , |u τ |) → π α (ρ, c 1 (α)) = T α ρ. Therefore, we recover an isothermal pressure equation-of-state in the limit. Then, since (5.49) is obtained by projecting (5.45) onto the hyperplane normal to u τ , the first term of the third line of (5.45), which is parallel to u τ , also vanishes in the limit. Finally, the last term of the third line of (5.45) combines with the second term at the left-hand side of (5.45) and yields the second term of the left-hand side of (5.47). Indeed, we get where P τ is defined like in the proof of Proposition 5.3. Therefore, c 2 (α) = (λ α − α 2 )c 1 (α) is given by (5.48). The remaining details are left to the reader. Now, system (5.46), (5.47) is in the form of the generic SOH model (4.21), (4.22) with c 2 < c 1 , like in [26]. We conclude that the inclusion of diffusive terms in the hydrodynamic model before passing to the fast relaxation limit was a successful approach. It has cured the deficiency of the SOH model derived in Section 4.2.2. In some sense, the purely hydrodynamic model obtained at Theorem 4.1 is too simple to describe the full complexity of the system but the inclusion of diffusive terms is enough to restore the adequate level of complexity.
It is intriguing that the only diffusive corrections that are kept in the fast relaxation limit are those coming from the self-propulsion force. It would be interesting to compute the O(τ ) corrections to this model. Then, the other diffusive terms (those arising from viscosity and non-locality of the interaction) would appear. The resulting model would then have more relevance as an approximation of the original kinetic model. It should also be compared to the diffusive SOH model obtained in [28], which is quite complex. It would be instructive to see if the present approach could help get a cleaner model.

Conclusion and perspectives
In this work, we have provided evidence of a phase transition from disordered to ordered motion in a hydrodynamic model of socially interacting agents with self-propulsion. The model we have investigated has been derived from a particle system combining a noisy Cucker-Smale consensus force and self-propulsion. We have shown that the phase transition appears in the limit of a large self-propulsion force and manifests itself as a change of type of the limit model (from hyperbolic to diffusive) at the crossing of a critical noise intensity. We have also shown that, in the hyperbolic regime, the resulting SOH (self-Organized Hydrodynamics) model suffers from unnecessary restrictions on the range of its coefficients. To remove these restrictions, we have computed diffusive correction to the model. With these diffusive corrections, the restrictions on the SOH model obtained in the limit of a large propulsion force disappear.
As pointed out in the core of the work, many points deserve further elucidation. A first one, currently under scrutiny, consists in performing the combined hydrodynamic and large self-propulsion force simultaneously at the level of the kinetic model. We anticipate that similar phase transitions will emerge and will be described by the same limit models, with possibly different coefficients. Other points would be worth being developed. The computation of diffusive corrections to the SOH model for instance would be of great practical use. Further investigations of quantities attached to the social force are also very promising, such as the role of a possible anisotropy of the observation kernels, or of a different scaling of its range. Finally, the understanding of the transition between the two phases and how the two different models can be matched at this transition is also crucial for applications.
for all f such that f /M u ∈ V u . We notice that for given u ∈ R d , R(u; f ) is a linear operator with respect to f and that Then, ψ is a CI of Q if and only if or equivalently, if and only if: ∀f such that f /M u ∈ V u and u f = u.
As a first step, we fix u ∈ R d and find all ψ ∈ V u which satisfy: Then, we will make u arbitrary. We note that any constant ψ is a solution of (A. 52).
Saying that u f = u is equivalent to saying that f (v − u)dv = 0. So ψ defined by (A.52) is such that the following implication holds: Since both sides of the implication of (A.53) are linear forms of f , by a standard theorem [13], there exists β u ∈ R d such that Using (A.51) and Green's formula (see (3.18)), and introducing the change of function g = f /M u , it follows that ψ is a solution of the following problem: with β ′ u = −β u /D. We will drop the primes in the following. Let φ = β u · (v − u). Then, problem (A.54) for ψ can be equivalently written according to the variational formulation: Since any constant ψ is a solution of (A.52), we can subtract ψM u dv from ψ and assume, without loss of generality, that ψ is such that ψ M u dv = 0. Next, we have: by the definition of M u . Therefore, if g is a constant, we have both (ψ, g) Vu = 0 (by the definition (3.16) of (·, ·) Vu ) and (φ, g) Hu = 0 (by (A.56)). Therefore, it is possible to restrict (A.55) to functions g such that (g, 1) Hu = 0. Consequently, we define the following space:V and variational formulation (A. 55) can be made precise as follows: Find ψ ∈V u such that (ψ, g) Vu = (φ, g) Hu , ∀g ∈V u . (A.57) The Poincaré inequality for Gaussian measures shows that for all φ ∈ V u , Then following from (A.58), the bilinear form (·, ·) Vu is coercive onV u . By the Lax-Milgram theorem, we deduce that there exists a unique ψ ∈V u such that (A.57) holds. But on the other hand, obvious calculation shows that ψ(v) = T β u · (v − u) is a particular solution. Since it belongs toV u , it is the only solution of the variational formulation (A.57).
Adding any constant, we have just shown that any solution of (A.52) is of the form ψ(v) = α + β · v, where α ∈ R and β ∈ R d are arbitrary. Since this form is independent of u, we deduce that such ψ's are solutions of (A.52) for all u ∈ R d and are therefore the only CI's. This ends the proof of Proposition 3.3.
B Appendix: Expansion of the force term Proof of Lemma 5.1. By the change of variables y = x − εR ξ and Taylor's formula, we have We have used the definition of k and the evenness of K with respect to ξ in order to cancel the odd order terms of the expansion. We have denoted by D 2 x f the Hessian matrix of f with respect to x (i.e. the matrix of the second order derivatives) and the symbol ':' refers to the contracted product of tensors. In a similar way, we have: Therefore, by expanding the ratio of (B.59) and (B.60) up to the fourth order and in view of (2.10), we find (5.24) with (5.26). Formula (5.25) immediately follows from (2.10).
C Appendix: proof of Theorem 5.2 Multiplying (5.27) by the collision invariants 1 and v, we are led to the conservation equations (in the same way as in Section 4): with ρ ε , j ε = ρ ε u ε , Σ ε , q ε given by (4.6) and (4.7) and . We also note that we can write The Chapman-Enskog expansion consists in closing the expressions of S ε , q ε and r ε by a first order expansion of f ε . To this end, we write the so-called macro-micro decomposition: By the definition of ρ ε and u ε , we have The first term at the right-hand side of (C.65) is the macroscopic part, as it is proportional to a local thermodynamical equilibrium and carries all information about the moments of the solution. The second term is the microscopic part. It carries no information about the macroscopic moments but instead carries information about the discrepancy between f ε and the local thermodynamical equilibrium. From Section 4, we know that, provided that ρ and u satisfy the Euler equations, the microscopic part is small of order ε. This is why this microscopic part is multiplied by ε in (C.65). We stress the fact that there is no approximation involved (at this step) in (C.65): it is a mere definition of f ε 1 . Inserting (C.65) into the formulas providing the expressions of the quantities involved in the moment equations (C.61), (C.62) (specifically equations (C.64), (4.7) and (C.63)), we can write the moment equations as follows: where we have used Green's formula for B ε 1 and (5.26), (C.63) for B ε 2 . Now, in order to compute B ε 1 and B ε 3 , we need to evaluate f ε 1 . But since we only look for O(ε) correction terms, and f ε 1 is multiplied by ε, we may compute f ε 1 up to terms of order O(ε). Inserting (C.65) into (3.17), the collision operator Q can be written: Inserting it in (5.27), we get: But we can neglect all terms of order ε or more in (C.73). Therefore, we are led to the following equation for f ε 1 : The inversion of (C.74) will give us f ε 1 . Once f ε 1 is obtained, we insert it into (C.69)-(C.70) and this leads us to the expressions of the Navier-Stokes terms. We first compute the right-hand side R ǫ : Lemma C.1. We have (dropping the superscript ε for the sake of clarity): For the third term (C.75), we compute Collecting the second term at the right-hand side of (C.79) together with (C.80) leads to an expression S which we can split in the following way: Collecting this expression with the first term at the right-hand side of (C.79) leads to (C.76).
In order to solve equation (C.74) for f ε 1 , we need to solve equations of the type where u is an arbitrary vector of R d and g is a given function. We refer the reader to Section 3.2 for the definitions of the spaces H u , V u andV u . We state the following lemma, whose proof is identical to that of proposition 3.3 and is left to the reader.
Proof. The proof that h(v − u)M u through g(v − u)M u satisfy (C.83) and (C.85) easily follows from classical formulas for moments of the Gaussian, which we leave to the reader. Then, we apply Lemma C.2, which gives the existence of L −1 u (h(v−u)M u ) through L −1 u (g(v−u)M u ) and the fact that they satisfy (C.84) and (C.86). Formulas (C.87), (C.88) follow from explicitly computing the action of L −1 u and using the uniqueness statement of Lemma C.2. Finally, equation (C.74) for f ε 1 , which can be written (up to order O(ε) terms) −L u f ε 1 = R, can be solved by f ε 1 = −L −1 u R since according to the first equation of (C.66), f ε 1 satisfies (C.84). By the linearity of L u and the decomposition (C.76) of R, we can write: Thanks to (C.87), (C.88), equation (C.89) follows.