# Convergence of the fully discrete JKO scheme

Anastasiia Hraivoronska, Filippo Santambrogio

► To cite this version:

Anastasiia Hraivoronska, Filippo Santambrogio. Convergence of the fully discrete JKO scheme. 2025. hal-05035393

HAL Id: hal-05035393

https://hal.science/hal-05035393v1

Preprint submitted on 15 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# CONVERGENCE OF THE FULLY DISCRETE JKO SCHEME

ANASTASIIA HRAIVORONSKA AND FILIPPO SANTAMBROGIO

**ABSTRACT.** The JKO scheme provides the discrete-in-time approximation for the solutions of evolutionary equations with Wasserstein gradient structure. We study a natural space-discretization of this scheme by restricting the minimization to the measures supported on the nodes of a regular grid. The study of the fully discrete JKO scheme is motivated by the applications to developing numerical schemes for the nonlinear diffusion equation with drift and the crowd motion model. The main result of this paper is the convergence of the scheme as both the time and space discretization parameters tend to zero in a suitable regime.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>1.1. The JKO scheme and the evolution PDEs that it approximates</td><td>2</td></tr><tr><td>1.2. Contribution</td><td>3</td></tr><tr><td>1.3. Related work</td><td>5</td></tr><tr><td>1.4. Numerical challenges and applications</td><td>6</td></tr><tr><td>2. Setting of the problem and statement of the main results</td><td>7</td></tr><tr><td>2.1. Continuous models</td><td>7</td></tr><tr><td>2.2. Discretization</td><td>10</td></tr><tr><td>2.3. Main results</td><td>12</td></tr><tr><td>3. Evolution driven by a potential</td><td>14</td></tr><tr><td>3.1. Relation between h and τ</td><td>14</td></tr><tr><td>3.2. Convergence of the fully discrete JKO with potential</td><td>15</td></tr><tr><td>4. Nonlinear diffusion</td><td>17</td></tr><tr><td>4.1. Interpolations and corresponding inequalities</td><td>18</td></tr><tr><td>4.2. Compactness results</td><td>27</td></tr><tr><td>4.3. Convergence results</td><td>32</td></tr><tr><td>5. Crowd motion model</td><td>36</td></tr><tr><td>References</td><td>44</td></tr></table>

## 1. INTRODUCTION

The theory of gradient flows in the space of probability measures provides powerful tools for studying a wide class of evolutionary equations [3]. The JKO scheme (or the minimizing movement scheme) plays an essential role in this theory, acting as a counterpart for the implicit Euler scheme for classical gradient-flow ODEs. It first appeared in Jordan, Kinderleherer, and Otto's paper [31], where they recognized the Fokker-Planck dynamics as a gradient flow of the Boltzmann entropy in the space of probability measures with an appropriate metric. The JKO scheme provides as an output a sequence of probability measures (see below), each one optimizing a variational problem involving a transport cost
to the previous one and an energy functional (the Boltzmann entropy in the case of Fokker-Planck), which converges to the solution of the desired PDE when a time step parameter $\tau$ tends to 0. As such, it is a means to study well-posedness and to develop numerical schemes for the equations of interest.

In this paper, we consider the most natural space-discretization of this scheme: we consider a uniform and regular grid and restrict our analysis to measures supported on the nodes of this grid. While the optimal transport problem and the Wasserstein distance that it defines are well-posed for atomic measures, our approach requires adapting the definition of some energy functionals to the same atomic measures, but this is very natural when we consider measures on a regular and fixed grid. We need hence to introduce another approximation parameter $h$, corresponding to the space step of the grid, and we study the convergence of the scheme when $\tau, h \to 0$. We identify for this sake a sharp condition, namely $h/\tau \to 0$, which is both necessary and sufficient to prove this convergence. Note that this condition is the opposite of a very standard one in numerical analysis, called the CFL condition, where instead $\tau/h$ is required to be small (called the CFL condition). This is not surprising as the CFL condition is standard in explicit methods, while the JKO scheme is implicit.

1.1. **The JKO scheme and the evolution PDEs that it approximates.** The JKO scheme is a discrete-in-time iterative scheme defined in the space of probability measures $P(\Omega)$ (in this paper $\Omega \subset \mathbb{R}^d$ is a bounded domain):

$$ (\text{JKO}) \qquad \rho_{k+1}^{\tau} \in \arg \min_{\rho \in \mathcal{P}(\Omega)} \left\{ \mathcal{F}(\rho) + \frac{1}{2\tau} W_2^2(\rho, \rho_k^{\tau}) \right\}, \quad \text{with given } \rho_0^{\tau} = \bar{\rho} \in \mathcal{P}(\Omega), $$

where $\tau > 0$ is a time step parameter, $\mathcal{F}: P(\Omega) \to \mathbb{R} \cup \{+\infty\}$ is an energy functional, and $W_2$ is the Wasserstein distance induced by the optimal transport problem with the quadratic cost

$$ (1.1) \qquad W_2^2(\mu, \nu) = \min_{\gamma \in \Gamma(\mu, \nu)} \iint |x - y|^2 \gamma(dx dy), $$

where $\Gamma(\mu, \nu) = \{\gamma \in P(\Omega \times \Omega) : P_{\#}^{1}\gamma = \mu, P_{\#}^{2}\gamma = \nu\}$ is the set of couplings between the marginals $\mu$ and $\nu$. The sequence of minimizers $\{\rho_k^{\tau}\}$ of (JKO) can be used to build a continuous-in-time curve $\rho^{\tau}$, for instance, with piecewise constant interpolation

$$ \rho^{\tau}(t) := \rho_{k+1}^{\tau} \quad \text{for } t \in (k\tau, (k+1)\tau]. $$

Under suitable assumption on the energy $\mathcal{F}$, one can prove that $\rho^{\tau}$ converges to a weak solution of

$$ (1.2) \qquad \partial_t \rho + \nabla \cdot \left( \rho \nabla \frac{\delta \mathcal{F}}{\delta \rho} \right) = 0 \quad \text{in } (0, T) \times \Omega, $$

with the initial datum $\rho_0 = \bar{\rho}$ and complemented with the non-flux boundary conditions on $\partial\Omega$, where $\delta\mathcal{F}/\delta\rho$ denotes the first variation of $\mathcal{F}$.

After the seminal work [31] for the linear Fokker-Planck equation, the JKO scheme has been proven to be a useful analytical tool for many classes of PDEs [3], including general aggregation-diffusion equations on the Keller-Segel model [9], cross-diffusion systems [22, 26, 32], and higher-order equations [19, 38].

In this paper, we consider the energy functionals that are the sum of internal and potential energy

$$ (1.3) \qquad \mathcal{F}(\rho) = \int_{\Omega} f(\rho) + \int_{\Omega} \rho V, $$

with a convex energy density $f: (0, \infty) \to \mathbb{R} \cup \{+\infty\}$ and an external potential $V: \Omega \to \mathbb{R}$. This energy gives rise to the drift-diffusion equations of the form

$$ (1.4) \qquad \partial_t \rho = \Delta L(\rho) + \nabla \cdot (\rho \nabla V), $$

where the diffusion term is defined by $L'(s) = sf''(s)$.
We provide a few examples that we have in mind as test cases for the results of this paper. The most well-known PDE of the type (1.4) is the linear Fokker-Planck equation corresponding to $f(s) = s \log s$ and $L(s) = s$. The second example is the porous media equation with drift, where the energy density is the power law $f_m(s) = s^m/(m-1)$ and $L_m(s) = s^m, m > 1$. We use these two cases as running examples, but we prove the results for more general convex superlinear $f$ (the exact assumptions are specified in Section 2.3) and $V \in \text{Lip}(\mathbb{R}^d)$.

The Hele-Shaw-type system is the third example that plays a special role in this paper:

$$ (1.5) \qquad \begin{cases} \partial_t \rho = \nabla p + \nabla \cdot (\rho \nabla V), \\ p \ge 0, \quad \rho \le 1, \quad p(1-\rho) = 0, \end{cases} $$

where we introduced the pressure variable $p$. It emerges in the incompressible limit of the porous media equation, i.e., $m \to \infty$ [29]. The gradient flow interpretation for (1.5) appears in [39, 40] with application to the crowd motion model. The idea is that the solutions of (1.5) can be obtained by means of the JKO scheme (JKO) choosing the energy functional

$$ \mathcal{F}_{\text{CM}}(\rho) = \begin{cases} \int_{\Omega} \rho V, & \rho \le 1, \\ +\infty, & \text{otherwise.} \end{cases} $$

One can see that $\mathcal{F}_{\text{CM}}$ has a similar structure as (1.3) with the degenerate energy density $f(s) = 0$ for $s \in [0, 1]$ and $f(s) = +\infty$ for $s > 1$, which is exactly the pointwise limit of the porous-media energy density $f_m$ as $m \to \infty$.

1.2. **Contribution.** In this paper, we explore a fully discrete JKO scheme which, to our knowledge, has not been studied before. We restrict the minimization at each step of the JKO scheme to the set of probability measures supported on the vertices of the grid $h\mathbb{Z}^d$ with a small space discretization parameter $h > 0$. Precisely, we define the discrete domain $\mathcal{T}^h := h\mathbb{Z}^d \cap \Omega$ and the set of corresponding probability measures

$$ \mathcal{P}(\mathcal{T}^h) := \{ \rho \in \mathcal{P}(\Omega) : \rho = \sum_{z \in \mathcal{T}^h} \rho_z \delta_z \}. $$

Then the fully discrete JKO scheme is defined as

$$ (\text{JKO}^{h, \tau}) \qquad \rho_{k+1}^{h, \tau} \in \arg \min_{\rho \in \mathcal{P}(\mathcal{T}^h)} \left\{ \mathcal{F}^h(\rho) + \frac{1}{2\tau} W_2^2(\rho, \rho_k^{h, \tau}) \right\}, \quad \text{with given } \rho_0^{h, \tau} \in \mathcal{P}(\mathcal{T}^h), $$

where $\mathcal{F}^h : \mathcal{P}(\mathcal{T}^h) \to \mathbb{R} \cup \{+\infty\}$ is the discrete counterpart of the energy $\mathcal{F}$ from (1.3) defined as

$$ (1.6) \qquad \mathcal{F}^h(\rho) = \sum_{z \in \mathcal{T}^h} f\left(\frac{\rho_z}{h^d}\right) h^d + \sum_{z \in \mathcal{T}^h} V(z) \rho_z. $$

For the crowd motion case, we similarly define

$$ (1.7) \qquad \mathcal{F}_{\text{CM}}^{h}(\rho^{h}) = \begin{cases} \sum_{z \in \mathcal{T}^{h}} V(z) \rho_{z}, & \frac{\rho_{z}}{h^{d}} \le 1 \text{ for every } z, \\ +\infty, & \text{otherwise.} \end{cases} $$

We stress that we use the standard Wasserstein distance in the discrete scheme $(\text{JKO}^{h, \tau})$ differently from other works where the optimal transport model itself is modified in order to take care of the discrete structure of the measures (see for instance [37], where the distance is defined in connection with a Markov chain on a graph).
As we will explain better below, we introduce the fully discrete JKO scheme having applications for numerical schemes in mind. In this sense, the crowd motion model is a particularly interesting case, because (JKO<sup>h,τ</sup>) becomes a linear optimization problem:

$$ \rho_{k+1}^{h,\tau} = P_{\#}^{1} \gamma_{k}^{h,\tau}, \quad \gamma_{k}^{h,\tau} \in \arg \min_{\gamma \in \mathcal{P}(\mathcal{T}^h \times \mathcal{T}^h)} \left\{ \iint \left\{ V(x) + \frac{|x-y|^2}{2\tau} \right\} \gamma(dx dy), \quad P_{\#}^{1} \gamma \le 1, \ P_{\#}^{2} \gamma = \rho_k^\tau \right\}. $$

In the case of the drift-diffusion equation, (JKO<sup>h,τ</sup>) is a nonlinear convex optimization problem. We will also mention throughout the paper how to consider the case of interaction potentials (corresponding to a drift given by $-\nabla W * \rho$, or to an energy functional including a term of the form $\int \int W(x-y) d\rho(x) d\rho(y)$), which is a small variant in terms of PDEs, but in terms of optimization it is less interesting, as it could give raise to non-convex terms in the optimization problem. For this reason, the interaction case is only mentioned as a remark and is not the main scope of our contribution.

The goal of this paper is to establish the convergence result for (JKO<sup>h,τ</sup>) in the joint limit $h, \tau \to 0$. The relation between $h$ and $\tau$ plays a crucial role in the convergence. To see how this relation affects the convergence, we consider a toy example with the movement driven only by an external potential. Specifically, let the energy functional with $V \in C^{1,1}(\mathbb{R}^d)$ be defined as:

$$ (1.8) \qquad \mathcal{F}_V^h(\rho^h) := \sum_{z \in \mathcal{T}^h} V(z) \rho_z^h. $$

In this case, the Dirac masses in the initial datum move independently of each other. Consider the movement of $\delta_{x_0}$ with $x_0 \in \mathcal{T}^h$. In the continuous setting, it simply means that particles move along the trajectories of the ODE $\dot{x} = -(\nabla V)(x)$ with $x(0) = x_0$. In the discrete setting, we look for a position $x_1 \in \mathcal{T}^h$ minimizing $V(x) + |x - x_0|^2/2\tau$. Using $x_0$ as a competitor for the minimizer, we obtain

$$ V(x_0) \geq V(x_1) + \frac{|x_1 - x_0|^2}{2\tau}. $$

Since the minimal distance that each particle has to travel is $h$, a necessary condition in order to observe some movement is $h/\tau \le 2\|\nabla V\|_{\text{Lip}}$. We see that if asymptotically $h/\tau > 2\|\nabla V\|_{\text{Lip}}$, then we cannot expect convergence to the continuous solution, because every subsequent minimizer is equal to $x_0$ and the evolution is "frozen".

We will further see that, for convergence, we need an even stronger assumption $h/\tau \to 0$, and this assumption is sharp. If $h \sim \tau$, we cannot guarantee the convergence to the solution because the discrete grid can introduce additional drift to the evolution. The example with the potential movement is detailed in Section 3.

Strikingly enough, the condition $h/\tau \to 0$ is indeed also sufficient to prove the convergence when the functional $\mathcal{F}$ also includes a term depending on the density, which means when the corresponding PDE has diffusion terms. Our main results are thus the convergence statements for (JKO<sup>h,τ</sup>) to the drift-diffusion equation and the crowd motion model. Let $\{\rho_k^{h,\tau}\}$ be the family of minimizers for (JKO<sup>h,τ</sup>) with $\mathcal{F}^h$ defined in (1.6) and let $\rho^{h,\tau}$ be the piecewise constant-in-time interpolation of $\{\rho_k^{h,\tau}\}$. Then the curve $\rho^{h,\tau}$ converges (up to a subsequence) as $h, \tau, h/\tau \to 0$ to a continuous curve $\rho \in C([0, T], \mathcal{P}(\Omega))$, which is a weak solution to (1.4). An analogous result holds for the crowd motion energy defined in (1.7) with the convergence to (1.5). The precise statements of the main results are given in Section 2.3. The strategy of the proof is based on the so-called *Energy Dissipation Inequality*, which characterizes the solutions of the PDE as those curves satisfying a certain (sharp) inequality on the rate of dissipation of the energy $\mathcal{F}$. This inequality will be obtained by passing to the limit similar inequalities obtained on the discrete problems thanks to the optimality conditions, and carefully mixing different forms of interpolations. Indeed, a variational interpolation à la De Giorgi will be used to make a slope term
appear (see Section 4.), and this interpolation will be composed of measures also supported on the grid, but a gedoesic interpolation (in the sense of geodesic in the Wasserstein space on the whole space, not only on the grid) will be used to handle a velocity term.

1.3. **Related work.** The JKO scheme has been actively employed to develop numerical schemes for the evolution equations with the gradient structure. From the theoretical point of view, it has advantages such as a variational formulation, providing gradient-structure preserving and stability properties. The difficulty is finding an efficient approach to discretizing the Wasserstein distance. The existing methods rely on various reformulations of the optimal transport problem, such as dynamical formulation by Benamou-Brenier [6], entropic regularization [23], and semi-discrete formulation [33, 42]. We briefly review the results in the literature on numerical schemes based on the JKO scheme.

The idea of the Benamou-Brenier formulation is that the Wasserstein distance $W_2$ admits an alternative formulation

$$ (1.9) \quad W_2^2(\mu, \nu) = \min \left\{ \frac{1}{2} \int_0^1 \int_{\Omega} \left| \frac{d\bar{j}_t}{d\rho_t} \right|^2 d\rho_t dt, \quad \partial_t \rho + \operatorname{div} j = 0, \quad \rho_0 = \mu, \rho_1 = \nu \right\}. $$

Here, to compute $W_2$, one considers all absolutely continuous curves $[0, 1] \ni t \mapsto \rho_t$ connecting $\mu$ and $\nu$ and minimizes the corresponding integral of the velocity field $v_t = d\bar{j}_t/d\rho_t$. Relying on this formulation, they implement the augmented Lagrangian method in [7] and the successive primal-dual method in [21] for the numerical approximation of (JKO). A disadvantage of this approach is that solving the system of optimality conditions requires inner time stepping. As an alternative, which avoids the sub-discretization, several works use the linearization of the Wasserstein distance [16, 36, 43]. The main idea is that the subsequent minimizers $\rho_k^\tau$ and $\rho_{k+1}^\tau$ of (JKO) are "close" to each other, and it is reasonable to replace the $W_2$ distance with a suitable weighted $H^{-1}$ norm.

Another approach to computing the Wasserstein distance relies on modifying it with the entropic regularization term because the regularized version can be computed efficiently with the Sinkhorn algorithm [23]. The entropic regularization is applied to the JKO scheme in [18], where the authors prove the convergence of the regularized scheme in the joint limit $\tau, \varepsilon \to 0$ with $\varepsilon/\tau \to 0$, where $\varepsilon > 0$ is the regularization parameter.

In [8], the authors treat the JKO scheme as a variational problem on the space of convex functions with the Monge-Ampère operator, for which they construct a discretization.

The crowd motion model corresponding to the constraint $\rho \le 1$ has also been treated with numerical methods inspired by the JKO scheme, often using variants of the same methods as for non-linear diffusion (see, in particular, [7] or [18]). In the original works on the gradient flow interpretation of this crowd motion model ([39, 40]) a different specific approach was used in order to compute the projection of a measure onto the set of admissible measures satisfying $\rho \le 1$ (and the drift was treated via a splitting method, i.e. first advecting the density ignoring the constraint, and then projecting). The projection was approximated via a stochastic procedure reminiscent of *Diffusion Limited Aggregation*, which was well efficient in practice but a convergence proof was out of reach.

Finally, we cite the works in [34] and [20], where both linear diffusion and crowd motion models are addressed. In spirit this is very different from what we do, since it is based on a discretization in “space” (replacing densities with moving particles) but not in time, while in our approach the starting point is the JKO scheme, which is a discretization in time, which we complete with a space discretization. Yet, a common point with our work is that they require a suitable extension of some energy functionals (such as the Boltzmann entropy or the constraint on the density) to atomic measures. Here we see how our work really chooses the simplest and most natural point of view: having a fixed grid allows us without ambiguity to spread the mass of an atom over the cell around it, while when particles move
this is more delicate. This is done via a Moreau-Yosida approximation of the functional in [34] and
via convolution in [20].

1.4. **Numerical challenges and applications.** Even if this paper is not devoted to numerical computations, but only to the proof of convergence of the method, we would like to explain here the interest of our approach for numerics. If the original JKO scheme is a way to approximate the solutions of a parabolic PDEs via a sequence of variational problems, the fully discrete JKO scheme that we propose is an effective way to approximate it via a sequence of finite-dimensional variational problems. In the examples that we analyze in detail, these optimization problems are convex. The unknown in the problems can be the new measure $\rho_{k+1}^{\tau}$ or the transport plan $\gamma$ between it and the previous one. If the number of points in the space discretization is $N = O(h^{-d})$, then the problem is set on a space $\mathbb{R}^N$ (with positivity and mass constraints) in the first case, in a space $\mathbb{R}^{N \times N}$ in the second. On the other hand, the function to be minimized is more explicit when expressed in terms of the transport plan as it includes a linear part and a non-linear one depending on its projection onto the space $\mathbb{R}^N$. In the particular case of the crowd motion model, the whole problem is linear in $\gamma$ and only subject to linear equality or inequality constraints, hence it is a typical linear programming problem (LP). This LP problem (or the one appearing as a part of the problem when $f$ is nonlinear) can be solved by the simplex method, but this can be very slow, or better by the network simplex, which has been implemented in a very efficient and fast way in [10, 11]. However, since our main result states that one should choose $h/\tau$ very small (because of the reverse CFL condition), the very large number of points in the space discretization can be an obstacle in the practical implementation of the method.

Luckily, there exist at least two ways to decrease the complexity of the problem by reducing the number of unknowns of the problem. The first one is a clever idea from [4], which exploits the separability of the quadratic cost in terms of costs on each coordinate direction and transform this into an LP problem with $N^{1+1/d} = O(h^{-(d+1)})$ unknowns instead of $N^2 = O(h^{-2d})$. Note that this requires the use of a regular grid for the discretization, which we do anyway because it is also needed in some of our proofs. If some parts of the proofs could be adapted to other discretizations, some others are really based on the crucial role played by the coordinate axis. We also mention the fact that this requirement to choose a regular grid prevents us from easily considering local refinements of the mesh, which is instead a very common feature of many efficient numerical approaches to solve parabolic PDEs.

The second improvement is based on the idea that, when $\tau$ is small, the displacement between $\rho_k^\tau$ and $\rho_{k+1}^\tau$ will be small. If one can prove an $L^\infty$ bound on the displacement $|x-y|$ in the support of the optimal plan $\gamma$, then $\gamma$ will be concentrated on a neighborhood of the diagonal and some points in the product space can be removed from the unknown. If a bound of the form $|x-y| \le C\tau^\alpha$ can be proven (ideally with $\alpha = 1$), the number of unknowns obtained thanks to [4] goes down to $\tau^\alpha h^{-(d+1)}$. A very general method to prove this kind of $L^\infty$ bounds is contained in [12], and the proof therein can be easily adapted to the case of a grid. This unfortunately provides an exponent $\alpha < 1$ while, with very different methods (PDE-based and difficult to adapt to the case of a grid), at least in the case of linear diffusion, one can obtain $\alpha = 1$ (see [28]). Notice however that in the crowd motion case the bound with $\alpha = 1$ is false, as one can see from the sharp estimate obtained in [24] where the velocity is proven to belong to $L^4$ in space and time, and not better, while $|x-y| \le C\tau$ would imply an $L^\infty$ bound on the velocity.

We also mention the fact that the discrete transport problem that we need to solve could also be solved via the Sinkhorn algorithm based on entropic regularization, which would insert an extra parameter $\epsilon$ to be handled. The recent result in [5] explains how to use Sinkhorn without taking $\epsilon/\tau \to 0$ if one wants to approximate an equation already including linear diffusion. However, we do not consider in this paper the combination of the two methods, i.e. Sinkhorn on a grid with a proof of convergence.
All in all, the effective implementation of the approximation suggested in the present paper opens a certain number of challenging questions both in terms of estimates and of complexity. The convergence rate in terms of h and τ should also be considered for a rigorous comparison to other methods. In this comparison, this fully discrete JKO scheme is for sure a costly approach due to the smallness of the parameter h, but it has the advantage of not suffering from instability or slower convergence in the regions where the density is small, and it works in the same way for linear or nonlinear diffusion. It is, of course, a bad choice to approximate the linear Fokker-Planck equation, but it can be competitive for nonlinear equations.

**Outline.** The rest of the paper is organized as follows. In Section 2, we present the setting of the problem in full detail, discuss the preliminaries, and formulate the main results. Section 3 is dedicated to the toy example with pure potential energy, where we motivate the necessity of the assumption $h/\tau \to 0$. The convergence result for the fully discrete JKO scheme for the nonlinear drift-diffusion equation is proven in Section 4. Finally, we adapt the proof of the convergence of the scheme to the crowd motion case in Section 5. Note that the precise assumptions on the integrand f and on the potential V defining the energy F will only be given at the beginning of Section 4.

**Acknowledgments.** The authors acknowledge the support of the European Union via the ERC AdG 101054420 EYAWKAJKOS.

## 2. SETTING OF THE PROBLEM AND STATEMENT OF THE MAIN RESULTS

### 2.1. Continuous models.
In this section, we discuss two PDEs that we expect to obtain in the discrete-to-continuum limit: the nonlinear drift-diffusion equation and the crowd motion model. In both cases, we need the characterization of the weak solutions by the energy-dissipation inequality.

Let $\Omega \subset \mathbb{R}^d$ be a bounded domain, either convex or with a smooth boundary.

#### 2.1.1. Nonlinear diffusion equation.
Consider the nonlinear diffusion equation with drift

$$ (2.1) \qquad \partial_t \rho + \text{div}(\rho \nabla(f'(\rho) + V)) = 0, \quad \text{in } (0, T) \times \Omega $$

complemented with the non-flux boundary condition

$$ \partial_{\nu} L(\rho) + \rho \partial_{\nu} V = 0 \quad \text{on } \partial\Omega, $$

where the function L is defined by $L'(s) = sf''(s)$ and $\nu$ is the outer normal vector on $\partial\Omega$. The free energy functional for this equation is

$$ (2.2) \qquad \mathcal{F}(\rho) = \begin{cases} \int_{\Omega} f(u) \, d\mathcal{L}^d + \int_{\Omega} V \, d\rho & \text{if } \rho = u\mathcal{L}^d \\ +\infty & \text{otherwise.} \end{cases} $$

An essential tool for the convergence results is the characterization of solutions by the energy-dissipation balance. We first define the Fisher information as follows.

**Definition 2.1 (Fisher information).** Let the energy be defined as in (2.2). We defined the corresponding Fisher information as the functional $S: \mathcal{P}(\Omega) \to [0, +\infty]$:

$$ (2.3) \qquad S(\rho) = \begin{cases} \frac{1}{2} \int_{\Omega} |\nabla \ell(u) + \sqrt{u} \nabla V|^2 \, d\mathcal{L}^d & \text{if } \rho = u\mathcal{L}^d \text{ and } \ell(u) \in H^1(\Omega), \\ +\infty & \text{otherwise,} \end{cases} $$
where $\ell : (0, \infty) \to \mathbb{R}$ is a function defined by $\ell'(s) = \sqrt{s}f''(s)$ for $s \in (0, \infty)$. Note that this definition is algebraically equivalent to $\frac{1}{2}\int_{\Omega}|\nabla(f'(u) + V)|^2 d\rho$ (which also equals $\frac{1}{2}\int_{\Omega}|\nabla\frac{\delta F}{\delta\rho}|^2 d\rho$), but we prefer for the rest of our analysis to use quantities which are not weighted with $\rho$.

**Example 2.2.** As mentioned in the introduction, we use the Fokker-Planck and porous media equations as running examples. For the Fokker-Planck equation, a simple calculation gives $\ell'(s) = s^{-1/2}$ and $\ell(s) = 2\sqrt{s}$. Thus, the Fisher information for $\rho = u\mathcal{L}^d$ is

$$S(\rho) = \int_{\Omega} |2\nabla\sqrt{u} + \sqrt{u}\nabla V|^2 d\mathcal{L}^d$$

For the porous media equation, where $f(s) = s^m/(m-1)$ with $m > 1$ and $V = 0$. In this case, $\ell'(s) = ms^{m-3/2}$ and $\ell(s) = ms^{m-1/2}/(m-1/2)$. The Fisher information for $\rho = u\mathcal{L}^d$ becomes

$$S(\rho) = \frac{m^2}{(m-1/2)^2} \int_{\Omega} |\nabla u^{m-1/2}|^2 d\mathcal{L}^d.$$

**Proposition 2.3 (Characterization of solutions).** Let a pair $(\rho, v)$ satisfy the continuity equation

$$ (2.4) \qquad \partial_t \rho + \operatorname{div}(\rho v) = 0 \quad \text{in } (0, T) \times \Omega $$

in the distributional sense with the no-flux boundary condition. Further, let $(\rho, v)$ satisfy the EDI:

$$ (2.5) \qquad F(\rho_T) - F(\rho_0) + \int_0^T \left\{ \frac{1}{2} \int_{\Omega} |v_t|^2 d\rho_t + S(\rho_t) \right\} dt \le 0. $$

If the chain rule holds, i.e., for any curve satisfying (2.4)

$$ (2.6) \qquad \frac{d}{dt} F(\rho_t) = \int \nabla(f'(u_t) + V) \cdot v_t d\rho_t, $$

then $\rho \in C((0, T); \mathcal{P}(\Omega))$ is a distributional solution for (2.1) with $v_t = -\nabla f'(\rho_t) - \nabla V$.

The point of Proposition 2.3 is that the EDI characterizes the weak solution given that the chain rule (2.6) holds. It can be quite technical to prove that the chain rule holds for an arbitrary $f$. The proof of the chain rule can be based on [3, Theorem 10.4.8, Theorem 10.4.13], where $f$ is assumed to be a convex differentiable function with superlinear growth and to satisfy the McCann condition [41], i.e., the map $s \mapsto f(s^{-d})s^d$ is convex and non increasing in $(0, +\infty)$, which yields the geodesic convexity of the internal energy. For example, the Fokker-Planck and porous media equations satisfy these assumptions.

Instead of the McCann condition, the chain rule can be proved assuming that the initial datum is bounded from above and below by a positive constant [1]. A relaxation of the McCann condition appears in [15], where $f$ is decomposed, after some approximation, into the difference of two functions satisfying the McCann condition, and a mild assumption on integrability of the initial datum is required. In this paper, we do not go into detail about proving the chain rule. This indeed only concerns the limit equation, and is actually quite standard for porous-medium-type equations. Hence, we will simply use the EDI condition (2.5) to characterize the solution of (2.1).

2.1.2. Crowd motion model. Here we explain the macroscopic crowd motion model [39, 40] and its formulation as a gradient flow in $(\mathcal{P}(\Omega), W_2)$ formalized using the JKO scheme and the energy-dissipation balance.

Consider a crowd described by an absolutely continuous measure $\rho = u\mathcal{L}^d \in \mathcal{P}(\Omega)$. To account for the hard constraint on the concentration of people, we introduce the congestion constraint requiring the density $u \in L^1(\Omega)$ to remain below 1. We assume the population wishes to move with the velocity $w =$
$-\nabla V$ given by a potential $V \in C^1(\mathbb{R}^d)$. The gradient form of the velocity is the necessary assumption for the model to have the gradient structure. The spontaneous velocity $w$ is not always possible to realize due to the congestion constraint. Thus, the actual velocity is defined as a $L^2$ projection of $w$ on the cone of admissible velocities $C_\rho$. Intuitively, the cone $C_\rho$ is the set of velocities that do not increase the density $u$ in the saturated region $\{u=1\}$. In the PDE terms, the time evolution of $\rho$ is defined by the advection equation with the projected velocity:

$$ (2.7) \qquad \begin{aligned} \partial_t \rho + \operatorname{div}(\rho v) &= 0 \\ v &= \operatorname{Proj}_{C_\rho} w. \end{aligned} $$

The cone of admissible velocities is defined in duality with the pressure variable $p$, which plays the role of the Lagrangian multiplier for the constraint $u \le 1$. Thus, the set of density-dependent pressure variables is

$$ H_u^1 = \{p \in H^1(\Omega) : p \ge 0 \text{ and } p(1-u) = 0 \text{ a.e. in } \Omega\} $$

and then

$$ C_\rho := \{v \in L^2(\Omega; \mathbb{R}^d) : \int_{\Omega} \nabla p \cdot v \, d\rho \le 0 \quad \forall p \in H_u^1\}. $$

The model has a gradient structure in $(\mathcal{P}(\Omega), W_2)$ with the driving energy

$$ (2.8) \qquad \mathcal{F}_{\text{CM}}(\rho) = \begin{cases} \int_{\Omega} V u \, d\mathcal{L}^d & \text{if } \rho = u \mathcal{L}^d \text{ and } u \le 1, \\ +\infty & \text{otherwise.} \end{cases} $$

One can define this gradient flow with the JKO scheme. As before, given an initial datum $\rho_0 \in \text{spt}(\mathcal{F}_{\text{CM}})$ and a time step $\tau > 0$, we set up an iterative scheme

$$ (\text{JKO}_{\text{cm}}^{\tau}) \qquad \rho_{k+1}^{\tau} \in \arg \min_{\rho \in \mathcal{P}(\Omega)} \left\{ \mathcal{F}_{\text{CM}}(\rho) + \frac{1}{2\tau} W_2^2(\rho, \rho_k^{\tau}) \right\}. $$

Again, we use the energy-dissipation balance to characterize solutions. The Fisher information in this case is defined as follows.

**Definition 2.4.** For a density-pressure pair $(u, p)$, we define the Fisher information as

$$ (2.9) \qquad S_{\text{CM}}(u, p) := \begin{cases} \frac{1}{2} \int_{\Omega} |\nabla p + \nabla V|^2 u \, d\mathcal{L}^d & \text{if } p \in H_u^1 \\ +\infty & \text{otherwise.} \end{cases} $$

**Proposition 2.5** (Characterization of solutions for crowd motion). Let a pair $(\rho, v)$ satisfy the continuity equation (2.4) in the distributional sense with the no-flux boundary condition. Let the initial datum $\rho_0 = u_0 \mathcal{L}^d$ be such that $\mathcal{F}_{\text{CM}}(\rho_0) < \infty$ and $p \in L^\infty([0, T]; H^1(\Omega))$ be the pressure variable such that $p(1-u) = 0$. Further, let $(\rho, p, v)$ satisfy the EDI:

$$ (2.10) \qquad \mathcal{F}_{\text{CM}}(\rho_T) - \mathcal{F}_{\text{CM}}(\rho_0) + \int_0^T \left\{ \frac{1}{2} \int_{\Omega} |v_t|^2 \, d\rho_t + S_{\text{CM}}(u_t, p_t) \right\} \, dt \le 0. $$

then $\rho \in C((0, T); \mathcal{P}(\Omega))$ is a distributional solution for the crowd motion model with $v_t = -\nabla p_t - \nabla V$.

*Proof.* Using that $(\rho, v)$ satisfies the continuity equation, we directly get that

$$ \mathcal{F}_{\text{CM}}(\rho_t) - \mathcal{F}_{\text{CM}}(\rho_0) = \int_{\Omega} V \, d\rho_t - \int_{\Omega} V \, d\rho_s = \int_0^t \int_{\Omega} \nabla V \cdot v_s \, d\rho_s \, dt. $$

To prove the chain rule, we need to show that

$$ \mathcal{F}_{\text{CM}}(\rho_t) - \mathcal{F}_{\text{CM}}(\rho_0) = \int_0^t \int_{\Omega} (\nabla V + \nabla p_s) \cdot v_s \, d\rho_s \, dt, $$
which means that we need to prove $\int \nabla p \cdot v d\rho = 0$. We proceed similarly as in [25, Proposition 4.7]. Fix a time $t_0 > 0$ and an admissible pressure $q \in H_{u_{t_0}}^1$. Consider the function $(t_0 - \varepsilon, t_0 + \varepsilon) \ni t \mapsto \int_{\Omega} q d\rho_t$. We have

$$ \int_{\Omega} q \, d\rho_t \le \int_{\Omega} q \, d\mathcal{L}^d \quad t \in (t_0 - \varepsilon, t_0 + \varepsilon) $$

and the equality is achieved when $t = t_0$. If the map is differentiable at $t_0$, then

$$ \left. \frac{d}{dt} \int_{\Omega} q \, d\rho_t \right|_{t=t_0} = 0, $$

and it implies

$$ \int_{0}^{t} \int_{\Omega} \nabla p_s \cdot v_s \, d\rho_s = 0. $$

The subtlety with the differentiability is that for any $q \in H^1(\Omega)$ we have for almost every $t$

$$ \frac{d}{dt} \int_{\Omega} q \, d\rho_t = \int_{\Omega} \nabla q \cdot v_t \, d\rho_t, $$

but it can happen that the derivative does not exist at $t_0$ for any $q$. Fortunately, we can show that the set of times, for which the derivative on the left-hand side exists for every $q \in H^1(\Omega)$ has full measure. By the continuity equation, this set of times can be chosen as all $t_0 \in [0, T]$ such that

$$ \int_{t_0}^{t_0+\varepsilon} v_t u_t \, dt \to v_{t_0} u_{t_0} \quad \text{weakly in } L^2(\Omega). $$

The latter convergence holds by the density argument and the bounds $\int_0^T \int_{\Omega} |v_t|^2 d\rho_t dt < \infty$ and $\|u_t\|_{\infty} \le 1$. $\square$

**2.2. Discretization.** We define a discretization of $\Omega$ as $\mathcal{T}^h = h\mathbb{Z}^d \cap \Omega$ for any $h > 0$. An arbitrary probability measure on $\mathcal{T}^h$ has the form

$$ \rho^h = \sum_{z \in \mathcal{T}^h} \rho_z^h \delta_z, \quad \sum_{z \in \mathcal{T}^h} \rho_z^h = 1. $$

It will be convenient for some calculations to use a piecewise constant counterpart of $\rho^h$. Instead of assigning the mass $\rho_z^h$ to one point $z \in \mathcal{T}^h$, we assign $\rho_z^h$ to the cell $Q_h(z) := z + (-h/2, h/2)^d$ around $z$ with the uniform density. Precisely, we define

$$ (2.11) \qquad \hat{\rho}^h := \sum_{z \in \mathcal{T}^h} \rho_z^h \mathcal{L}^d |_{Q_h(z)}. $$

The measure $\hat{\rho}^h$ is absolutely continuous with respect to the Lebesque measure and has the density:

$$ \hat{u}^h := \frac{d\hat{\rho}^h}{d\mathcal{L}^d} = \sum_{z \in \mathcal{T}^h} \frac{\rho_z^h}{h^d} 1_{Q_h(z)} = \sum_{z \in \mathcal{T}^h} u_z^h 1_{Q_h(z)}. $$

In the same way, we define the piecewise constant reconstruction for any bounded discrete function $f^h \in B(\mathcal{T}^h)$:

$$ (2.12) \qquad \hat{f}^h := \sum_{z \in \mathcal{T}^h} f_z^h 1_{Q_h(z)}. $$

For any $z \in \mathcal{T}^h$ (away from a boundary), the set of neighbors has $2d$ points of the form $z + h$ with $h \in \{he_1, -he_1, \dots, he_d, -he_d\} =: d_h$. We denote by $\Sigma^h \subset \mathcal{T}^h \times \mathcal{T}^h$ the set of pairs of neighboring
points, i.e., $\Sigma^h := \{(z, \zeta) : |z - \zeta| = h\}$. The interface between two neighboring cells $Q_h(z)$ and $Q_h(\zeta)$, $(z, \zeta) \in \Sigma^h$, is $I_h(z, \zeta) := \overline{Q_h(z)} \cap \overline{Q_h(\zeta)}$.

We define the fully discrete JKO scheme as a counterpart for (JKO) with the minimization restricted to the atomic measures on the regular lattice $\mathcal{T}^h$. We use the discrete energy functional defined as

$$ (2.13) \qquad \mathcal{F}^h(\rho^h) = \sum_{z \in \mathcal{T}^h} f\left(\frac{\rho_z^h}{h^d}\right)h^d + \sum_{z \in \mathcal{T}^h} V(z)\rho_z^h. $$

Note that the discrete version of the internal energy can be obtained by inserting the piecewise constant reconstruction $\hat{\rho}^h$ into the continuous internal energy

$$ \mathcal{F}(\hat{\rho}^h) = \int_{\Omega} f\left(\frac{d\hat{\rho}^h}{d\mathcal{L}^d}\right) d\mathcal{L}^d = \sum_{z \in \mathcal{T}^h} f\left(\frac{\rho_z^h}{h^d}\right)h^d = \mathcal{F}^h(\rho^h). $$

Instead, for the potential energy, we can use the atomic measure to obtain the discrete energy. The discrete energy for the crowd motion model takes the form

$$ (2.14) \qquad \mathcal{F}_{\text{CM}}(\rho^h) := \begin{cases} \sum_{z \in \mathcal{T}^h} V(z)\rho_z^h & \text{if } u^h = \frac{\rho^h}{h^d} \le 1, \\ +\infty & \text{otherwise,} \end{cases} $$

which is a natural counterpart for (2.8).

**Definition 2.6.** For a given $h > 0$ and $\tau > 0$ of the form $\tau = T/N$ for $T > 0$, $N \in \mathbb{N}$, and given initial data $\rho_0^h \in \mathcal{P}(\mathcal{T}^h)$, we define $\{\rho_k^{h,\tau}\}_{k=1,...,N}$ as

$$ (\text{JKO}^{h,\tau}) \qquad \rho_{k+1}^{h,\tau} \in \arg \min_{\rho^h \in \mathcal{P}(\mathcal{T}^h)} \left\{ \mathcal{F}^h(\rho^h) + \frac{1}{2\tau} W_2^2(\rho^h, \rho_k^h) \right\}. $$

It will be convenient to use the notation

$$ \mathcal{J}_{\rho_0^h}^{h,\tau}(\rho^h) := \mathcal{F}^h(\rho^h) + \frac{1}{2\tau} W_2^2(\rho^h, \rho_0^h). $$

**Remark 2.7 (Propagation of the $L^\infty$ bound).** Throughout the paper, we do not require additional assumptions on the initial data. Yet we add the remarks indicating where the proofs can be simplified if we assume a uniform $L^\infty$ bound on $\{\rho_0^h\}_{h>0}$. In those remarks, we rely on the propagation of the $L^\infty$ bound through iterations of the fully discrete JKO scheme. In particular, if we assume an $L^\infty$ bound on the initial data, then we do not need (f1) for the nonlinear diffusion equation without drift. Here we provide an argument for the propagation of the $L^\infty$ bound.

Consider $\mathcal{F}^h$ given in (2.13) with $V = 0$. Let $M := \max_{z \in \mathcal{T}^h} \rho_0^h(z)$ and suppose that $S := \{z \in \mathcal{T}^h : \rho_1^{h,\tau}(z) > M\}$ is not empty. Then

$$ \sum_{z \in S} \rho_1^{h,\tau}(z) > |S|M \ge \sum_{z \in S} \rho_0^h(z). $$

Therefore, the mass has to be transported from outside $S$ to inside $S$. More precisely, there exists $(x, y) \in \text{spt}(\gamma)$, with $\gamma$ being the optimal coupling between $\rho_1^{h,\tau}$ and $\rho_0^h$, such that $x \in S$, $y \notin S$, and $\gamma_{xy} > 0$.

We construct a competitor as $\tilde{\rho}_1^{h,\tau} := \rho_1^{h,\tau} - \gamma_{xy}\delta_x + \gamma_{xy}\delta_y$, together with a new transport plan $\tilde{\gamma} \in \mathcal{P}(\mathcal{T}^h \times \mathcal{T}^h)$ defined as $\tilde{\gamma} := \gamma + \gamma_{xy}\delta_{(y,y)} - \gamma_{xy}\delta_{(x,y)}$. This plan is a coupling between $\tilde{\rho}_1^{h,\tau}$ and $\rho_0^h$,
and using it we find

$$W_2^2(\rho_0^h, \tilde{\rho}_1^{h,\tau}) \le W_2^2(\rho_0^h, \rho_1^{h,\tau}).$$

Furthermore, by convexity of $f$, we get

$$ \begin{aligned} \mathcal{F}^h(\tilde{\rho}_1^{h,\tau}) &= \mathcal{F}^h(\rho_1^{h,\tau}) + f\left(\frac{\tilde{\rho}_1^{h,\tau}(x)}{h^d}\right)h^d - f\left(\frac{\rho_1^{h,\tau}(x)}{h^d}\right)h^d + f\left(\frac{\tilde{\rho}_1^{h,\tau}(y)}{h^d}\right)h^d - f\left(\frac{\rho_1^{h,\tau}(y)}{h^d}\right)h^d \\ &\le \mathcal{F}^h(\rho_1^{h,\tau}) - f'\left(\frac{\rho_1^{h,\tau}(x)}{h^d}\right)\gamma_{xy} + f'\left(\frac{\rho_1^{h,\tau}(y)}{h^d}\right)\gamma_{xy}, \end{aligned} $$

where we recall that $\rho_1^{h,\tau}(x) > M$ and $\rho_1^{h,\tau}(y) \le M$, thus, $\mathcal{F}^h(\tilde{\rho}_1^{h,\tau}) < \mathcal{F}^h(\rho_1^{h,\tau})$. Hence $\mathcal{J}_{\rho_0^h}^{h,\tau}(\tilde{\rho}_1^{h,\tau}) < \mathcal{J}_{\rho_0^h}^{h,\tau}(\rho_1^{h,\tau})$ we obtain a contradiction with $\rho_1^{h,\tau}$ being a minimizer.

2.3. **Main results.** Before we get to the main results, we study a toy model in Section 3. We consider evolution driven purely by a potential $V$, i.e., the energy functional is

$$\mathcal{F}_V(\rho) = \int_{\Omega} V \, d\rho.$$

This toy model is instructive because it illustrates the necessary relation between the time and space discretization parameters $\tau$ and $h$ for passing to the joint limit. Already for the toy model as well as the main results we need $h/\tau \to 0$.

The first main result is the convergence of $(\text{JKO}^{h,\tau})$ for the nonlinear diffusion. We require different assumptions on the energy density $f$ depending on the presence or absence of the external potential $V$. If $V = 0$, then we assume the following.

**Assumptions on the energy density.** The energy density $f : [0, +\infty) \to \mathbb{R}$ satisfies the following:

(f0) $f$ is convex and $C^2((0, +\infty))$

(f1) $f$ has superlinear growth at infinity.

The superlinearity assumption is quite standard in this setting. Removing it would require to re-define the functional $\mathcal{F}$ using the recession function $f^\infty$ (see, for instance, Chapter 7 in [44]). Moreover, this assumption guarantees equiintegrability of the densities provided by the scheme, which is crucial in some arguments throughout the paper. This requirement has been weakened in related papers, such as in [5] (but the arguments there are quite technical) or in [15] (but in such a paper equiintegrability is obtained by propagating in time the integrability properties of the initial density; the method to prove this strongly relied on displacement convexity and cannot be reproduced in our setting).

If the potential $V$ is present, then, in addition to (f0) and (f1), we need $f$ to satisfy a sort of quantifiable convexity assumption, which is made precise below.

**Assumptions on the energy density.** The energy density $f : [0, +\infty) \to \mathbb{R}$ satisfies (f0), (f1), and the following:

(i) $f'$ is strictly increasing;

(ii) there exist $C_f > 0$, $\theta < 1 + 1/d$, and $s_0 \ge 1$ such that for any $s \ge s_0$ : $f''(s) \ge C_f s^{-\theta}$

(iii) $V \in \text{Lip}(\mathbb{R}^d)$
**Remark 2.8.** For the energy densities $f(s) = s^m/(m-1)$, $m > 1$ corresponding to the porous media equation, (fV) follows directly from (f1) with $\theta = 1$. Yet, in general, there are examples of $f$ such as $f(s) = s^2/2 + \sin(s)$, which is convex and superlinear at infinity, but does not satisfy (fV).

**Theorem 2.9 (Main result for nonlinear diffusion).** Assume (f0), (f1), and either $V = 0$ or (fV). Let $\{\rho_k^{h,\tau}\}_{k=1,...,N}$ be the family of $(JKO^{h,\tau})$ minimizers as in Definition 2.6 with $\mathcal{F}^h$ defined in (2.13). Let the family of initial data $\{\rho_0^h\}_{h>0}$ be such that there exists $\rho_0 \in \mathcal{P}(\Omega)$ with finite energy $\mathcal{F}(\rho_0) < \infty$ such that

$$ \rho_0^h \to \rho_0 \quad \text{narrowly as } h \to 0 \quad \text{and} \quad \lim_{h \to 0} \mathcal{F}^h(\rho_0^h) = \mathcal{F}(\rho_0). $$

Then a suitable interpolation in time of $\{\rho_k^{h,\tau}\}_{k=1,...,N}$ converges as $h, \tau, h/\tau \to 0$ uniformly in time for the $W_2$ distance to an absolutely continuous curve $[0,T] \ni t \mapsto \rho_t \in \mathcal{P}(\Omega)$ such that $\rho_t$ is the distributional solution of (2.1).

**Remark 2.10.** Another common example of driving energy is interaction energy that takes the form

$$ \mathcal{F}_W(\rho) = \int_{\Omega} (W * \rho)(x)\rho(dx) $$

with the natural discrete counterpart

$$ \mathcal{F}_W^h(\rho^h) = \sum_{x,y \in \mathcal{T}^h \times \mathcal{T}^h} W(x-y)\rho_x^h \rho_y^h. $$

Assuming that the interaction potential $W: \mathbb{R}^d \to \mathbb{R}$ is regular enough, i.e. symmetric, non-negative, and $W \in \text{Lip}(\mathbb{R}^d) \cap C^1(\mathbb{R}^d \setminus \{0\})$, the proof of Theorem 2.9 can be adapted to deal with the sum of the internal and interaction energy instead of the internal and potential energy. This means that $(JKO^{h,\tau})$ approximates the solutions of the aggregation-diffusion equation.

We choose not to focus on the interaction energy in this paper, because the optimization problem in $(JKO^{h,\tau})$ becomes non-convex, which is less attractive from the numerical point of view presented in Section 1.4. Nevertheless, the adaptation of the proofs for $W * \rho$ instead of $V$ seems not to be too difficult, and we comment on them in Remarks 4.10 and 4.19.

For the crowd motion model, we need a smoother potential.

**Assumptions on the potential for the crowd motion model.**

$$ (\text{pV}) \qquad V \in C^1(\mathbb{R}^d) $$

The main convergence result in this case is as follows.

**Theorem 2.11 (Main result for crowd motion).** Assume (pV). Let $\{\rho_k^{h,\tau}\}_{k=1,...,N}$ be the family of $(JKO^{h,\tau})$ minimizers as in Definition 2.6 with $\mathcal{F}^h$ defined in (2.14). Let the family of initial data $\{\rho_0^h\}_{h>0}$ be such that there exists $\rho_0 = u_0 \mathcal{L}^d \in \mathcal{P}(\Omega)$ with $\|u_0\|_{L^\infty} \le 1$ such that

$$ \rho_0^h \to \rho_0 \quad \text{narrowly as } h \to 0 \quad \text{and} \quad \lim_{h \to 0} \mathcal{F}_{CM}(\rho_0^h) = \mathcal{F}_{CM}(\rho_0). $$

Then a suitable interpolation in time of $\{\rho_k^{h,\tau}\}_{k=1,...,N}$ converges as $h, \tau, h/\tau \to 0$ uniformly in time for the $W_2$ distance to an absolutely continuous curve $[0,T] \ni t \mapsto \rho_t \in \mathcal{P}(\Omega)$ such that $\rho_t$ is the distributional solution of (2.7).
### 3. EVOLUTION DRIVEN BY A POTENTIAL

In this section, we consider the fully discrete JKO scheme $(\text{JKO}^{h,\tau})$ with the potential energy

$$ \mathcal{F}^h(\rho^h) = \sum_{x \in \mathcal{T}^h} V(x) \rho_x^h. $$

We study this case as a toy example to understand the relation between discretization-in-space parameter $h$ and time step parameter $\tau$ that are needed to pass to the joint limit $h, \tau \to 0$. We see in Section 3.1 that the condition $h/\tau \to 0$ is a necessary condition for the convergence of $(\text{JKO}^{h,\tau})$ to hold. In Section 3.2, we show the convergence of $(\text{JKO}^{h,\tau})$ for the potential energy with $V \in C^{1,1}(\mathbb{R}^d)$.

**3.1. Relation between $h$ and $\tau$.** We observe that the one step of $(\text{JKO}^{h,\tau})$ can be rewritten as a linear programming problem:

$$ \begin{aligned} \min_{\rho^h} \mathcal{J}_{\rho_0^h}^{h,\tau}(\rho^h) &= \min_{\rho^h} \left\{ \sum_{x \in \mathcal{T}^h} V(x) \rho_x^h + \frac{1}{2\tau} \min_{\gamma \in \Gamma(\rho_0^h, \rho^h)} \left\{ \sum_{x,y \in \mathcal{T}^h \times \mathcal{T}^h} |x-y|^2 \gamma_{xy}^h \right\} \right\} \\ &= \min_{\gamma^h, P^1, P^2, \gamma^h = \rho_0^h} \left\{ \sum_{x,y \in \mathcal{T}^h \times \mathcal{T}^h} \left(V(x) + \frac{1}{2\tau} |x-y|^2\right) \gamma_{xy}^h \right\}. \end{aligned} $$

With this simple energy, one can consider the movement of each Dirac mass $\rho_0^h(x)\delta_{x(t)}$ independently of the others. Thus, to find the minimizer of $\mathcal{J}_{\rho_0^h}^{h,\tau}$ it is enough to solve for all $x_0 \in \text{spt}(\rho_0^h)$:

$$ (\text{Min-}\mathcal{T}_h) \qquad \min_{x \in \mathcal{T}^h} \left\{ V(x) + \frac{1}{2\tau} |x - x_0|^2 \right\}. $$

Here, we can make the first observation about the relation between $h$ and $\tau$. If $x_1$ is the minimizer of $(\text{Min-}\mathcal{T}_h)$, then

$$ V(x_1) + \frac{1}{2\tau} |x_1 - x_0|^2 \le V(x_0). $$

A simple rewriting gives

$$ \frac{1}{2\tau} |x_1 - x_0| \le \frac{V(x_0) - V(x_1)}{|x_1 - x_0|} \le \text{Lip}(V). $$

Since the distance between two different points in $\mathcal{T}^h$ is at least $h$, then a necessary condition for $x_1$ to be different from $x_0$ is

$$ (3.1) \qquad \frac{h}{\tau} \le 2\text{Lip}(V). $$

Consequently, if (3.1) does not hold, then the evolution is frozen.

To get error estimates, we compare the minimization problem $(\text{Min-}\mathcal{T}_h)$ constrained to the discrete set $\mathcal{T}^h$ with its counterpart without the constraint:

$$ (\text{Min-}\mathbb{R}^d) \qquad \min_{x \in \mathbb{R}^d} \left\{ V(x) + \frac{1}{2\tau} |x - x_0|^2 \right\}. $$

We begin with even simpler example of the linear potential.

**Example 3.1 (Linear potential).** Let $V(x) = -v \cdot x$ with some constant velocity $v \in \mathbb{R}^d$. The minimizer of $(\text{Min-}\mathbb{R}^d)$ is clearly $x_1 = x_0 + \tau v$. Let $x_1^h$ be the minimizer of $(\text{Min-}\mathcal{T}_h)$. Note that

$$ -v \cdot x + \frac{1}{2\tau} |x - x_0|^2 = \frac{1}{2\tau} \left( |x - (x_0 + \tau v)|^2 - |x_0 + \tau v|^2 + x_0^2 \right) = \frac{1}{2\tau} \left( |x - (x_0 + \tau v)|^2 - v \cdot x_0 - \frac{\tau}{2} |v|^2 \right), $$
therefore

$$x_1^h \in \arg \min_{x \in \mathcal{T}^h} \left\{ -v \cdot x + \frac{1}{2\tau} |x - x_0|^2 \right\} = \arg \min_{x \in \mathcal{T}^h} \left\{ |x - (x_0 + \tau v)|^2 \right\}.$$

This implies that $x_1^h = \text{Proj}_{\mathcal{T}^h}(x_1)$ and $|x_1 - x_1^h| \le \sqrt{dh}/2$ for all $x \in \mathbb{R}^d$. Furthermore, since the velocity does not depend on $x$, exactly the same error occurs at every iteration

$$|x_k - x_k^h| = k|x_1 - x_1^h| \quad k \in \mathbb{N}.$$

We then obtain $|x_1 - x_1^h| \sim h$ (except for some particular combinations of $v$ and $h$). Thus, for $t = k\tau$, we have

$$|x_k - x_k^h| = t \frac{|x_1 - x_1^h|}{\tau} \approx t \frac{h}{\tau}.$$

Therefore, $h/\tau \to 0$ is a necessary condition for the error to vanish as $h, \tau \to 0$.

Example 3.1 shows that the assumption $h/\tau \to 0$ is necessary for the convergence of $(\text{JKO}^{h,\tau})$ to the expected limiting minimizer. We assume $h/\tau \to 0$ in all the subsequent convergence results.

**3.2. Convergence of the fully discrete JKO with potential.** Let $V$ be an arbitrary $C^{1,1}(\mathbb{R}^d)$ potential. Denote by $x_1$ the minimizer of $(\text{Min-}\mathbb{R}^d)$. One can rewrite $(\text{Min-}\mathcal{T}_h)$ in the following way using the optimality condition $x_1 = x_0 - \tau \nabla V(x_1)$:

$$\begin{aligned} V(x) + \frac{1}{2\tau}(x - x_0)^2 &= \frac{1}{2\tau}|x - x_1|^2 + V(x) + \frac{1}{2\tau}|x - x_0|^2 - \frac{1}{2\tau}|x - x_0 + \tau \nabla V(x_1)|^2 \\ &= \frac{1}{2\tau}|x - x_1|^2 + V(x) - x \cdot \nabla V(x_1) + (\text{terms independent of } x) \\ &= \frac{1}{2\tau}|x - x_1|^2 + V(x) - V(x_1) - (x - x_1) \cdot \nabla V(x_1) + (\text{independent of } x). \end{aligned}$$

Therefore, $(\text{Min-}\mathcal{T}_h)$ is equivalent to

$$\min_{x \in \mathcal{T}^h} \left\{ \frac{1}{2\tau} |x - x_1|^2 + V(x) - V(x_1) - (x - x_1) \cdot \nabla V(x_1) \right\}.$$

We see the later formulation to have two parts: the distance to the minimizer of $(\text{Min-}\mathbb{R}^d) x_1$ and the part depending on the potential $V$. Since $V \in C^{1,1}(\mathbb{R}^d)$, it holds that

$$\begin{aligned} |V(x) - V(x_1) - (x - x_1) \cdot \nabla V(x_1)| &= \left| \int_0^1 \left( \nabla V(x_1 + \lambda(x - x_1)) - \nabla V(x_1) \right) d\lambda \cdot (x - x_1) \right| \\ &\le \frac{|\nabla V|_{\text{Lip}}}{2} |x - x_1|^2, \end{aligned}$$

where $|\cdot|_{\text{Lip}}$ is the Lipshitz constant of $\nabla V$.

We use two proximal operators defined as

$$\text{Prox}_{V}^{h,\tau}(z) := \arg \min_{x \in \mathcal{T}^h} \left\{ V(x) + \frac{|x-z|^2}{2\tau} \right\} \quad \text{and} \quad \text{Prox}_{V}^{\tau}(z) := \arg \min_{x \in \mathbb{R}^d} \left\{ V(x) + \frac{|x-z|^2}{2\tau} \right\},$$

and consider the corresponding sequences of minimizers. Given an arbitrary $x_0 \in \mathbb{R}^d$, we denote by $\{x_k^h\}_{k=1,...,N} \subset \mathcal{T}^h$ the sequence of the minimizers defined as $x_k^h := \text{Prox}_{V}^{h,\tau}(x_{k-1}^h)$ and by $\{x_k\}_{k=1,...,N} \subset \mathbb{R}^d$ the sequence of the minimizers defined as $x_k := \text{Prox}_{V}^{\tau}(x_{k-1})$.

**Lemma 3.2.** Consider the time interval $[0, T]$ with $T = N\tau$ for some $N \in \mathbb{N}$. Then the error estimate between the two sequences reads as

$$|x_N^h - x_N| \le C(T, |\nabla V|_{\text{Lip}}) \frac{h}{\tau}.$$
*Proof.* The difference between $x_{k+1}^h$ and $x_{k+1}$ can be decomposed into two parts:

$$|x_{k+1}^h - x_{k+1}| \le |\text{Prox}_V^{h,\tau}(x_k^h) - \text{Prox}_V^\tau(x_k^h)| + |\text{Prox}_V^\tau(x_k^h) - \text{Prox}_V^\tau(x_k)|.$$

The first term is the local error between minimization over $\mathbb{R}^d$ and minimizing over $\mathcal{T}^h$ with the same initial point $x_k^h$. The second term is the difference of minimizing over $\mathbb{R}^d$ for two different initial points $x_k^h$ and $x_k$.

To estimate the first type of error $|\text{Prox}_V^{h,\tau}(x_0) - \text{Prox}_V^\tau(x_0)|$ for an arbitrary $x_0 \in \mathbb{R}^d$, we use as a competitor a projection of $x_1$ on $\mathcal{T}^h$ denoted by $\tilde{x}_1$ ($|x_1 - \tilde{x}_1| \le h$):

$$\frac{1}{2\tau}(x_1^h - x_1)^2 + V(x_1^h) - V(x_1) - (x_1^h - x_1) \cdot \nabla V(x_1) \\ \le \frac{1}{2\tau}(\tilde{x}_1 - x_1)^2 + V(\tilde{x}_1) - V(x_1) - (\tilde{x}_1 - x_1) \cdot \nabla V(x_1).$$

Rearranging the terms gives

$$(x_1^h - x_1)^2 \le h^2 + 2\tau(V(\tilde{x}_1) - V(x_1^h) - (\tilde{x}_1 - x_1^h) \cdot \nabla V(x_1)) \\ = h^2 + 2\tau \int_0^1 (\nabla V(x_1^h + \lambda(\tilde{x}_1 - x_1^h)) - \nabla V(x_1)) \, d\lambda \cdot (\tilde{x}_1 - x_1^h) \\ \le h^2 + 2\tau \int_0^1 |\nabla V|_{\text{Lip}} |(\lambda - 1)(x_1 - x_1^h) + \lambda(\tilde{x}_1 - x_1)| \, d\lambda |\tilde{x}_1 - x_1^h| \\ \le h^2 + \tau |\nabla V|_{\text{Lip}} (|x_1 - x_1^h| + h)^2$$

Expanding the square, we obtain

$$(1 - |\nabla V|_{\text{Lip}} \tau) |x_1^h - x_1|^2 - 2 |\nabla V|_{\text{Lip}} \tau h |x_1^h - x_1| \le (1 + |\nabla V|_{\text{Lip}} \tau) h^2$$

Let $\tau > 0$ be such that $\tau \le (2|\nabla V|_{\text{Lip}})^{-1}$, then

$$\left( |x_1^h - x_1| - \frac{|\nabla V|_{\text{Lip}} \tau h}{1 - |\nabla V|_{\text{Lip}} \tau} \right)^2 \le \frac{h^2}{(1 - |\nabla V|_{\text{Lip}} \tau)^2},$$

and we get

$$|x_1^h - x_1| \le \frac{1 + |\nabla V|_{\text{Lip}} \tau}{1 - |\nabla V|_{\text{Lip}} \tau} h \le (1 + 4|\nabla V|_{\text{Lip}} \tau)h,$$

where we use the simple inequality $(1+x)/(1-x) \le 1+4x$ that holds for any $x \in [0, 1/2]$. Thus, for small enough $\tau$, the error introduced by the grid is bounded as

$$|\text{Prox}_V^{h,\tau}(x_0) - \text{Prox}_V^\tau(x_0)| \le (1 + 4|\nabla V|_{\text{Lip}} \tau)h := e_{h,\tau} \quad \forall x_0 \in \mathbb{R}^d.$$

To estimate the second part of the error, precisely, $|\text{Prox}_V^\tau(x_1^h) - \text{Prox}_V^\tau(x_1)|$, we write the optimality conditions:

$$\nabla V(x_2) + \frac{x_2 - x_1}{\tau} = 0 \quad \text{for } x_2 = \text{Prox}_V^\tau(x_1)$$

$$\nabla V(y_2) + \frac{y_2 - x_1^h}{\tau} = 0 \quad \text{for } y_2 = \text{Prox}_V^\tau(x_1^h).$$
Subtracting the second equation from the first one and multiplying the difference with $x_2^h - y_2$, we get

$$ (x_2 - y_2) \cdot (\nabla V(x_2) - \nabla V(y_2)) + (x_2 - y_2) \cdot \frac{x_2 - y_2 + x_1^h - x_1}{\tau} = 0. $$

We rearrange the terms as

$$ (1 - |\nabla V|_{\text{Lip}} \tau) |x_2 - y_2| \le |x_1 - x_1^h|. $$

A similar estimate can be obtained for any time step $k + 1$:

$$ (1 - |\nabla V|_{\text{Lip}} \tau) |x_{k+1} - y_{k+1}| \le |x_k - x_k^h|, $$

where $y_{k+1} = \text{Prox}_V^\tau(x_k^h)$. Then the total error at time step $k + 1$ is bounded as

$$ \begin{aligned} |x_{k+1} - x_{k+1}^h| & \le |x_{k+1} - y_{k+1}| + |y_{k+1} - x_{k+1}^h| \le \frac{1}{1 - |\nabla V|_{\text{Lip}} \tau} |x_k - x_k^h| + e_{h,\tau} \\ & \le e_{h,\tau} \sum_{j=0}^k \frac{1}{(1 - |\nabla V|_{\text{Lip}} \tau)^j} \le e_{h,\tau} \sum_{j=0}^k (1 + 2|\nabla V|_{\text{Lip}} \tau)^j \le e_{h,\tau} \frac{(1 + 2|\nabla V|_{\text{Lip}} \tau)^{k+1} - 1}{2|\nabla V|_{\text{Lip}} \tau} \\ & \le e_{h,\tau} \frac{1}{|\nabla V|_{\text{Lip}} \tau} e^{2|\nabla V|_{\text{Lip}} T}, \end{aligned} $$

where we use that $T = k\tau$. Recalling the estimate for $e_{h,\tau}$, we finally obtain

$$ |x_{k+1} - x_{k+1}^h| \le \frac{3}{|\nabla V|_{\text{Lip}}} e^{2|\nabla V|_{\text{Lip}} T} \frac{h}{\tau} $$

for $\tau \le (2|\nabla V|_{\text{Lip}})^{-1}$.

<figure><img src="image_2.png" /></figure>

## 4. NONLINEAR DIFFUSION

In this section, we consider the fully discrete JKO scheme (JKO$^{h,\tau}$) for the energy functional of the form

$$ (4.1) \qquad F^h(\rho^h) = \sum_{x \in \mathcal{T}^h} f\left(\frac{\rho_x^h}{h^d}\right) h^d + \sum_{x \in \mathcal{T}^h} V(x) \rho_x^h $$

with $f$ satisfying (f0), (f1), and either $V = 0$ or (fV). The goal is to prove the convergence of (JKO$^{h,\tau}$) stated in Theorem 2.9.

To explain the strategy that we will use to prove Theorem 2.9, we first recall the idea of the convergence proof for the standard JKO scheme and then indicate the adaptations required for the fully discrete case. The approach we have in mind relies on the variational (De Giorgi) interpolation, explained, for instance in [3, Chapter 3]. Given the minimizers $\{\rho_k^\tau\}$ of (JKO), the idea is to prove the inequality

$$ (4.2) \qquad F(\rho_{k+1}^\tau) - F(\rho_k^\tau) + \frac{W_2^2(\rho_k^\tau, \rho_{k+1}^\tau)}{2\tau} + \int_0^1 \frac{W_2^2(\rho_k^\tau, \rho_r^\tau)}{2\tau r^2} dr \le 0, $$

where $\rho_r^\tau$ is the variational interpolant between $\rho_k^\tau$ and $\rho_{k+1}^\tau$. Combining (4.2) with the lower bound on the Wasserstein distance by the slope of the energy [3, Lemma 3.1.1]

$$ (4.3) \qquad \frac{1}{2\tau^2 r^2} W_2^2(\rho_k^\tau, \rho_r^\tau) \ge |\partial F'|^2(\rho_r^\tau), $$
one gets

$$ (4.4) \qquad F(\rho_{k+1}^{\tau}) - F(\rho_{k}^{\tau}) + \frac{W_{2}^{2}(\rho_{k}^{\tau}, \rho_{k+1}^{\tau})}{2\tau} + \frac{\tau}{2} \int_{0}^{1} |\partial F'|^{2}(\rho_{r}^{\tau}) dr \le 0. $$

Summing up (4.4) for all discrete time steps gives

$$ (4.5) \qquad F(\rho_N^\tau) - F(\rho_0^\tau) + \sum_{k=0}^{N-1} \tau \frac{W_2^2(\rho_k^\tau, \rho_{k+1}^\tau)}{2\tau^2} + \frac{1}{2} \int_0^T \sum_{k=0}^{N-1} |\partial F'|^2(\rho_{r-k}^\tau) dr \le 0, $$

which is a discrete analog of EDI characterizing the solution of (2.1) (see Proposition 2.3 for the characterization by EDI). It remains to use an appropriate interpolation in time to obtain a sharp inequality which is convenient to pass to the limit to recover (2.5) (for more detail, we refer to [3]).

A challenging step in the fully discrete settings is to obtain a suitable replacement for (4.3). The definition of the local slope used in [3, Lemma 3.1.1] is unsuitable for passing to the discrete-to-continuum limit. Furthermore, we know from the consideration in Section 3 that the relation between the discretization parameters $h$ and $\tau$ plays a significant role in the convergence. Therefore, we expect the discrete counterpart for (4.3) to depend on the ratio $h/\tau$ explicitly.

We realize the following strategy to prove Theorem 2.9:

(1) We begin Section 4.1 by introducing the variational interpolation for the fully discrete JKO and obtaining an inequality resembling (4.2).

(2) We establish the crucial lower bound on the Wasserstein distance by the discrete Fisher information in Lemma 4.4 (for nonlinear diffusion without an external potential) and in Corollary 4.7 (with the potential). These inequalities provide the replacement for (4.3).

(3) We conclude Section 4.1 by showing the discrete-in-time-and-space version of EDI analogous to (4.5) (Lemma 4.12). This formulation of EDI provides the base to pass to the limit $h, \tau \to 0$.

(4) Section 4.2 is devoted to the required compactness results and the convergence result is wrapped up in Section 4.3.

### 4.1. Interpolations and corresponding inequalities.

**Definition 4.1 (Variational interpolation).** The variational interpolation $\tilde{\rho}^{h,\tau}$ for a sequence $\{\rho_k^{h,\tau}\}_{k \in \mathbb{N}}$ is defined for $t = (k+s)\tau$ with any $k \in 1, \dots, N$ and any $s \in (0, 1]$ as

$$ (4.6) \qquad \tilde{\rho}_{(k+s)\tau}^{h,\tau} = \arg \min_{\rho^h \in \mathcal{P}(\tau^h)} \left\{ F^h(\rho^h) + \frac{1}{2\tau s} W_2^2(\rho^h, \rho_k^h) \right\} $$

The standard way to get the inequality of the type (4.2) using the variational interpolation translates into the discrete setting without changes. We briefly present here the argument for completeness. Consider a function $G_k : [0, 1] \to \mathbb{R}$ defined as

$$ G_k(s) := \min_{\rho^h \in \mathcal{P}(\tau^h)} \left\{ F^h(\rho^h) + \frac{1}{2\tau s} W_2^2(\rho^h, \rho_k^{h,\tau}) \right\}. $$

The function $G_k$ has the following properties:

(i) $G_k(0^+) = F(\rho_k^{h,\tau})$ and $G_k(1) = F(\rho_{k+1}^{h,\tau}) + \frac{1}{2\tau} W_2^2(\rho_{k+1}^{h,\tau}, \rho_k^{h,\tau})$, where $\rho_{k+1}^{h,\tau}$ is the minimizer for (JKO);

(ii) $G_k$ is decreasing and, thus, differentiable almost everywhere;

(iii) $G'_k(s) = -\frac{1}{2\tau s^2} W_2^2(\tilde{\rho}_{(k+s)\tau}^{h,\tau}, \rho_k^{h,\tau})$ if $s \in (0, 1]$ is a point of differentiability.
Since $G_k$ is decreasing, we have

$$ \int_{0}^{1} G'_{k}(s) ds \geq G_{k}(1) - G_{k}(0^{+}), $$

therefore,

$$ (4.7) \qquad \mathcal{F}^{h}(\rho_{k+1}^{h,\tau}) - \mathcal{F}^{h}(\rho_{k}^{h,\tau}) + \frac{1}{2\tau} W_{2}^{2}(\rho_{k+1}^{h,\tau}, \rho_{k}^{h,\tau}) + \int_{0}^{1} \frac{1}{2\tau s^{2}} W_{2}^{2}(\tilde{\rho}_{(k+s)\tau}^{h,\tau}, \rho_{k}^{h,\tau}) ds \leq 0. $$

**Lemma 4.2 (Optimality condition).** An optimal measure $\tilde{\rho}_{(k+s)\tau}^{h,\tau}$ in (4.6) satisfies

$$ f'\left(\frac{\tilde{\rho}_{(k+s)\tau}^{h,\tau}(x)}{h^d}\right) + V(x) + \frac{\varphi^{h,\tau}(x)}{\tau s} = 0 \quad \text{for all } x \in \mathcal{T}^h, $$

where $\varphi^{h,\tau}$ is a Kantorovich potential corresponding to the optimal transport plan from $\tilde{\rho}_{(k+s)\tau}^{h,\tau}$ to $\rho_k^{h,\tau}$.

*Proof.* We use the notation

$$ \mathcal{J}_{\rho_0^h}^{h,\tau}(\rho^h) = \mathcal{F}^h(\rho^h) + \mathcal{W}_{\rho_0^h}^{h,\tau}(\rho^h), \quad \mathcal{W}_{\rho_0^h}^{h,\tau}(\rho^h) = \frac{1}{2\tau} W_2^2(\rho^h, \rho_0^h) $$

and the proof relies on the following dual formulation of the Wasserstein distance

$$ (4.8) \qquad W_2^2(\rho_1^h, \rho_0^h) = \sup_{\varphi^h \in L^1(\rho_1^h), \psi^h \in L^1(\rho_0^h)} \left\{ \sum_{x \in \mathcal{T}^h} \varphi^h(x) \rho_1^h(x) + \sum_{y \in \mathcal{T}^h} \psi^h(y) \rho_0^h(y) : \varphi^h \oplus \psi^h \le c \right\}, $$

with the cost function $c(x, y) = \frac{1}{2}|x - y|^2$.

We consider $\mathcal{J}_{\rho_0^h}^{h,\tau}$ as a functional on $L^\infty(\mathcal{T}^h)$. Given $\rho_0^h \in \mathcal{P}(\mathcal{T}^h)$, the functional $\mathcal{J}_{\rho_0^h}^{h,\tau}$ equals $+\infty$ on $L^\infty(\mathcal{T}^h) \setminus \mathcal{P}(\mathcal{T}^h)$, which one can show using (4.8) (or see [44, Proposition 7.17] for details).

Let $\rho_1^h$ be a minimizer of $\mathcal{J}_{\rho_0^h}^{h,\tau}$ that is $0 \in \partial \mathcal{J}_{\rho_0^h}^{h,\tau}(\rho_1^h)$. We claim that

$$ \partial \mathcal{J}_{\rho_0^h}^{h,\tau}(\rho_1^h) = \partial \mathcal{F}^h(\rho_1^h) + \partial \mathcal{W}_{\rho_0^h}^{h,\tau}(\rho_1^h). $$

To show this we use, for example, [17, Theorem 3.11], where the requirement is that there exists a point $v^h \in \text{dom}(\mathcal{F}^h) \cap \text{dom}(\mathcal{W}_{\rho_0^h}^{h,\tau})$ such that $\mathcal{F}^h$ is continuous in $v^h$. We choose $v^h$ to be the uniform normalized vector, i.e. $v^h(x) = 1/|\mathcal{T}^h|$ for all $x \in \mathcal{T}^h$. Indeed both $\mathcal{F}^h$ and $\mathcal{W}_{\rho_0^h}^{h,\tau}$ are finite in $v^h$ and $\mathcal{F}^h$ is continuous on $L^\infty(\mathcal{T}^h)$.

By [44, Proposition 7.17], the subdifferential of $W_2^2(\cdot, \rho_0^h)$ coincides with the set of Kantorovich potentials:

$$ \partial \mathcal{W}_{\rho_0^h}^{h,\tau}(\rho_1^h) = \left\{ \frac{\varphi^h}{\tau} \middle| \text{ for the Kantorovich potentials } \varphi^h \text{ from } \rho_1^h \text{ to } \rho_0^h \right\}. $$

For $\rho_1^h \in \mathcal{P}(\mathcal{T}^h)$, it holds that

$$ \left\{ f'\left(\frac{\rho_1^h}{h^d}\right) + V \right\} = \partial \mathcal{F}^h(\rho_1^h). $$

Thus, if $0 \in \partial \mathcal{J}_{\rho_0^h}^{h,\tau}(\rho_1^h)$, then we obtain the asserted optimality condition $\frac{\varphi^h}{\tau} + f'\left(\frac{\rho_1^h}{h^d}\right) + V = 0$ for some Kantorovich potentials $\varphi^h$ from $\rho_1^h$ to $\rho_0^h$. $\square$
We first consider the case where the driving energy is given by the internal energy without the external potential

$$ (4.9) \qquad \mathcal{F}^h(\rho^h) = \sum_{x \in \mathcal{T}^h} f\left(\frac{\rho_x^h}{h^d}\right) h^d. $$

We define a discrete counterpart for the Fisher information from Definition 2.1.

**Definition 4.3.** Let $\mathcal{F}^h$ be given as in (4.9). We define the discrete Fisher information $S^h : \mathcal{P}(\mathcal{T}^h) \to [0, +\infty)$ corresponding to $\mathcal{F}^h$ as

$$ S^h(\rho^h) = \frac{1}{4} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} \frac{|\ell(u_{x+\mathbf{h}}^h) - \ell(u_x^h)|^2}{h^2} h^d, \qquad u^h = \frac{\rho^h}{h^d}, $$

where $\ell$ is exactly the same as in Definition 2.1.

**Lemma 4.4.** Consider $\mathcal{F}^h$ defined in (4.9) and assume (f0). Let $\rho_0^h \in \mathcal{P}(\mathcal{T}^h)$ be given and $\rho^{h,\tau}$ be a minimizer of one step of (JKO$^{h,\tau}$). Then the following lower bound on the Wasserstein distance with the discrete Fisher information holds

$$ (4.10) \qquad \frac{1}{2\tau^2} W_2^2(\rho^{h,\tau}, \rho_0^h) \geq \left(1 - \frac{h}{2\tau}\right) S^h(\rho^{h,\tau}) - \frac{dh}{4\tau}. $$

*Proof.* Let $\gamma^{h,\tau}$ be an optimal coupling between the measures $\rho^{h,\tau}$ and $\rho_0^h$ and $(\varphi^{h,\tau}, \psi^{h,\tau})$ be the corresponding Kantorovich potentials. Let $(x, y) \in \text{spt}(\gamma^{h,\tau})$, then

$$ (4.11) \qquad \varphi^{h,\tau}(x) + \psi^{h,\tau}(y) = \frac{|x - y|^2}{2}. $$

For any $\mathbf{h} \in \mathbf{d}_h$, we have

$$ (4.12) \qquad \varphi^{h,\tau}(x + \mathbf{h}) + \psi^{h,\tau}(y) \leq \frac{|x + \mathbf{h} - y|^2}{2}. $$

Subtracting (4.11) from (4.12) and factorizing the difference of squares we obtain

$$ \varphi^{h,\tau}(x + \mathbf{h}) - \varphi^{h,\tau}(x) \leq \mathbf{h} \cdot (x - y + \mathbf{h}/2). $$

For what comes later, it is convenient to rearrange the terms as

$$ (4.13) \qquad (x - y) \cdot \frac{\mathbf{h}}{h} \geq \frac{\varphi^{h,\tau}(x + \mathbf{h}) - \varphi^{h,\tau}(x)}{h} - \frac{h}{2}. $$

If $(x - y) \cdot \mathbf{h} \geq 0$, it follows that

$$ (4.14) \qquad \left| (x - y) \cdot \frac{\mathbf{h}}{h} \right| \geq \frac{\left( \varphi^{h,\tau}(x + \mathbf{h}) - \varphi^{h,\tau}(x) \right)^{+}}{h} - \frac{h}{2}. $$

Choosing $\varphi^{h,\tau}$ so that we can apply the optimality condition from Lemma 4.2, we get

$$ \varphi^{h,\tau}(x + \mathbf{h}) - \varphi^{h,\tau}(x) = \tau \left( f'(u_x^{h,\tau}) - f'(u_{x+\mathbf{h}}^{h,\tau}) \right), $$

and it follows that

$$ (4.15) \qquad \left| (x - y) \cdot \frac{\mathbf{h}}{h\tau} \right| \geq \frac{\left( f'(u_x^{h,\tau}) - f'(u_{x+\mathbf{h}}^{h,\tau}) \right)^{+}}{h} - \frac{h}{2\tau}. $$
To deduce an inequality for square values from (4.15), we use that for any $a, b \ge 0$ such that $a \ge b - \epsilon$ for some $\epsilon > 0$, it holds that $a^2 \ge (1 - \epsilon)b^2 - \epsilon$ (see Lemma 4.5 below for the proof). Therefore, we obtain

$$ (4.16) \qquad \frac{|(x-y) \cdot \mathbf{h}|^2}{h^2 \tau^2} \ge \left(1 - \frac{h}{2\tau}\right) \left( \frac{\left(f'(u_x^{h,\tau}) - f'(u_{x+h}^{h,\tau})\right)^+}{h} \right)^2 - \frac{h}{2\tau}, \quad \text{if } (x-y) \cdot \mathbf{h} > 0. $$

In the other case, when $(x-y) \cdot \mathbf{h} < 0$, it means that $(x-y) \cdot (-\mathbf{h}) > 0$ and

$$ (x-y) \cdot \frac{-\mathbf{h}}{h} \ge \frac{(\varphi^{h,\tau}(x-\mathbf{h}) - \varphi^{h,\tau}(x))^+}{h} - \frac{h}{2}. $$

Following similar steps as for $(x-y) \cdot \mathbf{h} > 0$, we obtain counterparts of (4.14), (4.15), and (4.16) for $(x-y) \cdot \mathbf{h} < 0$. Eventually, we get

$$ (4.17) \qquad \frac{|(x-y) \cdot \mathbf{h}|^2}{h^2 \tau^2} \ge \left(1 - \frac{h}{2\tau}\right) \left( \frac{\left(f'(u_x^{h,\tau}) - f'(u_{x-h}^{h,\tau})\right)^+}{h} \right)^2 - \frac{h}{2\tau}, \quad \text{if } (x-y) \cdot \mathbf{h} < 0. $$

Summing (4.16) and (4.17) over $\mathbf{h} \in \mathbf{d}_h$, the left-hand side becomes

$$ (4.18) \qquad \sum_{\mathbf{h} \in \mathbf{d}_h} \frac{|(x-y) \cdot \mathbf{h}|^2}{h^2} = 2 \sum_{i=1}^{d} |(x-y) \cdot e_i|^2 = 2|x-y|^2 $$

and, consequently,

$$ \frac{|x-y|^2}{2\tau^2} \ge \frac{1}{4} \left(1 - \frac{h}{2\tau}\right) \sum_{\mathbf{h} \in \mathbf{d}_h} \left( \frac{\left(f'(u_x^{h,\tau}) - f'(u_{x+h}^{h,\tau})\right)^+}{h} \right)^2 - \frac{dh}{4\tau}. $$

Now we aim to use the estimate from the previous line to estimate the Wasserstein distance:

$$ \frac{1}{2\tau^2} W_2^2(\rho^{h,\tau}, \rho_0^h) = \frac{1}{2\tau^2} \sum_{(x,y) \in \mathcal{T}^h \times \mathcal{T}^h} |x-y|^2 \gamma^{h,\tau}(x,y). $$

Therefore,

$$ \frac{1}{2\tau^2} W_2^2(\rho^{h,\tau}, \rho_0^h) \ge \frac{1}{4} \left(1 - \frac{h}{2\tau}\right) \sum_{x \in \mathcal{T}^h} \rho_x^{h,\tau} \sum_{\mathbf{h} \in \mathbf{d}_h} \left( \frac{\left(f'(u_x^{h,\tau}) - f'(u_{x+h}^{h,\tau})\right)^+}{h} \right)^2 - \frac{dh}{4\tau}. $$

One can see the double sum over $\mathcal{T}^h$ and $\mathbf{d}_h$ as a sum over edges of the grid, where each edge is counted twice: as $(x, x + \mathbf{h})$ and $(x + \mathbf{h}, x)$. Hence the contribution to the sum of an undirected edge $(x, x + \mathbf{h})$ is

$$ \rho_x^{h,\tau} \left( \frac{\left(f'(u_x^{h,\tau}) - f'(u_{x+h}^{h,\tau})\right)^+}{h} \right)^2 + \rho_{x+h}^{h,\tau} \left( \frac{\left(f'(u_x^{h,\tau}) - f'(u_{x+h}^{h,\tau})\right)^+}{h} \right)^2 $$

Since we assumed that $f'$ is a monotonically increasing function, it means that the contribution of one edge equals to

$$ (4.19) \qquad \max(\rho_x^{h,\tau}, \rho_{x+h}^{h,\tau}) \frac{(f'(u_x^{h,\tau}) - f'(u_{x+h}^{h,\tau}))^2}{h^2} $$
and the lower bound on the sum can be rewritten as

$$ \sum_{x \in \mathcal{T}^h} \rho_x^{h, \tau} \sum_{\mathbf{h} \in \mathbf{d}_h} \left( \frac{(f'(u_x^{h, \tau}) - f'(u_{x+\mathbf{h}}^{h, \tau}))^+}{h} \right)^2 = \sum_{(x,y) \in \Sigma^h} \max(u_x^{h, \tau}, u_y^{h, \tau}) \frac{|f'(u_y^{h, \tau}) - f'(u_x^{h, \tau})|^2}{h^2} h^d. $$

Here we want to switch from dependence on $f'$ to $\ell$, where $\ell$ is defined up to an additive constant by the relation $\ell'(s) = \sqrt{s}f''(s)$. This change will be a discrete analog of the equality

$$ s|\nabla f'(s)|^2 = |\sqrt{s}f''(s)\nabla s|^2 = |\nabla \ell(s)|^2, $$

which is commonly used for different representations of the Fisher information. We assert that the following inequality holds true

$$ (4.20) \qquad \max(a, b)(f'(b) - f'(a))^2 \geq (\ell(b) - \ell(a))^2 \quad \text{for all } a, b \geq 0. $$

We can assume w.l.o.g. the inequality $b > a$: then, taking into account convexity of $f$, we get

$$ \sqrt{\max(a, b)(f'(b) - f'(a))} = \sqrt{b} \int_a^b f''(s) \, ds \geq \int_a^b \sqrt{s} f''(s) \, ds = \ell(b) - \ell(a). $$

Since both sides of the inequality are positive, we conclude

$$ \max(a, b)(f'(b) - f'(a))^2 \geq (\ell(b) - \ell(a))^2 $$

and we proved (4.20).

Applying (4.20) yields

$$ \sum_{(x,y) \in \Sigma^h} \max(u_x^{h, \tau}, u_y^{h, \tau}) \frac{|f'(u_y^{h, \tau}) - f'(u_x^{h, \tau})|^2}{h^2} h^d \geq \sum_{(x,y) \in \Sigma^h} \frac{|\ell(u_y^{h, \tau}) - \ell(u_x^{h, \tau})|^2}{h^2} h^d. $$

Summarizing the result, we obtain

$$ \frac{1}{2\tau^2} W_2^2(\rho^{h, \tau}, \rho_0^h) \geq \frac{1}{4} \left(1 - \frac{h}{2\tau}\right) \sum_{(x,y) \in \Sigma^h} \frac{|\ell(u_y^{h, \tau}) - \ell(u_x^{h, \tau})|^2}{h^2} h^d - \frac{dh}{4\tau}. \quad \square $$

**Lemma 4.5.** Let $a, b \geq 0$ be such that $a \geq b - \varepsilon$ for some $\varepsilon > 0$, then $a^2 \geq (1 - \varepsilon)b^2 - \varepsilon$.

*Proof.* If $b > \varepsilon$, then rearranging the terms gives

$$ a^2 \geq (b - \varepsilon)^2 = (1 - \varepsilon)b^2 + \varepsilon b^2 - 2\varepsilon b + \varepsilon^2 = (1 - \varepsilon)b^2 - \varepsilon + \varepsilon^2 + \varepsilon(b - 1)^2 $$

and the asserted inequality follows by dropping the positive term $\varepsilon^2 + \varepsilon(b - 1)^2$.

In the other case, if $b \in [0, \varepsilon]$, then

$$ (1 - \varepsilon)b^2 - \varepsilon \leq (1 - \varepsilon)\varepsilon^2 - \varepsilon = -\varepsilon\left[\left(\varepsilon - \frac{1}{2}\right)^2 + \frac{3}{4}\right] < 0 $$

and, trivially, $a^2 \geq 0 \geq 1 - \varepsilon$.

We now consider what changes if the driving energy is given by the internal energy plus an external potential energy as in (4.1). In this case, we define the discrete Fisher information as follows.
**Definition 4.6.** The discrete Fisher information $S_V^h : \mathcal{P}(\mathcal{T}^h) \to [0, +\infty)$ corresponding to $\mathcal{F}^h$ defined in (4.1) is given by

$$S_V^h(\rho^h) = \frac{1}{4} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} \frac{|\ell(u_{x+\mathbf{h}}^h) - \ell(u_x^h) + \sqrt{\max(u_x^h, u_{x+\mathbf{h}}^h)}(V(x+\mathbf{h}) - V(x))|^2}{h^2} h^d, \quad u^h = \frac{\rho^h}{h^d},$$

where $\ell$ is exactly the same as in Definition 2.1.

**Corollary 4.7.** Let the energy $\mathcal{F}^h$ be given by (4.1) with the potential $V \in \text{Lip}(\Omega)$ and let $f$ satisfy (f0) and (fV). Then, in the setting of Lemma 4.4, we have

$$(4.21) \qquad \frac{1}{2\tau^2} W_2^2(\rho^{h,\tau}, \rho_0^h) \ge \left(1 - \frac{h}{2\tau}\right) S_V^h(\rho^{h,\tau}) - e_h(\rho^{h,\tau}) - \frac{dh}{4\tau}.$$

with the Fisher information $S_h$ given in Definition 4.6 and the additional error term $e_h$ is such that $0 < e_h(\rho^{h,\tau}) \xrightarrow{h \to 0} 0$.

*Proof.* Notationwise, we omit the superscript $h, \tau$ in the proof and we distinguish the measure $\rho$ and its density $u_x = \rho_x/h^d$ for $x \in \mathcal{T}^h$.

The proof follows the similar steps as in the proof of Lemma 4.4. The main modification comes from the point that (4.19) does not have to hold anymore. The reason is that we cannot use the monotonicity of $f'$ directly because the potential $V$ can destroy the monotonicity property. Indeed, if we have a potential $V$, the contribution from one undirected edge $(x, y) \in \Sigma^h$ becomes

$$u_x \frac{[(f'(u_x) - f'(u_y) + V(x) - V(y))^+]^2}{h^2} + u_y \frac{[(f'(u_y) - f'(u_x) + V(y) - V(x))^+]^2}{h^2}.$$

Taking into account also the strict monotonicity of $f'$, meaning that $f'(u_y) - f'(u_x) = 0$ if and only if $u_x = u_y$, we can rewrite the contribution of an undirected edge as

$$\max(u_x, u_y) \frac{(f'(u_x) - f'(u_y) + V(x) - V(y))^2}{h^2} \mathbf{1}_{\{(f'(u_x) - f'(u_y) + V(x) - V(y))(f'(u_x) - f'(u_y)) \ge 0\}} \\ + \min(u_x, u_y) \frac{(f'(u_x) - f'(u_y) + V(x) - V(y))^2}{h^2} \mathbf{1}_{\{(f'(u_x) - f'(u_y) + V(x) - V(y))(f'(u_x) - f'(u_y)) < 0\}}.$$

We define the set of edges that violate the monotonicity property as

$$(4.22) \qquad \Sigma_{f,V}^h := \{(x, y) \in \Sigma^h : (f'(u_x) - f'(u_y) + V(x) - V(y))(f'(u_x) - f'(u_y)) < 0\}.$$

Summing up over all the edges and rearranging the terms give

$$\sum_{(x,y) \in \Sigma^h} \max(u_x, u_y) \frac{(f'(u_x) - f'(u_{x+\mathbf{h}}) + V(x) - V(y))^2}{h^2} h^d \\ + \sum_{(x,y) \in \Sigma^h_{f,V}} \left( \min(u_x, u_y) - \max(u_x, u_y) \right) \frac{(f'(u_x) - f'(u_y) + V(x) - V(y))^2}{h^2} h^d \\ = S^h(\rho^h) - e_h(u).$$

The error term is

$$e_h(u) := \sum_{(x,y) \in \Sigma^h_{f,V}} |u_x - u_y| \frac{(f'(u_x) - f'(u_y) + V(x) - V(y))^2}{h^2} h^d.$$
We take a closer look at the set $\Sigma_{f,V}^h$. If $f'(u_x) - f'(u_y) \ge 0$ and $f'(u_x) - f'(u_y) + V(x) - V(y) < 0$, then

$$0 \le f'(u_x) - f'(u_y) < V(y) - V(x) \le \mathrm{Lip}(V)h.$$

Consequently, for any $(x, y) \in \Sigma_{f,V}^h$, we can conclude that

$$(4.23) \qquad |f'(u_x) - f'(u_y)| \le \mathrm{Lip}(V)h$$

$$(4.24) \qquad |f'(u_x) - f'(u_y) + V(x) - V(y)| \le 2\mathrm{Lip}(V)h,$$

Therefore, the error term is bounded as

$$e_h(u) \le 2\mathrm{Lip}(V) \sum_{(x,y) \in \Sigma_{f,V}^h} |u_x - u_y| h^d.$$

Since $f'$ is strictly monotone, it has a well-defined inverse $g := (f')^{-1}$ and we can write

$$|u_x - u_y| = |g(f'(u_x)) - g(f'(u_y))| = \left| \int_{f'(u_x)}^{f'(u_y)} g'(s) \, ds \right|$$

We note that $f$ has to be defined only on $[0, \infty)$ and we extend it by $+\infty$ on $(-\infty, 0)$, then the Legendre dual $f^*$ satisfies the relation

$$(4.25) \qquad (f^*)'(s) = \begin{cases} g(s) & \text{for } s \ge f'(0), \\ 0 & \text{otherwise.} \end{cases}$$

We note that $f'(0)$ can equal $-\infty$, which means that $(f^*)'$ and $g$ coincide on the whole real line.

First we consider the case $\min(u_x, u_y) < s_0$. We claim that in this case $\max(u_x, u_y) < s_0+1$. Suppose the opposite, then by monotonicity of $f'$ we have

$$\mathrm{Lip}(V)h \ge |f'(u_x) - f'(u_y)| \ge f'(s_0 + 1) - f'(s_0) = f''(s_0 + \lambda_{s_0}),$$

with some $\lambda_{s_0} \in (0, 1)$. Employing (fV) again, we obtain

$$\mathrm{Lip}(V)h \ge f''(s_0 + \lambda_{s_0}) \ge \frac{1}{C_f(s_0 + 1)^\theta},$$

which cannot hold for small enough $h$, thus we have a contradiction and $|u_x - u_y| \le 1$. Therefore,

$$S_1^h := \sum_{(x,y) \in \Sigma_{f,V}^h} |u_x - u_y| h^d 1_{\{\min(u_x, u_y) < s_0\}} (x, y) \le \sum_{(x,y) \in \Sigma_{f,V}^h} |u_x - u_y| h^d 1_{\{\max(u_x, u_y) \le s_0+1\}} (x, y).$$

Since $(f^*)'$ is continuously differentiable and monotone on $[f'(0), f'(s_0 + 1)]$, we set $L_0 := \mathrm{Lip}((f^*)'_{|[f'(0), f'(s_0+1)]})$ and we obtain

$$\begin{align*}
S_1^h &= \sum_{(x,y) \in \Sigma_{f,V}^h} |(f^*)'(f'(u_x)) - (f^*)'(f'(u_y))| h^d 1_{\{\max(u_x, u_y) < s_0+1\}} (x, y) \\
&\le \sum_{(x,y) \in \Sigma_{f,V}^h} L_0 |f'(u_x) - f'(u_y)| h^d 1_{\{\max(u_x, u_y) < s_0+1\}} (x, y) \\
&\le L_0 \mathrm{Lip}(V) h 4d |\Omega| = o(1)|_{h \to 0}.
\end{align*}$$
If $\min(u_x, u_y) \ge s_0$, then $\min(f'(u_x), f'(u_y)) \ge f'(s_0) > -\infty$ and by the mean value theorem there exists $s_{xy} \in [f'(u_x), f'(u_y)]$ such that

$$|u_x - u_y| = \left| \int_{f'(u_x)}^{f'(u_y)} (f^*)''(s) ds \right| = (f^*)''(s_{xy}) |f'(u_x) - f'(u_y)|.$$

Since $f$ and $f^*$ are convex conjugates and twice continuously differentiable on the corresponding domains, they satisfy the relation

$$(f^*)''(f'(\sigma)) = \frac{1}{f''(\sigma)} \quad \text{for } \sigma \ge 0.$$

For any $s_{xy} \in [f'(u_x), f'(u_y)]$ there exists $\sigma_{xy} \in [u_x, u_y]$ such that $f'(\sigma_{xy}) = s_{xy}$, thus

$$|u_x - u_y| = \frac{1}{|f''(\sigma_{xy})|} |f'(u_x) - f'(u_y)| \le C_f \max(u_x, u_y)^\theta \mathrm{Lip}(V)h,$$

where we used (fV) to estimate $f''$. We notice that $\max_{x \in \mathcal{T}^h} u_x \le 1/h^d$, therefore,

$$\max(u_x, u_y)^\theta \le \max(u_x, u_y) \frac{1}{h^{d(\theta-1)}} \le \frac{u_x + u_y}{h^{d(\theta-1)}}$$

Combining the estimates for $(x, y) \in \Sigma_{f,V}^h$ such that $\min(u_x, u_y) \ge s_0$, we get

$$S_2^h := \sum_{(x,y) \in \Sigma_{f,V}^h} |u_x - u_y| h^d 1_{\{\min(u_x, u_y) \ge s_0\}} (x, y) \\ \le C_f \mathrm{Lip}(V) h \sum_{(x,y) \in \Sigma_{f,V}^h} \frac{u_x + u_y}{h^{d(\theta-1)}} h^d 1_{\{\min(u_x, u_y) \ge p\}} (x, y) \\ \le 2C_f \mathrm{Lip}(V) h^{1-d(\theta-1)} = o(1)|_{h \to 0},$$

where we take into account the assumption $1 - d(\theta - 1) > 0$ coming from (fV). This concludes the proof because we got

$$e_h(u) \le S_1^h + S_2^h = o(1)|_{h \to 0}.$$

<figure><img src="image_3.png" /></figure>

**Remark 4.8.** The linear Fokker-Planck case is easier to handle than the approach in the proof of Corollary 4.7. Indeed, the main point in the proof is the following: for two nearby points $x$ and $y$ in the grid, we assume $f'(u_x) + V(x) \ge f'(u_y) + V(y)$ and we need to bound from below $u_x$ with $\max\{u_x, u_y\}$. If $u_x \ge u_y$ there is nothing to prove, otherwise the computations of Corollary 4.7 exactly aim at obtaining a similar bound up to an error term. In the case where $f'(s) = \log s$, then we simply obtain $u_x \ge u_y e^{-\mathrm{Lip}(V)h}$, which allows to obtain the result up to a multiplicative factor which tends to 1 as $h \to 0$. Unfortunately, the computations are less explicit when we lose the logarithmic structure of the Fokker-Planck equation.

**Remark 4.9.** It might be possible to modify Lemma 4.4 for a more general finite-volume tessellation. Let $x$ and $z$ be generating points of the finite-volume cells $K_x$ and $K_z$ and $h_{zx} := z - x$. We also use a common assumption that $z - x$ is orthogonal to the interface $f_{zx} := \overline{K_x} \cap \overline{K_z}$. If we can define the transition kernel as $\kappa_{zx}$ in such a way that (4.18) becomes

$$\sum_{z \sim x} \kappa_{zx} \frac{|(x - y) \cdot \mathbf{h}_{zx}|^2}{|\mathbf{h}_{zx}|^2} = |x - y|^2,$$
then the conclusion should follow.

**Remark 4.10.** As mentioned in Remark 2.10, the proof of Corollary 4.7 can be adapted to deal with the interaction potential W instead of the external potential V. For this purpose, one needs to replace V(x) with the discrete convolution $\sum_{z \in \mathcal{T}^h} W(x-z)\rho_z^h$ in the definition of the Fisher information and, consequently, in the proof. The main point, why the proof works similarly for the interaction potential, is to notice that

$$ \left| \sum_{z \in \mathcal{T}^h} (W(x+h-z) - W(x-z)) \rho_z^h \right| \le \mathrm{Lip}(W)h \sum_{z \in \mathcal{T}^h} \rho_z^h = \mathrm{Lip}(W)h. $$

We are now in a good position to return to the inequality (4.7) and turn it into a form suitable for passing to the limit $h \to 0$. For this purpose, we need another type of interpolation in time, namely piecewise geodesic interpolation.

Let a sequence $\{\rho_k^{h,\tau}\}_{k=1,...,N}$ be the sequence of $(JKO^{h,\tau})$ minimizers as in Definition 2.6. Since $P_2(\Omega)$ is a geodesic space, all pairs $(\rho_k^{h,\tau}, \rho_{k+1}^{h,\tau})$, $k = 0, ..., N-1$ can be connected by a (possibly not unique) geodesic $[0, 1] \ni t \mapsto \tilde{\rho}_k^{h,\tau}(t)$. Then there exists a Borel time-dependent vector field $(0, 1) \ni t \mapsto \tilde{v}_k^{h,\tau}(t) \in L^2(\rho_k^{h,\tau}(t); \mathbb{R}^d)$ such that the continuity equation

$$ \partial_t \tilde{\rho}_k^{h,\tau} + \mathrm{div}(\tilde{\rho}_k^{h,\tau} \tilde{v}_k^{h,\tau}) = 0 $$

holds in the sense of distributions. These considerations allow us to introduce piecewise geodesic interpolation.

**Definition 4.11** (Piecewise geodesic interpolation). The piecewise geodesic interpolation is a curve $[0, T] \ni t \mapsto \tilde{\rho}_t^{h,\tau}$ defined as

$$ \tilde{\rho}_t^{h,\tau} := \rho_k^{h,\tau} \left( \frac{t-k\tau}{\tau} \right) \quad t \in [k\tau, (k+1)\tau), $$

where $\rho_k^{h,\tau}(t)$ is a geodesic connecting $\rho_k^{h,\tau}$ to $\rho_{k+1}^{h,\tau}$. Moreover, we denote by $\tilde{v}^{h,\tau}$ the velocity field defined as

$$ \tilde{v}_t^{h,\tau} := v_k^{h,\tau} \left( \frac{t-k\tau}{\tau} \right) \quad t \in (k\tau, (k+1)\tau). $$

To bring in the integration in time, we will also make use of the following measure

$$ (4.26) \qquad \lambda_{\epsilon}^{h,\tau}(dt) := \sum_{k=0}^{N-1} 1_{[(k+\epsilon)\tau, (k+1)\tau]}(t) \left(1 - \frac{h}{2(t-k\tau)}\right) \mathcal{L}(dt), $$

defined for an arbitrary $\epsilon \in (0, 1)$.

**Lemma 4.12.** Let $h > 0$, $\tau > 0$, and $\rho_0^h \in \mathcal{P}(\mathcal{T}^h)$ be given. Let also $\{\rho_k^{h,\tau}\}_k$ be the sequence of $(JKO^{h,\tau})$ minimizers from Definition 2.6. Then for any $\epsilon \in (0, 1)$, the following variational inequality holds

$$ 0 \ge \mathcal{F}^h(\rho_N^{h,\tau}) - \mathcal{F}^h(\rho_0^h) + \frac{1}{2} \int_0^T \|\tilde{v}_t^{h,\tau}\|_{L^2(\tilde{\rho}_t^{h,\tau})}^2 dt + \int_0^T S_V^h(\tilde{\rho}_t^{h,\tau}) \lambda_{\epsilon}^{h,\tau}(dt) + \frac{dT}{4\tau} \log \epsilon + o(1)_{h \to 0}. $$

*Proof.* A slight modification of the proof of Lemma 4.4 with $\tau s$ instead of $\tau$ shows that for the variational interpolant $\tilde{\rho}_{st}^{h,\tau}$ defined in Definition 4.1 we have

$$ (4.27) \qquad \frac{1}{2\tau^2 s^2} W_2^2(\tilde{\rho}_{s\tau}^{h,\tau}, \rho_0^h) \ge \left(1 - \frac{h}{2\tau}\right) S_V^h(\tilde{\rho}_{s\tau}^{h,\tau}) - \frac{dh}{4\tau s} + o(1)|_{h \to 0}. $$
Multiplying the latter inequality by $\tau$ and integrating over $s$, we get

$$ \frac{\tau}{2} \int_{\epsilon}^{1} \frac{1}{\tau^2 s^2} W_2^2(\tilde{\rho}_{(k+s)\tau}^{h,\tau}, \rho_k^{h,\tau}) ds \geq \tau \int_{\epsilon}^{1} \left(1 - \frac{h}{2s\tau}\right) S_V^h(\tilde{\rho}_{(k+s)\tau}^{h,\tau}) ds + \frac{dh}{4} \log \epsilon + \tau(1-\epsilon)\mathcal{o}(1)|_{h \to 0}, $$

for an arbitrary $\epsilon \in (0, 1)$. Then we obtain from (4.7)

$$ 0 \geq \mathcal{F}^h(\rho_{k+1}^{h,\tau}) - \mathcal{F}^h(\rho_k^{h,\tau}) + \frac{1}{2\tau} W_2^2(\rho_{k+1}^{h,\tau}, \rho_k^{h,\tau}) + \tau \int_{\epsilon}^{1} \left(1 - \frac{h}{2s\tau}\right) S_V^h(\tilde{\rho}_{(k+s)\tau}^{h,\tau}) ds \\ + \frac{dh}{4} \log \epsilon + \tau(1-\epsilon)\mathcal{o}(1)|_{h \to 0}. $$

Rescaling the time as

$$ \tau \int_{\epsilon}^{1} \left(1 - \frac{h}{2s\tau}\right) S_V^h(\tilde{\rho}_{(k+s)\tau}^{h,\tau}) ds = \int_{\epsilon\tau}^{\tau} \left(1 - \frac{h}{2s}\right) S_V^h(\tilde{\rho}_{k\tau+s}^{h,\tau}) ds, $$

and summing up from $k=0$ to $N-1$, we obtain that

$$ 0 \geq \mathcal{F}^h(\rho_N^{h,\tau}) - \mathcal{F}^h(\rho_0^{h}) + \frac{\tau}{2} \sum_{k=0}^{N-1} \frac{1}{\tau^2} W_2^2(\rho_{k+1}^{h,\tau}, \rho_k^{h,\tau}) + \int_{\epsilon\tau}^{\tau} \left(1 - \frac{h}{2s}\right) \sum_{k=0}^{N-1} S_V^h(\tilde{\rho}_{k\tau+s}^{h,\tau}) ds \\ + \frac{dT}{4\tau} \log \epsilon + o(1)|_{h \to 0}. $$

To switch to the time integral from 0 to $T$ in the term with the Fisher information, we use the measure $\lambda_{\epsilon}^{h,\tau}$ defined in (4.26):

$$ \int_{\epsilon\tau}^{\tau} \left(1 - \frac{h}{2s}\right) \sum_{k=0}^{N-1} S_V^h(\tilde{\rho}_{k\tau+s}^{h,\tau}) ds = \int_{0}^{T} S_V^h(\tilde{\rho}_t^{h,\tau}) \lambda_{\epsilon}^{h,\tau}(dt). $$

Now we need to rewrite the term $\frac{\tau}{2} \sum_{k=0}^{N-1} \frac{1}{\tau^2} W_2^2(\rho_{k+1}^{h,\tau}, \rho_k^{h,\tau})$ in an integral form. For this purpose, we will employ the piecewise geodesic interpolation. For $t \in (k\tau, (k+1)\tau)$, we have

$$ \|\tilde{v}_t^{h,\tau}\|_{L^2(\tilde{\rho}_t^{h,\tau})} = \frac{1}{\tau} W_2(\rho_{k+1}^{h,\tau}, \rho_k^{h,\tau}), $$

therefore, with the appropriate time rescaling, we obtain

$$ \frac{\tau}{2} \sum_{k=0}^{N-1} \frac{1}{\tau^2} W_2^2(\rho_{k+1}^{h,\tau}, \rho_k^{h,\tau}) = \frac{1}{2} \sum_{k=0}^{N-1} \int_{k\tau}^{(k+1)\tau} \int_{\Omega} |\tilde{v}_t^{h,\tau}|^2 d\tilde{\rho}_t^{h,\tau} dt = \frac{1}{2} \int_0^T \int_{\Omega} |\tilde{v}_t^{h,\tau}|^2 d\tilde{\rho}_t^{h,\tau} dt. $$

And the asserted variational inequality holds.

<figure><img src="image_4.png" /></figure>

**4.2. Compactness results.** We begin the lemma on the convergence of the discrete curves to the limit continuous curve. This step is standard for proving the convergence of various versions of the JKO schemes with the strategy explained, for instance, in [44, Section 8.3].

**Lemma 4.13.** The curves $\tilde{\rho}^{h,\tau}$ and $\tilde{\rho}^{h,\tau}$ defined in Definitions 4.1 and 4.11, respectively, converge up to a subsequence to the same limit curve $\rho \in C([0, T]; P(\Omega))$ as $h, \tau \to 0$ uniformly in time for the $W_2$ distance.

*Proof.* Using $\rho_k^{h,\tau}$ as a competitor for one step of (JKO$^{h,\tau}$), we have

$$ \mathcal{F}^h(\rho_{k+1}^{h,\tau}) + \frac{1}{\tau} W_2^2(\rho_k^{h,\tau}, \rho_{k+1}^{h,\tau}) \leq \mathcal{F}^h(\rho_k^{h,\tau}). $$
Summing up the inequalities for $k = 0, \dots, N-1$ gives

$$ \sum_{k=0}^{N-1} \frac{W_2^2(\rho_k^{h,\tau}, \rho_{k+1}^{h,\tau})}{\tau} \le \sum_{k=0}^{N-1} (\mathcal{F}^h(\rho_k^{h,\tau}) - \mathcal{F}^h(\rho_{k+1}^{h,\tau})) = \mathcal{F}^h(\rho_0^h) - \mathcal{F}(\hat{\rho}_N^h) \le \sup_{h>0} \mathcal{F}^h(\rho_0^h) - \inf_{\rho \in \mathcal{P}(\Omega)} \mathcal{F}(\rho) =: C_F. $$

We can use the latter inequality to bound the norm of the velocity $\tilde{v}^{h,\tau}$:

$$ \int_0^T \|\tilde{v}_t^{h,\tau}\|_{L^2(\tilde{\rho}_t^{h,\tau})}^2 dt = \sum_{k=0}^{N-1} \frac{1}{\tau} W_2^2(\rho_{k+1}^{h,\tau}, \rho_k^{h,\tau}) \le C_F. $$

Denoting by $\tilde{j}_t^{h,\tau}$ the corresponding flux $\tilde{\rho}_t^{h,\tau} = \tilde{\rho}_t^{h,\tau} \tilde{v}_t^{h,\tau}$, we also obtain

$$ \int_0^T |\tilde{j}_t^{h,\tau}|(\Omega) dt = \int_0^T \|\tilde{v}_t^{h,\tau}\|_{L^1(\tilde{\rho}_t^{h,\tau})} dt \le \sqrt{T} \left( \int_0^T \|\tilde{v}_t^{h,\tau}\|_{L^2(\tilde{\rho}_t^{h,\tau})}^2 dt \right)^{1/2} \le \sqrt{C_F T}. $$

Thus, the family $\{\tilde{j}_t^{h,\tau}\}_{h,\tau>0}$ is weakly-*$\delta$ compact in $\mathcal{M}([0,T] \times \Omega; \mathbb{R}^d)$. Since $(\tilde{\rho}^{h,\tau}, \tilde{j}^{h,\tau})$ satisfies the continuity equation, $t \mapsto \tilde{\rho}_t^{h,\tau}$ is an absolutely continuous curve satisfying for $s < t$

$$ \begin{aligned} W_2(\tilde{\rho}_t^{h,\tau}, \tilde{\rho}_s^{h,\tau}) &\le \int_s^t |(\tilde{\rho}_r^{h,\tau})'| dr \le \int_s^t \|\tilde{v}_r^{h,\tau}\|_{L^1(\tilde{\rho}_r^{h,\tau})} dr \\ &\le \sqrt{t-s} \left( \int_s^t \|\tilde{v}_r^{h,\tau}\|_{L^2(\tilde{\rho}_r^{h,\tau})}^2 dr \right)^{1/2} \le \sqrt{C_F(t-s)}. \end{aligned} $$

Notice that the estimates above hold for an arbitrary $h > 0$, so we can choose $h, \tau \to 0$. Therefore, we can apply the Ascoli–Arzelá theorem, and there exists a (not relabeled) subsequence $\{\tilde{\rho}^{h(\tau),\tau}\}$ and a limit curve $\rho \in C([0, T]; P(\Omega))$, such that

$$ \tilde{\rho}_t^{h,\tau} \to \rho_t \quad \text{uniformly for the } W_2 \text{ distance.} $$

Furthermore, the curves $\tilde{\rho}^{h,\tau}$ and $\tilde{\rho}^{h,\tau}$ converge to the same curve $\rho_t$. We use that the curves $\tilde{\rho}^{h,\tau}$ and $\tilde{\rho}^{h,\tau}$ coincide at $t = k\tau$ for all $k = 0, \dots, N-1$. For $t \in ((k-1)\tau, k\tau]$, we have

$$ \begin{aligned} W_2(\tilde{\rho}_t^{h,\tau}, \tilde{\rho}_t^{h,\tau}) &\le W_2(\tilde{\rho}_t^{h,\tau}, \tilde{\rho}_{k\tau}^{h,\tau}) + W_2(\tilde{\rho}_{k\tau}^{h,\tau}, \tilde{\rho}_t^{h,\tau}) \le \sqrt{C_F(t-k\tau)} + W_2(\rho_k^{h,\tau}, \rho_{k-1}^{h,\tau}) + W_2(\rho_{k-1}^{h,\tau}, \tilde{\rho}_t^{h,\tau}) \\ &\le 3\sqrt{C_F\tau}. \end{aligned} $$

Therefore,

$$ \tilde{\rho}_t^{h,\tau} \to \rho_t \quad \text{uniformly for the } W_2 \text{ distance.} \quad \square $$

**Remark 4.14.** There is no need to assume $h/\tau \to 0$ to obtain compactness. In fact, the compactness result holds for any relation between parameters $h$ and $\tau$. However, having $h/\tau \to 0$ will be crucial for proving that the limit curve $\rho_t$ is a solution of the corresponding PDE. If $h/\tau \to 0$ is not satisfied, it is possible to construct examples where $\rho_k^{h,\tau} = \rho_0^h$ for all $k = 1, \dots, N-1$ (the evolution is frozen) or where the error terms result in additional drift added to the original equation.

Now we aim to prove the strong compactness in space of the density. To simplify the notation, we define

$$ (4.28) \qquad \hat{\ell}^h := \sum_{x \in \mathcal{T}^h} \ell(u_x^h) 1_{Q_h(x)}. $$

To prove the crucial compactness result in Lemma 4.16, we need a version of the Poincaré-Wirtinger inequality presented in the following lemma. Several versions of discrete Poincaré-type inequalities are known for finite volume schemes [27]. This proof is a simple modification of the standard argument used for such inequalities (as in [27, Lemma 3.7]).
**Lemma 4.15** (Discrete Poincaré-Wirtinger inequality). Given $f^h \in B(\mathcal{T}^h)$, let $\hat{f}^h$ be the piecewise constant reconstruction defined in (2.12) and set

$$ \bar{f}^h := \int_{\Omega} \hat{f}^h(x) dx. $$

Then the following inequality holds

$$ \|\hat{f}^h - \bar{f}^h\|_2^2 \le C_{\Omega} \sum_{(p,q) \in \Sigma^h} |f_p^h - f_q^h|^2 h^{d-2}, $$

where $C_{\Omega} > 0$ is independent of $h$.

*Proof.* Fix an arbitrary $x \in \Omega$. There exists $p \in \mathcal{T}^h$ such that $x \in Q_h(p)$. For any $z \in \Omega$ such that $z \ni Q_h(p)$, there exists $q \in \mathcal{T}^h$ such that $x \in Q_h(q)$ and we can find a path between $p$ and $q$ consisting of the directed edges $(p_i, p_{i+1}) \in \Sigma^h$ such that $[x, z] \cap I_h(p_i, p_{i+1}) \neq \emptyset$, where $[x, z]$ is the line segment connecting $x$ and $z$. We denote by path$(x, z)$ this sequence of edges

$$ \text{path}(x, z) = ((p_0 = p, p_1), (p_1, p_2), \dots, (p_{K-1}, p_K = q)), \quad p_i \in \mathcal{T}^h, \quad (p_i, p_{i+1}) \in \Sigma^h. $$

To obtain a bound on the number of edges in the path, we use the following idea. A point $p' \in \mathcal{T}^h$ cannot be in the path if the distance from $p'$ to the line segment $[x, z]$ is more than $2h$. Thus, we can get $\text{bound } |\text{path}(x, z)| = K$ by comparing the volume of the set $H([x, z], 2h) := \{y \in \Omega : \text{dist}(y, [x, z]) < 2h\}$ with the volume of one cell $h^d$.

$$ (4.29) \qquad K \le \frac{|H([x, z], 2h)|}{h^d} \le C_d \frac{(|x-z| + 2h)h^{d-1}}{h^d} \le C_{\text{path}} \frac{|x-z|}{h}, $$

where the constant $C_d$ is the volume of $(d-1)$-dimensional ball with the radius 1.

Using first the triangular inequality and then the Hölder inequality, we get

$$ |\hat{f}^h(x) - \hat{f}^h(z)| \le \sum_{(p,q) \in \text{path}(x,z)} |f_p^h - f_q^h| \le \left( \sum_{(p,q) \in \text{path}(x,z)} |f_p^h - f_q^h|^2 \frac{1}{h} \right)^{1/2} \left( \sum_{(p,q) \in \text{path}(x,z)} h \right)^{1/2}. $$

By (4.29), we have

$$ |\hat{f}^h(x) - \hat{f}^h(z)|^2 \le C|x-z| \sum_{(p,q) \in \Sigma^h} |f_p^h - f_q^h|^2 \frac{1}{h} \{[x, z] \cap I_h(p, q) \ne \emptyset\}. $$

Fix an arbitrary $(p, q) \in \Sigma^h$. We denote by $R$ the diameter of $\Omega$ and by $B_R$ the ball with the radius $R$ centered at the origin. Integrating in $x$ and $z$ over $\Omega$ and using the change of variables $\eta = x - z$, we obtain:

$$ \int_{\Omega} \int_{\Omega} |x-z| \mathbf{1}\{[x, z] \cap I_h(p, q) \neq \emptyset\} dx dz \le \int_{B_R} \int_{\Omega} |\eta| \mathbf{1}\{[x, x+\eta] \cap I_h(p, q) \neq \emptyset\} dx d\eta \\ \le \int_{B_R} |\eta|^2 h^{d-1} d\eta \le C_R h^{d-1}, $$

where $C_R < \infty$ is a constant depending only on $R$. Therefore,

$$ \int_{\Omega} \int_{\Omega} |\hat{f}^h(x) - \hat{f}^h(z)|^2 dx dz \le C_{d,R} \sum_{(p,q) \in \Sigma^h} |f_p^h - f_q^h|^2 h^{d-2}. $$
To conclude the proof, it is just left to use the Jensen inequality:

$$
\begin{aligned}
\| \hat{f}^h - \bar{f}^h \|_2^2 &= \left| \int_{\Omega} |\hat{f}^h(x) - \int_{\Omega} \hat{f}^h(z) dz | \right|^2 dx \le \int_{\Omega} \left| \int_{\Omega} \hat{f}^h(x) - \hat{f}^h(z) \right|^2 dz dx \\
&\le C_{\Omega} \sum_{(p,q) \in \Sigma^h} |f_p^h - f_q^h|^2 h^{d-2}.
\end{aligned}
$$

<figure><img src="image_5.png" /></figure>

**Lemma 4.16.** Let a family ${\rho^h}_{h>0}$ with $\rho^h \in \mathcal{P}(\mathcal{T}^h)$ be such that $\hat{\rho}^h \to \rho$ narrowly in $\mathcal{P}(\Omega)$ and

$$
\sup_{h>0} S^h(\rho^h) := C_S < \infty.
$$

Then for any $\eta \in \mathbb{R}^d$ we have

$$
(4.30) \quad \sup_{h>0} \| \hat{\ell}^h \|_{L^2(\Omega)} := C_2 < \infty \quad \text{and} \quad \int_{\Omega_{|\eta|}} |\hat{\ell}^h(z-\eta) - \hat{\ell}^h(z)|^2 dz \le C|\eta| \max(|\eta|, h),
$$

where $\Omega_{\epsilon} = \{x \in \Omega : \operatorname{dist}(x, \partial\Omega) \ge \epsilon\}$.

Moreover, ${\{\hat{\ell}^h\}}_{h>0}$ is relatively compact in $L^2(\Omega)$, any subsequential limit $g$ of ${\{\hat{\ell}^h\}}_{h>0}$ belongs to $H^1(\Omega)$, and there exists a subsequence for which the limit has the form $g = \ell(\frac{d\rho}{d\mathcal{L}^d})$.

*Proof.* We first aim to establish the strong compactness in $L^2(\Omega)$ of the family ${\{\hat{\ell}^h - \bar{\ell}^h\}}_{h>0}$. By Lemma 4.15, we get the uniform bound for $L^2$ norm:

$$
\sup_{h>0} \| \hat{\ell}^h - \bar{\ell}^h \|_2^2 \le \sup_{h>0} C_{\Omega} S^h(\rho^h) \le C_{\Omega} C_S.
$$

To apply the Riesz-Fréchet-Kolmogorov theorem for the family ${\{\hat{\ell}^h - \bar{\ell}^h\}}_{h>0}$, we need to prove that for $\eta \in \mathbb{R}^d$ the value of $\|\hat{\ell}^h(\cdot - \eta) - \hat{\ell}^h\|_2$ is vanishing uniformly as $|\eta| \to 0$. The approach to this proof is similar to the one in Lemma 4.15. Notice that

$$
\int_{\Omega_{|\eta|}} |\hat{\ell}^h(z - \eta) - \hat{\ell}^h(z)|^2 dz = \sum_{x \in \mathcal{T}^h} \int_{Q_h(x) \cap \Omega_{|\eta|}} |\hat{\ell}^h(z - \eta) - \hat{\ell}^h(z)|^2 dz
$$

For a fixed $z \in \Omega$, there exists a unique pair $x, y \in \mathcal{T}^h$ such that $z \in Q_h(x)$ and $(z - \eta) \in Q_h(y)$.
Furthermore, we can find a path between $x$ and $y$ by following a chain of vertices $\{z_i\}_{i=0,\dots,K}$ such that
$Q_h(z_i) \cap [z, z - \eta]$, where $[z, z - \eta]$ is the line segment between $z$ and $z - h$. By construction, $z_0 = x$
and $z_K = y$. Therefore, using Jensen's inequality, we get

$$
\begin{aligned}
\sum_{x \in \mathcal{T}^h} \int_{Q_h(x) \cap \Omega_{|\eta|}} |\hat{\ell}^h(z - \eta) - \hat{\ell}^h(z)|^2 dz &= \sum_{x \in \mathcal{T}^h} \int_{Q_h(x) \cap \Omega_{|\eta|}} \left| \sum_{i=0}^{K-1} (\ell^h(z_{i+1}) - \ell^h(z_i)) \right|^2 dx \\
&\le \sum_{x \in \mathcal{T}^h} \int_{Q_h(x) \cap \Omega_{|\eta|}} K \sum_{i=0}^{K-1} |\ell^h(z_{i+1}) - \ell^h(z_i)|^2 dx \\
&\le \sum_{x \in \mathcal{T}^h} \int_{Q_h(x) \cap \Omega_{|\eta|}} K \sum_{z \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} |\ell^h(z + \mathbf{h}) - \ell^h(z)|^2 1_{H([x,x-\eta],2h)}(z) dx,
\end{aligned}
$$
where $H([x, x - \eta], 2h) = \{p \in \mathbb{R}^d : \text{dist}([x, x - \eta], p) < 2h\}$. The bound on the number of cells in the path was found in (4.29)

$$K \leq C_{\text{path}} \frac{|\eta| + h}{h}.$$

Therefore,

$$\begin{align*}
\sum_{x \in \mathcal{T}^h} \int_{Q_h(x)} K \sum_{z \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} |\ell^h(z + \mathbf{h}) - \ell^h(z)|^2 & 1_{H([x,x-\eta],3h)}(z) \,dx \\
&\leq C_{\text{path}} \frac{|\eta| + h}{h} \sum_{z \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} |\ell^h(z + \mathbf{h}) - \ell^h(z)|^2 \int_{\Omega_{|\eta|}} 1_{H([x,x-\eta],3h)}(z) \,dx \\
&\leq C_{\text{path}} \frac{|\eta| + h}{h} C_{d-1}(|\eta| + 2h)(2h)^{d-1} \sum_{z \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} |\ell^h(z + \mathbf{h}) - \ell^h(z)|^2 \\
&\leq C(|\eta| + h)^2 S^h(\rho^h).
\end{align*}$$

To sum up, we obtain

$$\int_{\Omega_{|\eta|}} |\ell^h(z - \eta) - \ell^h(z)|^2 dz \leq C(|\eta| + h)^2 S^h(\rho^h).$$

The Riesz-Fréchet-Kolmogorov theorem [14, Theorem 4.26] provides relative compactness of $\{\ell^h - \overline{\ell}^h\}_{h>0}$ in $L^2(\Omega)$. Let $g$ be an arbitrary subsequential limit of $\{\ell^h - \overline{\ell}^h\}_{h>0}$. An application of (4.30) gives

$$\int_{\Omega_{|\eta|}} |g(x - \eta) - g(x)|^2 dx = \lim_{h \to 0} \int_{\Omega_{|\eta|}} |\ell^h(x - \eta) - \ell^h(x)|^2 dx \leq C C_S |\eta|^2,$$

which implies that $g \in H^1(\Omega)$ by the characterization of $H^1$ in terms of difference quotients [35, Theorem 11.75].

We claim that $\{\overline{\ell}^h\}_{h>0}$ is uniformly bounded, i.e. $\inf_{h>0} \overline{\ell}^h > -\infty$ and $\sup_{h>0} \overline{\ell}^h < \infty$. First, we prove the upper bound. Suppose the opposite that there exists a (not relabeled) sequence such that $\overline{\ell}^h \to \infty$. The $L^2$ convergence $\ell^h - \overline{\ell}^h \to g$ implies that there exists a further subsequence such that $\ell^h(x) - \overline{\ell}^h \to g(x)$ pointwise a.e. on $\Omega$. Thus, together we obtain that $\ell^h(x) \to \infty$ for a.e. $x \in \Omega$. Since $\ell^h = \ell(\hat{u}^h)$ with $\ell$ being continuous and monotonically non-decreasing function, $\ell^h(x) \to \infty$ for a.e. $x \in \Omega$ implies $\hat{u}^h(x) \to \infty$ for a.e. $x \in \Omega$, which is a contradiction with $\int_{\Omega} \hat{u}^h(x) dx = 1$ for all $h > 0$.

For the lower bound, suppose that there exists a (not relabeled) subsequence such that $\overline{\ell}^h \to -\infty$. Arguing similarly to the upper bound, we arrive at the claim that there exists a subsequence such that $\ell^h(x) \to -\infty$ for a.e. $x \in \Omega$. Then we would deduce $\hat{u}^h(x) \to 0$ for a.e. $x \in \Omega$. Together with the equi-integrability of $\hat{u}^h$ (guaranteed by (f1)), this would be again a contradiction to the condition $\int_{\Omega} \hat{u}^h(x) dx = 1$.

Since $\{\overline{\ell}^h\}_{h>0}$ is uniformly bounded, the family $\{\ell^h\}_{h>0}$ is compact in $L^2(\Omega)$ with any subsequen- cial limit being in $H^1(\Omega)$. Again, one can get the pointwise a.e. convergence up to subsequence, and using that $\ell$ has a well-defined inverse, we conclude $\hat{u}^h(x) \to (\ell^{-1}\circ g)(x)$ a.e. on $\Omega$. On the other hand, we have $\rho^h \to \rho$ narrowly in $P(\Omega)$. Since the weak limit is unique, we obtain $\rho = (\ell^{-1}\circ g)S^d$.
Therefore, for a.e. $x \in \Omega$, it holds that

$$u(x) := \frac{d\rho}{d\mathcal{L}^d}(x) = \limsup_{r \to 0} \int_{B_r(x)} (\ell^{-1} \circ g)(z) dz = \ell^{-1}(g(x)).$$

Consequently, $g(x) = \ell(u(x))$ a.e. on $\Omega$.

<figure><img src="image_6.png" /></figure>

**Lemma 4.17.** Let $V \in \text{Lip}(\Omega)$ and $\{\rho^h\}_{h>0}$ be a family of probability measures on $\mathcal{P}(\mathcal{T}^h)$. The uniform bound on the Fisher information with the potential $V$ is equivalent to the uniform bound on the Fisher information without the potential, i.e.,

$$\sup_{h>0} S_V^h(\rho^h) < \infty \iff \sup_{h>0} S^h(\rho^h) < \infty.$$

*Proof.* Since the part of the Fisher information depending on the potential is uniformly bounded

$$\sum_{(x,y) \in \Sigma^h} (V(y) - V(x))^2 \max(u_x^h, u_y^h) h^{d-2} \le \text{Lip}^2(V) \sum_{(x,y) \in \Sigma^h} (u_x^h + u_y^h) h^d \le \text{Lip}^2(V),$$

the proof follows by the Young inequality.

<figure><img src="image_7.png" /></figure>

**4.3. Convergence results.** To complete the proof of Theorem 2.9, we need to consider the liminf inequality for the Fisher information. Various results on discrete-to-continuum $\Gamma$-convergence of quadratic functionals that apply to our setting are available in the literature, see for instance [2]. We present a simple proof of the liminf inequality for completeness.

We define an approximation for $\nabla \ell(u)$ in the following way

$$(4.31) \qquad \overline{\nabla}^h \hat{\ell}^h := \frac{1}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} 1_{Q_h(x+\mathbf{h}/2)} \frac{\mathbf{h}}{h} \frac{\hat{\ell}^h(x+\mathbf{h}) - \hat{\ell}^h(x)}{h}.$$

Note that we use the family of cells $\{Q_h(x+\mathbf{h}/2)\}_{x \in \mathcal{T}^h}$ which are dual to the family of cells $\{Q_h(x)\}_{x \in \mathcal{T}^h}$ used for the reconstruction of $\hat{\ell}^h$.

**Lemma 4.18 (Convergence of discrete gradients).** Let a family $\{\rho^h\}_{h>0}$ with $\rho^h \in \mathcal{P}(\mathcal{T}^h)$ be such that $\hat{\rho}^h \to \rho$ in narrowly $\mathcal{P}(\Omega)$ and

$$\sup_{h>0} S^h(\rho^h) := C_S < \infty.$$

Let $\{\hat{\ell}^h\}_{h>0}$ and $\{\overline{\nabla}^h \hat{\ell}^h\}_{h>0}$ be defined as in (4.28) and (4.31) respectively. Then, up to a subsequence,

$$\hat{\ell}^h \to \ell(u) \quad \text{strongly in } L^2(\Omega),$$

$$\overline{\nabla}^h \hat{\ell}^h \to \nabla \ell(u) \quad \text{weakly in } L^2(\Omega),$$

where $u = d\rho/d\mathcal{L}^d \in L^1(\Omega)$. Furthermore, we have

$$\liminf_{h \to 0} S_V^h(\rho^h) \ge S_V(\rho).$$

*Proof.* The asserted convergence $\hat{\ell}^h \to \ell(u)$ follows from Lemma 4.16.

For the convergence of the gradients $\overline{\nabla}^h \hat{\ell}^h \to \nabla \ell(u)$, we first show that $\overline{\nabla}^h \hat{\ell}^h \in L^2(\Omega)$ with the $L^2$-norm equal to the discrete Fisher information defined in Definition 4.3 up to a multiplicative constant.
Indeed, we calculate

$$
(4.32) \quad \begin{aligned}
\| \bar{\nabla}^h \hat{\ell}^h \|_{L^2}^2 &= \frac{1}{4} \int_{\Omega} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} 1_{Q_h(x+\mathbf{h}/2)}(z) \frac{|\hat{\ell}^h(x+\mathbf{h}) - \hat{\ell}^h(x)|^2}{h^2} \mathrm{d}z \\
&= \frac{1}{4} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} h^d \frac{|\hat{\ell}^h(x+\mathbf{h}) - \hat{\ell}^h(x)|^2}{h^2} = S^h(\rho^h).
\end{aligned}
$$

The assumed uniform bound on $S^h$ implies that there exists a (non-relabeled) subsequence and $\zeta \in L^2(\Omega)$ such that $\bar{\nabla}^h \hat{\ell}^h \to \zeta$ weakly in $L^2(\Omega)$.

To identify the limit $\zeta$, we show an approximate integration-by-parts formula for $\hat{\ell}^h$ and $\bar{\nabla}^h \hat{\ell}^h$. Let $\varphi \in C_c^\infty(\Omega; \mathbb{R}^d)$. On the one hand, we have

$$ \int (\nabla \cdot \varphi) \hat{\ell}^h = \sum_{x \in \mathcal{T}^h} \ell(u_x^h) \int_{Q_h(x)} (\nabla \cdot \varphi)(z) \, dz = \sum_{x \in \mathcal{T}^h} \ell(u_x^h) \sum_{\mathbf{h} \in \mathbf{d}_h} \int_{I(x, \mathbf{h})} \varphi(z) \cdot \frac{\mathbf{h}}{h} H^{d-1}(dz), $$

where we used the divergence theorem. We denote $\varphi_{x, \mathbf{h}}^h := f_{I(x, \mathbf{h})} \varphi(z) H^{d-1}(dz)$ and note that $\varphi_{x, \mathbf{h}}^h = \varphi_{x+\mathbf{h}, -\mathbf{h}}^h$ to obtain

$$
\begin{aligned}
\int (\nabla \cdot \varphi) \hat{\ell}^h &= \sum_{x \in \mathcal{T}^h} h^{d-1} \sum_{\mathbf{h} \in \mathbf{d}_h} \ell(u_x^h) \varphi_{x, \mathbf{h}}^h \cdot \frac{\mathbf{h}}{h} = \frac{1}{2} \sum_{x \in \mathcal{T}^h} h^{d-1} \sum_{\mathbf{h} \in \mathbf{d}_h} \left( \ell(u_x^h) \varphi_{x, \mathbf{h}}^h \cdot \frac{\mathbf{h}}{h} - \ell(u_{x+\mathbf{h}}^h) \varphi_{x, \mathbf{h}}^h \cdot \frac{\mathbf{h}}{h} \right) \\
&= -\frac{1}{2} \sum_{x \in \mathcal{T}^h} h^{d-1} \sum_{\mathbf{h} \in \mathbf{d}_h} \left( \ell(u_{x+\mathbf{h}}^h) - \ell(u_x^h) \right) \int_{Q_h(x+\mathbf{h}/2)} \varphi_{x, \mathbf{h}}^h \cdot \frac{\mathbf{h}}{h} \, dz.
\end{aligned}
$$

On the other hand, by the definition of $\bar{\nabla}^h \hat{\ell}^h$, we have

$$ \int \varphi \cdot \bar{\nabla}^h \hat{\ell}^h = \frac{1}{2} \sum_{x \in \mathcal{T}^h} h^{d-1} \sum_{\mathbf{h} \in \mathbf{d}_h} (\ell(u_{x+\mathbf{h}}^h) - \ell(u_x^h)) \int_{Q_h(x+\mathbf{h}/2)} \varphi(z) \cdot \frac{\mathbf{h}}{h} \, dz. $$

Since $\varphi$ is smooth, there exists $C_\varphi > 0$ such that $|\varphi_{x, \mathbf{h}}^h - \varphi(z)| \le C_\varphi h$ for all $z \in Q_h(x + \mathbf{h}/2)$ for any $x \in \mathcal{T}^h$ and $\mathbf{h} \in \mathbf{d}_h$. Therefore,

$$
\begin{aligned}
\left| \int (\nabla \cdot \varphi) \hat{\ell}^h + \int \varphi \cdot \bar{\nabla}^h \hat{\ell}^h \right| &\le \sum_{x \in \mathcal{T}^h} \frac{h^{d-1}}{2} \sum_{\mathbf{h} \in \mathbf{d}_h} |\ell(u_{x+\mathbf{h}}^h) - \ell(u_x^h)| \int_{Q_h(x+\mathbf{h}/2)} |\varphi(z) - \varphi_{x,\mathbf{h}}^h| \mathrm{d}z \\
&\le \frac{C_\varphi}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} |\ell(u_{x+\mathbf{h}}^h) - \ell(u_x^h)| h^d \\
&\le \frac{C_\varphi}{2} \left( \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} |\ell(u_{x+\mathbf{h}}^h) - \ell(u_x^h)|^2 h^{d-2} \right)^{1/2} \left( \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} h^{d+2} \right)^{1/2} \\
&\le \frac{C_\varphi}{2} \sqrt{S^h(\rho^h)} h \sqrt{2d|\Omega|}.
\end{aligned}
$$

Using again the uniform bound on $S^h(\rho^h)$, we obtain the approximate integration-by-parts formula:

$$ \int (\nabla \cdot \varphi) \hat{\ell}^h = - \int \varphi \cdot \bar{\nabla}^h \hat{\ell}^h + O(h). $$
Passing $h \to 0$ yields

$$ \int (\nabla \cdot \varphi) \ell(u) = - \int \varphi \cdot \zeta \quad \text{for any } \varphi \in C_c^\infty(\Omega; \mathbb{R}^d). $$

Since $\varphi$ is arbitrary, we obtain $\zeta = \nabla \ell(u)$.

We conclude that the liminf inequality for $S^h(\rho^h)$ holds due to (4.32) and the lower-semicontinuity of the $L^2$ norm:

$$ \liminf_{h \to 0} S^h(\rho^h) = \liminf_{h \to 0} \| \overline{\nabla}^h \hat{\ell}^h \|_{L^2}^2 \ge \| \nabla \ell(u) \|_{L^2}^2. $$

In the case where the external potential $V$ is present, we have to consider the other term that we denote by $\hat{\xi}^h \in L^2(\Omega; \mathbb{R}^d)$:

$$ \begin{aligned} \hat{\xi}_i^h &:= \hat{\mu}_i^h \overline{\nabla}^h \hat{V}^h \quad i = 1, \dots, d; \qquad \hat{\mu}_i^h := \sum_{x \in \mathcal{T}^h} 1_{Q_h(x+he_i/2)} \sqrt{\max(u_x^h, u_{x+he_i}^h)}, \\ \overline{\nabla}^h \hat{V}^h &:= \frac{1}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} 1_{Q_h(x+\mathbf{h}/2)} \frac{\mathbf{h}}{h} \frac{V(x+\mathbf{h}) - V(x)}{h}, \end{aligned} $$

that gives us $S_V^h = \| \overline{\nabla}^h \hat{\ell}^h + \hat{\xi}^h \|_{L^2}^2$.

The compactness result in Lemma 4.16 implies that up to a subsequence $\hat{u}^h$ converges a.e. to $u$. Together with the equi-integrability of $\hat{u}^h$ (which comes from the superlinearity assumption on $f$), this implies $\hat{u}^h \to u$ strongly in $L^1(\Omega)$. This in turns implies $\sqrt{\hat{u}^h} \to \sqrt{u}$ strongly in $L^2(\Omega)$. The, we denote $v_i^h := \hat{u}^h(\cdot + he_i) = \sum_{x \in \mathcal{T}^h} 1_{Q_h(x)} u_{x+he_i}^h$. We claim that $\max(\sqrt{\hat{u}^h}, \sqrt{v_i^h}) \to \sqrt{u}$ strongly in $L^2$ for all $i = 1, \dots, d$. Indeed,

$$ \begin{aligned} \| \max(\sqrt{\hat{u}^h}, \sqrt{v_i^h}) - \sqrt{u} \|_{L^2} &\le \| \sqrt{v_i^h} - \sqrt{\hat{u}^h} \|_{L^2} + \| \sqrt{\hat{u}^h} - \sqrt{u} \|_{L^2} \\ &\le \| \sqrt{\hat{u}^h}(\cdot + he_i) - \sqrt{u}(\cdot + he_i) \|_{L^2} + \| \sqrt{u}(\cdot + he_i) - \sqrt{u} \|_{L^2} + \| \sqrt{\hat{u}^h} - \sqrt{u} \|_{L^2} \\ &= \| \sqrt{u}(\cdot + he_i) - \sqrt{u} \|_{L^2} + \| \sqrt{\hat{u}^h} - \sqrt{u} \|_{L^2} \end{aligned} $$

where $\| \sqrt{u}(\cdot + he_i) - \sqrt{u} \|_{L^2} \to 0$ as $h \to 0$ holds by the density argument ([14, Lemma 4.3]). Thus, $\max(\sqrt{\hat{u}^h}, \sqrt{v_i^h}) \to \sqrt{u}$ strongly in $L^2$ and, consequently, $\hat{\mu}_i^h \to \sqrt{u}$ strongly in $L^2$ for all $i = 1, \dots, d$. Using a similar argument as in the first part of the proof, one can show that $\overline{\nabla}^h \hat{V}^h \to \nabla V$ weakly in $L^\infty$. Therefore, $\hat{\xi}^h \to \sqrt{u} \nabla V$ weakly in $L^2$ and we conclude that we have

$$ \liminf_{h \to 0} S_h^V \ge \| \nabla \ell(u) + \sqrt{u} \nabla V \|_{L^2}^2 = S^V. \quad \square $$

**Remark 4.19.** In continuation of Remarks 2.10 and 4.10 on including the interaction potential $W$ instead of the external potential $V$, we note that the proof of Lemma 4.18 can be adapted in the following way. One can show that

$$ (\widehat{W * \rho})^h := \sum_{x \in \mathcal{T}^h} 1_{Q_h(x)} \sum_{z \in \mathcal{T}^h} W(x-z) \rho_z^h \xrightarrow{h \to 0} W * \rho \quad \text{strongly in } L^1(\Omega), $$
using the regularity of $W$ and strong convergence $\hat{u}^h \to u$ in $L^1(\Omega)$. Secondly, one can also prove that the approximation for the gradient of the convolution

$$ \overline{\nabla}^h (\widehat{W * \rho})^h := \frac{1}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} 1_{Q_h(x+\mathbf{h}/2)} \frac{\mathbf{h}}{h} \frac{\sum_{z \in \mathcal{T}^h} (W(x+\mathbf{h}-z) - W(x-z)) \rho_z^h}{h} $$

converges weakly in $L^\infty$ to $\nabla(W * \rho)$. This can be done by means of the approximate integration-by-parts formula as in the first part of the proof, where the key estimate is

$$ \begin{aligned} \left| \int (\nabla \cdot \varphi)(\widehat{W * \rho})^h + \int \varphi \cdot \overline{\nabla}^h (\widehat{W * \rho})^h \right| & \le \frac{C_\varphi}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} \left| \sum_{z \in \mathcal{T}^h} (W(x+\mathbf{h}-z) - W(x-z)) \rho_z^h \right| h^d \\ & \le C_\varphi \mathrm{Lip}(W) |\Omega| dh \xrightarrow{h \to 0} 0. \end{aligned} $$

Finally, we have all the ingredients to summarize the proof of Theorem 2.9.

*Proof.* We start with the variational inequality derived in Lemma 4.12

$$ (4.33) \quad 0 \ge \mathcal{F}^h(\rho_N^{h,\tau}) - \mathcal{F}^h(\rho_0^h) + \frac{1}{2} \int_0^T \|\check{v}_t^{h,\tau}\|_{L^2(\rho_t^{h,\tau})}^2 dt + \int_0^T S_V^h(\tilde{\rho}_t^{h,\tau}) \lambda_\epsilon^{h,\tau}(dt) + \frac{dT}{\tau} \log \epsilon + o(1) \Big|_{h \to 0}. $$

We aim to pass $h, \tau, h/\tau \to 0$ (we will simply write $h, \tau \to 0$ implying that $h/\tau \to 0$ holds) in all components of (4.33).

We recall that the discrete energy can be written in the integral form as

$$ \mathcal{F}^h(\rho_T^{h,\tau}) = \int_{\Omega} f \left( \frac{d\hat{\rho}_T^{h,\tau}}{d\mathcal{L}^d} \right) dx + \int_{\Omega} \hat{V}^h(x) \hat{\rho}_T^{h,\tau}(dx). $$

Since $\rho_T^{h,\tau} \to \rho_T$ narrowly as $h, \tau \to 0$, then

$$ \liminf_{h, \tau \to 0} \int_{\Omega} f \left( \frac{d\hat{\rho}_T^{h, \tau}}{d\mathcal{L}^d} \right) dx \ge \int_{\Omega} f \left( \frac{d\rho_T}{d\mathcal{L}^d} \right) dx $$

and

$$ \begin{aligned} \lim_{h, \tau \to 0} \int_{\Omega} \hat{V}^h(x) \hat{\rho}_T^{h, \tau}(dx) &= \lim_{h, \tau \to 0} \int_{\Omega} V(x) \hat{\rho}_T^{h, \tau}(dx) + \lim_{h, \tau \to 0} \int_{\Omega} (\hat{V}^h(x) - V(x)) \hat{\rho}_T^{h, \tau}(dx) \\ &= \int_{\Omega} V(x) \rho_T(dx) + \lim_{h \to 0} \mathrm{Lip}(V)h = \int_{\Omega} V(x) \rho_T(dx). \end{aligned} $$

Using also the assumption on the initial data, we get

$$ \liminf_{h, \tau \to 0} (\mathcal{F}^h(\rho_T^{h,\tau}) - \mathcal{F}^h(\rho_0^{h,\tau})) \ge \mathcal{F}(\rho_T) - \mathcal{F}(\rho_0). $$

For the Fisher information, we claim that, for any fixed $\epsilon > 0$, we have

$$ \liminf_{h, \tau \to 0} \int_0^T S_V^h(\tilde{\rho}_t^{h,\tau}) \lambda_\epsilon^{h,\tau}(dt) \ge (1-\epsilon) \int_0^T S_V(\rho_t) dt. $$

This can be obtained by combining the two following facts:

* first, by Lemma 4.13 and Lemma 4.18, we have, for any sequence $t_n \to t$

$$ \liminf_{h, \tau \to 0, n \to \infty} S_V^h(\tilde{\rho}_{t_n}^{h,\tau}) \ge S(\rho_t); $$

* second, we have the narrow convergence $\lambda_\epsilon^{h,\tau} \to (1-\epsilon)\mathcal{L}$ as $h, \tau \to 0$.
these two facts together allow to conclude via a standard argument proven, for instance, in [13, Proposition 5.5].

By the lower semicontinuity of the Benamou-Brenier functional, we have

$$ \liminf_{h, \tau \to 0} \| \check{v}_t^{h, \tau} \|_{L^2(\rho_t^{h, \tau})}^2 \ge \|v_t\|_{L^2(\rho_t)}^2, $$

where we applied the convergence established in Lemma 4.13.

To summarize, in the limit $h, \tau \to 0$, we recover

$$ 0 \ge \mathcal{F}(\rho_T) - \mathcal{F}(\rho_0) + \int_0^T \left\{ \frac{1}{2} \int_{\Omega} |v_t|^2 d\rho_t + (1 - \varepsilon) \mathcal{S}(\rho_t) \right\} dt, $$

and the EDI inequality that we want to obtain as in Proposition 2.3 can be recovered taking $\varepsilon \to 0$. $\square$

## 5. CROWD MOTION MODEL

In this section, we study the convergence of the fully discrete JKO scheme for the crowd motion model presented in Section 2.1.2. We recall that the driving energy in this case becomes

$$ (5.1) \qquad \mathcal{F}_{\text{CM}}^{h}(\rho^{h}) = \begin{cases} \sum_{x \in \mathcal{T}^{h}} V(x) \rho_{x}^{h}, & \text{if } u^{h} \le 1 \\ +\infty, & \text{otherwise.} \end{cases} $$

for a probability measure $\rho^h = u^h h^d$ supported on the grid $\mathcal{T}^h$. The fully discrete JKO scheme reads

**Definition 5.1.** For a given $h > 0$ and $\tau > 0$ of the from $\tau = T/N$ for a fixed $T > 0$ and some $N \in \mathbb{N}$, we define

$$ (\text{JKO}_{\text{cm}}^{h, \tau}) \qquad \rho_{k+1}^{h, \tau} \in \arg \min_{\rho^h \in \mathcal{P}(\mathcal{T}^h)} \left\{ \mathcal{F}_{\text{CM}}^h(\rho^h) + \frac{1}{2\tau} W_2^2(\rho^h, \rho_k^{h, \tau}) \right\} $$

with a given initial datum $\rho_0^h \in \mathcal{P}(\mathcal{T}^h)$.

The main result of this section is the following convergence result for $(\text{JKO}_{\text{cm}}^{h, \tau})$.

**Theorem 5.2.** Assume $V \in C^1(\mathbb{R}^d)$. Let $\{\rho_k^{h, \tau}\}_{k=1,...,N}$ be the family of $(\text{JKO}_{\text{cm}}^{h, \tau})$ minimizers as in Definition 5.1. Let the family of initial data $\{\rho_0^h\}_{h>0}$ be such that there exists $\rho_0 = u_0 \mathcal{L}^d \in \mathcal{P}(\Omega)$ with $\|u_0\|_{L^\infty} \le 1$ such that

$$ \rho_0^h \to \rho_0 \quad \text{narrowly as } h \to 0 \quad \text{and} \quad \lim_{h \to 0} \mathcal{F}_{\text{CM}}(\rho_0^h) = \mathcal{F}_{\text{CM}}(\rho_0). $$

Then a suitable interpolation in time of $\{\rho_k^{h, \tau}\}_{k=1,...,N}$ converges as $h, \tau, h/\tau \to 0$ uniformly in time for the $W_2$ distance to an absolutely continuous curve $[0, T] \ni t \mapsto \rho_t \in \mathcal{P}(\Omega)$ such that $\rho_t$ is a distributional solution of (2.7).

For what follows, it is convenient to define upwind-type mean:

$$ (5.2) \qquad \Lambda_{p+V}(u_x, u_y) := \begin{cases} u_x, & p_x + V(x) \ge p_y + V(y) \\ u_y, & p_x + V(x) < p_y + V(y) \end{cases} \quad (x, y) \in \Sigma^h. $$

We use the notation $\Lambda_V(u_x, u_y)$ if the average depends only on the potential $V$ ($p=0$ in (5.2)).
**Definition 5.3.** Let $(u^h, p^h)$ be a density-pressure pair satisfying the relation $p_x^h(1-u_x^h)$ for all $x \in \mathcal{T}^h$. Given a smooth potential $V$, we define the discrete crowd-motion Fisher information as

$$S_{\text{CM}}^{h}(u^{h}, p^{h}) = \frac{1}{4} \sum_{(x,y) \in \Sigma^{h}} \left( (p_{x}^{h} - p_{y}^{h})^{2} + 2(p_{x}^{h} - p_{y}^{h})(V(x) - V(y)) + (V(x) - V(y))^{2} \Lambda_{V}(u_{x}^{h}, u_{y}^{h}) \right) h^{d-2}.$$

For $\varepsilon > 0$, we also introduce the truncation function $g_{\varepsilon}(s) := (s - \varepsilon)^{+}$, where $s^{+}$ denotes the positive part of $s$, e.i. $s^{+} = s1_{s>0}$.

**Lemma 5.4.** Let $\mathcal{F}_{\text{CM}}^h$ be defined as in (5.1) with $V \in C^1(\mathbb{R}^d)$. Let $\rho_0^h \in \mathcal{P}(\mathcal{T}^h)$ be given and $\rho^{h,\tau}$ be a minimizer of $(\text{JKO}_{\text{cm}}^{h,\tau})$ for some $h, \tau > 0$. Then the following lower bound holds

$$\frac{1}{2\tau^2} W_2^2(\rho^{h,\tau}, \rho_0^h) \geq \left(1 - \frac{h}{2\tau}\right) S_{\text{CM}}^h(u^{h,\tau}, g_{\varepsilon}(h)(p^{h,\tau})) - \frac{dh}{4\tau} + o(1)|_{h \to 0},$$
with $\varepsilon(h) = (\text{Lip}(V) + 1)h$.

*Proof.* There exists a Kantorovich potential $\varphi$ corresponding to the optimal transport from $\rho^{h,\tau}$ to $\rho_0^h$ such that

$$\sum_{x \in \mathcal{T}^h} \left( V(x) + \frac{\varphi_x}{\tau} \right) \rho_x^h \geq \sum_{x \in \mathcal{T}^h} \left( V(x) + \frac{\varphi_x}{\tau} \right) \rho_x^{h, \tau} \quad \text{for all } \rho^h \in \mathcal{P}(\mathcal{T}^h),$$

by an argument as in [39, Lemma 3.1]. Thus, for a fixed $\varphi$, $\rho^{h,\tau} = u^{h,\tau}h^d$ solves the minimization problem

$$\rho^{h,\tau} \in \arg \min_{\rho^h \in \mathcal{P}(\mathcal{T}^h), u^h \le 1} \sum_{x \in \mathcal{T}^h} \left( V(x) + \frac{\varphi_x}{\tau} \right) \rho_x^h.$$

The minimizer of this linear minimization problem with a linear constraint has to have the form

$$u_x^{h,\tau} = \begin{cases} 1, & V(x) + \varphi_x/\tau < c \\ \in [0, 1], & V(x) + \varphi_x/\tau = c \\ 0, & V(x) + \varphi_x/\tau > c, \end{cases}$$

where $c \in \mathbb{R}$ is chosen in the way to ensure $\sum_{x \in \mathcal{T}^h} \rho_x^{h,\tau} = 1$. Employing the Kantorovich potential $\varphi$, one can define a pressure-like function $p$ (as in [39, Lemma 3.3]):

$$(5.3) \qquad p := \left( c - V - \frac{\varphi}{\tau} \right)^{+}.$$

By a similar arguments as in Lemma 4.4, one can derive the inequality:

$$\frac{1}{2\tau^2} W_2^2(\rho^{h,\tau}, \rho_0^h) \geq \frac{1}{2} \left(1 - \frac{h}{2\tau}\right) \sum_{(x,y) \in \Sigma^h} \frac{(p_x - p_y + V(x) - V(y))^2}{h^2} \Lambda_{p+V}(u_x, u_y) h^d - \frac{dh}{4\tau},$$

where we use the notation $\Lambda_{p+V}(u_x, u_y)$ introduced in (5.2) and $(u, p)$ is the density-flux pair corresponding to $\rho^{h,\tau}$, i.e. $\rho^{h,\tau} = uh^d$ and $p(1-u) = 0$ on $\mathcal{T}^h$.

The approach of Corollary 4.7 is not applicable in this case, because the proof relies on the strict monotonicity of $f'$ that the pressure variable does not satisfy. We provide a different argument exploiting the properties of the pressure.

Consider an edge $(x, y) \in \Sigma^h$ and assume (without loss of generality) that $p_x + V(x) > p_y + V(y)$. The contribution of this edge to the sum above is

$$\frac{(p_x - p_y + V(x) - V(y))^2}{h^2} u_x h^d = \frac{(p_x - p_y)^2}{h^2} u_x h^d + 2 \frac{(p_x - p_y)(V(x) - V(y))}{h^2} u_x h^d + \frac{(V(x) - V(y))^2}{h^2} u_x h^d.$$
Fix $\varepsilon > 0$ such that $\varepsilon \ge (\mathrm{Lip}(V) + 1)h$. We claim that we have

$$ (5.4) \qquad \frac{(p_x - p_y)^2}{h^2} u_x h^d \ge \frac{(g_\varepsilon(p_x) - g_\varepsilon(p_y))^2}{h^2} h^d. $$

To show (5.4), we first use that $g_\varepsilon$ is 1-Lipshitz:

$$ (p_x - p_y)^2 u_x \ge (g_\varepsilon(p_x) - g_\varepsilon(p_y))^2 u_x, $$

and we can remove the coefficient $u_x$ from the right-hand side if we have $p_x > 0$ because of the density-pressure relation $p_x(1 - u_x) = 0$. If, instead, we assume $p_x = 0$, then, by the assumption on the edge $(x, y) \in \Sigma^h$, we have $p_y < V(x) - V(y) \le \mathrm{Lip}(V)h < \varepsilon$ and $g_\varepsilon(p_y) = 0$. Moreover, $p_x = 0$ implies $g_\varepsilon(p_x) = 0$ so that we have

$$ (p_x - p_y)^2 u_x \ge (g_\varepsilon(p_x) - g_\varepsilon(p_y))^2 = 0 $$

and we get (5.4). Thus, summing up over all the edges yields

$$ (5.5) \qquad \sum_{(x,y) \in \Sigma^h} \frac{(p_x - p_y)^2}{h^2} \Lambda_{p+V}(u_x, u_y) h^d \ge \frac{1}{2} \sum_{(x,y) \in \Sigma^h} \frac{(g_\varepsilon(p_x) - g_\varepsilon(p_y))^2}{h^2} h^d. $$

Consider the cross-term:

$$ (p_x - p_y)(V(x) - V(y)) u_x = (p_x - p_y)(V(x) - V(y)) + (u_x - 1)(p_x - p_y)(V(x) - V(y)). $$

If $p_x > 0$, then $u_x = 1$, and the last term equals 0. In the case when $p_x = p_y = 0$, the last term vanishes as well. It is left to consider the case when $p_x = 0$ and $p_y > 0$. Since we assume that $p_x + V(x) > p_y + V(y)$, then $V(x) > V(y)$, and it follows that the last term is positive. Thus, we get the inequality

$$ \begin{aligned} (p_x - p_y)(V(x) - V(y)) u_x &\ge (p_x - p_y)(V(x) - V(y)) \\ &= (g_\varepsilon(p_x) - g_\varepsilon(p_y))(V(x) - V(y)) + (p_x - g_\varepsilon(p_x) - p_y + g_\varepsilon(p_y))(V(x) - V(y)). \end{aligned} $$

Since the cross-term is symmetric in for edges $(x, y) \in \Sigma^h$, it will sum up as

$$ \begin{aligned} & \left| \sum_{(x,y) \in \Sigma^h} (p_x - g_\varepsilon(p_x) - p_y + g_\varepsilon(p_y))(V(x) - V(y)) h^{d-2} \right| \\ &= 2 \left| \sum_{x \in \mathcal{T}^h} (p_x - g_\varepsilon(p_x)) \sum_{h \in \mathbf{d}_h} (V(x) - V(x+h)) h^{d-2} \right| \\ &\le 2\varepsilon \sum_{x \in \mathcal{T}^h} \left| \sum_{h \in \mathbf{d}_h} (V(x) - V(x+h)) \right| h^{d-2}. \end{aligned} $$

Since $V \in C^1(\mathbb{R}^d)$, then for any $x \in \mathcal{T}^h$ away from the boundary the discrete divergence of the gradient is bounded as:

$$ \begin{aligned} \left| \sum_{h \in \mathbf{d}_h} (V(x+h) - V(x)) \right| &= \left| \sum_{i=1}^{d} (V(x+he_i) + V(x+he_i) - 2V(x)) \right| \\ &= \left| \sum_{i=1}^{d} \left( \int_{0}^{1} \nabla V(x + \lambda he_i) \, d\lambda \cdot (he_i) - \int_{0}^{1} \nabla V(x + (1-\lambda)he_i) \, d\lambda \cdot (he_i) \right) \right| \\ &\le \sum_{i=1}^{d} \int_{0}^{1} \left| \nabla V(x + \lambda he_i) - \nabla V(x + (1-\lambda)he_i) \right| \, d\lambda \, h \le d\omega(h)h, \end{aligned} $$
where $\omega$ is the modulus of continuity of $\nabla V$. Therefore, sum up, we obtain

$$ \left| \sum_{x \in \mathcal{T}^h} (p_x - g_\epsilon(p_x)) \sum_{h \in \mathcal{d}_h} (V(x) - V(x+h)) h^{d-2} \right| \le d\omega(h) \frac{\epsilon}{h} \sum_{x \in \mathcal{T}^h} h^d = d|\Omega|\omega(h) \frac{\epsilon}{h}. $$

In total, the cross-term enjoys the following lower bound:

$$ (5.6) \qquad 2 \sum_{(x,y) \in \Sigma^h} \frac{(p_x - p_y)(V(x) - V(y))}{h^2} \Lambda_{p+V}(u_x, u_y) h^d \\ \qquad \qquad \qquad \qquad \qquad \ge \sum_{(x,y) \in \Sigma^h} \frac{(g_\epsilon(p_x) - g_\epsilon(p_y))(V(x) - V(y))}{h^2} h^d + 2d|\Omega|\omega(h) \frac{\epsilon}{h}. $$

To rewrite the third term depending on the potential, we claim that $\Lambda_{p+V}(u_x, u_y) \ge \Lambda_V(u_x, u_y)$, where, analogously to (5.2), $\Lambda_V(u_x, u_y)$ is defined as

$$ \Lambda_V(u_x, u_y) := \begin{cases} u_x, & V(x) \ge V(y) \\ u_y, & V(x) < V(y). \end{cases} $$

To prove the claim, we consider the cases. It is clear that if $p_x = p_y = 0$, then $\Lambda_V(u_x, u_y) = \Lambda_{p+V}(u_x, u_y)$. The equality also holds if $p_x > 0$ and $p_y > 0$ because it follows that $u_x = u_y = 1$ and $\Lambda_V(u_x, u_y) = 1 = \Lambda_{p+V}(u_x, u_y)$. The last case is (w.l.o.g.) $p_x > 0$ and $p_y = 0$. In this case, two options are possible: either $p_x + V(x) \ge p_y + V_y$ and $\Lambda_{p+V}(u_x, u_y) = u_x = 1 \ge \Lambda_V(u_x, u_y)$, or $p_x + V(x) < p_y + V_y$ implying that $V(x) < V_y$ and $\Lambda_{p+V}(u_x, u_y) = u_y = \Lambda_V(u_x, u_y)$. Therefore,

$$ (5.7) \qquad \sum_{(x,y) \in \Sigma^h} \frac{(V(x) - V(y))^2}{h^2} \Lambda_{p+V}(u_x, u_y) h^d \ge \sum_{(x,y) \in \Sigma^h} \frac{(V(x) - V(y))^2}{h^2} \Lambda_V(u_x, u_y) h^d. $$

Combining (5.5), (5.6), and (5.7), we obtain for any $\epsilon \ge (\text{Lip}(V) + 1)h$

$$ \sum_{(x,y) \in \Sigma^h} \frac{(p_x - p_y + V(x) - V(y))^2}{h^2} \Lambda_{p+V}(u_x, u_y) h^d \ge S_{\text{CM}}^h(u, g_\epsilon(p)) + 2d|\Omega|\omega(h) \frac{\epsilon}{h}. $$

Choosing $\epsilon = (\text{Lip}(V) + 1)h$ provides that the error term $\sim \omega(h)$ converges to 0 as $h \to 0$. $\square$

**Theorem 5.5 (Compactness).** Let $\rho_0 \in P(\Omega)$. Let $\tilde{\rho}^{h,\tau} = \tilde{u}^{h,\tau} h^d$ be the variational interpolation from Definition 4.1 for the minimizers of $(\text{JKO}_{\text{cm}}^{h,\tau})$ and $\tilde{p}^{h,\tau}$ be the corresponding pressure variable constructed in (5.3).

Then there exists the density-pressure pair $(u, p)$ with $u \in L^1((0, T), L^1(\Omega))$ and $p \in L^2((0, T), H^1(\Omega))$ such that $p(1-u) = 0$ and

$$ \tilde{u}_t^{h,\tau} \to u_t \quad \text{weakly-} * \text{ in } L^\infty(\Omega) \text{ for } t \in (0, T), $$

$$ \tilde{p}^{h,\tau} \to p \quad \text{weakly in } L^2((0, T), L^2(\Omega)), $$

$$ \int_s^t \tilde{p}_r^{h,\tau} dr \to \int_s^t p_r dr \quad \text{strongly in } L^2(\Omega) \text{ for any } (s,t) \subset (0,T), $$

$$ \overline{\nabla}^h g_{\epsilon}(\tilde{p}^{h,\tau}) \to \nabla p \quad \text{weakly in } L^2((0, T), L^2(\Omega)), $$

where the discrete approximation for the gradient $\overline{\nabla}^h g_{\epsilon}(\tilde{p}^{h,\tau})$ is defined as in (4.31).
*Proof. Step 1. Convergence of the density.* By Lemma 4.13, $\tilde{\rho}^{h,\tau} \to \rho$ uniformly for the $W_2$ distance, where $\rho \in C([0,T]; P(\Omega))$. Since $\|\tilde{u}^{h,\tau}\|_{L^\infty} \le 1$, we have $\tilde{u}_t^{h,\tau} \to u_t$ weakly-* in $L^\infty$ for $t \in (0,T)$.

*Step 2. Convergence of the pressure.* From Lemma 5.4, we have

$$ \sup_{h,\tau>0, h/\tau<1} \int_0^T S_{CM}^h(\tilde{u}_t^{h,\tau}, g_{\epsilon}(h)(\tilde{p}_t^{h,\tau})) dt \le C $$

with $C > 0$ independent of $h$ and $\tau$ and $\epsilon(h) = (\text{Lip}(V)+1)h$. By a similar argument as in Lemma 4.17, one can deduce that the previous bound implies

$$ \sup_{h,\tau>0, h/\tau<1} \int_0^T \sum_{(x,y) \in \Sigma^h} (g_{\epsilon}(h)(p_t^{h,\tau}(y)) - g_{\epsilon}(h)(p_t^{h,\tau}(x)))^2 dt \le C, $$

with possibly different constant $C > 0$ but still independent of $h$ and $\tau$. We define

$$ g^{h,\tau} := \sum_{x \in \mathcal{T}^h} g_{\epsilon}(h)(\tilde{p}^{h,\tau}(x)) 1_{Q_h}(x). $$

Following the proof of Lemma 4.16 with $g_{\epsilon}(h)$ instead of $\hat{\epsilon}^h$, one can derive

$$ \sup_{h>0} \int_0^T \|g_t^{h,\tau} - \bar{g}^{h,\tau}\|_{L^2}^2 dt \le C \quad \text{and} \quad \int_0^T \int_{\Omega|\eta|} |g_t^{h,\tau}(z-\eta) - g_t^{h,\tau}(z)|^2 dz dt \le C|\eta| \max(|\eta|, h), $$

where $\bar{g}^{h,\tau} := \int_0^T \int_{\Omega} g_t^{h,\tau}(x) dx dt$. It is clear that $\bar{g}^{h,\tau} \ge 0$. We also claim that $\sup_{h,\tau>0} \bar{g}^{h,\tau} < \infty$.

Suppose that there exists a subsequence such that $\bar{g}^{h,\tau} \to \infty$, then there exists a subsequence such that $g_t^{h,\tau}(x) \to \infty$ pointwise a.e. on $(0,T) \times \Omega$. Consequently, $\tilde{p}_t^{h,\tau}(x) \to \infty$ pointwise a.e. on $(0,T) \times \Omega$ and $\tilde{u}_t^{h,\tau}(x) \to 1$ pointwise a.e. on $\Omega$ for a.e. $t \in (0,T)$, which is a contradiction with $\int_{\Omega} \hat{u}_t^{h,\tau}(x) dx = 1$ given that $|\Omega| > 1$.

Since $\bar{g}^{h,\tau}$ is uniformly bounded, we find that $g^{h,\tau}$ is uniformly bounded in $L^2((0,T), L^2(\Omega))$ and there exists $p \in L^2((0,T), L^2(\Omega))$ such that $g^{h,\tau} \to p$ weakly in $L^2((0,T), L^2(\Omega))$.

Consider $g_{s,t}^{h,\tau} := \int_t^s g_r^{h,\tau} dr$ and $p_{s,t} = \int_t^s p_r dr$ for some fixed $(t,s) \subset (0,T)$. Of course we have $g_{s,t}^{h,\tau} \to p_{s,t}$ weakly in $L^2(\Omega)$. Since $\|g_{s,t}^{h,\tau}(\cdot + \eta) - g_{s,t}^{h,\tau}\|_{L^2(\Omega|\eta|)} \le \frac{C}{t-s}|\eta| \max(|\eta|, h)$ for any $\eta \in \mathbb{R}^d$, the Riesz-Fréchet-Kolmogorov theorem implies $g_{s,t}^{h,\tau} \to p_{s,t}$ strongly in $L^2(\Omega)$. Furthermore, we have $\tilde{p}_{s,t}^{h,\tau} \to p_{s,t}$ strongly in $L^2(\Omega)$, because

$$ (5.8) \qquad \|\tilde{p}_{s,t}^{h,\tau} - p_{s,t}\|_{L^2} \le \|\tilde{p}_{s,t}^{h,\tau} - g_{s,t}^{h,\tau}\|_{L^2} + \|g_{s,t}^{h,\tau} - p_{s,t}\|_{L^2} \le \epsilon(h)|\Omega| + \|g_{s,t}^{h,\tau} - p_{s,t}\|_{L^2} \to 0. $$

*Step 3. Convergence of the density-pressure relation.* The proof of this step is an adaptation of a similar argument used in [39, Theorem 2.4]. Since the density-pressure relation holds for all the discrete pairs with $h > 0$, then

$$ \begin{align*} 0 &= \int_t^s \int_{\Omega} \tilde{p}_r^{h,\tau}(x)(1 - \tilde{u}_r^{h,\tau}(x)) dx dr \\ &= \int_t^s \int_{\Omega} (m_{\epsilon} * \tilde{p}_r^{h,\tau})(x)(1 - \tilde{u}_r^{h,\tau}(x)) dx dr + \int_t^s \int_{\Omega} (\tilde{p}_r^{h,\tau}(x) - (m_{\epsilon} * \tilde{p}_r^{h,\tau})(x))(1 - \tilde{u}_r^{h,\tau}(x)) dx dr \\ &=: I_{\epsilon}^{h,\tau}(t,s) + E_{\epsilon}^{h,\tau}(t,s), \end{align*} $$
where $m_\epsilon$ is a standard mollifier. The second term is bounded by the Hölder inequality

$$ egin{aligned} |E_{\epsilon}^{h, \tau}(t, s)| & \le \left( \int_{t}^{s} \|\tilde{p}_{r}^{h, \tau} - (m_{\epsilon} * \tilde{p}_{r}^{h, \tau})\|_{2}^{2} dr \right)^{1/2} \left( \int_{t}^{s} \|1 - \tilde{u}_{r}^{h, \tau}\|_{2}^{2} dr \right)^{1/2} \\ & \le \sqrt{|\Omega|} \left( \int_{t}^{s} \|\tilde{p}_{r}^{h, \tau} - (m_{\epsilon} * \tilde{p}_{r}^{h, \tau})\|_{2}^{2} dr \right)^{1/2}. \end{aligned} $$

We apply the change of variables $\eta = x - y$ and Jensen's inequality to get

$$ egin{aligned} \int_{t}^{s} \|\tilde{p}_{r}^{h, \tau} - (m_{\epsilon} * \tilde{p}_{r}^{h, \tau})\|_{2}^{2} dr &= \int_{t}^{s} \left| \int_{\Omega} m_{\epsilon}(x-y)(\tilde{p}_{r}^{h, \tau}(x) - \tilde{p}_{r}^{h, \tau}(y)) dy \right|^{2} dx \\ &\le \int_{t}^{s} \int_{\Omega} \int_{\{\eta \in \mathbb{R}^{d} : x-\eta \in \Omega\}} m_{\epsilon}(\eta) |\tilde{p}_{r}^{h, \tau}(x) - \tilde{p}_{r}^{h, \tau}(x-\eta)|^{2} d\eta dx dr \\ &= \int_{\mathbb{R}^{d}} m_{\epsilon}(\eta) \int_{t}^{s} \int_{\Omega_{|\eta|}} |\tilde{p}_{r}^{h, \tau}(x) - \tilde{p}_{r}^{h, \tau}(x-\eta)|^{2} dx dr d\eta \\ &\le \int_{\mathbb{R}^{d}} m_{\epsilon}(\eta) |\eta| \max(|\eta|, h) d\eta \le \epsilon \max(\epsilon, h), \end{aligned} $$

where the second-to-last inequality follows from a similar argument as in Lemma 4.16 applied to $\int_{t}^{s} \tilde{p}_{r}^{h, \tau} dr$. Thus, we obtain

$$ |E_{\epsilon}^{h, \tau}(t, s)| \le \sqrt{|\Omega| \epsilon \max(\epsilon, h)}. $$

For the first term, we have that

$$ egin{aligned} I_{\epsilon}^{h, \tau}(t, s) &= \int_{t}^{s} \int_{\Omega} (m_{\epsilon} * \tilde{p}_{r}^{h, \tau})(x)(1 - \tilde{u}_{t}^{h, \tau}(x)) dx dr + \int_{t}^{s} \int_{\Omega} (m_{\epsilon} * \tilde{p}_{r}^{h, \tau})(x)(\tilde{u}_{t}^{h, \tau}(x) - \tilde{u}_{r}^{h, \tau}(x)) dx dr \\ &=: A_{\epsilon}^{h, \tau}(t, s) + B_{\epsilon}^{h, \tau}(t, s) \end{aligned} $$

By [39, Lemma 3.4], we have the bound

$$ egin{aligned} |B_{\epsilon}^{h, \tau}(t, s)| & \le \int_{t}^{s} \|\nabla(m_{\epsilon} * \tilde{p}_{r}^{h, \tau})\|_{2} W_{2}(\rho_{t}^{h, \tau}, \rho_{r}^{h, \tau}) dr \le C \sqrt{s-t} \int_{t}^{s} \|\nabla(m_{\epsilon} * \tilde{p}_{r}^{h, \tau})\|_{2} dr \\ & \le C \left( \int_{t}^{s} \|\nabla(m_{\epsilon} * \tilde{p}_{r}^{h, \tau})\|_{2}^{2} dr \right)^{1/2}. \end{aligned} $$

We now aim to pass to the limit $A_{\epsilon}^{h, \tau}(t, s) + B_{\epsilon}^{h, \tau}(t, s) + E_{\epsilon}^{h, \tau}(t, s)$ as $h, \tau \to 0$. By (5.8), we get that

$$ \int_{t}^{s} (m_{\epsilon} * \tilde{p}_{r}^{h, \tau}) dr \to \int_{t}^{s} (m_{\epsilon} * p_{r}) dr \quad \text{as } h, \tau \to 0 \text{ strongly in } L^{2}(\Omega). $$

This convergence of convolution together with the weak-* convergence $u_t^{h, \tau} \to u_t$ in $L^\infty(\Omega)$ for a.e. $t \in (0, T)$ imply

$$ \lim_{h, \tau \to 0} A_{\epsilon}^{h, \tau}(t, s) = \int_{\Omega} \int_{t}^{s} (m_{\epsilon} * p_{r})(x) dr (1 - u_{t}(x)) dx := A_{\epsilon}(t, s). $$

Since $\int_{0}^{T} \|\tilde{p}_{t}^{h, \tau}\|_{L^{2}}^{2} dt$ is bounded, uniformly in $h$ and $\tau$, then $\int_{0}^{T} \|\nabla m_{\epsilon} * \tilde{p}_{t}^{h, \tau}\|_{L^{2}}^{2} dt$ is also bounded uniformly in $h$ and $\tau$ (with a constant depending on $\epsilon$). Therefore, the function $\|\nabla m_{\epsilon} * \tilde{p}_{t}^{h, \tau}\|_{L^{2}}^{2} \in L^{1}((0, T))$ converges weakly as $h, \tau \to 0$ to a measure $\mu \in \mathcal{M}^{+}((0, T))$ and

$$ \lim_{h \to 0, \tau \to 0} B_{\epsilon}^{h, \tau}(t, s) \le C_{\epsilon} \sqrt{\mu([t, s])} =: B_{\epsilon}(t, s). $$
For $E_{\epsilon}^{h,\tau}(t, s)$, we simply get $\lim_{h \to 0, \tau \to 0} |E_{\epsilon}^{h,\tau}(t, s)| \le \sqrt{|\Omega|}\epsilon := E_{\epsilon}$.

We now pass $s \to t$. For any Lebesgue point $t \in (0, T)$ of $p$, we get

$$ \lim_{s \to t} A_{\epsilon}(t, s) = \int_{\Omega} (m_{\epsilon} * p_{t})(x)(1 - u_{t}(x)) dx := A_{\epsilon} $$

and

$$ \lim_{s \to t} |B_{\epsilon}(t, s)| \le C_{\epsilon} \lim_{s \to t} \sqrt{\mu([t, s])} = 0 \quad \text{for a.e. } t \in (0, T). $$

Finally, as $\epsilon \to 0$, we obtain

$$ \lim_{\epsilon \to 0} A_{\epsilon} = \int_{\Omega} p_{t}(x)(1 - u_{t}(x)) dx $$

and the error term $E_{\epsilon}$ vanishes. In total, we obtain

$$ \int_{\Omega} p_{t}(x)(1 - u_{t}(x)) dx = 0 \quad \text{for almost all } t \in [0, T]. $$

**Step 4.** *Convergence of the pressure gradient.* A small modification of the proof of Lemma 4.18 yields

$$ \int_{0}^{T} \| \bar{\nabla}^{h} g_{t}^{h, \tau} \|_{L^{2}}^{2} dt = \int_{0}^{T} S_{\text{CM}}^{h}(\tilde{u}_{t}^{h, \tau}, \tilde{p}_{t}^{h, \tau}) dt. $$

Since the right-hand side is uniformly bounded by $\sup_{h>0} F_{\text{CM}}^{h}(\rho_{0}^{h}) < \infty$, there exists $\zeta \in L^{2}((0, T), L^{2}(\Omega))$ such that $\bar{\nabla}^{h} g^{h,\tau} \to \zeta$ weakly in $L^{2}((0, T), L^{2}(\Omega))$. Following similar lines as the proof of Lemma 4.18, we obtain for an arbitrary $\psi \in C_{c}^{\infty}((0, T) \times \Omega; \mathbb{R}^{d})$ that

$$ \left| \int_{0}^{T} \int_{\Omega} (\nabla_{x} \cdot \psi(x, t) g_{t}^{h, \tau}(x) + \psi(x, t) \cdot (\bar{\nabla}^{h} g_{t}^{h, \tau})(x)) dx dt \right| \le C h \left( \int_{0}^{T} S_{\text{CM}}^{h}(\tilde{u}_{t}^{h, \tau}, \tilde{p}_{t}^{h, \tau}) dt \right)^{1/2}. $$

This gives an approximate integration-by-parts formula

$$ \int_{0}^{T} \int_{\Omega} \nabla_{x} \cdot \psi(x, t) g_{t}^{h, \tau}(x) dx dt = \int_{0}^{T} \int_{\Omega} \psi(x, t) \cdot (\bar{\nabla}^{h} g_{t}^{h, \tau})(x) dx dt + O(h) $$

and passing to the limit $h, \tau \to 0$ we obtain

$$ \int_{0}^{T} \int_{\Omega} \nabla_{x} \cdot \psi(x, t) p_{t}(x) dx dt = \int_{0}^{T} \int_{\Omega} \psi(x, t) \cdot \zeta_{t}(x) dx dt. $$

Since $\psi$ is arbitrary, we have $\zeta = \nabla p$.

<figure><img src="image_8.png" /></figure>

**Lemma 5.6.** *Under the assumptions of Theorem 5.5, we get the following liminf inequality for the Fisher information:*

$$ \liminf_{h, \tau \to 0} \int_{0}^{T} S_{\text{CM}}^{h}(\tilde{u}_{t}^{h, \tau}, g_{\epsilon(h)}(\tilde{p}_{t}^{h, \tau})) dt \ge \int_{0}^{T} S_{\text{CM}}(u_{t}, p_{t}) dt. $$

*Proof.* We consider the three part of $S_{\text{CM}}^{h}$ one by one. For the first part, following the lines of Lemma 4.18, one gets that

$$ \begin{aligned} \liminf_{h, \tau \to 0} \int_{0}^{T} \frac{1}{4} \sum_{(x,y) \in \Sigma^{h}} (g_{\epsilon}(\tilde{p}_{x}^{h, \tau}(t)) - g_{\epsilon}(\tilde{p}_{y}^{h, \tau}(t)))^2 h^{d-2} dt &= \liminf_{h, \tau \to 0} \frac{1}{2} \int_{0}^{T} \|\bar{\nabla}^{h} g_{\epsilon}(\tilde{p}_{t}^{h, \tau})\|_{L^2}^2 dt \\ &\ge \frac{1}{2} \int_{0}^{T} \|\nabla p_{t}\|_{L^2}^2 dt. \end{aligned} $$
To handle the cross-term, we use that for all $x \in \mathcal{T}^h$ and $\mathbf{h} \in \mathbf{d}_h$

$$V(x + \mathbf{h}) - V(x) = (\nabla V)(x) \cdot \mathbf{h} + o(h) = \int_{Q_h(x+\mathbf{h}/2)} (\nabla V)(z) \, dz \cdot \mathbf{h} + o(h).$$

Thus, we have

$$\begin{aligned} & \frac{1}{2} \sum_{(x,y) \in \Sigma^h} (g_{\epsilon}(\tilde{p}_x^{h,\tau}) - g_{\epsilon}(\tilde{p}_y^{h,\tau}))(V(x) - V(y))h^{d-2} \\ &= \frac{1}{2} \sum_{x \in \mathcal{T}^h} h^d \sum_{\mathbf{h} \in \mathbf{d}_h} \int_{Q_h(x+\mathbf{h}/2)} (g_{\epsilon}(\tilde{p}_{x+\mathbf{h}}^{h,\tau}) - g_{\epsilon}(\tilde{p}_x^{h,\tau})) \frac{\mathbf{h}}{h^2} \cdot (\nabla V)(z) \, dz + o(1)|_{h \to 0} \\ &= \int_{\Omega} \tilde{\nabla}^h (g_{\epsilon}(\tilde{p}^{h,\tau}))(z) \cdot (\nabla V)(z) \, dz + o(1)|_{h \to 0}. \end{aligned}$$

Since $\tilde{\nabla}^h g_{\epsilon}(\tilde{p}^{h,\tau}) \to \nabla p$ weakly in $L^2((0,T), L^2(\Omega))$, we obtain

$$\lim_{h, \tau \to 0} \int_0^T \frac{1}{2} \sum_{(x,y) \in \Sigma^h} (g_{\epsilon}(\tilde{p}_x^{h,\tau})(t) - g_{\epsilon}(\tilde{p}_y^{h,\tau})(t))(V(x) - V(y))h^{d-2} \, dt = \int_0^T \int_{\Omega} \nabla p_t(z) \cdot \nabla V(z) \, dz \, dt.$$

We now turn to the last term in the Fisher information. The same idea was used in [30, Theorem 6.2], and we briefly repeat it here for completeness. A simple rewriting gives

$$\frac{1}{4} \sum_{(x,y) \in \Sigma^h} (V(x) - V(y))^2 \Lambda_V(u_x^h, u_y^h) h^{d-2} = \frac{1}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} u_x^h ((V(x) - V(x - \mathbf{h}))^+)^2 h^{d-2}.$$

Since the positive part is a Lipschitz continuous function and $V \in C^1$, we have

$$(V(x) - V(x - \mathbf{h}))^+ = (\nabla V(x) \cdot \mathbf{h})^+ + o(h),$$

therefore,

$$\begin{aligned} & \frac{1}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} u_x^h ((V(x) - V(x - \mathbf{h}))^+)^2 h^{d-2} = \frac{1}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} u_x^h ((\nabla V(x))^+ \cdot \mathbf{h})^2 h^{d-2} + o(1)|_{h \to 0} \\ &= \frac{1}{2} \sum_{x \in \mathcal{T}^h} u_x^h \langle \nabla V(x), \sum_{\mathbf{h} \in \mathbf{d}_h} \mathbf{h} \otimes \mathbf{h} \tilde{1}_{\nabla V(x) \cdot \mathbf{h}} \nabla V(x) \rangle h^{d-2} + o(1)|_{h \to 0}, \end{aligned}$$

where we use the notation $\tilde{1}_{\nabla V(x) \cdot \mathbf{h}} = 1\{\nabla V(x) \cdot \mathbf{h} > 0\} + \frac{1}{2}1\{\nabla V(x) \cdot \mathbf{h} = 0\}$. The tensor can be rewritten as

$$h^{d-2} \sum_{\mathbf{h} \in \mathbf{d}_h} \mathbf{h} \otimes \mathbf{h} \tilde{1}_{\nabla V(x) \cdot \mathbf{h}} = h^d \sum_{i=1}^d e_i \otimes e_i = h^d \mathrm{Id}.$$

Therefore,

$$\begin{aligned} & \frac{1}{2} \sum_{x \in \mathcal{T}^h} \sum_{\mathbf{h} \in \mathbf{d}_h} u_x^h ((V(x) - V(x - \mathbf{h}))^+)^2 h^{d-2} = \frac{1}{2} \sum_{x \in \mathcal{T}^h} u_x^h |\nabla V(x)|^2 h^d + o(1)|_{h \to 0} \\ &= \frac{1}{2} \sum_{x \in \mathcal{T}^h} \int_{Q_h(x)} |\nabla V(z)|^2 \hat{u}^h(z) \, dz + o(1)|_{h \to 0}, \end{aligned}$$
with $\hat{u}^h$ being piecewise constant reconstruction on $\{Q_h(x)\}_{x \in \mathcal{T}^h}$.

$$ \liminf_{h, \tau \to 0} \int_0^T \frac{1}{4} \sum_{(x,y) \in \Sigma^h} (V(x) - V(y))^2 \Lambda_V (\tilde{u}_x^{h,\tau}(t), \tilde{u}_y^{h,\tau}(t)) h^{d-2} dt \ge \liminf_{h, \tau \to 0} \int_0^T \int_0^T |\nabla V|^2 \tilde{u}_t^{h,\tau} dz dt \\ = \int_0^T \int_{\Omega} |\nabla V(z)|^2 u_t(z) dz dt. \quad \square $$

We conclude with the proof of Theorem 5.2.

*Proof.* Similarly to Lemma 4.12, we use one-step variational inequality for $(JKO_{cm}^{h,\tau})$ as in (4.7), the lower bound with the Fisher information from Lemma 5.4, geodesic interpolation $(\tilde{\rho}^{h,\tau}, \tilde{\nu}^{h,\tau})$, and time rescaling to arrive at

$$ 0 \ge F_{\text{CM}}^{h}(\rho_{N}^{h,\tau}) - F_{\text{CM}}^{h}(\rho_{0}^{h}) + \frac{1}{2} \int_{0}^{T} \|\tilde{\nu}_{t}^{h,\tau}\|_{L^{2}(\tilde{\rho}_{t}^{h,\tau})} dt + \int_{0}^{T} S_{\text{CM}}^{h}(\tilde{\rho}_{t}^{h,\tau}, g_{\epsilon}(h)(\tilde{\rho}_{t}^{h,\tau})) \lambda_{\epsilon}^{h,\tau}(dt) \\ + \frac{dT}{4} \frac{h}{\tau} \log \epsilon + o(1)_{h \to 0}. \quad \square $$

## REFERENCES

[1] M. Agueh. Existence of solutions to degenerate parabolic equations via the Monge-Kantorovich theory. *Advances in Differential Equations*, 10(3):309–360, 2005.

[2] R. Alicandro and M. Cicalese. A general integral representation result for continuum limits of discrete energies with superlinear growth. *SIAM journal on mathematical analysis*, 36(1):1–37, 2004.

[3] L. Ambrosio, N. Gigli, and G. Savaré. *Gradient flows: in metric spaces and in the space of probability measures*. Springer Science & Business Media, 2008.

[4] G. Auricchio, F. Bassetti, S. Gualandi, and M. Veneroni. Computing kantorovich-wasserstein distances on d-dimensional histograms using (d+1)-partite graphs. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, *Advances in Neural Information Processing Systems*, volume 31. Curran Associates, Inc., 2018.

[5] A. Baradat, A. Hraivoronska, and F. Santambrogio. Using Sinkhorn in the JKO scheme adds linear diffusion, 2025.

[6] J.-D. Benamou and Y. Brenier. A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. *Numerische Mathematik*, 84(3):375–393, 2000.

[7] J.-D. Benamou, G. Carlier, and M. Laborde. An augmented Lagrangian approach to Wasserstein gradient flows and applications. *ESAIM: Proceedings and surveys*, 54:1–17, 2016.

[8] J.-D. Benamou, G. Carlier, Q. Mérigot, and E. Oudet. Discretization of functionals involving the Monge–Ampère operator. *Numerische mathematik*, 134(3):611–636, 2016.

[9] A. Blanchet, V. Calvez, and J. A. Carrillo. Convergence of the mass-transport steepest descent scheme for the subcritical Patlak–Keller–Segel model. *SIAM Journal on Numerical Analysis*, 46(2):691–721, 2008.

[10] N. Bonneel. Lemon network simplex implementation. https://github.com/nbonneel/network_simplex. Accessed: 2025-04-11.

[11] N. Bonneel, M. van de Panne, S. Paris, and W. Heidrich. Displacement Interpolation Using Lagrangian Mass Transport. *ACM Transactions on Graphics* (SIGGRAPH ASIA 2011), 30(6), 2011.

[12] G. Bouchitté, C. Jimenez, and M. Rajesh. A new $L^\infty$ estimate in optimal mass transport. *Proc. Amer. Math. Soc.*, 135(11):3525–3535, 2007.

[13] A. Brancolini and G. Buttazzo. Optimal networks for mass transportation problems. *ESAIM: Control, Optimisation and Calculus of Variations*, 11(1):88–101, 2005.

[14] H. Brezis. *Functional analysis, Sobolev spaces and partial differential equations*. Springer Science & Business Media, 2010.

[15] T. Caillet and F. Santambrogio. Doubly nonlinear diffusive PDEs: new existence results via generalized Wasserstein gradient flows. *SIAM Journal on Mathematical Analysis*, 56(6):7043–7073, 2024.
[16] C. Cances, T. O. Gallouët, and G. Todeschi. A variational finite volume scheme for Wasserstein gradient flows. *Numerische Mathematik*, 146:437–480, 2020.

[17] G. Carlier. *Classical and modern optimization*. World Scientific, 2022.

[18] G. Carlier, V. Duval, G. Peyré, and B. Schmitzer. Convergence of entropic schemes for optimal transport and gradient flows. *SIAM Journal on Mathematical Analysis*, 49(2):1385–1418, 2017.

[19] G. Carlier and C. Poon. On the total variation Wasserstein gradient flow and the TV-JKO scheme. *ESAIM: Control, Optimisation and Calculus of Variations*, 25:42, 2019.

[20] J. A. Carrillo, K. Craig, and F. S. Patacchini. A blob method for diffusion. *Calculus of Variations and Partial Differential Equations*, 58:1–53, 2017.

[21] J. A. Carrillo, K. Craig, L. Wang, and C. Wei. Primal dual methods for Wasserstein gradient flows. *Foundations of Computational Mathematics*, pages 1–55, 2022.

[22] J. A. Carrillo, S. Fagioli, F. Santambrogio, and M. Schmidchen. Splitting schemes and segregation in reaction cross-diffusion systems. *SIAM Journal on Mathematical Analysis*, 50(5):5695–5718, 2018.

[23] M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. *Advances in neural information processing systems*, 26, 2013.

[24] N. David and B. Perthame. Free boundary limit of a tumor growth model with nutrient. *Journal de Mathématiques Pures et Appliquées*, 155:62–82, 2021.

[25] S. Di Marino, B. Maury, and F. Santambrogio. Measure sweeping processes. *Journal of Convex Analysis*, 23(2):567–601, 2016.

[26] R. Ducasse, F. Santambrogio, and H. Yoldaş. A cross-diffusion system obtained via (convex) relaxation in the JKO scheme. *Calculus of Variations and Partial Differential Equations*, 62(1):29, 2023.

[27] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. *Handbook of numerical analysis*, 7:713–1018, 2000.

[28] V. Ferrari and F. Santambrogio. Lipschitz estimates on the jko scheme for the fokker-plack equation on bounded convex domains. *Applied Mathematics Letters*, 112:106806, 2021.

[29] O. Gil and F. Quirós. Convergence of the porous media equation to Hele-Shaw. *Nonlinear Analysis: Theory, Methods & Applications*, 44(8):1111–1131, 2001.

[30] A. Hraivoronska, A. Schlichting, and O. Tse. Variational convergence of the Scharfetter–Gummel scheme to the aggregation-diffusion equation and vanishing diffusion limit. *Numerische Mathematik*, 156(6):2221–2292, 2024.

[31] R. Jordan, D. Kinderlehrer, and F. Otto. The variational formulation of the Fokker–Planck equation. *SIAM journal on mathematical analysis*, 29(1):1–17, 1998.

[32] I. Kim and A. R. Mészáros. On nonlinear cross-diffusion systems: an optimal transport approach. *Calculus of Variations and Partial Differential Equations*, 57(3):79, 2018.

[33] J. Kitagawa, Q. Mérigot, and B. Thibert. Convergence of a Newton algorithm for semi-discrete optimal transport. *Journal of the European Mathematical Society*, 21(9):2603–2651, 2019.

[34] H. Leclerc, Q. Mérigot, F. Santambrogio, and F. Stra. Lagrangian discretization of crowd motion and linear diffusion. *SIAM Journal on Numerical Analysis*, 58(4):2093–2118, 2020.

[35] G. Leoni. *A first course in Sobolev spaces*. American Mathematical Soc., 2017.

[36] W. Li, J. Lu, and L. Wang. Fisher information regularization schemes for Wasserstein gradient flows. *Journal of Computational Physics*, 416:109449, 2020.

[37] J. Maas. Gradient flows of the entropy for finite Markov chains. *Journal of Functional Analysis*, 261(8):2250–2292, 2011.

[38] D. Matthes, R. J. McCann, and G. Savaré. A family of nonlinear fourth order equations of gradient flow type. *Communications in Partial Differential Equations*, 34(11):1352–1397, 2009.

[39] B. Maury, A. Roudneff-Chupin, and F. Santambrogio. A macroscopic crowd motion model of gradient flow type. *Mathematical Models and Methods in Applied Sciences*, 20(10):1787–1821, 2010.

[40] B. Maury, A. Roudneff-Chupin, F. Santambrogio, and J. Venel. Handling congestion in crowd motion modeling. *Networks and Heterogeneous Media*, 6(3):485–519, 2011.

[41] R. J. McCann. A convexity principle for interacting gases. *Advances in mathematics*, 128(1):153–179, 1997.

[42] Q. Mérigot. A multiscale approach to optimal transport. In *Computer Graphics Forum*, volume 30, pages 1583–1592. Wiley Online Library, 2011.

[43] A. Natale and G. Todeschi. TPFA finite volume approximation of Wasserstein gradient flows. In *Finite Volumes for Complex Applications IX–Methods, Theoretical Aspects, Examples: FVCA 9, Bergen, Norway, June 2020 IX*, pages 193–201. Springer, 2020.

[44] F. Santambrogio. Optimal transport for applied mathematicians. *Birkhäuser, NY*, 55:58–63, 2015.