A blob method for diffusion

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 $\rho$ is a curve in the space of probability measures, are formally Wasserstein gradient flows of the energy

where $\mathcal{L}^{d}$ is $d$-dimensional Lebesgue measure. This implies that solutions $\rho(t,x)$ of (1) satisfy

for a generalized notion of gradient $\nabla_{W_{2}}$, which is formally given by

where $\delta\mathcal{E}/\delta\rho$ is the first variation density of $\mathcal{E}$ at $\rho$ (c.f. ).

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 aggregation-diffusion equations , stability of steady states of the Keller–Segel equation , and uniqueness of bounded solutions . The underlying gradient flow theory has been well developed in the case of convex (or, more generally, semiconvex) energies , and more recently, is being extended to consider energies with more general moduli of convexity .

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 . 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 . A wide variety of strategies have been developed for this approach: working with different discretizations of the space of Lagrangian maps , using alternative formulations of the variational structure , making use of convex analysis and computational geometry to solve the optimality conditions , and many others .

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 $\rho_{0}$ as a finite sum of $N$ Dirac masses, that is,

and solving the partial differential equation (1) reduces to solving a system of ordinary differential equations for the locations of the Dirac masses,

The particle solution $\rho^{N}(t)$ is the Wasserstein gradient flow of the energy (2) with initial data $\rho_{0}^{N}$, 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 $\rho^{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 . 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 . 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 $\varphi_{\varepsilon}$. In terms of the Wasserstein gradient flow perspective this translates into regularizing the interaction energy $(1/2)\int(W*\rho)\,d\rho$ as $(1/2)\int(W*\varphi_{\varepsilon}*\rho)\,d\rho$.

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 , and several recent works have considered analogous methods for equations of Wasserstein gradient flow type, including the Keller–Segel equation . 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 . Degond and Mustieles considered the case of linear diffusion ($m=1$) by interpreting the Laplacian as induced by a velocity field $v$, $\Delta\rho=\nabla\cdot(v\rho)$, $v=\nabla\rho/\rho$, and regularizing the numerator and denominator separately by convolution with a mollifier . 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 and extended to the case of the porous medium equation by Oelschläger and Lions and Mas-Gallic . 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. Oelschläger , on the other hand, succeeded in proving convergence of the deterministic particle method, as long as the corresponding solution of the porous medium equation is smooth and positive. 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 . 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 . Other approaches that respect the method’s variational structure have been recently proposed in one dimension by approximating particles by non-overlapping blobs . For further background on deterministic particle methods, we refer the reader to Chertock’s comprehensive review .

For more general nonlinear diffusion, we define

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 $\mathcal{U}_{\varepsilon}$ gives

which does not allow for a complete discretization of the integrals. On the contrary, in the second case, all convolutions involve $\rho$, so a similar computation (as it can be found in the proof of Corollary 5.5) shows that they 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 . 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 $\mathcal{F}_{\varepsilon}$ $\Gamma$-converges to $\mathcal{F}$ as $\varepsilon\to 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 $\Gamma$-convergence in hand, in Section 5 we then turn to the question of convergence of gradient flows, restricting to the case $m\geq 2$. Using the framework introduced by Sandier and Serfaty , we prove that, under sufficient regularity assumptions, gradient flows of the regularized energies converge as $\varepsilon\to 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(\varepsilon)$; 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^{\infty}$-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\sim 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=\left|\cdot\right|^{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\chi\log\left|\cdot\right|,\chi>0,m=1,2$) and the original two-dimensional equation ($V=0$, $W=\Delta^{-1}$, $m=1$). The one-dimensional equation has a critical mass $1$, and the two-dimensional equation has critical mass $8\pi$, at which point the concentration effects from the nonlocal interaction term balance with linear diffusion ($m=1$) . 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 $\varepsilon$ are refined.

There are several directions for future work. Our convergence theorem for $m\geq 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 $\rho^{m}$, which should remain bounded by the gradient flow structure; see equations (24) and (25). When $1\leq 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 $\varepsilon\to 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 $\varepsilon$. In particular, a theoretical result on the optimal scaling relation between $h$ and $\varepsilon$ remains open, though we observe good numerical performance for $\varepsilon=h^{1-p}$, $0<p\ll 1$. In a less technical direction, we foresee a use of the presented ideas in conjunction with splitting schemes for certain nonlinear kinetic equations , as well as in the fluids , since our numerical results demonstrate comparable rates of convergence to the particle strength exchange method, which has already gained attention in these contexts .

Preliminaries

2. 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.

whenever the integral converges. We consider mollifiers $\varphi$ satisfying the following assumption.

This assumption is satisfied by both Gaussians and smooth functions with compact support. Assumption 2.1 also ensures that $\varphi$ has finite first moment. For any $\varepsilon>0$, we write

Likewise, the technical assumption that $\varphi=\zeta*\zeta$, and therefore that $\varphi_{\varepsilon}=\zeta_{\varepsilon}*\zeta_{\varepsilon}$, allows us to regularize integrands involving the mollifier $\varphi_{\varepsilon}$; indeed, the following lemma provides sufficient conditions for moving functions in and out convolutions with mollifiers within integrals. (See also for a similar result.) This is an essential component in the proofs of both main results, Theorems 4.1 and 5.8, on the the $\Gamma$-convergence of the regularized energies and the convergence of the corresponding gradient flows. See Appendix A for the proof of this lemma.

3. Optimal transport, Wasserstein metric, and gradient flows

We now describe basic facts about optimal transport, including the Wasserstein metric and associated gradient flows. (See also for further background and more details on the definitions and remarks found in this section.)

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.

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.

The Wasserstein metric is formally Riemannian, and we may define the tangent space as follows.

We now turn to the definition of a gradient flow in the Wasserstein metric (c.f. [3, Proposition 8.3.1, Definition 11.1.1]).

We close this section with the following definition of the Wasserstein local slope.

When the functional $\mathcal{G}$ in Definition 2.9 is in addition semiconvex along geodesics the local slope $|\partial\mathcal{G}|$ is a strong upper gradient for $\mathcal{G}$. In this case a gradient flow of $\mathcal{G}$ is characterized as being a $2$-curve of maximal slope with respect to $|\partial\mathcal{G}|$; see [3, Theorem 11.1.3].

Regularized internal energies

The foundation of our blob method is the regularization of the internal energy $\mathcal{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 $\mathcal{F}_{\varepsilon}$, including convexity, lower semicontinuity, and differentiability. In what follows, we will suppose that our internal energies satisfy the following assumption.

Suppose $F\in C^{2}(0,+\infty)$ satisfies $\lim_{s\to+\infty}F(s)=+\infty$ and either $F$ is bounded below or $\liminf_{s\to 0}F(s)/s^{\beta}>-\infty$ for some $\beta>-2/(d+2)$. Suppose further that $U(s)=sF(s)$ is convex, bounded below, and $\lim_{s\to 0}U(s)=0$.

Thanks to this assumption we can define the internal energy corresponding to $F$ by

Assumption 3.1 implies that $F$ is nondecreasing. Indeed, by the convexity of $U(s)$ and the fact that $\lim_{s\to 0}sF(s)=0$,

which leads to $F^{\prime}(s)\geq 0$ for all $s\in(0,\infty)$.

Our assumption does not ensure that $\mathcal{F}$ is convex along Wasserstein geodesics, unless $F$ is convex.

McCann’s condition on the internal density $U$ for the convexity of the internal energy $\mathcal{F}$ can be stated on the function $F$ instead: the function $s\mapsto F(s^{-d})$ is nonincreasing and convex on $(0,\infty)$, i.e.,

which, by Remark 3.2, holds when for example $F$ is convex and satisfies Assumption 3.1.

We regularize the internal energies by convolution with a mollifier.

An important class of internal energies satisfying Assumption 3.1 are given by the (negative) entropy and Rényi entropies.

The entropy and Rényi entropies, and their regularizations, are given by

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\geq 1$, then we denote $\mathcal{E}$ by $\mathcal{E}^{m}$ and $\mathcal{E}_{\varepsilon}$ by $\mathcal{E}_{\varepsilon}^{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, $\varphi_{\varepsilon}$ 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.

Although the regularized energy in Definition 3.4 is of a novel form, it was noticed in [72, Proposition 6.9] that a previous particle method for diffusive gradient flows leads to a similar regularized internal energy after space discretization . 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 [31, Lemma 4.1].

where $C_{\varepsilon}=C_{\varepsilon}(m,\mu)\to 0$ as $\varepsilon\to 0$. Furthermore, for all $\delta>0$, we have

For all $\varepsilon>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 convergence of maps with respect to varying measures; see Proposition B.2. For $1<m\leq 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 $\varphi$ 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.

A key consequence of the preceding proposition is that the regularized energies are semiconvex along generalized geodesics, as we now show.

Suppose $F$ satisfies Assumption 3.1 and is convex. Then $\mathcal{F}_{\varepsilon}$ is $\lambda_{F}$-convex along generalized geodesics, where

We now use the previous results to characterize the subdifferential of the regularized internal energy. The structure of argument is classical (c.f. ), but due to the novel form of our regularized energies, we include the proof in Appendix A.

Suppose $F$ satisfies Assumption 3.1 and is convex. Let $\varepsilon>0$ and $\mu\in D(\mathcal{F}_{\varepsilon})$. Then

As a consequence of this characterization of the subdifferential, we obtain the analogous result for the full energy $\mathcal{E}_{\varepsilon}$, as in Definition 3.6. See Appendix A for the proof.

ΓΓ\Gamma-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 $\mathcal{F}_{\varepsilon}^{m}$ for $m\geq 1$. We begin by showing that $\mathcal{F}_{\varepsilon}^{m}$ $\Gamma$-converges to $\mathcal{F}$ as $\varepsilon\to 0$ with respect to the weak-∗ topology.

If $\mu_{\varepsilon}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu$, we have $\liminf_{\varepsilon\to 0}\mathcal{F}^{m}_{\varepsilon}(\mu_{\varepsilon})\geq\mathcal{F}^{m}(\mu)$.

We have $\limsup_{\varepsilon\to 0}\mathcal{F}^{m}_{\varepsilon}(\mu)\leq\mathcal{F}^{m}(\mu)$.

We begin by showing the result for $1\leq m\leq 2$, in which case the function $F$ is concave. We first show part (i). By Proposition 3.8, for all $\varepsilon>0$,

By Lemma 2.3, $\mu_{\varepsilon}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu$ implies $\zeta_{\varepsilon}*\mu_{\varepsilon}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu$. Therefore, by the lower semicontinuity of $\mathcal{F}^{m}$ with respect to weak-∗ convergence [3, Remark 9.3.8],

which gives the result. We now turn to part (ii). Again, by Proposition 3.8, for all $\varepsilon>0$,

where $C_{\varepsilon}\to 0$ as $\varepsilon\to 0$. Therefore, $\limsup_{\varepsilon\to 0}\mathcal{F}^{m}_{\varepsilon}(\mu)\leq\mathcal{F}^{m}(\mu)$.

so $(\zeta_{\varepsilon_{n}}*\mu_{\varepsilon_{n}})_{n}$ converges weakly in $L^{2}$ to $\mu$. By the Banach–Saks theorem (c.f. [74, Section 38]), up to taking a further subsequence of $(\zeta_{\varepsilon_{n}}*\mu_{\varepsilon_{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, Theorem 5.4.4(ii)],

Combining this with equation (18), we obtain

Now, we add a confining drift or interaction potential to our internal energies, so that energy minimizers exist and we may apply the previous $\Gamma$-convergence result to conclude that minimizers converge to minimizers. For the remainder of the section we consider energies of the form $\mathcal{E}_{\varepsilon}^{m}$ given in Definition 3.6, with the following additional assumptions on $V$ and $W$ to ensure that the energy is confining.

The potentials $V$ and $W$ are bounded below and one of the following additional assumptions holds:

Under these assumptions, the regularized energies $\mathcal{E}_{\varepsilon}^{m}$ are lower semicontinuous with respect to weak-∗ convergence ($m>1$) and Wasserstein convergence ($m=1$), where for the latter we assume $\varphi$ is a Gaussian (c.f. Proposition 3.9, and [3, Lemma 5.1.7], [66, Lemma 3.4] and [80, Lemma 2.2]).

We now prove existence of minimizers of $\mathcal{E}_{\varepsilon}^{m}$, for all $\varepsilon>0$.

First suppose $m>1$, so that $\mathcal{F}_{\varepsilon}\geq 0$ and $\mathcal{E}^{m}_{\varepsilon}$ is bounded below. By Remark 4.3, if (CV) holds, then any minimizing sequence of $\mathcal{E}^{m}_{\varepsilon}$ has a subsequence that converges in the weak-∗ topology. Likewise, if (CW) holds, then any minimizing sequence of $\mathcal{E}^{m}_{\varepsilon}$ has a subsequence that, up to translation, converges in the weak-∗ topology. By lower semicontinuity of $\mathcal{E}^{m}_{\varepsilon}$, the limits of minimizing sequences are minimizers of $\mathcal{E}^{m}_{\varepsilon}$.

Finally, we conclude that minimizers of the regularized energy converge to minimizers of the unregularized energy.

The proof is classical. We include it for completeness.

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 $\Gamma$-convergence of gradient flows, as we shall see in the next section.

ΓΓ\Gamma-convergence of gradient flows

We now consider gradient flows of the regularized energies $\mathcal{E}_{\varepsilon}^{m}$, as in Definition 3.6, for $m\geq 2$ and prove that, under sufficient regularity assumptions, gradient flows of the regularized energies converge to gradient flows of the unregularized energy as $\varepsilon\to 0$. For simplicity of notation, we often write $\mathcal{E}_{\varepsilon}^{m}$ and $\mathcal{F}_{\varepsilon}^{m}$ for $\varepsilon\geq 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.

More generally, our results naturally extend to drift and interaction energies that are merely $\omega$-convex; see . 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.

Let $\varepsilon\geq 0$ and $m\geq 2$. Suppose $\mathcal{E}^{m}_{\varepsilon}$ is as in Definition 3.6 and $V$ and $W$ satisfy Assumption 5.1. Then, for any $\mu_{0}\in\overline{D(\mathcal{E}_{\varepsilon}^{m})}$, there exists a unique gradient flow of $\mathcal{E}_{\varepsilon}^{m}$ with initial datum $\mu_{0}$.

For $\varepsilon>0$, Proposition 3.9 ensures that $\mathcal{F}_{\varepsilon}^{m}$ is lower semicontinuous with respect to weak-∗ convergence, hence also $2$-Wasserstein convergence. For $\varepsilon=0$, the unregularized internal energy $\mathcal{F}^{m}$ is also lower semicontinuous with respect to weak-∗ and $2$-Wasserstein convergence [66, 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, $\mathcal{E}_{\varepsilon}^{m}$ is lower semicontinuous for all $\varepsilon\geq 0$.

In the case $\varepsilon=0$, gradient flows of the energies $\mathcal{E}^{m}$ are characterized as solutions of the partial differential equation (1); c.f. [3, Theorems 10.4.13 and 11.2.1], [25, Theorem 2.12]. Now, we show that gradient flows of the regularized energies $\mathcal{E}_{\varepsilon}^{m}$ can also be characterized as solutions of a partial differential equation.

A consequence of the previous proposition is that, for the regularized energies $\mathcal{E}^{m}_{\varepsilon}$, 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.

To see that (23) is well-posed, first note that the function

is Lipschitz. Likewise, Assumption 5.1 ensures $y_{i}\mapsto\nabla V(y_{i})$ and $y_{i}\mapsto\sum_{j\in I}\nabla W(y_{i}-y_{j})$ are continuous and one-sided Lipschitz. Therefore, the ODE system (23) is well-posed forward in time.

Multiplying both sides by $m_{i}$, summing over $i$, and taking $\mu_{\varepsilon}=\sum_{i\in I}\delta_{X_{i}(\cdot)}m_{i}$ for $t\in[0,T]$ gives

for $v$ as in (22). Therefore, $\mu_{\varepsilon}$ is a weak solution of the continuity equation with velocity field $v$. Furthermore, for all $T>0$

We now turn to the $\Gamma$-convergence of the gradient flows of the regularized energies, using the scheme introduced by Sandier–Serfaty and then generalized by Serfaty , 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.)

$\displaystyle\liminf_{\varepsilon\to 0}\int_{0}^{t}|\mu_{\varepsilon}^{\prime}|(s)^{2}\,ds\geq\int_{0}^{t}|\mu^{\prime}|(s)^{2}\,ds$,

$\displaystyle\liminf_{\varepsilon\to 0}\mathcal{E}^{m}_{\varepsilon}(\mu_{\varepsilon}(t))\geq\displaystyle\mathcal{E}^{m}(\mu(t))$,

$\displaystyle\liminf_{\varepsilon\to 0}\int_{0}^{t}|\partial\mathcal{E}^{m}_{\varepsilon}|^{2}(\mu_{\varepsilon}(s))\,ds\geq\int_{0}^{t}|\partial\mathcal{E}^{m}|^{2}(\mu(s))\,ds$.

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 $\Gamma$-convergence of the gradient flows, a key assumption is that the following quantity is bounded uniformly in $\varepsilon>0$ along the gradient flows $\mu_{\varepsilon}$ of the regularized energies $\mathcal{E}^{m}_{\varepsilon}$:

Still, $\|\mu_{\varepsilon}\|_{BV_{\varepsilon}^{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 $\|\mu_{\varepsilon}\|_{BV_{\varepsilon}^{m}}$ as a nonlocal approximation of the $L^{1}$-norm of the gradient of $\mu^{m}$.

We begin with a technical lemma we shall use to prove the convergence of the gradient flows.

where $C_{\zeta}>0$ is as in Assumption 2.1.

First, we consider $I_{1}$. Since, in the integral, $|x-y|<\varepsilon^{\bar{r}}$, we obtain

Since $0\leq(m-2)/(m-1)<1$, Jensen’s inequality gives

With this technical lemma in hand, we now turn to the $\Gamma$-convergence of the gradient flows.

for some $\mu(0)\in D(\mathcal{E}^{m})$. Furthermore, suppose that the following hold:

$\sup_{\varepsilon>0}\int_{0}^{T}\|\mu_{\varepsilon}(t)\|_{BV_{\varepsilon}^{m}}dt<\infty$;

It remains to verify conditions (S0), (S1), (S2), and (S3) from Theorem 5.6. Item (S0) holds by assumption (A0). Item (S1) follows by the same argument as in [38, Theorem 5.6]. Item (S2) is an immediate consequence of the fact that $\mu_{\varepsilon}(t)\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu(t)$ for almost every $t\in[0,T]$, our main $\Gamma$-convergence Theorem 4.1, and the lower semicontinuity of the potential and interaction energies with respect to weak-∗ convergence [3, Lemma 5.1.7].

We devote the remainder of the proof to showing Condition (S3). We shall use the following fact throughout: combining Assumption (A2) with Proposition 3.8 implies that

To prove (S3) we may assume, without loss of generality, that $\liminf_{\varepsilon\to 0}\int_{0}^{T}|\partial\mathcal{E}^{m}_{\varepsilon}|(\mu_{\varepsilon}(t))^{2}\,dt$ is finite, so by Fatou’s lemma

so $\liminf_{\varepsilon\to 0}|\partial\mathcal{E}^{m}_{\varepsilon}|(\mu_{\varepsilon}(t))<\infty$ for almost every $t\in[0,T]$. In particular, up to taking subsequences, we may assume that, for almost every $t\in[0,T]$, $\{|\partial\mathcal{E}^{m}_{\varepsilon}|(\mu_{\varepsilon}(t))\}_{\varepsilon}$ is bounded uniformly in $\varepsilon>0$. By Corollary 3.13,

when (31) holds for almost every $t\in[0,T]$. Furthermore, the inequality in (32) is, by Proposition B.2(ii), a consequence of

First, we address the terms with the drift and interaction potentials $V$ and $W$. Combining Assumption 5.1 on $V$ and $W$ with Assumption (A5.8) on $\mu_{\varepsilon}$ ensures that $|\nabla V|$ is uniformly integrable in $d\mu_{\varepsilon}\otimes d\mathcal{L}^{d}$ and $(x,y)\mapsto|\nabla W(x-y)|$ is uniformly integrable $d\mu_{\varepsilon}\otimes d\mu_{\varepsilon}\otimes d\mathcal{L}^{d}$.Therefore, by [3, Lemma 5.1.7], $(\mu_{\varepsilon})_{\varepsilon}$ converging weakly-∗ to $\mu$ ensures that

Recalling the abbreviation $p_{\varepsilon}:=(\varphi_{\varepsilon}*\mu_{\varepsilon})^{m-2}\mu_{\varepsilon}$, we rewrite the inner integral on the left-hand side of (34) as

Applying Lemma 5.7 together with (29) and (A3), and integrating by parts, we obtain

Now we move $\nabla f$ out of the convolution. By Lemma 2.2, there exists $p>0$ so

where we again use (27). Using the inequality in (28) and that $\{\int_{0}^{T}\mathcal{F}^{m}_{\varepsilon}(\mu_{\varepsilon}(t))\,dt\}_{\varepsilon}$ is uniformly bounded in $\varepsilon$,

Taking $g=\nabla f$, choosing $R>1$ so that $\kappa_{R}\equiv 1$ on the support of $\nabla f$, and combining the above equation with equation (35), we obtain

for all $\varepsilon>0$. Combining this with (38) gives

By another application of Hölder’s inequality, this guarantees

and $\mu$ is the gradient flow of $\mathcal{E}^{m}$ with initial data $\mu(0)$.

The above theorem generalizes a result by Lions and Mas-Gallic on a numerical scheme for the porous medium equation $\partial_{t}\mu=\Delta\mu^{2}$ on a bounded domain with periodic boundary conditions to equations of the form (1) on Euclidean space.

Furthermore, since the energy $\mathcal{F}^{2}_{\varepsilon}$ decreases along solutions to the gradient flow, we have

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é [65, 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_{\varepsilon>0}\mathcal{F}^{1}(\mu_{\varepsilon}(0))$ so that, for all $t\in[0,T]$,

Since $D^{2}V$ and $D^{2}W$ are bounded, $|\nabla V|$ and $|\nabla W|$ grow at most linearly. Consequently, there exists $C^{\prime}>0$, depending on $V$, $W$, and $C_{\kappa}$ so that

Likewise, by Lemma 5.7, there exists $r>0$ so that, for all $t\in[0,T]$,

Therefore, there exists $C^{\prime\prime}>0$ so that, for all $t\in[0,T]$,

As the right-hand side is independent of $R>1$, by sending $R\to+\infty$ by the dominated convergence theorem we obtain that for $\varepsilon^{r}<1/(2C^{\prime\prime})$,

We may combine this with the inequality in (43) to obtain, for all $t\in[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 (A5.8) is a consequence of the inequality in (45). Assumption (A1) is a consequence of the inequalities in (42), (44) and (46).

Numerical results

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 $\varepsilon>0$, define the evolution of these measures by

where $\{X_{i}(t)\}_{i\in Q_{R}^{h}}$ are solutions to the ODE system (23) on a time interval $[0,T]$ with initial data $X_{i}(0)=ih$. If $h=o(\varepsilon)$ as $\varepsilon\to 0$ and Assumptions (A5.8)–(A2) from Theorem 5.8 hold, then $(\mu_{\varepsilon}(t))_{\varepsilon}$ converges in the weak-∗ topology to $\mu(t)$ as $\varepsilon\to 0$ for almost every $t\in[0,T]$, where $\mu(t)$ is the unique solution of (1) with initial datum $\mu(0)$.

so $\mu_{\varepsilon}(0)\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu(0)$ as $\varepsilon\to 0$ (and so, as $h\to 0$). Likewise, for all $\varepsilon,h>0$, ${\mathop{\rm supp\ }}\mu_{\varepsilon}(0)\subseteq B_{R}(0)$. Consequently, since $V$ and $W$ are continuous,

By Theorem 4.1, we have that $\liminf_{\varepsilon\to 0}\mathcal{F}^{m}_{\varepsilon}(\mu_{\varepsilon}(0))\geq\mathcal{F}^{m}(\mu_{\varepsilon}(0)).$ By Proposition 3.8, for all $\varepsilon>0$ 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,\|\nabla\zeta\|_{\infty}$, and the space dimension. Therefore, provided that $h=o(\varepsilon)$ as $\varepsilon\to 0$, we obtain that $\zeta_{\varepsilon}*\mu_{\varepsilon}-\zeta_{\varepsilon}*\mu\to 0$ in $L^{m}$. ∎

In Theorem 6.1, we proved that, as long as Assumptions (A5.8)–(A2) from Theorem 5.8 hold along the particle solutions $\{\mu_{\varepsilon}\}_{\varepsilon}$, 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 (A5.8) is that the $(m-1)$th moment of the particle solution

is bounded uniformly in $t,\varepsilon$, and $h$. In particular, this is satisfied if the particles remain compactly supported in a ball.

A sufficient condition for Assumption (A1) is that

with $p_{\varepsilon}=(\varphi_{\varepsilon}*\mu_{\varepsilon})^{m-2}\mu_{\varepsilon}$, remains bounded uniformly in $t$, $\varepsilon$, and $h$. In fact, for purely diffusive problems, we observe that this quantity is not only bounded uniformly in $\varepsilon$ 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 $\varepsilon$ and $h$ and converges to the corresponding norm of the steady state as $t\to\infty$; see Figure 6 below.

A sufficient condition for Assumption (A2) is that the blob solution converges to a limit in $L^{1}$ and $L^{\infty}$, 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 (A5.8)–(A2) may be verified numerically in order to give confidence that the limit of any blob method solution is, in fact, the correct exact solution.

2. Numerical implementation

We now describe the details of our numerical implementation. In all of the numerical examples which follow, our mollifiers $\zeta_{\varepsilon}$ and $\varphi_{\varepsilon}$ 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 $\mu_{\varepsilon}(t)$ as described in Theorem 6.1. As a mild simplification, we consider the mass of each particle to be given by $m_{i}=\mu(0,ih)h^{d}$, where $\mu(0,ih)$ is the value of the initial datum $\mu(0)$ at the grid point $ih$. For the numerical examples we consider, in which $\mu(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 , 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 . In particular, we compute the evolution of the particle trajectories via the SciPy implementation of the Fortran VODE solver , 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(\varepsilon)$ as $\varepsilon\to 0$. Numerically, we observe the fastest rate of convergence with $\varepsilon=h^{1-p}$, for $0<p\ll 1$, as $h\to 0$. Since computational speed decreases as $p$ approaches 0, we take $\varepsilon=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 .

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,

We measure the accuracy of our numerical method with respect to the $L^{1}$-, $L^{\infty}$-, and Wasserstein metrics. To compute the $L^{1}$- and $L^{\infty}$-errors, we take the difference between the exact solution and the blob solution (50) and evaluate discrete $L^{1}$- and $L^{\infty}$-norms using the following formulas:

We compute the Wasserstein distance between our particle solution $\mu_{\varepsilon}$ in (48) and the exact solution $\mu$ in one dimension using the formula

where $F_{\mu_{\varepsilon}}^{-1}$ and $F_{\mu}^{-1}$ are the generalized inverses of the cumulative distribution functions of $\mu$ and $\mu_{\varepsilon}$, respectively; c.f. [3, Theorem 6.0.2]. We evaluate the integral in (51) numerically using the SciPy implementation of the Fortran library QUADPACK . 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 .

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 $K=K(m,d)$ chosen so that $\int\psi_{m}(\tau,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 $\tau=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, $\rho_{0}(\cdot)=0.3\psi_{1}(\cdot+1,0.0225)+0.7\psi_{1}(\cdot-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\to 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 $m$th 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 (A1) 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^{\infty}$-norm and examine the scaling of the error with the grid spacing $h$. (Recall that $\varepsilon=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^{\infty}$-error scales quadratically for $m=1$, superlinearly for $m=2$, and sublinearly for $m=3$. This deterioration of the rate of $L^{\infty}$-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 7, we consider the one-dimensional variant of the Keller–Segel equation ($V=0$, $W(\cdot)=2\chi\log\left|\cdot\right|$, $m=1$) studied in . Its interest is that it has a defined critical value $\chi$ for unit mass leading to the dichotomy of blow-up versus global existence. For $\chi=1.5$ and initial data of mass one, solutions blow up in finite time. We consider initial data given by a Gaussian $\psi_{1}(\tau,\cdot)$, $\tau=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 . 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.

In Figure 8, we consider a nonlinear variant of the Keller–Segel equation ($V=0$, $W(\cdot)=2\chi\log\left|\cdot\right|$, $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(\cdot)=1/(2\pi)\log\left|\cdot\right|$, $m=1$) in two dimensions. In Figures 9, 10, and 11, the initial data is given by a Gaussian $\psi_{1}(\tau,\cdot)$, $\tau=0.16$, scaled to have mass that is either supercritical ($>8\pi$), critical ($=8\pi$), or subcritical ($<8\pi$) with respect to blowup behavior. In particular, for supercritical initial data, solutions blow up in finite time . 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.0\bar{3}$ and varying mass $7\pi,8\pi,$ and $9\pi$, 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.0\bar{3}$, 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\to 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\pi,8\pi,$ and $9\pi$. The slopes of the second moment agree well with the theoretically predicted slopes given in Figure 10.

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.

Thus, we conclude our result by estimating the above quantity by

We now give the proof that if $\mu_{\varepsilon}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu$, then $\varphi_{\varepsilon}*\mu_{\varepsilon}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu$.

Since $\mu_{\varepsilon}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu$, the second term goes to zero. We bound the first term as follows:

which goes to zero as $\varepsilon\to 0$. ∎

Next, we prove the inequalities relating the regularized internal energies to the unregularized internal energies.

We begin with (11). To prove the left inequality, we may assume without loss of generality that $\mu\in D(\mathcal{F})$. First, we show the result for the entropy ($m=1$). Note that

By Jensen’s inequality for the convex function $s\mapsto s\log s$, the relative entropy is nonnegative, which gives the result. Now, we show the left inequality in (11) for $1<m\leq 2$. By the above-the-tangent property of the concave function $F_{m}$ and Hölder’s inequality, we get

Now, we show (12). Since $F_{m}$ is convex for $m\geq 2$, this is simply a consequence of reversing the inequalities in the last two inequalities.

Finally, we consider the lower bounds (13). When $m=1$, these follow from the right inequality in (11), a Carleman-type estimate [31, Lemma 4.1] ensuring that $\mathcal{F}_{\varepsilon}^{m}(\zeta_{\varepsilon}*\mu)\geq-(2\pi/\delta)^{d/2}-\delta M_{2}(\zeta_{\varepsilon}*\mu)$ for all $\delta>0$, and the fact that

When $m>1$, we simply use that $F_{m}\geq 0$. ∎

We now give the proof that, for all $\varepsilon>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 $\varphi$ to be a Gaussian.

We now show (i). Suppose $\mu_{n}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\mu$. By Lemma B.3, we have

Combining the two previous inequalities, we obtain $\liminf_{n\to\infty}\mathcal{F}^{m}_{\varepsilon}(\mu_{n})\geq\mathcal{F}^{m}_{\varepsilon}(\mu)$, giving the result.

Define $f_{n}:=\log(\varphi_{\varepsilon}*\mu_{n})$ and $q(\cdot):=C_{0}|\cdot-x_{0}|^{2}+C_{1}$. Then, by Lemma B.3, we have

Since $\mu_{n}\to\mu$ in the Wasserstein metric,

Furthermore, by (53) and the fact that $\log(\cdot)$ is continuous on $(0,+\infty)$,

Thus, combining (55), (56), and (57), we obtain,

Now we turn to the proof that the regularized energies are differentiable along generalized geodesics.

where $c_{s,\alpha}(y,z)=(1-s)\varphi_{\varepsilon}*\mu_{1}(y)+s\varphi_{\varepsilon}*\mu_{\alpha}^{2\to 3}((1-\alpha)y+\alpha x)$. Using Taylor’s theorem compute

where $D_{\alpha}(y,z,v,w)$ is a term depending on the Hessian of $\varphi_{\varepsilon}$ satisfying

Hence, since $F^{\prime}$ is nondecreasing,

Thus, to complete the result, it suffices to show that there exists $g\in L^{1}(\gamma\otimes\gamma)$ so that

since the result then follows by the dominated convergence theorem. Since $F^{\prime}$ is nondecreasing we may take

Next, we apply the result of the previous proof to characterize the subdifferential of the regularized energies.

which shows that $f$ is nondecreasing, and so $f(1)\geq\lim_{\alpha\to 0}f(\alpha)$, which implies (after integrating against $d\gamma(x,y)$)

Then, by (58) and antisymmetry of $\nabla\varphi_{\varepsilon}$, compute

Now compute, using the antisymmetry of $\nabla\varphi_{\varepsilon}$,

where passing the limit $\alpha\to 0$ inside the integral in the first line is justified by the fact that $H^{\prime}$ is bounded. Then, by the definition of the local slope of $\mathcal{F}_{\varepsilon}$,

since, by definition of the 2-Wasserstein distance,

which shows the desired result. Since $|\partial\mathcal{F}_{\varepsilon}|(\mu)$ is the unique minimal norm element of $\partial\mathcal{F}_{\varepsilon}$, this also shows that we actually have equality in the right-hand side above.

Combining this with Proposition 3.10, we obtain

Rewriting the expression from equation (14) gives

Finally, we prove the characterization of the subdifferential of the full regularized energies $\mathcal{E}^{m}_{\varepsilon}$.

where $\mu_{0}$, $\mu_{1}$, $\lambda$ and $\xi$ are as in the proof of Proposition 3.12. ∎

Appendix B Weak convergence of measures

In this appendix, we recall several fundamental results on the weak convergence of measures. We begin with a result due to Ambrosio, Gigli, and Savaré on convergence of maps with respect to varying probability measures. This plays a key role in our proofs of both the $\Gamma$-convergence of the energies and the $\Gamma$ convergence of the gradient flows.

Furthermore, we say that $(v_{n})_{n}$ converges strongly to $v$ in $L^{p}$, $p>1$, if

We close by recalling a generalization of Fatou’s lemma, for varying measures.

Acknowledgments: The authors thank Andrew Bernoff, Andrea Bertozzi, Eric Carlen, Yanghong Huang, Inwon Kim, Dejan Slepčev, and Fangbo Zhang for many helpful discussions.

References