A blob method for diffusion

As a counterpoint to classical stochastic particle methods for diffusion, we develop a deterministic particle method for linear and nonlinear diffusion. At first glance, deterministic particle methods are incompatible with diffusive partial differential equations since initial data given by sums of Dirac masses would be smoothed instantaneously: particles do not remain particles. Inspired by classical vortex blob methods, we introduce a nonlocal regularization of our velocity field that ensures particles do remain particles, and we apply this to develop a numerical blob method for a range of diffusive partial differential equations of Wasserstein gradient flow type, including the heat equation, the porous medium equation, the Fokker-Planck equation, the Keller-Segel equation, and its variants. Our choice of regularization is guided by the Wasserstein gradient flow structure, and the corresponding energy has a novel form, combining aspects of the well-known interaction and potential energies. In the presence of a confining drift or interaction potential, we prove that minimizers of the regularized energy exist and, as the regularization is removed, converge to the minimizers of the unregularized energy. We then restrict our attention to nonlinear diffusion of porous medium type with at least quadratic exponent. Under sufficient regularity assumptions, we prove that gradient flows of the regularized energies converge to solutions of the porous medium equation. As a corollary, we obtain convergence of our numerical blob method, again under sufficient regularity assumptions. We conclude by considering a range of numerical examples to demonstrate our method's rate of convergence to exact solutions and to illustrate key qualitative properties preserved by the method, including asymptotic behavior of the Fokker-Planck equation and critical mass of the two-dimensional Keller-Segel equation.


Introduction
For a range of partial differential equations, from the heat and porous medium equations to the Fokker-Planck and Keller-Segel equations, solutions can be characterized as gradient flows with respect to the quadratic Wasserstein distance. In particular, solutions of the equation where ρ is a curve in the space of probability measures, are formally Wasserstein gradient flows of the energy where L d is d-dimensional Lebesgue measure. This implies that solutions ρ(t, x) of (1) satisfy for a generalized notion of gradient ∇ W 2 , which is formally given by where δE/δρ is the first variation density of E at ρ (c.f. [3,27,28,81]). Over the past twenty years, the Wasserstein gradient flow perspective has led to several new theoretical results, including asymptotic behavior of solutions of nonlinear diffusion and aggregationdiffusion equations [27,28,68], stability of steady states of the Keller-Segel equation [9,11], and uniqueness of bounded solutions [26]. The underlying gradient flow theory has been well developed in the case of convex (or, more generally, semiconvex) energies [2-4, 24, 54, 76, 81, 82], and more recently, is being extended to consider energies with more general moduli of convexity [5,26,28,35].
Wasserstein gradient flow theory has also inspired new numerical methods, with a common goal of maintaining the gradient flow structure at the discrete level, albeit in different ways. Recent work has considered finite volume, finite element, and discontinuous Galerkin methods [8,16,21,60,79]. Such methods are energy decreasing, positivity preserving, and mass conserving at the semidiscrete level, leading to high-order approximations. They naturally preserve stationary states, since dissipation of the free energy provides inherent stability, and often also capture the rate of asymptotic decay. Another common strategy for preserving the gradient flow structure at the discrete level is to leverage the discrete-time variational scheme introduced by Jordan, Kinderlehrer, and Otto [54]. A wide variety of strategies have been developed for this approach: working with different discretizations of the space of Lagrangian maps [42,55,[65][66][67], using alternative formulations of the variational structure [43], making use of convex analysis and computational geometry to solve the optimality conditions [7], and many others [10,17,23,29,31,46,47,83].
In this work, we develop a deterministic particle method for Wasserstein gradient flows. The simplest implementation of a particle method for equation (1), in the absence of diffusion, begins by first discretizing the initial datum ρ 0 as a finite sum of N Dirac masses, that is, where δ x i is a Dirac mass centered at x i ∈ R d . Without diffusion and provided sufficient regularity of V and W , the solution ρ N of (1) with initial datum ρ N 0 remains a sum of Dirac masses at all times t, so that and solving the partial differential equation (1) reduces to solving a system of ordinary differential equations for the locations of the Dirac masses, ∇W (x i − x j )m j , i ∈ {1, . . . , N }.
The particle solution ρ N (t) is the Wasserstein gradient flow of the energy (2) with initial data ρ N 0 , so in particular the energy decreases in time along this spatially discrete solution. The ODE system (5) can be solved using range of fast numerical methods, and the resulting discretized solution ρ N (t) can be interpolated in a variety of ways for graphical visualization.
This simple particle method converges to exact solutions of equation (1) under suitable assumptions on V and W , as has been shown in the rigorous derivation of this equation as the mean-field limit of particle systems [22,24,51]. Recent work, aimed at capturing competing effects in repulsive-attractive systems and developing methods with higher-order accuracy, has considered enhancements of standard particle methods inspired by techniques from classical fluid dynamics, including vortex blob methods and linearly transformed particle methods [19,36,45,48]. Bertozzi and the second author's blob method for the aggregation equation obtained improved rates of convergence to exact solutions for singular interaction potentials W by convolving W with a mollifier ϕ ε . In terms of the Wasserstein gradient flow perspective this translates into regularizing the interaction energy (1/2) (W * ρ) dρ as (1/2) (W * ϕ ε * ρ) dρ.
When diffusion is present in equation (1), the fundamental assumption underlying basic particle methods breaks down: particles do not remain particles, or in other words, the solution of (1) with initial datum (3) is not of the form (4). A natural way to circumvent this difficulty, at least in the case of linear diffusion (m = 1), is to consider a stochastic particle method, in which the particles evolve via Brownian motion. Such approaches were originally developed in the classical fluids case [33], and several recent works have considered analogous methods for equations of Wasserstein gradient flow type, including the Keller-Segel equation [49,51,52,61]. The main practical disadvantage of these stochastic methods is that their results must be averaged over a large number of runs to compensate for the inherent randomness of the approximation. Furthermore, to the authors' knowledge, such methods have not been extended to the case of degenerate diffusion m > 1.
Alternatives to stochastic methods have been explored for similar equations, motivated by particle-in-cell methods in classical fluid, kinetic, and plasma physics equations. These alternatives proceed by introducing a suitable regularization of the flux of the continuity equation [34,74]. Degond and Mustieles considered the case of linear diffusion (m = 1) by interpreting the Laplacian as induced by a velocity field v, ∆ρ = ∇ · (vρ), v = ∇ρ/ρ, and regularizing the numerator and denominator separately by convolution with a mollifier [40,73]. For this regularized equation, particles do remain particles, and a standard particle method can be applied. Well-posedness of the resulting system of ordinary differential equations and a priori estimates relevant to the method were studied by Lacombe and Mas-Gallic [57] and extended to the case of the porous medium equation by Lions and Mas-Gallic [59,62]. In the case m = 2 on bounded domains, Lions and Mas-Gallic succeeded in showing that solutions to the regularized equation converge to solutions of the unregularized equation, as long as the initial data has uniformly bounded entropy. Unfortunately, this assumption fails to hold when the initial datum is given by a particle approximation (3), and consequently Lions and Mas-Gallic's result doesn't guarantee convergence of the particle method. An alternative approach, now known as the particle strength exchange method, incorporates instead the effects of diffusion by allowing the weights of the particles m i to vary in time. Degond and Mas-Gallic developed such a method for linear diffusion (m = 1) and proved second order convergence with respect to the initial particle spacing [38,39]. The main disadvantage of these existing deterministic particle methods is that, with the exception of Lions and MasGallic's work when m = 2, they do not preserve the gradient flow structure [59]. Other approaches that respect the method's variational structure have been recently proposed in one dimension by approximating particles by non-overlapping blobs [25,30]. For further background on deterministic particle methods, we refer the reader to Chertock's comprehensive review [32].
The goal of the present paper is to introduce a new deterministic particle method for equations of the form (1), with linear and nonlinear diffusion (m ≥ 1), that respects the problem's underlying gradient flow structure and naturally extends to all dimensions. In contrast to the above described work, which began by regularizing the flux of the continuity equation, we follow an approach analogous to Bertozzi and the second author's blob method for the aggregation equation and regularize the associated internal energy F. For a mollifier ϕ ε (x) = ϕ(x/ε)/ε d , x ∈ R d , ε > 0, we define For more general nonlinear diffusion, we define F ε (ρ) = F (ϕ ε * ρ) dρ, F : (0, ∞) → R.
Under sufficient regularity conditions, we prove that solutions of the regularized gradient flows converge to solutions of equation (1); see Theorem 5.8. When m = 2 and the initial datum has bounded entropy, we show that these regularity conditions automatically hold, thus generalizing Lions and Mas-Gallic's result for the porous medium equation on bounded domains to the full equation (1) on all of R d ; see Corollary 5.9 and [59, Theorem 2]. For this regularized equation (8), particles do remain particles; see Corollary 5.5. Consequently, our numerical blob method for diffusion consists of taking a particle approximation for (8). We conclude by showing that, under sufficient regularity conditions, our blob method's particle solutions converge to exact solutions of (1); see Theorem 6.1. We then give several numerical examples illustrating the rate of convergence of our method and its qualitative properties.
A key advantage of our approach is that, by regularizing the energy functional and not the flux, we preserve the problem's gradient flow structure. Still, at first glance, our regularization of the energy (6) may seem less natural than other potential choices. For example, one could instead consider the following more symmetric regularization Although studying the above regularization is not without interest, we focus our attention on the regularization in (6) and (7) for numerical reasons. Indeed, computing the first variation density of U ε gives δU ε δρ = ϕ ε * (U • (ϕ ε * ρ)), for F ε . In the first case, one can see that replacing ρ by a sum of Dirac masses still requires the computation of an integral convolution with ϕ ε , whereas in the second case, all the convolutions reduce to finite sums, which are numerically less costly. Another advantage of our approach, in the m = 2 case, is that our regularization of the energy can naturally be interpreted as an approximation of the porous medium equation by a very localized nonlocal interaction potential. In this way, our proof of the convergence of the associated particle method provides a theoretical underpinning to approximations of this kind in the computational math and swarming literature [56,58]. Further advantages our blob method include the ease with which it may be combined with particle methods for interaction and drift potentials, its simplicity in any dimension, and the good numerical performance we observe for a wide choice of interaction and drift potentials.
Our paper is organized as follows. In Section 2, we collect preliminary results concerning the regularization of measures via convolution with a mollifier, including a mollifier exchange lemma (Lemma 2.2), and relevant background on Wasserstein gradient flow and weak convergence of measures. In Section 3, we prove several results on the general regularized energies (7), which are of a novel form from the perspective of Wasserstein gradient flow theory, combining aspects of the well-known interaction and internal energies. We show that these regularized energies are semiconvex and differentiable in the Wasserstein metric and characterize their subdifferential with respect to this structure; see Propositions 3.10-3.12. In Section 4, we prove that F ε Γ-converges to F as ε → 0 and that minimizers converge to minimizers, when in the presence of a confining drift or interaction term; see Theorems 4.1 and 4.5. With this Γ-convergence in hand, in Section 5 we then turn to the question of convergence of gradient flows, restricting to the case m ≥ 2. Using the framework introduced by Sandier and Serfaty [75,77], we prove that, under sufficient regularity assumptions, gradient flows of the regularized energies converge as ε → 0 to gradient flows of the unregularized energy, recovering a generalization of Lions and Mas-Gallic's results when m = 2; see Theorem 5.8 and Corollary 5.9. Finally, in Section 6, we prove the convergence of our numerical blob method, under sufficient regularity assumptions, when the initial particle spacing h scales with the regularization like h = o(ε); see Theorem 6.1.
We close with several numerical examples, in one and two dimensions, analyzing the rate of convergence to exact solutions with respect to the 2-Wasserstein metric, L 1 -norm, and L ∞ -norm and illustrating qualitative properties of the method, including asymptotic behavior of the Fokker-Planck equation and critical mass of the two-dimensional Keller-Segel equation; see Section 6.3. In particular, for the heat equation and porous medium equations (V = W = 0, m = 1, 2, 3), we observe that the 2-Wasserstein error depends linearly on the grid spacing h ∼ N −1/d for m = 1, 2, 3, while the L 1 -norm depends quadratically on the grid spacing for m = 1, 2 and superlinearly for m = 3. We apply our method to study long time behavior of the nonlinear Fokker-Planck equation (V = |·| 2 /2, W = 0, m = 2), showing that the blob method accurately captures convergence to the unique steady state. Finally, we conduct a detailed numerical study of equations of Keller-Segel type, including a one-dimensional variant (V = 0, W = 2χ log |·| , χ > 0, m = 1, 2) and the original two-dimensional equation (V = 0, W = ∆ −1 , m = 1). The one-dimensional equation has a critical mass 1, and the two-dimensional equation has critical mass 8π, at which point the concentration effects from the nonlocal interaction term balance with linear diffusion (m = 1) [12,41]. We show that the same notion of criticality is present in our numerical solutions and demonstrate convergence of the critical mass as the grid spacing h and regularization ε are refined.
There are several directions for future work. Our convergence theorem for m ≥ 2 requires additional regularity assumptions, which we are only able to remove in the case m = 2 when the initial data has bounded entropy. In the case of m > 2 or more general initial data, it remains an open question how to control certain nonlocal norms of the regularized energies, which play an important role in our convergence result; see Theorem 5.8. Formally, we expect these to behave as approximations of the BV -norm of ρ m , which should remain bounded by the gradient flow structure; see equations (24) and (25). When 1 ≤ m < 2, it is not clear how to use these nonlocal norms to get the desired convergence result or whether an entirely different approach is needed. Perhaps related to these questions is the fact that our estimate on the semiconvexity of the regularized energies (6) deteriorates as ε → 0, while we expect that the semiconvexity should not deteriorate along smooth geodesics; see Proposition 3.11. Finally, while our results show convergence of the blob method for diffusive Wasserstein gradient flows, they do not quantify the rate of convergence in terms of h and ε. In particular, a theoretical result on the optimal scaling relation between h and ε remains open, though we observe good numerical performance for ε = h 1−p , 0 < p 1. In a less technical direction, we foresee a use of the presented ideas in conjunction with splitting schemes for certain nonlinear kinetic equations [1,20], as well as in the fluids [48], since our numerical results demonstrate comparable rates of convergence to the particle strength exchange method, which has already gained attention in these contexts [40].

Preliminaries
2.1. Basic notation. For any r > 0 and x ∈ R d we denote the open ball of center x and radius r by B r (x). Given a set S ⊂ R d , we write 1 S : R d → {0, 1} for the indicator function of S, i.e., 1 S (x) = 1 for x ∈ S and 1 S (x) = 0 otherwise. We say a function A : R d → R has at most quadratic growth if there exist c 0 , c 1 > 0 so that |A(x)| ≤ c 0 + c 1 |x| 2 for all x ∈ R d .
Let P(R d ) denote the set of Borel probability measures on R d , and for, any p ∈ N, P p (R d ) denotes elements of P(R d ) with finite pth moment, M p (R d ) := R d |x| p dµ(x) < +∞. We write L d for the d-dimensional Lebesgue measure, and for given µ ∈ P(R d ), we write µ L d if µ is absolutely continuous with respect to the Lebesgue measure. Often we use the same symbol for both a probability measure and its Lebesgue density, whenever the latter exists. We let L p (µ; R d ) denote the Lebesgue space of functions with pth power integrable against µ.
Given σ a finite, signed Borel measure on R d , we denote its variation by |σ|. For a Borel set E ⊂ R d we write σ(E) for the σ-measure of set E. For a Borel map T : R d → R d and µ ∈ P(R d ), we write T # µ for the push-forward of µ through T . We let id : For a sequence (µ n ) n ⊂ P(R d ) and some µ ∈ P(R d ), we write µ n * µ if (µ n ) n converges to µ in the weak- * topology of probability measures, i.e., in the duality with bounded continuous functions.

Convolution of measures.
A key aspect of our approach is the regularization of the energy (2) via convolution with a mollifier. In this section, we collect some elementary results on the convolution of probability measures, including a mollifier exchange lemma, Lemma 2.2. For any µ ∈ P(R d ) and measurable function φ, the convolution of φ with µ is given by whenever the integral converges. We consider mollifiers satisfying the following assumption.
This assumption is satisfied by both Gaussians and smooth functions with compact support. Assumption 2.1 also ensures that ϕ has finite first moment. For any ε > 0, we write Throughout, we use the fact that the definition of convolution allows us to move mollifiers from the measure to the integrand. In particular, for any φ bounded below and ψ ∈ L 1 (R d ) even, Likewise, we often use the following lemma which provides sufficient conditions for moving functions in and out convolutions with mollifiers within integrals. (See also [59] for a similar result.) The proof is contained in Appendix A.
LEMMA 2.2 (mollifier exchange lemma). Let f : R d → R be Lipschitz continuous with constant L f > 0, and let σ and ν be finite, signed Borel measures on R d . There is p = p(q, d) > 0 so that We conclude this section with a proposition stating that if a sequence of measures converges in the weak- * topology of P(R d ), then the mollified sequence converges to the same limit. We refer the reader to Appendix A for the proof. PROPOSITION 2.3. Let µ ε be a sequence in P(R d ) such that µ ε * µ as ε → 0 for some µ ∈ P(R d ). Then ϕ ε * µ ε * µ.

2.3.
Wasserstein metric. For µ, ν ∈ P(R d ), we denote the set of transport plans from µ to ν by Γ(µ, ν) := {γ ∈ P(R d × R d ) | π 1 # γ = µ, π 2 # γ = ν}, where π 1 , π 2 : R d × R d → R d are the projections of R d × R d onto the first and second copy of R d , respectively. The Wasserstein distance W 2 (µ, ν) between two probability measures µ, ν ∈ P 2 (R d ) is given by and a transport plan γ o is optimal if it attains the minimum in (9). We denote the set of optimal transport plans by Γ o (µ, ν). If either µ is absolutely continuous with respect to Lebesgue measure, then there is a unique optimal transport plan γ o , and T o is unique up to sets of µ-measure zero is known as the optimal transport map from µ to ν. Convergence with respect to the Wasserstein metric is stronger than weak- * convergence. In particular, if (µ n ) n ⊂ P 2 (R d ) and µ ∈ P 2 (R d ), then W 2 (µ n , µ) → 0 as n → ∞ ⇐⇒ µ n * µ and M 2 (µ n ) → M 2 (µ) as n → ∞ .
In order to define Wasserstein gradient flows, we will require the following notion of regularity in time with respect to the Wasserstein metric.
Along such curves, we have a notion of metric derivative. DEFINITION 2.5 (metric derivative). Given µ ∈ AC 2 loc ((0, ∞); P 2 (R d )), its metric derivative is An important class of curves in the Wasserstein metric are the (constant speed) geodesics. Given µ 0 , µ 1 ∈ P 2 (R d ), geodesics connecting µ 0 to µ 1 are of the form If γ o is induced by a map T o , then More generally, given µ 1 , µ 2 , µ 3 ∈ P 2 (R d ), a generalized geodesic connecting µ 2 to µ 3 with base µ 1 is given by such that π 1,2 # γ ∈ Γ o (µ 1 , µ 2 ) and π 1, the projection of onto the first and ith copies of R d . When the base µ 1 coincides with one of the endpoints µ 2 or µ 3 , generalized geodesics are geodesics.
A key property for the uniqueness and stability of Wasserstein gradient flows is that the energies are convex, or more generally semiconvex, along generalized geodesics. DEFINITION 2.6 (semiconvexity along generalized geodesics). We say a functional G : P 2 (R d ) → (−∞, ∞] is semiconvex along generalized geodesics if there is λ ∈ R such that for all µ 1 , µ 2 , µ 3 ∈ P 2 (R d ) there exists a generalized geodesic connecting µ 2 to µ 3 with base µ 1 such that For any subset X ⊂ P(R d ) and functional G : X → (−∞, ∞], we denote the domain of G by D(G) = {µ ∈ X | G(µ) < +∞}, and we say that G is proper if D(G) = ∅. As soon as a functional is proper and lower semicontinuous with respect to the weak-* topology, we may define its subdifferential ; see [3, Definition 10.3.1 and Equation 10.3.12]. Following the approach in [24], the notion of subdifferential we use in this paper is, in fact, the following reduced one. DEFINITION 2.7 (subdifferential). Given G : P 2 (R d ) → (−∞, ∞] proper and lower semicontinuous, µ ∈ D(G), and ξ : R d → R d with ξ ∈ L 2 (µ; R d ), then ξ belongs to the subdifferential of G at µ, written ξ ∈ ∂G(µ), if as ν The Wasserstein metric is formally Riemannian, and we may define the tangent space as follows. 8 DEFINITION 2.8. Let µ ∈ P 2 (R d ). The tangent space at µ is where the closure is taken in L 2 (µ; R d ).
We now turn to the definition of a gradient flow in the Wasserstein metric (c.f.
for almost every t > 0 such that µ is a weak solution of the continuity equation i.e., µ is a solution to the continuity equation in duality with C ∞ c (R d ). We close this section with the following definition of the Wasserstein local slope.
where the subscript + denotes the positive part.
REMARK 2.11. When the functional G in Definition 2.9 is in addition semiconvex along geodesics the local slope |∂G| is a strong upper gradient for G. In this case a gradient flow of G is characterized as being a 2-curve of maximal slope with respect to |∂G|; see [3, Theorem 11.1.3].

Regularized internal energies
The foundation of our blob method is the regularization of the internal energy F via convolution with a mollifier. This allows us to preserve the gradient flow structure and approximate our original partial differential equation (1) by a sequence of equations for which particles do remain particles. In this section, we consider several fundamental properties of the regularized internal energies F ε , including convexity, lower semicontinuity, and differentiability. In what follows, we will suppose that our internal energies satisfy the following assumption. ASSUMPTION 3.1 (internal energies). Suppose F ∈ C 2 (0, +∞) satisfies lim s→+∞ F (s) = +∞ and either F is bounded below or lim inf s→0 F (s)/s β > −∞ for some β > −2/(d + 2). Suppose further that U (s) = sF (s) is convex, bounded below, and lim s→0 U (s) = 0.
Thanks to this assumption we can define the internal energy corresponding to F by otherwise.
If F is bounded below, this is well-defined on all of P(R d ). If lim inf s→0 F (s)/s β > −∞ for some β > −2/(d + 2), this is well-defined on P 2 (R d ); see [3, Example 9.3.6]. which leads to F (s) ≥ 0 for all s ∈ (0, ∞). 9 Our assumption does not ensure that F is convex along Wasserstein geodesics, unless F is convex. which, by Remark 3.2, holds when for example F is convex.
We regularize the internal energies by convolution with a mollifier.
DEFINITION 3.4 (regularized internal energies). Given F : (0, ∞) → R satisfying Assumption 3.1, we define, for all µ ∈ P(R d ), the regularized internal energies by Note that, for all µ ∈ P(R d ) and An important class of internal energies satisfying Assumption 3.1 are given by the (negative) entropy and Rényi entropies.
Note that, as per our observation just below the definition of F, the entropy F 1 is well-defined on P 2 (R d ) and the Rényi entropies (F m , m > 1) are well-defined on all of P(R d ). Also note that the regularized entropies (F m ε , m ≥ 1, ε > 0) are well-defined on all of P(R d ). In order to approximate solutions of equation (1), we will consider combinations of the above regularized internal energies with potential and interaction energies.
When F = F m for some m ≥ 1, then we denote E by E m and E ε by E m ε . The regularized internal energy in Definition 3.4 incorporate a blend of interaction and internal phenomena, through the convolution with the mollifier, or potential, ϕ ε and the composition with the function F . To our knowledge, this is a novel form of functional on the space of probability measures. We now describe some of its basic properties: energy bounds and lower semicontinuity, when F is the logarithm or a power, and differentiability, convexity and subdifferential characterization when F is convex. For the existence and uniqueness of gradient flows associated to this regularized energy, see Section 5. REMARK 3.7. Although the regularized energy in Definition 3.4 is of a novel form, it was noticed in [69, Proposition 6.9] that a previous particle method for diffusive gradient flows leads to a similar regularized internal energy after space discretization [25,30]. The essential difference between these two methods stands in the choice of the mollifier, which, instead of satisfying 2.1, is a very singular potential.
We begin with inequalities relating the regularized internal energies to the unregularized energies. See Appendix A for the proof, which is a consequence of Jensen's inequality and a Carleman-type estimate on the lower bound of the entropy [30, Lemma 4.1].
For all ε > 0, the regularized entropies are lower semicontinuous with respect to weak-* convergence (m > 1) and Wasserstein convergence (m = 1). For m > 2, we prove this using a theorem of Ambrosio, Gigli, and Savaré on the converge of maps with respect to varying measures; see Proposition B.2. For 1 < m ≤ 2, this is a consequence of Jensen's inequality. For m = 1, we apply both Jensen's inequality and a version of Fatou's lemma for varying measures; see Lemma B.3. In this case, we also require that the mollifier ϕ is a Gaussian, so that we can get the bound from below required by Fatou's lemma. We refer the reader to Appendix A for the proof.
When F is convex, the regularized internal energies are differentiable along generalized geodesics. The proof relies on the fact that F is differentiable and ϕ ε ∈ C 2 (R d ), with bounded Hessian; see Appendix A. PROPOSITION 3.10 (differentiability). Let F satisfy Assumption 3.1 and be convex. Given A key consequence of the preceding proposition is that the regularized energies are semiconvex along generalized geodesics, as we now show. PROPOSITION 3.11 (convexity). Suppose F satisfies Assumption 3.1 and is convex. Then F ε is λ F -convex along generalized geodesics, where Proof. Let (µ 2→3 α ) α∈[0,1] be a generalized geodesic connecting two probability measures µ 2 , µ 3 ∈ P 2 (R d ) with base µ 1 ∈ P 2 (R d ); see (10). We have, using the above-the-tangent inequality for convex functions, Therefore, by Proposition 3.10, , which gives the result.
We now use the previous results to characterize the subdifferential of the regularized internal energy. The structure of argument is classical (c.f. [3,24,54]), but due to the novel form of our regularized energies, we include the proof in Appendix A.
As a consequence of this characterization of the subdifferential, we obtain the analogous result for the full energy E ε , as in Definition 3.6. See Appendix A for the proof. COROLLARY 3.13. Suppose F satisfies Assumption 3.1 and is convex. Let ε > 0 and µ ∈ D(E ε ). Suppose V, W ∈ C 1 (R d ) are semiconvex, with at most quadratic growth, and W is even. Then .

Γ-convergence of regularized internal energies
We now turn to the convergence of the regularized energies and, when in the presence of confining drift or interaction terms, the corresponding convergence of their minimizers. In this section, and for the remainder of the work, we consider regularized entropies and Rényi entropies of the form F m ε for m ≥ 1. We begin by showing that F m ε Γ-converges to F as ε → 0. THEOREM 4.1 (F ε Γ-converges to F). If m = 1, consider (µ ε ) ε ⊂ P 2 (R d ) and µ ∈ P 2 (R d ), and if m > 1, consider (µ ε ) ε ⊂ P(R d ) and µ ∈ P(R d ).
We begin by showing the result for 1 ≤ m ≤ 2, in which case the function F is concave. We first show part (i). By Proposition 3.8, for all ε > 0, which gives the result. We now turn to part (ii). Again, by Proposition 3.8, for all ε > 0, . We now consider the case when m > 2. Part (ii) follows quickly: by Proposition 3.8, Young's convolution inequality, and the fact that . Taking the supremum limit as ε → 0 then gives the result. Let us prove part (i). Without loss of generality, we may suppose that lim inf ε→0 F m ε (µ ε ) is finite. Furthermore, there exists a positive sequence (ε n ) n such that ε n → 0 and lim n→+∞ F m εn (µ εn ) = lim inf ε→0 F m ε (µ ε ). In particular, there exists C > 0 for which F m εn (µ εn ) < C for all n ∈ N. By Jensen's inequality for the convex function x → x m−1 and the fact that ζ ε * ζ ε = ϕ ε for all ε > 0, Thus, since F m εn (µ εn ) < C for all n ∈ N, we have ζ εn * µ εn L 2 (R d ) < C := (C(m − 1)) 1/2(m−1) . We now use this bound on the L 2 -norm of ζ εn * µ εn to deduce a stronger notion of convergence of ζ εn * µ εn to µ. First, since (µ εn ) n converges weakly- * to µ as n → ∞, Proposition 2.3 ensures that (ζ εn * µ εn − µ εn ) n converges weakly- * to 0. Since the L 2 -norm is lower semicontinuous with respect to weak- * convergence [64, Lemma 3.4], we have so that µ ∈ L 2 (R d ). Furthermore, up to another subsequence, we may assume that (ζ εn * µ εn ) n converges weakly in so (ζ εn * µ εn ) n converges to µ weakly in L 2 . By the Banach-Saks theorem (c.f. [71, Section 38]), up to taking a further subsequence of (ζ εn * µ εn ) n , the Cesàro mean (v k ) k defined by We now use this stronger notion convergence to conclude our proof of part (i). Since m > 2 and Furthermore, recalling the definition of the regularized energy and applying [3, Therefore, to finish the proof, it suffices to show that w(x) ≥ µ(x) for µ-almost every x ∈ R d . By Lemma 2.2 and the fact that ζ εn * ζ εn = ϕ εn for all n ∈ N, there exists p > 0 and C ζ > 0 so that Combining this with equation (18), we obtain Finally, using equation (17) and the definition of (v k ) k as a sequence of convex combinations of the family Since the limit in (19) exists, it coincides with its Cesàro mean on the right-hand side of the above This gives w(x) ≥ µ(x) for µ-almost every x ∈ R d , which completes the proof. Now, we add a confining drift or interaction potential to our internal energies, so that energy minimizers exist and we may apply the previous Γ-convergence result to conclude that minimizers converge to minimizers. For the remainder of the section we consider energies of the form E m ε given in Definition 3.6, with the following additional assumptions on V and W to ensure that the energy is confining. 14 ASSUMPTION 4.2 (confining assumptions). The potentials V and W are bounded below and one of the following additional assumptions holds: V has compact sublevel sets; (CV) Under these assumptions, the regularized energies E m ε are lower semicontinuous with respect to weak- * convergence (m > 1) and Wasserstein convergence (m = 1), where for the latter we assume ϕ is a Gaussian (c.f. Proposition 3.9, and [3, Lemma 5.
We now prove existence of minimizers of E m ε , for all ε > 0.
ε is bounded below. By Remark 4.3, if (CV) holds, then any minimizing sequence of E m ε has a subsequence that converges in the weak- * topology. Likewise, if (CW) holds, then any minimizing sequence of E m ε has a subsequence that, up to translation, converges in the weak- * topology. By lower semicontinuity of E m ε , the limits of minimizing sequences are minimizers of E m ε . Now, suppose m = 1. By Proposition 3.8, for all δ > 0 and µ ∈ P 2 (R d ), , Consequently, by the assumption in (CV ) and the fact that W is bounded below by, say,C ∈ R, we can choose δ = C 0 /2 and obtaiñ Hence any minimizing sequence (µ n ) n ⊂ P 2 (R d ) has bounded second moment. Thus, (µ n ) n has a subsequence that converges in the weak- * topology to a limit with bounded second moment. By the lower semicontinuity of E m ε the limit must be a minimizer of E m ε . Finally, we conclude that minimizers of the regularized energy converge to minimizers of the unregularized energy.
holds, then we have µ ε * µ, up to a subsequence and translation, where again µ minimizes E m . Now suppose m = 1. If Assumption (CV ) holds and ϕ is a Gaussian, then for any sequence Proof. The proof is classical. We include it for completeness. We only prove the result under Assumptions (CV)/(CV ) since the argument for (CW) is analogous. For any ε > 0, since µ ε is a minimizer of E m ε we have that E m ε (µ ε ) ≤ E m ε (ν) for all ν ∈ P(R d ) 15 if m > 1, and for all ν ∈ P 2 (R d ) if m = 1. Taking the infimum limit of the left-hand side and the supremum limit of the right-hand side, Theorem 4.1(ii) ensures that Since E m is proper there exists ν ∈ P(R d ) if m > 1 and ν ∈ P 2 (R d ) if m = 1 so that the right-hand side is finite. Thus, up to a subsequence, we may assume that {E m ε (µ ε )} ε is uniformly bounded. When m > 1, F ε (µ) ≥ 0 for all ε ≥ 0, and this implies that { V dµ ε } ε is uniformly bounded, so {µ ε } ε is tight. When m = 1, the inequality in (20) ensures that {M 2 (µ ε )} ε is uniformly bounded, so again {µ ε } ε is tight. Thus, up to a subsequence, (µ ε ) ε converges weakly- * to a limit µ ∈ P(R d ) if m > 1 and µ ∈ P 2 (R d ) if m = 1. By Theorem 4.1(i) and the inequality in (21), we obtain REMARK 4.6 (convergence of minimizers). One the main difficulties for improving the topology in which the convergence of the minimizers happen is that we do not control L m -norms of the regularized minimizing sequences due to the special form of our regularized energy. This is the main reason we only get weak- * convergence in the previous result and the main obstacle to improve results for the Γ-convergence of gradient flows, as we shall see in the next section.

Γ-convergence of gradient flows
We now consider gradient flows of the regularized energies E m ε , as in Definition 3.6, for m ≥ 2 and prove that, under sufficient regularity assumptions, gradient flows of the regularized energies converge to gradient flows of the unregularized energy as ε → 0. For simplicity of notation, we often write E m ε and F m ε for ε ≥ 0 when we refer jointly to the regularized and unregularized energies. We begin by showing that the gradient flows of the regularized energies are well-posed, provided that V and W satisfy the following convexity and regularity assumptions.
ASSUMPTION 5.1 (convexity and regularity of V and W ). The potentials V, W ∈ C 1 (R d ) are semiconvex, with at most quadratic growth, and W is even. Furthermore, there exist C 0 , C 1 > 0 so REMARK 5.2 (ω-convexity). More generally, our results naturally extend to drift and interaction energies that are merely ω-convex; see [35]. However, given that the main interest of the present work is approximation of diffusion, we prefer the simplicity of Assumption (5.1), as it allows us to focus our attention on the regularized internal energy.
ε is as in Definition 3.6 and V and W satisfy Assumption 5.1. Then, for any µ 0 ∈ D(E m ε ), there exists a unique gradient flow of E m ε with initial datum µ 0 .
Proof. It suffices to verify that E m ε is proper, coercive, lower semicontinuous with respect to 2-Wasserstein convergence, and semiconvex along generalized geodesics; c.f. [3,Theorem 11.2.1]. (See also [3, Equation (2.1.2b)] for the definition of coercive.) If ε > 0, then F m ε is finite on all of P 2 (R d ), and if ε = 0, then F m is proper. Thus, our assumptions on V and W ensure that E m ε is proper. Clearly F m ε is bounded below. Hence, since the semiconvexity of V and W ensures that their negative parts have at most quadratic growth, E m ε is coercive. For ε > 0, Proposition 3.9 ensures that F m ε is lower semicontinuous with respect to weak- * convergence, hence also 2-Wasserstein convergence. For ε = 0, the unregularized internal energy F m is also lower semicontinuous with respect to weak- * and 2-Wasserstein convergence [64,Lemma 3.4]. Since V and W are lower semicontinuous and their negative parts have at most quadratic growth, the associated potential and interaction energies are lower semicontinuous with respect to 2-Wasserstein convergence [3, Lemma 5.1.7, Example 9.3.4]. Therefore, E m ε is lower semicontinuous for all ε ≥ 0.
In the case ε = 0, gradient flows of the energies E m are characterized as solutions of the partial differential equation (1); c.f. [3, Theorems 10.4.13 and 11.2.1], [24,Theorem 2.12]. Now, we show that gradient flows of the regularized energies E m ε can also be characterized as solutions of a partial differential equation. Moreover, Then, by Definition 2.9 and Corollary 3.13, µ ε is a weak solution to the continuity equation with velocity field (22). Conversely, suppose µ ε is a weak solution to the continuity equation with velocity field (22). By A consequence of the previous proposition is that, for the regularized energies E m ε , particles remain particles, i.e. a solution of the gradient flow with initial datum given by a finite sum of Dirac masses remains a sum of Dirac masses, and the evolution of the trajectories of the particles is given by a system of ordinary differential equations.
is well-posed for all T > 0. Furthermore, µ ε = i∈I δ X i (·) m i belongs to AC 2 ([0, T ]; P 2 (R d )) and is the gradient flow of E m ε with initial conditions µ ε (0) : Proof. To see that (23) is well-posed, first note that the function is Lipschitz. Likewise, Assumption 5.1 ensures y i → ∇V (y i ) and y i → j∈I ∇W (y i − y j ) are continuous and one-sided Lipschitz. Therefore, the ODE system (23) is well-posed forward in time.
on an interval [0, T ], for some fixed T . We abbreviate by v i = v i (X 1 , X 2 , . . . , X N ) the velocity field for X i in (23). For any test function ϕ ∈ C ∞ c (R d × (0, T )), the fundamental theorem of calculus ensures that, for all i ∈ I, Combining this with (23), we obtain Multiplying both sides by m i , summing over i, and taking for v as in (22). Therefore, µ ε is a weak solution of the continuity equation with velocity field v.
We now turn to the Γ-convergence of the gradient flows of the regularized energies, using the scheme introduced by Sandier-Serfaty [75] and then generalized by Serfaty [77], which provides three sufficient conditions for concluding convergence. We will use the following variant of Serfaty's result, which allows for slightly weaker assumptions on the gradient flows of the regularized energies, but follows from the same argument as Serfaty's original result. (See also Remark 2.11 on the correspondence between Wasserstein gradient flows and curves of maximal slope.) THEOREM 5.6 (c.f. [77,Theorem 2]). Let m ≥ 2. Suppose that, for all ε > 0, µ ε belongs to AC 2 ([0, T ]; P 2 (R d )) and is a gradient flow of E m ε with well-prepared initial data, i.e., Suppose further that there exists a curve µ in P 2 (R d ) such that, for almost every t ∈ [0, T ], µ ε (t) * µ(t) and (S1) lim inf Then µ ∈ AC 2 ([0, T ]; P 2 (R d )), and µ is a gradient flow of E m .
For simplicity of notation, in what follows we shall at times omit dependence on time when referring to curves in the space of probability measures.
In order to apply Serfaty's scheme in the present setting to obtain Γ-convergence of the gradient flows, a key assumption is that the following quantity is bounded uniformly in ε > 0 along the gradient flows µ ε of the regularized energies E m ε : where we use the abbreviation p ε := (ϕ ε * µ ε ) m−2 µ ε . This quantity differs from ∇δF m ε /δµ ε L 1 (µε;R d ) merely by the placement of the absolute value sign: Serfaty's scheme allows one to assume, without loss of generality, that |F m ε |(µ ε ) is bounded uniformly in ε > 0 for almost everywhere t ∈ [0, T ], and Hölder's inequality ensures that |F m ε |(µ ε ) = ∇δF m ε /δµ ε L 2 (µε;R d ) ≥ ∇δF m ε /δµ ε L 1 (µε;R d ) ; see Proposition 3.12. Consequently, we miss the bound we require on µ ε BV m ε merely by placement of the absolute value sign in inequality (24). Still, µ ε BV m ε has a useful heuristic interpretation. Through the proof of Theorem 5.8, we obtain see the inequality (33) and Proposition B.2. Consequently, one may think of µ ε BV m ε as a nonlocal approximation of the L 1 -norm of the gradient of µ m .
We begin with a technical lemma we shall use to prove the convergence of the gradient flows.
Proof. We argue similarly as in Lemma 2.2. Let f : By Assumption 2.1, C ζ is so that ζ(x) ≤ C ζ |x| −q for q > d + 1 for all x ∈ R d . Chooser so that Now, we break the integral with respect to dµ ε (y) above into integrals over the domain B εr (x) and R d \ B εr (x), bounding the above quantity by First, we consider I 1 . Since, in the integral, |x − y| < εr, we obtain We apply the inequality in (52) to obtain ζ ε (x − y)|x − y| ≤ C ζ εr with r :=r(1 − q) + q − d in the integral-the inequality in (26) ensuresr > 1. Consequently, where, in the last inequality, we use that ∇ζ ε L 1 (R d ) = ∇ζ L 1 (R d ) /ε and, by Jensen's inequality for the concave function s (m−2)/(m−1) , Since 0 ≤ (m − 2)/(m − 1) < 1, Jensen's inequality gives This gives the result by taking r := min(r,r − 1).
With this technical lemma in hand, we now turn to the Γ-convergence of the gradient flows.
Furthermore, note that if µ m ∈ W 1,1 (R d ) and ∇µ m + ∇V µ + (∇W * µ)µ = ξµ for some ξ ∈ L 2 (µ; R d ), Observe that Proposition B.2 is stated for probability measures-we can easily rescale dµ ε ⊗ dL d to be a probability measure by diving the above equations by T > 0.
We conclude this section by showing that, in the case when m = 2 and for V, W ∈ C 2 (R d ) with bounded Hessians, whenever the initial data of the gradient flows have bounded second moments and internal energies, we automatically obtain Assumptions (A1)-(A3). Consequently, in this special case, we are able to conclude the convergence of the gradient flows without these additional assumptions.

REMARK 5.10 (Previous work, m = 2). The above theorem generalizes a result by Lions and Mas-Gallic [59] on a numerical scheme for the porous medium equation ∂ t µ = ∆µ 2 on a bounded domain with periodic boundary conditions to equations of the form (1) on Euclidean space.
Proof of Corollary 5.9. First, we show that sup ε>0 ζ ε * (µ ε (0)) L 2 (R d ) < ∞. The fact that D 2 V and D 2 W are bounded ensures |V | and |W | grow at most quadratically. Combining this with equations (40)- (41), which ensure {E m ε (µ(0))} ε and {M 2 (µ ε (0))} ε are bounded uniformly in ε > 0, we obtain Furthermore, since the energy F 2 ε decreases along solutions to the gradient flow, we have sup Next, we show that our assumption that the initial data has bounded entropy (41) ensures , for some C > 0 depending on d, V , W and sup ε>0 log µ ε (0) dµ ε (0). Formally differentiating the entropy F 1 (µ) = log(µ) dµ along the gradient flows µ ε , we expect that, for all t ∈ [0, T ], Hence, for any t ∈ [0, T ], This computation can be made rigorous by first proving the analogous inequality along discrete time gradient flows using the flow interchange method of Matthes, McCann, and Savaré [63, Theorem 3.2] and then sending the timestep to zero to recover the above inequality in continuous time. Thus, there exists K 0 > 0 depending on V, W and sup ε>0 F 1 (µ ε (0)) so that, for all t ∈ [0, T ], Finally, by a Carleman-type estimate [30,Lemma 4 Now, we use this estimate to show that {M 2 (µ ε (t))} ε is uniformly bounded in ε for all t ∈ [0, T ]. Let κ be a smooth cutoff function with 0 ≤ κ ≤ 1, for all |x| < 1/2 and κ(x) = 0 for all |x| > 2. Given R > 0, define κ R (x) = κ(x/R), so that ∇κ R L ∞ (R d ) ≤ 1/R and D 2 κ R L ∞ (R d ) ≤ 4/R 2 . Then there exists C κ > 0 so that for all R > 1, |∇(κ R (x)x 2 )| ≤ C κ |x| and |D 2 (κ R (x)x 2 )| ≤ C κ for all x ∈ R d . By Proposition 5.4, µ ε is a weak solution of the continuity equation. Therefore choosing κ R (x)|x| 2 as our test function, we obtain, for all t ∈ [0, T ], Since D 2 V and D 2 W are bounded, |∇V | and |∇W | grow at most linearly. Consequently, there exists C > 0, depending on V , W , and C κ so that Likewise, by Lemma 5.7, there exists r > 0 so that, for all t ∈ [0, T ], for C depending on C κ , sup ε>0 F 2 ε (µ ε (0)), and ∇ζ L 1 (R d ) . In the second inequality, we use that Therefore, there exists C > 0 so that, for all t ∈ [0, T ], As the right-hand side is independent of R > 1, by sending R → +∞ by the dominated convergence theorem we obtain that for ε r < 1/(2C ), Therefore, by Gronwall's inequality, there existsC depending on C , C and T (and independent of ε) so that 26 We may combine this with the inequality in (43) to obtain, for all t ∈ [0, T ], We now use these results to verify the assumptions of Theorem 5.8 hold, so that we may apply this result to conclude convergence of the gradient flows. Assumption (A1) is a consequence of the inequality in (45). Assumption (A2) is a consequence of the inequalities in (42), (44) and (46).

Numerical results
6.1. Numerical method and convergence. We now apply the theory of regularized gradient flows developed in the previous sections to develop a blob method for diffusion, allowing us to numerically simulate solutions to partial differential equations of Wasserstein gradient flow type (1). We begin by describing the details of our numerical scheme and applying Theorem 5.8 to prove its convergence, under suitable regularity assumptions.
where Q i is the cube centered at ih of side length h. Next, for ε > 0, define the evolution of these measures by where {X i (t)} i∈Q h R are solutions to the ODE system (23) on a time interval [0, T ] with initial data X i (0) = ih. If h = o(ε) as ε → 0 and Assumptions (A1)-(A3) from Theorem 5.8 hold, then (µ ε (t)) ε converges in the weak- * topology to µ(t) as ε → 0 for almost every t ∈ [0, T ], where µ(t) is the unique solution of (1) with initial datum µ(0).
For simplicity of notation, we suppress the dependence on time and show ζ ε * µ ε → µ in L m as ε → 0. By the assumptions that µ ∈ D(E m ) with compact support and V and W are continuous, we have µ ∈ L m (R d ). Consequently ζ ε * µ → µ in L m as ε → 0, and it is enough to show that ζ ε * µ ε − ζ ε * µ → 0 in L m . Using that T is a transport map from µ ε to µ, Combining the decay of ∇ζ from Assumption 2.1 with the fact that ∇ζ is continuous, there exists C > 0 so that |∇ζ(x)| ≤ C(1 B (x)+|x| −q 1 R d \B (x)), where B = B 1 (0) is the unit ball centered at the origin. Note that if |x−y| ≥ 2h, then for all α ∈ [0, 1], |x−(1−α)T (y)−αy| ≥ |x−y|−h ≥ |x−y|/2 and |x − (1 − α)T (y) − αy| ≤ 3|x − y|/2. Thus, by the assumptions on our mollifier, we have Thus, taking the L m -norm with respect to x, doing a change of variables, and applying Minkowski's inequality, we obtain where c > 0 depends on C, ∇ζ ∞ , and the space dimension. Therefore, provided that h = o(ε) as ε → 0, we obtain that ζ ε * µ ε − ζ ε * µ → 0 in L m . REMARK 6.2 (compact support of initial data). In Theorem 6.1, we assume that the initial datum of the exact solution µ(0) ∈ D(E m ) is compactly supported. More generally, under the same assumptions on V , W , and m, given any ν 0 ∈ D(E m ) ∩ P 2 (R d ) without compact support, there existsν 0 ∈ D(E m ) with compact support such that ν 0 andν 0 are arbitrarily close in the Wasserstein distance. Furthermore, by the contraction inequality for gradient flows of E m , the solution ν with initial data ν 0 and the solutionν with initial dataν 0 satisfy where C > 0 depends on T and the semiconvexity of V and W [3, Theorem 11.2.1]. In this way, any solution of (1) with initial datum in D(E m ) ∩ P 2 (R d ) can be approximated by a solution with compactly supported initial datum, so that our assumption of compact support in Theorem 6.1 is not restrictive. REMARK 6.3 (Assumptions (A1)-(A3)). In Theorem 6.1, we prove that, as long as Assumptions (A1)-(A3) from Theorem 5.8 hold along the particle solutions {µ ε } ε , then any limit of these particle solutions must be the corresponding gradient flow of the unregularized energy. Verifying these conditions analytically can be challenging; see Theorem 5.9. However, numerical results can provide confidence that these conditions hold along a given particle approximation. A sufficient condition for Assumption (A1) is that the (m − 1)th moment of the particle solution i∈Q h R |X i (t)| m−1 m i is bounded uniformly in t, ε, and h. In particular, this is satisfied if the particles remain compactly supported in a ball. A sufficient condition for Assumption (A2) is that with p ε = (ϕ ε * µ ε ) m−2 µ ε , remains bounded uniformly in t, ε, and h. In fact, for purely diffusive problems, we observe that this quantity not only bounded uniformly in ε and h, but decreases in time along our numerical solutions; see Figure 3 below. For the nonlinear Fokker-Planck equation, we observe that this quantity is bounded uniformly in ε and h and converges to the corresponding norm of the steady state as t → ∞; see Figure 6 below. A sufficient condition for Assumption (A3) is that the blob solution converges to a limit in L 1 and L ∞ , uniformly on bounded time intervals. Again, we observe this numerically, in both one and two dimensions, and both for purely diffusive equations and the nonlinear Fokker-Planck equation; see Figures 4-6 below. In this way, Assumptions (A1)-(A3) may be verified numerically in order to give confidence that the limit of any blob method solution is, in fact, the correct exact solution.
6.2. Numerical implementation. We now describe the details of our numerical implementation. In all of the numerical examples which follow, our mollifiers ζ ε and ϕ ε are given by Gaussians, In addition to Gaussian mollifiers, we also performed numerical experiments with a range of compactly supported and oscillatory mollifiers and observed similar results. In practice, Gaussian mollifiers provided the best balance between speed of computation and speed of convergence. We construct our numerical particle solutions µ ε (t) as described in Theorem 6.1. As a mild simplification, we consider the mass of each particle to be given by m i = µ(0, ih)h d , where µ(0, ih) is the value of the initial datum µ(0) at the grid point ih. For the numerical examples we consider, in which µ(0) is a continuous function, the rate of convergence is indistinguishable from defining m i as in (47).
The system of ordinary differential equations that prescribes the evolution of the particle locations (c.f. (23) and (48)) can be solved numerically in a variety of ways, and we observe nearly identical results independent of our choice of ODE solver. In analogy with previous work on blob methods in the fluids case [6], we find that the numerical error due to the choice of time discretization is of lower order than the error due to the regularization and spatial discretization. We implement the blob method in Python, using the Numpy, SciPy, and Matplotlib libraries [50,53,80]. In particular, we compute the evolution of the particle trajectories via the SciPy implementation of the Fortran VODE solver [15], which uses either a backward differentiation formula (BDF) method or an implicit Adams method, depending on the stiffness of the problem.
Our convergence result, Theorem 6.1, requires that h = o(ε) as ε → 0. Numerically, we observe the fastest rate of convergence with ε = h 1−p , for 0 < p 1, as h → 0. Since computational speed decreases as p approaches 0, we take ε = h 0.99 in the following simulations. In these examples, we discretize the initial data on a line (d = 1) or square of sidelength 5.0 (d = 2), centered at 0.
Finally, to visualize our particle solution (48) and compare it to the exact solutions in L p -norms, we construct a blob solution obtained by convolving the particle solution with a mollifier, By Proposition 2.3, if µ ε * µ as ε → 0, where µ is the exact solution, then we also haveμ ε * µ.
Consequently our convergence result, Theorem 6.1, also applies to this blob solution.
We measure the accuracy of our numerical method with respect to the L 1 -, L ∞ -, and Wasserstein metrics. To compute the L 1 -and L ∞ -errors, we take the difference between the exact solution and the blob solution (50) and evaluate discrete L 1 -and L ∞ -norms using the following formulas: We compute the Wasserstein distance between our particle solution µ ε in (48) and the exact solution µ in one dimension using the formula where F −1 µε and F −1 µ are the generalized inverses of the cumulative distribution functions of µ and µ ε , respectively; c.f. [3, Theorem 6.0.2]. We evaluate the integral in (51) numerically using the SciPy implementation of the Fortran library QUADPACK [70]. In two dimensions, we compute the Wasserstein error by discretizing the exact and blob solutions as piecewise constant functions on a fine grid and then using the Python Optimal Transport library to compute the discrete Wasserstein distance between them. In particular, we use the Earth Mover's Distance function in this library, which is based on the network simplex algorithm introduced by Bonneel, van de Panne, Paris, and Heidrich [14].
6.3. Simulations. Using the method described in the previous section, we now give several examples of numerical simulations. We consider initial data given by linear combinations of Gaussian and Barenblatt profiles, which we denote as follows: and κ = β 2 m − 1 m , and K = K(m, d) chosen so that ψ m (τ, x)dx = 1.
In Figure 1, we compare exact and numerical solutions to the heat and porous medium equations (V = W = 0, m = 1, 2, 3), with initial data given by a Gaussian (m = 1) or Barenblatt (m = 2, 3) function with scaling τ = 0.0625. The top row shows the evolution of the density on a large spatial scale, at which the exact and numerical solutions are visually indistinguishable for m = 1 and m = 2. However, for m = 3 the fat tails of the numerical simulation peel away from the exact solution at small times. The second row depicts the numerical simulations for m = 3 on a smaller spatial scale, illustrating how the tails of the numerical simulation converge to the exact solution as the spacing of the computational grid is refined.
In Figure 2, we compute solutions of the one-dimensional heat and porous medium equations (V = W = 0, m = 1, 2, 3), illustrating the role of the diffusion exponent m. The initial data is given by a linear combination of Gaussians, ρ 0 (·) = 0.3ψ 1 (· + 1, 0.0225) + 0.7ψ 1 (· − 1, 0.0225), and the grid spacing is h = 0.01. For m = 1, the infinite speed of propagation of support of solutions to the heat equation is reflected both at the level of the density, for which the gap between the two bumps fills quickly, and also in the particle trajectories, which quickly spread to fill in areas of low mass. In contrast, for m = 2 and m = 3, we observe finite speed of propagation of support, as well as the emergence of Barenblatt profiles as time advances. In Figure 3, we compute the evolution of the nonlocal Sobolev norm (49) along the numerical solutions from Figures 1 and 2. In both cases, we observe that the quantity converges as h → 0 and decreases in time. This gives further credence to the heuristic that the nonlocal Sobolev norm is an approximation of the L 1 -norm of the gradient of the mth power of the exact solution, which does decrease in time along the exact solution; see (24) and (25). In particular, this provides numerical evidence that Assumption (A2) from our main convergence theorem, Theorem 5.8, is satisfied.
In Figure 4, we analyze the rate of convergence of our numerical scheme in one dimension. We compute the error between numerical and exact solutions of the heat and porous medium equations (m = 1, 2, 3) in Figure 1 at time t = 0.05, with respect to the 2-Wasserstein distance, L 1 -norm, and L ∞ -norm and examine the scaling of the error with the grid spacing h. (Recall that ε = h 0.99 throughout.) Plotting the errors on a logarithmic scale, we observe that the Wasserstein error depends linearly on the grid spacing for all values of m. The L 1 -norm scales quadratically for m = 1 and 2 and superlinearly for m = 3. Finally, the L ∞ -error scales quadratically for m = 1, superlinearly for m = 2, and sublinearly for m = 3. This deterioration of the rate of L ∞ -convergence for m = 3 is due to the sharp transition at the boundary of the exact solution; see the second row of Figure 1. In Figure 5, we perform the same analysis on the rate of convergence of our method in two dimensions and observe similar rates of convergence as in the one-dimensional case.    In Figure 6, we simulate solutions to the nonlinear Fokker-Planck equation (V (·) = |·| 2 /2, W = 0, m = 2) and consider the rate of convergence to the steady state of the equation, ψ 2 (0.25, x). In the top row, we compute the error between the numerical solution at time t = 1.2 and the steady state with respect to the Wasserstein, L 1 -, and L ∞ -norms for various choices of grid spacing h. We consider solutions with Barenblatt initial data (m = 2, τ = 0.15). We plot the error's dependence on h with a logarithmic scale and compute the slope of the line of best fit to determine the scaling relationship between the error and h. We observe similar rates of convergence as in the case of the heat and porous medium equations; see Figure 5. In the middle rows, we give snapshots of the evolution of the blob method solution, as it converges to the steady state. We consider Barenblatt initial data (m = 2, τ = 0.15) and double bump initial data given by a linear combination of Barenblatts, ρ 0 (x) = 0.7ψ 2 (x−(1.25, 0), 0.1)+0.3ψ 2 (x+(1.25, 0), 0.1). The grid spacing is h = 0.02. In the bottom row, we compute the evolution of the nonlocal Sobolev norm (49) along the numerical solutions from the middle rows. In both cases, we observe that this quantity converges for h small. For Barenblatt initial data, it decreases in time along the numerical solution and agrees well with the value of ∇µ 2 L 1 (R d ) along the exact solution µ. For the double bump initial data, it remains bounded in time, converging asymptotically to ∇(ψ 2 (0.25, ·)) 2 L 1 , where ψ 2 (0.25, ·) is the steady state. Again, this supports the interpretation of the nonlocal Sobolev norm as an approximation of the L 1 -norm of the gradient of the mth power of the exact solution and provides numerical evidence for Assumption (A2) from Theorem 5.8.

34
Convergence Analysis: Two-Dimensional Diffusion Figure 5. Rate of convergence of blob method for two-dimensional heat and porous medium equations.
In the remaining numerical examples, we apply our method to simulate solutions of Keller-Segel type equations, with the interaction potential W given by 2χ log |·| for χ > 0. In one dimension, the derivative of this potential is not integrable, and we remove its singularity it setting it equal to 2χ/ε for all x ∈ R d such that |x| < ε. In two dimensions, the gradient of this potential is integrable, and we regularize it by convolving it with a mollifier ϕ ε , as done in previous work by the second author on a blob method for the aggregation equation [36].
In Figure 7, we consider the one-dimensional variant of the Keller-Segel equation (V = 0, W (·) = 2χ log |·|, m = 1) studied in [18]. Its interest is that it has a defined critical value χ for unit mass leading to the dichotomy of blow-up versus global existence. For χ = 1.5 and initial data of mass one, solutions blow up in finite time. We consider initial data given by a Gaussian ψ 1 (τ, ·), τ = 0.25, discretized on the interval [−4.5, 4.5] with grid spacing h = 0.009. We compare the evolution of the second moment of our blob method solutions with the second moment of the exact solution. We also compare our results with those obtained in previous work via a one-dimensional Discrete Gradient Flow (DGF) particle method [25,30]. By refining our spatial grid with respect to the DGF particle method, we observe modest improvements. (Alternative simulations, with similar spatial and time discretizations as used in the DGF method, yielded similar results as obtained by DGF.) The blow-up of solution is not only evident in the second moment, which converges to zero linearly in time, but also in the evolution of the particle trajectories. In particular, we observe particle trajectories merging on several occasions as time advances.   Figure 7. Left: Comparison of the evolution of the second moment along exact solutions (solid line) with blob method solutiosn (dashed line) and previous numerical results by the DGF particle method [25]. Right: Evolution of particle trajectories, with colors indicating the relative mass of each particle.  Figure 8. Left: Evolution of the second moment. Right: Evolution of particle trajectories, with colors indicating the relative mass of each particle.
In Figure 8, we consider a nonlinear variant of the Keller-Segel equation (V = 0, W (·) = 2χ log |·|, m = 2) in one dimension, with initial data and discretization as in Figure 7. We observe the convergence to a steady state both at the level of the second moment and the particle trajectories.
In Figures 9-12 we consider the classical Keller-Segel equation (V = 0, W (·) = 1/(2π) log |·|, m = 1) in two dimensions. In Figures 9, 10, and 11, the initial data is given by a Gaussian ψ 1 (τ, ·), τ = 0.16, scaled to have mass that is either supercritical (> 8π), critical (= 8π), or subcritical (< 8π) with respect to blowup behavior. In particular, for supercritical initial data, solutions blow up in finite time [13,41]. In Figure 9, we analyze the blow-up behavior. We compute the evolution of the second moment of solutions for fixed grid spacing h = 0.03 and varying mass 7π, 8π, and 9π, illustrating how initial data with larger mass aggregates more quickly at the origin. In Figure 10, we consider the evolution of the second moment for the solutions from Figure 9. For fixed grid spacing h = 0.03, we observe that the second moment depends linearly on time, and we compute its slope using the line of best fit. We then analyze how the slope of this line converges to the theoretically predicted slope as the grid spacing h → 0.
In Figure 11, we consider the evolution of the second moment for the supercritical mass solution from Figure 9 on a longer time interval. As in the one-dimensional case (see Figure 7), we are able to get approximately halfway to the time when the second moment becomes zero before the second moment of our numerical solution begins to peel away from the second moment of the exact solution. Indeed, one of the benefits of our blob method approach is that the numerical method naturally extends to two and more dimensions, and we observe similar numerical performance independent of the dimension. We also plot the evolution of particle trajectories, observing the tendency of trajectories in regions of larger mass to be driven largely by pairwise attraction, while trajectories in regions of lower mass feel more strongly the effects of diffusion.
Finally, in Figure 12, we consider the evolution of the density and second moment for double bump initial data, with initial mass 7π, 8π, and 9π. The slopes of the second moment agree well with the theoretically predicted slopes given in Figure 10.     Figure 6. In particular, we consider constant multiples M = 7π, 8π, and 9π and again observe that larger values of M correspond to faster aggregation at the origin. Bottom: We consider the evolution of the second moment along particle solutions, for each choice of M . We estimate the slope of the line using the line of best fit.

Appendix A. Proofs of preliminary results
We now turn to the proofs of some of the elementary lemmas and propositions from Sections 2 and 3. We begin with the proof of the mollifier exchange lemma. By Jensen's inequality for the convex function s → s log s, the relative entropy is nonnegative, which gives the result. Now, we show the left inequality in (11) for 1 < m ≤ 2. By the above-the-tangent property of the concave function F m and Hölder's inequality, we get , the first term goes to zero as ε → 0 and the second term remains bounded. This gives the result.
We now turn to the right inequality in (11) in the case 1 ≤ m ≤ 2. By the fact that ϕ ε = ζ ε * ζ ε and Jensen's inequality for the concave function F m , for all x ∈ R d we have Now, we show (12). Since F m is convex for m ≥ 2, this is simply a consequence of reversing the inequalities in the last two inequalities.
When m > 1, we simply use that F m ≥ 0.
We now give the proof that, for all ε > 0, the regularized energies are lower semicontinuous with respect to weak-* convergence (m > 1) and Wasserstein convergence (m = 1), where in the latter case, we require ϕ to be a Gaussian.
Proof of Proposition 3.9. We first show (i). Let (µ n ) n ⊂ P(R d ) and µ ∈ P(R d ) be such that µ n * µ; we must show lim inf n→∞ F m ε (µ n ) ≥ F m ε (µ). Without loss of generality, we replace (µ n ) n by a subsequence which attains the limit on the left-hand side, and we may assume that this limit is finite. Consequently, there exists C > 0 so that We now consider the case when 1 < m ≤ 2. For any a, b > 0, |a m−1 −b m−1 | ≤ |a−b| m−1 . Combining this with Jensen's inequality for the concave function s → s m−1 , Since ϕ ε ∈ C b (R d ), µ n * µ ensures that ζ ε * µ n → ζ ε * µ pointwise. The integrand of the first term is bounded above by 2 ζ ε 2 L ∞ (R d ) , so by the dominated convergence theorem, the first integral converges to zero. Since the integrand of the second term is continuous and bounded by ϕ ε m−1 L ∞ (R d ) , the fact that (µ n ) n converges weakly- * to µ ensures this second term converges to F ε (µ). This gives the result for 1 < m ≤ 2. We now deal with the case m > 2. Inequality (53) ensures that F m ε (µ n ) = ϕ ε * µ n L m−1 (µn;R d ) < C for all n ∈ N; so by Proposition B.2 there exists w ∈ L m−1 (µ; R d ) such that, up to another subsequence, for all It suffices to show that w ≥ ϕ ε * µ holds µ-almost everywhere. Then, the inequality in (54) gives Since µ n * µ and f, ζ ε ∈ C b (R d ), we have ζ ε * (f µ n ) → ζ ε * (f µ) and ζ ε * µ n → ζ ε * µ pointwise.
Now we turn to the proof that the regularized energies are differentiable along generalized geodesics.
Next, we apply the result of the previous proof to characterize the subdifferential of the regularized energies.