J A Carrillo
email: carrillo@imperial.ac.uk

F James

F Lagoutière

N Vauchelet
email: nicolas.vauchelet@upmc.fr

The Filippov characteristic flow for the aggregation equation with mildly singular potentials

Keywords: aggregation equation, nonlocal conservation equations, measure-valued solutions, gradient flow, Filippov's flow, finite volume schemes. 2010 AMS subject classifications: 35B40, 35D30, 35L60, 35Q92, 49K20

Existence and uniqueness of global in time measure solution for the multidimensional aggregation equation is analyzed. Such a system can be written as a continuity equation with a velocity field computed through a self-consistent interaction potential. In Carrillo et al. (Duke Math J (2011)) [16], a well-posedness theory based on the geometric approach of gradient flows in measure metric spaces has been developed for mildly singular potentials at the origin under the basic assumption of being λ-convex. We propose here an alternative method using classical tools from PDEs. We show the existence of a characteristic flow based on Filippov's theory of discontinuous dynamical systems such that the weak measure solution is the pushforward measure with this flow. Uniqueness is obtained thanks to a contraction argument in transport distances using the λ-convexity of the potential. Moreover, we show the equivalence of this solution with the gradient flow solution. Finally, we show the convergence of a numerical scheme for general measure solutions in this framework allowing for the simulation of solutions for initial smooth densities after their first blow-up time in L p -norms.

Introduction

This paper is devoted to the so-called aggregation equation in d space dimension

∂ t ρ = div (∇ x W * ρ)ρ , t > 0, x ∈ R d , (1.1)
complemented with the initial condition ρ(0, x) = ρ ini . Here, W plays the role of an interaction potential whose gradient ∇ x W (xy) measures the relative force exerted by an infinitesimal mass localized at a point y onto an infinitesimal mass located at a point x. This system appears in many applications in physics and population dynamics. In the framework of granular media, equation (1.1) is used to describe the large time dynamics of inhomogeneous kinetic models (see [START_REF] Benedetto | A kinetic equation for granular media[END_REF][START_REF] Carrillo | Contractions in the 2-Wasserstein length space and thermalization of granular media[END_REF][START_REF] Toscani | Kinetic and hydrodynamic models of nearly elastic granular flows[END_REF]). Model of crowd motion with a nonlinear dependancy of the term ∇ x W * ρ are also encountered in [START_REF] Colombo | A class of nonlocal models for pedestrian traffic[END_REF][START_REF] Crippa | Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow[END_REF]. In population dynamics, (1.1) provides a biologically meaningful description of aggregative phenomena. The description of the collective migration of cells by swarming leads to such non-local interaction PDEs (see e.g. [START_REF] Morale | An interacting particle system modelling aggregation behavior: from individuals to populations[END_REF][START_REF] Okubo | Diffusion and Ecological Problems: Modern Perspectives[END_REF][START_REF] Topaz | Swarming patterns in a two-dimensional kinematic model for biological groups[END_REF]). Another example is the modelling of bacterial chemotaxis. In this framework, the quantity S = W * ρ is the chemoattractant concentration which is a substance emitted by bacteria allowing them to interact with each others. The dynamics can be macroscopically modelled by the Patlak-Keller-Segel system [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF][START_REF] Patlak | Random walk with persistence and external bias[END_REF]. In the kinetic framework, the Othmer-Dunbar-Alt model is usually used, its hydrodynamic limit leads to the aggregation equation (1.1) [START_REF] Dolak | Kinetic models for chemotaxis: Hydrodynamic limits and spatiotemporal mechanisms[END_REF][START_REF] Filbet | Derivation of hyperbolic models for chemosensitive movement[END_REF][START_REF] James | Chemotaxis: from kinetic equations to aggregation dynamics[END_REF]. In many of these examples, the potential W is usually mildly singular, i.e. W has a weak singularity at the origin. Due to this weak regularity, finite time blow-up of regular solutions has been observed for such systems and has gained the attention of several authors (see e.g. [START_REF] Li | Long time asymptotics of kinetic models of granular flows[END_REF][START_REF] Bodnar | An integro-differential equation arising as a limit of individual cell-based models[END_REF][START_REF] Bertozzi | Blow-up in multidimensional aggregation equation with mildly singular interaction kernels[END_REF][START_REF] Bertozzi | L p theory for the multidimensional aggregation equation[END_REF][START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]). Finite time concentration is sometimes considered as a very simple mathematical way to mimick aggregation of individuals, as opposed to diffusion. Finally, attraction-repulsion potentials have been recently proposed as very simple models of pattern formation due to the rich structure of the set of stationary solutions, see [START_REF] Perthame | Traveling plateaus for a hyperbolic Keller-Segel system with attraction and repulsion: existence and branching instabilities[END_REF][START_REF] Brecht | On soccer balls and linearized inverse statistical mechanics[END_REF][START_REF] Brecht | Predicting pattern formation in particle interactions[END_REF][START_REF] Balagué | Dimensionality of local minimizers of the interaction energy[END_REF][START_REF] Bertozzi | Ring Patterns and their Bifurcations in a Nonlocal Model of Biological Swarms[END_REF] for instance.

Since finite time blow-up of regular solutions occurs, a natural framework to study the existence of global in time solutions is to work in the space of probability measures. However, several difficulties appear due to the weak regularity of the potential. In fact, the definition of the product of ∇W * ρ with ρ is a priori not well defined. This fact has already been noticed in one dimension in [START_REF] James | Chemotaxis: from kinetic equations to aggregation dynamics[END_REF][START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF]. Using defect measures in a two-dimensional framework, existence of weak measure solutions for parabolic-elliptic coupled system has been obtained in [START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF][START_REF] Dolbeault | The two-dimensional Keller-Segel model after blow-up[END_REF]. However, uniqueness is lacking. Measured valued solutions for the 2D Keller-Segel system have been considered in [START_REF] Luckhaus | Measure valued solutions of the 2D Keller-Segel system[END_REF] as limit of solutions of a regularized problem.

For the aggregation equation (1.1), a well-posedness theory for measure valued solutions has been considered using the geometrical approach of gradient flows in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]. This technique has been extended to the case with two species in [START_REF] Di Francesco | Measure solutions for non-local interaction PDEs with two species[END_REF]. The assumptions on the potential in order to get this well-posedness theory of measure valued solutions use certain convexity of the potential that allows for mild singularity of the potential at the origin.

In this paper, we assume that the interaction potential W : R d → R satisfies the following properties:

(A0) W is Lipschitz continuous, W (x) = W (-x) and W (0) = 0.

(A1) W is λ-convex for some λ ≤ 0, i.e. W (x) -λ 2 |x| 2 is convex.

(A2) W ∈ C 1 (R d \ {0}).
Such a potential will be named as a pointy potential. A typical example is a fully attractive Morse type potential, W (x) = 1e -|x| , which is -1-convex.

Let us emphasize that we only consider Lipschitz potentials which allows to bound the velocity field, whereas in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF], linearly growing at infinity potentials are allowed. In other words, we assume that there exists a nonnegative constant w ∞ such that for all x = 0,

|∇W (x)| ≤ w ∞ .
(1.

2)

The main reason for this restriction is to be able to work with suitable characteristics for this velocity field as explained below.

Denoting a = -∇W * ρ the macroscopic velocity, equation (1.1) can be considered as a conservative transport equation with velocity field a. Then a traditional definition for solutions is the one defined thanks to the characteristics corresponding to this macroscopic velocity. However, the velocity a is not Lipschitz and therefore we cannot defined classical solutions to the characteristics equation. To overcome this difficulty, Filippov [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Side[END_REF] has proposed a notion of solution which extend the classical one. Using this so-called Filippov flow, Poupaud & Rascle [START_REF] Poupaud | Measure solutions to the linear multidimensional transport equation with discontinuous coefficients[END_REF] have proposed a notion of solution to the conservative linear transport equation defined by X # ρ ini where X is the Filippov flow corresponding to the macroscopic velocity. However a stability result of the flow was still lacking until recently [START_REF] Bianchini | An estimate on the flow generated by monotone operators[END_REF], and thus there are no results with this technique for nonlinear equations of the form (1.1). We notice that in one dimension and for linear equations, these solutions are equivalent to the duality solutions defined in [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF][START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF], which have been successfully used in [START_REF] James | Chemotaxis: from kinetic equations to aggregation dynamics[END_REF][START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF] to tackle (1.1) in the one dimensional case.

On the other hand, although the geometric approach of gradient flows furnishes a wellposedness framework for well-posedness, this approach does not allow to define a characteristic flow corresponding to the macroscopic velocity a = -∇W * ρ. In this work, we focus on improving the understanding of these solutions by showing that under assumptions (A0)-(A2) on the potential, the solutions can be understood as their initial data pushed forward by suitable characteristic flows.

In order to achieve this goal, we first generalize the theory developed in [START_REF] Poupaud | Measure solutions to the linear multidimensional transport equation with discontinuous coefficients[END_REF] to the nonlinear aggregation equation (1.1). The first difficulty is, as it was in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF], to identify the right definition of the nonlinear term and the nonlinear product. This was solved in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF] by identifying the element of minimal norm by subdifferential calculus. We revisit this issue by clarifying that this is the right definition of the nonlinear term if we approximate a pointy potential by smooth symmetric potentials. Once the identification of the right velocity field has been done, we use the crucial stability results of Filippov's flows in [START_REF] Bianchini | An estimate on the flow generated by monotone operators[END_REF] to pass to the limit in the nonlinear terms. This leads to the construction of global measure solutions of the form X # ρ ini , where X is the Filippov flow associated to the velocity vector field a. This is the point where we need globally bounded velocity vector fields since Filippov's theory [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Side[END_REF] was only developed under these assumptions. In this way, we extend to the muti-dimensional case the results in [START_REF] James | Chemotaxis: from kinetic equations to aggregation dynamics[END_REF] (for a particular choice of the potential W ) and in [START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF].

Moreover, we are able to adapt arguments for uniqueness already used for the aggregation equation and for nonlinear continuity equations as in [START_REF] Loeper | Uniqueness of the solution to the Vlasov-Poisson system with bounded density[END_REF][START_REF] Carrillo | Uniqueness of bounded solutions to aggregation equations by optimal transport methods[END_REF][START_REF] Bertozzi | L p theory for the multidimensional aggregation equation[END_REF] to show the contraction property of the Wasserstein distance for our constructed solutions. This leads to a uniqueness result for our constructed solutions and to show the equivalence between the notion of gradient flow solutions and these Filippov's flow characteristics solutions. Let us further comment that in the one dimensional case it has been noticed that there is a link between solutions to (1.1) and entropy solutions to scalar conservation law for an antiderivative of ρ (see [START_REF] Bodnar | An integro-differential equation arising as a limit of individual cell-based models[END_REF][START_REF] Bonaschi | Equivalence of gradient flows and entropy solutions for singular nonlocal interaction equations in 1D[END_REF][START_REF] James | Chemotaxis: from kinetic equations to aggregation dynamics[END_REF][START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF]). This link has allowed to consider extensions of the model (1.1) with a nonlinear dependency of the term ∇W * ρ.

Finally, let us mention that apart from particle methods to the aggregation equations, very few numerical schemes have been proposed to simulate solutions of the aggregation equation after blow-up. The so-called sticky particle method was shown to be convergent in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF] and used to obtain qualitative properties of the solutions such as the finite time total collapse. However, this method is not that practical to deal with finite time blow-up and the behavior of solutions after blow-up in dimensions larger than one. In one dimension, such numerical simulations thanks to a particle scheme have been obtained by part of the authors in [START_REF] James | Chemotaxis: from kinetic equations to aggregation dynamics[END_REF]. Moreover, in the one dimensional case and with a nonlinear dependency of the term ∇W * ρ, they propose in [START_REF] James | Numerical method for one-dimensional aggregation equations[END_REF] a finite volume scheme allowing to simulate the behaviour after blow up and prove its convergence. Finally, extremely accurate numerical schemes have been developed to study the blow-up profile for smooth solutions, see [START_REF] Huang | Asymptotics of blowup solutions for the aggregation equation[END_REF][START_REF] Huang | Self-similar blowup solutions to an aggregation equation in R n[END_REF]. In fact, part of the authors recently proposed an energy decreasing finite volume method [START_REF] Carrillo | A Finite-Volume Method for Nonlinear Nonlocal Equations with a Gradient Flow Structure[END_REF] for a large class of PDEs including in particular (1.1) but no convergence result was given. Here, we give a convergence result for a finite volume scheme and for general measures as initial data. This allows for numerical simulations of solutions in dimension greater than one allowing to observe the behaviour after blow-up occurs.

The outline of the paper is the following. Next section is devoted to the definition of our notion of weak measure solutions for the aggregation equation. After introducing some notations, we first recall the basic results as obtained by Poupaud & Rascle [START_REF] Poupaud | Measure solutions to the linear multidimensional transport equation with discontinuous coefficients[END_REF] on measure solutions for conservative linear transport equations. Then we define the notion of solutions defined by a flow and state the main result of this paper in Theorem 2.5. Finally, we recall the existence result of gradient flow solutions in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF] and state their equivalence with solutions defined by a flow. Section 3 is devoted to the proof of the existence and uniqueness result. The main ingredient of the proof of existence is a one-sided Lipschitz property of the macroscopic velocity and an atomization strategy by approximating with finite Dirac Deltas. A contraction argument in Wasserstein distance for these solutions allows to recover the uniqueness. In Section 4, we investigate the numerical approximation of such solutions. A finite volume scheme is proposed and its convergence is established for general measure valued solutions. An illustration thanks to numerical simulations is also provided showing the ability of the scheme to capture the finite time total collapse and the qualitative interaction between different aggregates after the first blow-up in L p -norms. Finally, an Appendix is devoted to some technical Lemmas useful throughout the paper.

Weak measure solutions for the aggregation equation

All along the paper, we will make use of the following notations. We denote M loc (R d ) the space of local Borel measures on R d . For ρ ∈ M loc (R d ), we denote by |ρ|(R d ) its total variation. We denote M b (R d ) the space of measures in M loc (R d ) with finite total variation. From now on, the space of measures M b (R d ) is always endowed with the weak topology σ(M b , C 0 ). For T > 0, we denote S M := C([0, T ]; M b (R d )σ(M b , C 0 )). Finally, we define the space of probability measures with finite second order moment:

P 2 (R d ) = µ nonnegative Borel measure, µ(R d ) = 1, |x| 2 µ(dx) < ∞ .
This space is endowed with the Wasserstein distance d W defined by (see e.g. [START_REF] Villani | Optimal transport, old and new[END_REF][START_REF] Villani | Topics in optimal transportation[END_REF])

d W (µ, ν) = inf γ∈Γ(µ,ν) |y -x| 2 γ(dx, dy) 1/2 (2.1)
where Γ(µ, ν) is the set of measures on R d × R d with marginals µ and ν, i.e.

Γ(µ, ν) = γ ∈ P 2 (R d × R d ); ∀ ξ ∈ C 0 (R d ),
ξ(y 0 )γ(dy 1 , dy 2 ) = ξ(y 1 )µ(dy 1 ),

ξ(y 2 )γ(dy 1 , dy 2 ) = ξ(y 2 )ν(dy 2 ) .
From a minimization argument, we know that in the definition of d W the infimum is actually a minimum. A map that realizes the minimum in the definition (2.1) of d W is called an optimal plan, the set of which is denoted by Γ 0 (µ, ν). Then for all γ 0 ∈ Γ 0 (µ, ν), we have

d 2 W (µ, ν) = |y -x| 2 γ 0 (dx, dy).

Weak measure solutions for conservative transport equation

We recall in this Section some useful results on weak measure solutions to the conservative transport equation

∂ t u + div(bu) = 0; u(t = 0) = u 0 . (2.2)
We assume here that the vector field b is given. We start by the following definition of characteristics [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Side[END_REF] :

Definition 2.1 Let us assume that b = b(t, x) ∈ R d is a vector field defined on [0, T ] × R d with T > 0. A Filippov characteristic X(t; s, x) stems from x ∈ R d at time s is a continuous function X(•; s, x) ∈ C([0, T ], R d ) such that ∂ ∂t X(t; s, x) exists a.e. t ∈ [0, T ] satisfying ∂ ∂t X(t; s, x) ∈ Convess(b)(t, •) (X(t; s, x)) a.e. t ∈ [0, T ]; X(s; s, x) = x.
From now on, we will use the notation X(t, x) = X(t; 0, x).

In this definition Convess(E) denotes the essential convex hull of a set E. We remind the reader the definition for the sake of completeness, see [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Side[END_REF][START_REF] Aubin | Differential inclusions. Set-valued maps and viability theory[END_REF] for more details. We denote by Conv(E) the classical convex hull of E, i.e., the smallest closed convex set containing E. Given the vector field b(t, •) : R d -→ R d , the essential convex hull at point x is defined as

{Convess(b)(t, •)}(x) = r>0 N ∈N 0 Conv [b (t, B(x, r) \ N )] ,
where N 0 is the set of zero Lebesgue measure sets. Then, we have the following existence and uniqueness result of Filippov characteristics under the mere assumption that the vector field b is one-sided Lipschitz. R d )) satisfies the OSL condition, that is for all x and y in R d , for all t ∈ [0, T ],

(b(t, x) -b(t, y)) • (x -y) ≤ α(t) x -y 2 , for α ∈ L 1 (0, T ). (2.3)
Then there exists an unique Filippov characteristic X associated to this vector field.

An important consequence of this result is the existence and uniqueness of weak measure solutions for the conservative linear transport equation. This result has been proved by Poupaud and Rascle [START_REF] Poupaud | Measure solutions to the linear multidimensional transport equation with discontinuous coefficients[END_REF].

Theorem 2.3 ([42]) Let T > 0. Let b ∈ L 1 ([0, T ], L ∞ (R d
)) be a vector field satisfying the OSL condition (2.3). Then for any u 0 ∈ M b (R d ), there exists a unique measure solution u in S M to the conservative transport equation (2.2) such that u(t) = X(t) # u 0 , i.e. for any φ ∈ C 0 (R d ), we have

R d φ(x)u(t, dx) = R d φ(X(t, x))u 0 (dx), for t ∈ [0, T ].
Finally, we recall the following stability result for the Filippov characteristics which has been established by Bianchini

Solutions defined by Filippov's flow

We state in this Section the main result of this paper dealing with the existence and uniqueness of measure solutions defined thanks to the Filippov characteristics for the aggregation equation (1.1). For ρ ∈ C([0, T ], P 2 (R d )), we define the velocity field a ρ by a ρ (t, x) =y =x ∇W (xy)ρ(t, dy).

(2.4)

This choice of macroscopic velocity will be justified by the convergence result of Lemma 3.1 below. We remark that this definition of the velocity field coincides with the one based on subdifferential calculus done in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF], see next subsection. Due to the λ-convexity of W (A1), we deduce that for all x, y in R d \ {0} we have

(∇W (x) -∇W (y)) • (x -y) ≥ λ x -y 2 . (2.5)
For the sake of simplicity of the notations, we introduce

∇W (x) = ∇W (x), for x = 0; 0, for x = 0,
such that, by definition of the velocity (2.4), we have

a ρ (t, x) = - R d ∇W (x -y)ρ(t, dy) . (2.6)
Moreover, since W is even, ∇W is odd and by taking y = -x in (2.5), we deduce that inequality (2.5) is true even when x or y vanishes for ∇W :

∀ x, y ∈ R d , ( ∇W (x) -∇W (y)) • (x -y) ≥ λ x -y 2 .
(2.7)

We are now ready to state the main result of this paper. Its proof is postponed until Section 3 below. Theorem 2.5 Let W satisfy assumptions (A0)-(A2) and let ρ ini be given in P 2 (R d ). Given T > 0, there exists a unique Filippov characteristic flow X such that the pushforward measure ρ := X # ρ ini is a distributional solution of the aggregation equation

∂ t ρ + div( a ρ ρ) = 0, ρ(0, •) = ρ ini , (2.8)
where a ρ is defined by (2.4).

Besides, if ρ ini and µ ini are two given nonnegative measure in P 2 (R d ), then the corresponding pushforward measures ρ and µ satisfy for all t ∈ [0, T ]

d W (ρ(t), µ(t)) ≤ e -2λt d W (ρ ini , µ ini ).
(2.9)

Remark 2.6 Let us point out that the exponent in the stability estimate in d W in (2.9) can be improved to -λt if both initial measures ρ ini and µ ini have the same center of mass.

Gradient flow solutions

We recall the definition of gradient flow solutions as defined in [START_REF] Ambrosio | Gradient flows in metric space of probability measures[END_REF][START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]. Let W be the energy of the system defined by

W(ρ) = 1 2 R d ×R d W (x -y) ρ(dx)ρ(dy).
(2.10)

We say that µ ∈ AC 2 loc ([0, +∞);

P 2 (R d )) if µ is locally Hölder continuous of exponent 1/2 in time with respect to the distance d W in P 2 (R d ).
Definition 2.7 (Gradient flows) Let W satisfy assumptions (A0)-(A2). We say that a map µ ∈ AC 2 loc ([0, +∞); P 2 (R d )) is a solution of a gradient flow equation associated to the functional W, defined in (2.10), if there exists a Borel vector field

v such that v(t) ∈ T an µ(t) P 2 (R d ) for a.e. t > 0, v(t) L 2 (µ) ∈ L 2
loc (0, +∞), the continuity equation

∂ t µ + div vµ = 0,
holds in the sense of distributions, and v(t) = -∂ 0 W(µ(t)) for a.e. t > 0. Here ∂ 0 W(µ) denotes the element of minimal norm in ∂W(µ), which is the subdifferential of W at the point µ.

We refer to [START_REF] Ambrosio | Gradient flows in metric space of probability measures[END_REF][START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF] for details about the definition of the subdifferential since we will not make use of them in the sequel. The existence and uniqueness result of [16, Theorem 2.12 and 2.13 ] can now be synthetized as follows.

Theorem 2.8 ( [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]) Let W satisfy assumptions (A0)-(A2). Given ρ ini ∈ P 2 (R d ), there exists a unique gradient flow solution of (1.1), i.e. a curve

ρ GF ∈ AC 2 loc ([0, ∞); P 2 (R d )) satisfying ∂ρ GF (t) ∂t + div(v(t)ρ GF (t)) = 0, in D ′ ([0, ∞) × R d ), v(t, x) = -∂ 0 W(ρ GF )(t, x) = - y =x ∇W (x -y) ρ GF (t, dy),
with ρ GF (0) = ρ ini . Moreover, the following energy identity holds for all 0 ≤ t 0 ≤ t 1 < ∞:

t 1 t 0 R d |∂ 0 W * ρ GF | 2 ρ GF (t, dx)dt + W(ρ GF (t 1 )) = W(ρ GF (t 0 )).
Theorems 2.5 and 2.8 furnish two notions of solutions to (1.1) which are solutions in the sense of distributions. Then we should wonder on the link between this two notions. The following result states their equivalence. Theorem 2.9 Let W satisfy assumptions (A0)-(A2). Let ρ ini ∈ P 2 (R d ) be given. Let us denote ρ the solution of Theorem 2.5 and by ρ GF the solution of Theorem 2.8. Then we have ρ ∈ AC 2 loc ([0, ∞); P 2 (R d )) and ρ = ρ GF . As a consequence of this equivalence result, there exists a unique solution ρ which satisfies in the sense of distribution (2.8) with a ρ defined in (2.4). This solution is a pushforward measure by a characteristic flow: ρ = X # ρ ini .

Existence and uniqueness

Macroscopic velocity and one-sided estimate

In order to justify the choice of the expression of the macroscopic velocity in (2.4), we prove a stability result for symmetric potentials. Moreover, we state in Lemma 3.3 the important one-sided Lipschitz property for this macroscopic velocity. Lemma 3.1 Let us assume that W satisfies assumptions (A0)-(A2). Let (W n ) n∈N * be a sequence of even functions in C 1 (R d ) satisfying (A1) and (1.2) with the same constants λ and w ∞ not depending on n and such that

sup x∈R d \B(0, 1 n ) ∇W n (x) -∇W (x) ≤ 1 n , for all n ∈ N * . (3.1)
If the sequence ρ n ⇀ ρ weakly as measures, then for every continuous compactly supported φ, we have

lim n→+∞ R d ×R d φ(x)∇W n (x -y)ρ n (dx)ρ n (dy) = R d ×R d \D φ(x)∇W (x -y)ρ(dx)ρ(dy),
where

D is the diagonal in R d : D = {(x, x), x ∈ R d }.
Proof. The construction of such an approximating sequence of potentials can be obtained for instance using the Moreau-Yosida regularization, see [START_REF] Ambrosio | Gradient flows in metric space of probability measures[END_REF] and [START_REF] Carrillo | Gradient flows for non-smooth interaction potentials[END_REF]Proposition 3.5]. Let us focus on the last property. We first notice that by symmetry of W n , we have for all φ ∈ Lip(R d ),

R d φ(x)a n (x)ρ n (dx) = 1 2 R d ×R d (φ(x) -φ(y))∇W n (x -y)ρ n (dx)ρ n (dy).
We recall that since ρ n ⇀ ρ weakly as measures, we have that ρ n ⊗ ρ n ⇀ ρ ⊗ ρ weakly as measures. Let ε > 0. Since φ is continuous on a compact set, it is uniformly continuous therefore there exists α > 0 such that |φ(x)φ(y)| ≤ ε for |x -y| ≤ α. Then, defining

D α = {(x, y) ∈ R d × R d , |x -y| < α}
for any α > 0, we split the latter integral into :

R d ×R d (φ(x) -φ(y)) ∇W n (x -y)ρ n (dx)ρ n (dy) -∇W (x -y)ρ(dx)ρ(dy) = R d ×R d \Dα (φ(x) -φ(y)) ∇W n (x -y)ρ n (dx)ρ n (dy) -∇W (x -y)ρ(dx)ρ(dy) + Dα (φ(x) -φ(y)) ∇W n (x -y)ρ n (dx)ρ n (dy) -∇W (x -y)ρ(dx)ρ(dy) .
For the last term of the right hand side, we use the fact that φ is uniformly continuous and (1.2) for W and W n to prove that

Dα (φ(x) -φ(y)) ∇W n (x -y)ρ n (dx)ρ n (dy) -∇W (x -y)ρ(dx)ρ(dy) ≤ Cε.
For the first term, we have

R d ×R d \Dα (φ(x) -φ(y)) ∇W n (x -y)ρ n (dx)ρ n (dy) -∇W (x -y)ρ(dx)ρ(dy) = R d ×R d \Dα (φ(x) -φ(y)) ∇W n (x -y) -∇W (x -y) ρ n (dx)ρ n (dy) + R d ×R d \Dα (φ(x) -φ(y)) ∇W (x -y) ρ n (dx)ρ n (dy) -ρ(dx)ρ(dy) .
Using (3.1) we deduce that the first term of the right hand side is bounded by ε for n large enough. For the second term, we use the fact that (x, y) → (φ(x)φ(y)) ∇W (xy) is continuous and compactly supported and the tight convergence of ρ n towards ρ to prove it is bounded by ε when n is large enough. This concludes the proof. Remark 3.2 In other words, this Lemma states that if W n is an approximating smooth and even sequence for W and for any sequence ρ n converging to ρ in S M , then, denoting a n = ∇W n * ρ n , we have the convergence of the flux a n ρ n ⇀ a ρ ρ in the weak topology S M with a ρ defined in (2.4). A similar convergence result has been proved in [START_REF] Poupaud | Diagonal defect measures, adhesion dynamics and Euler equation[END_REF], although in this paper, the potential is less regular and in particular it does not satisfies (A0) neither the bound (1.2). Then at the limit the author recovers a defect measure which vanishes in our case. Such result has also been used in [START_REF] Dolbeault | The two-dimensional Keller-Segel model after blow-up[END_REF] to define weak solution for the two-dimensional Keller-Segel system for chemotaxis.

Lemma 3.3 Let ρ(t) ∈ M b (R d ) be nonnegative such that |ρ(t, •)|(R d ) ≤ c
for all t ≥ 0. Then under assumptions (A0) -(A2) the function (t, x) → a ρ (t, x) defined in (2.4) or equivalently in (2.6) satisfies the one-sided Lipschitz (OSL) estimate

( a ρ (t, x) -a ρ (t, y)) • (x -y) ≤ -λ|ρ|(R d ) x -y 2 . (3.2)
Proof. This result is an easy consequence of the λ-convexity of the potential. In fact, by definition (2.6), we have

a ρ (x) -a ρ (y) = - R d ∇W (x -z) -∇W (y -z) ρ(dz).
Using inequality (2.7) and the nonnegativity of ρ, we readily obtain (3.2).

From Lemma 3.3, we deduce that if ρ ∈ C([0, T ], P 2 (R)) and a ρ is defined as in (2.4), we can define the Filippov characteristic flow, denoted X, associated to the velocity field a ρ (see [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Side[END_REF]). Then we consider the push-forward measure

ρ P R := X # ρ ini .
Poupaud & Rascle [START_REF] Poupaud | Measure solutions to the linear multidimensional transport equation with discontinuous coefficients[END_REF] have shown that this measure is the unique measure solution of the conservative linear transport equation

∂ t ρ P R + div( a ρ ρ P R ) = 0.
The difficulty here is that the measure ρ used in the definition of the macroscopic velocity a ρ is a priori not the same as ρ P R . Actually, the whole aim of the next subsection is to prove that they are equal.

Existence

In this subsection, we prove the existence part of Theorem 2.5. We follow the idea of atomization consisting in approximating the solution by a finite sum of Dirac masses or particles, and then passing to the limit. This approach has been very successful for the aggregation equation, see [START_REF] Bertozzi | Blow-up in multidimensional aggregation equation with mildly singular interaction kernels[END_REF][START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF][START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF][START_REF] Bonaschi | Equivalence of gradient flows and entropy solutions for singular nonlocal interaction equations in 1D[END_REF] for instance.

Approximation with Dirac masses. Let us assume that the initial density is given by ρ ini,N (x) = N i=1 m i δ(xx 0 i ), with x 0 i = x 0 j for i = j, for a finite integer N and belongs to P 2 (R d ), i.e. we have

N i=1 m i = 1, M 2 (0) := N i=1 m i |x 0 i | 2 < +∞. (3.3)
Then we look for a solution of the aggregation equation given by

ρ N (t, x) = N i=1 m i δ(x -x i (t)).
By definition (2.4) we have

a ρ N (t, x) =          - N j=1 m j ∇W (x -x j (t)) , if x = x i , i = 1, . . . , N, - j =i m j ∇W (x i (t) -x j (t)) , otherwise.
For such a macroscopic velocity, we can define the Filippov characteristic X N as in Definition 2.1. In fact, from Lemma 3.3, a ρ N satisfies the OSL condition, which allows to define uniquely the Filippov characteristic. It is obvious from the essential convex hull definition that

- j =i m j ∇W (x i (t) -x j (t)) ∈ {Convess( a ρ N )(t, •)}(x i (t)).
Then setting the classical ODE system x ′ i (t) =j =i m j ∇W (x i (t)x j (t)), the solution will be defined up to the time t c of the first collision between two or more particles. By uniqueness of the Filippov characteristic, X N (t, x 0 i ) = x i (t) until that time. At time t c , one has to recompute the velocity field, since the colliding particles will stick together for later times according to the rule given by

∂ ∂t X N (t; s, x) ∈ {Convess( a ρ N )(t, •)}( X N (t; s, x)) a.e. t ∈ [0, T ]; X N (s; s, x) = x.
This construction of the characteristics coincides with the one done in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]Remark 2.10]. In other words, the Filippov flow coincides with this time evolution+collision+gluing of particles procedure.

Next, we define ρ N P R = X N # ρ ini,N . By construction, this measure satisfies in the sense of distributions ∂ t ρ N P R + div a ρ N ρ N P R = 0. Moreover, from the definition of the pushforward measure, we can write

a ρ N P R = - R d ∇W (x -y)ρ N P R (dy) = - R d ∇W (x -X N (t, y))ρ ini,N (dy).
By definition of ρ ini,N , we deduce

a ρ N P R (t, x) = - N i=1 m i R d ∇W (x -X N (t, y))δ(y -x 0 i ) = - N i=1 m i ∇W (x -X N (t, x 0 i )) = a ρ N (t, x).
Thus we conclude that ρ N P R = ρ N . Let us consider now the bound on the second moment. We define

M N 2 (t) := N i=1 m i |x i (t)| 2 . Differentiating, we have d dt M N 2 (t) = 2 N i=1 N j=1 m i m j x i ∇W (x i -x j ).
Using (1.2), we deduce that

d dt M N 2 (t) ≤ 2C N i=1 N j=1 m i m j |x i |.
From the Cauchy-Schwarz inequality and the fact that i m i = 1, we deduce

d dt M N 2 (t) ≤ K(1 + M N 2 (t)). (3.4) Since M N 2 (0) is finite from (3.
3), we deduce from a Gronwall Lemma that for all t ∈ [0, T ] we have M N 2 (t) < +∞. By continuity of the Filippov flow, we have that ρ N ∈ C([0, T ], P 2 (R d )). Moreover, using (1.2), we deduce that

| a ρ N (t, x)| ≤ C.
(3.5)

Passing to the limit N → +∞. Let us assume that ρ ini ∈ P 2 (R d ) and consider an approximation ρ ini,N ∈ P 2 (R d ) given by a finite sum of Dirac masses such that ρ ini,N → ρ ini tightly in M b (R) as N → +∞ with a uniform in N bound of the second moments, or equivalently, d W (ρ ini,N , ρ ini ) → 0 as N → ∞. We have proved above that we can construct a Filippov flow X N and a measure

ρ N = X N # ρ ini,N ∈ C([0, T ], P 2 (R d ))
such that in the distributional sense

∂ t ρ N + div( a ρ N ρ N ) = 0,
where a ρ N is defined by (2.4). From (3.5), we have that

a ρ N is bounded in L ∞ ([0, T ] × R d ).
Thus a ρ N converges up to a subsequence towards b in L ∞ t,xweak * . We can pass to the limit in the distributional sense in the one-sided Lipschitz inequality (3.2) satisfied by a ρ N , since the right hand side of this inequality does not depend on N . Then b satisfies the OSL condition and we can define X b the Filippov flow corresponding to b. From the L ∞ t,xweak * convergence above, it is obvious that a ρ N converges weakly to b in L 1 ([0, T ]; L 1 loc (R d )). Therefore, we can apply Theorem 2.4, and deduce that

X N → X b locally in C([0, T ] × R d ) as N → +∞.
Moreover, for every φ ∈ C 0 (R d ), we have

R d φ(x)ρ N (t, dx) = R d φ( X N (t, x))ρ ini,N (dx).
Since ρ ini,N ⇀ ρ ini tightly and

X N (t, x) → X b (t, x) locally in C([0, T ] × R d ), we deduce that for any R > 0, lim N →+∞ B(0,R) φ X N (t, x) ρ ini,N (dx) = B(0,R) φ(X b (t, x)) ρ ini (dx).
Denoting as above M N 2 (0) (resp. M 2 (0)) the second order moment of ρ ini,N (resp. ρ ini ), we infer that

R d \B(0,R) ρ ini,N (dx) ≤ M N 2 (0) R 2 ≤ d W (ρ ini,N , ρ ini ) + M 2 (0) R 2 .
This implies that for all φ ∈ C 0 (R d ),

R d φ( X N (t, x))ρ ini,N (dx) -→ N →+∞ R d φ(X b (t, x))ρ ini (dx) = R d φ(x) X b # ρ ini (dx).
We deduce that ρ N ⇀ ρ := X b # ρ ini in S M as N → +∞. Finally, from this latter convergence, we deduce by applying Lemma A.1 that a ρ N → a ρ a.e. By uniqueness of the limit, we conclude that b = a ρ a.e.

Bound on the second moment. Finally, we recover the bound in P 2 (R d ). We first notice that due to the approximation of the initial data done in the previous step, we know that M N 2 (0) is bounded uniformly in N . Taking into account this fact together with (3.4), there exists a nonnegative constant C T depending only on T and the initial data ρ ini such that

M N 2 (t) = R d | X N (t, x)| 2 ρ ini,N (dx) ≤ C T .
Then |x| 2 ρ N (t) is a bounded sequence of nonnegative measures that converges weakly as measures to |x| 2 ρ(t). Therefore, by the Banach-Alaoglu theorem, we get

M 2 (t) = R d |X b (t, x)| 2 ρ ini (dx) ≤ lim inf N →∞ M N 2 (t) ≤ C T .
This ends the proof of existence.

Uniqueness

The proof of the uniqueness relies on a contraction property with respect to the Wasserstein distance d W . In the framework of general gradient flows, this property has been established using the λ-geodesically convexity of the energy in [START_REF] Ambrosio | Gradient flows in metric space of probability measures[END_REF]Theorem 11.1.4], see also [START_REF] Carrillo | Contractions in the 2-Wasserstein length space and thermalization of granular media[END_REF][START_REF] Loeper | Uniqueness of the solution to the Vlasov-Poisson system with bounded density[END_REF][START_REF] Carrillo | Uniqueness of bounded solutions to aggregation equations by optimal transport methods[END_REF][START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF] for related results. We show here an equivalent result for our notion of solution. The proof relies strongly on the definition of the solution as a pushforward measure associated to a flow and of the λ-convexity of W .

Proposition 3.4 Let us assume that W satisfies assumptions (A0) -(A2). Let ρ 0 and ρ 0 be two nonnegative measure in P 2 (R d ). Let ρ and ρ in C([0, T ], P 2 (R d )) be solutions of the aggregation equation as in Theorem 2.5 with initial data ρ 0 and ρ 0 respectively. Then for all t > 0,

d W (ρ(t), ρ(t)) ≤ e -2λt d W (ρ 0 , ρ 0 ) .
Moreover if ρ 0 and ρ 0 have the same center of mass, then for all t > 0,

d W (ρ(t), ρ(t)) ≤ e -λt d W (ρ 0 , ρ 0 ) .
Proof. Let ρ 0 and ρ 0 be two nonnegative measure in P 2 (R d ). We first choose an optimal plan γ 0 ∈ Γ 0 (ρ 0 , ρ 0 ) such that we have

d 2 W (ρ 0 , ρ 0 ) = R d ×R d |x 1 -x 2 | 2 γ 0 (dx 1 , dx 2 ).
We regularize the potential W , as in Lemma 3.1, by

W ε ∈ C 1 (R d ) such that W ε is λ-convex, W ε (-x) = W ε (x), |∇W ε | ≤ |∇W | and sup x∈R d \B(0,ε) ∇W ε (x) -∇W (x)| ≤ ε.
As in subsection 3.2, we construct a Filippov flow X ε associated to the velocity field

a ε := -R d ∇W ε (x -y)ρ ε (t, dy) such that ρ ε = X ε # ρ 0 ∈ C([0, T ], P 2 (R d )
) is a measure solution to the aggregation equation ∂ t ρ ε + div(a ε ρ ε ) = 0 with initial data ρ 0 . For this flow we have

d dt X ε (t, x) = - R d ∇W ε (x -y)ρ ε (t, dy); X ε (0, x) = x.
Similarly we construct

ρ ε = X ε # ρ 0 ∈ C([0, T ], P 2 (R d )) associated to the velocity field a ε := -R d ∇W ε (x -y) ρ ε (t, dy).
By definition of the pushforward measure, we have that

a ε (t, x) = - R d ∇W ε (x -X ε (t, y))ρ 0 (dy), a ε (t, x) = - R d ∇W ε (x -X ε (t, y)) ρ 0 (dy).
Moreover, from the definition of the optimal plan γ 0 we can rewrite

a ε (t, x) = - R d ×R d ∇W ε (x -X ε (t, y 1 )) γ 0 (dy 1 , dy 2 ), (3.6)
a ε (t, x) = - R d ×R d ∇W ε (x -X ε (t, y 2 )) γ 0 (dy 1 , dy 2 ). (3.7)
Since ρ ε (t) belongs to P 2 (R d ), we have that

R d |x| 2 ρ ε (t, dx) = R d |X ε (t, x)| 2 dρ 0 (x) < ∞.
The same estimate holds true for ρ ε . Then we can consider the quantity

I ε (t) = R d ×R d X ε (t, x 1 ) -X ε (t, x 2 ) 2 γ 0 (dx 1 , dx 2 ).
We notice that for t = 0, we have I(0) = d 2 W (ρ 0 , ρ 0 ). We have

d dt I ε = 2 R d ×R d a ε (t, X ε (t, x 1 )) -a ε (t, X ε (t, x 2 )) • (X ε (t, x 1 ) -X ε (t, x 2 )) γ 0 (dx 1 , dx 2 ).
From the definition of the velocity field (3.6)-(3.7), we have

d dt I ε = -2 (R d ) 4 ∇W ε (X ε (t, x 1 ) -X ε (t, y 1 )) -∇W ε ( X ε (t, x 2 ) -X ε (t, y 2 )) • (X ε (t, x 1 ) -X ε (t, x 2 )) γ 0 (dx 1 , dx 2 )γ 0 (dy 1 , dy 2 ).
From assumption W ε (-x) = W ε (x), we deduce that ∇W ε is odd. Then ∇W ε (X ε (t, x) -X ε (t, y)) = -∇W ε (X ε (t, y) -X ε (t, x)) for all x, y. By exchanging the role of (x 1 , x 2 ) and (y 1 , y 2 ) in this latter equality and using the symmetry of ∇W ε we deduce that

d dt I ε = 2 (R d ) 4 ∇W ε (X ε (t, x 1 ) -X ε (t, y 1 )) -∇W ε ( X ε (t, x 2 ) -X ε (t, y 2 )) • (X ε (t, y 1 ) -X ε (t, y 2 )) γ 0 (dx 1 , dx 2 )γ 0 (dy 1 , dy 2 ).
Summing these two latter equalities, we obtain

d dt I ε = - (R d ) 4 ∇W ε (X ε (t, x 1 ) -X ε (t, y 1 )) -∇W ε ( X ε (t, x 2 ) -X ε (t, y 2 )) • (X ε (t, x 1 ) -X ε (t, y 1 ) -X ε (t, x 2 ) + X ε (t, y 2 )) γ 0 (dx 1 , dx 2 )γ 0 (dy 1 , dy 2 ).
From the λ-convexity of W , we deduce from (2.5) that

d dt I ε ≤ -λ (R d ) 4 X ε (t, x 1 ) -X ε (t, y 1 ) -X ε (t, x 2 ) + X ε (t, y 2 ) 2 γ 0 (dx 1 , dx 2 )γ 0 (dy 1 , dy 2 ). (3.8) We recall that λ ≤ 0 and |ρ 0 |(R d ) = | ρ 0 |(R d ) = 1. A direct Young inequality leads to d dt I ε ≤ -4λI ε . (3.9)
Applying the Gronwall lemma, we deduce that

I ε (t) ≤ e -4λt I(0) = e -4λt d 2 W (ρ 0 , ρ 0 ). (3.10)
If the initial data have the same center of mass, then it is easy to check that the center of mass remains the same for both solutions for all times, that is, for all t ≥ 0

M 1 = R d x ρ ε (t, dx) = R d X ε (t, x) ρ 0 (dx) = (R d ) 2 X ε (t, x 1 ) γ 0 (dx 1 , dx 2 ) = R d y ρ ε (t, dy) = R d X ε (t, y) ρ 0 (dy) = (R d ) 2 X ε (t, x 2 ) γ 0 (dx 1 , dx 2 ) .
Thus, one can check that

(R d ) 4 (X ε (t, x 1 ) -X ε (t, x 2 )) • (X ε (t, y 1 ) -X ε (t, y 2 )) γ 0 (dx 1 , dx 2 )γ 0 (dy 1 , dy 2 ) = 0 .
Now, expanding the square in (3.8), we improve the decay by a factor of 2 in (3.9) getting

d dt I ε ≤ -2λI ε .
From now on, we stick to the general case to pass to the limit ε → 0 in (3.10). Since ρ ε (t) is bounded in P 2 (R d ) independently on ε, we deduce from the Prokhorov theorem that we can extract a subsequence such that ρ ε (t) ⇀ ρ(t) tightly. Then applying Lemma A.2 in the Appendix we deduce that a ε → a ρ for a.e. t ∈ [0, T ], x ∈ R d , where a ρ is defined in (2.4). Then we have shown that we can construct a Filippov characteristic flow X associated to the velocity field a ρ . Applying the stability result of [START_REF] Bianchini | An estimate on the flow generated by monotone operators[END_REF], recalled in Theorem 2.4, we deduce that X ε → X locally in C([0, T ] × R d ). We can proceed analogously for ρ(t), and thus, for any R > 0 we have lim ε→0 B(0,R)

|X(t, x) -X ε (t, x)| 2 ρ 0 (dx) = lim ε→0 B(0,R) | X(t, x) -X ε (t, x)| 2 ρ 0 (dx) = 0 .
We conclude that

B(0,R)×B(0,R) X ε (t, x 1 ) -X ε (t, x 2 ) 2 -X(t, x 1 ) -X(t, x 2 ) 2 γ 0 (dx 1 , dx 2 ) -→ 0 (3.11)
as ε → 0. Now, using (3.10) together with (3.11), we deduce

B(0,R)×B(0,R) X(t, x 1 ) -X(t, x 2 ) 2 γ 0 (dx 1 , dx 2 ) ≤ e -4λt d 2 W (ρ 0 , ρ 0 ) ,
for all R > 0, leading to our final desired estimate

I(t) := R d ×R d |X(t, x 1 ) -X(t, x 2 )| 2 γ 0 (dx 1 , dx 2 ) ≤ e -4λt d 2 W (ρ 0 , ρ 0 ). (3.12)
Finally, by definition of the Wasserstein distance (2.1), we deduce d 2 W (ρ, ρ) ≤ I(t) and the contraction inequality (2.9) follows directly.

The uniqueness of solution in Theorem 2.5 is then a trivial consequence of this contraction property. In fact, applying Proposition 3.4 for two solutions ρ and ρ with the same initial data ρ ini , we deduce from (3.12) that X = X on supp(ρ ini ) which implies that ρ = ρ.

Equivalence with gradient flow solutions

This subsection is devoted to the proof of the equivalence of solution defined by the Filippov flow with the gradient flow solution as stated in Theorem 2.9. For ρ ini given in P 2 (R d ), we denote ρ the solution of Theorem 2.5 and ρ GF the solution of Theorem 2.8. We have proved above the existence of a Filippov characteristic flow X such that ρ = X # ρ ini and ρ satisfies in the sense of distributions

∂ t ρ + div( a ρ ρ) = 0.
From the bound on a ρ in (3.5), we deduce since ρ belongs to C([0, T ],

P 2 (R d )) that a ρ is bounded in L 2 ([0, T ], L 2 (ρ(t))). Thus using Theorem 8.3.1 of [1], we deduce that ρ ∈ AC 2 ([0, T ], P 2 (R d )).
We can conclude that ρ is a gradient flow solution, see [START_REF] Ambrosio | Gradient flows in metric space of probability measures[END_REF]Sections 8.3 and 8.4] and [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]. We conclude the proof using the uniqueness of gradient flow solutions. As a consequence, the solutions constructed in Theorem 2.5 satisfy the energy identity, for all 0 ≤ t 0 ≤ t 1 < ∞,

t 1 t 0 R d | a ρ (t, x)| 2 ρ(t, dx) + W(ρ(t 1 )) = W(ρ(t 0 )).

Numerical approximation

This Section is devoted to the convergence of a numerical scheme for simulating solutions given by Theorem 2.5. The theory of existence developed in the previous section will allow to prove convergence of standard finite volume schemes, provided the discretized macroscopic velocity is accurately defined. Before that, we would like to comment on particle schemes.

Particle scheme

The contraction estimate in d W for solutions leads to a theoretical estimate of the convergence error of the particle scheme used in the first step of the proof of Theorem 2.5. This was already pointed out in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF] in the framework of gradient flow solutions and used for qualitative behavior properties. We just remind the main result here for completeness. Let us consider an initial distribution given by a finite sum of N Dirac masses ρ ini,N = N i=1 m i δ(xx 0 i ). We consider the sticky particles dynamics given by

x ′ i (t) = - j =i m j ∇W (x i (t) -x j (t)), x i (0) = x 0 i , i = 1, . . . , N.
These dynamics are well defined provided x i (t) = x j (t). When two or more particles meet, we stick them and the resulting system follows the same dynamics with one or more particle less. This system of ODEs plus the collision+gluing particle procedure gives the solution ρ N (t) = N i=1 m i δ(xx i (t)) of Theorem 2.5 at time t ≥ 0 with initial data ρ ini,N as explained in the first step of its proof. i ) an approximation such that d W (ρ ini , ρ ini,N ) → 0 as N → +∞. Given T > 0, then the corresponding solution ρ N with initial data ρ ini,N defined above verifies

sup t∈[0,T ] d W (ρ(t), ρ N (t)) -→ N →+∞ 0.
The previous corollary is a direct consequence of the stability property in Theorem 2.5. Although this result is very nice from the theoretical viewpoint, it is not that useful for simulating the evolution of equation (1.1) for fully attractive potentials in practice. The reason is twofold. On one hand, to get a good control on the error after a long time one needs a very large number of particles. On the other hand, the treatment of the collision between particles and the gluing procedure is not too difficult in one dimension but it is very cumbersome (and difficult to control its error) in more dimensions. Nevertheless, particle simulations lead to a very good understanding of qualitative properties of solutions for attractive-repulsive potentials where collisions do not happen, see [START_REF] Brecht | On soccer balls and linearized inverse statistical mechanics[END_REF][START_REF] Brecht | Predicting pattern formation in particle interactions[END_REF][START_REF] Balagué | Dimensionality of local minimizers of the interaction energy[END_REF][START_REF] Bertozzi | Ring Patterns and their Bifurcations in a Nonlocal Model of Biological Swarms[END_REF] for instance. We finally mention the recent result of convergence of smooth particle schemes toward smooth solutions of the aggregation equation before blow-up in [START_REF] Craig | A blob method for the aggregation equation[END_REF].

Finite volume discretization

In the next three subsections, we will concentrate on the convergence of a finite volume scheme for the solutions constructed in Theorem 2.5 with general measures as initial data. The one dimensional case has been considered in [START_REF] James | Numerical method for one-dimensional aggregation equations[END_REF]. This case is particular since we can define an antiderivative of the measure solution ρ which is then a BV function solution of an equation obtained by integrating the aggregation equation. This fact is very much connected to the relation of the one dimensional case with conservation laws as in [START_REF] Bodnar | An integro-differential equation arising as a limit of individual cell-based models[END_REF][START_REF] James | Chemotaxis: from kinetic equations to aggregation dynamics[END_REF][START_REF] Bonaschi | Equivalence of gradient flows and entropy solutions for singular nonlocal interaction equations in 1D[END_REF][START_REF] James | Equivalence between duality and gradient flow solutions for one-dimensional aggregation equations[END_REF]. Then the convergence of the numerical scheme relies on a TVD property. We refer the reader to [START_REF] James | Numerical method for one-dimensional aggregation equations[END_REF] for more details such as the importance of a good choice of the macroscopic velocity which is emphasized with some numerical examples.

We focus in this work to higher dimensions where such techniques cannot be applied. For the sake of clarity, we restrict ourselves to the case d = 2. We consider a cartesian grid x i = i∆x and y j = j∆y, for i ∈ Z and j ∈ Z. We denote by C ij the cells C ij = [x i , x i+1 ) × [y j , y j+1 ). The time discretization is given by t n = n∆t, n ∈ N. As usual, we denote ρ n ij an approximation of ρ(t n , x i , y j ). We consider that the potential W is given and satisfies assumptions (A0)-(A2).

Following the idea in [START_REF] James | Numerical method for one-dimensional aggregation equations[END_REF], we propose the following discretization. For a given nonnegative measure ρ ini ∈ P 2 (R 2 ), we define for i, j ∈ Z 2 ,

ρ 0 ij = 1 ∆x∆y C ij ρ ini (dx, dy) ≥ 0. ( 4.1)
Since ρ ini is a probability measure, the total mass of the system is i,j ρ 0 ij ∆x∆y = 1. Assuming that an approximating sequence (ρ n ij ) i,j is known at time n, then we compute the approximation at time t n+1 by :

ρ n+1 ij = ρ n ij - ∆t ∆x a x n i+1/2j ρ n i+1/2j -a x n i-1/2j ρ n i-1/2j - ∆t ∆y a y n ij+1/2 ρ n ij+1/2 -a y n ij-1/2 ρ n ij-1/2 + ∆t 2∆x w ∞ ρ n i+1j -2ρ n ij + ρ n i-1j + ∆t 2∆y w ∞ ρ n ij+1 -2ρ n ij + ρ n ij-1 , (4.2
) where w ∞ is defined in (1.2). We have used the notation

ρ i+1/2j = ρ ij + ρ i+1j 2 , ρ ij+1/2 = ρ ij + ρ ij+1 2 , a xi+1/2j = a xij + a xi+1j 2 , a y ij+1/2 = a y ij + a y ij+1 2 .
The macroscopic velocity is defined by

a xij = 1 ∆x∆y k,ℓ ρ kℓ D x W kℓ ij , a y ij = 1 ∆x∆y k,ℓ ρ kℓ D y W kℓ ij , (4.3)
where

D x W kℓ ij := C kℓ C ij ∂ x W x -x ′ , y -y ′ dxdy dx ′ dy ′ , D y W kℓ ij := C kℓ C ij ∂ y W x -x ′ , y -y ′ dxdy dx ′ dy ′ .
We notice after a straightforward change of variable that we have also

a xi+1/2j = 1 ∆x∆y k,ℓ ρ k+1/2ℓ D x W kℓ ij , a y ij+1/2 = 1 ∆x∆y k,ℓ ρ kℓ+1/2 D y W kℓ ij . (4.4)
Let us finally remark that this scheme is close to the Lax-Friedrichs flux formula for conservation laws. Therefore, it introduces some numerical viscosity in the simulations. This will be clear in the error terms obtained in the convergence proof since we will have error estimates depending on second order derivatives, see subsection 4.4.

Properties of the scheme

The following Lemma states a CFL-like condition for the scheme : Lemma 4.2 Let us assume that W satisfies (A0)-(A2) and consider ρ ini ∈ P 2 (R 2 ). We define ρ 0 ij by (4.1). Let us assume that the condition

w ∞ 1 ∆x + 1 ∆y ∆t ≤ 1 2 , (4.5)
is satisfied. Then the sequences computed thanks to the scheme defined in (4.2)-(4.3) satisfy for all i, j and n,

ρ n ij ≥ 0, |a x n ij | ≤ w ∞ , |a y n ij | ≤ w ∞ . Proof.
The total initial mass of the system is ∆x∆y i,j ρ 0 ij = 1. Since the scheme (4.2) is conservative, we have for all n ∈ N, ∆x∆y i,j ρ n ij = 1. We can rewrite equation (4.2) as

ρ n+1 ij = ρ n ij 1 - ∆t ∆x a x n i+1/2j -a x n i-1/2j 2 - ∆t ∆y a y n ij+1/2 -a y n ij-1/2 2 - ∆t ∆x w ∞ - ∆t ∆y w ∞ + ρ n i+1j ∆t 2∆x w ∞ -a x n i+1/2j + ρ n i-1j ∆t 2∆x w ∞ + a x n i-1/2j + ρ n ij+1 ∆t 2∆y w ∞ -a y n ij+1/2 + ρ n ij-1 ∆t 2∆y w ∞ + a y n ij-1/2 . (4.6)
Let us prove by induction on n that for all i, j, n we have ρ n ij ≥ 0. Let us assume that for a given n ∈ N we have ρ n ij ≥ 0 for all i, j. Then, from definition (4.3) and assumption (1.2) we clearly have that

|a x n ij | ≤ w ∞ ∆x∆y i,j ρ n ij = w ∞ ; |a y n ij | ≤ w ∞ .
Then assuming that the condition (4.5) holds, we deduce that in the scheme (4.6) all the coefficients in front of ρ n ij , ρ n i-1j , ρ n i+1j , ρ n ij-1 , and ρ n ij+1 are nonnegative. Thus, using the induction assumption, we deduce that ρ n+1 ij ≥ 0 for all i, j.

In the following Lemma, we gather some properties of the scheme: mass conservation, center of mass conservation and finite second order moment. Lemma 4.3 Let us assume that W satisfies (A0)-(A2) and consider ρ 0 ij defined by (4.1) for some ρ ini ∈ P 2 (R 2 ). Let us assume that (4.5) is satisfied. Then the sequence (ρ n ij ) constructed thanks to the numerical scheme (4.2)-(4.3) satisfies:

(i) Mass conservation and conservation of the center of mass: for all n ∈ N * , we have

i,j∈Z 2 ρ n ij ∆x∆y = i,j∈Z 2 ρ 0 ij ∆x∆y = 1 , i,j∈Z 2 x i ρ n ij = i,j∈Z 2 x i ρ 0 ij , i,j∈Z 2 y j ρ n ij = i,j∈Z 2 y j ρ 0 ij .
(ii) Bound on the second moment: there exists a constant C > 0 such that for all n ∈ N * , we have

M n 2 := i,j∈Z 2 (x 2 i + y 2 j )ρ n ij ∆x∆y ≤ e Ctn M 0 2 + 1 -1, (4.7)
where we recall that t n = n∆t.

Proof. We first notice that due to Lemma 4.2, we have that for all n, i, j the sequence (ρ n ij ) is nonnegative.

(i) The mass conservation is directly obtained by summing over i and j equation (4.2). For the center of mass, we have from (4.2) after using a discrete integration by parts :

i,j∈Z 2 x i ρ n+1 ij = i,j∈Z 2 x i ρ n ij - ∆t ∆x i,j∈Z 2 a x n i+1/2j ρ n i+1/2j x i -x i+1 + ∆t 2∆x w ∞ i,j∈Z 2 ρ n ij x i-1 -2x i + x i+1 .
From the definition x i = i∆x, we deduce i,j∈Z 2

x i ρ n+1 ij = i,j∈Z 2 x i ρ n ij -∆t i,j∈Z 2 a x n i+1/2j ρ n i+1/2j .
By definition of the macroscopic velocity (4.4), we have

i,j∈Z 2 a x n i+1/2j ρ n i+1/2j = 1 ∆x∆y i,j k,ℓ D x W kℓ ij ρ n k+1/2ℓ ρ n i+1/2j . Since the function ∂ x W is odd, we deduce that D x W kℓ ij = -D x W ij kℓ .
Then by exchanging the role of i, j and k, ℓ is the latter sum, we deduce that it vanishes. Thus, i,j∈Z 2

x i ρ n+1 ij = i,j∈Z 2
x i ρ n ij and we proceed in the same way with y j instead of x i .

(ii) For the second moment, still using (4.2) and a discrete integration by parts, we get i,j∈Z 2

x 2 i ρ n+1 ij = i,j∈Z 2 x 2 i ρ n ij - ∆t ∆x i,j∈Z 2 a x n i+1/2j ρ n i+1/2j x 2 i -x 2 i+1 + ∆t 2∆x w ∞ i,j∈Z 2 ρ n ij x 2 i-1 -2x 2 i + x 2 i+1 .
By definition x i = i∆x, we have (

x 2 i -x 2 i+1 ) = -2x i+1/2 ∆x and (x 2 i-1 -2x 2 i + x 2 i+1 ) = 2∆x 2 . Thus, i,j∈Z 2 x 2 i ρ n+1 ij = i,j∈Z 2 x 2 i ρ n ij + 2∆t i,j∈Z 2 a x n i+1/2j ρ n i+1/2j x i+1/2 + w ∞ ∆t∆x,
where we have used the conservation of the mass. From Lemma 4.2, we deduce that |a x n i+1/2j | ≤ w ∞ . Thus, after applying a Cauchy-Schwarz inequality and using the mass conservation, we get

i,j∈Z 2 a x n i+1/2j ρ n i+1/2j x i+1/2 ∆x∆y ≤ w ∞ 2 1 + i,j∈Z 2
x 2 i+1/2 ρ n i+1/2j ∆x∆y .

We deduce then that there exists a nonnegative constant C such that i,j∈Z 2

x 2 i ρ n+1 ij ∆x∆y ≤ 1 + C∆t i,j∈Z 2 x 2 i ρ n ij ∆x∆y + C∆t.
Doing the same with the term i,j∈Z 2 y 2 j ρ n+1 ij , we deduce that there exists a nonnegative constant C such that M n+1 2 ≤ 1 + C∆t M n 2 + C∆t. We conclude the proof using a discrete Gronwall Lemma.

Convergence of the numerical approximation

Let us denote by ∆ = max{∆x, ∆y}. We define the reconstruction

ρ ∆ (t, x, y) = n∈N i∈Z j∈Z ρ n ij 1 [n∆t,(n+1)∆t)×C ij (t, x, y), (4.8)
Therefore, we have by definition of

a n ij = (a x n ij , a y n ij ) in (4.3) that a n ij = 1 ∆x∆y C ij ∇W * ρ ∆ (t n , x, y) dxdy.
In the same manner, we define

a ∆ (t, x, y) = n∈N i∈Z j∈Z a n ij 1 [n∆t,(n+1)∆t)×C ij (t, x, y).
Then we have the following convergence result:

Theorem 4.4 Let us assume that W satisfies (A0)-(A2) and consider ρ ini ∈ P 2 (R 2 ). We define ρ 0 ij by (4.1). Let T > 0 be fixed. Then, if (4.5) is satisfied, the discretization ρ ∆ converges weakly in M b ([0, T ] × R 2 ) towards the solution ρ of Theorem 2.5 as ∆ := max{∆x, ∆y} goes to 0 with ∆t satisfying (4.5).

Proof. From Lemma 4.2, we have that ρ n ij ≥ 0 provided the condition (4.5) is satisfied. Moreover, by conservation of the mass we deduce that the sequence nonnegative bounded measures (ρ ∆ ) ∆ satisfies for all t ∈ [0, T ], |ρ ∆ (t)|(R 2 ) = 1. Therefore, we can extract a subsequence, still denoted (ρ ∆ ) ∆ , converging for the weak topology towards ρ as ∆t, ∆x and ∆y go to 0 satisfying (4.5), i.e. We choose ∆t > 0 and N T ∈ N * such that condition (4.5) holds and T = ∆tN T . Let φ ∈ D([0, T ] × R 2 ) be smooth and compactly supported. We denote

∀ φ ∈ C 0 ([0, T ] × R 2 ), T 0 R 2 φ(t, x, y)ρ ∆ (t, x, y) dxdydt -→ T 0 R 2 φ(t,
ψ n i,j = t n+1 tn C ij φ(t, x, y) dtdxdy, such that T 0 R 2 ρ ∆ (t, x, y)φ(t, x, y) dtdxdy = N T n=0 i∈Z j∈Z ρ n ij ψ n i,j .
In particular, we have n,i,j

1 ∆t ρ ∆ (t n+1 , x i , y j ) -ρ ∆ (t n , x i , y j ) ψ n i,j = - n,i,j ρ n i,j ψ n i,j -ψ n-1 i,j ∆t = - T 0 R 2
ρ ∆ (t, x, y) φ(t, x, y)φ(t -∆t, x, y) ∆t dtdxdy.

We have φ(t, x, y)φ(t -∆t, x, y) = ∂ t φ(t, x, y)∆t + O(∆t 2 ). From the weak convergence of ρ ∆ and the fact that ρ ∆ is a bounded measure with a bound not depending on the mesh, we deduce that the latter integral converges to

- T 0 R 2 ∂ t φ(t, x, y)ρ(t, dx, dy) dt.
By the same token, we have

n,i,j 1 2∆x ρ ∆ (t n , x i+1 , y j ) -2ρ ∆ (t n , x i , y j ) + ρ ∆ (t n , x i-1 , y j ) ψ n i,j = T 0 R 2 ρ ∆ (t, x, y) φ(t, x + ∆x, y) -2φ(t, x, y) -φ(t, x -∆x, y) 2∆x dtdxdy,
Using the fact that |φ(t, x + ∆x, y) -2φ(t, x, y)φ(t, x -∆x, y)| ≤ ∂ xx φ ∞ ∆x 2 , we deduce that this latter integral converges towards 0 as ∆t, ∆x and ∆y go to 0. Futhermore, we have

n,i,j 1 ∆x a x n i+1/2j ρ n i+1/2j -a x n i-1/2j ρ n i-1/2j ψ n i,j = = - 1 4∆x n,i,j a x n ij ρ n ij ψ n i+1,j -ψ n i-1,j + a x n i+1j ρ n ij ψ n i+1,j -ψ n i,j + a x n i-1j ρ n ij ψ n i,j -ψ n i-1,j = - 1 4∆x T 0 R 2 a x∆ (t,
R 2 a x∆ ρ ∆ ξ(x 1 , y 1 ) dx 1 dy 1 = = i,j,k,ℓ 1 ∆x∆y C ij C kℓ ρ n kℓ ρ n ij ∂ x W (x -x ′ , y -y ′ ) dxdydx ′ dy ′ C ij ξ(x 1 , y 1 ) dx 1 dy 1 ,
on the other hand,

R 2 ∂ x W * ρ ∆ ρ ∆ ξ(x, y) dxdy = i,j,k,ℓ C ij C kℓ ρ n kℓ ρ n ij ∂ x W (x -x ′ , y -y ′ )ξ(x, y) dx ′ dy ′ dxdy.
Moreover, for any test function ξ smooth and compactly supported, we have for all x, y ∈ C ij , 1 ∆x∆y

C ij ξ(x 1 , y 1 ) dx 1 dy 1 = ξ(x, y) + O(∆x) + O(∆y). Thus we have R 2 a x∆ ρ ∆ ξ(x, y) dxdy = R 2 ∂ x W * ρ ∆ ρ ∆ ξ(x, y) dxdy + O(∆x) + O(∆y).
Finally, we deduce from (4.10)

n,i,j 1 ∆x a x n i+1/2j ρ n i+1/2j -a x n i-1/2j ρ n i-1/2j ψ n i,j = - 1 2 I 1 + I 2 + O(∆x) + O(∆y), (4.11)
where

I 1 = T 0 R 4 ∂ x W (x -x ′ , y -y ′ )ρ ∆ (t, x ′ , y ′ )ρ ∆ (t, x, y)∂ x φ(t, x, y) dxdydt, I 2 = 1 2 T 0 R 4 ∂ x W (x + ∆x -x ′ , y -y ′ ) + ∂ x W (x -∆x -x ′ , y -y ′ ) ρ ∆ (t, x ′ , y ′ )ρ ∆ (t, x, y)∂ x φ(t, x, y) dx ′ dy ′ dxdydt.
As a direct consequence of Lemma 3.1, we have

I 1 -→ ∆→0 T 0 R 2 ∂ x W * ρ(t, x, y)ρ(t, x, y)∂ x φ(t, x, y) dtdxdy.
For the term I 2 , we proceed as in the proof of Lemma 3.1. We recall the main ingredients of this proof. First, using the symmetry of W , we write

I 2 = 1 4 T 0 R 4 ∂ x W (x + ∆x -x ′ , y -y ′ ) + ∂ x W (x -∆x -x ′ , y -y ′ ) ρ ∆ (t, x ′ , y ′ )ρ ∆ (t, x, y) ∂ x φ(t, x, y) -∂ x φ(t, x ′ , y ′ ) dx ′ dy ′ dxdydt.
We introduce the set D α = {(x, y, x ′ , y ′ ) s.t. |xx ′ | + |yy ′ | < α} for some positive coefficient α < ∆x and split the latter integral into the sum of the integral over R 4 \ D α and over D α .

Using the uniform continuity of ∂ x φ, we deduce that the integral over D α is small for small α.

Then, using the fact that by continuity of ∂ x W on R 4 \{0}, we have for all (x, x ′ , y, y

′ ) ∈ R 4 \D α lim ∆x→0 ∂ x W (x + ∆x -x ′ , y -y ′ ) + ∂ x W (x -∆x -x ′ , y -y ′ ) = ∂ x W (x -x ′ , y -y ′ ),
we deduce

I 2 -→ ∆→0 T 0 R 2 ∂ x W * ρ(t,
x, y)ρ(t, x, y)∂ x φ(t, x, y) dtdxdy.

Therefore, we conclude from (4.11)

lim ∆→0 n,i,j 1 ∆x a x n i+1/2j ρ n i+1/2j -a x n i-1/2j ρ n i-1/2j ψ n i,j = = - T 0 R 2 ∂ x φ(t, x, y) ∂ x W * ρ(t, x, y)ρ(t, x, y) dtdxdy.
Finally, multiplying equation (4.2) by ψ n i,j , summing over n, i, j and taking the limit ∆t, ∆x, ∆y to 0, we obtain

T 0 R 2 ∂ t φ(t, x, y) + ∇W * ρ(t, x, y) • ∇φ(t, x, y) ρ(t, dx, dy) = 0.
Thus ρ is a solution in the sense of distributions of the aggregation equation (1.1). We proceed now as in the proof of Theorem 2.5. Due to the assumptions on the potential, we have that -∇W * ρ ∈ L 2 ((0, T ), L 2 (ρ(t))). Then, we deduce that ρ ∈ AC 2 ([0, T ], P 2 (R 2 )) using [1, Theorem 8.3.1]. By uniqueness of the gradient flow solution and the equivalence Theorem 2.9, we conclude that ρ is the solution of Theorem 2.5. Since the limit is unique, we deduce that the whole sequence is converging towards this limit.

Numerical simulations

We perform in this subsection some numerical simulations obtained by implementing the scheme described above (4.2)-(4.3). We will consider two examples of potential which fit the assumptions (A0)-(A2), that is

W 1 (x) = 1 -e -5|x| ; W 2 (x) = |x|.
For such potentials it is known, see [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF]Section 4], that finite time collapse occurs. More precisely, for any compactly suported initial data, there exists a finite time beyond which the solution is given by a single Dirac Delta mass located at the center of mass. We verify here that we can observe such phenomena thanks to the numerical scheme introduced above.

In our numerical simulations, we consider an initial data given by the sum of three regular bumps:

ρ 0 (x, y) = exp(-C x (x -1/4) 2 -C x (y -1/3) 2 ) + exp(-C x (x -0.8) 2 -C x (y -0.6) 2 ) +0.9 exp(-C x (x -0.4) 2 -C x (y -0.6) 2 ), with c x = 100.
Due to the finite time collapse result, we expect the convergence in finite time of the solution towards a single Dirac Delta. In fact, this is what we observe in Figure 1 for W 1 and in Figure 2 for W 2 . However, comparing the two Figures, the qualitative properties of the convergence towards a single Dirac Delta are not the same depending on the choice of the potential.

In fact, within the dynamics given in Figure 1, we can distinguish two phases in the simulation. In a first phase, we notice the concentration of the density into small masses : we can consider that the numerical solution for time t = 1.8 s is a sum of three numerical Dirac masses with small numerical diffusion. Then these three masses aggregate into two and finally one single mass. On the contrary, for the potential W 2 , we observe in Figure 2 that the numerical solution stays regular and bounded until it forms one single bump and then it collapses.

This tends to indicate the existence of two different time scales: the one corresponding to a radial self-similar collapse onto a single Dirac, and the one corresponding to the interactions between different Dirac Deltas. In the case of the potential W 1 , we observe a faster time scale for the self-similar blow-up of regular solutions into several Dirac Deltas, then the trajectories are given by the sticky particle dynamics for these aggregates. Whereas for the potential W 2 the time scale of the self-similar blow-up is slower compared to the dynamics of the attraction of the aggregates, and then the blow up occurs after all regular bumps aggregate into a single regular bump before the final fate of total collapse.

A very nice feature of this numerical scheme is that it allows for simulations after the first blow-up happens with seemingly good approximation in the measure sense by comparison to the particle simulations, see the one dimensional case [START_REF] James | Numerical method for one-dimensional aggregation equations[END_REF]. The regularization induced on the Dirac Deltas by the numerical diffusion of the scheme does not seem to change the qualitative properties of the solution. Given ε > 0, we use the property (A.1) to get an estimate on the second term in (A.3)

R d \B(x, 1 k ) ∇(W k -W )(x -y)ρ n (dy) ≤ 1 k ≤ ε (A.4) for k ≥ K 1 . Now, we fix K 2 ≥ K 1 such that ρ B(x, 2 K 2 ) \ {x} ≤ ε 4 . (A.5)
We choose a continuous function 0 ≤ ξ ≤ 1 such that ξ(x) = 1 on B(x, 1 K 2 ) and ξ(x) = 0 on R d \ B(x, 2 K 2 ). Then ξ ∈ C c (R d ) and for all k ≥ K 2 , we have

ρ n B(x, 1 k ) \ {x} ≤ ρ n B(x, 1 K 2 ) \ {x} ≤ R d ξ(x)ρ n (dx) ≤ R d ξ(x)(ρ n -ρ)(dx) + ε 4 ,
where we use (A.5) for the last inequality. From the weak convergence as measures of ρ n towards ρ, we have that for n ≥ N 1 large enough

R d ξ(x)(ρ n -ρ)(dx) ≤ ε 4 .
Thus, for k ≥ K 2 we obtain

ρ n B(x, 1 k ) \ {x} ≤ ε 2 ,
uniform in n ≥ N 1 . Therefore, we can bound the first term of the right hand side in (A.3) as

B(x, 1 k )\{x} ∇(W k -W )(x -y)ρ n (dy) ≤ 2w ∞ ρ n (B(x, 1 k ) \ {x}) ≤ w ∞ ε .
Collecting the last inequality with (A.4), we deduce that lim k→∞ x =y ∇W k (xy)ρ n (dy) =

x =y ∇W (xy)ρ n (dy),

We choose a continuous function 0 ≤ ξ ≤ 1 such that ξ(x) = 1 on B(x, 1 N 2 ) and ξ(x) = 0 on R d \ B(x, 2 N 2 ). Then ξ ∈ C c (R d ) and for all n ≥ N 2 , we have

ρ n B(x, 1 n ) \ {x} ≤ ρ n B(x, 1 N 2 ) \ {x} ≤ R d ξ(x)ρ n (dx) ≤ R d ξ(x)(ρ n -ρ)(dx) + ε 4 ,
where we use (A.10) for the last inequality. From the tight convergence of ρ n towards ρ, we have that for n ≥ N 3 large enough (and

N 3 > N 2 ), R d ξ(x)(ρ n -ρ)(dx) ≤ ε 4 .
Thus, for n ≥ N 3

ρ n B(x, 1 n ) \ {x} ≤ ε 2 .
Plugging this latter inequality into (A.9) and from (A.8), we deduce that for n ≥ N 3 ,

x =y

∇W n (xy) -∇W (xy) ρ n (dy) ≤

1 n + C ε 2
Finally, combining this latter inequality with (A.7), we deduce that for n ≥ N 3 ,

a n (x) -a(x)| ≤ ε 4 + 1 n + C ε 2 ,
for a.e. x ∈ R d .

Theorem 2 . 2 (

 22 [START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Side[END_REF]) Let T > 0. Let us assume that the vector field b ∈ L 1 loc (R; L ∞ (

  and Gloyer [8, Theorem 1.2] Theorem 2.4 Let T > 0. Assume that the sequence of vector fields b n converges weakly to b in L 1 ([0, T ], L 1 loc (R d )). Then the Filippov flow X n generated by b n converges locally in C([0, T ]×R d ) to the Filippov flow X generated by b.

Corollary 4 . 1

 41 Let ρ ini ∈ P 2 (R d ), we denote ρ ∈ C([0, T ]; P 2 (R d )) the corresponding solution in Theorem 2.5 with initial data ρ ini . Let ρ ini,N be given in P 2 (R d ) by ρ ini,N (x) = N i=1 m i δ(xx 0

Figure 1 :

 1 Figure 1: Dynamics of the cell density ρ with intial data given by the sum of three bumps in the case W 1 (x) = 1e -5|x| .

Figure 2 :

 2 Figure 2: Dynamics of the cell density ρ with intial data given by the sum of three bumps in the case W 2 (x) = |x|.

  Then, from (4.3), we deduce that for any test function ξ we have on the one hand

										(4.9)
				x, y)ρ ∆ (t, x, y) φ(t, x + ∆x, y) -φ(t, x -∆x, y)
				+a x∆ (t, x + ∆x, y)ρ ∆ (t, x, y) φ(t, x + ∆x, y) -φ(t, x, y)
				+a x∆ (t, x -∆x, y)ρ ∆ (t, x, y) φ(t, x, y) -φ(t, x -∆x, y) dtdxdy.
	Using a Taylor expansion, the mass conservation and the bound (1.2), we deduce from (4.9)
	that
	n,i,j	1 ∆x	a x	n i+1/2j ρ n i+1/2j -a x	n i-1/2j ρ n i-1/2j ψ n i,j = -	1 4	0	T	R 2	2a

x∆ (t, x, y)ρ ∆ (t, x, y) ∂ x φ(t, x, y)

+ a x∆ (t, x + ∆x, y) + a x∆ (t, x -∆x, y) ρ ∆ (t, x, y)∂ x φ(t, x, y) dtdxdy + O(∆x).

(4.10)

Acknowledgements. JAC acknowledges support from projects MTM2011-27739-C04-02, 2009-SGR-345 from Agència de Gestió d'Ajuts Universitaris i de Recerca-Generalitat de Catalunya, the Royal Society through a Wolfson Research Merit Award, and the Engineering and Physical Sciences Research Council (UK) grant number EP/K008404/1. NV acknowledges partial support from the french "ANR blanche" project Kibord : ANR-13-BS01-0004.

Appendix Technical Lemmas

In this appendix we state some technical lemmas which are used in the paper.

Lemma A.1 Let us assume that W satisfies assumptions (A0)-(A2). Let (ρ n ) n∈N be a se-uniformly for n ≥ N 1 . The same argument shows in particular that lim k→∞ x =y ∇W k (xy)ρ(dy) =

x =y ∇W (xy)ρ(dy).

We conclude by passing into the limit k → ∞ in (A.2).

Lemma A.2 Let us assume that W satisfies assumptions (A0)-(A2). Let (W n ) n∈N * be a sequence of even functions in C 1 (R d ) satisfying (A1) and (1.2) with constants λ and w ∞ not depending on n and such that

Let (ρ n ) n∈N be a sequence of measures in P 2 (R d ) such that ρ n ⇀ ρ tightly. Then we have

x =y ∇W (xy)ρ(dy), for a.e. x ∈ R d .

Proof. Let us denote by

∇W n (xy)ρ n (dy), and a(x) := -

x =y ∇W (xy)ρ(dy).

We notice that since W n is even, we have ∇W n (0) = 0, then

∇W n (xy)ρ n (dy).

Let ε > 0, from Lemma A.1 we deduce that there exists N 1 ∈ N * such that for all n ≥ N 1 ,

Then using (A.6), we deduce that

∇W n (xy) -∇W (xy) ρ n (dy).

(A.8) Now, we proceed as in the proof of Lemma A.1. From assumptions on W n and W , we deduce (see (1.2)) that there exists a constant C such that