Continuum limit of self-driven particles with orientation interaction
aa r X i v : . [ m a t h - ph ] O c t Continuum limit of self-driven particles withorientation interaction
P. Degond (1) , S. Motsch (1) (1) Institute of Mathematics of Toulouse UMR 5219 (CNRS-UPS-INSA-UT1-UT2),Universit´e Paul Sabatier, 118, route de Narbonne, 31062 Toulouse cedex, Franceemail: [email protected], [email protected]
Abstract
We consider the discrete Couzin-Vicsek algorithm (CVA) [1, 9, 19, 36], whichdescribes the interactions of individuals among animal societies such as fish schools.In this article, we propose a kinetic (mean-field) version of the CVA model andprovide its formal macroscopic limit. The final macroscopic model involves a con-servation equation for the density of the individuals and a non conservative equationfor the director of the mean velocity and is proved to be hyperbolic. The derivationis based on the introduction of a non-conventional concept of a collisional invariantof a collision operator.
Acknowledgements:
The first author wishes to thank E. Carlen and M. Carvalho fortheir interest in this work and their helpful suggestions.Preprint of an article submitted for consideration in Mathematical Models and Methodsin Applied Sciences (M3AS) c (cid:13)
Key words:
Individual based model, Fish behavior, Couzin-Vicsek algorithm, asymp-totic analysis, orientation interaction, hydrodynamic limit, collision invariants
AMS Subject classification:
Introduction
The discrete Couzin-Vicsek algorithm (CVA) [1, 9, 19, 36] has been proposed as a modelfor the interactions of individuals among animal societies such as fish schools. The in-dividuals move with a velocity of constant magnitude. The CVA model describes in adiscrete way the time evolution of the positions of the individuals and of their velocityangles measured from a reference direction. At each time step, the angle is updated to anew value given by the director of the average velocity of the neighbouring particles, withaddition of noise. The positions are updated by adding the distance travelled during thetime step by the fish in the direction speficied by its velocity angle.For the modeling of large fish schools which can reach up to several milion individuals,it may be more efficient to look for continuum like models, which describe the fish societyby macroscopic variables (e.g. mean density, mean velocity and so on). Several suchphenomenological models exist (see e.g. [25, 34, 35]). Several attemps to derive continuummodels from the CVA model are also reported in the literature [23, 29, 30], but thederivation and the mathematical ’qualities’ of the resulting models have not been fullyanalyzed yet. One can also refer to [16, 26] for related models. An alternate model, thePersistent Turning Walker model, has been proposed in [18] on the basis of experimentalmeasurements. Its large-scale dynamics is studied in [12]. Additional references on swarmaggregation and fish schooling can be found in [7]. Among other types of animal societies,ants have been the subject of numerous studies and the reader can refer (among otherreferences) to [21, 33], and references therein.In this work, we propose a derivation of a continuum model from a kinetic versionof the CVA algorithm. For that purpose, we first rephrase the CVA model as a timecontinuous dynamical system (see section 2). Then, we pass to a mean-field version ofthis dynamical system (section 3). This mean field model consists of a kinetic equation ofFokker-Planck type with a force term resulting from the alignement interactions betweenthe particles. More precisely, the mean-field model is written: ε ( ∂ t f ε + ω · ∇ x f ε ) = −∇ ω · ( F ε f ε ) + d ∆ ω f ε + O ( ε ) , (1.1) F ε ( x, ω, t ) = ν (Id − ω ⊗ ω )Ω ε ( x, t ) , (1.2)Ω ε ( x, t ) = j ε ( x, t ) | j ε ( x, t ) | , and j ε ( x, t ) = Z υ ∈ S υ f ε ( x, υ, t ) dυ . (1.3)Here f ε ( x, ω, t ) is the particle distribution function depending on the space variable x ∈ R , the velocity direction ω ∈ S and the time t . d is a scaled diffusion constant and F ε ( x, ω, t ) is the mean-field interaction force between the particles which depends on aninteraction frequency ν . This force tends to align the particles to the direction Ω ε whichis the director of the particle flux j ε . the operators ∇ ω · and ∆ ω are respectively thegradient and the Laplace-Beltrami operators on the sphere. The matrix (Id − ω ⊗ ω ) isthe projection matrix onto the normal plane to ω . ε ≪ ε also equals the ratio of the microscopic time scale (the mean timebetween interactions) to the macroscopic observation time.The ’hydrodynamic limit’ ε → Theorem 1.1 (formal) The limit ε → of f ε is given by f = ρM Ω where ρ = ρ ( x, t ) ≥ is the total mass of f and Ω = Ω( x, t ) ∈ S is the director of its flux: ρ ( x, t ) = Z ω ∈ S f ( x, ω, t ) dω, (1.4)Ω = j | j | , j ( x, t ) = Z ω ∈ S f ( x, ω, t ) ω dω, (1.5) and M Ω is a given function of ω · Ω only depending on ν and d which will be specifiedlater on (see (4.16)). Furthermore, ρ ( x, t ) and Ω( x, t ) satisfy the following system of firstorder partial differential equations: ∂ t ρ + ∇ x · ( c ρ Ω) = 0 . (1.6) ρ ( ∂ t Ω + c (Ω · ∇ )Ω) + λ ( Id − Ω ⊗ Ω) ∇ x ρ = 0 , (1.7) where the convection speeds c , c and the interaction constant λ will be specified in thecourse of the paper (see (4.41) and (4.63)). Hydrodynamic limits have first been developed in the framework of the Boltzmanntheory of rarefied gases. The reader can refer to [8, 11, 31] for recent viewpoints as wellas to [6, 15, 37] for major landmarks in its mathematical theory. Hydrodynamic limitshave been recently investigated in traffic flow modeling [4, 20] as well as in supply chainresearch [3, 14].From the viewpoint of hydrodynamic limits, the originality of theorem 1.1 lies in thefact that the collision operator (i.e. the right-hand side of (1.1)) has a three dimensionalmanifold of equilibria (parametrized by the density ρ and the velocity director Ω) but hasonly a one-dimensional set of collisional invariants (corresponding to mass conservation).Indeed, the interaction does not conserve momentum and one should not expect anycollisional invariant related to that conservation. The problem is solved by introducing abroader class of collisional invariants, such that their integral (with respect to ω ) againstthe collision operator cancels only when the collision operator is applied to a subclass offunctions. Here, a generalized class of collision invariants is associated with each directionΩ on the sphere and the corresponding subclass of functions have their flux in the directionof Ω. We show that such generalized collision invariants exist and that they lead to (1.7).In section 4.4, we show that this system is hyperbolic. The detailed qualitative study ofthe system as well as numerical simulations will be the subject of future work. A summaryof this work can be found in [13]. 3n important consequence of this result is that the large-scale dynamics of the CVAmodel does not present any phase transition, in contrast with the observations of [36].Indeed, the equilibrium is unique (for given density and velocity director). Therefore,the model cannot exhibit any bi-stable behavior where shifts between two competingequilibria would trigger abrupt phase transitions, like in rod-like polymers (see e.g. [24]and references therein). Instead, the equilibrium gradually shifts from a collective onewhere all particles point in the same direction to an isotropic one as the diffusion constant d increases from 0 to infinity. Additionally, the hyperbolicity of the model does not allowlines of faults across unstable elliptic regions, like in the case of multi-phase mixtures orphase transitions in fluids or solids (see e.g. the review in [22] and references therein).With these considerations in mind, a phenomenon qualitatively resembling a phasetransition could occur if the coefficients c , c and λ have sharp variations in some smallrange of values of the diffusion coefficient d . In this case, the model could undergo a rapidchange of its qualitative features which would be reminiscent of a phase transition. Oneof our future goals is to verify or discard this possibility by numerically computing theseconstants.There are many questions which are left open. For instance, one question is about thepossible influence of a limited range of vision in the backwards direction. In this case, theasymetry of the observation will bring more terms in the limit model. Similarly, one couldargue that the angular diffusion should produce some spatial dissipation. Indeed, suchdissipation phenomena are likely to occur if we retain the first order correction in the seriesexpansion in terms of the small parameter ε (the so-called Hilbert or Chapman-Enskogexpansions, see e.g. [11]). A deeper analysis is needed in order to determine the preciseform of these diffusion terms. Another question concerns the possibility of retaining someof the non-local effects in the macroscopic model. It is likely that the absence of phasetransition in the present model is related to the fact that the large-scale limit cancelsmost of the non-local effects (at least at leading order). The question whether retainingsome non-locality effects in the macroscopic limit would allow the appearence of phasetransitions at large scales would indeed reconcile the analytical result with the numericalobservations. A result in this direction obtained with methods from matrix recursionscan be found in [10]. Finally, the alignement interaction is only one of the aspects of theCouzin model, which also involves repulsion at short scales and attraction at large scales.The incorporation of these effects would allow to build a complete continuum model whichwould account for all the important features of this kind of social interaction. The Couzin-Vicsek algorithm considers N point particles in R labeled by k ∈ { , . . . N } with positions X nk at the discrete times t n = n ∆ t . The magnitude of the velocity is thesame for all particles and is constant in time denoted by c >
0. The velocity vector iswritten c ω nk where ω nk belongs to the unit sphere S = { ω s.t. | ω | = 1 } of R .4he Couzin-Vicsek algorithm is a time-discrete algorithm that updates the velocitiesand positions of the particles at every time step ∆ t according to the following rules.(i) The particle position of the k -th particle at time n is evolved according to: X n +1 k = X nk + c ω nk ∆ t. (2.1)(ii) The velocity director of the k -th particle, ω nk , is changed to the director ¯ ω nk of theaverage velocity of the neighboring particles with addition of noise. This algorithm triesto mimic the behaviour of some animal species like fish, which tend to align with theirneighbors. Noise accounts for the inaccuracies of the animal perception and cognitivesystems. The neighborhood of the k -th particle is the ball centered at X nk with radius R > ω nk is given by: ¯ ω nk = J nk | J nk | , J nk = X j, | X nj − X nk |≤ R ω nj . (2.2)In the Couzin-Vicsek algorithm, the space is 2-dimensional and the orientations are vectorsbelonging to the unit sphere S in R . One can write ω nk = e iθ nk with θ nk defined modulo2 π , and similarly ¯ ω nk = e i ¯ θ nk . In the original version of the algorithm, a uniform noise ina small interval of angles [ − α, α ] is added, where α is a measure of the intensity of thenoise. This leads to the following update for the phases: θ n +1 k = ¯ θ nk + ˆ θ nk , (2.3)where ˆ θ nk are independent identically distributed random variables with uniform distri-bution in [ − α, α ]. Then, ω n +1 k = e iθ n +1 k . In [36], Vicsek et al analyze the dynamics ofthis algorithm and experimentally demonstrate the existence of a threshold value α ∗ . For α < α ∗ , a coherent dynamics appears after some time where all the particles are nearlyaligned. On the other hand, if α > α ∗ , disorder prevails at all times.Here, we consider a three dimensional version of the Couzin-Vicsek algorithm, of whichthe two-dimensional original version is a particular case. Of course, formula (2.2) for theaverage remains the same in any dimension. For simplicity, we also consider a Gaussiannoise rather than a uniformly distributed noise as in the original version of the algorithm.Therefore, our algorithm updates the velocity directors according to: ω n +1 k = ˆ ω nk , (2.4)where ˆ ω nk are random variables on the sphere centered at ¯ ω nk with Gaussian distributionsof variance √ D ∆ t where D is a supposed given coefficient. If the Gaussian noise isdiscarded, the evolution of the orientations is given by ω n +1 k = ¯ ω nk , (2.5)where ¯ ω nk is the average defined at (2.2). 5ow, we would like to take the limit ∆ t → | ω nk | = | ω n +1 k | , we have ( ω n +1 k − ω nk )( ω n +1 k + ω nk ) =0. Therefore, defining ω n +1 / k = ( ω n +1 k + ω nk ) / ω n +1 k − ω nk ∆ t = 1∆ t (Id − ω n +1 / k ⊗ ω n +1 / k )(¯ ω nk − ω nk ) , (2.6)where Id denotes the Identity matrix and the symbol ⊗ denotes the tensor product ofvectors. The matrix Id − ω n +1 / k ⊗ ω n +1 / k is the orthogonal projector onto the planeorthogonal to ω n +1 / k . Relation (2.6) simply expresses that the vector ω n +1 k − ω nk belongsto that plane.Now, we let ∆ t →
0. Then, the positions X k ( t ) and the orientations ω k ( t ) becomecontinuous functions of time. If we let ∆ t → ∂ω k /∂t . The right hand side, however, does not seem to have an obvious limit. Thisis due to an improper choice of time scale. Indeed, if we run the clock twice as fast, theparticles will interact twice as frequently. In the limit ∆ t →
0, the number of interactionsper unit of time is infinite and we should not expect to find anything interesting if we donot rescale the time. In order to have the proper time-scale for the model, we need toreplace the tick of the clock ∆ t by a typical interaction frequency ν of the particles underconsideration. For instance, in the case of fish, ν − is the typical time-interval betweentwo successive changes in the fish trajectory to accomodate the presence of other fish inthe neighbourhood. Therefore, we start from a discrete algorithm defined by ω n +1 k − ω nk ∆ t = ν (Id − ω n +1 / k ⊗ ω n +1 / k )(¯ ω nk − ω nk ) , (2.7)together with (2.1) and in the limit ∆ t →
0, we find the following continuous dynamicalsystem: dX k dt = c ω k , (2.8) dω k dt = ν (Id − ω k ⊗ ω k )¯ ω k , (2.9)where we have used that (Id − ω k ⊗ ω k ) ω k = 0. If the Gaussian noise is retained, then, thelimit ∆ t → dX k dt = c ω k , (2.10) dω k = (Id − ω k ⊗ ω k )( ν ¯ ω k dt + √ D dB t ) , (2.11)where dB t is a Brownian motion with intensity √ D . Of course, this ∆ t → ν may depend on the angle between ω k and ¯ ω k , namely ν = ν (cos θ k ), with cos θ k = ω k · ¯ ω k . Indeed, it is legitimate to thinkthat the ability to turn is dependent on the target direction. If we are considering fish,the ability to turn a large angle is likely to be reduced compared to small angles. We willassume that ν (cos θ ) is a smooth and bounded function of cos θ . We now consider the limit of a large number of particles N → ∞ . We first considerthe case without Gaussian noise. For this derivation, we proceed e.g. like in [32]. Weintroduce the so-called empirical distribution f N ( x, ω, t ) defined by: f N ( x, ω, t ) = 1 N N X k =1 δ ( x − X k ( t )) δ ( ω, ω k ( t )) . (3.1)Here, the distribution ω ∈ S → δ ( ω, ω ′ ) is defined by duality against a smooth function ϕ by the relation: h δ ( ω, ω ′ ) , ϕ ( ω ) i = ϕ ( ω ′ ) . We note that δ ( ω, ω ′ ) = δ ( ω − ω ′ ) because the sphere S is not left invariant by thesubtraction operation. However, we have similar properties of δ such as δ ( ω, ω ′ ) = δ ( ω ′ , ω )where this relation is interpreted as concerning a distribution on the product S × S .Then, it is an easy matter to see that f N satisfies the following kinetic equation ∂ t f N + cω · ∇ x f N + ∇ ω · ( F N f N ) = 0 , (3.2)where F N ( x, ω, t ) is an interaction force defined by: F N ( x, ω, t ) = ν (cos θ N ) (Id − ω ⊗ ω )¯ ω N , (3.3)with cos θ N = ω · ¯ ω N and ¯ ω N ( x, ω, t ) is the average orientation around x , given by:¯ ω N ( x, ω, t ) = J N ( x, t ) | J N ( x, t ) | , J N ( x, t ) = X j, | X nj − x |≤ R ω nj . (3.4)If, by any chance, J N is equal to zero, we decide to assign to ¯ ω N ( x, ω, t ) the value ω (whichis the only way by which ¯ ω N ( x, ω, t ) can depend on ω ). In the sequel, this convention willnot be recalled but will be marked by showing the dependence of ¯ ω upon ω We recall the expressions of the gradient and divergence operator on the sphere. Let x = ( x , x , x ) be a cartesian coordinate system associated with an orthonormal basis( e , e , e ) and let ( θ, φ ) be a spherical coordinate system associated with this basis, i.e. x = sin θ cos φ , x = sin θ sin φ , x = cos θ . Let also ( e θ , e φ ) be the local basis associated7ith the spherical coordinate system ; the vectors e θ and e φ have the following coordinatesin the cartesian basis: e θ = (cos θ cos φ, cos θ sin φ, − sin θ ), e φ = ( − sin φ, cos φ, f ( ω ) be a scalar function and A = A θ e θ + A φ e φ be a tangent vector field. Then: ∇ ω f = ∂ θ f e θ + 1sin θ ∂ φ f e φ , ∇ ω · A = 1sin θ ∂ θ ( A θ sin θ ) + 1sin θ ∂ φ A φ . If the cartesian coordinate system is such that e = ¯ ω N , then F N = − ν (cos θ ) sin θ e θ . (3.5)Back to system (3.2)-(3.4), we note that relation (3.4) can be written¯ ω N ( x, ω, t ) = J N ( x, t ) | J N ( x, t ) | , J N ( x, t ) = Z | y − x |≤ R, υ ∈ S υf N ( y, υ, t ) dy dυ . (3.6)We will slightly generalize this formula and consider ¯ ω N ( x, ω, t ) defined by the followingrelation:¯ ω N ( x, ω, t ) = J N ( x, t ) | J N ( x, t ) | , J N ( x, t ) = Z y ∈ R , υ ∈ S K ( | x − y | ) υ f N ( y, υ, t ) dy dυ , (3.7)where K ( | x | ) is the ’observation kernel’ around each particle. Typically, in formula (3.6), K ( | x | ) is the indicator function of the ball centered at the origin and of radius R butwe can imagine more general kernels modeling the fact that the influence of the particlesfades away with distance. We will assume that this function is smooth, bounded andtends to zero at infinity.Clearly, the formal mean-field limit of the particle system modeled by the kineticsystem (3.2), (3.3), (3.7) is given by the following system: ∂ t f + cω · ∇ x f + ∇ ω · ( F f ) = 0 , (3.8) F ( x, ω, t ) = ν (cos ¯ θ ) (Id − ω ⊗ ω )¯ ω ( x, ω, t ) , (3.9)¯ ω ( x, ω, t ) = J ( x, t ) | J ( x, t ) | , J ( x, t ) = Z y ∈ R , υ ∈ S K ( | x − y | ) υ f ( y, υ, t ) dy dυ , (3.10)with cos ¯ θ = ω · ¯ ω . It is an open problem to rigorously show that this convergence holds.For interacting particle system, a typical result is as follows (see e.g. [32]). Suppose thatthe empirical measure at time t = 0 converges in the weak star topology of boundedmeasures towards a smooth function f I ( x, ω ). Then, f N ( x, ω, t ) converges to the solution f of (3.8)-(3.10) with initial datum f I , in the topology of continous functions of timeon [0 , T ] (for arbitrary T >
0) with values in the space of bounded measures endowedwith the weak star topology. We will admit that such a result is true (may be with somemodifid functinal setting) and leave a rigorous convergence proof to future work.We will also admit that the mean-field limit of the stochastic particle system (2.10),(2.11) consists of the following Kolmogorov-Fokker-Planck equation ∂ t f + cω · ∇ x f + ∇ ω · ( F f ) = D ∆ ω f, (3.11)8gain coupled with (3.9), (3.10) for the definition of F and ¯ ω , and where ∆ ω denotes theLaplace-Belltrami operator on the sphere:∆ ω f = ∇ ω · ∇ ω f = 1sin θ ∂ θ (sin θ∂ θ f ) + 1sin θ ∂ φφ f. We are interested in the large time and space dynamics of the mean-field Fokker-Planckequation (3.11), coupled with (3.9), (3.10).So far, the various quantities appearing in the system have physical dimensions. Wefirst introduce the characteristic physical units associated with the problem and scalethe system to dimensionless variables. Let ν the typical interaction frequency scale.This means that ν (cos θ ) = ν ν ′ (cos θ ) with ν ′ (cos θ ) = O (1) in most of the range ofcos θ . We now introduce related time and space scales t and x as follows: t = ν − and x = ct = c/ν . This choice means that the time unit is the mean time betweeninteractions and the space unit is the mean distance traveled by the particles betweeninteractions. We introduce the dimensionless diffusion coefficient d = D/ν . Note that D has also the dimension of a frequency so that d is actually dimensionless. We alsointroduce a scaled observation kernel K ′ such that K ( x | x ′ | ) = K ′ ( | x ′ | ). Typically, if K is the indicator function of the ball of radius R , K ′ is the indicator of the ball of radius R ′ = R/x . The assumption that the interaction is non local means that R ′ = O (1). Inother words, the observation radius is of the same order as the mean distance travelledby the particles between two interactions. This appears consistent with the behaviour ofa fish, but would probably require more justifications. In the present work, we shall takethis fact for granted.Let now t ′ = t/t , x ′ = x/x the associated dimensionless time and space variables.Then, system (3.11), coupled with (3.9), (3.10) is written in this new system of units(after dropping the primes for the sake of clarity): ∂ t f + ω · ∇ x f + ∇ ω · ( F f ) = d ∆ ω f, (4.1) F ( x, ω, t ) = ν (cos ¯ θ ) (Id − ω ⊗ ω )¯ ω ( x, ω, t ) , with cos ¯ θ = ω · ¯ ω, (4.2)¯ ω ( x, ω, t ) = J ( x, t ) | J ( x, t ) | , J ( x, t ) = Z y ∈ R , υ ∈ S K ( | x − y | ) υ f ( y, υ, t ) dy dυ , (4.3)The system now depends on only one dimensionless parameter d and two dimensionlessfunctions which describe the behaviour of the fish: ν (cos ¯ θ ) and K ( x ), which are allsupposed to be of order 1.Up to now, the system has been written at the microscopic level, i.e. at time andlength scales which are characteristic of the dynamics of the individual particles. Our9oal is now to investigate the dynamics of the system at large time and length scalescompared with the scales of the individuals. For this purpose, we adopt new time andspace units ˜ t = t /ε , ˜ x = x /ε with ε ≪
1. Then, a set of new dimensionless variablesis introduced ˜ x = εx , ˜ t = εt . In this new set of variables, the system is written (again,dropping the tildes for clarity): ε ( ∂ t f ε + ω · ∇ x f ε ) = −∇ ω · ( F ε f ε ) + d ∆ ω f ε , (4.4) F ε ( x, ω, t ) = ν ( ω · ¯ ω ε ) (Id − ω ⊗ ω )¯ ω ε ( x, ω, t ) , (4.5)¯ ω ε ( x, ω, t ) = J ε ( x, t ) | J ε ( x, t ) | , J ε ( x, t ) = Z y ∈ R , υ ∈ S K (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) x − yε (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) υ f ε ( y, υ, t ) dy dυ , (4.6)Our goal in this paper is to investigate the formal limit ε → ω ε interms of ε . Lemma 4.1
We have the expansion: ¯ ω ε ( x, ω, t ) = Ω ε ( x, t ) + O ( ε ) , (4.7) where Ω ε ( x, t ) = j ε ( x, t ) | j ε ( x, t ) | , and j ε ( x, t ) = Z υ ∈ S υ f ε ( x, υ, t ) dυ . (4.8)The proof of this lemma is elementary, and is omitted. That the remainder in (4.7) isof order ε is linked with the fact that the observation kernel is isotropic. If an anisotropickernel had been chosen, such as one favouring observations in the forward direction, thena term of order ε would have been obtained. This additional term would substantiallychange the dynamics. We leave this point to future work.The quantity j ε ( x, t ) is the particle flux. We will also use the density, which is definedas a moment of f as well: ρ ε ( x, t ) = Z υ ∈ S f ε ( x, υ, t ) dυ . (4.9)Thanks to lemma 4.1, system (4.4)-(4.6) is written ε ( ∂ t f ε + ω · ∇ x f ε ) = −∇ ω · ( F ε f ε ) + d ∆ ω f ε + O ( ε ) , (4.10) F ε ( x, ω, t ) = ν ( ω · Ω ε ) (Id − ω ⊗ ω )Ω ε ( x, t ) , (4.11)Ω ε ( x, t ) = j ε ( x, t ) | j ε ( x, t ) | , and j ε ( x, t ) = Z υ ∈ S υ f ε ( x, υ, t ) dυ . (4.12)We note that observing the system at large scales makes the interaction local and thatthis interaction tends to align the particle velocity to the direction of the local particle10ux. This interaction term is balanced at leading order by the diffusion term which tendsto spread the particles isotropically on the sphere. Obviously, an equilibrium distributionresults from the balance of these two antogonist phenomena.In the remainder of the paper, we write F [ f ε ] for F ε . We introduce the operator Q ( f ) = −∇ ω · ( F [ f ] f ) + d ∆ ω f, (4.13) F [ f ] = ν (Id − ω ⊗ ω )Ω[ f ] , (4.14)Ω[ f ] = j [ f ] | j [ f ] | , and j [ f ] = Z ω ∈ S ω f dω . (4.15)We note that Ω[ f ] is a non linear operator of f , and so are F [ f ] and Q ( f ). In theremainder, we will always suppose that f is as smooth and integrable as necessary. Weleave the question of finding the appropriate functional framework to forthcoming work.The operator Q acts on the angle variable ω only and leaves the other variables x and t as parameters. Therefore, it is legitimate to study the properties of Q as an operatoracting on functions of ω only. This is the task performed in the following section. Q We begin by looking for the equilibrium solutions, i.e. the functions f which cancel Q .Let µ = cos θ . We denote by σ ( µ ) an antiderivative of ν ( µ ), i.e. ( dσ/dµ )( µ ) = ν ( µ ). Wedefine M Ω ( ω ) = C exp( σ ( ω · Ω) d ) , Z M Ω ( ω ) dω = 1 . (4.16)The constant C is set by the normalization condition (second equality of (4.16)) ; itdepends only on d and on the function σ but not on Ω.We have the following: Lemma 4.2 (i) The operator Q can be written as Q ( f ) = d ∇ ω · (cid:20) M Ω[ f ] ∇ ω (cid:18) fM Ω[ f ] (cid:19)(cid:21) , (4.17) and we have H ( f ) := Z ω ∈ S Q ( f ) fM Ω[ f ] dω = − d Z ω ∈ S M Ω[ f ] (cid:12)(cid:12)(cid:12)(cid:12) ∇ ω (cid:18) fM Ω[ f ] (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) dω ≤ . (4.18) (ii) The equilibria, i.e. the functions f ( ω ) such that Q ( f ) = 0 form a three-dimensionalmanifold E given by E = { ρM Ω ( ω ) | ρ ∈ R + , Ω ∈ S } , (4.19)11 nd ρ is the total mass while Ω is the director of the flux of ρM Ω ( ω ) , i.e. Z ω ∈ S ρ M Ω ( ω ) dω = ρ (4.20)Ω = j [ ρM Ω ] | j [ ρM Ω ] | , j [ ρM Ω ] = Z ω ∈ S ρM Ω ( ω ) ω dω. (4.21) Furthermore, H ( f ) = 0 if and only if f = ρM Ω for arbitrary ρ ∈ R + and Ω ∈ S . The function σ being an increasing function of µ (since ν > M Ω is maximal for ω · Ω = 1, i.e. for ω pointing in the direction of Ω. Therefore, Ω plays the same role as theaverage velocity of the classical Maxwellian of gas dynamics. The role of the temperatureis played by the normalized diffusion constant d : it measures the ’spreading’ of theequilibrium about the average direction Ω. Here the temperature is fixed by the value ofthe diffusion constant, in contrast with classical gas dynamics where the temperature is athermodynamical variable whose evolution is determined by the energy balance equation.An elementary computation shows that the flux can be written j [ ρM Ω ] = h cos θ i M ρ Ω , (4.22)where for any function g (cos θ ), the symbol h g (cos θ ) i M denotes the average of g over theprobability distribution M Ω , i.e. h g (cos θ ) i M = Z M Ω ( ω ) g ( ω · Ω) dω = R π g (cos θ ) exp( σ (cos θ ) d ) sin θ dθ R π exp( σ (cos θ ) d ) sin θ dθ . (4.23)We note that h g (cos θ ) i M does not depend on Ω but depends on d . In particular, h g (cos θ ) i M → g (1) when d → h g (cos θ ) i M → ¯ g , the arithmetic average of g over the sphere, when d → ∞ (with ¯ g = R g ( ω · Ω) dω = R π g (cos θ ) sin θ dθ ). Therefore, h cos θ i M → d → h cos θ i M → d → ∞ . For a large diffusion, theequilibrium is almost isotropic and the magnitude of the velocity tends to zero while fora small diffusion, the distribution is strongly peaked in the forward direction and themagnitude of the velocity tends to 1, which is the velocity of the individual particles. Proof of lemma 4.2:
To prove (i), we introduce a reference frame such that e = Ω[ f ].In spherical coordinates, we have M Ω[ f ] ( ω ( θ, φ )) = C exp( d − σ (cos θ )) . (4.24)Therefore, ∇ ω (ln M Ω[ f ] ) = ∇ ω [ ln { C exp( d − σ (cos θ )) } ]= d − ∇ ω ( σ (cos θ ))= − d − ν (cos θ ) sin θ e θ = d − F [ f ] , (4.25)12here ln denotes the logarithm and the last equality results from (3.5). Then, we deducethat d ∇ ω · (cid:20) M Ω[ f ] ∇ ω (cid:18) fM Ω[ f ] (cid:19)(cid:21) = d ∇ ω · (cid:2) ∇ ω f − f ∇ ω (ln M Ω[ f ] ) (cid:3) = d ∆ ω f − ∇ ω · ( F [ f ] f ) = Q ( f ) . (4.26)(4.18) follows directly from (4.17) and Stokes theorem.(ii) follows directly from (i). If Q ( f ) = 0, then H ( f ) = 0. But H ( f ) is the integralof a non-negative quantity and can be zero only if this quantity is identically zero, whichmeans f = ρM Ω[ f ] for a conveniently chosen ρ . Since Ω[ f ] can be arbitrary, the resultfollows. The remaining statements are obvious.Our task now is to determine the collision invariants of Q , i.e. the functions ψ ( ω ) suchthat Z ω ∈ S Q ( f ) ψ dω = 0 , ∀ f. (4.27)Using (4.17), this equation can be rewritten as Z ω ∈ S fM Ω[ f ] ∇ ω · ( M Ω[ f ] ∇ ω ψ ) dω = 0 , ∀ f. (4.28)Clearly, if ψ = Constant, ψ is a collisional invariant. On the other hand, there is noother obvious conservation relation, since momentum is not conserved by the interactionoperator. The constants span a one-dimensional function space, while the set of equilibriais a three-dimensional manifold. So, we need to find some substitute to the notion ofcollisional invariant, otherwise, in the limit ε →
0, the problem will be under-determined,and in particular, we will lack an equation for Ω (appearing in the expression of theequilibrium).To solve the problem, we slightly change the viewpoint. We fix Ω ∈ S arbitrarily,and we ask the problem of finding all ψ ’s which are collisional invariants of Q ( f ) for all f with director Ω[ f ] = Ω. Such a function ψ is not a collisional invariant in the strict sense,because (4.27) is valid for all f but only for a subclass of f . But this weaker concept ofa collisional invariant is going to suffice for our purpose. So, for fixed Ω, we want to findall ψ ’s such that Z ω ∈ S fM Ω ∇ ω · ( M Ω ∇ ω ψ ) dω = 0 , ∀ f such that Ω[ f ] = Ω . (4.29)Now, saying that Ω[ f ] = Ω is equivalent to saying that j [ f ] is aligned with Ω[ f ], or againto 0 = Ω × j [ f ] = Z ω ∈ S f (Ω × ω ) dω. (4.30)13his last formula can be viewed as a linear constraint and, introducing the Lagrangemultiplier β of this constraint, β being a vector normal to Ω, we can restate the problemof finding the ’generalized’ collisional invariants (4.29) as follows: Given Ω ∈ S , find all ψ ’s such that there exist β ∈ R with Ω · β = 0, and Z ω ∈ S fM Ω {∇ ω · ( M Ω ∇ ω ψ ) − β · (Ω × ω ) M Ω } dω = 0 , ∀ f. (4.31)Now, (4.31) holds for all f without constraint and immediately leads to the followingproblem for ψ : ∇ ω · ( M Ω ∇ ω ψ ) = β · (Ω × ω ) M Ω . (4.32)The problem defining ψ is obviously linear, so that the set C Ω of generalized collisionalinvariants associated with the vector Ω is a vector space. It is convenient to introducea cartesian basis ( e , e , Ω) and the associated spherical coordinates ( θ, φ ). Then β =( β , β ,
0) and β · (Ω × ω ) = ( − β sin φ + β cos φ ) sin θ . Therefore, we can successivelysolve for ψ and ψ , the solutions of (4.32) with right-hand sides respectively equal to − sin φ sin θM Ω and cos φ sin θM Ω .We are naturally looking for solutions in an L ( S ) framework, since ψ is aimed atconstructing marcroscopic quantities by integration against f with respect to ω . There-fore, one possible framework is to look for both f and ψ in L ( S ) to give a meaning tothese macroscopic quantities. We state the following lemma: Lemma 4.3
Let χ ∈ L ( S ) such that R χ dω = 0 . The problem ∇ ω · ( M Ω ∇ ω ψ ) = χ, (4.33) has a unique weak solution in the space ◦ H ( S ) , the quotient of the space H ( S ) by thespace spanned by the constant functions, endowed with the quotient norm. Proof:
We apply the Lax-Milgram theorem to the following variational formulation of(4.33): Z ω ∈ S M Ω ∇ ω ψ · ∇ ω ϕ dω = Z ω ∈ S χϕ dω, (4.34)for all ϕ ∈ ◦ H ( S ). The function M Ω is bounded from above and below on S , so thebilinear form at the left-hand side is continous on ◦ H ( S ). The fact that the average of χ over S is zero ensures that the right-hand side is a continuous linear form on ◦ H ( S ).The coercivity of the bilinear form is a consequence of the Poincare inequality: ∃ C > ∀ ψ ∈ ◦ H ( S ): | ψ | H ≥ C || ψ || ◦ L := C min K ∈ R || ψ + K || L , (4.35)14here | ψ | H is the H semi-norm. We note that the Poincare inequality would not holdwithout taking the quotient.So, to each of the right-hand sides χ = − sin φ sin θM Ω or χ = cos φ sin θM Ω whichhave zero average on the sphere, there exist solutions ψ and ψ respectively (unique upto constants) of problem (4.33). We single out unique solutions by requesting that ψ and ψ have zero average on the sphere: R ψ k dω = 0, k = 1 ,
2. We can state the followingcorollary to lemma 4.3:
Proposition 4.4
The set C Ω of generalized collisional invariants associated with the vec-tor Ω which belong to H ( S ) is a three dimensional vector space C Ω = Span { , ψ , ψ } More explicit forms for ψ and ψ can be found. By expanding in Fourier series withrespect to φ , we easily see that ψ = − g (cos θ ) sin φ, ψ = g (cos θ ) cos φ, (4.36)where g ( µ ) is the unique solution of the elliptic problem on [ − , − (1 − µ ) ∂ µ ( e σ ( µ ) /d (1 − µ ) ∂ µ g ) + e σ ( µ ) /d g = − (1 − µ ) / e σ ( µ ) /d . (4.37)We note that no boundary condition is needed to specify g uniquely since the operatorat the left-hand side of (4.37) is degenerate at the boundaries µ = ±
1. Indeed, it is aneasy matter, using again Lax-Milgram theorem, to prove that problem (4.37) has a uniquesolution in the weighted H space V defined by V = { g | (1 − µ ) − / g ∈ L ( − , , (1 − µ ) / ∂ µ g ∈ L ( − , } . Furthermore, the Maximum Principle shows that g is non-positive.For convenience, we introduce h ( µ ) = (1 − µ ) − / g ∈ L ( − ,
1) or equivalently h (cos θ ) = g (cos θ ) / sin θ . We then define ~ψ ( ω ) = (Ω × ω ) h (Ω · ω ) = ψ e + ψ e . (4.38) ~ψ is the vector generalized collisional invariant associated with the direction Ω. ε → The goal of this section is to prove theorem 1.1.Again, we suppose that all functions are as regular as needed and that all convergencesare as strong as needed. The rigorous proof of this convergence result is outside the scopeof this article.We start with eq. (4.10) which can be written ε ( ∂ t f ε + ω · ∇ x f ε ) = Q ( f ε ) + O ( ε ) . (4.39)15e suppose that f ε → f when ε →
0. Then, from the previous equation, Q ( f ε ) = O ( ε )and we deduce that Q ( f ) = 0. By lemma 4.2, f = ρM Ω , with ρ ≥ ∈ S . Now,since Q operates on the variable ω only, this limit does not specify the dependence of f on ( x, t ), and consequently, ρ and Ω are functions of ( x, t ).To find this dependence, we use the generalized collisional invariants. First, we con-sider the constant collisional invariants, which merely means that we integrate (4.39) withrespect to ω . We find the continuity equation ∂ t ρ ε + ∇ x · j ε = 0 , (4.40)where ρ ε and j ε are the density and flux as defined above. It is an easy matter to realizethat the right-hand side is exactly zero (and not O ( ε )). In the limit ε → ρ ε → ρ and j ε → j = c ρ Ω with c = h cos θ i M , (4.41)and we get ∂ t ρ + ∇ x · ( c ρ Ω) = 0 . (4.42)Now, we multiply (4.39) by ~ψ ε = h ( ω · Ω[ f ε ]) (Ω[ f ε ] × ω ), integrate with respect to ω and take the limit ε →
0. We note that Ω[ f ε ] → Ω and that ~ψ ǫ is smooth enough(given the functional spaces used for the existence theory), and consequently, ~ψ ε → ~ψ = h ( ω · Ω) (Ω × ω ). Therefore, in the limit ε →
0, we get:Ω × X = 0 , X := Z ω ∈ S ( ∂ t ( ρM Ω ) + ω · ∇ x ( ρM Ω )) h ( ω · Ω) ω dω. (4.43)Saying that Ω × X = 0 is equivalent to saying that the projection of X onto the planenormal to Ω vanishes or in other words, that(Id − Ω ⊗ Ω) X = 0 . (4.44)This is the equation that we need to make explicit in order to find the evolution equationfor Ω.Elementary differential geometry gives the derivative of M Ω with respect to Ω actingon a tangent vector d Ω to the sphere as follows: ∂M Ω ∂ Ω ( d Ω) = d − ν ( ω · Ω) ( ω · d Ω) M Ω . (4.45)We deduce that ∂ t ( ρM Ω ) = M Ω ( ∂ t ρ + d − ν ρ ( ω · ∂ t Ω)) , (4.46)( ω · ∇ x )( ρM Ω ) = M Ω (( ω · ∇ x ) ρ + d − ν ρ ω · (( ω · ∇ x )Ω)) . (4.47)16ombining these two identities, we get: ∂ t ( ρM Ω ) + ω · ∇ x ( ρM Ω ) == M Ω (cid:2) ∂ t ρ + ω · ∇ x ρ + d − νρ ( ω · ∂ t Ω + ( ω ⊗ ω ) : ∇ x Ω ) (cid:3) , (4.48)where the symbol ’:’ denotes the contracted product of two tensors (if A = ( A ij ) i,j =1 ,..., and B = ( B ij ) i,j =1 ,..., are two tensors, then A : B = P i,j =1 ,..., A ij B ij ) and ∇ x Ω is thegradient tensor of the vector Ω: ( ∇ x Ω) ij = ∂ x i Ω j . Therefore, the vector X , is given by: X = Z ω ∈ S (cid:2) ∂ t ρ + ω · ∇ x ρ + d − νρ ( ω · ∂ t Ω + ( ω ⊗ ω ) : ∇ x Ω ) (cid:3) ω h M Ω dω (4.49)The four terms in this formula, denoted by X to X , are computed successively usingspherical coordinates ( θ, φ ) associated with a cartesian basis ( e , e , Ω) where e and e aretwo vectors normal to Ω. In the integral (4.49), the functions h = h (cos θ ), ν = ν (cos θ )and M Ω = C exp( σ (cos θ ) d ) only depend on θ . Therefore, the integrals with respect to φ only concern the repeated tensor products of ω .We first have that R π ω dφ = 2 π cos θ Ω, so that X = Z ω ∈ S ∂ t ρ ω h M Ω dω = 2 π ∂ t ρ Z π cos θ h (cos θ ) M Ω (cos θ ) sin θ dθ Ω , (4.50)and (Id − Ω ⊗ Ω) X = 0.Now, an easy computation shows that Z π ω ⊗ ω dφ = π sin θ (Id − Ω ⊗ Ω) + 2 π cos θ Ω ⊗ Ω . (4.51)We deduce that X = Z ω ∈ S (( ω ⊗ ω ) ∇ x ρ ) h M Ω dω == π Z π sin θ h M Ω sin θ dθ (Id − Ω ⊗ Ω) ∇ x ρ ++ 2 π Z π cos θ h M Ω sin θ dθ (Ω · ∇ x ρ ) Ω , (4.52)which leads to:(Id − Ω ⊗ Ω) X = π Z π sin θ h M Ω sin θ dθ (Id − Ω ⊗ Ω) ∇ x ρ, (4.53)Using (4.51) again, we find: X = d − ρ Z ω ∈ S (( ω ⊗ ω ) ∂ t Ω) ν h M Ω dω == πd − ρ Z π sin θ ν h M Ω sin θ dθ (Id − Ω ⊗ Ω) ∂ t Ω ++ 2 πd − ρ Z π cos θ ν h M Ω sin θ dθ (Ω · ∂ t Ω) Ω . (4.54)17he second term at the r.h.s. of (4.54) vanishes since ∂ t Ω is normal to Ω (Ω being a unitvector). For the same reason, (Id − Ω ⊗ Ω) ∂ t Ω = ∂ t Ω and we are left with:(Id − Ω ⊗ Ω) X = πd − ρ Z π sin θ ν h M Ω sin θ dθ ∂ t Ω . (4.55)We now need to compute the integral with respect to φ of the third tensor power of ω . After some computations, we are left with Z π ω ⊗ ω ⊗ ω dφ = π sin θ cos θ ((Id − Ω ⊗ Ω) ⊗ Ω + Ω ⊗ (Id − Ω ⊗ Ω) ++[(Id − Ω ⊗ Ω) ⊗ Ω ⊗ (Id − Ω ⊗ Ω)] :24 )+2 π cos θ Ω ⊗ Ω ⊗ Ω , (4.56)where the index ’: 24’ indicates contraction of the indices 2 and 4. In other words, thetensor element ( R π ω ⊗ ω ⊗ ω dφ ) ijk equals π sin θ cos θ when ( ij, k ) equals any of thetriples (1 , , , , , , , , , , , , π cos θ when ( ij, k ) =(3 , ,
3) and is equal to 0 otherwise. Using Einstein’s summation convention, the followingformula follows:( Z π ω ⊗ ω ⊗ ω dφ ) ∇ x Ω = (cid:18)Z π ω ⊗ ω ⊗ ω dφ (cid:19) ijk ∂ x j Ω k == π sin θ cos θ ((Id − Ω ⊗ Ω) ij Ω k ∂ x j Ω k + Ω i (Id − Ω ⊗ Ω) jk ∂ x j Ω k ++(Id − Ω ⊗ Ω) ik Ω j ∂ x j Ω k )+2 π cos θ Ω i Ω j Ω k ∂ x j Ω k , (4.57)But since Ω is a unit vector, Ω k ∂ x j Ω k = ∂ x j ( | Ω | ) = 0 and the first and fourth terms inthe sum vanish. The expression simplifies into:( Z π ω ⊗ ω ⊗ ω dφ ) ∇ x Ω = π sin θ cos θ ((Id − Ω ⊗ Ω) : ( ∇ x Ω)) Ω ++ π sin θ cos θ (Id − Ω ⊗ Ω)((Ω · ∇ )Ω) , (4.58)The first term is parallel to Ω. Besides, since Ω is a unit vector, (Ω · ∇ )Ω is normal to Ω.So, we finally get(Id − Ω ⊗ Ω)(( Z π ω ⊗ ω ⊗ ω dφ ) ∇ x Ω) = π sin θ cos θ (Ω · ∇ )Ω , (4.59)This leads to the following formula for X :(Id − Ω ⊗ Ω) X = d − ρ (Id − Ω ⊗ Ω) (cid:18)Z ω ∈ S ( ω ⊗ ω ⊗ ω )( ∇ x Ω) ν h M Ω dω (cid:19) = πd − ρ Z π sin θ cos θ ν h M Ω sin θ dθ (Ω · ∇ )Ω (4.60)18ow, we insert the expressions of X to X into (4.44). Using notation (4.23), wefinally find the evolution equation for Ω: d − ρ h sin θ ν h i M ∂ t Ω + d − ρ h sin θ cos θ ν h i M (Ω · ∇ )Ω ++ h sin θ h i M (Id − Ω ⊗ Ω) ∇ x ρ = 0 . (4.61)By the maximum principle, the function h is non-positive. Therefore, we can definesimilar averages as (4.23), substituting M Ω with sin θ ν h M Ω and we denote such averagesas h g i (sin θ ) νhM . With such a notation, (4.61) becomes: ρ ( ∂ t Ω + c (Ω · ∇ )Ω) + λ (Id − Ω ⊗ Ω) ∇ x ρ = 0 , (4.62)with c = h cos θ i (sin θ ) νhM , λ = d (cid:28) ν (cid:29) (sin θ ) νhM (4.63)Collecting the mass and momentum eqs (4.42) and (4.62), we find the final macroscopicmodel of the Couzin-Vicsek algorithm: ∂ t ρ + ∇ x · ( c ρ Ω) = 0 . (4.64) ρ ( ∂ t Ω + c (Ω · ∇ )Ω) + λ (Id − Ω ⊗ Ω) ∇ x ρ = 0 , (4.65)with the coefficients c , c and λ given by (4.41) and (4.63). This ends the proof oftheorem 1.1. The detailed study (both theoretical and numerical) of the properties of the continuummodel (1.4), (1.5), will be the subject of future work. As a preliminary step, we look atthe hyperbolicity of the model.First, thanks to a temporal rescaling, t = t ′ /c , we can replace c by 1, c by c := c /c and λ by λ ′ = λ/c . We will omit the primes for simplicity. Then, the system reads: ∂ t ρ + ∇ x · ( ρ Ω) = 0 . (4.66) ρ ( ∂ t Ω + c (Ω · ∇ )Ω) + λ (Id − Ω ⊗ Ω) ∇ x ρ = 0 , (4.67)This rescaling amounts to saying that the magnitude of the velocity of the individualparticles is equal to 1 /c in the chosen system of units.We choose an arbitrary fixed cartesian coordinate system (Ω , Ω , Ω ) and use sphericalcoordinates ( θ, φ ) in this system (see section 3). Then, Ω = (sin θ cos φ, sin θ sin φ, cos θ ).A simple algebra shows that ( ρ, θ, φ ) satisfy the system ∂ t ρ + ∂ x ( ρ sin θ cos φ ) + ∂ y ( ρ sin θ sin φ ) + ∂ z ( ρ cos θ ) = 0 . (4.68) ∂ t θ + c (sin θ cos φ ∂ x θ + sin θ sin φ ∂ y θ + cos θ∂ z θ ) ++ λ (cos θ cos φ ∂ x ln ρ + cos θ sin φ ∂ y ln ρ − sin θ ∂ z ln ρ ) = 0 . (4.69) ∂ t φ + c (sin θ cos φ ∂ x φ + sin θ sin φ ∂ y φ + cos θ∂ z φ ) ++ λ ( − sin θ sin φ ∂ x ln ρ + sin θ cos φ ∂ y ln ρ ) = 0 . (4.70)19upposing that ρ, θ, φ are independent of x and y amounts to looking at waves whichpropagate in the z direction at a solid angle ( θ, φ ) with the velocity director Ω. In thisgeometry, the system reads: ∂ t ρ + cos θ ∂ z ρ − ρ sin θ ∂ z θ = 0 . (4.71) ∂ t θ + c cos θ ∂ z θ − λ sin θ ∂ z ln ρ = 0 . (4.72) ∂ t φ + c cos θ ∂ z φ = 0 . (4.73)This is a first order system of the form ∂ t ρ∂ t θ∂ t φ + A ( ρ, θ, φ ) ∂ z ρ∂ z θ∂ z φ = 0 , (4.74)with A ( ρ, θ, φ ) = cos θ − ρ sin θ − λ sin θρ c cos θ
00 0 c cos θ , (4.75)The eigenvalues γ ± and γ of the matrix A ( ρ, θ, φ ) are readily computed and are givenby γ = c cos θ, γ ± = 12 h ( c + 1) cos θ ± (cid:0) ( c − cos θ + 4 λ sin θ (cid:1) / i . (4.76)Two special cases are noteworthy. The case θ = 0 (modulo π ) corresponds to waveswhich propagate parallel to the velocity director. In this case, two eigenvalues are equal: γ = γ + = c and γ − = 1. The eigenvectors corresponding to these three eigenvalues arerespectively the density ρ , and the angles θ and φ . So far, the relative magnitude of c and 1 are not known. But, whatever the situation ( c bigger or smaller or even equal to1), the matrix is diagonalizable and therefore the system is hyperbolic.The case θ = π/ π ) corresponds to waves propagating normally to thevelocity director. In this case, γ ± = ± √ λ are opposite and γ = 0. The system for( ρ, θ ) reduces to a special form of the nonlinear wave equation. The sound speed whichpropagates in the medium due to the interactions between the particles has magnitudeequal to 2 √ λ .If θ has an arbitrary value, then, a combination of these two phenomena occurs. Forthe two waves associated with γ ± , there is a net drift at velocity ( c + 1) cos θ and twosound waves with velocities (( c − cos θ +4 λ sin θ (cid:1) / . However, the speed of thewave associated with γ , is not equal to the drift of the two sound waves. A disymmetryappears which is not present in the usual gas dynamics equations. The resolution of theRiemann problem is left to future work. 20 Conclusion
In this paper, we have studied the large-scale dynamics of the Couzin-Vicsek algorithm.For that purpose, we have rephrased the dynamics as a time-continuous one and haveformulated it in terms of a kinetic Fokker-Planck equation. Then, a hydrodynamic scalingof this kinetic equation is introduced with small parameter ε and the limit when ε → ε diffusive corrections, the incorpo-ration of more non-locality effects in the asymptotics and finally, the accounting of theother types of interactions, being of repulsive or attractive type. References [1] M. Aldana and C. Huepe,
Phase transitions in self-driven many-particle systemsand related non-equilibrium models: a network approach , J. Stat. Phys., 112, no 1/2(2003), pp. 135–153.[2] I. Aoki,
A simulation study on the schooling mechanism in fish , Bulletin of the JapanSociety of Scientific Fisheries, 48 (1982), pp. 1081–1088.[3] D. Armbruster, P. Degond and C. Ringhofer,
A model for the dynamics of largequeuing networks and supply chains , SIAM J. Appl. Math., 66 (2006), pp. 896–920.[4] A. Aw, A. Klar, M. Rascle and T. Materne,
Derivation of continuum traffic flowmodels from microscopic follow-the-leader models , SIAM J. Appl. Math., 63 (2002),pp. 259–278.[5] D. R. Brillinger, H. K. Preisler, A. A. Ager, J. G. Kie and B. S. Stewart,
Employingstochastic differential equations to model wildlife motion , Bull Braz Math Soc, 33(2002), pp. 385–408.[6] R. E. Caflisch,
The fluid dynamic limit of the nonlinear Boltzmann equation , Comm.Pure Appl. Math., 33 (1980), pp. 651-666.[7] S. Camazine, J-L. Deneubourg, N. R. Franks, J. Sneyd, G. Theraulaz andE. Bonabeau, Self-Organization in Biological Systems, Princeton University Press,2002.[8] C. Cercignani, R. Illner, M. Pulvirenti, The mathematical theory of dilute gases,Springer-Verlag, New-York, 1991. 219] I. D. Couzin, J. Krause, R. James, G. D. Ruxton and N. R. Franks,
Collective Memoryand Spatial Sorting in Animal Groups , J. theor. Biol., 218 (2002), pp. 1–11.[10] F. Cucker, S. Smale,
Emergent Behavior in Flocks , IEEE Transactions on AutomaticControl, 52 (2007), pp. 852–862.[11] P. Degond,
Macroscopic limits of the Boltzmann equation: a review , in Modeling andcomputational methods for kinetic equations, P. Degond, L. Pareschi, G. Russo (eds),Modeling and Simulation in Science, Engineering and Technology Series, Birkhauser,2003, pp. 3–57.[12] P. Degond, S. Motsch, Large-scale dynamics of the Persistent Turning Walker modelof fish behavior, preprint[13] P. Degond and S. Motsch,
Macroscopic limit of self-driven particles with orientationinteraction , note, to be published.[14] P. Degond, C. Ringhofer,
Stochastic dynamics of long supply chains with randombreakdowns , `a paraˆıtre dans SIAM J. Appl. Math.[15] R. J. DiPerna, P. L. Lions,
On the Cauchy Problem for Boltzmann Equations: GlobalExistence and Weak Stability , The Annals of Mathematics, 130, (1989), pp. 321-366.[16] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi and L. Chayes,
Self-propelled particleswith soft-core interactions: patterns, stability and collapse , Phys. Rev. Lett., 2006.[17] L. Edelstein-Keshet,
Mathematical models of swarming and social aggregation , in-vited lecture, The 2001 International Symposium on Nonlinear Theory and its Ap-plications, (NOLTA 2001) Miyagi, Japan (Oct 28-Nov 1, 2001).[18] J. Gautrais, S. Motsch, C. Jost, M. Soria, A. Campo, R. Fournier, S. Bianco and G.Th´eraulaz,
Analyzing fish movement as a persistent turning walker , in preparation.[19] G. Gr´egoire, and H. Chat´e,
Onset of collective and cohesive motion , Phys. Rev. Lett.,92 (2004) 025702.[20] D. Helbing,
Traffic and related self-driven many-particle systems , Reviews of modernphysics, 73 (2001), pp. 1067–1141.[21] C. Jost et al.,
From individual to collective ant displacements in heterogenous envi-ronments , preprint, 2007.[22] B. L. Keyfitz,
A geometric theory of conservation laws which change type , Zeitschriftfur Angewandte Mathematik und Mechanik, 75, (1995), 571-581.[23] V. L. Kulinskii, V. I. Ratushnaya, A. V. Zvelindovsky, D. Bedeaux,
Hydrodynamicmodel for a system of self-propelling particles with conservative kinematic constraints ,Europhys. Lett., 71 (2005), pp. 207–213.2224] H. Liu, H. Zhang and P.W. Zhang
Axial Symmetry and Classification of StationarySolutions of Doi-Onsager Equation on the Sphere with Maier-Saupe Potential , Comm.Math. Sci. 3 (2), (2005), 201-218.[25] A. Mogilner and L. Edelstein-Keshet,
A non-local model for a swarm , J. Math. Biol.,38 (1999), pp. 534–570.[26] A. Mogilner, L. Edelstein-Keshet, L. Bent and A. Spiros,
Mutual interactions, poten-tials, and individual distance in a social aggregation , J. Math. Biol., 47 (2003), pp.353–389.[27] J. K. Parrish and S. V. Viscido,
Traffic rules of fish schools: a review of agent-basedapproaches , in ’Self-Organization and Complexity’, CK Hemelrijk (ed.), CambridgeUniversity Press, 2003.[28] J. K. Parrish, S. V. Viscido and D. Gr¨unbaum,
Self-organized fish schools: an exam-ination of emergent properties , The biological bulletin, 202 (2002), pp. 296–305.[29] V. I. Ratushnaya, D. Bedeaux, V. L. Kulinskii and A. V. Zvelindovsky,
Collectivebehaviour of self propelling particles with kinematic constraints ; the relations betweenthe discrete and the continuous description , Physica A, to appear.[30] V. I. Ratushnaya, V. L. Kulinskii, A. V. Zvelindovsky, D. Bedeaux,
Hydrodynamicmodel for the system of self propelling particles with conservative kinematic con-straints; two dimensional stationary solutions
Physica A, 366, (2006), pp. 107–114.[31] Y. Sone, Kinetic Theory and Fluid Dynamics, Birkhauser, 2002.[32] H. Spohn, Large scale dynamics of interacting particles, Springer, Berlin, 1991.[33] Theraulaz et al.,
Spatial patterns in ant colonies , Proceedings of the NationalAcademy of Sciences, 99 (2002), pp. 9645–9649.[34] C. M. Topaz and A. L. Bertozzi,
Swarming patterns in a two-dimensional kinematicmodel for biological groups , SIAM J. Appl. Math, 65 (2004), pp. 152–174.[35] C. M. Topaz, A. L. Bertozzi, M. A. Lewis,
A nonlocal continuum model for biologicalaggregation , Bull. Math. Biol., 68 (2006), pp. 1601–1623.[36] T. Vicsek, A. Czir´ok, E. Ben-Jacob, I. Cohen and O. Shochet,
Novel type of phasetransition in a system of self-driven particles , Phys. Rev. Lett., 75 (1995), pp. 1226–1229.[37] S-H. Yu,