Convergence analysis of upwind type schemes for the aggregation equation with pointy potential

A numerical analysis of upwind type schemes for the nonlinear nonlocal aggregation equation is provided. In this approach, the aggregation equation is interpreted as a conservative transport equation driven by a nonlocal nonlinear velocity ﬁeld with low regularity. In particular, we allow the interacting potential to be pointy, in which case the velocity ﬁeld may have discontinuities. Based on recent results of existence and uniqueness of a Filippov ﬂow for this type of equations, we study an upwind ﬁnite volume numerical scheme and we prove that it is convergent at order 1 / 2 in Wasserstein distance. The paper is illustrated by numerical simulations that indicate that this convergence order should be optimal.


F R A N Ç OIS D E L A R U E F R É D É R IC L A G O U TIÈ R E NIC O L A S VA U C H E L E T C O N V E R G E N C E A N A LY SIS O F U P W IN D T Y P E S C H E M E S F O R T H E A G G R E G ATIO N E Q U ATIO N W IT H P OIN T Y P O T E N TI A L A N A LY S E D E C O N V E R G E N C E D E S C H É M A S D E T Y P E U P W IN D P O U R L ' É Q U ATIO N D ' A G R É G ATIO N AV E C U N P O T E N TI A L P OIN T U
Abstract. -A numerical analysis of upwind type schemes for the nonlinear nonlocal aggregation equation is provided. In this approach, the aggregation equation is interpreted as a conservative transport equation driven by a nonlocal nonlinear velocity field with low regularity. In particular, we allow the interacting potential to be pointy, in which case the velocity field may have discontinuities. Based on recent results of existence and uniqueness of Keywords: Aggregation equation, upwind finite volume scheme, convergence order, measure-valued solution. 2010 Mathematics Subject Classification: 35B40, 35D30, 35L60, 35Q92, 49K20. DOI: https://doi.org/10.5802/ahl.30 (*) The authors acknowledge partial support from the french "ANR blanche" project Kibord: ANR-13-BS01-0004, as well as from the "BQR Acceuil EC 2017" grant from Université Lyon 1. a Filippov flow for this type of equations, we study an upwind finite volume numerical scheme and we prove that it is convergent at order 1/2 in Wasserstein distance. The paper is illustrated by numerical simulations that indicate that this convergence order should be optimal.

Introduction
This paper is devoted to the numerical approximation of measure valued solutions to the so-called aggregation equation in space dimension d. This equation reads with the initial condition ρ(0, · ) = ρ ini . Here, W plays the role of an interaction potential whose gradient ∇ x W (x − y) measures the relative force exerted by a unit mass localized at a point y onto a unit 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 [BCP97,CMV06,Tos04]. Models of crowd motion with a nonlinear term of the form ∇ x W * ρ are also addressed in [CGLM12,CLM13]. In population dynamics, (1.1) provides a biologically meaningful description of aggregative phenomena. For instance, the description of the collective migration of cells by swarming leads to such a kind of PDEs with non-local interaction, see e.g. [MCO05,OL02,TB04]. 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 one another. The dynamics can be macroscopically modelled by the Patlak-Keller-Segel system [KS70,Pat53]. In the kinetic framework, the most frequently used model is the Othmer-Dunbar-Alt system, the hydrodynamic limit of which leads to the aggregation equation (1.1), see [DS05,FLP05,JV13]. In many of these examples, the potential W is usually mildly singular, i.e. W has a weak singularity at the origin. Because of this low regularity, smooth solutions of such systems may blow-up in finite time, see e.g. [BV06,BLR11,CDF + 11,LT04]. In the latter case, finite time concentration may be regarded as a very simple mathematical way to account for aggregation of individuals, as opposed to diffusion.
Since finite time blow-up of smooth solutions may occur and since equation (1.1) conserves mass, a natural framework to study the existence of global in time solutions is to work in the space of probability measures. In this regard, two strategies have been proposed in the literature. In [CDF + 11], the aggregation equation is seen as a gradient flow taking values in the Wasserstein space and minimizing the interaction energy.
In [CJLV16,JV13,JV16,LV16], this system is considered as a conservative transport equation with velocity field ∇ x W * ρ. Then a unique flow, say Z = (Z(t, · )) t 0 , can be constructed, hence allowing to define the solution as a pushforward measure by the flow, namely ρ = (ρ(t) = Z(t, · ) # ρ ini ) t 0 . When the singularity of the potential is stronger than the mild form described above, such a construction has been achieved in the radially symmetric case in [BGL12], but uniqueness is then lacking. Actually, the assumptions on the potential W that are needed to ensure the well-posedness of the equation in the space of measure valued solutions require a certain convexity property of the potential that allows only for a mild singularity at the origin. More precisely, we assume that the interaction potential W : R d → R satisfies the following properties: (A1) W (x) = W (−x) and W (0) = 0; (A2) W is λ-convex for some λ ∈ R, i.e. W (x) − λ 2 |x| 2 is convex; (A3) W ∈ C 1 (R d \ {0}); (A4) W is Lipschitz-continuous. Such a potential will be referred to as a pointy potential. Typical examples are fully attractive potentials W (x) = 1 − e −|x| , which is −1-convex, and W (x) = |x|, which is 0-convex. Notice that the Lipschitz-continuity of the potential allows to bound the velocity field: there exists a nonnegative constant w ∞ such that for all x = 0, Observe also that (A4) forces λ in (A2) to be non-positive, as otherwise W would be at least of quadratic growth, whilst (A4) forces it to be at most of linear growth. However, we shall sometimes discard (A4), when the initial datum is compactly supported. In this case, as W − λ|x| 2 /2 is convex, it is locally Lipschitz-continuous, so that W is locally Lipschitz-continuous, what will be sufficient for compactly supported initial data. In that case it perfectly makes sense to assume λ > 0 in (A2). For numerical analysis, we will assume in this case that the potential is radial, that is to say that W is a function of the sole scalar |x|, W (x) = W(|x|). Although very accurate numerical schemes have been developed to study the blowup profile for smooth solutions, see [HB10,HB12], very few numerical schemes have been proposed to simulate the behavior of solutions to the aggregation equation after blow-up. The so-called sticky particles method was shown to be convergent in [CDF + 11] and used to obtain qualitative properties of the solutions such as the time of total collapse. However, this method is not so practical to catch the behavior of the solutions after blow-up in dimension d larger than one. In dimension d = 1, this question has been addressed in [JV13]. In higher dimension, particle methods have been recently proposed and studied in [CB16,CPCCC15], but only the convergence of smooth solutions, before the blowup time, has been proved. Finite volume schemes have also been developed. In [JV15], the authors propose a finite volume scheme to approximate the behavior of the solution to the aggregation equation (1.1) after blow-up and prove that it is convergent. A finite volume method for a large class of PDEs including in particular (1.1) has been also proposed in [CCH15], but no convergence result has been given. Finally, a finite volume scheme of Lax-Friedrichs type for general measures as initial data has been introduced and investigated in [CJLV16]. Numerical simulations of solutions in dimension greater than one have been obtained, allowing to observe the behavior after blow-up. Moreover, convergence towards measure valued solutions has been proved. However, no estimate on the order of convergence has been established so far. In the current work, we provide a precise estimate of the order of convergence in Wasserstein distance for an upwind type scheme. This scheme is based on an idea introduced in [JV13] and used later on in [CJLV16,JV15]. It consists in discretizing properly the macroscopic velocity so that its product with the measure solution ρ is well-defined. In this paper, we introduce an upwind scheme for which this product is treated accurately, and we prove its convergence at order 1/2 in Wasserstein distance (the definition of which is recalled below).
For a given velocity field, the study of the order of convergence for the finite volume upwind scheme for the transport equation has received a lot of attention. This scheme is known to be first order convergent in L ∞ norm for any smooth initial data in C 2 (R d ) and for well-suited meshes, provided a standard stability condition (Courant-Friedrichs-Lewy condition) holds, see [BGP05]. However, this order of convergence falls down to 1/2 in L p norm when considering non-smooth initial data or more general meshes. This result has been first proved in the Cartesian framework by Kuznetsov in [Kuz76]. In [Des04], a 1/2 order estimate in the L ∞ ([0, T ], L 2 (R d )) norm for H 2 (R d ) initial data has been established. Finally in [DL11,MV07], a 1/2 order estimate in L 1 has been proved for initial data in L 1 (R d ) ∩ BV (R d ), whilst, for Lipschitz-continuous initial data, an estimate of order 1/2−ε in L ∞ for any ε > 0 has been obtained in [DL11,Mer07]. We emphasize that the techniques used in [Mer07,MV07] and [DL11] are totally different. In the former, the strategy of proof is based on entropy estimates, whereas in the latter, the proof relies on the construction and the analysis of stochastic characteristics for the numerical scheme. Finally, when the velocity field is only L ∞ and one-sided Lipschitz-continuous, solutions of the conservative transport equation are defined only in the sense of measures. In this regard, Poupaud and Rascle [PR97] have proved that solutions of the conservative transport equation could be defined as the pushforward of the initial condition by a flow of characteristics. A stability estimate for such solutions has been stated later in [BG11]. In dimension d = 1, these solutions, as introduced in [PR97], are equivalent to duality solutions, as defined in [BJ98]. Numerical investigations may be found in [GJ00]. In such a framework with a low regularity, numerical analysis requires to work with a sufficiently weak topology, which is precisely what has been done in [DLV17]. Therein, the convergence at order 1/2 of a finite volume upwind scheme has been shown in Wasserstein distance by means of a stochastic characteristic method, as done in [DL11]. Observe also that, recently, such an approach has been successfully used in [SS17] for the numerical analysis of the upwind scheme for the transport equation with rough coefficients. In the current work, we adapt the strategy initiated in [DLV17] to prove the convergence at order 1/2 of an upwind scheme for the aggregation equation for which the velocity field depends on the solution in a nonlinear way. We will strongly use the fact that, as mentioned above, measure valued solutions of (1.1) are constructed by pushing forward the initial condition by an R d -valued flow. Noticeably, we entirely reformulate the stochastic approach used in [DLV17] by means of analytical tools. In the end, our proof is completely deterministic. Although using analytical instead of probabilistic arguments do not change the final result (neither nor the general philosophy of the proof), it certainly makes the whole more accessible for the reader. As we pointed out, the key fact in [DLV17] is to represent the scheme through a Markov chain; here, the main idea is to use the sole transition kernel of the latter Markov chain to couple the measurevalued numerical solution at two consecutive times (and hence to bypass any use of the Markov chain itself). We refer to Remark 4.2 below for more details.
The outline of the paper is the following. In the next section, we introduce the notations and recall the theory for the existence of a measure solution to (1.1). Then we present the upwind scheme and state the main result: the scheme is convergent at order 1/2. In case when the potential W is strictly convex and radially symmetric and the initial condition has a bounded support, the rate is claimed to be uniform in time. Section 3 is devoted to the properties of the scheme. The proof of the main result for a Cartesian grid mesh is presented in Section 4. In Section 5, we explain briefly how to extend our result to simplicial meshes. Finally, numerical illustrations are given in Section 6. In particular, we show that the order of convergence is optimal and we provide several numerical simulations in which we recover the behavior of the solutions after blow-up time.

Notations
Throughout the paper, we will make use of the following notations. We denote by C 0 (R d ) the space of continuous functions from R d to R that tend to 0 at ∞. We denote by M b (R d ) the space of Borel signed measures whose total variation is finite. For ). For ρ a measure in M b (R d ) and Z a measurable map, we denote Z # ρ the pushforward measure of ρ by Z; it satisfies, for any continuous function φ, We call P(R d ) the subset of M b (R d ) of probability measures. We define the space of probability measures with finite second order moment by Here and in the following, | · | 2 stands for the square Euclidean norm, and · , · for the Euclidean inner product. where Γ(µ, ν) is the set of measures on R d × R d with marginals µ and ν, i.e.
By a minimization argument, we know that the infimum in the definition of d W is actually a minimum. A measure 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 We will make use of the following properties of the Wasserstein distance. Given µ ∈ P 2 (R d ) and two µ-square integrable Borel measurable maps X, Y :

Existence of a unique flow
In this section, we recall the existence and uniqueness result for the aggregation equation (1.1) obtained in [CJLV16] (and extend it a bit for non-globally Lipschitzcontinuous potentials). For ρ ∈ C([0, T ]; P 2 (R d )), we define the velocity field a ρ by where we have used the notation Due to the λ-convexity of W , see (A3), we deduce that, for all x, y in R d \ {0}, Moreover, since W is even, ∇W is odd and by taking y = −x in (2.3), we deduce that inequality (2.3) still holds for ∇W , even when x or y vanishes: This latter inequality provides a one-sided Lipschitz-continuity (OSL) estimate for the velocity field a ρ defined in (2.2), i.e. we have We recall that, for a velocity field b ∈ L ∞ ([0, +∞); L ∞ (R d )) d satisfying an OSL estimate, i.e.
In this definition, {Convess( a ρ )(t, · )}(x) denotes the essential convex hull of the vector field a ρ (t, · ) at x. We remind briefly the definition for the sake of completeness (see [AC84,Fil64] for more details). We denote by Conv(E) the classical convex hull of a set E ⊂ R d , i.e., the smallest closed convex set containing E. Given the vector field a ρ (t, · ) : R d → R d , its essential convex hull at point x is defined as where N 0 is the set of zero Lebesgue measure sets. Moreover, we have the semi-group property: for any t, τ, s ∈ [0, +∞) such that t τ s and x ∈ R d , From now on, we will make use of the notation Z(t, x) = Z(t; 0, x). Using this characteristic, it has been established in [PR97] that solutions to the conservative transport equation with a given bounded and one-sided Lipschitz-continuous velocity field could be defined as the pushforward of the initial condition by the Filippov characteristic flow. Based on this approach, existence and uniqueness of solutions to (1.1) defined by a Filippov flow has been established in [CJLV16]. More precisely the statement reads: Theorem 2.1 ([CJLV16, Theorem 2.5 and 2.9]). -(1) Let W satisfy assumptions (A1)-(A4) and let ρ ini be given in P 2 (R d ). Then, there exists a unique solution ρ ∈ C([0, +∞); P 2 (R d )) satisfying, in the sense of distributions, the aggregation equation where a ρ is defined by (2.2). This solution may be represented as the family of pushforward measures (ρ(t) := Z ρ (t, · ) # ρ ini ) t 0 where (Z ρ (t, · )) t 0 is the unique Filippov characteristic flow associated to the velocity field a ρ .
Moreover, the flow Z ρ is Lipschitz-continuous and we have At last, if ρ and ρ are the respective solutions of (2.6) with ρ ini and ρ ini, as initial conditions in P 2 (R d ), then (2) Let W satisfy (A1)-(A3) and be radial, λ be (strictly) positive and let ρ ini be given in P 2 (R d ) with compact support included in B ∞ (M 1 , R), where M 1 is the first moment of ρ ini (i.e. its center of mass) and B ∞ (M 1 , R) the closed ball for the infinite norm on R d centered at M 1 with radius R. Then, there exists a unique solution ρ ∈ C([0, +∞); P 2 (R d )) with support included in B ∞ (M 1 , R) satisfying, in the sense of distributions, the aggregation equation (2.6) where a ρ is defined by (2.2). Moreover, the flow Z ρ is Lipschitz-continuous and we have At last, if ρ ini and ρ ini, have a bounded support, then, The stability estimates that are present in this result are Dobrushin type estimates in the quadratic Wasserstein distance, in the case where the kernel is not Lipschitzcontinuous but only one-sided Lipschitz-continuous. See [Dob79] and [Gol16].
We mention that the solution, which is here represented by the Filippov characteristic flow, may be also constructed as a gradient flow solution in the Wasserstein space P 2 (R d ), see [CDF + 11]. Here it is also important to remark that (2.7) is true under the sole assumptions (A1)-(A3) whenever λ > 0 (which is a mere consequence of (2.9) and (2.10) below). In that case, it ensures that B 2 (M 1 , R) (the closed Euclidean ball) is preserved by the flow without the assumption that W is radial. As a result, it may be tempting to address the analysis below without requiring the potential to be radial. Nevertheless, the problem is that the numerical scheme does not satisfy a similar property. Indeed, the Euclidean ball B 2 (M 1 , R) is not convex from a numerical point of view, that is to say, if we regard the mesh underpinning the scheme, then the union of the square cells whose center is included in B 2 (M 1 , R) is not convex. Due to this drawback, the flow associated to the scheme does not preserve the ball B 2 (M 1 , R). This is in contrast with Lemma 3.3 below, which shows that, in the radial setting, the ball B ∞ (M 1 , R + ∆x) is kept stable by the scheme, where ∆x is the step of the spatial mesh. This latter fact is the reason why we here assume that the potential is radial.
Proof. -For the first two statements of the Theorem, existence of a unique solution and Lipschitz-continuity of the flow, we refer to [CJLV16]. These statements remain true whenever the sole (A1)-(A3) hold true, W is radial, λ is (strictly) positive and the support of ρ ini is bounded, provided that the notion of solution is limited to collections (ρ(t, · )) t 0 that have a compact support, uniformly in t in compact subsets. Indeed, if we denote by M 1 (t) the center of mass of the solution at time t, namely M 1 (t) := R d ρ(t, dx), then this center of mass is known to be preserved: M 1 (t) = M 1 (0) =: M 1 (see [CJLV16] or Lemma 3.2 below for the discrete counterpart). Now, if λ 0 and if W is radial, is preserved by the flow and guarantees that ρ(t) has its support included in B ∞ (M 1 , R) for any time t 0, if it is the case for t = 0. Given the fact that the support of ρ(t) remains bounded in B ∞ (M 1 , R), everything works as if W was globally Lipschitzcontinuous. Existence and uniqueness of a solution to the aggregation equation can thus be proved by a straightforward localization argument. Indeed, observe that from the very definition of the velocity a, the Lipschitz-continuity constant of W that is involved in the existence and uniqueness theory is the local one of W on the compact Now it only remains to prove the two inequalities regarding the Wasserstein distance between solutions starting from different data. Under assumptions (A1)-(A4) on the potential, it was proven in [CJLV16], but with a constant 2|λ| instead of |λ| in the exponential (as in [Dob79] and [Gol16], where the convolution operator is however replaced with a slightly more general integral operator), thus we here provide a proof of the present better estimate.
We consider the two Filippov flows (Z ρ (t, · )) t 0 and (Z ρ (t, · )) t 0 as defined in the statement of Theorem 2.1. We recall that To simplify, we just write Z(t, · ) = Z ρ (t, · ) and Z (t, · ) = Z ρ (t, · ). Then, for any x, y ∈ R d and for a.e. t 0, Call π ∈ Γ 0 (ρ ini , ρ ini, ) an optimal plan between ρ ini and ρ ini, . Then, This identity must be understood as an expression for |Z(t, x) − Z (t, y)| 2 in terms of the time integral of the right-hand side. Integrating the latter in (x, y) with respect to π, we get Thanks to the fact that ∇W is odd, see (A1), we can write, by a symmetry argument, Using (2.4), we obtain Observe that the above right-hand side is equal to 1st case. -If λ 0, we deduce from (2.9) and (2.10) that which suffices to complete the proof of the first claim by noting that 2nd case. -If λ 0, we just use the fact that the right-hand side in (2.9) is non-positive. Proceeding as above, this permits to complete the proof of the second claim.

Main result
The aim of this paper is to prove the convergence at order 1/2 of an upwind type scheme in distance d W for the aggregation equation. The numerical scheme is defined as follows. We denote by ∆t the time step and consider a Cartesian grid with step ∆x i in the ith direction, i = 1, . . . , d; we then let ∆x := max i ∆x i . We also introduce the following notations. For a multi-index For a given nonnegative measure ρ ini ∈ P 2 (R d ), we put, for any J ∈ Z d , Since ρ ini is a probability measure, the total mass of the system is J∈Z d ρ 0 J = 1. We then construct iteratively the collection ((ρ n J ) J∈Z d ) n∈N , each ρ n J being intended to provide an approximation of the value ρ(t n , x J ), for J ∈ Z d . Assuming that the approximating sequence (ρ n J ) J∈Z d is already given at time t n := n∆t, we compute the approximation at time t n+1 by: The notation (a) + = max{0, a} stands for the positive part of the real a and respectively (a) − = max{0, −a} for the negative part. The macroscopic velocity is defined by Since W is even, we also have: (2.14) The main result of this paper is the proof of the convergence at order 1/2 of the above upwind scheme. More precisely the statement reads: (1) Assume that W satisfies hypotheses (A1)-(A4) and that the so-called strict 1 2 -CFL condition holds: with w ∞ as in (1.2).
For ρ ini ∈ P 2 (R d ), let ρ = (ρ(t)) t 0 be the unique measure solution to the aggregation equation with initial data ρ ini , as given by Theorem 2.1. Define Then, there exists a nonnegative constant C, only depending on λ, w ∞ and d, such that, for all n ∈ N * , (2) Assume that W is radial and satisfies hypotheses (A1)-(A3) with λ (strictly) positive, that ρ ini is compactly supported in B ∞ (M 1 , R) where M 1 is the center of mass of ρ ini , and that the CFL condition (2.15) holds, with w ∞ defined as Assume also that ∆t 1/2 and 2λ∆t < 1. Then, there exists a nonnegative constant C, only depending on λ, w ∞ , d and R such that, for all n ∈ N * , (2.16) is valid, as well as which proves that the error can be uniformly controlled in time.
We stress the fact that, under the setting defined in (2), (2.16) is valid. In small time, it provides a better estimate than (2.18). As indicated in the statement, the constant C in (2.18) may depend on the value of R in the assumption Supp(ρ ini ) ⊂ B ∞ (M 1 , R).
We also point out that, although the computations below are performed for the sole upwind scheme, the first part of the statement, which holds true under the full set of hypotheses (A1)-(A4), can be straightforwardly adapted to other diffusive schemes, see for instance our previous article [DLV17]. As for (2), the statement remains true provided that the supports of the approximating measures (ρ n ) n 0 remain bounded as n grows up. It must be stressed that there are some schemes for which the latter property fails (e.g. Lax-Friedrichs' scheme).
Moreover, as already mentioned in Introduction, the convergence rate is optimal; this latter fact will be illustrated by numerical examples in Section 6.
Example 2.3. -In one dimension, the scheme (2.12) reads where i is just taken in Z. The scheme has then the following interpretation. Given ρ n ∆x = j∈Z ρ n j δ x j , we construct the approximation at time t n+1 by implementing the following two steps: • The Delta mass ρ n i located at position x i moves with velocity a n i to the position x i + a n i ∆t. Under the CFL condition w ∞ ∆t ∆x (which is obviously weaker than what we require in (2.15)), the point x i + a n i ∆t belongs to the interval [x i , x i+1 ] if a n i 0, and to the interval [x i−1 , x i ] if a n i 0.

ANNALES HENRI LEBESGUE
• Then the mass ρ n i is split into two parts; if a n i 0, a fraction a n i ∆t/∆x of it is transported to the cell i + 1, while the remaining fraction is left in cell i; if a n i 0, the same fraction |a n i |∆t/∆x of the mass is not transported to the cell i + 1 but to the cell i − 1. This procedure may be regarded as a linear interpolation of the mass ρ n i between the points x i and x i+1 if a n i 0 and between the points x i and x i−1 if a n i 0. This interpretation holds only in the one dimensional case. However thanks to this interpretation, we can define a forward semi-Lagrangian scheme in any dimension on (unstructured) simplicial meshes, which is then different from (2.12). Such a scheme is introduced in Section 5.
Finally, we emphasize that this scheme differs from the standard finite volume upwind scheme in which the velocity is computed at the interface a n i+1/2 . This subtlety is due to the particular structure of the equation, as the latter requires the product a ρ ρ to be defined properly. A convenient way to make it proper is to compute, in the discretization, the velocity and the density at the same grid points. This fact has already been noticed in [JV15, GV16] and is also illustrated numerically in Section 6.

Properties of the scheme
The following lemma explains why we called CFL the condition on the ratios (∆t/∆x i ) i=1,...,d that we formulated in the statement of Theorem 2.2. Proof. -The total initial mass of the system is J ρ 0 J = 1. By summing equation (2.12) over J, we can show that the mass is conservative, namely, for all n ∈ N * , We prove by induction on n that ρ n J 0 for all J ∈ Z d and for all n ∈ N. Indeed, if, for some n ∈ N, it holds ρ n J 0 for all J ∈ Z d , then, by definition (2.13) and assumption (1.2), we clearly have Then, assuming that the condition (2.15) holds, we deduce that, in the relationship (3.1), all the coefficients in front of ρ n J , ρ n J−e i and ρ n J+e i , i = 1, . . . , d, are nonnegative. Thus, using the induction assumption, we deduce that ρ n+1 In the following lemma, we collect two additional properties of the scheme: the conservation of the center of mass and the finiteness of the second order moment. 11). Then, the sequence (ρ n J ) J∈Z d given by the numerical scheme (2.12)-(2.13) satisfies: (1) Conservation of the center of mass. For all n ∈ N * , We conclude the proof using a discrete version of Gronwall's lemma.
In case when W is radial and satisfies (A1)-(A3), λ is (strictly) positive and ρ ini has a bounded support, Lemmas 3.1 and 3.2 become: Lemma 3.3. -Assume that W is radial and satisfies (A1)-(A3), λ is (strictly positive) and ρ ini has a bounded support, then the conclusions of Lemmas 3.1 and 3.2 remain true provided that w ∞ is defined as in (2.17).
Moreover, for any R 0 such that Supp(ρ ini ) ⊂ B ∞ (M 1 , R), it holds, for any n ∈ N, The meaning of Lemma 3.3 is pretty clear. For R as in the statement, the mass, as defined by the numerical scheme, cannot leave the ball B ∞ (M 1,∆x , R + ∆x). We here recover the same idea as in Theorem 2.1.
Proof. -As long as we can prove that the mass, as defined by the numerical scheme, cannot leave the ball B ∞ (M 1,∆x , R + ∆x), the proof is similar to that of Lemmas 3.1 and 3.2. So, we focus on the second part of the statement.
We first recall that ρ 0 Below, we prove by induction that the same holds true for any n ∈ N. To do so, we assume that there exists an integer n ∈ N such that, for all J ∈ Z d , ρ n The goal is then to prove that, for any J satisfying (3.2), ρ n+1 J = 0. By (3.1), it suffices to prove that, for any coordinate i ∈ {1, . . . , d} and any J as in (3.2),  ∆x ) i 0 and the argument below is the same). Hence, (x J+e i 0 ) i 0 > R+∆x+(M 1,∆x ) i 0 and, by the induction hypothesis, ρ n J+e i 0 = 0, which proves the first equality in (3.3) when i = i 0 . In order to prove the second equality when i = i 0 , we notice from (2.13) that As W is radial and λ > 0, ∇W (x − y) is positively proportional to x − y. Hence, Remark 3.4. -Lemma 3.3 is the main rationale for requiring W to be radial. Indeed, the counter-example below shows that the growth of the support of ρ ini can be hardly controlled whenever λ > 0 and W is just assumed to satisfy (A1)-(A3). Consider for instance the following potential in dimension d = 2: where q is a free integer whose value will be fixed later on. It is well checked that Standard computations show that the smallest eigenvalue of the Hessian matrix (which is independent of (x 1 , x 2 )) is so that W is λ-convex with λ converging to 1/2 as q tends to ∞.

ANNALES HENRI LEBESGUE
Take now a centered probability measure ρ and compute the first coordinate of the velocity field a ρ . By centering, In particular, if x 2 = 1, then ( a ρ ) 1 (x 1 , 1) = q − x 1 , which is non-negative as long as x 1 < q. Therefore, if the numerical scheme is initialized with some centered ρ 0 ∆x supported by the unit square [−1, 1] 2 , it holds ( a ρ 0 ∆x ) 1 (1, 1) > 0, if q > 1. Hence, provided that condition (2.15) holds true, ρ 1 ∆x charges the point (1 + ∆x, 1). Since the numerical scheme preserves the centering, we also have ( a ρ 1 ∆x ) 1 (1 + ∆x, 1) > 0, if q > 1 + ∆x, and then ρ 2 ∆x also charges the point (1 + 2∆x, 1), and so on up until (∆x q/∆x , 1). This says that there is no way to control the growth of the support of the numerical solution in terms of the sole lower bound of the Hessian matrix. Somehow, the growth of ∇W plays a key role. This is in stark contrast with the support of the real solution, which may be bounded independently of q, as emphasized in the proof of Theorem 2.1.
A possible way to overcome the fact that the numerical scheme does not preserve any ball containing the initial support in the general case when W is not radial would be to truncate the scheme. We feel more reasonable not to address this question in this paper, as it would require to revisit in deep the arguments used to tackle the case λ 0.

Comparison with a potential non-increasing scheme
It must be stressed that the scheme could be defined differently in order to force the potential (or total energy: R d ×R d W (x − y) ρ(dx) ρ(dy)) to be non-increasing. Basically, this requires the velocity a to be defined as a discrete derivative.
For simplicity, we provide the construction of the scheme in dimension 1 only. For a probability measure ∈ P(Z) and a cell I ∈ Z, we consider the following two discrete convolutions of finite differences: where, as before, ∆x is obtained by pushing forward by the mapping y → ∆x y. The two terms above define velocities at the interfaces of the cell I. Namely, we call the first term −a . Of course, the sign − in the former term guarantees the consistency of the notation, that is a Following (2.12), the scheme is defined by: for n ∈ N and J ∈ Z. It is shown in [CCH15] that the potential is non-increasing for the semi-discretized version of this scheme, which is to say that, up to a remainder of order 2 in ∆t (the value of ∆x being fixed), the potential of the fully discretized scheme does not increase from one step to another. The proof of the latter claim follows from a direct expansion of the quantity by using the updating rule for ρ n+1 J in terms of ρ n J , ρ n J−1 and ρ n J+1 . The numerical scheme investigated in this paper does not satisfy the same property. Indeed, we provide a counter example, which shows that the potential may increase when W is convex, as a consequence of the numerical diffusion. However, the same example, but in dimension 1, shows that the scheme (3.4) may not be convergent for certain forms of potential for which Theorem 2.2 applies, see Subsection 6.3.
Then, denoting by ρ 1 the distribution at time 1 obtained by implementing the upwind scheme, it holds that: where the Landau symbol O( · ) may depend upon p.
Choosing p > 1/2 in (3.5), we see that the potential may increase at the same rate as the time step.
Proof. -We first compute the potential at time 0. To do so, we compute R 2 |x − y|ρ 0 (dy), for x ∈ {0, e 1 , e 2 }: In order to compute the potential at time 1, we compute the velocity at each of the above points. Observing that the velocity at point x is given by the formula: with the convention 0 0 = 0, we get: We then compute the new masses at time 1. There is one additional point which is charged: e 1 + e 2 = (1, 1). We have: We now have all the required data to compute the potential at time 1.
Finally, the potential at time 1 is given by: We expand the above right-hand side in powers of ∆t. The zero-order term is exactly equal to R 2 R 2 |x − y|ρ 0 (dx)ρ 0 (dy). So, we just compute the terms in ∆t. It is equal to which completes the proof.

Order of convergence
This section is devoted to the proof of Theorem 2.2.

Preliminaries
Before presenting the proof, we introduce some notations and establish some useful properties. We first define the following interpolation weights: for J ∈ Z d and y ∈ R d , we let The terminology interpolation weights is justified by the following straightforward observation. Given a collection of reals (h J ) J∈Z d indexed by the cells of the mesh, which we may regard as a real-valued function h : x J → h J defined at the nodes of the mesh, we may define an interpolation of h = (h J ) J∈Z d by letting Obviously, the sum in the right-hand side makes sense since only a finite number of weights are non-zero for a given value of y. Clearly, the functional I is an interpolation operator. As explained below, I makes the connection between the analysis we perform in this paper and the one we performed in our previous work [DLV17]. Several crucial facts must be noticed. The first one is that, contrary to what one could guess at first sight, the weights are not necessarily non-negative. For a given J ∈ Z d , take for instance y = (y i = (J i − 1 2 )∆x i ) i=1,...,d ∈ C J . Then α J (y) = 1 − d 2 , which is obviously negative if d 3. However, the second point is that, for useful values of y, the weights are indeed non-negative provided that the CFL condition (2.15) is in force. For a given J ∈ Z d , call indeed U J the subset of C J of so-called useful values that are in C J , as given by Then, for any J, L ∈ Z d and any y ∈ U L , α J (y) is non-negative, which is a direct consequence of the CFL condition (2.15). In fact, the CFL condition (2.15) says more, and this is the rationale for the additional factor 1 2 in (2.15): U J is included in C J . Of course, the consequence is that, under the CFL condition (2.15), we have, for any J ∈ Z d , x J + a n J ∆t ∈ C J , where a n J is the d-dimensional vector with entries Another key fact is that the definition of α J (y) in (4.1) is closely related to the definition of the numerical scheme (2.12). Indeed, we have the following formula, for any J, L ∈ Z d , In particular, we may rewrite (2.12) as which is the core of our analysis below. In this regard, The following lemma gathers some useful properties.
Lemma 4.1. -Let (α L (y)) L∈Z d ,y∈R d be defined as in (4.1). Then, for any y ∈ R d , we have Proof. -There exists a unique J ∈ Z d such that y ∈ C J . Then, we compute Then, using the fact that which completes the proof.
Remark 4.2. -Lemma 4.1 prompts us to draw a comparison with our previous paper [DLV17]. For a given y ∈ R d in the set of useful values U := ∪ J∈Z d U J , namely y ∈ U J for some J ∈ Z d , the collection of weights (α L (y)) L∈Z d forms a probability measure, as the weights are non-negative and their sum is 1! In particular, I(h)(y) in (4.2), for y ∈ U , may be interpreted as an expectation.
Using the same terminology as in [DLV17] (which is in fact the terminology of the theory of Markov chains), those weights should be regarded as transition probabilities: For a given y in the set of useful values, α L (y) reads as the probability of jumping from a certain state depending on the sole value of y to the node x L . Of course, the interpretation of the so-called certain state depending on the sole value of y is better understood from (4.3). In (4.3), if we fix a cell L ∈ Z d (or equivalently a node x L ), then α J (x L +∆ta n L ) should read as the probability of passing from the node x L to the node x J (or from the cell L to the cell J) at the n th step of a (time inhomogeneous) Markov chain having the collection of nodes (or of cells) as state space. In this regard, (4.4) is nothing but the Kolmogorov equation for the corresponding Markov chain, as (ρ n J ) J∈Z d can be interpreted as the law at time n of the Markov chain driven by the latter transition probabilities. The reader can easily check that the so-called stochastic characteristic used in [DLV17] is in fact this Markov chain.
Below, we do not make use of the Markov chain explicitly. Still, we use the weights (α J (y)) J∈Z d ,y∈R d to construct a coupling between the two measures ρ n ∆x and ρ n+1 ∆x , that is to construct a specific element of Γ(ρ n ∆x , ρ n+1 ∆x ). In [DLV17], this coupling does not explicitly show up but it is in fact implicitly used, as it coincides with the joint law of two consecutive states of the aforementioned Markov chain.
In a nutshell, the reader can reformulate the whole analysis below in a probabilistic fashion. The only (conceptual) difficulty to do so is that, in contrast with [DLV17], the Markov chain is here nonlinear: as a n in (2.13) depends on ρ n , the transition probabilities of the Markov do depend upon the marginal law of the Markov chain itself, which fact gives rise to a so-called nonlinear Markov chain!

Proof of Theorem 2.2
1st step. -We first consider the case where the initial datum is given by where we recall that ρ 0 J is defined in (2.11). For n ∈ N * , let us define D n := d W (ρ(t n ), ρ n ∆x ). Clearly, with our choice of initial datum, we have D 0 = 0.
Let γ be an optimal plan in Γ 0 (ρ(t n ), ρ n ∆x ), we have Let us introduce a n ∆x , the piecewise affine in each direction reconstruction of the velocity such that for all J ∈ Z d , a n ∆x (x J ) = a n J Denote also by Z := Z ρ the flow given by Theorem 2.1, when ρ ini is prescribed as above. Recalling the definition of α J (y) from (4.1), we then consider a new measure γ , defined as the image of γ by the kernel K that associates with a point (x, y) ∈ R d × R d the point (Z(t n+1 ; t n , x), x L ) with measure α L (y + ∆ta n ∆x (y)), namely, for any two Borel subsets A and B of R d , where δ z denotes the Dirac mass at point z, and then Equivalently, for any bounded Borel-measurable function θ : ; t n , x), x L )α L (y + ∆ta n ∆x (y)) γ(dx, dy).
(1) The probabilistic reader will easily recognize the standard computation of the L 2 norm of a random variable in terms of its variance and its expectation, which indeed plays, but under a conditional form, a key role in [DLV17].

Thus, equation (4.7) rewrites
Injecting into (4.6), we deduce where we used the fact that ρ n ∆x is the second marginal of γ. By definition, ρ n ∆x (y) = J∈Z d ρ n J δ J (y), so that Moreover using the definition of α L in (4.1), we compute where we used, for the last inequality, the CFL condition (2.15) and the fact that the velocity (a n J ) J is uniformly bounded (see Lemma 3.1 or Lemma 3.3). Then, (4.8) gives 2nd step. -We estimate the error between the exact characteristic Z(t n+1 ; t n , x) and the forward Euler discretization y + ∆ta n ∆x (y). By definition of the characteristics (2.5), we have ∇W (Z(s; t n , x) − Z(s; t n , ξ))ρ(t n , dξ)ds.
By conservation of the center of mass, see Lemma 3.2 (1), we deduce that since we have chosen the initial data such that ρ ini = ρ 0 ∆x . Using also the bound |Z(s; t n , x) − x| w ∞ (s − t n ), we may bound the last term of (4.10) by w 2 ∞ ∆t 2 . Moreover, using again Young's inequality and the estimate |Z(s; t n , x) − x| w ∞ (s − t n ), we have, for any ε > 0, It implies, for any ε ∈ (0, 1), Thus we deduce that −2λ Injecting this latter inequality into (4.10) and taking ε = ∆t, we deduce Hence, since 2λ(1 − ∆t)∆t < 1, we have by induction, recalling that D 0 = 0, Using the assumption ∆t 1/2, we conclude the proof of Theorem 2.2 (2) in the case ρ ini = ρ 0 ∆x . 4th step. -We are left with the case ρ ini = ρ 0 ∆x . Let us define ρ (t) = Z (t) # ρ 0 ∆x , the exact solution with initial data ρ 0 ∆x . From the triangle inequality, we have The last term in the right hand side may be estimated thanks to the above computations. For the first term in the right hand side, we use the estimates in Theorem 2.1 (we apply (1) if λ 0 and (2) if λ > 0): where (λ) − = max(−λ, 0) is the negative part of λ. Let us define τ : σ)x, for x ∈ C J . We have that τ (0, · ) = id and τ (1, · ) # ρ ini = ρ 0 ∆x . Then (4.11) We deduce d W (ρ ini , ρ 0 ∆x ) ∆x. Then, we get

Unstructured mesh
We can extend our convergence result to more general meshes. For the sake of simplicity of the notation, we present the case of a triangular mesh in two dimensions. This approach can be easily extended to meshes made of simplices, in any dimension.

Forward semi-Lagrangian scheme
Let us consider a triangular mesh T = (T k ) k∈Z with nodes (x i ) i∈Z . We assume this mesh to be conformal: A summit cannot belong to an open edge of the grid. The triangles (T k ) k∈Z are assumed to satisfy k∈Z T k = R 2 and T k ∩ T l = ∅ if k = l (in particular, the cells are here not assumed to be closed nor open). For any triangle T with summits x, y, z, we will use also the notation (x, y, z) = T . We denote by V(T ) = V(x, y, z) the area of this triangle, and h(T ) its height (defined as the minimum of the three heights of the triangle T ). We make the assumption that the mesh satisfies := inf k∈Z h(T k ) > 0.
For any node x i , i ∈ Z, we denote by K(i) the set of indices indexing triangles that have x i as a summit, and we denote by T i the set of all triangles of T that have x i as a summit: thus For any triangle T k , k ∈ Z, we denote by the set of indices indexing the summits of T k (for some arbitrary order, whose choice has no importance for the sequel). We consider the following scheme, which may be seen as a forward semi-Lagrangian scheme on the triangular mesh.
• For an initial distribution ρ ini of the PDE (1.1), define the probability weights (ρ 0 i ) i∈Z through the following procedure: Consider the one-to-one mapping ι : Z k → ι(k) ∈ Z such that, for each k ∈ Z, x ι(k) is a node of the triangle T k ; ι is thus a way to associate a node with a cell; then, for all i ∈ Z, let ρ 0 i = k:ι(k)=i ρ ini (T k ). Observe from (4.11) that ρ 0 ∆x = j∈Z ρ 0 j δ x j is an approximation of ρ ini . • Assume that, for a given n ∈ N, we already have probability weights (ρ n i ) i∈Z such that ρ n ∆x = j∈Z ρ n j δ x j is an approximation of ρ(t n , · ), where ρ is the solution to (1.1) with ρ ini as initial condition. For i ∈ Z, we let , and y n i := x i + a n i ∆t.
Under the CFL-like condition (5.1) w ∞ ∆t , y n i belongs to one (and only one) of the elements of T i . We denote by k n i the index of this triangle, namely y n i ∈ T k n i . • We use a linear splitting rule between the summits of the triangle T k n i : the mass ρ n i is sent to the three points x I 1 (k n i ) , x I 2 (k n i ) , x I 3 (k n i ) according to the barycentric coordinates of y n i in the triangle.
Let us make more precise the latter point. Let T = (x, y, z) ∈ T , and ξ ∈ T . We define the barycentric coordinates of ξ with respect to x, y and z, λ T x , λ T y and λ T z : and then have ξ = λ T x (ξ)x+λ T y (ξ)y+λ T z (ξ)z. Note also that λ T x (ξ)+λ T y (ξ)+λ T z (ξ) = 1. Therefore, we have the following fundamental property, which will be used in the sequel: In the same spirit as in Section 4, we here define the interpolation weights by: For j ∈ Z, and y ∈ R 2 , Then, the numerical scheme reads (5.5) ρ n+1 j = i∈Z ρ n i α j (x i + a n i ∆t), j ∈ Z, n ∈ N.

Convergence result
By the same token as in Section 4, we can use Lemma 5.1 and Theorem 2.1 to prove that the numerical scheme (5.5) is of order 1/2: Theorem 5.2. -Assume that W satisfies hypotheses (A1)-(A4). For ρ ini ∈ P 2 (R d ), let (ρ(t)) t 0 be the unique measure solution to the aggregation equation with initial data ρ ini , as given by Theorem 2.1. Let us also consider a triangular conformal mesh (T k ) k∈Z with nodes (x j ) j∈Z such that = inf k∈Z h(T k ) > 0 and the CFL condition (5.1) holds true. We denote by ∆x the longest edge in the mesh.
Define ((ρ n j ) j∈Z ) n∈N as in (5.5) and let ρ n ∆x := j∈Z ρ n j δ x j , n ∈ N.
Then, there exists a nonnegative constant C, independent of the discretization parameters, such that, for all n ∈ N * , Importantly, we do not claim that (2) in the statement of Theorem 2.2 remains true in the framework of Theorem 5.2. Indeed, it would require to prove that the support of the numerical solution remains included in a ball when the support of the initial condition is bounded. As made clear by the proof of Lemma 3.3, this latter fact depends on the geometry of the mesh.

Numerical illustrations
We now address several numerical examples. In Subsection 6.2, we show that the rate of convergence established in Theorem 2.2 is optimal in a one-dimensional example. This prompts us to start with a short reminder on the Wasserstein distance in dimension d = 1. In Subsection 6.3, we provide several numerical examples in dimension d = 1 for the Newtonian potential, whilst examples in dimension d = 2 are handled in Subsection 6.4.

Wasserstein distance in one dimension
The numerical computation of the Wasserstein distance between two probability measures in any dimension is generally quite difficult. However, in dimension d = 1, there is an explicit expression of the Wasserstein distance and this allows for direct computations, including numerical purposes, as shown in the pioneering work [GT06]. Indeed, any probability measure µ on the real line R can be described thanks to its cumulative distribution function F (x) = µ((−∞, x]), which is a right-continuous and non-decreasing function with F (−∞) = 0 and F (+∞) = 1. Then we can define the generalized inverse Q µ of F (or monotone rearrangement of µ) by Q µ (z) = F −1 (z) := inf{x ∈ R : F (x) > z}; it is a right-continuous and non-decreasing function, defined on [0, 1). For every non-negative Borel-measurable map ξ : In particular, µ ∈ P 2 (R) if and only if Q µ ∈ L 2 ((0, 1)). Moreover, in the onedimensional setting, there exists a unique optimal transport plan realizing the minimum in (2.1). More precisely, if µ and ν belong to P p (R), with monotone rearrangements Q µ and Q ν , then Γ 0 (µ, ν) = {(Q µ , Q ν ) # L (0,1) } where L (0,1) is the restriction to (0, 1) of the Lebesgue measure. Then we have the explicit expression of the Wasserstein distance (see [RR98,Vil03]) , and the map µ → Q µ is an isometry between P 2 (R) and the convex subset of (essentially) non-decreasing functions of L 2 ([0, 1)). We will take advantage of this expression (6.1) of the Wasserstein distance in dimension 1 in our numerical simulations to estimate the numerical error of the upwind scheme (2.12). This scheme in dimension 1 on a Cartesian mesh reads, with time step ∆t and cell size ∆x: ∆t ∆x (a n j ) + ρ n j − (a n j+1 ) − ρ n j+1 − (a n j−1 ) + ρ n j−1 + (a n j ) − ρ n j .
With this scheme, we define the probability measure ρ n ∆x = j∈Z ρ n j δ x j . Then the generalized inverse of ρ n ∆x , denoted by Q n ∆x , is given by ρ n k .

Optimality of the order of convergence
Thanks to formula (6.1) in dimension d = 1, we can verify numerically the optimality of our result. Let us consider the potential W (x) = 2x 2 for |x| 1 and W (x) = 4|x| − 2 for |x| > 1; such a potential verifies our assumptions (A1)-(A4) with λ = 0. We choose the initial datum ρ ini = 1 2 δ −x 0 + 1 2 δ x 0 with x 0 = 0.25. Then the solution to the aggregation equation (1.1) is given by . Therefore, letting u n j := k j ρ n k for j ∈ Z, we can easily compute the error at time t n = n∆t by means of the two formulas (6.1)-(6.3): We then define the numerical error as e := max n T /∆t e n . We display in Figure 6.1 the numerical error with respect to the number of nodes in logarithmic scale, as computed with the above procedure (the time steps being chosen in a such a way that the ratio (2.15) in the CFL condition is kept constant). We observe that the computed numerical error is of order 1/2.

Newtonian potential in one dimension
An interesting and illustrative example is the Newtonian potential in dimension d = 1. Let us indeed consider the case W (x) = |x| and an initial datum given by the sum of two masses located at points x i 1 and x i 2 of the grid mesh, namely ρ ini = 1 2 δ x i 1 + 1 2 δ x i 2 , with say x i 1 < x i 2 . The solution of the aggregation equation in Theorem 2.1 is given by ρ(t) = 1 2 δ x 1 (t) + 1 2 δ x 2 (t) , where Indeed, recalling definition (2.2), we have, for t < x i 2 − x i 1 : At t = x i 2 − x i 1 , the two particles collapse, then for t x i 2 − x i 1 , we have ρ(t) = δ 1 2 (x i 1 +x i 2 ) . Standard finite volume upwind scheme. This simple example explains why we have chosen the scheme (6.2) instead of the standard finite volume upwind scheme introduced in Subsection 3.2. In dimension d = 1 and on a Cartesian grid, this latter one reads (6.4) ρ n+1 i = ρ n i − ∆t ∆x (a n i+1/2 ) + ρ n i − (a n i+1/2 ) − ρ n i+1 − (a n i−1/2 ) + ρ n i−1 + (a n i−1/2 ) − ρ n i , where a n i+1/2 = − k∈Z ρ n k sign(x i+1/2 − x k ). Assume indeed that, at time t n , for some n ∈ N, we have obtained the approximation ρ n i = 0 for i ∈ Z \ {i 1 , i 2 }, and ρ n i 1 = ρ n i 2 = 1/2. We then compute So, when applying the upwind scheme for i ∈ {i 1 − 1, i 1 , i 1 + 1}, we get Doing the same computation for i ∈ {i 2 − 1, i 2 , i 2 + 1}, we deduce that ρ n+1 = ρ n . Thus the above upwind scheme may be not able to capture the correct dynamics of Dirac deltas. The above computation is illustrated by the numerical results in Figure 6.2, where a comparison between the numerical results obtained with (6.4) (left) and with (6.2) (right) is displayed. We observe that the Dirac deltas are stationary when using the scheme (6.4), whereas the scheme (6.2) permits to catch the right dynamics. Another interesting numerical illustration of this phenomenon is provided by Figure 6.3. In this example, we choose the potential W (x) = 1 − e −2|x| , which is −4-convex, and a smooth initial datum given by the sum of two Gaussian functions: ρ ini (x) = 1 M (e −20(x−0.5) 2 + e −20(x+0.5) 2 ), where M = ρ ini L 1 is a normalization coefficient. With this choice, we observe that the solution blows-up quickly. Dirac deltas appear in finite time and, as observed above, the scheme (6.4) ( Fig. 6.3-left) does not allow to capture the dynamics after blow-up time, whilst the scheme (6.2) (Fig. 6.3-right) succeeds to do so. For these numerical simulations, the numerical spatial domain is [−1.25, 1.25]; it is discretized with a uniform Cartesian grid of 800 nodes, and the ratio in the CFL condition (2.15) is 1/2.

Comparison with Burgers-Hopf equation.
Considering the potential W (x) = 1 2 |x|, it has been proved in [JV16] (see also [BCDFP15]) that the following equivalence holds true: ρ is the solution in Theorem 2.1 if and only if u = −W * ρ is the entropy solution of the Burgers-Hopf equation ∂ t u + 1 2 ∂ x u 2 = 0. Let (ρ n i ) i∈Z,n∈N be given by the scheme (2.12)-(2.13). By conservation of the total mass, see Lemma 3.2, we have k∈Z ρ n k = 1. Introducing we deduce, by summing (2.12) and by using the fact that ρ n i = −(u n i − u n i−1 ), that the family (u n i ) i∈Z,n∈N satisfies: As in Fig. 6.2, the upwind scheme (6.4) does not capture the right dynamics of the Dirac deltas after blow-up time.
where, with (2.13), we have Then a n i = − Moreover, as ρ n i remains nonnegative under the CFL condition (see Lemma 3.1), Similarly, we get (a n i+1 ) − (u n i+1 − u n i ) = − a n i+1 (u n i+1 − u n i ) so that the scheme (6.5) for u finally rewrites Then we may apply the main result of this paper and deduce the convergence at order 1/2 of the above scheme: Lemma 6.1. -Let u ini be given in BV (R) such that ∂ x u ini 0 and T V (u ini ) = 1. Define the family (u n i ) i∈Z,n∈N by means of (6.6), with the initial data u 0 i := 1 2 + ∂ x u ini (−∞, x i+ 1 2 ), and let u n ∆x := i∈Z u n i 1 [x i ,x i+1 ) . Let u be the entropy solution to the Burgers equation ∂ t u + 1 2 ∂ x u 2 = 0 with u ini as initial condition. Then, there exists C 0, independent of the discretization parameters, such that if the CFL condition ∆t < ∆x is satisfied, one has u(t n ) − u n ∆x L 1 C( √ t n ∆x + ∆x).
Remark 6.2. -We do not claim that the scheme converges for any initial datum of the Cauchy problem for the Burgers equation (and actually it does not). The convergence result above only applies to a non-increasing initial condition belonging to [−1/2, 1/2].
Note that this scheme is not conservative, but, surprisingly (see [HLF94]) this does not prevent it from converging toward the right solution.
Proof. -First remark that the CFL condition that is here required is w ∞ ∆t < 1 2 ∆x, with w ∞ = 1/2 as W (x) = 1 2 |x|. The entropy solution u of the Burgers equation with a nonincreasing BV initial datum is a nonincreasing BV function. By Cauchy-Schwarz inequality, we have where (ρ(t)) t 0 is the solution of (1.1), with W (x) = 1 2 |x| as before and ρ ini = −∂ x u ini as initial condition, and (ρ n ∆x ) n 0 is the numerical solution obtained by Scheme (2.12) with d = 1 together with initial condition (2.11) (numerical solution whose convergence at order 1/2 is stated in Theorem 2.2).
The claim follows provided we prove that (6.7) R |u(t n , x) − u n ∆x (x)| dx = 1 0 |Q ρ(t n ) (z) − Q ρ n ∆x (z)| dz. In order to prove (6.7), we notice that, from a geometrical point of view, the left hand side of equality (6.7) corresponds to the area between the curves x → u(t n , x) and x → u n ∆x (x). Also, the right hand side is a measure of the area between their generalized inverses. However, the graph of the pseudo-inverse of a function may be obtained by flipping the graph of the function with respect to the diagonal. Since this operation conserves the area, we deduce that both areas are equal, that is (6.7) holds.
Another way to prove the identity (6.7) is to observe that the solution u of the Burgers-Hopf equation reads: where ρ is the solution in Theorem 2.1. In fact, as the number of points x for which ρ(t, {x}) > 0 is at most countable for any given t > 0, we have the almost everywhere equality: Similarly, So, to complete the proof, it suffices to use the fact that, for any two probability measures µ and µ on R, R |µ((x, +∞)) − µ ((x, +∞))|dx = 1 0 |Q µ (z) − Q µ (z)|dz, see [BL19, Theorems 2.9 and 2.10], noticing that the function Q µ we use here is the right continuous version of the quantile function used in [BL19].

Numerical simulation in two dimensions
As an illustration, we propose now a numerical example in two dimensions. The spatial domain is the square [0, 1] × [0, 1]; it is discretized with N x = 70 nodes in the x-direction and N y = 70 nodes in the y-direction; we take a time step ∆t = 10 −3 . We consider two different initial data: the sum of three bumps (as in [CJLV16]) where M is a normalization constant such that ρ ini L 1 = 1; and an initial density with a square shape With these numerical data, we compare the numerical results between the two potentials W 1 (x) = 1 − e −5|x| and W 2 (x) = 5|x|. For |x| close to 0, we have that ∇W 1 ∼ ∇W 2 . Then the short range interaction is similar between both potentials, but the long range interaction is different. The numerical results are displayed in Figures 6.4 and 6.6 for the potential W 1 (x) = 1 − e −5|x| and in Figures 6.5 and 6.7 for the potential W 2 (x) = 5|x|.
In each case, we observe, as expected, the aggregation in finite time of ρ towards a Dirac delta. Indeed it has been proved in [CDF + 11] that when the initial data is compactly supported, solutions converge towards a Dirac delta in finite time. We also observe that the time dynamics during this step of concentration is different between potentials W 1 and W 2 .
The case with an initial datum with three bumps has been implemented in [CJLV16] with a Lax-Friedrichs scheme. We obtain here similar results but we observe a smaller numerical diffusion. Then we can make similar comments for the comparison between the two potentials W 1 and W 2 . For the potential W 1 , we observe that each bump coalesces into a Dirac delta, then the three remaining Dirac deltas merge into a single Dirac delta (see Fig. 6.4). For the potential W 2 , the solution seems to be more regular and Dirac deltas seems to appear for larger time (see Fig. 6.5).
For the initial data with a square shape, the density ρ keeps, for both potentials, a shape similar to the initial square shape which tightens as time increases. However with the potential W 1 (Fig. 6.6), we notice a strong concentration at the corners of the square, whereas in the case of the potential W 2 (Fig. 6.7) the density is