# Diffusion in a continuum model of self-propelled particles with alignment interaction

###### Abstract

In this paper, we provide the corrections to the hydrodynamic model derived by Degond & Motsch from a kinetic version of the model by Vicsek & coauthors describing flocking biological agents. The parameter stands for the ratio of the microscopic to the macroscopic scales. The corrected model involves diffusion terms in both the mass and velocity equations as well as terms which are quadratic functions of the first order derivatives of the density and velocity. The derivation method is based on the standard Chapman-Enskog theory, but is significantly more complex than usual due to both the non-isotropy of the fluid and the lack of momentum conservation.

1-Universit de Toulouse; UPS, INSA, UT1, UTM ;

Institut de Math matiques de Toulouse ;

F-31062 Toulouse, France.

2-CNRS; Institut de Math matiques de Toulouse UMR 5219 ;

F-31062 Toulouse, France.

email:

3- Department of Mathematics,

City University of Hong Kong,

83 Tat Chee Avenue, Kowloon, Hong Kong

email:

Acknowledgements: This work was supported by the Marie Curie Actions of the European Commission in the frame of the DEASE project (MEST-CT-2005-021122) and by the french ’Agence Nationale pour la Recherche (ANR)’ in the frame of the contract ’Panurge’ (ANR-07-BLAN-0208-03).

Key words: Flocking, Vicsek model, alignment interaction, asymptotic analysis, hydrodynamic limit, Chapman-Enskog expansion.

AMS Subject classification: 35Q80, 35L60, 35K99, 82C22, 82C31, 82C44, 82C70, 92D50

## 1 Introduction

This paper is a development of a previous work [26, 27] about continuum models of self-propelled particles subject to alignment interaction. This class of models describes swarming behaviour among biological species and attempts at providing a simplified theoretical framework to experimental observations (see recent observations in Refs. [3, 5, 8, 32, 40]).

The starting point is a particle model (or Individual-Based Model (IBM)) discussed in e.g. Refs. [2, 16, 39, 41], where the interaction between biological agents such as fish or birds is described by three interaction ranges: a close range where repulsion occurs, a long range where attraction prevails and a medium range where agents tend to align with each other. In the Vicsek model,[53] the alignment interaction is singled out and analyzed. More precisely, each particle moves with a constant and uniform speed and aligns with the average direction of all neighbours within an interaction distance , up to some angular fluctuation. Vicsek and co-authors[53] show that phase transitions from disorder to order appear as the noise intensity decreases or the density increases. This model has triggered a wealth of publications[1, 4, 13, 22, 33, 34] and given rise to various variants[14, 50].

The task of deriving kinetic (Boltzmann-like) or continuum (fluid-like) models from this model has been undertaken by various approaches (see Refs. [6, 46] for kinetic models and Refs. [6, 21, 23, 46, 47] for fluid models). However, these models are based on physical arguments and the mathematical approach of Ref. [26, 27] leads to a different type of model: indeed, this model is not of diffusive nature like in Refs. [6, 21, 23, 46] and is local, by contrast to Ref. [47]. The differences arise because the model of Refs. [26, 27] is derived in the large time and space scale limit in which, at leading order, the interactions are local and the diffusivities, negligible.

The goal of this paper is to investigate how the model of Ref. [26, 27] must be adapted to account for small but finite nonlocality and diffusivity. We will see that the resulting model involves complex effects due to the non-isotropy of the fluid (by contrast to standard fluids). This point will be further developed in section 2. Similar to the rarefied gas dynamics case (see e.g. Refs. [9, 12, 24, 49]), the derivation is based on the Chapman-Enskog expansion method. However, in the present case, the computations are significantly more complex because of the non-isotropy of the fluid and of the lack of momentum conservation.

The model derived in Ref. [26, 27] has been extended in Ref. [31] to account for anisotropic vision and density-dependent interaction frequency. Its numerical resolution is performed in Ref. [44].

Other Individual-Based Models involving attraction-repulsion interactions can be found in Refs. [28, 30, 43, 45] and continuum models, in Refs. [7, 10, 15, 20, 25, 42, 51, 52]. The existence of flocking or non-flocking behaviour for the Cucker-Smale model[17, 18, 19] (which is similar to the Vicsek model but without noise nor speed constraint) has received a great deal of attention[11, 29, 35, 36, 37, 38, 48].

## 2 Position of the problem and main result

In Refs. [26, 27] the following continuum model of self-driven particles with alignment interaction has been derived:

(2.1) | |||

(2.2) |

where is the number density of the particles and is the direction of their average velocity, which satisfies . , and are constants which are computed from the underlying microscopic dynamics, where and , , . The matrix denotes the orthogonal projection matrix onto the plane orthogonal to . The notations Id and respectively refer to the identity matrix and to the tensor product.

The derivation of this model, according to Ref. [27], proceeds through several (formal) asymptotic limits. The starting point is a time-discrete particle model proposed by Vicsek and co-authors[53]. In the Vicsek model, the particles move with a constant speed and, at discrete times, align their velocities to the mean velocity of their neighbours, up to some small noise. In Ref. [27], a continuous in time version of this particle model is first proposed in which the alignment interaction is modeled through a relaxation term and the noise, by a Brownian motion on the particle velocities. The formal mean-field limit of this time-continuous particle model (as the number of particles tends to infinity) leads to the following nonlinear Fokker-Planck model:

(2.3) | |||

(2.4) | |||

(2.5) |

Here is the particle distribution function depending on the space variable , the velocity direction and the time . is a scaled diffusion constant associated with the Brownian noise, and is the mean-field interaction force between the particles which depends on an interaction frequency . This force tends to align the particles to the direction which is the direction of the particle flux in a neighbourhood of weighted by the kernel . Typically, if is the indicator function of the ball centered at and of radius then is the particle flux integrated over a ball centered at and of radius . The matrix is the projection matrix onto the normal plane to and we assume that the collision frequency may depend on , i.e. on the cosine of the angle between and .

Notation convention: The (and below, ) symbols indicate the nabla and Laplacian operators with respect to while and denote the nabla and Laplace-Beltrami operators with respect to . Expressions of and in spherical coordinates will be recalled later.

System (2.3)-(2.5) is written in a scaled form: the time and space scales have been chosen such that the particle speed is exactly and that both and are of order unity. With these so-called microscopic scales, the typical time and distance between two particle interactions are both . We refer to Ref. [27] for a discussion of this point.

By contrast, model (2.1), (2.2) is designed to capture the large scale dynamics only, while averaging out the microscopic scales. Therefore, the passage from model (2.3)-(2.5) to model (2.1), (2.2) requires a change of scales. Let be a measure of the ratio of the microscopic length scale to the size of the observation domain. Here, the relevant scaling is a hydrodynamic scaling, which means that is also equal to the ratio of the microscopic time scale to the macroscopic observation time. To express model (2.3)-(2.5) in terms of the macroscopic time and space scales, we perform the change of variables , . In doing so, we must assign a scale to the interaction kernel (i.e. to the interaction range ). A key assumption in the present work, following Ref. [27], is that this interaction range is microscopic, and therefore of order when using the macroscopic scales. This change of variables leads to (dropping the tildes for clarity):

(2.6) | |||

(2.7) | |||

(2.8) |

It is an easy matter to see that has the following expansion:

(2.9) | |||

(2.10) | |||

(2.11) |

is the direction of the local flux and the constant depends on the interaction kernel through:

Accordingly, can be expanded:

(2.12) | |||

(2.13) | |||

(2.14) |

where is the derivative of with respect to . Because of the dependence of upon the distance only, all odd powers of vanish in the expansion. This would not be the case if we considered more general kernels such as those[31] depending on the angle between and . The consideration of more general kernels is left to future work.

Consequently, we will consider the following expanded Fokker-Planck model

(2.15) |

where the terms and are defined by (2.13) and (2.14) We note that, at leading order, the interaction force only depends on the local flux and that the corrections due to the nonlocality of the interaction force appear in the terms only. This is due to the assumption that the radius of the interaction region is very small (of order ) in the macroscopic variables.

In Ref. [27], it has been proved that model (2.1), (2.2) is the formal hydrodynamic limit of the mean-field model (2.15). Additionally, Ref. [27] provides the connection between the coefficients , and of the macroscopic model to the coefficients and of the microscopic one. The goal of this paper is to investigate what diffusive corrections are obtained when keeping the corrections in the Chapman-Enskog expansion of . These terms describe the response of the system to the appearence of gradients of the state variables and . Here, because the fluid is anisotropic, and has only invariance through rotations about , these gradients must be split into their components parallel and perpendicular to .

To this aim, we denote by

(2.16) |

the orthogonal projection matrices onto the plane normal to and onto the line spanned by respectively. For a given vector , we recall that

Using these projections, any vector field and tensor field can be decomposed into parallel and transverse components according to:

(2.17) |

defined by

Now, we decompose gradient fields according to their parallel and normal components to . For a scalar function , we define the normal and parallel gradients as

Similarly, we may decompose the gradient of a vector field into

using the tensor decomposition (2.17). Applying this decomposition to itself, we find:

(2.18) | |||

(2.19) |

The last line is a consequence of , which is found by taking the derivative of the relation .

In compressible Navier-Stokes equations[24], the diffusion terms can be expressed as functions of only two quantities constructed with the gradient of the velocity field : the traceless rate of strain tensor , and the divergence field (the exponent denotes the matrix transpose). Here, the anisotropy of the problem gives rise to different diffusivities in the directions parallel or normal to and we need to split the matrix into a larger number of separate entities. To this aim, we note that

(2.20) |

where ’Tr’ denotes the trace of a tensor. The traceless tensor is decomposed in its symmetric and anti-symmetric parts and :

(2.21) | |||||

(2.22) |

These relations will be used in the form:

(2.23) | |||||

(2.24) |

The non-zero block of is a 2 by 2 symmetric traceless tensor and the non-zero block of is a 2 by 2 anti-symmetric tensor. We note that, for a given vector :

(2.25) |

Similarly, through (2.18), depends only on and we have:

(2.26) |

where denotes the curl of . Physically, describes the rate of tilt of as one moves along the flow lines (see figure 1). The other quantities describe elementary flow patterns in the plane normal to : refers to convergent or divergent flows in the direction normal to while refers to swirling patterns around (see figure 2) and to shear patterns with one converging and one diverging orthogonal directions. Restricted to the plane normal to , is a symmetric traceless matrix. Therefore, it can be expressed as a linear combination of the two elementary matrices:

each corresponding to an elementary flow pattern (see figure 3).

In the compressible Navier-Stokes equations, there are no diffusion terms in the mass equation and diffusion terms in the energy equation are expressed in terms of the temperature gradient as a whole. Here, the temperature is constant (it is fixed by the noise level) but, because of the lack of momentum conservation, the diffusion terms in the mass conservation equation are not zero. Additionally, both the mass and velocity diffusions depend on as well as on . Therefore, similar to , we decompose into its parallel and normal components (respectively and ).

Finally, the diffusion terms are composed of two parts: the first one is a quadratic form of the gradients in the set

###### Theorem 2.1

(formal) The following model

(2.27) | |||

(2.28) |

where and given below, provides a second order approximation of the moments of the solution and of the initial model (2.15). The right-hand-sides are given by:

(2.29) | |||

(2.30) | |||

(2.31) | |||

(2.32) |

where , , for , for are coefficients, possibly depending on , which are given below. Additionally, we have .

The structure of is as announced: it is decomposed into a term which is a quadratic function of the gradients and a term which consists of derivatives of these gradients. Both terms involve coefficients which may depend on . The quadratic part combines products of the parallel gradient of , , with the perpendicular gradient of , , and similarly for (the parallel gradient of being , the perpendicular ones being defined as any of gradients in the list ) and products of the perpendicular gradient of with the perpendicular gradients of or parallel gradients of with parallel gradients of . The diffusive part involves only parallel gradients of the parallel gradient of , or perpendicular gradients of perpendicular gradients of , and finally a parallel gradient of a perpendicular gradient of . In spite of its complex expression, has a lot of structure, since a general function of this form would have 21 different terms in the quadratic part and 10 in the diffusion part, instead of respectively 8 and 5.

The main objective of this paper is the proof of this theorem. In future work, the properties of this system will be analyzed. In particular, the question of the well-posedness of the system under the sole condition will be investigated, at least on a simpler model. Indeed, the other diffusion coefficients have no definite signs, but the constraint should prevent the formation of singularities, by contrast to the usual backwards heat equation.

The organization of the proof is as follows. In section 3, we recall the main properties of the collision operator that are proved in Ref. [27] and provide additional properties of the linearized operator about an equilibrium. This prepares the terrain for the Chapman-Enskog expansion which is performed in section 4. We start section 4 by an exposition of the main steps to be accomplished. Then, we successively examine the solvability condition for the existence of the first order correction to the equilibrium, the finding of an analytical expression of it as a function of elementary solutions of the linearized collision operator, and finally, the computations of the moments of this correction which precisely give rise to the expressions of and . We begin with the properties of in the next section.

## 3 Properties of the collision operator and of its linearization

### 3.1 Preliminaries

We recall the expressions of the gradient and divergence operator on the sphere. Let be a cartesian coordinate system associated with an orthonormal basis and let be a spherical coordinate system associated with this basis, i.e. , , . Let also be the local basis associated with the spherical coordinate system ; the vectors and have the following coordinates in the cartesian basis:

denotes the Laplace-Belltrami operator on the sphere:

We write for . We introduce the ’collision’ operator, which corresponds to the leading order term of (2.15):

(3.1) | |||

(3.2) | |||

(3.3) |

We note that is a non linear operator. From now on, we assume that is as smooth and integrable as necessary. We note that acts as an operator on functions of only and that the possible dependence of these functions on can be ignored. The properties of have been demonstrated in Ref. [27] are developed in the next section.

### 3.2 Properties of

#### 3.2.1 Null-space of

For , let . We denote by an antiderivative of , i.e. . We define

(3.4) |

The constant is set by the normalization condition (second equality of (3.4)) ; it depends only on and on the function but not on . We note that, when is constant, and that is the so-called Von-Mises distribution. The Von-Mises distribution extends the notion of Gaussian for functions defined on the sphere and is also known as the circular Gaussian. In the present case, the Von-Mises distribution is centered at (or peaked at) .

The following lemma states what are the elements of the null-space of , i.e. what are the equilibria of the problem (see Ref. [27] for the proof):

###### Lemma 3.1

(i) The operator can be written as

(3.5) |

and we have

(3.6) |

(ii) The equilibria, i.e. the functions such that form a three-dimensional manifold given by

(3.7) |

and is the total mass while is the direction of the flux of , i.e.

(3.8) | |||

(3.9) |

Furthermore, if and only if .

#### 3.2.2 Generalized collision invariants

The second set of lemmas state what are the generalized collision invariants of . Indeed, we recall that the collision invariants are classically defined as the functions such that

(3.12) |

However, it is readily seen that the linear vector space of collision invariants is of dimension one, while the hydrodynamic limit requires that its dimension be equal to the dimension of , which is in the present case. To find the missing collision invariants, we slightly weaken the definition. We fix arbitrarily, and we define a Generalized Collision Invariant (or GCI) associated to as a function which satisfies (3.12) only for functions with direction . This constraint is linear and can be resolved by the introduction of a Lagrange multiplier. This leads to the following definition:

###### Definition 3.2

Let be given. is a Generalized Collision Invariant (or GCI) associated to if and only if

(3.13) |

It is shown in Ref. [27] that, using (3.5) and Green’s formula, (3.13) leads to the following problem defining the GCI’s associated to a direction : , such that and

(3.14) |

This problem is obviously linear, so that the set of GCI’s associated to is a linear vector space. In a cartesian basis and the associated spherical coordinates , we have with , . Therefore, we can successively solve for and , the solutions of (3.14) with right-hand sides respectively equal to and . The following lemma provides the framework for solving (3.14). It is based on Lax-Milgram theorem and is proved in Ref. [27]:

###### Lemma 3.3

Let such that . The problem

(3.15) |

has a unique weak solution in the space , the quotient of the space by the space spanned by the constant functions, endowed with the quotient norm.

So, to each of the right-hand sides or which have zero average on the sphere, there exist solutions and respectively (unique up to constants) of problem (3.15). We single out unique solutions by requesting that and have zero average on the sphere: , . Then, we have, as a consequence of Lemma 3.3:

###### Proposition 3.4

The set of generalized collisional invariants associated with the vector which belong to is a three dimensional vector space .

More explicit forms for and can be found. By expanding in Fourier series with respect to , we easily see that

(3.16) |

where is a solution of the elliptic problem on :

(3.17) |

To solve this problem, we apply the following lemma:

###### Lemma 3.5

Let , . Let belong to such that there exists and . Then, for any , there exists a unique solution of the problem

(3.18) |

Additionally, the maximum principle holds: if is non-positive (respectively non-negative), then, so is .

The next lemma, in the spirit of the previous one, will prove useful in the sequel:

###### Lemma 3.6

Let , , , . Then, for any , there exists a unique solution of the problem

(3.19) |

Both lemma are direct consequences of Lax-Milgram’s theorem. The first one is proved in Ref. [27]. For the second Lemma, thanks to a Poincaré inequality we note that the semi-norm of is equivalent to the norm of on the quotient . For both problems, no boundary conditions at need to be prescribed. This is due to the degeneracy of the elliptic operator at these points.

Applying Lemma 3.5 shows that the function , solution of (3.17), is uniquely defined in the space . For convenience, we introduce or equivalently . We then define

(3.20) |

is the vector generalized collisional invariant associated to the direction and is uniquely defined by the problem

(3.21) |

We note that, by the maximum principle, .

### 3.3 Linearization about the equilibrium state

We introduce a macro-micro decomposition of :

(3.22) |

and being the density and mean velocity direction of , i.e.

(3.23) |

and given by (3.3). These definitions of and are equivalent to saying that belongs to the space:

(3.24) |

For any given , we define the following linear operator operating on :

(3.25) |

Inserting (3.22) into (3.5) yields that

We now precise the functional setting. Let

###### Lemma 3.7

Consider as an operator from into its dual , defined by the bilinear form on :

(3.26) |

Then,

(i) The null-space of is the linear space spanned by .

(ii) Let . The equation has a solution in if and only if satisfies the solvability condition:

(3.27) |

And is unique if it additionally satisfies (3.27). This unique solution is written and is called the pseudo-inverse of .

Proof: Statement (i) is obvious. To prove (ii), we note that the bilinear form (3.26) defines a bilinear form on the quotient space . The solvability condition (3.27) is the necessary and sufficient condition for an element of to belong to the dual space . Then, the proof of (ii) follows from the application of Lax-Milgram theorem. Indeed, by Poincare’s lemma, the bilinear form (3.26) defines a norm on which is equivalent to the norm of . Therefore, there exists a unique solution of the problem

(3.29) |

which uniquely defines a solution in provided the cancellation condition (3.27) is imposed on . Finally, inserting in (3.29) and using (3.21), we deduce that

Property (iii) is an immediate consequence of this identity.

### 3.4 Coefficients of the hydrodynamic model

We finish these preliminaries by recalling the expressions of the coefficients , and as they were derived in Ref. [27]. For any functions , with , we denote by the average of over the probability distribution defined by :

(3.30) |

Then, the constants , and of model (2.1), (2.2) are given by:

(3.31) |

With simple computations, one can check that these definitions are equivalently expressed by the following relations:

(3.32) | |||

(3.33) | |||

(3.34) |

## 4 The Chapman-Enskog expansion

### 4.1 Setting up the expansion

We introduce the macro-micro decomposition (3.22) in a scaled form:

(4.1) |

where and are the density and velocity direction of the solution of the kinetic model (2.15) and where . This leads to:

(4.2) |

This justifies the scaling (4.1) because appears as an quantity, showing that the correction to equilibrium is .