Thierry Goudon
email: thierry.goudon@inria.fr

Boniface Nkonga
email: boniface.nkonga@unice.fr

Michel Rascle
email: michel.rascle@unice.fr

Magali Ribot
email: magali.ribot@unice.fr

Self-organized populations interacting under pursuit-evasion dynamics

Keywords: Collective behavior. Self-propelling particles. Self-organization. Kinetic models. Hydrodynamic models. Math. Subject Classification. 92D25, 92C17, 74A25, 76N10

We discuss the modelling of interacting populations through pursuit-evasion --or attraction-repulsion --principles : preys try to escape chasers, chasers are attracted by the presence of preys. We construct a hierarchy of models, ranging from ODEs systems with finite numbers of individuals of each population, to hydrodynamic systems. First order macroscopic models look like generalized "two-species Keller-Segel equations". But, due to cross-interactions, we can show that the system does not exhibit any blow up phenomena in finite time. We also obtain second order models, that have the form of systems of balance laws, derived from kinetic models. We bring out a few remarkable features of the models based either on mathematical analysis or numerical simulations.

Introduction

The formation of space-time heterogeneous patterns is a universal feature of living organisms. Many attempts can be found to model the behavioral mechanisms that lead to the observed self-organization of interacting populations. We refer the reader to [START_REF] Parrish | Self-organized fish schools: an examination of emergent properties[END_REF] for a detailed introduction to such phenomena, in the specific case of fish schools, and to the review in [START_REF] Vicsek | Collective motion[END_REF]. The question can be rephrased as to describe the information exchanges between self-propelled individuals: the individuals can use the information contained in a certain subdomain of their environment, so that, according to a set of basic rules, the motion of the whole population organizes with a remarkable pattern. This phenomenon is referred to as "flocking" or "swarming". The mathematical modelling of these natural behaviors has motivated an intense research activity, after the seminal works of Vicsek et al. [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF]. Flocking can be represented by hydrodynamical models [START_REF] Toner | Flocks, herds, and schools: A quantitative theory of flocking[END_REF] as well as many particles systems, where the interaction between individuals is embodied into some potential [START_REF] Cucker | Emergent behavior in flocks[END_REF][START_REF] Cucker | On the mathematics of emergence[END_REF]. This potential encodes how the motion of an individual adapts to the others, e. g. by adjusting the relative velocities. The advantage of hydrodynamical models is to describe the dynamic through a reduced set of macroscopic quantities, like the concentration and bulk velocity. Individual-based and continuum models, stochastic or not, have led to original problems for mathematical analysis and fascinating numerical simulations that reproduce certain features of natural phenomena [START_REF] Chuang | State transition and the continuum limit for the 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]. Hierarchies of models, à la BBGKY, have been derived, which offer a complete picture, ranging from N -particles systems to hydrodynamic models, through kinetic equations where individuals are described according to the principles of statistical physics, see [START_REF] Carrillo | Particle, kinetic, and hydrodynamic models of swarming[END_REF][START_REF] Degond | Macroscopic limits and phase transition in a system of self-propelled particles[END_REF][START_REF] Ha | From particle to kinetic and hydrodynamic descriptions of flocking[END_REF][START_REF] Vecil | A numerical study of attraction/repulsion collective behavior models: 3D particle analyses and 1D kinetic simulations[END_REF]. Roughly speaking, in these models, the motion of the individuals is driven by the combination of self-propulsion, friction and an attractive/repulsive potential. The latter has the general shape of the Morse potential, describing the tendency to pack individuals together, up to a certain critical distance where a repelling effect dominates. Interestingly, these models lead to a large variety of possible behaviors, with complex selection mechanisms driven by certain thresholds on the parameters of the equations.

Here, we are considering a different situation since we deal with two populations and we address the question of constructing a mathematical model that produces such self-organized patterns through pursuit-evasion --or attraction-repulsion --simple principles. We shall discuss "toy-models", certainly (and intentionally) too rough to capture quantitative features; nevertheless, the interactions of the idealized populations we are dealing with are sufficient to bring out relevant behaviors. We neglect direct interactions between individuals of the same specie, the motion is only determined by the potential created by the other population. The dynamic can be seen as the interaction between preys and chasers, described by their respective concentrations: preys are repelled by the chasers while chasers are attracted by the presence of preys. The model can equally be interpreted as a simple "cops and robbers game" (we warn the reader not to confuse with a more complex problem referred to with this name in computer science and graph theory). We refer the reader to [START_REF] Fasano | Equilibrium of two populations subject of chemotaxis[END_REF][START_REF] Tyutyunov | A minimal model of pursuit-evasion in a predator-prey system[END_REF] for similar attempts, but with a different definition of the attraction-repulsion mechanisms. The modelling also echoes to certain applications in robotics. In this work, the potential is obtained through convolution formulae with the densities. A possible model for designing the interaction potential can be obtained by mimicking the repulsive/attractive effects used when dealing with charged or gravitational particles respectively. The latter principle has been adapted from astrophysics, see [START_REF] Chandrasekhar | Brownian motion, dynamical friction and stellar dynamics[END_REF], to population dynamics and it leads to remarkable aggregation phenomena, typical of the behavior of certain populations of bacteria, see [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF]. At first sight, (one of) our hydrodynamic models share(s) the structure of the Keller-Segel system. This system has motivated a huge amount of mathematical works because it exhibits interesting singularity formations: we refer the reader for instance to [START_REF] Greenberg | Stability results for a diffusion equation with functional drift approximating a chemotaxis model[END_REF][START_REF] Herrero | Singularity patterns in a chemotaxis model[END_REF][START_REF] Jäger | On explosions of solutions to a system of partial differential equations modelling chemotaxis[END_REF] for the analysis of such phenomena. In the classical Keller-Segel model the individuals (cells, bacteria) move according to the gradient of the concentration of a substance they emit themselves: the higher the concentration of individuals, the higher the production of the attracting chemical signal. This principle is the basis of chemotaxis. This aggregation process is counter-balanced by spacial diffusion: the two effects compete to determine whether or not the solution blows up in finite time. Thresholds on the initial mass can be discussed accordingly. At least in dimension 2, the situation is quite well understood now. The Keller-Segel system can be obtained through hydrodynamic regimes from a kinetic model for chemotaxis, based on run and tumbling responses to the chemoattractant: this approach is proposed and analyzed in [START_REF] Bournaveas | Global existence for a kinetic model of chemotaxis via dispersion and Strichartz estimates[END_REF] and [START_REF] Perthame | Transport Equations in Biology[END_REF]Section. 5.7]. Several modifications of the model have been introduced in order to prevent the overcrowding: concentration-dependent chemotactic sensitivity and diffusion coefficient, reaction terms, cross-diffusion effects, etc. [START_REF] Bedrossian | Local and global well-posedness for aggregation equations and Patlak-Keller-Segel models with degenerate diffusion[END_REF][START_REF] Carrillo | Cross diffusion and nonlinear diffusion preventing blow up in the Keller-Segel model[END_REF] We refer the reader to the overviews on chemotaxis models in [START_REF] Hillen | A user's guide to PDE models for chemotaxis[END_REF][START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences[END_REF][START_REF] Horstmann | Generalizing the Keller-Segel model: Lyapunov functionals, steady state analysis, and blow-up results for multi-species chemotaxis models in the presence of attraction and repulsion between competitive interacting species[END_REF] and [START_REF] Perthame | Transport Equations in Biology[END_REF]Chap. 5] for further details and results. It is also worth mentioning that similar ideas are also at the basis of PDEs systems proposed to model criminal behavior [START_REF] Short | A statistical model of criminal behavior[END_REF]: these models are intended to reproduce the formation of "hotspots" of criminal activity.

The paper is organized as follows. We start by introducing first order models: the concentrations of preys and chasers obey transport equations, the velocities of which are gradients of potentials satisfying Poisson equations. The right hand side of the Poisson equation is proportional to the concentration of the opposite population, the sign depending whether the effect is attractive or repulsive. Coming back to a single specie, we obtain a diffusionless Keller-Segel equation, as analyzed in [START_REF] Nieto | High-field limit of the Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF], but we shall see that in the present context the crossed effects between preys and chasers prevent the formation of blow up. We discuss N -particles versions of the model. Within this interpretation, it turns out that it might be relevant to replace the Poisson kernel by convolution kernels that take into account further distance effects. Next, we turn to second order models where the presence of preys and chasers is interpreted as creating attractive and repulsive forces. Hence, we obtain individualbased models that have the form of non-linear ODEs systems derived from the standard principles of classical mechanics. We propose a kinetic version of such models. Finally, based on asymptotic arguments we set up a hierarchy of hydrodynamic-like systems for the interacting populations, see e. g. [START_REF] Carrillo | Particle, kinetic, and hydrodynamic models of swarming[END_REF][START_REF] Degond | Macroscopic limits and phase transition in a system of self-propelled particles[END_REF][START_REF] Ha | From particle to kinetic and hydrodynamic descriptions of flocking[END_REF][START_REF] Goudon | Hydrodynamic limit for the Vlasov-Poisson-Fokker-Planck system: analysis of the two-dimensional case[END_REF][START_REF] Nieto | High-field limit of the Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Poupaud | Parabolic limit and stability of the Vlasov-Poisson-Fokker-Planck system[END_REF] for similar arguments in different contexts. Section 3 is concerned with numerical simulations. We investigate the behavior of individual-based and continuum models in 1D and 2D, discussing the role of the parameters entering into the models. The simulations exhibit an interesting variety of behaviors, that could be valuable sources for further mathematical analysis. Eventually, Section 4 is devoted to the analysis of the first order continuum model. In contrast with the single specie problem analyzed in [START_REF] Nieto | High-field limit of the Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Goudon | Multidimensional high-field limit of the electrostatic Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF], concentration in finite time cannot occur as far as initial data are bounded functions. Therefore, we establish the existence and uniqueness of bounded weak solutions, for a wide variety of interaction kernels.

A hierarchy of models for pursuit-evasion dynamics

In what follows we will discuss models based either on ODEs or PDEs. The two viewpoints are intimately connected, as a consequence of the following basic remark about the transport equation

∂ t ρ + ∇ x • (ρu) = 0. ( 1
)
We start by assuming that u : R × R N → R N is smooth enough, so that characteristics curves are well defined by the ODE

d dt X(t; s, x) = u(t, X(t; s, x)), X(s; s, x) = x.
Namely X(t; s, x) is the position at time t of a particle driven by the velocity field u, knowing that it starts from position x at time s. Our discussion will use the following claim (for the sake of self-containedness the proof is detailed in Appendix A, see also [START_REF] Bouchut | Introduction to the mathematical theory of kinetic equations[END_REF]).

Proposition 2.1

i) The measure I i=1 δ(x = X(t; 0, x 0,i )) is a solution of (1) associated to the initial data I i=1 δ(x = x 0,i ). ii) For ρ init ∈ L p (R N ), 1 ≤ p ≤ ∞, the unique solution of (1) having ρ init as initial data is given by

ρ(t, x) = ρ init (X(0; t, x)) exp - t 0 (∇ x • u)(s, X(s; t, x)) ds . The solution lies in C 0 ([0, T ]; L p (R N )) when p is finite, or in C 0 ([0, T ]; L ∞ (R N )- weak-⋆) otherwise.
We shall deal with non linear models where, roughly speaking, the velocity u depends on the concentration ρ through non-local definitions. It leads to mathematical difficulties since the regularity necessary to define properly the characteristic curves is not directly guaranteed in this context (see for instance [START_REF] Goudon | Multidimensional high-field limit of the electrostatic Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Nieto | High-field limit of the Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF] for the analysis of a similar problem). Nevertheless it is worth bearing in mind Proposition 2.1 to make connection, at least formally, between individual-based modelling and PDEs description.

First order models: coupling determines velocities

We consider two interacting populations: the preys (robbers) and the chasers (cops), described by their concentrations ρ p (t, x) and ρ c (t, x), respectively. We start with a simple model where the dynamics is driven by the interactions between the two species, according to the following basic rule: the presence of the other population directly influences the velocity of the individuals, the chasers being attracted by the preys, the preys trying to escape to the chasers. The concentrations obey the following transport equations

∂ t ρ p + ∇ x • (ρ p V p ) = 0, ∂ t ρ c + ∇ x • (ρ c V c ) = 0, (2)
where we assume that the velocities V p , V c derive from the corresponding potentials Φ p , Φ c :

V p = ∇ x Φ p , V c = ∇ x Φ c . ( 3
)
Directly inspired from chemotaxis theory, see e. g. [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF], we define them through the following Poisson equations

∆ x Φ p = αρ c (with α > 0), -∆ x Φ c = ρ p . ( 4
)
Actually, this definition should be thought of as a convolution relation

Φ p (t, x) = -αE ⋆ ρ c (t, x), Φ c (t, x) = E ⋆ ρ p (t, x)
where E stands for the elementary solution of (-∆), mind the minus sign, on R N , whose definition depends on the space dimension:

E(x) =                    - |x| 2 if N = 1, 1 2π ln(|x|) if N = 2, C N |x| N -2 if N > 2 with C N > 0.

Remark 2.2 (Bounded domains, boundary conditions)

The definition (4) should be adapted when the motion is considered in a bounded domain Ω ⊂ R N . In such a situation it is natural to assume that the flux vanishes at the boundary: denoting by ν(x) the unit outward normal vector on x ∈ ∂Ω, we have

V p • ν = 0, V c • ν = 0,

which amounts to impose Neumann boundary conditions on the potentials

∂ ν Φ p = 0, ∂ ν Φ c = 0.
The Poisson equation should be adapted to this case. We assume instead that individuals are sensitive to fluctuations about the mean value of the "opposite" concentration, namely (4) is replaced by

∆ x Φ p (t, x) = α ρ c (t, x) - 1 |Ω| Ω ρ c (t, y) dy with α > 0, -∆ x Φ c (t, x) = ρ p (t, x) - 1 |Ω| Ω ρ p (t, y) dy .
(5)

Equations ( 5) have equally to be used dealing with periodic domains.

The coupling between the two species makes the dynamic completely different from the Keller-Segel system. In particular, we are able to justify the existence-uniqueness of globally defined weak solutions, as well as the preservation of regularity.

Theorem 2.3 Let ρ p,Init , ρ c,Init be non negative functions in L 1 ∩ L ∞ (R N ). Let 0 < T < ∞. Then, there exists a unique weak solution (ρ p , ρ c ) ∈ L ∞ (0, T ; L 1 ∩ L ∞ (R N )) of (2)-(4) with initial data ρ p,Init , ρ c,Init . The solution lies in C 0 ([0, T ]; L q (R N )) for any 1 ≤ q < ∞. If, moreover, ∇ x ρ p,Init , ∇ x ρ c,Init belong to L ∞ (R N ), then the solution lies in L ∞ (0, T ; W 1,∞ (R N )).
This statement will be discussed in Section 4. Of course with velocities defined by (3)-( 4), it is far from clear how the result in Proposition 2.1 can be applied: regularity of the velocity fields can be lacking to define correctly the characteristics curves. We shall go back to the technical issues later on. However, the formula can be helpful to guide the intuition and we are going to use it formally. For instance, it is instructive to study the mere one-dimensional framework where x ∈ R. We get

∂ x Φ p (t, x) = α 2 sgn(x -y)ρ c (t, y) dy = α 2 x -∞ ρ c (t, y) dy - α 2 +∞ x ρ c (t, y) dy, ∂ x Φ c (t, x) = - 1 2 sgn(x -y)ρ p (t, y) dy = - 1 2 x -∞ ρ p (t, y) dy + 1 2 +∞ x ρ p (t, y) dy.
It corresponds to the physical intuition in the discrete 1D case: consider a chaser and a prey. If the chaser is on his left side, the prey tends to escape moving to the right, and so does the chaser. In this situation, with p 0 , c 0 ∈ R the initial positions of the prey and the chaser, their trajectories are defined by

X p (t) = α 2 sgn(p 0 -c 0 )t + p 0 1 0≤t≤T 0 + Q 0 1 t>T 0 , X c (t) = 1 2 sgn(p 0 -c 0 )t + c 0 1 0≤t≤T 0 + Q 0 1 t>T 0 , T 0 =          +∞ if α ≥ 1, 2|p 0 -c 0 | 1 -α if α < 1, Q 0 = X p (T 0 ) = X c (T 0 ) = p 0 -c 0 1 -α + c 0 when α < 1.
Figure 1 shows a numerical simulation of this specific case. Up to the regularity issue, by Proposition 2.1, we can check that the pair (ρ p (t,

x) = δ(x = X p (t)), ρ c (t, x) = δ(x = X c (t)
) is a solution of (2)-(4) (see Remark 2.4 below). On Figure 2, we exhibit the same numerical test in a periodic setting. More precisely, the unknowns X p and X c are constrained to remain in the domain [0, L] (prey and chaser can be thought of as running on a ring of length L). The kernel which appears in the definition of the velocity is adapted:

min k∈Z sgn(X p -(X c -kL)) replaces sgn(X p -X c ).
When V p < V c (corresponding to α < 1), the behavior is the same as in the free-space situation: the chaser reaches the prey in finite time. When V p > V c , the velocity of the prey is larger than the one of the chaser, an equilibrium situation establishes where the two individuals remain at a constant distance.

Remark 2. [START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF] It is clear that Proposition 2.1 does not apply when the concentrations are (sum of) Dirac masses. This case is not covered by Theorem 2.3 neither. In this situation the velocity defined by the convolution formula is a bounded measurable function and the product ρV is not well defined. The difficulty consists in defining the product when jumps of V coincide with atoms of the measure ρ, see [START_REF] Poupaud | Measure solutions to the linear multi-dimensional transport equation with non smooth coefficients[END_REF] and the references therein for related issues. This is precisely the case when the two individuals meet (see on Figure 1, with the encouter time 0.3 ≤ T 0 ≤ 0.4). The observed numerical solution where the velocity of the two individuals vanishes after the encounter might look "natural" on a physical viewpoint, but the mathematical basis to decide how to select this specific solution might be quite subtle. In particular, it is interesting to compare the result with exactly two individuals, one prey and one chaser, in Figure 1, with the case in Figure 3 where we consider a set of several preys (resp. chasers) located randomly around a given position according to a peaked normal law. At the meeting point, we observe that the population of preys splits into two parts, one block escaping to the chasers, the other being pursued by the whole population of chasers, which adopts a reduced speed. Therefore, the stability of measure valued solutions of (2)-( 4) is certainly a complicated issue. We shall go back to these issues elsewhere. on the relative position of the two individuals, which is natural, but the modulus is constant and does not depend on the distance between the individuals. It might be relevant to use instead a velocity profile as in Figure 4. Such a profile depends on 3 parameters:

• A maximal speed: the largest speed that individuals may reach when imposed by the circumstances,

• A security radius: when individuals of the other category are at a distance smaller than this threshold, the motion holds at the maximal speed,

• A vigilance radius: the speed of pursuit/escape decays smoothly as the distance to other individuals becomes larger and there is no reason to move when they are far enough.

Of course, the value of the parameters might differ for the two populations and it will be interesting to investigate how their variations influence the behavior of the populations.

In what follows, we denote by Γ a function with such a profile. For instance we can set Adopting such a modelling (and still reasoning in 1D for the time being) we consider a set of N p preys and N c chasers. They are described by their positions

x ∈ R → Γ(x) = V × 1 |x|≤R + e -(|x|-R) 2 /ǫ 1 |x|>R ∈ [0, V ]. (6)
X p,1 , ...X p,Np , X c,1 , ...X c,Nc ∈ R. The dynamics is embodied into the ODE system d dt X p,i = Nc j=1 sgn(X p,i -X c,j ) Γ p (X p,i -X c,j ), d dt X c,i = Np j=1 sgn(X p,j -X c,i ) Γ c (X p,j -X c,i ). (7)
Let us set [START_REF] Bouchut | Introduction to the mathematical theory of kinetic equations[END_REF] with velocities defined by this formula instead of ( 3)-( 4) (that would correspond to Γ p = α/2, Γ c = 1/2). Of course it makes sense to consider (2), coupled to [START_REF] Chuang | State transition and the continuum limit for the 2D interacting, self-propelled particle system[END_REF], for initial data which are not necessarily sums of Dirac masses.

V p (t, x) = sgn(x -y)Γ p (x -y) ρ c (t, y) dy, V c (t, x) = -sgn(x -y)Γ c (x -y) ρ p (t, y) dy. (8) By Proposition 2.1, ρ p (t, x) = Np i=1 δ(x = X p,i (t)), ρ c (t, x) = Nc i=1 δ(x = X c,i (t)) satisfies the PDE
It is straightforward to extend this approach to any dimension. With positions in R N , we consider scalar functions G p and G c : (0, ∞) → R and we define the potentials through the convolution with the radially symmetric functions

x ∈ R N → G p,c (|x|) ∈ R,
where |•| stands for the usual Euclidean norm on R N . For instance G p,c can be obtained from a primitive of a function like [START_REF] Carrillo | Cross diffusion and nonlinear diffusion preventing blow up in the Keller-Segel model[END_REF]. In the multidimensional framework, the ODE system becomes

d dt X p,i = Nc j=1 X p,i -X c,j |X p,i -X c,j | G ′ p (|X p,i -X c,j |), d dt X c,i = Np j=1 X p,j -X c,i |X p,j -X c,i | G ′ c (|X p,j -X c,i |). (9)
Likewise, for the PDE, the velocities in ( 2) are given by

V p (t, x) = x -y |x -y| G ′ p (|x -y|)ρ c (t, y) dy, V c (t, x) = - x -y |x -y| G ′ c (|x -y|)ρ p (t, y) dy. ( 10
)
It generalizes the case (3)-( 4) where velocities are obtained by convolution with the elementary solution of -∆. Note however the following fundamental differences between (4) and the convolution with x → x |x| G ′ p,c (|x|): • On the one hand, the singularity of the convolution kernel is removed, which likely simplifies the mathematical analysis,

• On the other hand, with the convolution model, the divergence of the velocity field does not have a definite sign, while (4) imposes

∇ x • V p ≥ 0, ∇ x • V c ≤ 0. Indeed the jacobian matrix of x |x| G ′ p,c (|x|) is x ⊗ x |x| 2 G ′′ p,c (|x|) + I - x ⊗ x |x| 2 G ′ p,c (|x|) |x| . Its trace is G ′′ p,c (|x|) + (N -1) G ′ p,c (|x|)
|x| , which, in general, might change sign (see Appendix B for further details).

Second-order models: repulsive/attractive dynamics 2.2.1 Discrete and kinetic models

In the next step, we think of the attractive/repulsive effect as the action of forces. Such a description naturally arises dealing with gravitational/charged particles [START_REF] Chandrasekhar | Brownian motion, dynamical friction and stellar dynamics[END_REF][START_REF] Nieto | High-field limit of the Vlasov-Poisson-Fokker-Planck system[END_REF]. Let us examine this viewpoint at the discrete level. Now, the preys and chasers dynamics is described by the evolution of the position/velocity pairs (X p,j , Ξ p,j ), and (X c,k , Ξ c,k ), respectively. These quantities are functions of t ≥ 0, valued in R N . The motion of the individuals is then obtained through the fundamental principle of dynamics:

d dt X p,j = Ξ p,j , d dt X c,k = Ξ c,k , d dt Ξ p,j = - 1 τ (Ξ p,j -∇Φ p (t, X p,j )), d dt Ξ c,k = - 1 τ (Ξ c,k -∇Φ c (t, X c,k )). ( 11
)
This ODE system is completed by initial conditions that prescribe the position/velocity pairs at t = 0. In [START_REF] Cucker | Emergent behavior in flocks[END_REF], τ is interpreted as a relaxation time. Individuals are subject to a drag force, which is proportional with opposite direction to their velocity and an interaction force embodied into the potentials Φ p and Φ c . We use the generic definition devised above with a convolution kernel:

Φ p (t, x) = Nc k=1 G p (|x -X c,k |), Φ c (t, x) = - Np j=1 G c (|x -X p,j |).
It means that chasers have a repulsive effect on the preys, which conversely attract the chasers. As said above the definition of the potential by the Poisson equation ( 4) is contained in this framework. Furthermore, we can take into account possible velocity fluctuations. Then the ODEs are replaced by stochastic differential equations, involving independent copies of the Brownian motion

dX p,j = Ξ p,j dt, dΞ p,j = - 1 τ (Ξ p,j -∇Φ p (t, X p,j )) dt + 2θ/τ dB j , dX c,k = Ξ c,k dt, dΞ c,k = - 1 τ (Ξ c,k -∇Φ c (t, X c,k )) dt + 2θ/τ dB k . ( 12
)
Coming back to PDEs, the description of the population now involves distribution functions in phase space that we denote p(t, x, ξ) and c(t, x, ξ) respectively. The integral Ω×O p(t, x, ξ) dξ dx (resp. Ω×O c(t, x, ξ) dξ dx) gives the number of preys (resp. chasers) having at time t ≥ 0 a position x ∈ Ω with a velocity ξ ∈ O. These quantities are driven by the following Vlasov-type equations

∂ t p + ξ • ∇ x p + 1 τ ∇ x Φ p • ∇ ξ p = 1 τ ∇ ξ • (ξp + θ∇ ξ p) ( 13
)
and

∂ t c + ξ • ∇ x c + 1 τ ∇ x Φ c • ∇ ξ c = 1 τ ∇ ξ • (ξc + θ∇ ξ c) ( 14
)
where the potentials are given by

Φ p (t, x) = + G p (x -y)c(t, y, ξ) dξ dy, Φ c (t, x) = - G c (x -y)p(t, y, ξ) dξ dy.
We pass from the individual based description to this kinetic framework by a direct application of Proposition 2.1 in phase space: 13)-( 14). However ( 13)-( 14) can be considered with more general initial distribution functions. When θ = 0, ( 13) and ( 14) are pure Vlasov equations, for θ > 0 they involve the Fokker-Planck operator

p(t, x, ξ) = Np i=1 δ(x = X p,i )⊗δ(ξ = Ξ p,i ) and c(t, x, ξ) = Nc i=1 δ(x = X c,i ) ⊗ δ(ξ = Ξ c,i ) satisfy (
Lf = ∇ ξ • (ξf + θ∇ ξ f ) = ∇ ξ • θM ∇ ξ f M , M (ξ) = e -ξ 2 /(2θ) (2πθ) N/2 ,
which induces some relaxation effects towards the Maxwellian state M (ξ).

For further purposes let us rewrite the equations in dimensionless form

∂ t p + βξ • ∇ x p + γ∇ x Φ p • ∇ ξ p = 1 ǫ ∇ ξ • (ξp + θ∇ ξ p), (15)
and

∂ t c + βξ • ∇ x c + γ∇ x Φ c • ∇ ξ c = 1 ǫ ∇ ξ • (ξc + θ∇ ξ c). ( 16
)
Here, given time and length units T, L, as well as a velocity unit V, the dimensionless parameters are defined by

β = TV L , γ = TV M τ V , ǫ = τ T , θ = θ V 2 . ( 17
)
(We remind that V M is the maximal speed involved in the definition of the kernel G .) Next, we associate to the distribution functions b, c the following "macroscopic" quantities

ρ p (t, x) = p(t, x, ξ) dξ, ρ c (t, x) = c(t, x, ξ) dξ, J p (t, x) = β ξ p(t, x, ξ) dξ, J c (t, x) = β ξ c(t, x, ξ) dξ, P p (t, x) = ξ ⊗ ξ p(t, x, ξ) dξ, P c (t, x) = ξ ⊗ ξ c(t, x, ξ) dξ.
We consider the zeroth and first moments of ( 13) and ( 14). We get

∂ t ρ p + ∇ x • J p = 0, ∂ t ρ c + ∇ x • J c = 0, (18)
and

∂ t J p + β 2 ∇ x • P p -βγρ p ∇ x Φ p = - 1 ǫ J p , ∂ t J c + β 2 ∇ x • P c -βγρ c ∇ x Φ c = - 1 ǫ J c .
The system is not closed because P p and P c cannot be expressed in general by means of ρ p , J p and ρ c , J c , respectively. Nevertheless, it can be used to discuss asymptotic regimes, depending on the respective behavior of the dimensionless parameters β, γ, ǫ, θ.

Macroscopic models: Entropy Minimization Closure

Before dealing with asymptotics issues, let us propose a possible closure. We seek a relevant expression of P p,c depending on ρ p , J p , respectively ρ c , J c . Bearing in mind the dissipation property of the Fokker-Planck operator

ln(f )L(f ) dξ = -4θ 2 M ∇ ξ f M 2 dξ ≤ 0,
we define approximate distribution functions p, c as minimizers of the entropy functional f → f ln(f /M ) dξ under the constraint of fixed zeroth and first order moments, namely (1, βξ)f dξ = (ρ p , J p ) or (ρ c , J c ) respectively. The minimizers have the general form M (ξ) exp(λ + µξ), with λ(t, x) ∈ R, µ(t, x) ∈ R N the Lagrange multipliers associated to the constraints. We close the macroscopic system by setting

P p (t, x) = ξ ⊗ ξ p(t, x, ξ) dξ, P c (t, x) = ξ ⊗ ξ c(t, x, ξ) dξ.
We refer the reader to [START_REF] Coulombel | Diffusion approximation and entropybased moment closure for kinetic equations[END_REF] for analysis of such a closure based on entropy minimization principles [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Levermore | Entropy-based moment closures for kinetic equations[END_REF]. In the present context, the obtained system reads

∂ t J p + β 2 ∇ x • J p ⊗ J p ρ p + β 2 ∇ x ( θρ p ) -γρ p ∇ x Φ p = - 1 ǫ J p , ∂ t J c + β 2 ∇ x • J c ⊗ J c ρ c + β 2 ∇ x ( θρ c ) -γρ c ∇ x Φ c = - 1 ǫ J c , (19)
togteher with [START_REF] Fasano | Equilibrium of two populations subject of chemotaxis[END_REF]. We get nothing but isothermal Euler systems for ρ p , J p and ρ c , J c , coupled through the force terms with the potentials

Φ p (t, x) = G p (x -y)ρ c (t, y) dy, Φ c (t, x) = -G c (x -y)ρ p (t, y) dy.
We can expect that the solutions of this macroscopic system have a behavior close to moments of the original microscopic model, at least for small ǫ's.

Macroscopic models: low and high field limits

In the same spirit as in plasma physics or for attractive dynamics [START_REF] Nieto | High-field limit of the Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Poupaud | Parabolic limit and stability of the Vlasov-Poisson-Fokker-Planck system[END_REF], two asymptotic regimes are relevant for ( 15)-( 16) and make the connection with the models presented in Section 2.1. We are interested in the low field regime, where the coefficients in ( 17) scale as follows

β = γ = 1 √ ǫ , ǫ ≪ 1,
and the high field regime where

β = 1, γ = 1 ǫ , ǫ ≪ 1.
In the former situation we guess that b and c relax toward centered Maxwellians ρ p (t, x)M (ξ), ρ c (t, x)M (ξ), respectively. Plugging this ansatz into the moment equations we arrive at

J p = -θ∇ x ρ p + ρ p ∇ x Φ p , J c = -θ∇ x ρ c + ρ c ∇ x Φ c ,
so that we are led to the following system of convection-diffusion equations

∂ t ρ p + ∇ x • (ρ p ∇ x Φ p -θ∇ x ρ p ) = 0, ∂ t ρ c + ∇ x • (ρ c ∇ x Φ c -θ∇ x ρ c ) = 0. (20)
With only one specie, we obtain by this way the Keller-Segel system [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF], and actually the derivation of this equation from the Vlasov-Fokker-Planck equation dates back to [START_REF] Chandrasekhar | Brownian motion, dynamical friction and stellar dynamics[END_REF] in astrophysics. For the analysis of the low field regime, for both the attractive and repulsive dynamics, we refer the reader to [START_REF] Ghani | Diffusion limit of the Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Goudon | Hydrodynamic limit for the Vlasov-Poisson-Fokker-Planck system: analysis of the two-dimensional case[END_REF][START_REF] Nieto | High-field limit of the Vlasov-Poisson-Fokker-Planck system[END_REF]. The high field regime is much more delicate. Formally the distribution functions are supposed to tend to shifted Maxwellians, and we directly obtain the convection equations [START_REF] Bouchut | Introduction to the mathematical theory of kinetic equations[END_REF], see e. g. [START_REF] Poupaud | Parabolic limit and stability of the Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Goudon | Multidimensional high-field limit of the electrostatic Vlasov-Poisson-Fokker-Planck system[END_REF]. It is likely that these techniques can be adapted to study the coupled problem with two species under attraction/repulsion dynamics; this is beyond the scope of the present paper.

Numerical experiments

We discuss on numerical grounds a few features of the model hierarchy, pointing out very different behaviors depending on the model. The numerical schemes we are going to use are not particularly elaborate, and they do not deserve too much details. They are based on time-splitting strategies: knowing the state of the population, we determine the potential, which is used to update the population on a time step by using standard ODE or PDE solvers, etc. We also discuss the role of the boundary conditions, comparing solutions of the equations set on unbounded domains, or with periodic or Neumann conditions.

Standard and modified 1D models

We start with the ODE system (7) in 1D. Initial positions of the individuals are chosen randomly in [0, 1]. The case where Γ is constant is referred to as the standard model (it corresponds to the convolution with the kernel of the Poisson equation). For the simulation, this constant is .6 for the preys, 1.15 for the chasers. Figures 5 on the top and in the middle present a quite generic organization where the preys split into two groups, one going up, the other going down, in which the individuals have all the same velocity; the chasers all go in the same direction, with the same velocity. This organization appears after a short transient state. When increasing the number of individuals the self-organization establishes faster, as it can be seen by comparing the two figures. Dealing with an even number of preys, we can obtain a situation with the whole population of chasers just stuck between the two groups of preys, see Figure 5 on the bottom. In the case of a periodic setting, we can observe the formation of different types of quasi-periodic patterns as shown on Figure 6. It is worth pointing that the numerical patterns are quite sensitive to the time step: having a fair representation of the solution requires small time steps. Patterns are much more complex when dealing with another convolution kernel, and we can see support effects. In Figure 7 on top the convolution kernel is given by [START_REF] Carrillo | Cross diffusion and nonlinear diffusion preventing blow up in the Keller-Segel model[END_REF] Preys maximal velocity is 1, and the chasers are sensibly faster with a maximal velocity 1.5. However the safety radius of the preys, set to 1, is larger than the safety radius of the chasers, .3. We have ǫ = .4 for the preys, ǫ = 2 for the chasers. For Figure 7 on bottom, the kernels are purely Gaussian (which means that the radius R vanishes in [START_REF] Carrillo | Cross diffusion and nonlinear diffusion preventing blow up in the Keller-Segel model[END_REF]). The other parameters here are as follows: preys' velocity .8, chasers' velocity 1.1, ǫ is .6 for the preys, .5 for the chasers. In both cases we observe bifurcations, with changes either on the amplitude or the direction of the displacement, and possible transient states without motion. We also remark the formation of a cluster of chasers who do not feel the preys, and thus stand at the same place. Figure 8 is concerned with the same situation with periodic boundary conditions. Again, we observe the formation of quasi-periodic patterns. In particular, notice on the Figure on the bottom that a self-organization establishes after a short transient time. Anyway, a wide variety of patterns is possible with this model, depending on the initial configuration. We end this section with simulations of the second-order model [START_REF] Cucker | Emergent behavior in flocks[END_REF]. Typical shapes of the solutions are given in Figure 9 for the standard kernel or the modified kernel. The trajectories of the chasers can be highly involved in the earlier stages. For larger time we observe again segregation effects. In the periodic case, the results are much different from what we observed with the first-order model; they are all of the same type, as displayed in Figure 10

Standard and modified 2D models

For the standard model in 2D, we face the difficulty of the singularity of the convolution kernel: as X c,i becomes close to X p,j then X c,i -X p,j |X c,i -X p,j | 2 blows up. It makes the problem stiff, a difficulty which is well known for similar problems in chemotaxis and astrophysics. Numerically, we regularize a little bit the kernel, and we stop the simulation when the distance between two particles is less than a given (small enough) threshold. Figure 11 shows two typical situations. The two tests are realized with the same conditions: V p = .4, V c = .6; only the initial conditions differ (they are still chosen randomly). Initial positions are marked by the diamond symbol. On the left, one of the prey does not have time to escape; it is rapidly caught by the chaser and the simulation stops. Conversely, on the right we observe that the spreading of the preys is the strategy that allows them to escape to the chasers: they are all far away from the chasers, the motion of the two species becomes very slow. In a periodic box, we can observe similar situations. We see on the left of Figure 12 the situation of preys being caught by the chaser when the velocities are such that V p < V c . However, on the right of Figure 12, when V c < V p , we can see the preys escaping and being pursued in a chaotic way by the chasers. In Figure 13, we find simulations of the 2D modified model. There are initially 34 preys and 20 chasers. In each case the reference velocities are 0.1 for the preys, .5 for the chasers. The vigilance parameters for the preys are R p = 1, ǫ p = .1. For the first test (left) we set R c = .2, ǫ c = .4 for the chasers. These parameters are reduced for the second run (right): R c = .1, ǫ c = .2. Hence chasers are only influenced by preys in their close neighbourhood. We indeed observe that some of the chasers do not move after a while, once the preys are too far. In the periodic case, similar behaviors are shown at Figures 14 and we note that the preys and then the chasers tend to concentrate on the domain boundaries, and to the corners after a while.

Increasing the vigilance parameters could be unefficient from the chasers viewpoint, as shown in Figure 15 (in the classical setting on the left, in the periodic setting on the right). Here we have set V p = .02, R p = .3, ǫ p = .2, V c = .5, R c = 1, ǫ c = 1. In Figure 12: Standard model in the periodic case: Evolution with 18 preys, 7 chasers. On the left, in the case when V p < V c ; on the right when V c < V p this case, chasers are not able to choose between all the preys and they all converge to the 1-median of the set of prey positions (the conjecture can be numerically checked for simple examples, with 3 or 4 preys).

First order PDE models

This Section is concerned with the numerical simulation of the PDE system (2). The equation is completed by initial conditions

ρ p (t = 0, x) = ρ p,Init (x), ρ c (t = 0, x) = ρ c,Init (x).
We shall work in bounded domains, and we also need boundary conditions on the potential. We use either the Neumann conditions

∇ x Φ k • ν = ∂ ν Φ k = 0,
Φ k (x + ℓ) = Φ k (x) for any ℓ ∈ N N , x ∈ [0, 1] N .
As pointed out in Remark 2.2, with such boundary conditions, the Poisson equation is well posed under suitable compatibility conditions on the data. Note also that the important quantity is actually the gradient of the potential Φ k . Therefore, we enforce the existence-uniqueness of the solution of the Poisson equation by the introduction of a global condition and an augmented variable λ k ∈ R:

Ω Φ k dx = 0 and Ω λ k dx = Ω b k dx.
Equations for the potentials Φ k are then modified as follows

   ∆Φ k -λ k = b k Ω Φ k dx = 0
with the appropriate boundary condition (Neumann or periodic) where b p = αρ c , b c = -ρ p and λ k ∈ R is an augmented variable.

The numerical scheme is designed as follows. We consider numerical approximations over simple structured meshes (quadrangles in 2D). The computational domain is denoted Ω h , h > 0 being a measure of the mesh refinement. For given densities at the nth time step, ρ n p and ρ n c , the numerical time step can be decomposed as follows: Computation of the potentials Φ n p and Φ n c by inverting the modified elliptic equations, understood in the following weak formulation

         Ω h ∇ h Φ n c • ∇ h ϕ j dx + Ω h λ c ϕ j dx = Ω h ρ n p ϕ j dx, Ω h Φ n c dx = 0,
and

         Ω h ∇ h Φ n p • ∇ h ϕ j dx + Ω h λ p ϕ j dx = - Ω h αρ n c ϕ j dx, Ω h Φ n p dx = 0,
where ϕ j are piecewise linear test functions, Φ n p = i Φ n p,i ϕ i and Φ n c = i Φ n c,i ϕ i , and ∇ h is a discrete difference operator. The problem can be cast in matrix form

A e e T 0 Φ λ = b 0
with e a certain vector in R P h , and A a P h × P h matrix, P h ∈ N being defined by the refinement of the grid. It can be checked that the linear system is invertible.

Evaluation of the velocities fields V n p and V n c by applying the discrete gradient on potentials computed at the previous step

V n p = ∇ h Φ n p and V n c = ∇ h Φ n c .

Updates of the densities ρ n+1

p and ρ n+1 c by using a mere vertex centered upwind scheme in explicit form:

ρ n+1 p -ρ n p δt + ∇ h • ρ n p V n p = 0 and ρ n+1 c -ρ n c δt + ∇ h • (ρ n c V n c ) = 0.
The definition of the numerical flux is embodied into the discrete divergence operator (∇ h •).

At each time step, boundary conditions are applied, consistently with the numerical test case under consideration.

One-dimensional simulations, Neumann and periodic boundary conditions

We perform a series of numerical tests on the domain [0, 1], as summarised in the following table:

α ρ c,Init (x) ρ p,Init (x) Boundary Case 1D_01 1 √ 10 exp -10 x -1 3 2 √ 5 exp -5 x -2 3 2 Neumann Case 1D_11 1 √ 10 exp -10 x -1 3 2 √ 5 exp -5 x -2 3 2 Periodic Case 1D_12 0.5 2 (1.1 + sin (4πx) ) 2 (1.1 + cos (2πx) ) Periodic
Results are reported in Figures 1617181920for the 1D case. We represent the trajectories of (fictitious) individual particles of the two species. For both species, we pick 6 points in the time-space plane (represented by the bullets) and we follow, backward and forward, the trajectories emanating from these recording points: they appear as bold curves in the picture. (Note that the reconstruction algorithm might fail when the trajectory reaches the boundary.) The first remarkable feature on the picture is the formation of quasi-periodic structures, This is particularly visible with Neumann boundary conditions, see Figures 1617. We do not represent the densities, but we can infer the evolution as follows. With Neumann boundary conditions, the two populations start by moving to the right (boundary x = 1). When peaks form at the far right end of the domain, preys move back to the left (boundary x = 0). Then new bumps of preys are created on the left part of the domain, which, in turn, attract the chasers. The populations move back and forth in this way. The concentration profiles are sensibly different working with periodic boundary conditions, see Figures 181920like if individuals were moving on a ring. In this case, of course, the populations are influenced from both sides. Secondly, we note that trajectories issued from neighboring locations can separate rapidly. This observation is in agreement with the ODEs simulations. We clearly observe the pursuit phenomena, the shape of the streamlines of the preys being more or less reproduced by the chasers, with some delay (compare Figures 16 the shape of the solutions: we do not see a universal profile to establish. Finally, we report the simulations obtained with different grids, the time step being adapted accordingly to guarantee the CFL stability condition of the transport equations. In agreement with what has been observed with ODEs simulation, we point out that time and space scales have to be resolved quite finely: the effect of numerical diffusion can be sensitive on individual trajectories, especially for large times of simulation. It is very likely that the numerical investigation of the asymptotic behavior would require a more involved numerical scheme.

Two-dimensional simulations, periodic boundary conditions

In dimension 2, we work with the following initial data reported1 in Figures 2122. We show the streamlines of the preys and the concentration levels of the two species for several times (the warmer the contour, the higher the concentration). We clearly observe the attraction/repulsion mechanism (compare the left and right columns of Figure 21) and the formation of periodic structures (see Figure 22). To be more specific, the dynamic of the system can be decomposed, in this case, into four main phases of motion: preys move a) toward the point (0.5, 0.5), see Figure 22-left at time 37; b) toward the point (1, 0.5), see Figure 21-top at time 47; c) toward the point (0.5, 1), see Figure 21-middle at time 55; d) and finally toward the point (1, 1) Figure 21-bottom at time 64, before returning to the initial position (0.5, 0.5), see Figure 22-right at time 71. We remind the reader that boundary conditions are periodic such that, for instance, (1, 1) is a translation of the points with coordinates (0, 0), (1, 0) and (0, 1). The motion of the chasers follow this pattern, up to some delay, see Figure 21-right column. Again, we point out that fine enough grids and small enough time steps should be used, otherwise the numerical diffusion may alter the pattern formation.

Analysis of the first order models: Proof of Theorem 2.3

In this section, we focus on the system (2)-( 4): we shall prove the global existenceuniqueness of weak solutions associated to bounded and integrable initial data, and establish that smooth data produce smooth solutions. Derivation of a priori estimate is the necessary preliminary step for proving the existence of weak solutions. The difficulty relies on the lack of regularity of the potentials which are actually defined by convolution formula involving the singular kernel of the Poisson equation. However, in contrast with the attractive equation for a single specie, here the coupling does not induce the formation of singular solution with apparition of Dirac mass in finite time, see [START_REF] Greenberg | Stability results for a diffusion equation with functional drift approximating a chemotaxis model[END_REF][START_REF] Jäger | On explosions of solutions to a system of partial differential equations modelling chemotaxis[END_REF][START_REF] Herrero | Singularity patterns in a chemotaxis model[END_REF][START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences[END_REF]. Instead, the solutions remain bounded at least on any finite time interval. It allows us to construct weak solutions in L 1 ∩L ∞ defined on any time interval [0, T ], 0 < T < ∞. These weak solutions can be shown to be unique. Furthermore, we also show that initial data with bounded gradients produce solutions with the same regularity, in contrast with the situation met with non linear conservation laws of fluid mechanics. The analysis of the coupling (3)-( 4) relies on the following standard statements from harmonic analysis. We refer the reader to [START_REF] Lieb | Analysis[END_REF]Th. 4.3]. We use this claim with λ = N -1, g = ρ given in L r and f spans L p ; it leads to estimate the velocity V = ∇E ⋆ ρ in L p ′ , 1/p ′ = 1/r -1/N , as soon as the concentration ρ belongs to L r with 1 < r < N . This can be improved in the present situation, since we are dealing with bounded and integrable densities.

Lemma 4.1 (Hardy-Littlewood-Sobolev inequality)

Let 1 < p, r < ∞ and 0 < λ < N . Assume 1/p + 1/r = 2 -λ/N . There exists a constant C > 0 such that for any f ∈ L p (R N ) and g ∈ L r (R N ). we have f (x)g(y) |x -y| λ dy dx ≤ C f L p g L r .
Lemma 4.2 Let g ∈ L 1 ∩ L ∞ (R N ). Then V (x) = g(y) |x-y| N -1 dy belongs to L ∞ (R N ) with V L ∞ ≤ N (N -1) 1/N -1 |S N -1 | 1-1/N g 1-1/N L ∞ g 1/N L 1 .
Proof. We estimate by splitting the convolution integral, with some M > 0:

|V (x)| ≤ |x-y|≥M ... dy + |x-y|<M ... dy ≤ M 1-N g L 1 + g L ∞ |S N -1 | M 0 r N -1 dr r N -1 = M 1-N g L 1 + |S N -1 |M g L ∞ .
We conclude by optimizing with respect to M . Finally, derivatives of the velocity field are estimated with the following statement, see e. g. [41, Chap. VI], [START_REF] Duoandikoetxea | Fourier Analysis[END_REF]Th. 4.12].

Lemma 4.3 (Calderon-Zygmung inequality) There exists

M ⋆ > 0 such that for any 1 < p < ∞ and any g ∈ L p (R N ), the function V (x) = g(y) |x-y| N -1 dy satisfies ∇V L p ≤ M ⋆ p g L p .
The keypoint of the analysis relies on the following a priori estimate.

Proposition 4.4 (A priori estimates) Assume that the initial data satisfy:

ρ p,Init ≥ 0, ρ c,Init ≥ 0, ρ p,Init , ρ c,Init ∈ L 1 ∩ L ∞ (R N ).
Then smooth solutions of (2)-( 4) satisfy

0 ≤ ρ p (t, x) ≤ ρ p,Init L ∞ , 0 ≤ ρ c (t, x) ≤ ρ c,Init L ∞ e t ρ p,Init L ∞ .

29

Proof. We write the equations in the non-conservative form

∂ t ρ p + ∇ x Φ p • ∇ x ρ p = -ρ p ∆ x Φ p = -αρ p ρ c , ∂ t ρ c + ∇ x Φ c • ∇ x ρ c = -ρ c ∆ x Φ c = +ρ p ρ c
Let us denote by X p (resp. X c ) the characteristic curves associated to ∇ x Φ p (resp. ∇ x Φ c ). This assumes the regularity of each potential. We obtain

ρ p (t, x) = ρ p,Init (X p (0; t, x)) exp -α t 0 ρ c (s, X p (s; t, x)) ds and ρ c (t, x) = ρ c,Init (X c (0; t, x)) exp t 0 ρ p (s, X c (s; t, x)) ds .
Clearly, that implies ρ p , ρ c ≥ 0. Then we deduce also that ρ p is dominated by ρ p,Init L ∞ , and finally we obtain the estimate for ρ c .

We turn to the proof of Theorem 2.3. Let 0 < T < ∞ be fixed once for all. With the information in Proposition 4.4, and owing to the conservative nature of the equations, we seek solutions of ( 2)-( 4) in the following closed and convex set

C T = (ρ p , ρ c ) ∈ L ∞ (0, T ; L 1 ∩ L ∞ (R N )), 0 ≤ ρ p (t, x) ≤ ρ p,Init L ∞ , 0 ≤ ρ c (t, x) ≤ ρ c,Init L ∞ e T ρ p,Init L ∞ , ρ p (t, x) dx = ρ p,Init (x) dx, ρ c (t, x) dx = ρ c,Init (x) dx .
We are going to use a fixed point argument, based on the Schauder theorem. The proof relies heavily on the analysis of transport equations, and renormalization techniques [START_REF] Di Perna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF]. The proof splits into several steps:

• To obtain the stability and compactness of sequences of solutions, we need to consider first a regularized (non linear) problem, the solutions of which are obtained as fixed points of an application defined through linear transport equations.

• Then, we let the regularizing parameter go to 0.

• Finally, we establish the uniqueness of the weak solutions of the non linear problem. We shall also investigate the regularity.

We introduce a sequence of mollifiers

ζ ǫ (x) = 1 ǫ N ζ(x/ǫ), ζ ∈ C ∞ c , 0 ≤ ζ(x) ≤ 1, ζ(x) dx = 1.
Given (ρ p , ρ c ) ∈ C T , we define (ρ p , ρc ) as to be the solution of the following linear PDEs system

∂ t ρp + V ǫ p • ∇ x ρp = -α(ζ ǫ ⋆ ρ c )ρ p , ∂ t ρc + V ǫ c • ∇ x ρc = (ζ ǫ ⋆ ρ p )ρ c , ρp (0, x) = ρ p,Init (x), ρc (0, x) = ρ c,Init (x),
where the velocities are given by

V ǫ p = ∇ x Φ ǫ p , V ǫ c = ∇ x Φ ǫ c , ∆Φ ǫ p = αζ ǫ ⋆ ρ c , -∆Φ ǫ c = ζ ǫ ⋆ ρ p .
Note that

ζ ǫ ⋆ρ p L ∞ ≤ ρ p L ∞ , ζ ǫ ⋆ρ p L 1 ≤ ρ p L 1 , ζ ǫ ⋆ρ c L ∞ ≤ ρ c L ∞ , ζ ǫ ⋆ρ c L 1 ≤ ρ c L 1 .
According to Lemma 4.2, we infer that V ǫ p and V ǫ c are bounded in L ∞ ((0, T ) × R N ), uniformly with respect to ǫ. Furthermore ∇

x • V ǫ p = αζ ǫ ⋆ ρ c and ∇ x • V ǫ c = -ζ ǫ ⋆ ρ p are bounded in L ∞ (0, T ; L 1 ∩L ∞ (R N ))
, and by Lemma 4.3, the components of the jacobian matrices ∇ x V ǫ p and ∇ x V ǫ c are bounded in L ∞ (0, T ; L q (R N )), for any 1 < q < ∞, uniformly with respect to ǫ. This is not enough to make use of the classical theory of characteristic curves; nevertheless we can apply the results of [START_REF] Di Perna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] which guarantee that ρp , ρc are well defined. We also refer the reader to the fully detailed analysis of the conservative transport equation in [4, Chap. VI, Sect. 1]. Let us denote T : (ρ p , ρ c ) → (ρ p , ρc ). This mapping satisfies T (C T ) ⊂ C T (the L ∞ estimate are obtained by adapting readily the arguments in Proposition 4.4). Furthermore, still as a direct application of results in [START_REF] Di Perna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF], see also [4, Th. VI.1.9], T is continuous, with C T endowed with the L 1 ((0, T ) × R N ) norm. Let us now check that T is a compact mapping. Let (ρ p,n , ρ c,n ) n∈N be a sequence in C T . It is thus bounded in L ∞ (0, T ; L 1 ∩ L ∞ (R N )). We can extract a subsequence, still labelled ρ p,n and ρ c,n , which converges weakly-⋆ in L ∞ ((0, T ) × R N ), and weakly in L 1 ((0, T ) × R N ). We denote by ρ p and ρ c the limits. For any fixed ǫ, we have

V ǫ p,n ⇀ V ǫ p and V ǫ c,n ⇀ V ǫ c
weakly in any L q ((0, T ) × R N ), for 1 ≤ q < ∞, and weakly-⋆ in L ∞ ((0, T ) × R N ). By virtue of Lemma 4.2, and 4.3, V ǫ p,n and V ǫ c,n are bounded in L ∞ ((0, T ) × R N )), with the jacobian matrices ∇ x V ǫ p,n and ∇ x V ǫ c,n bounded in L ∞ (0, T ; L q (R N )), for any 1 < q < ∞, uniformly with respect to ǫ. Consequently the following "space-compactness" property holds

lim h→0 sup n,ǫ |V ǫ p,n (t, x + h) -V ǫ p,n (t, x)| p + |V ǫ c,n (t, x + h) -V ǫ c,n (t, x)| p dx dt = 0.
Besides, ρp,n and ρc,n are also bounded in L ∞ (0, T ; L 1 ∩ L ∞ (R N )). Therefore, we can assume that ρp,n and ρc,n converge weakly-⋆ in L ∞ ((0, T ) × R N ), and weakly in L 1 ((0, T ) × R N ) to ρp and ρc , respectively. Moreover, ∂ t ρp,n = -∇ x • (V ǫ p,n ρp,n ), and

∂ t ρc,n = -∇ x • (V ǫ c,n ρcn ) are bounded in L ∞ (0, T ; W -1,∞ (R N ))
. This allows to pass to the limit in the products

V ǫ p,n ρp,n ⇀ V ǫ p ρp , V ǫ c,n ρc,n ⇀ V ǫ c ρc , at least in D ′ ((0, T ) × R N ) as n → ∞.
We refer the reader to [14, Th. II.7] or [START_REF] Lions | Compressible models[END_REF]Lemma 5.1] for details on such compactness arguments. Hence the limits satisfy the transport equations ∂ t ρp +∇ x •(V ǫ p ρp ) = 0, and ∂ t ρc +∇ x •(V ǫ c ρc ) = 0. For applying the fixed point argument, we need to establish that ρp,n and ρc,n converge strongly. This property requires further compactness on the divergence of the velocity fields. This is where the regularization plays a role. Indeed, for any fixed ǫ > 0,

∇ x • V ǫ p,n , ∇ x • V ǫ c,n satisfy lim h→0 sup n |∇ x • V ǫ p,n (t, x + h) -∇ x • V ǫ p,n (t, x)| p dx dt + |∇ x • V ǫ c,n (t, x + h) -∇ x • V ǫ c,n (t, x)| p dx dt = 0,
and we justify that ρp,n and ρc,n converge strongly in L 1 ((0, T ) × R N ) to ρ p and ρ c by applying again the renormalization techniques of [START_REF] Di Perna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] and proving that ρ 2 p,n , ρ 2 c,n converge weakly to ρ 2 p , ρ 2 c , see the second step of the proofs of [START_REF] Di Perna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF]Th. II.4], or [4, Th.VI.1.9]. We have thus proved that T satisfies the hypothesis of the Schauder theorem on the set C T . Hence, T admits a fixed point (ρ ǫ p , ρ ǫ c ) = T (ρ ǫ p , ρ ǫ c ), which is precisely a solution of the nonlinear system

∂ t ρ ǫ p + ∇ x • (V ǫ p ρ ǫ p ) = 0, ∂ t ρ ǫ c + ∇ x • (V ǫ c ρ ǫ c ) = 0 with V ǫ p = ∇ x Φ ǫ p , V ǫ c = ∇ x Φ ǫ c , ∆Φ ǫ p = αζ ǫ ⋆ ρ ǫ c , -∆Φ ǫ c = ζ ǫ ⋆ ρ ǫ p .
It remains to let ǫ go to 0. But ρ ǫ p and ρ ǫ c lie in C T , thus are bounded in L ∞ (0, T ; L 1 ∩ L ∞ (R N )), and it suffices to repeat the compactness argument used above to remove the regularization. It proves the existence of weak solutions for data in L 1 ∩ L ∞ . Continuity with respect to time is a consequence of renormalization techniques, see [START_REF] Di Perna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF]Th. II.3] or [START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF]Th. VI.1.3]. Now we turn to uniqueness. We adapt the tricky argument introduced in [START_REF] Robert | Unicité de la solution faible à support compact de l'équation de Vlasov-Poisson[END_REF] for the Vlasov-Poisson equation. Let (ρ p i , ρ c i ), i = 1, 2, be two solutions in L ∞ (0, T ; L 1 ∩ L ∞ (R N )), corresponding to the same initial data. Here we slightly modify the notations.

Let us set R = R 1 -R 2 := ρ p 1 -ρ p 2 , S = C 1 -C 2 := ρ c 1 -ρ c 2 , W = Φ p 1 -Φ p 2 and Z = Φ c 1 -Φ c 2 We thus have ∂ t R + ∇ x • (V p 1 R + R 2 ∇ x W ) = 0 = ∂ t S + ∇ x • (V c 1 S + C 2 ∇ x Z), with ∆ x W = αS, ∆ x Z = -R, R(0, •) = S(0, •) = 0.
We compute

1 2 d dt |∇ x W | 2 + |∇ x Z| 2 dx = -αW ∂ t S + Z∂ t R dx = - α∇ x W • (V c 1 S + C 2 ∇ x Z) -∇ x Z • (V p 1 R + R 2 ∇ x W ) dx = -I -(αC 2 -R 2 ) ∇ x W • ∇ x Z dx where I = α∇ x W • V c 1 S -∇ x Z • V p 1 R dx = ∇ x W • V c 1 ∆ x W dx + ∇ x Z • V p 1 ∆ x Z dx.
Observe that

∇ x W • V c 1 ∆ x W dx = - i,j ∂ 2 x i ,x j W V (i) c 1 ∂ x j W dx + ∂ x i W ∂ x j V (i) c 1 ∂ x j W dx .

32

Integrating by part again, we obtain i,j

∂ 2 x i ,x j W V (i) c 1 ∂ x j W dx = - i,j ∂ x j W ∂ x i V (i) c 1 ∂ x j W dx - i,j ∂ x j W V (i) c 1 ∂ 2 x i ,x j W dx = - 1 2 i,j ∂ x j W ∂ x i V (i) c 1 ∂ x j W dx = + 1 2 R 1 |∇ x W | 2 dx.
A similar manipulation applies for the last term in the expression of I. Therefore, we arrive at

I = - 1 2 R 1 |∇ x W | 2 dx + α 2 C 1 |∇ x Z| 2 dx -D 2 Φ c 1 ∇ x W • ∇ x W dx -D 2 Φ p 1 ∇ x Z • ∇ x Z dx
where D 2 Φ stands for the Hessian matrix of the potential function Φ. Let us set

E (t) = |∇ x W | 2 + |∇ x Z| 2 dx.
Bear in mind that there exists C T > 0 such that 0 ≤ E (t) ≤ C T , by virtue of the known estimates on the R i 's and C i 's, and the Hardy-Littlewood-Sobolev inequality.

We wish to establish that the differential inequality

d dt E ≤ M pE 1-1/p holds for any p ≥ 2, 0 ≤ t ≤ T , with a constant M which depends on T , α, N , R 1,2 L 1 , R 1,2 L ∞ , C 1,2 L 1 , C 1,2 L ∞ ,
but that remains uniform with respect to p ≥ 2. We shall use the shorthand notation A B when the inequality A ≤ M B holds with such a constant M , without trying to make the definition of the constant M , which might vary from a line to another, precise. By using the Cauchy-Schwarz inequality together with |ab| ≤ 1 2 (a 2 + b 2 ), we obtain

d dt E E + D 2 Φ c 1 ∇ x W • ∇ x W dx + D 2 Φ p 1 ∇ x Z • ∇ x Z dx II .
By using the Hölder and Calderon-Zygmund inequalities, for any 1 < p < ∞, 1/p ′ + 1/p = 1, the last two terms can be dominated by

II ≤ D 2 Φ c 1 L p |∇ x W | 2p ′ dx 1/p ′ + D 2 Φ p 1 L p |∇ x Z| 2p ′ dx 1/p ′ ≤ M ⋆ p R 1 L p ∇W 2/p L ∞ |∇ x W | 2 dx 1-1/p + C 1 L p ∇Z 2/p L ∞ |∇ x Z| 2 dx 1-1/p . Note that ρ L p ≤ ρ 1-1/p L ∞ ρ 1/p L 1 ≤ ρ L ∞ + ρ L 1 can
L ∞ ≤ (1 + ∇ x Φ p 1 L ∞ + ∇ x Φ p 2 L ∞ )
1, by virtue of Lemma 4.2. Of course, a similar estimate holds with ∇Z. Hence, we get

II p |∇ x W | 2 dx 1-1/p + |∇ x Z| 2 dx 1-1/p .
The elementary inequality a r + b r ≤ 2 1-r (a + b) r ≤ 2(a + b) r which holds for any a, b ≥ 0, 0 < r < 1 yields

d dt E E + pE 1-1/p pE 1-1/p ,
for any p ≥ 2, since we already know that E is bounded. We remind that E (0) = 0. For given ǫ > 0, we compare E with the solution of the ODE

d dt z ǫ (t) = M p(ǫ + z ǫ (t)) 1-1/p , with z ǫ (0) = ǫ, that is z ǫ (t) = ((2ǫ) 1/p + M t) p -ǫ.
We have E (t) ≤ z ǫ (t) for any ǫ > 0 and letting ǫ go to 0 we deduce that E ≤ (M t) p holds for any 2 ≤ p < ∞. By letting p go to ∞ we infer that E (t) = 0 for 0 ≤ t ≤ 1/M . We repeat the argument on successive time intervals of length 1/M , and we conclude that E vanishes on the whole interval [0, T ]. It implies ∇W = ∇Z = 0, and thus R = S = 0.

We end with the regularity analysis. Lemma 4.3 does not provide uniform bound on derivatives of the velocity fields. To this end, we need the following classical statement [2, Lemma 4.2].

Lemma 4.5 Let

ρ ∈ L 1 ∩ L ∞ (R N ) with ∇ x ρ ∈ L ∞ (R N ). Let us set Φ = E ⋆ ρ.
Then there exists M > 0, which does not depend on ρ, such that

D 2 Φ L ∞ ≤ M 1 + ρ L 1 + ρ L ∞ (1 + ln(1 + ∇ x ρ L ∞ )) Let us denote β i = ∂ x i ρ p , γ i = ∂ x i ρ c and β = (β 1 , ..., β N ), γ = (γ 1 , ..., γ N ). By differentiating (2) we obtain (∂ t + V p • ∇ x )β i + ∂ x i V p • β = -α(ρ c β i + ρ p γ i ), (∂ t + V c • ∇ x )γ i + ∂ x i V c • γ = (ρ p γ i + ρ c β i )
Applying usual estimates for transport equation leads to

β i L ∞ (t) β i L ∞ (0) + t 0 β L ∞ (s)(1 + ∂ x i V p L ∞ (s)) + γ i L ∞ (s) ds, γ i L ∞ (t) γ i L ∞ (0) + t 0 γ L ∞ (s)(1 + ∂ x i V c L ∞ (s)) + β i L ∞ (s) ds.
Here and below indicates an estimate involving a constant depending only on T , α and the L 1 and L ∞ norms of the initial data. Since

∂ x i V p L ∞ (1 + ln(1 + ∇ρ c L ∞ )), ∂ x i V c L ∞ (1 + ln(1 + ∇ρ p L ∞ )) (21)
and using a ln(1 + b)

+ b ln(1 + a) ≤ (a + b) ln(1 + a + b), we deduce that G = β L ∞ (t) + γ L ∞ (t) satisfies G (t) G (0) + t 0 G (s)(1 + ln(1 + G (s))) ds. Let t → z(t) be the solution of the ODE d dt z = M (1 + z)(1 + ln(1 + z)), z(0) = z Init > G (0)
. As a matter of fact, we find z(t) = exp(Ae M t -1) -1, with A = 1 + ln(1 + z Init ). Then, by comparizon principle, G (t) ≤ z(t) holds for any 0 ≤ t ≤ T < ∞. Hence ∇ x ρ p and ∇ x ρ c are bounded in L ∞ , and so is the jacobian of the velocities, using [START_REF] Greenberg | Stability results for a diffusion equation with functional drift approximating a chemotaxis model[END_REF].

The proof can be readily adapted to the case where ( 2) is coupled to [START_REF] Crandall | Some relations between non-expansive and order preserving mappings[END_REF]. [START_REF] Crandall | Some relations between non-expansive and order preserving mappings[END_REF] and initial data

• G p,c , G ′ p,c and G ′′ p,c belong to L ∞ (R), • G ′ p,c ≥ 0. We further assume that r → r N -1 G ′′ (r) ∈ L 1 (R). Let ρ p,Init , ρ c,Init be non nega- tive functions in L 1 ∩ L ∞ (R N ). Let 0 < T < ∞. Then, there exists a unique weak solution (ρ p , ρ c ) ∈ L ∞ (0, T ; L 1 ∩ L ∞ (R N )) of (2), with
(ρ p,Init , ρ c,Init ). If, moreover, ∇ x ρ p,Init , ∇ x ρ c,Init belong to L ∞ (R N ), then the solution lies in L ∞ (0, T ; W 1,∞ (R N )).
According to the computation in Appendix B, we have the following estimate for the velocity field. The proof follows directly from the formulae detailed in Appendix B, and standard convolution inequalities.

Lemma 4.7 Let

ρ ∈ L 1 (R N ), ρ ≥ 0 and set Φ(x) = G (|x -y|)ρ(y) dy, where G ∈ C 1 (R N
) is an non decreasing function with bounded first and second derivatives. Then,

V = ∇Φ belongs to L ∞ (R N ), with V L ∞ ≤ G ′ L ∞ ρ L 1 . Furthermore, we have • If N = 1, then ∂ x V = d -2G ′ ρ, with d L ∞ ≤ G ′′ L ∞ ρ L 1 , • If N > 1, then ∇ x • V = d 1 + d 2 , with d 1 L ∞ ≤ G ′′ L ∞ ρ L 1 and d 2 ∈ L 1 (R N ), d 2 ≥ 0. If we assume moreover ρ ∈ L ∞ (R N ), then, V lies in W 1,∞ .
The analog of Proposition 4.4 for (10) casts as follows.

ρ p,Init ≥ 0, ρ c,Init ≥ 0, ρ p,Init , ρ c,Init ∈ L 1 ∩ L ∞ (R N ).
Smooth solutions of (2) with (10) satisfy

0 ≤ ρ p (t, x) ≤ ρ p,Init L ∞ e t G ′′ p L ∞ ρ c,Init L 1 , 0 ≤ ρ c (t, x) ≤ ρ c,Init L ∞ e t( G ′′ c L ∞ ρ p,Init L 1 +Λ ρp L ∞ ((0,t)×R N ) ) , with Λ the L 1 norm of x → G ′ c (|x|) |x| when N > 1, or Λ = 1 when N = 1.
Proof. We consider the case where N > 1. We write the equations in the nonconservative form

∂ t ρ p + V p • ∇ x ρ p = -ρ p ∇ x • V p ≤ -pd p , ∂ t ρ c + V c • ∇ x ρ c = -ρ c ∇ x • V c where d p L ∞ ≤ G ′′ p L ∞ ρ c (t, •) L 1 = G ′′ p L ∞ ρ c,Init L 1 . It provides the bound on ρ p (t, x). Then, we get |∇ x • V c | ≤ G ′′ c L ∞ ρ p,Init L 1 + Λ ρ p L ∞ ((0,t)×R N )
, which ends the proof. The adaptation to the case N = 1 is straightforward.

The existence proof reproduces the argument of the model coupled with the Poisson equation, up to the modification of the set C T which has to account for the estimates detailed in Proposition 4.8. Note also that in this context, the velocity field is Lipschitz with respect to the space variable, which allows to define characteristics. Consequently solutions associated to data in W 1,∞ (R N ) data belong to L ∞ (0, T ; W 1,∞ (R N )) too. Let us consider the situation N > 1. Note that space-compactness of the velocity fields and their divergence is guaranteed without requiring any further regularization recipe, since x → G ′ (|x|) |x| and x → G ′′ (|x|) are integrable. This is a consequence of the following estimate U (x + h -y)ρ(y) dy -U (x -y)ρ(y) dy ≤ ρ L ∞ |U (z + h) -U (z)| dz which tends to 0 as |h| → 0, by virtue of the continuity of translation in L 1 , uniformly over functions ρ lying in a given ball of L ∞ (R N ).

When the velocity V is given, the solution mapping ρ Init → ρ, the solution of ∂ t ρ+∇ x •(V ρ) = 0 with ρ(t = 0, x) = ρ Init (x) is both order preserving --ρ Init,1 ≥ ρ Init,2 implies ρ 1 ≥ ρ 2 --and integral preserving --ρ(t, x) dx = ρ Init (x) dx. Therefore, the Crandall-Tartar Lemma [START_REF] Crandall | Some relations between non-expansive and order preserving mappings[END_REF] directly implies the contraction property |ρ 1 (t, x)ρ 2 (t, x)| dx ≤ |ρ Init,1 (x) -ρ Init,2 (x)| dx. In particular, this inequality applies to the fixed point of the mapping T , which leads to Consequently, uniqueness of the solution of the non linear problem is guaranteed. For N = 1, the proof of the Poisson case can be readily adapted; details are left to the reader.

Conclusion

We have introduced a hierarchy of models describing the motion of two populations the interactions of which are driven by pursuit-evasion principles. It leads to define crossed self-consistent interaction potentials. Namely, we consider the concentration of preys and chasers, ρ p (t, x) and ρ c (t, x) respectively. For individual-based models they can be thought of as sums of Dirac masses. We associate to the concentrations, the potentials Φ p (t, •) = G p * ρ c (t, •) and Φ c (t, •) = -G c * ρ p (t, •), for certain radially symmetric kernels G p,c . Roughly speaking the action of a population on the other one is embodied into the gradients ∇Φ p,c . The sign in the definition of the potential determines whether the effect is repulsive (chasers are repellent for the preys) or attractive (preys are attractive for the chasers). The definition of the potentials from Poisson equations, as it comes in the standard theory of chemotaxis, is a particular case of this situation. We design a hierarchy of models that can be summarized in the following As the number of individuals increases, the computational effort to reproduce the behavior of individual-based models becomes non affordable. It motivates the use of more macroscopic equations, which are also able to capture self-organization mechanisms relevant in life sciences. We have analyzed in detail the first order model ( 2) either with the Poisson kernel or a more general convolution relation. In contrast with similar models for a single specie, blow up do not occur in finite time in the present situation. The numerical simulations confirm the ability of the models in producing remarkable pattern formation. While it is likely that the models can be enriched by taking into account further interaction effects, the analysis of typical patterns, their stability and the influence of the modelling parameters should motivate further mathematical works. It proves i). We obtain ii) by integrating between s = 0 and s = t the following identity: d ds ρ(s, X(s; t, x)) = -ρ(s, X(s; t, x)) × (∇ x • u)(s, X(s; t, x)).

A Proof of

B Computation of the divergence of the velocity field

The velocity field is defined by the convolution

V (x) = E ⋆ ρ(x), E (x) = x |x| G ′ (|x|) = ∇ x (G (|x|).
The function G : R → R is supposed to belong to W 2,∞ , with bounded first and second derivatives. As a matter of fact, E ∈ L ∞ (R N ). Furthermore its derivation far from the origin does not lead to any difficulty and the derivative is given by the matrix valued function

Ψ(x) = x ⊗ x |x| 2 G ′′ (|x|) + I - x ⊗ x |x| 2 G ′ (|x|) |x| .
Observe that Ψ belongs to L ∞ (R N ) too. Hence ∇E is computed as follows: for any We conclude that

∇E = Ψ if N > 1, Ψ -2G ′ δ(x = 0) = G ′′ -2G ′ δ(x = 0) if N = 1.
For N > 1, the divergence of the velocity field is therefore given by the convolution with Tr(Ψ) = G ′′ (|x|) + (N -1) G ′ (|x|) |x| .

Figure 1 :Figure 2 :

 12 Figure 1: The chaser running after the prey: numerical simulation with one individual of each specie.

Figure 3 :

 3 Figure 3: A set of chasers (left) running after a set of preys (right) with V p < V c : the population of preys splits at the meeting point.

Figure 4 :

 4 Figure 4: Typical profile of the velocity function (6).

  preys, 12 chasers, final time 0.1.

Figure 5 :

 5 Figure 5: Standard model : evolution with different numbers of individuals. The chasers are displayed on the left and the preys on the right

Figure 6 :

 6 Figure 6: Standard model in the periodic setting: evolution with different numbers of individuals. The chasers are displayed on the left and the preys on the right

  kernel. Evolution with 10 preys, 14 chasers, final time 1.

Figure 7 :

 7 Figure 7: Modified model. The chasers are displayed on the left and the preys on the right

Figure 8 :

 8 Figure 8: Modified model in the periodic setting. The chasers are displayed on the left and the preys on the right.

Figure 9 :

 9 Figure 9: Second order model with the standard (top) and the modified (bottom) kernel.

Figure 10 :

 10 Figure 10: Second order model with the standard kernel in the periodic case.

Figure 11 :

 11 Figure 11: Standard model: Evolution with 18 preys, 7 chasers.

Figure 13 :

 13 Figure 13: Modified model: Evolution with 34 preys, 20 chasers, final time 2

Figure 14 : 13 -

 1413 Figure 14: Modified model in the periodic case: Evolution with 34 preys, 20 chasers. The parameters are the same as the ones of Figure 13

Figure 15 :

 15 Figure 15: Modified model: Evolution with 34 preys, 11 chasers with a large vigilance radius, final time 5. On the left, on the unbounded domain and on the right, on the periodic domain

Figure 16 :

 16 Figure 16: Case 1D_01: Neumann boundary condition. Trajectories of the preys for N h = 400 (left) and N h = 800 (right).

Figure 17 :

 17 Figure 17: Case 1D_01: Neumann boundary condition. Trajectories of the chasers for N h = 400 (left) and N h = 800 (right).

Figures

  Figures 20) clearly influence the shape of the solutions: we do not see a universal profile to establish. Finally, we report the simulations obtained with different grids, the time step being adapted accordingly to guarantee the CFL stability condition of the transport equations. In agreement with what has been observed with ODEs simulation, we point out that time and space scales have to be resolved quite finely: the effect of numerical diffusion can be sensitive on individual trajectories, especially for large times of simulation. It is very likely that the numerical investigation of the asymptotic behavior would require a more involved numerical scheme.

  Figures 20) clearly influence the shape of the solutions: we do not see a universal profile to establish. Finally, we report the simulations obtained with different grids, the time step being adapted accordingly to guarantee the CFL stability condition of the transport equations. In agreement with what has been observed with ODEs simulation, we point out that time and space scales have to be resolved quite finely: the effect of numerical diffusion can be sensitive on individual trajectories, especially for large times of simulation. It is very likely that the numerical investigation of the asymptotic behavior would require a more involved numerical scheme.

  01 + 0.8 sin (πx) sin (πy) ) .The domain is the square [0, 1] × [0, 1]. We set α = 0.5. The number of grid points in each space direction is N h = 100. Results with periodic boundary conditions are

Figure 18 :

 18 Figure 18: Case 1D_11: Periodic boundary condition. Trajectories of the preys for N h = 400 (left) and N h = 800 (right).

Figure 19 :

 19 Figure 19: Case 1D_11: Periodic boundary condition. Trajectories of the chasers for N h = 400 (left) and N h = 800 (right).

Figure 20 :

 20 Figure 20: Case 1D_12: Periodic boundary condition. Trajectories of the preys (left) and the chasers (right). N h = 500.

Figure 21 :

 21 Figure 21: 2D simulations with periodic boundary conditions. Streamlines of the preys (left), concentration of the preys (middle), concentration of the chasers (right) at t = 47 (top), t = 55 (middle) and t = 63 (bottom).

Figure 22 :

 22 Figure 22: 2D simulations with periodic boundary conditions. Streamlines of the preys at t = 37 and t = 71.

Theorem 4 . 6

 46 Let G p,c : R → R fulfill the following conditions

Proposition 4 . 8 (

 48 A priori estimates) Let G p , G c satisfy the assumptions of Lemma 4.7. The initial data are required to satisfy

|ρ p 1

 1 (t, x) -ρ p 2 (t, x)| dx + |ρ c 1 (t, x) -ρ c 2 (t, x)| dx ≤ |ρ p,Init,1 (x) -ρ p,Init,2 (x)| dx + |ρ c,Init,1 (x) -ρ c,Init,2 (x)| dx.

Proposition 2. 1 .

 1 Let ϕ ∈ C ∞ c ([0, ∞) × R N ). By definition of distributional derivatives, we compute∂ t I i=1 δ(x = X(t; 0, x 0,i )) + ∂ x u I i=1 δ(x = X(t; 0, x 0,i )) , ϕ = -I i=1 δ(x = X(t; 0, x 0,i )), ∂ t ϕ + u∂ x ϕϕ + u∂ x ϕ)(t, X(t; 0, x 0,i )) dt -, X(t; 0, x 0,i )) dt -I i=1ϕ(0, x 0,i ) = 0.

S N - 1 G

 1 ϕ ∈ C ∞ c (R N ), ∇E , ϕ = -lim ǫ→0 |x|>ǫ E (x) ∇ϕ(x) dx = lim ǫ→0 |x|>ǫ ∇E (x) ϕ(x) dx + B(ǫ)with B(ǫ) the boundary term, given by an integral over {|x| = ǫ}, the expression of which depends on the space dimension. The dominated convergence theorem yields limǫ→0 |x|>ǫ ∇E (x) ϕ(x) dx = Ψ(x) ϕ(x) dx.Consider the boundary term in dimension N > 1; it reads (with ν(x) the outward unit vector on the sphere {|x| = ǫ}, that isx |x| ) B(ǫ) = -|x|=ǫ E (x)ϕ(x)ν(x) dσ(x) = ′ (ǫ)ϕ(ǫω) ǫ N -1 dω --→ ǫ→0 0.In dimension N = 1, the boundary term becomesB(ǫ) = -G ′ (ǫ)(ϕ(ǫ) + ϕ(-ǫ)) --→ ǫ→0 -2G ′ (0)ϕ(0).

  .

								Predators												Preys
		0.9												1.5
		0.8
														1
		0.7
		0.6												0.5
	position	0.5											position
		0.4												0
		0.3
														-0.5
		0.2
		0.1	0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2	-1	0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2
								time												time
							Evolution with 9 preys, 5 chasers, final time 0.2.
		0.1	0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1

  and do not change significantly with the kernel.

							Predators												Preys
	2												15
	1
													10
	0
													5
	-1
	-2												0
	-3
													-5
	-4
													-10
	-5
	-6	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	-15	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8
							Predators												Preys
	8												12
	7												10
													8
	6
													6
	5
													4
	4
													2
	3
													0
	2
													-2
	1												-4
	0	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	-6	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8

  on ∂Ω.

						Preys and Predators										Preys and Predators
	1												1
	0.9												0.9
	0.8												0.8
	0.7												0.7
	0.6												0.6
	0.5												0.5
	0.4												0.4
	0.3												0.3
	0.2												0.2
	0.1												0.1
	0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

  and 17, Figures 18 and 19, Figures 20 left and right). The boundary conditions (compare Figures 16-17 with Figures 18-19) and the initial data (compare Figures 18-19 with

  be dominated independently of the exponent p. Similarly, restricting to p ≥ 2, we have the rough estimate ∇W

	2/p

  table:

	2nd order models ∇Φ is a force field ODEs in phase space	-→ Large friction	1st order models ∇Φ is a velocity field ODEs for positions
	(11) or (12)		(9)
	Kinetic models
	(13)-(14)
	Macroscopic model		Drift Diffusion models
	for concentrations/momentum		for concentrations
	Entropy Minimization Principle		("crossed Keller-Segel system")
	(19)		(20)
			Purely convective model
			(2)

This simple example shows that the modelling is questionable, already in 1D. In this model, the velocity is only determined by sgn(X p -X c ): the direction depends

Movies can be consulted at the URL http://www-sop.inria.fr/members/Thierry.Goudon/movieBC.mpeg (periodic boundary conditions) and http://www-sop.inria.fr/members/Thierry.Goudon/movieN.mpeg (Neumann boundary conditions).