Emmanuel Boissard
email: emmanuel.boissard@math.univ-toulouse.fr

Pierre Degond
email: pierre.degond@math.univ-toulouse.fr

Sebastien Motsch
email: smotsch@cscamm.umd.edu

Trail formation based on directed pheromone deposition

Keywords: Self-propelled particles, pheromone deposition, directed pheromones, alignment interaction, Individual-Based Model, trail detection, pattern formation, kinetic models, fluid models AMS Subject classification: 35Q80, 35L60, 82C22, 82C31, 82C70, 82C80, 92D50

We propose an Individual-Based Model of ant-trail formation. The ants are modeled as self-propelled particles which deposit directed pheromones and interact with them through alignment interaction. The directed pheromones intend to model pieces of trails, while the alignment interaction translates the tendency for an ant to follow a trail when it meets it. Thanks to adequate quantitative descriptors of the trail patterns, the existence of a phase transition as the ant-pheromone interaction frequency is increased can be evidenced. Finally, we propose both kinetic and fluid descriptions of this model and analyze the capabilities of the fluid model to develop trail patterns. We observe that the development of patterns by fluid models require extra trail amplification mechanisms that are not needed at the Individual-Based Model level.

Introduction

One of the many features displayed by self-organized collective motion of animals or individuals is the formation of trails. For instance, ant displacements are characterized by their organization into lanes consisting of a large number of individuals, for the purpose of exploring the environment or exploiting its resources. Another example involving species with higher cognitive capacities is the formation of mountain trails by hikers or herds of animals. In both cases, the main feature is that the interaction between the individuals is not direct, but instead, is mediated by a chemical substance or by the environment. Indeed, ants lay down pheromones as chemical markers. These pheromones are sensed by other individuals which use them to adjust their path. In the case of mountain trails, the modification of the soil by walkers facilitates the passage of the next group of individuals and attracts them. This phenomenon is well-known to biologists under the name of stigmergy, a concept first forged by Pierre-Paul Grassé [START_REF] Grassé | Termitologia: Comportement, socialité, écologie, evolution, systématique[END_REF] to describe the coordination of social insects in nest building.

The formation of trails by ants has been widely studied in the biological literature [START_REF] Beckers | Trail laying behaviour during food recruitment in the antLasius niger (L.)[END_REF][START_REF] Couzin | Self-organized lane formation and optimized traffic flow in army ants[END_REF][START_REF] Detrain | The influence of the physical environment on the self-organised foraging patterns of ants[END_REF][START_REF] Edelstein-Keshet | Simple models for trail-following behaviour; trunk trails versus individual foragers[END_REF][START_REF] Edelstein-Keshet | Trail following in ants: individual properties determine population behaviour[END_REF][START_REF] Watmough | Modelling the formation of trail networks by foraging ants[END_REF]. One general observation is the fact that trail formation is a self-organized phenomenon and expresses the emergence of a large-scale order stemming from simple rules at the individual level. Indeed, ant colonies in the numbers of thousands of individuals or more arrange into lines without resorting to long-range signaling or hierarchical organization. Another striking feature is the variability of the trail patterns, which may range from densely woven networks to a few large trails. This flexibility may result from the ability of the individuals to adapt their activity to variable external conditions such as food availability, temperature, terrain conditions, the presence of predators, etc. Trail plasticity derives from internal and external factors: for example it may vary according to the species of ant under consideration or depends on the properties of the soil. Our goal is to quantities).

All these previously cited two-dimensional models assume a pre-existing lattice structure. One questions which is seldom addressed is whether this pre-existing lattice may influence the formation of the trails. For instance, it is well-known that lattice Boltzmann models with too few velocities have incorrect behavior. One may wonder if similar effects could be encountered with spatially discrete ant trail formation models. For this reason, in the present work, we will depart from a latticelike spatial organization and treat the motion of the ants in the two-dimensional continuous space.

In this work, we propose a time and space continuous Individual-Based Model for self-organized trail formation. In this model, self-propelled particles interact by laying down pheromone trails that indicate both the position and direction of the trails. Ants adapt their course by following trails deposited by others, therefore reinforcing existing trails while evaporation of pheromones allows weaker trails to disappear. The ant dynamics is time-discrete and is a succession of free flights and velocity jumps occuring at time intervals ∆t. Velocity jumps occur with an exponential probability. Two kinds of velocity jumps are considered: purely random jumps which translate the ability of ants to explore a new environment, and trailrecruitment jumps. In order to perform the latter, ants look for trails in a disk around themselves, pick up one of these trails with uniform probability and adopt the direction of the chosen trail.

This model bears analogies with chemotaxis models. Chemotaxis is the name given to remote attraction interaction through chemical signaling in colonies of bacteria. Mathematical modeling of chemotaxis has been largely studied. Macroscopic models were first introduced by Keller and Segel in the form of a set of parabolic equations [START_REF] Keller | Model for chemotaxis[END_REF]. These equations can be obtained as macroscopic limits of kinetic models [START_REF] Erban | From individual to collective behavior in bacterial chemotaxis[END_REF][START_REF] Filbet | Derivation of hyperbolic models for chemosensitive movement[END_REF][START_REF] Hillen | The diffusion limit of transport equations derived from velocity-jump processes[END_REF][START_REF] Othmer | The diffusion limit of transport equations II: chemotaxis equations[END_REF][START_REF] Othmer | Aggregation, blowup, and collapse: The ABC's of taxis in reinforced random walks[END_REF]. Kinetic models describe the evolution of the population density in position-velocity space. In [START_REF] Stevens | The derivation of chemotaxis equations as limit dynamics of moderately interacting stochastic many-particle systems[END_REF] a direct derivation of the Keller-Segel model from a stochastic many-particle model is given. The common feature of most chemotaxis models is the appearance of blow-up, which corresponds to the fast aggregation of individuals at a specific point in space (see e.g. [START_REF] Blanchet | Two-dimensional Keller-Segel model: optimal critical mass and qualitative properties of the solutions[END_REF][START_REF] Calvez | Volume effects in the Keller-Segel model: energy estimates preventing blow-up[END_REF]). By contrast, in the present paper, the dynamics gives rise to the spontaneous organization of lane-like spatial patterns, much alike to the observed behavior of ants. The reason for this different morphogenetic behavior is the directed nature of the mediator of the interaction.

The model also bears analogies with the kinetic model of cell migration developed in [START_REF] Painter | Modelling cell migration strategies in the extracellular matrix[END_REF]. In this model, cells move in a medium consisting of interwoven extra-cellular fibers in the direction of one, randomly chosen fiber direction. As they move, cells specifically destroy the extracellular fibers which are transverse to their motion. The induced trail reinforcement mechanism produces a network.

There are two differences with the trail formation mechanism that we present here. The first one is that the dynamics starts from a prescribed set of motion directions and gradually reduces this set, while our algorithms builds up the set of available directions gradually and new directions are created through random velocity jumps. The second difference is the role of trail evaporation in our algorithm, which has no equivalent in the cell motion model. Indeed, trail evaporation is a major ingredient for network plasticity.

In the last part of the present work, we derive a kinetic formulation of the proposed ant trail formation model in the spirit of [START_REF] Painter | Modelling cell migration strategies in the extracellular matrix[END_REF]. Then, the fluid limit of this kinetic model is considered. We show that the resulting fluid model can exhibit trail formation only if some concentration mechanism is involved, while numerical simulations indicate that trail formation may occur without such a mechanism. Therefore, the appearance of trails is enhanced when the model provides more information about the ant velocity distribution function.

The outline of this paper is as follows. In section 2, we provide the model description. Section 3 is devoted to the analysis of the numerical simulations. We establish a methodology for the detection of trail patterns from a simulation outcome, and analyze the dependency of the observed features on the model parameters. In section 4, we formally establish a set of kinetic equations that describes the dynamics and we investigate their fluid limit. A conclusion in section 5 draws some perspectives of this work.

An Individual-Based-Model of ant behavior based on directed pheromone deposition

We consider N "ants" in a flat (2-dimensional) domain: each ant is described by its position x i ∈ R 2 and the direction of its motion ω i . The vector ω i is supposed to be of unit-length, i.e. ω i ∈ S 1 , where S 1 denotes the unit circle. We also consider pieces of trails described by pairs (y p , ω p ) where y p ∈ R 2 is the trail piece position and ω p is a unit vector describing the trail direction (see figure 1). In the case of ants, the marking of the trail is realized by a chemical marker, namely a trail pheromone. We assume that the ants can distinguishably perceive the direction of the trail of this chemical marker and that they are able to follow, not the line of steepest gradient, like in chemotaxis, but the direction of this trail. Note that in the case of walkers or sheep in an outdoor terrain, the marking of the trail is realized by the modification of the terrain consecutive to the passage of the walkers, such as flattened grass. For wild white bears, this modification is realized by the trail left in the snow by the animals. In the sequel, we concentrate on the modeling of ant trail formation, and we will indistinguishably refer to these pieces of trails as 'trails' or 'trail pheromones', or simply, 'pheromones'. The set of pheromones varies with time, since new pheromones are created by the deposition process and pheromones disappear after some time in order to model the evaporation process.

We will denote by P(t) the set of pheromones at time t.

The simulated ants follow a random walk process. During free flights, Ant i moves in direction ω i at a constant speed c, i.e. is subject to the differential equation: ẋi = cω i , ωi = 0, (

where the dots stand for time derivatives. This free motion is randomly interrupted by velocity jumps. When Ant i undergoes a velocity jump at time t, its velocity direction before the jump ω i (t-0) is suddenly changed into a different one ω i (t+ 0). The jump times are drawn according to Poisson distributions. In practice, a time discrete algorithm is used, with time steps ∆t. With such a discretization, a Poisson process of frequency λ is represented by an event occuring with probability 1-e -λ∆t over this time step. There are two kinds of jumps: random velocity jumps and trail recruitment jumps.

Random velocity jumps. In this case, ω i (t+0) differs from ω i (t-0) by a random angle ε, i.e. (ω i (t -0), ω i (t + 0)) = ε, where ε is drawn out of a Gaussian distribution p(ε) with zero mean and variance σ 2 , periodized over [0, 2π], i.e.

p(ε) dε

= n∈Z 1 (2πσ 2 ) 1/2 exp - (ε + 2nπ) 2 2σ 2 dε, ε ∈ [0, 2π].
The frequency of the Poisson process is constant in time and denoted by λ r .

Trail recruitment jump. In this case, ω i (t+0), is picked up with uniform probability among the directions ω p of the trail pheromones located in the ball B R (x i (t)) of radius R centered at x i (t) (see figure 1). B R (x i (t)) is the ant detection region and R its detection radius. More precisely, defining the set

S i (t) = {p ∈ P(t) , |x p -x i (t)| ≤ R},
Ant i chooses an index p in S i (t) with uniform probability and sets

ω i (t + 0) = ω p . (2.2)
A variant of this mechanism involves a nematic interaction (i.e. the deposited trails have no specific orientation). In this case, the new direction is defined by

ω i (t + 0) = ±ω p , such that (ω i (t + 0), ω p ) is acute. (2.
3)

The nematic interaction makes more biological sense, since it seems difficult to envision a mechanism which would allow the ants to detect the orientation of a given

(x i , ω i ) R (x p , ω p )
Figure 1: Ants follow a random walk process. Each ant is moving on a straight line until it undergoes a random velocity jump (left) or a trail recruitment jump (right).

In this picture, a trail pheromone is located in the disk centered at the ant location when the jump occurs, and of given radius R.

trail. The use of a uniform probability to select the interacting pheromone can be questioned. For instance, the choice of the trail pheromone p could be dependent on the angle (ω i (t), ω p ) like considered in [START_REF] Vincent | The effect of a non-uniform turning kernel on ant trail morphology[END_REF], but in the present work we will discard this effect. Ants may also preferably choose the largest trails indicated by a large concentration of pheromones in one given direction. This 'preferential choice' will be discussed in connection to the kinetic and fluid models in section 4 but discarded in the numerical simulations of the Individual-Based Model. The frequency of the Poisson process is given by λ p M i (t) where M i (t) is the number of pheromones in the detection region of Ant i: M i (t) = Card(S i (t)), and λ p is the trail recruitment frequency per unit pheromone. The dependency of the jump frequency upon the number of detected pheromones accounts for the observed increase of the alignment probability with the pheromone density. Of course, nonlinear functions of the pheromone density could be chosen as well. For instance, some saturation of the detection capability occurs at large pheromone densities, such as investigated in [START_REF] Rauch | Pattern formation and functionality in swarm models[END_REF]. This effect will also be discarded here. We also discard any consideration of the detection mechanism, such as discussed e.g. in [START_REF] Calenbuhr | A model for osmotropotactic orientation (II)*[END_REF][START_REF] Calenbuhr | A model for osmotropotactic orientation (I)*[END_REF][START_REF] Couzin | Self-organized lane formation and optimized traffic flow in army ants[END_REF].

During their walk, ants leave trail pheromones at a certain deposition rate ν d . If at time t, ant i deposits a pheromone, a new pheromone particle is created at position x i (t) with the direction ω i (t). Hence, we postulate that: At deposition times t, a pheromone p is created with (x p , ω p ) = (x i (t), ω i (t)).

Pheromones have a life-time T p and remain immobile during their lifetime. In this work, pheromone diffusion is neglected. Pheromone deposition and evaporation times are modeled by Poisson processes: each ant has a probability ν d per unit of time to lay down a pheromone and each pheromone has a probability 1/T p per unit of time to disappear.

Pheromone deposition mediates the interactions between the ants. This interaction is nonlocal in both space and time (because the ant which has deposited a pheromone may have moved away quite far before another ant interacts with it). Random velocity jumps and trail recruitment jumps have opposite effects. Random velocity jumps generate diffusion at large scales whereas trail recruitment jumps tend to produce concentrations of the ants trajectories on the pheromone trails. Therefore, the pheromone-meditated interaction induces correlations of the ants motions and these correlations result in trail formation.

The trail recruitment process together with pheromone evaporation result in network plasticity. To illustrate this mechanism, let us consider the following simplified situation. Suppose an ant reaches a "crossroads" of trails, meaning a spot where pheromones point in two different directions denoted by 1 and 2. Suppose there are n 1 pheromones in one direction and n 2 in the other one. The probability for the ant to choose to orient in direction i; i = 1, 2, is equal to the ratio n i n 1 +n 2 . When the ant turns to its new direction, it may release a pheromone which will serve to reinforce this branch. Eventually, one under-selected branch of the crossroads will vanish due to evaporation of the pheromones. Note that the choice of the surviving branch depends on random fluctuations of this process: therefore, the outcome of this situation is non-deterministic and even an initially strongly populated branch has a non-zero probability to vanish away.

Simulations and results

Choice of the modeling and numerical parameters

We use experimentally determined parameter values as often as possible. Since parameters are species-dependent, we focus on the species 'Lasius Niger'.

In our model, the motion of a single ant is described by three quantities: the speed c, the frequency of random velocity jumps λ r and their amplitude σ. These three parameters have been estimated in different studies [START_REF] Bernadou | Does substrate coarseness matter for foraging ants? an experiment with lasius niger (Hymenoptera; formicidae)[END_REF][START_REF] Casellas | From individual to collective displacements in heterogeneous environments[END_REF] which give us a range of possible values. We choose rather low estimations of λ r and σ (see Table 1) since in real experiments the estimation of these coefficients counts both for random jumps and recruitment by trails. The deposition rate of pheromones ν d and their life time T p have also been measured experimentally for 'Lasius Niger' [START_REF] Beckers | Trail laying behaviour during food recruitment in the antLasius niger (L.)[END_REF]. After leaving a food source, an ant drops on the average .5 pheromone per second. This experimental value gives us an upper bound for ν d because it corresponds to an estimation of ν d in a very specific situation where the ant activity level is high. In our simulations, we use pheromones, we also take a low estimation of the pheromone lifetime (T p = 100 s) otherwise the domain becomes saturated with pheromones. By contrast, the interaction between ants and pheromones has not been quantified experimentally. For this reason, we do not have experimental values for the pheromone detection radius R and the alignment probability per unit of time λ p . In our simulations, we fix the radius of perception R equal to 1 cm (corresponding roughly to 2 body lengths). The alignment probability λ p remains a free parameter in our model. By changing the value of λ p , we can tune the influence of the pheromone-mediated interaction between the ants. A low value of λ p corresponds to a weak interaction, whereas, for large values of λ p , the ant velocities become controlled by the pheromone directions. In our simulations, λ p varies from 0 to 3 s -1 .

ν d = .2 s -1 .
For simplicity, all simulations are carried out in a square domain of size L = 100 cm with periodic boundary conditions. For the initial condition, 200 ants are randomly distributed in the domain. Their velocity ω i is chosen uniformly on the circle S 1 . The ant-pheromone interaction is always taken nematic unless otherwise stated.

We can estimate the average number of pheromones M at equilibrium, when the average is taken over realizations. The evolution of M (t) is given by the following differential equation:

d M (t) dt = ν d N - 1 T p M (t),
where N is the number of ants, ν d and T p are (resp.) the deposition rate and pheromone lifetime. Thus, at equilibrium (i.e. d M /dt = 0), the average number of pheromones is given by ν d T p N. For our choice of parameters (Table 1), this corresponds to 4000 pheromones.

Detection of trails

Evidence of trail formation

The typical outcomes of the model are shown in figure 2. After an initial transient, we observe the formation of a network of trails. This network is not static, as we observe in the two graphics: the network at time t = 2000 s is significantly different from the network observed at time t = 1000 s. Here, the goal is to provide statistical descriptors of this trail formation phenomenon and to analyze it. , ∆t = .05 (see also Table 1).

Definition of a trail

To quantify the amount of particles that are organized into trails at a given time, we consider the collection of all particles, that is to say, the union of the sets of ants and of pheromones. Indeed collecting the pheromones allows us to trace back the recent history of the individuals. To define a trail, we fix two parameters: a distance r max and an angle θ max . We say that particle P i = (x i , ω i ) (P i being either an ant or a pheromone) is linked to particle P j = (x j , ω j ) if the distance between the two particles is less than r max and the angle between ω i and ω j is either less than θ max or greater than πθ max . In other words, we define a relationship (see figure 3):

P i ∼ P j if and only if |x i -x j | < r max and sin(ω j -ω i ) < sin θ max . (3.1)
Using this relationship, particles can be sorted into different trails: we say that P and Q belong to the same trail if there exists particles P 1 , . . . , P k (a path) such that P ∼ P 1 , P 1 ∼ P 2 , . . ., P k ∼ Q. Thus, a trail is defined as the connected components of the particles under the relationship (3.1). A trail is approximately a slowly turning lane of particles. As a first example, in figure 4, we display the partitioning into trails of the previous simulations (figure 2) at time t = 2000 s with r max = 2 cm and θ max = 45 • . For these values, the largest trail (drawn in red in figure 4) consists of 2670 particles and the second largest (drawn in orange in figure 4) is made of 254 particles.

Statistics of the trails

We expect that trail formation results in the development of a small number of large trails, while unorganized states are characterized by a large number of small trails, most of them being reduced to single elements. Therefore, trail formation can be detected by observing the trail sizes. With this aim, we denote by S i (t) the size of the trail to which particle i belongs at time t. Let N (t) be the total number of particles, i.e. N (t) = N + P(t) where N is the number of ants and P(t) is the number of pheromones at time t. We form

p t (S) = Card({i | S i (t) = S}) N (t) , S ∈ N.
p t (S) is the probability that a particle belongs to a trail of size S at time t. An unorganized state is therefore characterized by a quickly decaying p t (S) as a function of S while a state where the particles are highly organized into trails displays a bimodal p t (S) with high values for large values of S. To display the distribution of p t (S) is easy: it is nothing but the histogram of the trail sizes S i , collected from several independent simulations with identical parameters.

As an illustration, we provide the distribution p t (S) for the set parameters used to generate Figs. 2 and4, (i.e. λ p = 2 s -1 , r max = 2 cm and θ max = 45 • ), with 1000 realizations. We clearly observe in Fig. 5 (left) that the distribution p t (S) is bimodal: a first maximum is observed near the minimal value of S, i.e. S = 1, and a second maximum is observed near the values S ≈ 2500. This indicates that a particle (i.e. an ant or a pheromone) belongs to either a small-size trail (S < 100) or to a large-size trail (S ≈ 2500). As a control sample for our statistical measurement, we run the same simulations but cutting off the ant-pheromone interaction (i.e. λ p = 0) and proceed to the same analysis. In Fig. 5 (right), we observe that without the influence of the pheromones (blue histogram) the probability p t (s) is only concentrated near the value S = 1 and decays very fast to almost vanish for S > 500. (Right) Histograms of S with and without trail-pheromone interaction (λ p = 2 and λ p = 0 resp.). The parameters for this simulation are the same as in Figs. 2 and4.

Trail size

As observed in figure 5, ant-pheromone interactions lead to the formation of trails which are evidenced by the transformation of the shape of the distribution p t (S).

To perform a systematic parametric analysis of the trail formation phenomenon, we use the mean S of the distribution S:

S = S∈N S p t (S).
The quantity S quantifies the level of organization of the system into trails. Indeed, large values of S indicate a high level of organization into trails while smaller values of S are the signature of a disordered system. For example, in Fig. 5, we have S = 1333.7 when the ant-pheromone interaction is on with interaction frequency λ p = 2. By contrast, its value falls down to S = 76.8 when the antpheromone interaction is turned off (i.e. λ p = 0) and the system is in a fully disordered state.

Our first use of the mean trail size S is to show that it stabilizes to a fixed value after an initial transient. Fig 6 shows the mean trail size S (t) for one simulation (dashed line) and averaged over 1000 different simulations (solid line). It appears that, after some transient, S (t) presents a lot of fluctuations about an averaged value. If the simulation is reproduced a large number of times and the mean trail size S (t) is averaged over all these realizations, the convergence towards a constant value becomes apparent.

Therefore, statistical analysis of the trail patterns using the mean trail size S become significant only once this constant value has been reached. In the forthcoming sections, analysis will be performed for simulation times equal to 2000 s, which is significantly larger than the time needed for the stabilization of S (about 800 s).

Evidence of a phase transition

In Fig. 7, we display S as a function of λ p . For each value of λ p , we estimate S by averaging it over 1000 independent simulations. We observe an abrupt increase of S when λ p varies from 0 to 1 which means that a sharp transition from an unorganized system to a system organized into trails arises. For larger values of λ p , the influence of the ant-pheromone interaction saturates and S reaches a plateau at the approximate value S ≈ 1300. The transition from disorder to trails also depends on the other parameters of the model. For example, if we increase the noise by increasing the random jump frequency λ r , the corresponding value of S decreases. In order to restore the previous value of S the ant-pheromone interaction frequency λ p must be increased simultaneously. In Fig. 8, we plot S as a function of both the random jump frequency λ r and the ant-pheromone interaction frequency λ p . We estimate S by averaging it over 100 realizations for each value of the pair (λ p , λ r ). We still observe a fast transition from an unorganized state ( S < 100) to a state organized into trails ( S ≥ 1000) when λ p increases. However, as the noise λ r increases, this transition becomes smoother. Moreover, the plateau reached by S when λ p is large is still comprised between 1300 and 1500 for all values of λ r , but reaching this plateau for large values of λ r requires larger value of λ p .

At the value λ r = 0, the transition from disorder ( S ≤ 100) to trail-like organization ( S ≥ 1000) is the fastest. However, the plateau reached by S when λ p is large is significantly lower than for larger values of λ r ( S ≈ 1000 instead of 1300). This could be attributed to the fact that, without random jumps, the level of diffusion is too low, the ants do not mix enough, and trails have little opportunities to merge.

On the other hand, we can look for another explanation of this paradoxical lower value of S when λ r is very small. Indeed, we notice that, in this case, the formed trails are much narrower than for larger values of λ r . Fig. 9 (left) shows a simulation result using a quite small random jump frequency of λ r = .2. We observe that the trails are narrower and more straight than those obtained with the larger value λ r = 2 (figure 2). We can quantify statistically this feature by changing the parameters of trail detection r max and θ max . We reduce the maximum distance (r max = 1.5 cm) and the maximum angle (θ max = 35 • ). With these smaller values, two particles are less likely to be connected. Then we proceed to the same analysis as in figure 8, by estimating the mean size of the trails S as a function of λ r and λ p , averaged over 100 realizations. As we observe in figure 9 (right), the mean size of the trails S is much larger for smaller values of λ r and we recover the same behavior as that observed for larger values of λ r . This discussion illustrates the difficulty of working with an estimator which depends on arbitrary choices of scales (here the space and angular threshold of trail detection). A discussion of the dependence of the trail width upon the biological parameters is developed in the next section.

Trail width

A way to highlight the dependence of the trail width upon the model parameters is to compute a two-particle correlation distribution. Let a particle (ant or pheromone) i be located at position x i and velocity ω i . Denote by ω ⊥ i the orthogonal vector to ω i in the direct orientation. For all particles j = i, we form the vector where δ is the Dirac delta, provides the probability that, given a first particle (located at say x 0 with orientation ω 0 ), a second particle lies at location x 0 +ω ⊥ 0 X x + ω 0 X y (see figure 10 (left) for an illustration of the construction of f 2 ). Looking at this 2-particle density, trails appear as concentrations near a line passing through the origin and directed in the y-direction. Figure 10 (right) provides a histogram of the two-particle density f 2 for the simulation corresponding to the right picture of fig. 2. The above mentioned concentration is clearly visible. Additionally, the typical width of this concentration gives access to the typical width of the trails.

X ij = (x j -x i ) • ω ⊥ i (x j -x i ) • ω i . The distribution 2 N (N -1) f 2 (X), with f 2 (X) = f 2 (X x , X y ) = (i,j), i =j δ(X -X ij ),
In order to better estimate the typical width of the trails, we plot cuts of the two-particle density f 2 along the line {y = 0} (see figure 10). In practice, these cuts are determined by computing the following density

f2 (r) = (i,j), i =j, |(X ij )y|≤ξ δ(r -|(X ij ) x |).
where ξ is suitable chosen (of the order of 1 cm). Figure 11 (left) displays f2 (r) as a function of r for different values of the trail recruitment frequency λ p and a fixed value of the random jump frequency equal to λ r = 2 s -1 . It appears that f2 is higher decreases faster for larger values of λ p . The decay of f2 can give an estimate of the width of the trail: if we approximate the decay of f2 by an exponential, f2 (r) ≈ f 0 exp -r r 0 for r ≈ 0, then r 0 measures the typical width of the trail. This quantity can be estimated using the formula:

r 0 = 1 |(ln f2 ) ′ (0)| . (3.2)
As we observe in Table 2, the width r 0 increases as λ p decreases. Therefore, increasing the trail recruitment frequency increases the intensity of the particles interactions and produces trails with smaller width. Figure 11 (right) displays f2 (r) as a function of r for different values of the random jump frequency λ r and a fixed value of the trail recruitment frequency equal to λ p = 2 s -1 . Here, the trail width r 0 estimated from f2 is larger for large values of λ r (see Table 2), indicating that the typical width of the trails increases with increasing λ r , as it should.

We also observe a discontinuity at r = 0 for all the functions f2 (figure 11). These jumps are easily explained by the deposit process: each time an ant drops a pheromone, the new pheromone and the ant are located at the same position exactly. This results in a peak of concentration of f2 at r = 0. 11.

λ p r 0 (
We estimate the derivative of ln f2 (r) near 0 using the values of r between .1 and 2.

Kinetic and continuum descriptions 4.1 Framework

In this section, we propose meso-and macro-scopic descriptions of the previously discussed ant dynamics. We first propose a kinetic model, i.e. a model for the probability distributions of ants and pheromones. The derivation of this kinetic model is formal and based on analogies with the underlying discrete dynamics. A rigorous derivation of the kinetic model from the discrete dynamics is up to now beyond reach. Issues such as the validity of the chaos propagation property [START_REF] Cercignani | The mathematical theory of dilute gases[END_REF], which is the key for proving such results, may be quite difficult to solve. Then, fluid limits of this kinetic model will be considered. We will notice that the resulting fluid models can only exhibit the development of trails if some concentration mechanism is added, while the numerical simulations above indicate that such a mechanism is not needed at the level of the Individual-Based Model.

Kinetic model

In this section, we introduce the kinetic model of the discrete ant-pheromone interaction on a purely formal basis. We introduce the ant distribution function F (x, ω, t) and the pheromone distribution function G(x, ω, t), for x ∈ R 2 , ω ∈ S 1 and t ≥ 0. They are respectively the number density in phase-space (x, ω) of the ants (respectively of the pheromones), i.e. the number of such particles located at position x with orientation ω at time t. Here, we remark that the consideration of directed pheromones requires the introduction of a pheromone density in (position, orientation) phase-space in the same manner as for the ants.

Trail dynamics. The trail dynamics is described by the ordinary differential equation:

∂ t G(x, ω, t) = ν d F (x, ω, t) -ν e G(x, ω, t), ( 4.1)
This equation can be easily deduced from the evolution of the probability density of the underlying stochastic Poisson process. The first term describes deposition by the ants according to a Poisson process of frequency ν d while the second term results from the finite lifetime expectancy T p = ν -1 e of the pheromones. Pheromones are supposed immobile, which explains the absence of any convection or diffusion operator in this model. Discarding pheromone diffusion is done for simplicity only and can be easily added. It would add a term ∆ x G or ∆ ω G at the right-hand side of (4.1) according to whether one considers spatial or orientational diffusion.

Ant dynamics. The evolution of the ant distribution function is ruled by the following kinetic equation:

∂ t F + c ω • ∇ x F = Q(F ). (4.2)
The left-hand side describes the ant motion with constant speed c in the direction ω.

The right hand side is a Boltzmann-type operator which describes the rate of change of the distribution function due to the velocity jump processes.

Q is decomposed into Q = Q r + Q p ,
where Q r and Q p respectively describe the random velocity jumps and the trail recruitment jumps. Both operators Q k (x, ω, t), k = p or r express the balance between gain and loss due to velocity jumps, i.e.

Q k (x, ω, t) = Q + k -Q - k .
The gain term Q + k describes the rate of increase of F (x, ω, t) due to particles which have post-jump velocity ω and pre-jump velocity ω ′ . Similarly, the loss term Q - k describes the rate of decay of F (x, ω, t) due to particles jumping from ω to another velocity ω ′ . The jump probability P k (ω → ω ′ )dω ′ is the probability per unit time that a particle with velocity ω jumps to the neighborhood dω ′ of ω ′ due to jump process k. Therefore, the expression of Q k is:

Q k (F )(x, ω, t) = S 1 P k (ω ′ → ω)F (x, ω ′ , t) -P k (ω → ω ′ )F (x, ω, t) dω ′ . (4.3)
where the positive term corresponds to gain and the second term, to loss. By symmetry, we note that

S 1 Q k (F )(x, ω, t) dω = 0, (4.4)
for any distribution F . This expresses that the local number density of particles is preserved by the velocity jump process. Now we describe the expressions of the jump probabilities P k . For both processes, we postulate the existence of a detailed balance principle, which means that the ratio of the direct and inverse collision probabilities are equal to the ratios of the corresponding equilibrium probabilities

P k (ω ′ → ω) P k (ω → ω ′ ) = h k (ω) h k (ω ′ ) , ( 4.5)
where h k is the equilibrium probability of the process k (k = r or k = p). Using (4.5), we can define:

Φ k (ω ′ , ω) = 1 h k (ω) P k (ω ′ → ω) = Φ k (ω, ω ′ ),
which is symmetric by exchange of ω and ω ′ and write

Q k (F )(x, ω, t) = S 1 Φ k (ω, ω ′ ) h k (ω) F (x, ω ′ , t) -h k (ω ′ ) F (x, ω, t) dω ′ . ( 4.6)
From this equation, it is classically deduced that the equilibria, i.e. the solutions of Q k (F ) = 0 are given by F (x, ω, t) = ρ(x, t)h k (ω) with arbitrary ρ. We recall the argument here for the sake of completeness. Indeed, such F are clearly equilibria. Reciprocally, if F is an equilibrium, then, using the symmetry of Φ k leads to

0 = S 1 Q k (F ) F h k dω = - 1 2 (S 1 ) 2 Φ k (ω, ω ′ ) h k (ω) h k (ω ′ ) F (x, ω ′ , t) h k (ω ′ ) - F (x, ω, t) h k (ω) 2 dω dω ′ .
The last expression is the integral of a non-negative function which therefore must be identically zero for any choice of (ω, ω ′ ). It follows that the only equilibria are functions of the form ρh k with ρ only depending on (x, t). It is not clear if the biological processes actually do satisfy the detailed balance property but this hypothesis simplifies the discussion. Indeed, with this assumption, the equilibria h k and the jump probabilities Φ k can be specified independently.

Trail recruitment jumps.

For trail recruitment, we first need to specify the equilibrium distribution as a function of the pheromone distribution. Several options are possible: non-local interactions, local ones, preferential choice, nematic interactions.

1. Non-local interaction. We first introduce the sensing application:

S R (x, ω, t) = 1 πR 2 |x-y|<R G(y, ω, t)dy,
where R represents the perception radius of the particle, i.e. the maximal distance at which it can feel a deposited pheromones. The quantity S R (x, ω, t) represents the density of pheromones pointing towards ω which can be perceived by an ant at point x in its perception area. We also define

T R (x, t) = S 1 S R (x, ω, t) dω,
the pheromone total density within the perception radius, regardless of orientation. Then, we let the equilibrium distribution of the trail recruitment process as follows:

h p (ω) = g R (x, ω, t) := S R (x, ω, t) T R (x, t) , ( 4.7)
which, by construction, is a probability density. Now, The expression for the transition probability reads:

Φ p (ω → ω ′ ; x, t) = λ p γ(T R (x, t))φ p (ω • ω ′ ), (4.8)
where λ p is the trail-recruitment frequency and γ is a dimensionless increasing function of T which accounts for the fact that recruitment by trails increases with pheromone density (in the discrete particle dynamics, we have taken γ(T ) = πR 2 T , the total number of pheromones in the sensing region). The function φ p (ω • ω ′ ) represents the angular dependence of the interaction process and is such that

1 2π S 1 φ p (ω • ω ′ ) dω ′ = 1. (4.9)
We assume that it is independent of the pheromone distribution for simplicity. Inserting (4.8) into (4.6), the trail-recruitment operator is written:

Q p (F )(x, ω, t) = λ p γ(T R (x, t)) S 1 φ p (ω • ω ′ )(g R (ω) F (x, ω ′ , t) -g R (ω ′ ) F (x, ω, t))dω ′ . (4.10)
The choice of φ p which corresponds to the discrete dynamics discussed in the previous sections is φ p (ω • ω ′ ) = 1. Inserting this prescription into (4.10) and using (4.7) leads to the simplified operator

Q p (F )(x, ω, t) = λ p γ(T R (x, t)) ρ(x, t) S R (x, ω, t) T R (x, t) -F (x, ω, t) , with ρ(x, t) = F (x, ω, t)dω,
the local ant density at x.

2.

Local interaction. This corresponds to taking the limit of the sensing radius to zero: R → 0 which leads to

h p (ω) = g(ω) := G(x, ω, t) T (x, t) , T (x, t) = S 1 G(x, ω, t) dω.
T is the local trail density. Then, the expression of the collision operator is easily deduced from (4.10) by changing g R into g. In the case where φ p = 1, we get the expression:

Q p (F )(x, ω, t) = λ p γ(T (x, t)) (ρ(x, t)g(x, ω, t) -F (x, ω, t)) .
3. Preferential choice. We can envision a mechanism by which the ants can sense and choose the most frequently used trails. A possible way to model this preferential choice is by postulating an equilibrium distribution of the form

h p (ω) = g [k] R (ω) = g k R (ω) S 1 g k R (ω) dω , ( 4.11)
with a power k > 1. Indeed, it can be shown [5] that the maxima of g

[k]
R are larger than those of g R and similarly, the minima are lower. Additionally, the monotony is preserved, i.e.

g R (ω) ≤ g R (ω ′ ) =⇒ g [k] R (ω) ≤ g [k] R (ω ′ ), ∀(ω, ω ′ ) ∈ (S 1 ) 2 .
Therefore, taking g

[k]
R (ω) as equilibrium distribution of the ant-pheromone interaction means that the ants choose the trails ω with a higher probability when the trail density in direction ω is high and with lower probability when the trail density is low. The expression of the collision operator is easily deduced from (4.10) by changing g R into g

[k] R . In the case where φ p = 1, we get:

Q p (F )(x, ω, t) = λ p γ(T R (x, t)) ρ(x, t)g [k] R (x, ω, t) -F (x, ω, t) .
This mechanism can also be combined with a local interaction, by replacing g R by the local angular pheromone probability g. We note that this mechanism is not implementable in the discrete dynamics because the operation g → g [k] is only defined for measures g which belong to the Lebesgue space L k (S 1 ). However, sums of Dirac deltas, which correspond to the measure g in the Individual-Based Model, do not belong to this space. Therefore, a smoothing procedure must be applied to such measures beforehand. Since, it is not possible to obtain experimental data about the smoothing procedure and the power k, the preferential choice model has not been used in the numerical experiments of the previous sections.

4. Nematic interaction. The above described ant-pheromone interactions are polar ones, i.e. the pheromones are supposed to have both a direction and an orientation. However, we can easily propose a nematic interaction, for which an ant of velocity ω chooses ω ′ among the pheromone directions and their opposite in such a way that the angle (ω, ω ′ ) is acute, i.e. such that ω • ω ′ > 0. For this purpose, we modify the equilibria of the trail recruitment operator as follows:

h p (x, ω, t) = g (sym) R := S R (x, ω, t) + S R (x, -ω, t) 2T R (x, t) ,
and suppose that

φ p (ω • ω ′ ) = 0, when ω • ω ′ ≤ 0. (4.
12)

The expression of the collision operator is easily deduced from (4.10) by making the change of g R into g

(sym) R and imposing the restriction (4.12). In the case where

φ p (ω • ω ′ ) = 2H(ω • ω ′ ),
where H is the Heaviside function (i.e. the indicator function of the positive real line), we find

Q p (F )(x, ω, t) = λ p γ(T R (x, t)) 1 T R (x, t) ρ + ω (x, t) S R (x, ω, t)+ +ρ - ω (x, t) S R (x, -ω, t) -F (x, ω, t) , with ρ ± ω (x, t) = F (x, ω ′ , t) H(±ω • ω ′ ) dω ′ ,
is the local density of ants pointing in a direction making respectively an acute angle (for ρ + ω ) or obtuse angle (for ρ - ω ) with ω at x. Random velocity jumps. For random velocity jumps, we assume a uniform equilibrium

h r (ω) = 1 2π , into ∂ t G = ν d F -ν e G, (4.19) ∂ t F + c ω • ∇ x F = Q r (F ) + Q p (F ), (4.20) Q p (F ) = λ p γ(T ) [ρ h -F ] , (4.21) Q r (F ) = λ r ρ 2π -F , (4.22) h = g [k] , g = G T , T = S 1 G dω, ρ = S 1
F dω, (4.23) where the meaning of the power [k] operation has been defined at (4.11). We now change to dimensionless variables. We let t 0 , x 0 , ρ 0 , T 0 , be respectively units of time, space, ant density and pheromone density and we introduce x ′ = x/x 0 , t ′ = t/t 0 , ρ ′ = ρ/ρ 0 , T ′ = T /T 0 , F ′ = F/ρ 0 , G ′ = G/T 0 as new variables and unknowns. Specifically, t 0 is chosen to be the macroscopic time scale (e.g. the observation time scale). Similarly, x 0 is the macroscopic length scale (e.g. the size of the experimental arena). We impose x 0 = ct 0 , so that the time and space derivatives in (4.20) are of the same orders of magnitude. This scaling allows us to observe the system at the convection scale where the convection speed of the ant density is finite. We introduce the following dimensionless parameters:

νd = ν d t 0 , νe = ν e t 0 T 0 ρ 0 , λp = λ p t 0 , λr = λ r t 0 .
We make the assumption that the macroscopic time scale t 0 is very large compared to the microscopic time scales λ -1 r and λ -1 p which are both supposed to be of the same orders of magnitude. Indeed, during the time needed for patterns to develop, ants make a large number of jumps of either kind. Following this assumption, we introduce:

ε = 1 λp = 1 λ p t 0 ≪ 1, σ = λr λp = λ r λ p = O(1).
Concerning the pheromone dynamics, we assume that νd and νe are of the same orders of magnitude, which amounts to supposing that pheromone deposition and evaporation balance each other. Indeed, if one of these two antagonist phenomena predominates, then, after some transient the pheromone density will become either too low or too large and we cannot expect any interesting patterns to emerge in this case. We introduce

η = 1 νd = 1 ν d t 0 , κ = νe νd = ν e ν d = O(1).
In what follows, we will assume that η = O(1) i.e. that the pheromone dynamics occurs at the macroscopic time scale. After rescaling, system (4.19)-(4.23) becomes (dropping the primes for the sake of clarity):

η ∂ t G ε = F ε -κ G ε , (4.24) ε (∂ t F ε + ω • ∇ x F ε ) = Q(F ε ), (4.25)
with the collision operator Here, we suppose that η = O(1) i.e. we assume that the pheromone dynamics occurs at the macroscopic scale. We show that the limit ε → 0 of (4.24)-(4.28) consists of the following system for the ant density ρ(x, t), pheromone density T (x, t) and pheromone distribution function g(x, ω, t):

Q = Q r + Q p given by Q(F ) = (Q r + Q p )(F ) = (γ(T ) + σ) (µ ρ -F ), (4.26) µ = γ(T )h + σ 2π γ(T ) + σ , h = g [k] , g = G T , (4.27) T = S 1 G dω, ρ = S 1 F dω. ( 4
∂ t ρ + ∇ x • γ(T ) γ(T ) + σ j h = 0., (4.29)
η∂ t T = ρ -κT, (4.30)

η∂ t g = ρ T γ(T )g [k] + σ 2π γ(T ) + σ -g , ( 4.31)
with h = g [k] and where j ϕ = S 1 ϕ(ω) ω dω, denotes the flux of any function ϕ(ω).

Eq. (4.31) is a closed equation for g. Once g is determined and inserted into (4.29) the evolution of the ant density ρ can be computed. The ant distribution function f is equal to µ at any time, with µ given by (4.27). Indeed, in this limit, supposing that F ε → F , we get Q(F ) = 0 from (4.25). Therefore, from (4.26), we obtain

F = ρµ, or f = µ. (4.32)
The equation for ρ(x, t) is obtained by integrating (4.25) with respect to ω and using (4.4). We find:

∂ t ρ + ∇ x • j F = 0.
Remarking that j F = ρj µ and that the flux of the isotropic distribution vanishes, we finally get from (4.27): Some comments are now in order. In the limit ε → 0 the ant distribution function instantaneously relaxes to the distribution µ. This distribution reflects the antagonist effects of trail recruitment and random velocity jumps. Indeed, µ is the convex combination of the equilibrium distributions h and 1 2π of the two processes respectively. The weights, respectively equal to γ(T )/(γ(T ) + σ) and σ/(γ(T ) + σ) show that the influence of the trail recruitment process is more pronounced at large pheromone densities, since γ increases with T . On the other hand, if the frequency of random jump σ is increased, the trail recruitment process is comparatively less important.

j µ = γ(T ) γ(T ) + σ j h , ( 4
Case k = 1: no preferential choice. If the ants do not implement a preferential choice of the largest trails, i.e. if k = 1, eq. (4.31) simplifies into

η∂ t g = ρ T σ γ(T ) + σ 1 2π -g . (4.35)
This is a classical relaxation equation of g towards the isotropic distribution 1 2π . As a consequence, in this case, there is no trail formation and the large time behavior of the system leads to a homogeneous steady state. This description can be complemented by looking at the pheromone flux j g . Indeed, (4.35) leads to η∂ t j g = -ρ T σ γ(T ) + σ j g .

As a consequence, the direction of the local pheromone flux never changes and its intensity decays to 0 as t → ∞. Additionally, eq. (4.33) which in the case k = 1 gives j µ = γ(T ) γ(T )+σ j g shows that the ant flux is always proportional to and smaller than the pheromone flux. Therefore, it also converges to 0 for large times. Note that this direction may not correspond to the maximum of the pheromone distribution g. Therefore, the ant flux may not be aligned with any particular trail, defined as such a maximum.

The ant distribution µ is just the convex combination of the pheromone distribution g and of the isotropic distribution. Therefore, the ant distribution is always smoother than the pheromone distribution. The random velocity jump process, even if very weak, seems to prevent a positive feedback between the ant and pheromone distributions which could lead to the formation of trails. Of course, these conclusions hold only when ε → 0, i.e. if the equilibrium of the ant jump operator is instantaneously reached. The fact that the simulations do indeed show the formation of trails without any implementation of a preferential choice seems to indicate that the fast microscopic dynamics plays an important role in the formation of trails which the macroscopic model is unable to capture. We also note that if σ = 0, the pheromone distribution is constant in time. This is due to the fact that, in the absence of random velocity jumps, newly created pheromones are deposited according to a distribution which coincides exactly with the current pheromone distribution, resulting in an exact zero balance for this distribution. Therefore, even if σ = 0, no trails can develop.

Case k > 1: existence of a preferential choice. In this case, Eq. (4.31) is a non-local equation due to the operator g → g [k] . No analysis is available yet (to our knowledge) for such an equation (some preliminary results can be found in [5]). The large-time behavior of the system depends on the limit as t → ∞ of eq. (4.31). We note that (4.31) may produce concentrations [5]. Indeed, the contribution of the largest trails is amplified and the ant flux becomes more strongly correlated to the direction of the largest trails. Therefore if the ants choose preferably the largest trail, the resulting concentration dynamics may counterbalance the effect of the random velocity jumps and a positive feedback between the ants and the pheromones is more likely to occur. The study of this case is deferred to future work.

Conclusion on macroscopic models.

We have shown that macroscopic models are unable to develop trail formation without some mechanism allowing to amplify the variations of the pheromone distribution function. We have provided an example of such a mechanism, referred to as the preferential choice and which consists for the ants to choose the strong trails with higher probability than the weak ones. However, the need for such an amplification mechanism is not observed on the simulations of the microscopic model (see section 3). This difference may indicate that the use of such macroscopic models is not fully justified for this dynamics. In particular, the chaos property (see e.g. [START_REF] Cercignani | The mathematical theory of dilute gases[END_REF]), which is the corner stone of the derivation of macroscopic models, may not be valid. Further rigorous mathematical studies are needed to make this point clearer.

Conclusion

In this article, we have introduced an Individual-Based Model of ant-trail formation. The ants are modeled as self-propelled particles which deposit directed pheromones (or pieces of trails) and interact with them through alignment interaction. We have introduced a trail detection technique which provides numerical evidence for the formation of trail patterns, and allowed us to quantify the effects of the biological parameters on the pattern formation. Finally, we have proposed both kinetic and fluid descriptions of this model and analyzed the capabilities of the fluid model to develop trail patterns. From the biological viewpoint, the model can be further improved. The ant and pheromone dynamics can be complexified for instance by adding extra pheromone diffusion, anisotropy or saturation in the pheromone detection mechanism, or by investigating the effect of a non-homogeneous medium. From the mathematical viewpoint, a rigorous derivation of the kinetic and fluid equations are still open problems.

Figure 2 :

 2 Figure2: A typical output of the model at two different times. The ants are represented in blue and the pheromones in green. We clearly observe the formation of trails. Parameters of the simulation: λ p = 2s -1 , ∆t = .05 (see also Table1).

Figure 3 :

 3 Figure 3: In this example, Ants 1, 2 and 3 are linked together: they form a trail. Ant 4 is not linked to Ant 2 since their directions are too different.

  Ø t = 2000 s

Figure 4 :Figure 5 :

 45 Figure 4: The two largest trails (drawn in red and orange) for the simulation depicted in figure 2 at time t = 2000 s. Parameters for the estimation of the trails: r max = 2 cm and θ max = 45 • .

Figure 6 :

 6 Figure 6: Mean trail size S as a function of time for one simulation (dashed line) and averaged over 1000 different simulations (solid line).

Figure 7 :

 7 Figure 7: The mean S of the distribution p t (s) as a function of the ant-pheromone interaction frequency λ p for a fixed value of the random jump frequency λ r .

Figure 8 :Figure 9 :

 89 Figure 8: The mean S of the distribution p t (s) as a function of the pair (λ r , λ p ), The cuts of this surface at a fixed value of λ r shows the same behavior as in figure 7.

Figure 10 : 3 Figure 11 :

 10311 Figure 10: (Left) Construction of the two-particle distribution f 2 . (Right) Histogram of the two-particle density f 2 for the test-case corresponding to the right picture of fig. 2. f 2 (X) is represented via a color scale as a function of the two components of X.

. 28 )

 28 Macroscopic limit ε → 0 of the kinetic model (4.24)-(4.28).

Table 1 :

 1 Since in our model, all the ants lay down Table of the parameters used in the simulations.

	Parameters	Value
	L Box size	100 cm
	N Number of ants	200
	c Ant speed	2 cm/s
	λ r Random jump frequency	2 s -1
	σ Random jump standard deviation .1
	ν d Pheromone deposition rate	.2 s -1
	T p Pheromone lifetime	100 s
	R Detection radius	1 cm
	λ p Trail recruitment frequency	0-3 s -1

Table 2 :

 2 Estimations of the width r 0 (3.2) of the trails using f2 given in figure

		cm)	λ r r 0 (cm)
	3	2.796	0	2.532
	2	3.181	1	2.964
	1	4.350	2	3.181
			3	3.345

  .33) and consequently, ρ satisfies (4.29). To compute the pheromone distribution function g, we integrate (4.24) with respect to ω and get (4.30).

				Then, combining (4.30)
	with (4.24), we deduce that	η∂ t g =	ρ T	(f -g).	(4.34)
	But, with (4.32) and (4.27), we deduce that g satisfies (4.31).

Theraulaz for fruitful discussions. P. Degond wishes to acknowledge the hospitality of the Mathematics Department of Tsinghua University where this research was completed. The work of S. Motsch is partially supported by NSF grants DMS07-07949, DMS10-08397 and FRG07-57227.

with a given jump probability Φ r (ω, ω ′ ) = λ r φ r (ω, ω ′ ).

Here, φ r satisfies the same normalization condition (4.9) as the trail recruitment jump transition probability and λ r is the random velocity jump frequency. With (4.6), we find the expression of Q r :

Summary of the kinetic model. Below, we collect all equations of the kinetic model. We have written the model in the framework of non-local interaction, preferential choice and nematic interaction. The restriction to simpler rules is easily deduced.

S R (x, ω, t) dω. (4.18)

In the following section, we consider fluid limits of the present kinetic model.

Macroscopic model

Scaling. In order to study the macroscopic limit of the kinetic model (4.13)-(4.18), we use the local interaction approximation R = 0, with non-nematic interaction and uniform transition probabilities φ r = 1, φ p = 1. In this case, the model simplifies


P Degond

J-G Liu
email: jliu@phy.duke.edu

S Motsch
email: smotsch@cscamm.umd.edu

Hydrodynamic models of self-organized dynamics: derivation and existence theory

Keywords: Self-propelled particles, alignment dynamics, hydrodynamic limit, diffusion correction, weakly non-local interaction, symmetrizable hyperbolic system, energy method, local well-posedness, capillary force, attraction-repulsion potential AMS Subject classification: 35L60, 35K55, 35Q80, 82C05, 82C22, 82C70, 92D50

 with the addition of an attraction-repulsion interaction potential. Introducing different scalings than in [10], the non-local effects of the alignment and attraction-repulsion interactions can be kept in the hydrodynamic limit and result in extra pressure, viscosity terms and capillary force. The systems are shown to be symmetrizable hyperbolic systems with viscosity terms. A local-in-time existence result is proved in the 2D case for the viscous model and in the 3D case for the inviscid model. The proof relies on the energy method.

Introduction

The context of this paper is the hydrodynamic limit of a kinetic model for self-propelled particles. The self-propulsion speed is supposed to be constant and identical for all the particles. Therefore, the velocity variable reduces to its orientation. The particle interactions consist in two parts: an alignment rule which tends to relax the particle velocity to the local average orientation and an attraction-repulsion rule which makes the particles move closer or farther away from each other. This model is inspired both by the Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] and the Couzin model [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF][START_REF] Couzin | Collective Memory and Spatial Sorting in Animal Groups[END_REF].

The model studied in this paper is a generalization of the model of [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] with the addition of an attraction-repulsion interaction potential. More importantly, a different scaling is investigated. In this scaling, the non-local effects of the alignment and attractionrepulsion interactions are kept in the hydrodynamic limit and result in extra pressure and viscosity terms. Beyond the statement of the model, the main result of the present paper is a local-in-time existence theorem in the 2D case for the viscous model (when the nonlocal effects are retained) and in the 3D case for the inviscid model (when the non-local effects are omitted). Both proofs rely on a suitable symmetrization of the system and on the energy method.

There has been an intense literature about the modeling of interactions between individuals among animal societies such as fish schools, bird flocks, herds of mammalians, etc. We refer e.g. to [START_REF] Aldana | Phase transitions in self-driven many-particle systems and related non-equilibrium models: a network approach[END_REF][START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF][START_REF] Couzin | Collective Memory and Spatial Sorting in Animal Groups[END_REF][START_REF] Grégoire | Onset of collective and cohesive motion[END_REF] but an exhaustive bibliography is out of reach. Among these models, the Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] has received particular attention due to its simplicity and the universality of its qualitative features. This model is a discrete particle model (or 'Individual-Based Model' or 'Agent-Based model') which consists of a timediscretized set of Ordinary Differential Equations for the particle positions and velocities. A time-continuous version of this model and its kinetic formulation are available in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]. A rigorous derivation of this kinetic model from the time-continuous Vicsek model can be found in [START_REF] Bolley | Mean-field limit for the stochastic Vicsek model[END_REF]. In the present paper, we extend this model by adding an attraction-repulsion force.

Hydrodynamic models are attractive over particle ones due to their computational efficiency. For this reason, many such models have been proposed in the literature [START_REF] Carrillo | Double milling in self-propelled swarms from kinetic theory[END_REF][START_REF] Carrillo | Self-propelled interacting particle systems with roosting force[END_REF][START_REF] Chuang | State transitions and the continuum limit for a 2D interacting, self-propelled particle system[END_REF][START_REF] D'orsogna | Self-propelled particles with soft-core interactions: patterns, stability and collapse[END_REF][START_REF] Mogilner | A non-local model for a swarm[END_REF][START_REF] Mogilner | Mutual interactions, potentials, and individual distance in a social aggregation[END_REF][START_REF] Topaz | Swarming patterns in a two-dimensional kinematic model for biological groups[END_REF][START_REF] Topaz | A nonlocal continuum model for biological aggregation[END_REF]. However, most of them are phenomenological. [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] proposes one of the first rigorous derivations of a hydrodynamic version of the Vicsek model (see also [START_REF] Kulinskii | Hydrodynamic model for a system of self-propelling particles with conservative kinematic constraints[END_REF][START_REF] Ratushnaya | Collective behaviour of self propelling particles with kinematic constraints ; the relations between the discrete and the continuous description[END_REF][START_REF] Ratushnaya | Hydrodynamic model for the system of self propelling particles with conservative kinematic constraints; two dimensional stationary solutions[END_REF] for phenomenological derivations). It has been expanded in [START_REF] Degond | A Macroscopic Model for a System of Swarming Agents Using Curvature Control[END_REF] to account for a model of fish behavior where particles interact through curvature control, and in [START_REF] Degond | Diffusion in a continuum model of self-propelled particles with alignment interaction[END_REF] to include diffusive corrections. Other variants have also been investigated. For instance, [START_REF] Frouvelle | A continuous model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF] studies the influence of a vision angle and of the dependency of the alignment frequency upon the local density. [START_REF] Degond | Macroscopic limits and phase transition in a system of self-propelled particles[END_REF][START_REF] Frouvelle | Dynamics, in a kinetic model of oriented particles with phase transition[END_REF] propose a modification of the model which results in phase transitions from disordered to ordered equilibria as the density increases and reaches a threshold, in a way similar to polymer models [START_REF] Doi | The theory of polymer dynamics[END_REF][START_REF] Onsager | The effects of shape on the interaction of colloidal particles[END_REF].

The organization of the paper is as follows. In section 2, we introduce the model of self-propelled particles and set up the associated kinetic equation. We then discuss various scalings which lead to the derivation of the studied hydrodynamic models. We introduce four dimensionless parameters in the problem: the scaled interaction mean-free path ε, the radius of the interaction region η, the noise intensity δ/ε and the relative strength between the attraction-repulsion and the alignment forces s. The scaling considered in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] ignores the attraction-repulsion force and supposes that ε = η → 0, δ = O(1). Here, we investigate four different scaling relations.

1. The weakly non-local interaction scaling without noise:

η = √ ε, δ = 0, s = η 2 .
The resulting model is a viscous hydrodynamic model with constrained velocity on the unit sphere. For this case, we assume that the solutions of the kinetic equation are monokinetic. We justify this assumption by studying the space homogeneous kinetic model and prove that the solutions converge on the fast ε time scale to the monokinetic distribution. We also highlight the variational structure of this space homogeneous kinetic model. Note that the scaling assumption η = √ ε is different from the one used in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]. It corresponds to increasing the size of the interaction region in the microscopic variables by a factor 1/ √ ε, as ε → 0. Therefore, more and more non-local effects are picked up in the hydrodynamic limit. These non-local effects give rise to the viscosity terms in the macroscopic models which make an original addition from previous work.

2. The local interaction scaling with noise. This is the scaling proposed in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] which is recalled here just for the sake of comparisons. It consists in η ≪ ε, δ = 0(1), s ≤ η 2 . The resulting model is the inviscid hydrodynamic model with constrained velocity on the unit sphere.

However, we make the additional assumption that the zero-th order moment of the potential is zero, which expresses some kind of balance between the attraction and repulsion effects. This results in a model like in the previous scaling, but with the addition of a term analog to the capillary force, induced from the attractive part of the potential.

In section 3, we prove local well-posedness for all the models derived in section 2, except the last one (capillary force scaling). All the remaining systems have the same form of a symmetrizable hyperbolic system with additional viscosity. In section 3.1, we prove the local-in-time existence of solutions for the viscous system in 2D and in section 3.2, we show the same result for the inviscid system in 3D based on the energy method. Finally, a conclusion is drawn in section 4.

Derivation of hydrodynamic models 2.1 Individual-Based Model of self-alignment with attractionrepulsion

The starting point of this study is an Individual-Based Model of particles interacting through self-alignment [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] and attraction-repulsion [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF][START_REF] Couzin | Collective Memory and Spatial Sorting in Animal Groups[END_REF]. Specifically, we consider N particles x k ∈ R d moving at a constant speed v k ∈ S d-1 . Each particle adjusts its velocity to align with its neighbors and to get closer or further away. The evolution of each particle is modeled by the following dynamics:

dx k dt = v k (2.1)
dv k = P v ⊥ k v k dt + √ 2d dB k t . (2.2)
Here, P v ⊥ k is the projection matrix onto the normal plane to v k :

P v ⊥ = Id -v ⊗ v.
It ensures that v k stays of norm 1. B k t is a Brownian motion and d represents the noise intensity. Both the alignment and attraction-repulsion rules are encoded in the vector v k :

v k = j k + r k |j k + r k | ,
where j k counts for the alignment and r k for the attraction-repulsion:

j k = j K(|x j -x k |)v j , r k = j Φ ′ (|x j -x k |) x j -x k |x j -x k | . (2.3)
The kernel K is a positive function, Φ ′ can be both negative (repulsion) and positive (attraction). In figure 1, we give an example of functions K and Φ ′ modeling the popular "zone-based" model for fish behavior [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF][START_REF] Couzin | Collective Memory and Spatial Sorting in Animal Groups[END_REF][START_REF] Parrish | Self-organized fish schools: an examination of emergent properties[END_REF]. When the number of particles becomes large (i.e. N → ∞), one can formally derive the equation satisfied by the particle distribution function f (x, v, t) (i.e. the probability distribution of the particles in phase-space (x, v)). Under suitable assumptions [START_REF] Bolley | Mean-field limit for the stochastic Vicsek model[END_REF][START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Sznitman | Topics in propagation of chaos, École d'été de probabilités de Saint-Flour XIX-1989[END_REF], f satisfies:

f t + v • ∇ x f = -∇ v • [(P v ⊥ v f )f ] + d∆ v f, (2.4)
where

v f = j f + r f |j f + r f | , j f = x ′ ,v ′ K(|x ′ -x|)v ′ f (x ′ , v ′ , t) dx ′ dv ′ , r f = -∇ x x ′ ,v ′ Φ(|x ′ -x|) f (x ′ , v ′ , t) dx ′ dv ′ .
The function Φ is the antiderivative of Φ ′ which vanishes at infinity (i.e. Φ(r) r→∞ -→ 0). Using the distribution f , we want to identify the asymptotic behavior of the model in different regimes. This is the purpose of the next section.

Scaling parameters

Introducing two dimensionless parameters ε and η, the starting point is the following scaled version of the previous kinetic model for the distribution function f (x, v, t):

f t + v • ∇ x f = - 1 ε ∇ v • (P v ⊥ v η f )f + d∆ v f, (2.5)
where

v η f = j η f + η 2 r η f |j η f + η 2 r η f | , (2.6)
j η f = (x ′ ,v ′ )∈R n ×S n-1 K |x ′ -x| η v ′ f (x ′ , v ′ , t) dx ′ dv ′ , (2.7
)

r η f = -∇ x (x ′ ,v ′ )∈R n ×S n-1 Φ |x ′ -x| η f (x ′ , v ′ , t) dx ′ dv ′ . (2.8)
The first term (given by j η f ) expresses the alignment interaction (like in the Vicsek dynamics [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF]) while the second term (given by r η f ) expresses the repulsion interaction (like in e.g. [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF]). The expression of v η f means that the alignment term will prevail in the limit η → 0. In this paper, we will consider various possible assumptions concerning the relative speeds of convergence of ε and η to 0.

We denote:

ξ∈R n K(|ξ|) dξ = k 0 , 1 2n ξ∈R n K(|ξ|) |ξ| 2 dξ = k, ξ∈R n Φ(|ξ|) dξ = φ.
We can always assume that k 0 = 1. The potential Φ is said to be repulsive if φ ≥ 0. Defining the moments ρ f and ρ f u f of f by

ρ f = v ′ ∈S n-1 f (v ′ ) dv ′ , ρ f u f = v ′ ∈S n-1 f (v ′ ) v ′ dv ′ ,
we have the following Taylor expansion of v η f :

v η f = Ω f + η 2 1 ρ f |u f | ℓ f + o(η 2 ), Ω f = u f |u f | , ℓ f = P Ω ⊥ f (k∆(ρ f u f ) -φ∇ x ρ f ).
Inserting this expression into the kinetic equation (2.5), we get

f t + v • ∇ x f = - 1 ε ∇ v • [(P v ⊥ Ω f )f ] - η 2 ε 1 ρ f |u f | ∇ v • [(P v ⊥ ℓ f )f ] + d∆ v f + o( η 2 ε ).
(2.9)

We now consider three different scaling limits which lead to models for which we will prove local existence of classical solutions.

Weakly non-local interaction scaling without noise

In this scaling limit, we assume no noise d = 0 and the following ordering between the two parameters ε and η:

ε → 0, η → 0, η 2 ε → 1.
f ε satisfies (keeping only the O(1) terms in ε or larger):

f ε t + v • ∇ x f ε = - 1 ε ∇ v • [(P v ⊥ Ω f ε )f ε ] - 1 ρ f ε |u f ε | ∇ v • [(P v ⊥ ℓ f ε )f ε ] . (2.10)
For f ε to converge, we need to assume that the leading order term at the right-hand side of (2.10) vanishes, i.e. that f ε satisfies:

∇ v • [(P v ⊥ Ω f ε )f ε ] = 0, ∀ε > 0.
This is equivalent to assuming that f ε is a monokinetic distribution, i.e.

f ε (x, v, t) = ρ ε (x, t) δ(v, Ω ε (x, t)), (2.11)
where δ(v, v) is the delta distribution on the sphere at the point v. The assumption of monokinetic distribution requires, in order to be consistent, that there is no noise, which is the reason for assuming d = 0.

Proposition 2.1 For monokinetic solutions (2.11), ρ and Ω are independent of ε and satisfy the following system:

∂ t ρ + ∇ x • (ρΩ) = 0, (2.12) ∂ t (ρΩ) + ∇ x • (ρΩ ⊗ Ω) + φP Ω ⊥ ∇ x ρ = kP u ⊥ ∆(ρu).
(2.13)

Proof. The result follows from multiplying (2.10) by 1 and v and using Green's formula.

Remark 2.1 The repulsive force contributes for a pressure term at the left-hand side of the momentum equation, which otherwise would not be strictly hyperbolic, and would fall in the class of Pressureless Gas Dynamics models [START_REF] Bouchut | On zero pressure gas dynamics[END_REF].

In order to justify the monokinetic assumption, we consider the spatial homogeneous equation:

f t = - 1 ε ∇ v • [(P v ⊥ Ω f )f ] , (2.14)
and show that its solution relaxes to a monokinetic distribution (2.11) at the fast ε time scale. More precisely, we have Proposition 2.2 We assume that f | t=0 dv = 1 so that we have f (v, t) dv = 1 for all times. Therefore, u = v f (v) dv (we omit the index f when the context is clear). We also assume that u t=0 = 0, otherwise, the dynamics is not defined. Then,

f (t) → f (∞) where f (∞) is of the form f (∞) = δ(v, Ω).
Proof. We introduce the variance:

F (f ) = |v -u| 2 f (v) dv = (1 -(u • v)) f (v) dv.
We have

F (f ) = 1 -|u| 2 .
We note that |u| 2 is the classical order parameter [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF]. In Lemma 2.3 below, we prove that F (f ) satisfies the following dissipation equation:

F (f ) t + 2 ε |u| 2 -(v • u) 2 |u| f (v) dv = 0. (2.15) Since |u| 2 -(v • u) 2 ≥ 0, we have F (f ) t ≤ 0. Therefore, F (f ) is a decreasing function of time. Furthermore, if f is a distribution such that |u| 2 -(v • u) 2 |u| f (v) dv = 0,
then, either f = δ(v, Ω) or f is of the form:

f = αδ(v, Ω) + (1 -α)δ(v, -Ω),
where α ∈ [0, 1] and Ω ∈ S n-1 . The second form is called a dipole. So, unless f is a delta or a dipole, F (f ) is strictly decaying. Since a dipole is unstable, it is never reached in the course of the dynamics. Therefore, F (f ) → 0 as t → ∞. Since F (f ) = 0 is equivalent to f = δ(v, Ω), this shows that f → δ(v, Ω) as t → ∞ and the typical convergence time is ε.

Lemma 2.3 Any solution f of (2.14) satisfies (2.15).

Proof. We write

εF (f ) t = (1 -(u • v)) ∂ t f (v) dv -u t • v f (v) dv = (1 -(u • v)) ∂ t f (v) dv -v ∂ t f (v) dv • u = ∂ t ( f (v) dv) -2 (u • v) ∂ t f (v) dv = -2 (u • v) ∂ t f (v) dv.
In the last equality, we have used that f (v, t) dv = 1. Now, multiplying (2.14) by -2(u • v), integrating with respect to v and using Green's formula, we get:

εF (f ) t + 2 (u • P v ⊥ Ω) f (v) dv = 0. But u • P v ⊥ Ω = 1 |u| (|u| 2 -(u • v) 2 ),
which leads to the result.

F (f ) is a free energy for the problem (2.14) and provides a variational structure. First, let us denote by ν = ∇ f F the gradient of F with respect to f . It is defined by

∇ f F , g = δF δf (g)
where g is an increment of f , i.e. a function g(v) satisfying g(v) dv = 0 (so that f + g satisfies the admissibility condition (f + g) dv = 1).

Proposition 2.4 Eq. (2.14) can be recast as

f t - 1 ε|u| ∇ v • {[∇ v (∇ f F )]f } = 0, (2.16)
which shows that the flow of (2.14) has a gradient flow structure in the Wasserstein metrics [START_REF] Villani | Topics in optimal transportation[END_REF]. We have:

F (f ) t + 1 ε|u| |∇ v (∇ f F )| 2 f dv = 0.
(2.17)

Remark 2.2 Eq. (2.17) provides another proof of the decay of F (f ) with time.

Proof. We compute:

δF δf (g) = -2u • g(v) v dv = -2(u • v), g , which yields ν = ∇ f F = -2u • v.
A simple computation shows that

∇ v ν = -P v ⊥ u = -|u|P v ⊥ Ω.
Therefore, eq. (2.14) can be written as (2.16). Now, multiplying by ∇ f F , integrating over v and using Green's formula, we get (2.17), which ends the proof.

Local interaction scaling with noise

In this scaling we assume that d = O( 1 ε ). More precisely, we let:

d = δ ε ,
with δ a given constant. We also assume that ε and η are such that:

ε → 0, η → 0, η 2 ε → 0.
With this last assumption, the O( η 2 ε ) term in (2.9), which results from the non-locality of the average alignment direction, vanishes. Therefore, this scaling keeps only the local contribution of the alignment interaction. The resulting asymptotic problem, keeping only terms of order O(1) or larger, is written:

f ε t + v • ∇ x f ε = - 1 ε {∇ v • [(P v ⊥ Ωε)f ε ] + δ∆ v f ε } . (2.18)
The limit of (2.18) as ε → 0 has been studied in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] in dimension 3 and in [START_REF] Frouvelle | A continuous model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF] in any dimensions. The result is stated in the following theorem.

Theorem 2.5 We have f ε → ρM Ω where M Ω (v) is the Von Mises-Fischer distribution:

M Ω (v) = exp(β(v • Ω)) dv v∈S n-1 exp(β(v • Ω)) dv , β = 1 δ , (2.19)
and ρ and Ω satisfy the following system:

∂ t ρ + c 1 ∇ x • (ρΩ) = 0, (2.20) ρ(∂ t Ω + c 2 Ω • ∇ x Ω) + δP Ω ⊥ ∇ x ρ = 0.
(2.21)

The constants c 1 and c 2 are defined by

c 1 = v∈S n-1 M Ω (v) (v • Ω) dv, c 2 = v∈S n-1 M Ω (v) h(v • Ω) (1 -(v • Ω) 2 ) (v • Ω) dv v∈S n-1 M Ω (v) h(v • Ω) (1 -(v • Ω) 2 ) dv ,
where h(v • Ω) is the Generalized Collision Invariants (GCI) [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] and is defined as follows in the n-dimensional case [START_REF] Frouvelle | A continuous model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF]. Set

ψ a (v) = h(Ω • v) (a • v)
where a ∈ R n is any vector such that a • Ω = 0. Then, ψ a is the unique solution in the Sobolev space H 1 (S n-1 ) with zero mean, of the following elliptic problem:

-∆ v ψ -β(Ω • ∇ v )ψ = a • v.
Proof. We refer to [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] in the three dimensional case and [START_REF] Frouvelle | A continuous model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF] in the general n-dimensional case.

Weakly non-local interaction scaling with noise

In this section, we propose a scaling which unifies the two previous ones. In this scaling we assume that d = δ ε , with δ is a given O(1) constant and that:

ε → 0, η → 0, η 2 ε → 1.
Here η 2 ε → 1 instead of 0 like in the previous section. Inserting these assumptions into (2.9), and keeping terms of order O(1) or larger, we get

f ε t + v • ∇ x f ε + 1 ρ f ε |u f ε | ∇ v • [(P v ⊥ ℓ f ε )f ε ] = - 1 ε {∇ v • [(P v ⊥ Ω f ε )f ε ] + δ∆ v f ε } . (2.22)
The limit ε → 0 is investigated in the following theorem:

Theorem 2.6 We have f ε → ρM Ω where M Ω (v) is the Von Mises-Fischer distribution (2.19
). ρ and Ω satisfy the following system:

∂ t ρ + c 1 ∇ x • (ρΩ) = 0, (2.23) ρ(∂ t Ω + c 2 Ω • ∇ x Ω) + (δ + c 3 φ)P Ω ⊥ ∇ x ρ = c 3 kc 1 P Ω ⊥ ∆(ρΩ), (2.24)
where the constants c 1 and c 2 are defined as in Theorem 2.5 and

c 3 = (n -1)δ + c 2 c 1 .
(2.25)

Remark 2.3 We notice that c k → 1 as δ → 0 for k = 1, 2, 3 and we recover the noiseless system (2.12), (2.13) when δ → 0.

Proof. We write (2.22) as

(T 1 + T 2 )f ε = 1 ε Q(f ε ),
where T 1 + T 2 and Q are respectively the operators appearing at the left and right hand sides of (2.22).

T 1 = ∂ t + v • ∇
x and T 2 is the remaining part of the left-hand side.

Integrating over v and letting f ε → ρM Ω leads to the mass conservation equation (2.20) unchanged, since T 2 is in divergence form and vanishes through integration with respect to v. Now, to get the momentum equation, we proceed like in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]. From the Generalized Collision Invariant property [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], it follows that

P Ω ⊥ v∈S n-1 T (ρM Ω ) h v dv = 0.
Now, the term

P 1 := P Ω ⊥ v∈S n-1 T 1 (ρM Ω ) h v dv,
gives rise to the same expression as in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]. This expression is

P 1 = βαρ∂ t Ω + γΩ • ∇ x Ω + αP Ω ⊥ ∇ x ρ, with α = 1 n -1 v∈S n-1 M Ω (v) h (1 -(v • Ω) 2 ) dv, γ = 1 (n -1)δ v∈S n-1 M Ω (v) h (1 -(v • Ω) 2 ) (v • Ω) dv.
Dividing by αβ, we find the coefficients c 2 and δ of (2.21) (we recall that βδ = 1). We introduce the notation

ℓ := ℓ ρM Ω = P Ω ⊥ (kc 1 ∆(ρΩ) -φ∇ x ρ),
and consider

P 2 := P Ω ⊥ v∈S n-1 T 2 (ρM Ω ) h v dv = 1 c 1 P Ω ⊥ v∈S n-1 ∇ v • [(P v ⊥ ℓ)M Ω ] h v dv.
Using Green's formula, we get

P 2 = - 1 c 1 P Ω ⊥ v∈S n-1 [P v ⊥ ℓ] • ∇ v (h v) M Ω dv.
We note that (P v ⊥ ℓ) • ∇ v φ = ℓ • ∇ v φ, with φ being any component of h v. We deduce that

P 2 = - 1 c 1 P Ω ⊥ v∈S n-1 (ℓ • ∇ v )(h v) M Ω dv = - 1 c 1 P Ω ⊥ v∈S n-1 ∇ v (h v) M Ω dv T ℓ. (2.26)
Now, we use the formulas:

S n-1 ∇ v g dv = (n -1) S 2
vg dv

S n-1 (∇ v g)h dv = (n -1) S 2 vgh dv - S 2 (∇ v h)g dv
for any pair of scalar functions g, h on S n-1 . We recall that ∇ v M Ω = βP v ⊥ ΩM Ω . Since ℓ • Ω = 0, we compute the matrix

D := v∈S n-1 ∇ v (h v) M Ω dv T P Ω ⊥ = (n -1) v∈S n-1 (v ⊗ v)h M Ω dv P Ω ⊥ -β v∈S n-1 (v ⊗ P v ⊥ Ω)h M Ω dv P Ω ⊥ := (n -1)D 1 -D 2 .
We decompose

v = v ⊥ + v , v ⊥ = P Ω ⊥ v, v = (v • Ω)Ω.
Using this decomposition and the fact that integrals of odd degree polynomials of v ⊥ over S n-1 vanish, we have:

D 1 = v∈S n-1 (v ⊥ ⊗ v ⊥ )h M Ω dv P Ω ⊥ = αP Ω ⊥ ,
and

D 2 = β v∈S n-1 ((v + v ⊥ ) ⊗ (Ω -(Ω • v))(v + v ⊥ ))h M Ω dv P Ω ⊥ .
Owing to the fact that any term of the form (A ⊗ v )P Ω ⊥ = 0 for any vector A, we have since v is parallel to Ω:

D 2 = -β v∈S n-1 (v ⊥ ⊗ v ⊥ )(Ω • v) h M Ω dv P Ω ⊥ = -γP Ω ⊥ .
Inserting these results into (2.26), we get

P 2 = - (n -1)α + γ c 1 ℓ.
Collecting all the results and dividing by αβ, we are led to the momentum equation (2.24), which ends the proof.

Attraction-repulsion potential: induced capillary force

In this section, we investigate the case where the attraction-repulsion force term is of the same order as the alignment term in the expression of the alignment direction v η f , i.e. we assume that

v η f = j η f + r η f |j η f + r η f | , (2.27)
|ξ| ξ * ξ 0 Φ(|ξ|)
Figure 2: The attraction-repulsion potential Φ where j η f and r η f are respectively given by (2.7) and (2.8). Note that, by contrast to (2.6), there is no η 2 in front of r η f in (2.27). The Taylor expansion of v η f is now given by

v η f = Ωf + η 2 1 ρ f |u f | ℓ f + o(η 2 ), Ωf = u f -φ∇ x ρ f |u f -φ∇ x ρ f | , ℓ f = P Ω⊥ f (k∆(ρ f u f ) -φ 2 ∇ x ∆ρ f ),
where

1 2n x ′ ∈R n Φ(|ξ|) |ξ| 2 dξ = φ 2 .
Here, we suppose like in [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF], that the potential is repulsive at short scales and attractive at large scales (See Fig. 2). Therefore, Φ(|ξ|) is supposed to decrease for |ξ| ∈ [0, ξ * ] and to increase for |ξ| ∈ [ξ * , +∞). Furthermore, since Φ(|ξ|) is supposed integrable on R n , we have Φ(|ξ|) → 0 as |ξ| → ∞. It results that Φ(ξ * ) < 0 and that Φ ≥ 0 for |ξ| ∈ [0, ξ 0 ] and Φ ≤ 0 for |ξ| ∈ [ξ 0 , +∞) where ξ 0 < ξ * . We make the additional assumption that the zero-th order moment vanishes: φ = 0, which expresses the balance between the attractive and repulsive parts of Φ. Given the above assumptions, the second moment is negative:

φ 2 < 0.
With these assumptions, the Taylor expansion of v η f simplifies and becomes:

v η f = Ω f + η 2 1 ρ f |u f | ℓ f + o(η 2 ), Ω f = u f |u f | , ℓ f = P Ω ⊥ f (k∆(ρ f u f ) -φ 2 ∇ x ∆ρ f ). (2.28)
Now, we can develop the same theory as before, assuming that

ε → 0, η → 0, η 2 ε → 1.
Inserting these assumptions into (2.9), and keeping terms of order O(1) or larger, we get (2.22) but with ℓ given by (2.28). The limit ε → 0 can be performed like in section 2.5 and we obtain the following theorem:

Theorem 2.7 We have f ε → ρM Ω where M Ω (v) is the Von Mises-Fischer distribution (2.19
). ρ and Ω satisfy the following system:

∂ t ρ + c 1 ∇ x • (ρΩ) = 0, (2.29) ρ(∂ t Ω + c 2 Ω • ∇ x Ω) + δP Ω ⊥ ∇ x ρ = c 3 kc 1 P Ω ⊥ ∆(ρΩ) + c 3 |φ 2 |P Ω ⊥ ∇ x ∆ρ, (2.30)
where the constants c 1 and c 2 are defined as in Theorem 2.5 and c 3 as in Theorem 2.6.

Remark 2.4

The last term at the right-hand side of (2.30) has the same expression as the capillary force in fluid dynamics, except for the projection operator P Ω ⊥ . This capillary force is induced from the attractive part of the potential Φ.

3 Existence theory

Existence in 2D with viscosity

This section is concerned with a local existence result in 2D for a system of the general form

∂ t ρ + ∇ x • (ρΩ) = 0, (3.1) ρ(∂ t Ω + cΩ • ∇ x Ω) + P Ω ⊥ ∇ x (p(ρ)) = µP Ω ⊥ ∆(ρΩ), (3.2)
where the constants c ∈ R and µ ≥ 0 are given and the pressure relation p(ρ) satisfies p ′ (ρ) > 0. All systems derived in the previous section can by recast in this form, with a particular choice of c, µ and p(ρ), after time rescaling, except for the last one (section 2.6) involving the capillary force. The system is supplemented with initial data ρ 0 > 0 and Ω 0 such that |Ω 0 | = 1. We assume that the domain is the square box Π 2 = [0, 1] 2 with periodic boundary conditions.

Theorem 3.1 We assume that the initial data belong to H m (Π 2 ) with m > 2. Then, there exists a time T > 0 and a unique solution

(ρ, ϕ) in L ∞ ([0, T ], H m (Π 2 )) ∩ H 1 ([0, T ], H m-1 (Π 2 )
) such that ρ remains positive. If, in addition, µ > 0, then, the solution also belongs to L 2 ([0, T ], H m+1 (Π 2 )).

Proof. In 2D, we can set Ω = (cos ϕ, sin ϕ). We recall that

∂ t Ω = Ω ⊥ ∂ t ϕ, ∇ x • Ω = (Ω ⊥ • ∇ x )ϕ, P Ω ⊥ = Ω ⊥ ⊗ Ω ⊥ ,
with Ω ⊥ = (-sin ϕ, cos ϕ). Then, we have

∆(ρΩ) = ∆ρ Ω + 2Ω ⊥ (∇ x ρ • ∇ x ϕ) -2ρΩ |∇ x ϕ| 2 + ρΩ ⊥ ∆ϕ, Ω ⊥ • ∆(ρΩ) = ρ∆ϕ + 2 (∇ x ρ • ∇ x ϕ).
Therefore, system (3.1), (3.2) is written:

(∂ t + Ω • ∇ x )ρ + ρ (Ω ⊥ • ∇ x )ϕ = 0, (3.3)
(∂ t + cΩ • ∇ x )ϕ + p ′ (ρ) ρ (Ω ⊥ • ∇ x )ρ = µ ∆ϕ + 2 ∇ x ρ • ∇ x ϕ ρ . (3.4)
Introduce ρ = a(ρ) and λ(ρ) such that

a ′ (ρ) = p ′ (ρ) ρ , λ(ρ) = a ′ (ρ)ρ, h(ρ) = 2 lnρ. (3.5)
Then, system (3.3), (3.4) becomes:

(∂ t + Ω • ∇ x )ρ + λ(ρ) (Ω ⊥ • ∇ x )ϕ = 0, (3.6) (∂ t + cΩ • ∇ x )ϕ + λ(ρ) (Ω ⊥ • ∇ x )ρ = µ (∆ϕ + ∇ x h(ρ) • ∇ x ϕ) .
(3.7)

From (3.6), we have the following a priori estimate (maximum principle):

ρ min exp(- t 0 ∇ x ϕ(•, s) L ∞ (Π 2 ) ds) ≤ ρ ≤ ρ max exp( t 0 ∇ x ϕ(•, s) L ∞ (Π 2 ) ds), (3.8)
where ρ min = min x∈Π 2 ρ 0 (x), ρ max = max x∈Π 2 ρ 0 (x). We remind the following lemmas [START_REF] Taylor | Partial Differential Equations III[END_REF]:

Lemma 3.2 For any pair of functions f , g in H m (R n ) ∩ L ∞ (R n ), we have: f g H m ≤ C ( f H m g L ∞ + f L ∞ g H m ) . If additionally, we suppose that ∇f ∈ L ∞ (R n ), we have, for any α ∈ N n , with |α| = n i=1 α i = m : D α (f g) -f D α g H m ≤ C ( f H m g L ∞ + ∇ x f L ∞ g H m-1 ) , where D α = ∂ x α 1 1 ...x αn n .
Now, with |α| ≤ m, we take the D α derivative of (3.6) and multiply it by D α ρ and integrate it with respect to x. Similarly, we take the D α derivative of (3.7) and multiply it by D α ϕ and integrate it with respect to x. We sum up the resulting identities. Using the notation f, g = Π 2

f g dx, we find:

0 = D α ρ, D α ρt + D α ϕ, D α ϕ t + D α ρ, D α ((Ω • ∇ x )ρ) + c D α ϕ, D α ((Ω • ∇ x )ϕ) + D α ρ, D α (λ(ρ) (Ω ⊥ • ∇ x )ϕ) + D α ϕ, D α (λ(ρ) (Ω ⊥ • ∇ x )ρ) -µ D α ϕ, D α ∆ϕ -µ D α ϕ, D α (∇ x h(ρ) • ∇ x ϕ) = I 1 + . . . + I 5
Then:

I 1 = 1 2 d dt ( D α ρ 2 + D α ϕ 2 ),
and

I 4 = µ D α ∇ϕ 2 ,
where • just indicates an L 2 norm. Now, for the remaining terms, we have the following lemma Lemma 3.3 We have:

|I k | ≤ C( ρ W 1,∞ + ϕ W 1,∞ ) ( ρ 2 H m + ϕ 2 H m ), k = 2, 3, |I 5 | ≤ µ 2 ∇D α ϕ 2 + C( ρ 2 W 1,∞ + ϕ 2 W 1,∞ ) ( ρ 2 H m + ϕ 2 H m ),
where C denote generic constants depending on the parameters of the problem.

The proof of the lemma is postponed at the end. Adding all these terms together for all possible indices α such that |α| ≤ m, we have,

1 2 d dt ( ρ 2 H m + ϕ 2 H m ) + µ ∇ϕ 2 H m ≤ µ 2 ∇ϕ 2 H m + +C( ρ 2 W 1,∞ + ϕ 2 W 1,∞ + 1) ( ρ 2 H m + ϕ 2 H m ). For m ≥ n 2 + 1, we have ρ W 1,∞ + ϕ W 1,∞ ≤ C( ρ H m + ϕ H m ),
and get 1 2

d dt ( ρ 2 H m + ϕ 2 H m ) + µ 2 ∇ϕ 2 H m ≤ C ( ρ 2 H m + ϕ 2 H m + 1) 2 .
Gronwall's inequality leads to the local existence of a solution (ρ, ϕ) in L ∞ ([0, T ], H m (Π 2 )) which, if µ > 0, also belongs to L 2 ([0, T ], H m+1 (Π 2 )) and which satisfies the a priori bound (3.8). To get time regularity, we directly use eqs. (3.6), (3.7), take the H m-1 norm, apply Lemma 3.2, and find

ρt H m-1 + ϕ t H m-1 ≤ C ϕ H m+1 + C( ρ W 1,∞ + ϕ W 1,∞ ) ( ρ H m + ϕ H m ).
Using the previous estimates, we deduce that (ρ, ϕ) also belongs to

H 1 ([0, T ], H m-1 (Π 2 )).
The estimates on ρ immediately transfer to ρ since a(ρ) is smooth and invertible for ρ > 0.

Proof of Lemma 3.3. Estimate of I 5 : Using Green's formula and Cauchy-Schwartz inequality, we have:

|I 5 | ≤ µ ∇D α ϕ ∇ x h(ρ) • ∇ x ϕ H m-1 ≤ C ∇D α ϕ ( ρ H m ∇ϕ L ∞ + ∇ρ L ∞ ϕ H m ) ≤ µ 2 ∇D α ϕ 2 + C ( ρ H m ∇ϕ L ∞ + ∇ρ L ∞ ϕ H m ) 2 ≤ µ 2 ∇D α ϕ 2 + C( ρ 2 W 1,∞ + ϕ 2 W 1,∞ ) ( ρ 2 H m + ϕ 2 H m ).
The second inequality uses Lemma 3.2 and the third one uses Young's inequality.

Estimate of I 3 : We write

I 3 = D α ρ, λ(ρ) (Ω ⊥ • ∇ x )D α ϕ + D α ϕ, (λ(ρ) (Ω ⊥ • ∇ x )D α ρ) + D α ρ, D α (λ(ρ) (Ω ⊥ • ∇ x )ϕ) -λ(ρ) (Ω ⊥ • ∇ x )D α ϕ + D α ϕ, D α (λ(ρ) (Ω ⊥ • ∇ x )ρ) -λ(ρ) (Ω ⊥ • ∇ x )D α ρ = J 1 + J 2 + J 3 .
Using Green's formula, we find

|J 1 | = | ∇ • (λ(ρ) Ω ⊥ ) D α ρ, D α ϕ | ≤ C( ρ W 1,∞ + ϕ W 1,∞ ) ( ρ 2 H m + ϕ 2 H m ).
Now, using Cauchy-Schwartz inequality and applying Lemma 3.2, we find that J 2 and J 3 satisfy the same inequality.

Estimate of I 2 : The proof is similar as for I 3 and is omitted. If sin θ > 0, then this system is a symmetrizable hyperbolic system. We can apply proposition 2.1 p. 425 of [START_REF] Taylor | Partial Differential Equations III[END_REF] and the following theorem follows immediately: Theorem 3. [START_REF] Bouchut | On zero pressure gas dynamics[END_REF] We assume that the initial data (ρ 0 , θ 0 , ϕ 0 ) belong to H m (Π 3 ) with m > 5/2 with ρ 0 > 0, sin θ 0 > 0. Then, there exists a time T > 0 and a unique solution (ρ, θ, ϕ) in L ∞ ([0, T ], H m (Π 3 )) ∩ H 1 ([0, T ], H m-1 (Π 3 )) such that ρ remains positive.

Conclusion

In this paper, we have derived hydrodynamic systems from kinetic models of self-propelled particles with alignment interaction and attraction-repulsion force. We have particularly focused on the inclusion of diffusion terms under the assumption of weakly non-local interactions. Then, we have proved the local-in-time existence of solutions for the viscous system in 2D and a similar result for the inviscid system in 3D. The methods rely on a suitable symmetrization and on the energy method. Future works in this direction will consist in continuing the exploration of the mathematical structure of the system and particularly, trying to prove local existence of the viscous system in 3D and the treatment of the geometric singularity near sin θ = 0. Another direction of work will consist of the numerical quantification of the viscosity as a consequence of the non-locality of the interaction.

Figure 1 :

 1 Figure 1: An example of functions K and Φ ′ modeling the "zone-based" model: repulsionalignment-attraction.

in

  Cartesian coordinates x = (x, y, z), where A k (U), k = 0, . . . , 4 are all symmetric matrices and

The weakly non-local interaction scaling with noise. This scaling unifies the two previous scalings. It consists in η = √ ε, δ = 0(1), s = η 2 . Again, the resulting model is a viscous hydrodynamic model with constrained velocity on the unit sphere, but with modified coefficients as compared to the first scaling. We note however, that in the zero noise limit δ → 0, we recover the system obtained with the first scaling, which provides another justification of the monokinetic assumption in the derivation of the model.

Capillary force scaling. This corresponds to η = √ ε, δ = 0(1), s = 1. Therefore, here, the attraction repulsion force is of the same order as the alignment force.

The authors wish to acknowledge the hospitality of Mathematical Sciences Center and Mathematics Department of Tsinghua University where this research was performed. The research of J.-G. L. was partially supported by NSF grant DMS 10-11738.

Existence in 3D without viscosity

In this section, we investigate the local existence for the inviscid problem in 3 dimensions:

where the parameters and data have the same meaning as in section 3.1. We consider the system in the domain Π 3 = [0, 1] 3 with periodic boundary conditions. For this purpose, we use the spherical coordinates associated to a fixed Cartesian basis. In this basis, denoting by θ ∈ [0, π] the latitude and ϕ ∈ [0, 2π] the longitude, we have Ω = (sin θ cos ϕ, sin θ sin ϕ, cos θ) T , and we let Ω θ and Ω ϕ be the derivatives of Ω with respect to θ and ϕ. We note that

We will use the formulas

where a is an arbitrary vector. Introduce ρ and λ(ρ) as in (3.5). Then, system (3.9), (3.10) becomes:

this system is written