Convergence of a linearly transformed particle method for aggregation equations

We study a linearly transformed particle method for the aggregation equation with smooth or singular interaction forces. For the smooth interaction forces, we provide convergence estimates in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^1$$\end{document}L1 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^\infty $$\end{document}L∞ norms depending on the regularity of the initial data. Moreover, we give convergence estimates in bounded Lipschitz distance for measure valued solutions. For singular interaction forces, we establish the convergence of the error between the approximated and exact flows up to the existence time of the solutions in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^1 \cap L^p$$\end{document}L1∩Lp norm.


Introduction
In this work, we are interested in showing the convergence of approximated particle schemes to the Cauchy problem for the so-called aggregation equation.This equation determines the evolution of a probability density ρ(t, x) defined by (1.1) Here, −∇W (x − y) measures the interaction force that an infinitesimal particle located at y ∈ R d will exert on a particle located at x ∈ R d .As a result, we will call W the interaction potential.Since the total mass is preserved, without loss of generality, we assume ˆRd ρ(t, x) dx = ˆRd ρ 0 (x) dx = 1 ∀t ≥ 0.
The microscopic dynamics of N particles X i , i = 1, . . ., N , interacting through the potential W are given by where the inertia term is assumed to be negligible compared to friction [45,46].The macroscopic dynamics (1.1) consists of a continuity equation where the velocity field is given by u(t, x) = −(∇W * ρ(t))(x), which is the mean-field limit of the microscopic system when N → ∞ under certain conditions on the potential [32,37,18,20].
We will focus the rest of the introduction on the well-posedness of the continuous equation (1.1) and the numerical methods proposed for its approximation.The equation (1.1) has the formal structure of being a gradient flow of a functional in the set of probability measures.Indeed, defining the interaction energy as for any probability measure µ, we find u = −∇ δF δµ where δF δµ is the formal variation of the functional F[µ].This observation leads to a natural formal Lyapunov functional for the solutions of equation (1.1).In fact, we expect solutions to satisfy the identity for all t ≥ 0. This structure can be rendered fully rigorous for C 1 -potentials [1] and it allows for mildly singular potentials at the origin [20,21,25] provided the interaction potential has some convexity property called λ-convexity.
On the other hand, global in time unique weak measure solutions can be constructed for any probability measure as initial data under suitable smoothness assumptions on the interaction potential.In this work, whenever we refer to smooth potentials, we mean that the interaction potential satisfies ∇W ∈ W 1,∞ (R d ).For smooth potentials, the approach introduced by Dobrushin for the Vlasov equation [32] using the bounded Lipschitz distance between probability measures, see [37,18,14] for further details, gives a well-posedness theory of weak measure solutions.
However, many of the interesting features related to blow-up dynamics and stationary states happen for potentials that are singular at the origin.Typical examples to bear in mind are combinations of repulsive attractive power-law potentials of the form W (x) = |x| a a − |x| b b with −d ≤ b < a and the convention |x| 0 0 = log |x|, or fully attractive potentials W (x) = |x| a a with a > −d, suitably cut-off at infinity.In this work, whenever we refer to singular potentials we mean that the interaction potential is not smooth but satisfies for some constant C > 0, and in addition we assume that ∇W is bounded away from the origin if α < 0. These conditions allow for singularities at the origin up to Newtonian but not including it.In particular, our singular potentials are such that ∇W ∈ W 1,q  loc (R d ) with a range depending on α: 1 ≤ q < d α+1 .Note that the power-law potentials satisfy locally the conditions of being a singular potential in the range 2 − d < b < 2 for repulsive-attractive and in the range 2 − d < a < 2 for fully attractive.The various well-posedness theories for measure solutions fail as soon as the potential becomes singular at the origin.However, weak solutions in Lebesgue spaces can be obtained.A local-in-time well-posedness theory was obtained in [10,18] for initial data in (L 1 ∩ L p )(R d ) with p = q the conjugate exponent of q, and in [7,9] a local-in-time well-posedness theory for initial data in (L 1 ∩ L ∞ )(R d ) was developed for singularities up to and including a Newtonian singularity at the origin, corresponding to α = d − 1.In this work, we will use the setting introduced in [18].The Newtonian case is very specific because of the relation between the divergence of the velocity field and the density becomes local.
Under the above assumptions of either smooth or singular potentials, the proofs of the global-in-time well-posedness of weak measure solutions and the local-in-time well-posedness of weak solutions for initial data in (L 1 ∩ L p )(R d ) spaces are essentially based on the fact that the velocity field is regular enough to have meaningful characteristics.It is proved in [32,37,10,18] that the velocity field of the constructed solutions is continuous in time and Lipschitz continuous in space.Then, the flow map Φ t (x), defined by the unique solution of the characteristic system is a diffeomorphism for all times t ≥ 0. In all cases, the solution built in [32,37,10,18] is obtained by characteristics and given by ρ(t) = Φ t #ρ 0 .Here, T #µ denotes the push-forward of a measure through a measurable map T : R A very interesting question is the rigorous derivation of the continuum description (1.1) starting from the microscopic dynamics (1.2) for both regular and singular potentials.This is the so-called mean-field limit problem.The mean-field limit results contain as a by-product convergence results for the classical particle method.More precisely, proving that (1.1) is the mean-field limit of the system (1.2) as N → ∞ is equivalent to show that the empirical measure converges weakly in measure sense to the solution of (1.1) provided that this weak convergence holds initially.Even if the particle method is proved to be convergent of order 1 N , the convergence error is only controlled in the bounded Lipschitz or Wasserstein-type distances between measures [32,37,18,20].
Vortex-blob methods, originally introduced for the 2D Euler equations for incompressible fluids, see [44] and the references therein, have also been adapted recently to the aggregation equation [8] with fixed shapes, where the approximate densities are shown to converge with arbitrary orders but only in negative Sobolev norms.Particle methods were also used in plasma physics for the Vlasov-Poisson system [30], where they are usually called smooth Particle In Cell (PIC) methods.
In the Linearly Transform Particle (LTP) method, introduced by Campos Pinto in [12] following an idea of Cohen and Perthame [29], particles are pushed on to discrete times according to an approximation of the exact flow as in standard particle methods.Moreover, particles have their own shape, which is transformed in the discrete evolution in order to better approach the local flow using a linearization of the exact flow.To our knowledge the LTP method has only been used for a linear transport equation [29] or for a Vlasov-Poisson system [13] involving measurepreserving characteristic flows.The technical difficulties posed by the deformation of the flows in our present case have been overcome by detailed estimates of the Jacobian matrices and determinants.These estimates have allowed us to control the error on the densities via the errors of the flows to finally obtain the convergence results.Certain Sobolev regularity is needed on the initial data to obtain convergence of the LTP method in Lebesgue spaces for both smooth and singular potentials.However, a general result of convergence for weak measure solutions is obtained in an appropriate distance for measures.
Let us finally mention that other numerical methods have been proposed in the literature for the aggregation equation.In [16], the authors proposed a finite volume scheme which is shown to be energy preserving, i.e., it keeps the property that the energy functional is dissipated along the semidiscrete flow.Finite volume and finite difference schemes have been shown to be convergent to weak measure solutions of the aggregation equations for mildly singular potentials in [41,25].
In this work, we extend the LTP method to the aggregation equation seen as one of the most important representatives of a class of nonlinear continuity equations with non divergence free velocity fields in any dimensions.We start by summarizing the basic ideas of the numerical LTP method in Section 2 together with the preliminaries and notations used in this work.Section 3 is devoted to give convergence results for smooth potentials in Lebesgue spaces.Depending on the regularity of the initial data, we will be able for smooth potentials to control errors in L 1 and L ∞ .For initial data just being a probability measure, we will show in Section 4 the convergence in bounded Lipschitz distance.In the case of singular potentials, we will control in Section 5 the error upto the existence time of the solution of (1.1) in L 1 and L p with p suitably chosen.We finally show in Section 6 the performance of this method in one dimension validating the numerical implementation with explicit solutions and making use of it to study certain not well-known qualitative features of the evolution of (1.1) with several smooth and singular potentials.

Basic properties of the exact flow.
In the setting of our main results, the velocity field of the exact solution to (1.1) is always continuous in t and Lipschitz continuous in x.The solution of the characteristic system is well-defined and it has unique global in time solutions for all initial data x ∈ R d .Moreover, the general solution of the characteristic system is a diffeomorphism in R d .The general flow map will be denoted by F s,t (x) for all t, s ∈ R and x ∈ R d .
The flow map satisfies (2.1) and the Jacobian matrix and its determinant satisfy the differential equations Using u(τ, y) = −(∇W * ρ(τ ))(y), this yields Estimates are then easily derived when u ∈ L ∞ (0, ∞; W 1,∞ (R d )).We will write L := sup t∈[0,∞) u(t, •) W 1,∞ .For instance, using (2.2) and J s,s (x) = I d we find (2.5) sup and in particular the characteristic flow is Lipschitz (relative to any norm in R d ), Furthermore, we derive from (2.3) and (2.5) that (2.7) sup and using (2.4) we also find Let us remark that the previous estimates (2.5)-(2.9)can also be obtained in a time interval [0, T ] for locally Lipschitz velocity fields u ∈ L ∞ (0, T ; W 1,∞ (R d )) for some T > 0, with constant L T := sup t∈[0,T ] u(t, •) W 1,∞ .These estimates will be used in Section 5, where the dependence on T of the Lipschitz constant will be omitted for the sake of simplicity.
2.2.Linearly Transformed Particles.As in standard particle methods, the density ρ is represented with weighted macro-particles, and as in smooth particle methods, particles have here a finite and smooth shape.Thus, we approximate the initial density ρ 0 on a Cartesian grid of size h > 0 by (2.10) with particle shapes obtained by scaling and translating a reference function, i.e., (2.11) Here the reference shape is assumed to have a compact support supp(ϕ) ⊂ B(0, R o ), be bounded and satisfy In this work we will require that the shape functions are Lipschitz, and we can either consider for the reference shape the tensor-product hat function or the B3-spline As for the weights ω k = ω k (h, ρ 0 ), they are usually defined as (2.14) d ρ 0 (x)dx, however this will not be sufficient to prove the convergence of our particle scheme without additional smoothness assumptions on the initial density ρ 0 .Indeed, using standard arguments (see e.g.[12,28]) based on the fact that the approximation ρ 0 → ρ 0 h = k∈Z d ω k ϕ 0 h,k is local, bounded in any L p space and preserves the affine functions, one easily verifies the following estimate.

Proposition 1. If ρ 0
h is initialized as in (2.10) with weights and shape function given by (2.14) and (2.11), respectively, then we have In our analysis we will need second-order estimates which are then available for ρ 0 ∈ W 2,p (R d ).However, if we allow negative weights then second-order estimates are also available in a dual norm, as follows.Consider weights defined as (2.16) with integration kernels bi-orthogonal to the shape functions in the sense that (2.17) holds with δ k,k the Kronecker symbol.Similar to the shape functions, they can be obtained by scaling and translating a reference φ (assumed again compactly supported, bounded and satisfying (2.2)) with a different normalization, namely For instance if ϕ is the above tensor-product hat function (2.12) then for the integration kernel we may take φ(x) = i≤d Notice that estimate (2.15) still holds with these weights.Now, from the duality (2.17) we can derive a convenient second-order estimate which only relies on the first-order smoothness of ρ 0 .It is expressed in the dual norm where q is the conjugate exponent of p and w, v is the duality pair that coincides with the integral of the product wv as soon as the latter is integrable.
Proof.It follows from the duality relation (2.17 h,k and standard arguments show that the approximation v → ṽh satisfies an error estimate similar to (2.15) for s = 1.Using the Hölder inequality this gives and the proof is completed due to the definition of the W −1,p (R d ) norm.
We observe that both the above initializations yield (2.20) sup and since the shape functions are assumed to be non-negative, (2.2) gives (2.21) with a constant depending only on φ.We now describe the LTP method.As mentioned in the introduction, compared to standard particle methods, the LTP method follows the shape of smooth particles.Therefore we need to track not only the particle positions but also their deformations given by the Jacobian matrices.Given discrete trajectories x n k approximating the exact ones F 0,tn (x 0 k ) on the discrete times t n := n∆t, n = 0, 1, . . ., N := T /∆t, and non singular approximations J n k of the forward Jacobian matrices J tn,tn+1 (x n k ), the particle shapes ϕ n+1 h,k are recursively defined as the push-forward of ϕ n h,k along the affine flow ) , which approximates the exact flow F tn,tn+1 around x n k .Here x n+1 k can also be seen as an approximation to F tn,tn+1 (x n k ), as will be specified below.In short, we define where j n k := det(J n k ) > 0. Starting from ϕ 0 h,k defined as in (2.11), this gives particles of the form (2.23) where the deformation matrix D n k and the particle volume h n k are defined by (2.24) It follows from the above process that D n k is an approximation to the backward Jacobian matrix J tn,0 (x n k ), whereas h n k approximates the elementary volume h d multiplied by the Jacobian determinant of the forward flow F 0,tn at x 0 k .Moreover, the particle shape ϕ n h,k is the push-forward of ϕ 0 h,k along the integrated flow which can be seen as a linearization of F 0,tn around x 0 k (for n = 0 we set Indeed, it follows from the above definitions that (2.26) , and we easily verify that Finally, the LTP approximation of the density at time t n is defined as with weights ω k constant in time and computed as in (2.14) or (2.16).According to (2.26), we have ´ϕn h,k = ´ϕ0 h,k = ´ϕ, and thus the conservation of mass ( ´ρn h = ´ρ0 h ) holds at the discrete level.Moreover, using the fact that the particle shapes are non-negative, we find as in (2.21) 2.3.Approximated Jacobian matrices and particle positions.To complete the description of the numerical method (2.23)-(2.24),(2.27), we are left to specify how to compute the particle center x n+1 k and the discrete Jacobian matrix J n k involved in the affine flow (2.22).Before doing so we observe that if the matrices D 2 W (x) and D 2 W (y) commute for all x and y, then the exact solution to the ODE (2.2) takes an exponential form.However, in the general case the matrix J tn,tn+1 (x) is not equal to but the difference is small, as shown next.
with a constant C independent of n ≤ N − 1 and ∆t.
Proof.Given n ≤ N − 1 and x ∈ R d , we denote for simplicity and we observe that |B(τ . From the above bound for B we readily find (a) where we have used (2.5) in the last inequality.The result follows.
Then we can define x n+1 k as the approximation given by a numerical scheme discretizing (2.30) when replacing the exact density ρ at discrete times in [t n , t n+1 ] by its LTP approximation ρ n h .In the convergence analysis, we consider particle trajectories x n k and approached Jacobian matrices J n k defined by an explicit Euler scheme: (2.31) Note that this expression can be seen as an approximation to (2.29) using a rectangular rule in the time integral (here will not take into account the approximation error of convolution products).Accordingly, we set (2.32) ) .Using (3.1) and the L 1 bound (2.28) on ρ n h , we see that this approximation yields Clearly, higher-order time discretizations are also possible.

It is easily verified that the difference between these approximations satisfies
as long as we have ∇W ∈ W 1,q (R d ) and sup 0≤n≤ T ∆t ρ n h L p ≤ C with p = q .2.4.General strategy of the convergence proofs.In order to establish error estimates for the approximation of the density ρ(t n ) by ρ n h we will use Gronwall arguments that involve errors on the flows and on the Jacobian determinants.Since the velocity fields depend nonlinearly on the densities, we need to couple these errors with the density approximation error, and since the k-th particle is pushed forward by the approximated flow F n h,k during the time interval [t n , t n+1 ], we need to control the local error between this approximation and the exact flow F tn,tn+1 .To this end we define a first error term on the support of the smooth particles, (2.33) e n F := sup In our analysis, we shall also need to track the error on an extended domain which accounts for the deformation of the particle support by the exact flow, namely The error corresponding to the integrated flow (2.25) is then defined as Using the fact that the exact flow is Lipschitz, see (2.6), it is easy to bound this term by accumulating the local flow errors, e ), but in the analysis we will need a finer control, see Lemma 4 below.We will also need to control the error of the Jacobian determinants for each particle, thus we define (2.35) e n j := sup .
Finally we will need to track carefully the particles that affect the local value of the approximated density.For this purpose, we let 3. L 1 and L ∞ convergence for smooth potentials In this section we assume that the potential is smooth, as defined in the introduction.This means that ∇W ∈ W 1,∞ (R d ).In this case, the Lipschitz norm of u is bounded by ∇W W 1,∞ : indeed letting |•| denote the Euclidean norm in R d as well as its associated matrix norm, we have for all and similarly for u, so that estimates (2.5)-(2.9)hold with L = C ∇W W 1,∞ .However, to obtain convergence rates in L p -spaces we need more regularity on the solutions.In turn we assume that ρ 0 ∈ W 1,1 + (R d ) in this section and we compute the weights with the formula (2.16) involving the dual kernels.According to the propagation of regularity of solutions to (1.1) in Proposition 9 in the Appendix, this ensures that the unique solution to (1.1) satisfies ρ ∈ L ∞ (0, T ; W 1,1 (R d )) for all T > 0.
Given the solution ρ to (1.1), we will use the shortcut notation, ρ n (x) := ρ(t n , x) for x ∈ R d .From now on, C denotes a generic constant independent of h and ∆t, depending only on L = sup t≤T |u(t)| W 1,∞ , d and the exact solution.
Moreover, we assume that both h and ∆t are bounded by an absolute constant.We denote by (3.2) the errors in L 1 and L ∞ norms.
3.1.Estimates on the flows and related terms.We first control the particle overlapping from the approximation error on the flow.
Lemma 1.There exists a constant C independent on h and ∆t such that From (2.25) we see that z k ∈ S 0 h,k .Using the Lipschitz bound (2.6) we then write , and the result follows.
Using the formulas (2.31), (2.32) and the a priori L 1 bound (2.28) on the approximated densities ρ n h we easily derive uniform estimates for the approximated Jacobian matrices and the particle supports.

Lemma 2. The approximated Jacobian determinants satisfy
∞ ∆t for a constant uniform in k and n ≤ N .In particular, J n k is always invertible and As for the deformation matrices We next show that the support of the particle approximation is of order h.
with constants C independent of ∆t and h.
To control the approximation errors for the velocity and the Jacobian matrices, we next introduce the generic error for τ ∈ [t n , t n+1 ], 0 ≤ n ≤ N − 1 and with a constant C independent of ∆t and h.
Proof.Using that u(τ, y) = −(∇W * ρ(τ ))(y) = −(∇W * (F 0,τ #ρ 0 ))(y), we write For the first term we write Using the expression (2.1) for the exact flow, estimate (3.7) and the equality ρ(s) L 1 = 1 gives then For (b), using the Lipschitz regularity of the flow (2.6) and the error bound (2.19) on the initial data we find Finally, we observe that for y ∈ S n h,l we have (F n h,l ) −1 (y) ∈ S 0 h,l from (2.25), and and arguing as in (2.28) this gives By gathering the above estimates, we complete the proof.

Proposition 5. If the initial density satisfies ρ
holds with a constant C depending only on d, T , L, and ρ 0 W 1,1 .Moreover, at x = x n k , we have sup For (c) we use the one-to-one mapping Φ : where we have used (2.8) in the first inequality.Using (2.9) this allows to bound Turning next to the (d) term, we introduce where in the last inequality we have used (see (2.1) and Lemma 3) To end the proof we will show that up to a sign and a translation, Ξ α is uniformly close to the identity mapping.Let G(y) := (F τ,tn − I)(F tn,τ (x) − y) so that Ξ α (y) = −y + αG(y) + (1 − α)x n k + αF tn,τ (x).From (2.7) we infer The desired bound |d| ≤ C(h + ∆t) follows by gathering the above steps.
We can now compute an estimate for the error of the Jacobian determinants.
where C is a positive constant depending only on T , L, and ρ L ∞ (0,T :W 1,1 ) .
Proof.According to (2.4) and (2.32), we have ) , so that Proposition 5 gives the desired result.
From Proposition 5 we also derive an estimate for the error between Jacobian matrices.
Proof.Using the matrix Jtn,tn+1 (x) defined by (2.29), Proposition 3 gives and to bound the remaining error we proceed as in the proof of Corollary 1: denoting )dτ , we use the exponential matrix expressions (2.31) and (2.29) to compute ) and using (2.28) this yields so that the desired result follows again from Proposition 5.
Remark 2. If ρ 0 is only assumed to be an L 1 (R d ) function (or a Radon measure), then ξ n (D 2 W ) can be bounded by a constant using the L 1 bound on ρ n h , see (2.28), and the W 1,∞ (R d ) smoothness of ∇W .Arguing as in the proof above we then find an error estimate for the Jacobian matrices on the order of ∆t.
We next turn to the approximation errors involving the forward characteristic flows and we establish a series of estimates.Proof.Given x ∈ S 0 k,h we write y = F 0,tn (x) and ỹk = F n h,k (x) ∈ S n h,k .We have with a constant C independent of ∆t and h.
Proof.Given x ∈ Sn h,k , we rewrite the linearized flow (2.22) as follows, . Using (2.31) and the expression (2.1) for the exact flow, we then compute where the inequality follows from (3.9) (note that here C depends on ρ 0 W 1,1 ).For (b), we easily get using estimate (3.11) in Corollary 2 and Lemma 3 that Turning to (c) we next differentiate (2.3) and obtain for 1 This yields where we used that ρ ∈ L ∞ (0, T ; for some C, see (2.5).Invoking the Gronwall Lemma, we then obtain where C only depends on d, T , L and ρ 0 W 1,1 .With a Taylor expansion this gives for some η n k between x and x n k and a constant C that only depends on d, T , L and ρ 0 W 1,1 .Combining the above estimates yields the desired result.
We finally provide estimates for e F n and ẽn F .
Proof.Using (3.12), (3.13) and the fact that e C∆t + C∆t ≤ e 2C∆t , we find follows by a summation using e F 0 = 0.The bound on ẽn F is obtained with (3.13).

Proof of L 1 and L
holds with a constant C depending only on d, T , L, and ρ 0 W 1,1 .
Proof.Let y ∈ R d .Using the relation ρ(t n ) = F tn,tn−1 #ρ(t n−1 ) and the form (2.27) of the approximate solution together with the fact that , we decompose the error ρ(t n , y) − ρ n h (y) into three parts as (3.14) F tn,tn−1 (y) j tn,tn−1 (y) .
Estimate of A n L 1 : Using the one-to-one change of variable x = F tn,tn−1 (y), we easily find that ˆRd Estimate of B n L 1 : By means of the same change of variable and the relation j tn,tn−1 (y) = (j tn−1,tn (x)) −1 , we obtain ˆRd |B n (y)|dy ≤ ˆRd due to (2.5), (2.28) and (2.35), indeed x can be taken in S n−1 h,k in the k-th term.Estimate of C n L 1 : Writing again x = F tn,tn−1 (y), we observe that in the kth term, we must consider the cases where y ∈ S n h,k and those where x ∈ S n−1 h,k .Thus, x must be taken in the extended particle support Sn−1 h,k , see (2.34).Using the incremental relation (2.24) we then estimate (2.33).To obtain a global bound we next observe that the measure of Sn−1 h,k is of order (h + ∆t) d ≤ Ch d according to Lemma 3 and the assumption ∆t ≤ Ch, as well as that of F tn−1,tn ( Sn−1 h,k ) according to (2.8).Using the above observations and the fact that the reference shape ϕ is assumed to be Lipschitz, we find where the last inequality follows from the uniform bounds on the matrices D n k and their determinants (Lemma 2), and from the estimates inside (2.28).

Conclusion:
We now combine all the estimates above and (3.10) in Corollary 1 to obtain Using Corollary 3 to estimate the flow error yields We next derive L ∞ -estimates.Here the required regularity propagates in time.As proved in the Appendix, Proposition 9, the unique solution to (1.1) belongs to holds with a constant independent of h and ∆t.
Proof.Given y ∈ R d , we decompose ρ(t n , y) − ρ n h (y) into three terms as in (3.14).
Estimate of A n L ∞ : Using the bound (2.8) on the exact Jacobian determinant, we find A n L ∞ ≤ e C∆t ε n−1 .Estimate of B n L ∞ : Writing again x = F tn,tn−1 (y), we observe that the k-th term vanishes if x ∈ S n−1 h,k .In particular, the sum can be restricted to the indices k in the set K n−1 (x).Gathering the bounds (3.4) on h n k , (2.20) on ω k and (3.3) on Estimate of C n L ∞ : Similarly as in the proof of Theorem 1, we observe that the k-th summand in C n (y) must be considered when y ∈ S n h,k or when x ∈ S n−1 h,k (or both).Clearly the cardinality of the corresponding index set satisfies Using the Lipschitz smoothness of the reference shape function ϕ as in (3.15), and again the bounds (3.4) on h n k , (2.20) on ω k and (3.3) on κ n , we write Conclusion: Combining the estimates above, we have Now, with the assumptions made here Theorem 1 applies, hence Corollaries 1 and 3 provide error estimates for the Jacobian and flow errors.Specifically, we have Plugging these estimates into (3.16)yields then due to ∆t h 1 and θ 0 ≤ 2. We again conclude with the discrete Gronwall Lemma.
Remark 3.Under the condition ∆t ≤ Ch in the convergence theorems, we have obtained the convergence estimates in L 1 and L ∞ with the terms of the form ∆t/h.We obviously need the assumption ∆t = o(h) to get the convergence results.

Convergence for measure solutions with smooth potentials
In this part, we consider measure valued solutions to the system (1.1) using the bounded Lipschitz distance.More precisely, let ρ 1 , ρ 2 ∈ M(R d ) be two Radon measures.Then the bounded Lipschitz distance d BL (ρ 1 , ρ 2 ) between ρ 1 and ρ 2 is given by Since the interaction potential W satisfies ∇W ∈ W 1,∞ (R d ), a well-posedness theory for measure valued solutions to (1.1) can be developed by using the classical results of Dobrushin [32], see [37,18] for related results.
To estimate the error between the exact flow and its local linearizations we now revisit some results from the previous section, namely Proposition 6, given the low regularity of the solutions.As in the previous section, we denote ρ n = ρ(t n ).Proposition 7. Let ρ 0 be an initial Radon measure on R d , and ρ n h be the approximation constructed in (2.27).If W satisfies ∇W ∈ W 1,∞ (R d ), then the flow error defined on the particles support (2.33) holds for 0 ≤ n ≤ N with a constant C independent of h and ∆t.
Proof.Let x ∈ S n h,k .We decompose the linearized flow as in Proposition 6, We next rewrite (a) = ´tn+1 tn ) dτ using (2.31) and (2.1), and estimate the integrand by From Using next a change of variable and the relation ρ(τ ) = F tn,τ #ρ n we get Combining the estimates above, we obtain For the estimate of (b), we easily get from Remark 2 that |(b)| ≤ Ch∆t.Finally, we observe that (c) cannot be estimated as in the proof of Proposition 6, due to the lesser regularity of the densities.We then proceed as follows, where we used estimate (3.6) for x ∈ S n h,k , and the estimates (2.6) and (2.7).Theorem 3. Let ρ 0 be an initial probability measure on R d , and ρ n h be the approximation constructed in (2.27).Assume that the interaction potential W satisfies holds, where C depends only on d and L.

Remark 4.
Observe that a convergence condition on the approximation of the initial data in Theorem 3 such as d BL (ρ 0 , ρ 0 h ) h is easily achieved by using a uniform quadrangular mesh of size h d and approximating the initial data ρ 0 by a sum of Dirac deltas via transporting the mass of ρ 0 inside each d-dimensional cube to its center.A cut-off procedure to leave small mass outside a large ball allows us to reduce to a finite number of Dirac deltas in this approximation.Finally, the error produced between smoothed particles and Dirac deltas is obviously of order h in the d BL distance.

Proof of Theorem 3. Since ρ
)e L∆t and we estimate (b) with where the last inequality uses the estimates inside (2.28).This leads to and using Lemma 7 we obtain with constants independent of ∆t and h.The proof is then completed using Gronwall's inequality as in Theorem 1.

L 1 and L p convergence for singular potentials
In this part, we are interested in L p -convergence between the solution and its approximation allowing for more singular potentials.With this aim, we consider the solutions of the equation ( 1 to be determined depending on the singularity of the potential.Since we are dealing with both attractive and repulsive potentials, we can only expect local in time existence and uniqueness of solutions as in [10,18].In those references, a local in time well-posedness theory in L ∞ (0, T ; L 1 (R d ) ∩ L p (R d )) was developed under suitable assumptions on the potentials.The solutions are constructed by characteristics since the velocity fields are still Lipschitz continuous in x.However, to prove convergence rates we need more regularity on the solutions.For the existence of solutions to (1.1) , we provide a priori estimates in Appendix A, Proposition 10.These estimates combined with the existing literature [18,10] show the well-posedness of solutions in the desired class.In our presentation we will follow the setting of local existence introduced in [18].
Let us remind the set of hypotheses on the interaction potential called singular potentials in the introduction.We assume that there exists L > 0 such that (5.1) and for −1 ≤ α < 0 and In particular, singular potentials satisfy ∇W ∈ W 1,q loc (R d ) for all 1 ≤ q < d α+1 .Note that (5.1) implies (see [18,39]) We remind the reader that these assumptions are enough to guarantee that the velocity fields are bounded and Lipschitz continuous with respect to x locally in time for densities in (L 1 ∩ L p )(R d ) where p is the conjugate exponent of q.Note that q = p < d α+1 is equivalent to α < −1 + d p , giving us the condition on the initial data for the well-posedness theory.Indeed, it follows from (5.1) that for some constant C depending on L, q and d, and a similar estimate holds for u using (5.2) and the fact that ∇W is bounded away from the origin.
Let T * be the maximal time of existence of weak solutions ρ ∈ L ∞ (0, T ; (L 1 ∩ L p )(R d )) with T < T * constructed in [18].Additional regularity will be needed on these solutions ensured by Proposition 10 of Appendix A under suitable initial data assumptions.In this section we consider T < T * , and we denote again t n = n∆t with 0 ≤ n ≤ N and ∆t = T /N for some given positive integer N .We introduce the following notations: As for the convergence analysis, we point out that the proof of Section 3 cannot be directly applied.Indeed, it is not obvious to obtain an a piori bound on sup 0≤n≤N ρ n h L p uniformly in h and ∆t, which we need to estimate (∇W * ρ n h ) and (D 2 W * ρ n h ).In order to do that, we will prove by induction that there is some h * > 0 for which We remind the reader that our error analysis between exact and approximated solutions for singular potentials requires non-negative weights for the particles, and this imposes us to give higher regularity on the initial data ρ 0 ∈ W 2,p (R d ), see Proposition 1.Using the results in [18] and Appendix A, we can obtain the existence and uniqueness of a solution ρ ∈ L ∞ (0, T ; (L 1 ∩ W 2,p )(R d )).However, in the next results we need less regularity in the solutions than on the initial data.Therefore, we prefer to keep both the assumptions stating the needed properties on the solution ρ and the initial data ρ 0 to emphasize this fact.Under the (induction) assumption that Γ n h is bounded uniformly in h and ∆t, we can derive the following estimates.Lemma 5.If M > 0 and n ≤ N are such that Γ n h ≤ M , and if the solution to with a constant C M depending on M but not on h and ∆t.
Proof.A straightforward computation yields In a similar way to (5.4), we also bound products like , from which we derive estimates similar to those of Lemma 2. In particular, following the proof of Lemma 3 we find that for x ∈ S m h,k , ( and for x ∈ Sm h,k we find |x − x m k | ≤ C M (h + ∆t).Note that this latter estimate involves bounding (5.5) on S m+1 h,k which only requires the norm ρ l h for l ≤ m, so that the resulting estimate indeed involves a constant depending on M .
We next give the estimates of u(τ, F tm,τ ) − u m k for τ ∈ [t m , t m+1 ] and ξm (D 2 W ) for 0 ≤ m ≤ n − 1.The proof can be obtained by using similar arguments as in Proposition 5 with the help of Lemma 5 and a second-order estimate provided either by Proposition 2 or by a standard L p error estimate as described in Proposition 1.We omit its proof, but point out that the crucial point is the smoothness assumptions (5.3) on the singular potential and the Lipschitz bound (5.4) on the velocity field.
We can also adapt the proof of Corollary 1, Lemma 6, and Proposition 6 to obtain the following result.
with constants C M depending on M but not on h and ∆t.
We finally connect the errors to the L 1 ∩ L p bounds on the densities.
Using this together with (5.6) completes the proof.
We are now in a position to show the uniform L 1 ∩ L p bounds on the density.
Proposition 8. Assume that the interaction potential W is singular in the sense of (5.1)-(5.2),and let ρ be a solution to the equation Then for all M > 0, there exists h * (M ) > 0 such that Proof.We use an induction argument on n.Since Γ 0 h = Γ 0 h h 2 , clearly there exists h 0 (M ) such that Γ 0 h ≤ M for all h < h 0 (M ).We then assume that n < N and h n (M ) > 0 are such that sup 0<h≤hn(M ) For the remaining of the proof we then consider m ≤ n and h ≤ h n (M ).In particular, we observe that the Lemmas above can be used with this value of M .Decomposing the error as in Theorem 1, we write .
Using arguments similar than in Theorem 1 we find For the estimate of C m+1 (y), we use the interpolation inequality and the estimates in Theorems 1 and 2 to get Using Lemma 8 and the fact that Γ m h ≤ M and ∆t h we find that both e F m and e F m+1 are bounded by C M h, thus and the above estimates yield We also observe that in the proof of Theorem 1, all the steps leading to the estimate (where we remind that θ m = ρ m − ρ m h L 1 ) are valid in the case of singular potentials.This yields On the other hand, it follows from Lemmas 7 and 8 that and ẽm where we used the assumption ∆t h 2 .Thus we find Since this is valid for all m ≤ n, it follows from Gronwall's lemma that Γ n+1 h ≤ C M h holds for some constant C M > 0. We remind the reader that C M is the generic constant depending on M but independent of h and ∆t.In particular, setting h n+1 (M ) := min(h n (M ), M/C M ) allows to write sup This ends the induction argument and the proof, by taking h * (M ) = h N (M ).
Putting together all the results in this section, we obtain the main convergence result in (L 1 ∩ L p )(R d ).
Theorem 4. Assume that the interaction potential W is singular in the sense of (5.1)-(5.2),and let ρ be a solution to the equation (1.1) given by Proposition 8 and a constant C independent of h and ∆t.

Numerical Results
We will present in this Section some numerical examples in one dimension, with different interaction potentials and initial densities to showcase some of the features already observed in numerical and theoretical analysis of the aggregation equation (1.1) in [35,36,6,40,9,3].In this way, we first validate our numerical implementation in order to explore some less-known properties about the behavior of its solutions in one dimension.A further more complete numerical study in 2D of this method will be reported elsewhere.These examples already show the wide range of different behaviors of solutions to the aggregation equation.
6.1.Numerical method: validation and implementation.We have implemented the numerical method described in Section 2.2 using Python.We use different initial conditions depending on the behaviors we would like to show.Specifically, we consider as initial densities (6.1) in order to have asymmetric, discontinuous symmetric and compactly supported smooth initial data respectively.These initial densities have been normalized to have unit mass.Shape functions for the particle method are here B3-splines given by (2.13).We first examine the validation of our code by comparison of the numerical solution and the exact solution of (1.1) with W (x) = x 2 .Due to the conservation of the center of mass, ∀t ≥ 0, ˆR xρ(t, x)dx = ˆR xρ 0 (x)dx := λ , the solution is explicitly given by (6.4) ρ(t, x) = ρ 0 (x − λ)e 2t + λ e 2t , using the method of characteristics.Figure 2 (left) shows the exact solution of (1.1) with initial data (6.1) and the numerical solution computed with the LTP method, together with the L 1 and L ∞ errors with respect to h.
Let us now compare the results with classical particle methods.One of the drawbacks of classical particle methods in which the density is reconstructed with shape function of same size is the need to choose adequate values of ε.Indeed, if ε is too small compared to the distance between two particles, the reconstructed density will vanish between particles and is thus irrelevant; and if ε is too large the reconstructed density will be too spread out and the results lack accuracy, as it is demonstrated in Figure 2 (right).Figure 3 presents the L 1 and L ∞ errors between a standard Smooth Particle (SP) method (with different values of ε) and our LTP method for the potential W (x) = x 2 for which solutions are explicit by the method of characteristics, see (6.4).On the left picture, we observe that the optimal ε, at this instant t, for a classical particle method is well captured with the LTP method.In both cases, W (x) = x 2 and ρ 0 is given by (6.1) with errors computed at t = 0.5.
One could object that the gain is not significantly better with the LTP method.However, since particles aggregate, the average distance between two particles decreases exponentially in time, and consequently the optimal size ε for reconstruction in classical particle method is not the same during the whole simulation.Therefore, an evolution in time of ε is much better adapted.Notice that the case of potential W (x) = x 2 is not particularly the best example to show the higher accuracy of the LTP method with respect to the classical particle method since all particles have the same size at time t because j 0,t ex (x) = e −2t , see (6.4).Moreover, the gain of accuracy with the LTP method is even clearer while using B1-spline as shape functions instead of B3-spline, as it is demonstrated on Figure 4. 6.2.Numerical Simulations.We now take advantage of the method to explore the behavior for other attractive potentials of type W (x) = |x| a a , a > 1.Notice that for a ≥ 2 the potential is smooth while for 1 < a < 2 is singular once W is cut-off at infinity or if the initial data is compactly supported since the effective values of the potential lie on a bounded set and W can be cut-off at infinity without changing the solution.Figure 5 presents the numerical results obtained by the LTP method in the case of a = 1.5 and a = 2.5.We represent the approximation of the density ρ n h , and also the reconstructed velocity u n h and the reconstructed size of particles h n by piecewise linear interpolation such that Potentials and their derivatives are also represented.In both cases, we observe that the density converges to a Dirac mass. Figure 5 also shows that for a = 2.5, W ∈ L ∞ loc , no finite-time blow-up in L ∞ appears, opposite to the case a = 1.5 in agreement with the results proved in [7].Notice also the different qualitative behavior in their trend to blow-up as studied in [40].
Notice again that for b ≥ 2 the potential is smooth while for 1 < b < 2 is singular once W is cut-off at infinity or if the initial data is compactly supported as discussed above.Figure 6 presents the approximated density ρ n h , reconstructed velocity u n h and size of particles h n obtained by the LTP method in the case of the attractive-repulsive potentials with (a, b) = (3, 1.5) and (a, b) = (3, 2.5).In this case ρ 0 is given by (6.2).
We observe that the long time asymptotics for b = 2.5 are characterized by the concentration of mass equally onto Dirac deltas at two points in infinite time, while for b = 1.5 we obtain a convergence in time towards a steady L 1 density profile seemingly diverging at the boundary of the support.This last behavior has been reported in several simulations and related problems [6].However, it has not been rigorously proven yet.Let us point out that the set of stationary states when the interaction potential is analytic in 1D consists of a finite number of Dirac deltas as proven in [35,36].This result also holds for W (x) = |x| a a − |x| b b , 2 < b < a, as it will be reported in [23].
Figure 7 also represents the time evolution of the approximated density for (a, b) = (3, 2.5), with ρ 0 given by (6.3).Solutions in the range 2 < b < a for initial data in L 1 ∩ L ∞ exist globally in time, see [37].The numerical evidence shows that all solutions converge towards stationary states consisting of finite number of Dirac Deltas as t → ∞ in this range.
Finally, we show in Figure 8 the results of the stationary state of the SP method versus the LTP method for the potential (a, b) = (4, 2.5) with N = 100.We observe how the good local adaption of the size of the particles makes our approximation much better with no oscillations with respect to the SP method showing the good performance of the LTP method in this case and its good properties at work.As mentioned in the introduction, vortex-blob type methods have been shown to converge for the aggregation equation (1.1).They obtained convergence estimates in suitable L p norms for the velocity fields and the associated characteristics fields  while the error for the densities was controlled in suitable W −1,p -norms in [8,Th. 3.8].The error estimates for vortex-blob and SP methods depend as usual on the regularization of particles and the fixed particle size related in a suitable way to get convergence.We have proven that the LTP method has in contrast direct error estimates for the densities in L p depending on the initial mesh size showing that the local adaptation of the shape has this benefit on the error estimates too.

Appendix A. A priori estimates on the regularity of solutions
In this part, we deduce a priori estimates on the regularity of equation (1.1) that combined with the global/local in time well posedness theory obtained in [32,37,10,18], leads to the existence of solutions with the desired properties to apply the convergence results of previous sections.
As we remind the reader in the introduction and in several places along the text, there are two different well-posedness settings: for smooth and for singular potentials.In both cases under the assumptions on the initial data the velocity fields are continuous in time and Lipschitiz continuous in space.In the smooth  potential case, this property holds globally in time leading to unique global in time measure solutions [32,37].In the singular potential case, this property holds locally in time only since there exist blowing-up of solutions for fully atractive potentials, see [7,10,18].In both cases, the flow map associated to the velocity field is welldefined and solutions are obtained by pushing forward the initial data through the flow map.
In this Appendix, we present first a global-in-time propagation of regularity result in the smooth potential case adapted to our hypotheses on the convergence results.On the other hand, we show a local-in-time propagation of regularity result in the singular potential case.Proposition 9. Assume that the interaction potential W satisfies ∇W ∈ W 1,∞ (R d ).
Since u ∈ W 1,∞ (R d ), we obtain This completes the proof.
This implies where we used u ∈ W 1,∞ (R d ) and the estimate (A.4).On the other hand, D 2 u(s, •) L ∞ can be estimated by Hence, we have and by applying Gronwall's inequality to conclude the desired result.Similar arguments were used in [2] to construct classical solutions.

Lemma 4 .
For 0 ≤ n ≤ N − 1, the following estimate holds (3.12) e F n+1 ≤ e C∆t e F n + ẽn F with a constant C independent of ∆t and h.

Lemma 7 .
If M > 0 and n ≤ N are such that Γ n h ≤ M , and if the solution ρ with constants C M depending on M but not on h and ∆t.Proof.Since Lemma 4 only relies on the Lipschitz smoothness of the exact flow, we have e F m+1 ≤ e C∆t e F m + ẽm F for all m.Then from (5.6) we derive e F m+1 ≤ e (C+C M )∆t e F m + C M ∆t h 2 + ∆t + hΓ m h for m ≤ n, so that Gronwall's inequality (together with Γ m h = max m ≤m Γ m h ) yields

Figure 2 .
Figure 2. Comparisons for W (x) = x 2 , at t = 0.5 with initial data (6.1) and ∆t = 10 −4 .Left Figure: Comparison between the exact solution ρ(t, x) given by (6.4) and ρ n h computed by LTP method for various h.Right Figure: Approximated values ρ n h obtained with a classical Smooth Particle (SP) method for h = 0.01 and different fixed particle sizes ε.

Figure 3 .
Figure 3. Left Figure: Log-Log Plot of the L 1 and L ∞ errors of the SP method using various values for the particle radius ε, versus those of the LTP method (both using h = 0.01 and ∆t = 0.01).Right Figure: Log-Log Plot of the L 1 and L ∞ errors for the SP and LTP methods with ∆t = 0.0001 and different values of h = ε.In both cases, W (x) = x 2 and ρ 0 is given by (6.1) with errors computed at t = 0.5.

Figure 4 .
Figure 4. Comparisons between exact solution (at t = 0.5) and approximated solutions ρ n h with LTP method and SP method, with W (x) = x 2 /2, ρ 0 given by (6.1), h = 1/25 and ∆t = 10 −3 .On the left Figure, the shape function ϕ is a hat function, whereas on the right Figure, ϕ is a B3-spline.

Figure 5 .
Figure 5. Approximated density and reconstructed velocity and size of particles computed by the LTP method with h = 0.01 for W (x) = |x| a a , with a = 1.5 or a = 2.5 and ρ 0 given by (6.2) with the number of time-steps N = 200.

Figure 6 .
Figure 6.Approximated density and reconstructed velocity and size of particles computed by the LTP method with h = 0.01 for W (x) = |x| a a − |x| b b , with a = 3 and b = 1.5 or b = 2.5 and ρ 0 given by (6.2) with the number of time-steps N = 200.