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

HAL
open science

# Discretizing the Fokker-Planck equation with second-order accuracy: a dissipation driven approach

Clément Cancès, Léonard Monsaingeon, Andrea Natale

► To cite this version:

Clément Cancès, Léonard Monsaingeon, Andrea Natale. Discretizing the Fokker-Planck equation with second-order accuracy: a dissipation driven approach. 2024. hal-04720129

HAL Id: hal-04720129

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

Preprint submitted on 3 Oct 2024

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.
# DISCRETIZING THE FOKKER-PLANCK EQUATION WITH SECOND-ORDER ACCURACY: A DISSIPATION DRIVEN APPROACH

CLÉMENT CANCÈS, LÉONARD MONSAINGEON, AND ANDREA NATALE

**ABSTRACT.** We propose a fully discrete finite volume scheme for the standard Fokker-Planck equation. The space discretization relies on the well-known square-root approximation, which falls into the framework of two-point flux approximations. Our time discretization is novel and relies on a tailored nonlinear mid-point rule, designed to accurately capture the dissipative structure of the model. We establish well-posedness for the scheme, positivity of the solutions, as well as a fully discrete energy-dissipation inequality mimicking the continuous one. We then prove the rigorous convergence of the scheme under mildly restrictive conditions on the unstructured grids, which can be easily satisfied in practice. Numerical simulations show that our scheme is second order accurate both in time and space, and that one can solve the discrete nonlinear systems arising at each time step using Newton's method with low computational cost.

## 1. INTRODUCTION

1.1. **Fokker-Planck equation and Wasserstein gradient flows.** Because of their broad interest in physics [2, 10, 36, 56], biology [6, 16, 22] or social sciences [27, 47], Wasserstein gradient flows have been the object of strong interest by the mathematical community in the last decades. A prototypical example of such Wasserstein gradient flows is the Fokker-Planck equation

$$ (1.1a) \qquad \partial_t \rho + \operatorname{div} F = 0, $$

$$ (1.1b) \qquad F + \rho \nabla V + \nabla \rho = 0, $$

set in a space time domain $Q_T = (0, T) \times \Omega$, where $\Omega$ is a convex and bounded open subset of $\mathbb{R}^d$ that we further assume to be polyhedral for meshing purposes, and where $T$ is an arbitrary finite time horizon. The background potential $V \in C^2(\bar{\Omega})$ is always assumed to be given and smooth. The kinematics is complemented by a nontrivial initial condition $\rho^0$ and no-flux boundary conditions

$$ (1.1c) \qquad F \cdot n_{\partial\Omega} = 0 \quad \text{on } (0, T) \times \partial\Omega. $$

We always assume that

$$ (1.1d) \qquad \rho(0, \cdot) = \rho^0 \ge 0 \quad \text{with} \quad \mathcal{H}(\rho^0) < +\infty \quad \text{and} \quad \int_{\Omega} \rho^0 > 0. $$

Here the (negative) entropy $\mathcal{H}(\rho)$ and free energy $\mathcal{E}(\rho)$ are defined as

$$ \mathcal{H}(\rho) := \int_{\Omega} (\rho \log \rho - \rho + 1), \quad \mathcal{E}(\rho) := \mathcal{H}(\rho) + \int_{\Omega} \rho V, $$

<footnote n="2010">2010 Mathematics Subject Classification. 65M12, 65M08, 35K10.</footnote>

<footnote>Key words and phrases. Fokker Planck equation, finite volumes, energy dissipation, second order time and space discretization, convergence.</footnote>
and the associated stationary Gibbs measure reads

$$ \pi := e^{-V}. $$

A formal multiplication of (1.1a) by $\log \frac{\rho}{\pi}$ provides that

$$ (1.2) \qquad \frac{d}{dt} \mathcal{E}(\rho) + \int_{\Omega} \rho \left| \nabla \log \frac{\rho}{\pi} \right|^2 = 0, $$

meaning that the energy is dissipated along time at a precise rate. Since the seminal work of Otto [34, 53] and the book of Ambrosio, Gigli and Savaré [1], the minimizing movement scheme (also referred to as the JKO scheme in this setting) has been playing a central role for the analysis and numerical discretization of gradient flows. It indeed enjoys strong stability properties as well as a variational structure, in the sense that it amounts to a minimization problem at each time step. Being the limiting curve obtained by the convergence of a JKO scheme is even one of the possible characterizations of abstract metric gradient flows, see [1]. However, although very attractive from a theoretical point of view, the JKO scheme suffers from several downsides when it comes to practical implementation. First, it is merely first order accurate in time. Second, the optimality condition of the problem to be solved at each time-step amounts to a continuous in time mean field game, that further needs to be discretized, either with inner time stepping [5, 14] or by linearizing the Wasserstein distance [12, 40, 52]. Moreover, full space discretization of the JKO scheme on fixed grids creates difficulties which do not arise in semi-discrete problems [44, 19, 33, 46], and as a consequence the JKO time discretization is often replaced by the computationally cheaper Backward Euler scheme [13]. Lagrangian particle schemes [38, 28, 51] and moving meshes [45, 15] on the other hand, are generally only first order accurate in space.

Our new approach here will circumvent these numerical bottlenecks and relies instead on another by-now classical characterization of gradient flows (still formal at this stage): Any smooth density $\rho = \rho(t, x)$ solving the continuity equation (1.1a) with no-flux boundary conditions satisfies

$$ (1.3) \qquad \frac{d}{dt} \mathcal{E}(\rho) = \int_{\Omega} \log \frac{\rho}{\pi} \partial_t \rho = \int_{\Omega} F \cdot \nabla \log \frac{\rho}{\pi} \geq -\frac{1}{2} \int_{\Omega} \frac{|F|^2}{\rho} - \frac{1}{2} \int_{\Omega} \rho \left| \nabla \log \frac{\rho}{\pi} \right|^2. $$

The key observation is that equality holds in Young's inequality above if and only if (1.1b) is fulfilled. Therefore if $\rho$ is such that (1.1a) holds and satisfies in addition the reverse Energy-Dissipation Inequality (EDI):

$$ (1.4) \qquad \frac{d}{dt} \mathcal{E}(\rho) \leq -\frac{1}{2} \int_{\Omega} \frac{|F|^2}{\rho} - \frac{1}{2} \int_{\Omega} \rho \left| \nabla \log \frac{\rho}{\pi} \right|^2, $$

then $\rho$ must also satisfy (1.1b). The first term in the right-hand side dissipation only depends on the kinematics through the continuity equation (1.1a), while the second part, known as the Fisher information functional, is related to the specific choice of an energy through the first variation $\frac{\delta \mathcal{E}}{\delta \rho} = \log \frac{\rho}{\pi}$. In order to make the EDI formulation rigorous, and following ideas of [4], one introduces the one-homogeneous, convex, and lower semi-continuous Benamou-Brenier function $B: [0, +\infty) \times \mathbb{R}^d \to [0, +\infty]$ defined by

$$ (1.5) \qquad B(a, b) := \begin{cases} \frac{|b|^2}{2a} & \text{if } a > 0, \\ 0 & \text{if } b = 0 \text{ and } a = 0, \\ +\infty & \text{otherwise,} \end{cases} $$
and observes that the Fisher information rewrites as a convex function of $\rho$ under the form

$$ \mathcal{R}(\rho) := \frac{1}{2} \int_{\Omega} \rho \left| \nabla \log \frac{\rho}{\pi} \right|^2 = 2 \int_{\Omega} \pi \left| \nabla \sqrt{\frac{\rho}{\pi}} \right|^2 . $$

Then, integrating (1.4) over time leads to the following notion of EDI solution to (1.1), thoroughly developed in [1].

**Definition 1.** A curve $\rho \in C([0, T]; L^1(\Omega))$ is an EDI solution to (1.1) corresponding to the initial solution $\rho^0$ if, denoting by $\rho^T = \rho(T, \cdot)$, there holds

$$ (1.6) \qquad \mathcal{E}(\rho^T) + \int_0^T \mathcal{R}(\rho) + \inf_F \left\{ \int_{Q_T} B(\rho, F) \right\} \le \mathcal{E}(\rho^0), $$

where the infimum is taken among vector fields $F \in L^1(Q_T; \mathbb{R}^d)$ satisfying the continuity equation $\partial_t \rho + \text{div} F = 0$ with initial/terminal data $\rho^0, \rho^T$ and no-flux boundary conditions:

$$ (1.7) \qquad \int_{Q_T} \rho \partial_t \varphi + \int_{Q_T} F \cdot \nabla \varphi - \int_{\Omega} \rho^T \varphi(T, \cdot) + \int_{\Omega} \rho^0 \varphi(0, \cdot) = 0, \quad \forall \varphi \in C^1(\bar{Q}_T). $$

In the previous definition $L^1(\Omega)$ denotes the space of integrable functions equipped with its weak topology. It is not difficult to check [1] that densities $\rho \in L^1(Q_T)$ satisfying (1.7) with finite kinetic energy $\iint_{Q_T} B(\rho, F)$ are $L^1(\Omega)$-continuous in time, and satisfy $\rho(0) = \rho^0, \rho(T) = \rho^T$.

For the sake of self-completeness we collect basic properties of EDI solutions in Appendix A. Let us only mention at this stage that a) EDI solutions are unique, and b) they are automatically distributional solutions of the Fokker-Planck equation (1.1).

**1.2. Our contribution and organization of the paper.** Our goal here is to propose a fully discrete finite volume scheme based on a two-point flux approximation (TPFA) which satisfies a discrete counterpart of the EDI formulation (1.4) while being second order accurate in both space and time.

The space discretization we adopt relies on the well-established square-root approximation (SQRA) scheme [41, 32]. As exploited in [11], this space discretization enjoys a dissipative structure involving hyperbolic cosine dissipation potentials, cf. Section 3.2, very much related to [49, 54, 55]. Our main contribution here concerns the time discretization, and our strategy consists in capturing the energy dissipation

$$ (1.8) \qquad \int_{n\tau}^{(n+1)\tau} \mathcal{R}(\rho) + \int_{n\tau}^{(n+1)\tau} \int_{\Omega} B(\rho, F) $$

over each time step in a sufficiently accurate way to be discussed shortly. In [35, 58] this is achieved by recursively minimizing a discrete version of the full dissipation functional appearing in the left-hand side of the EDI (1.6), resulting in a variational – but first order scheme. In this work our approximation of the above quantity will rely instead on a mid-point rule in order to recover second order accuracy. More precisely, given a density $\rho^n$ at time $t^n = n\tau$, we first introduce a density $\theta^{n+1/2}$ at the intermediate time $t^{n+1/2} = (n+1/2)\tau$. The next time step will be obtained as a particular extrapolation

$$ \rho^{n+1} = \Xi(\rho^n, \theta^{n+1/2}) $$

to be detailed in Section 2.2 below. Defining the Fokker-Planck flux at time $t^{n+1/2}$

$$ F^{n+1/2} = \theta^{n+1/2} \nabla \log \frac{\theta^{n+1/2}}{\pi} = \pi \nabla \frac{\theta^{n+1/2}}{\pi} = \nabla \theta^{n+1/2} + \theta^{n+1/2} \nabla V $$
(or rather, its discrete SQRA counterpart (2.10) later on), we further impose the discrete continuity equation

$$ \frac{\rho^{n+1} - \rho^n}{\tau} + \operatorname{div} F^{n+1/2} = 0, $$

and we finally approximate the dissipation (1.8) by

$$ \tau \left\{ \int_{\Omega} B(\theta^{n+1/2}, F^{n+1/2}) + \mathcal{R}(\theta^{n+1/2}) \right\}. $$

In addition to the desirable preservation of mass and positivity, our specific choice of extrapolation combined with the variational structure of the SQRA flux will entail a discrete (upper) chain rule for the energy density (see Lemma 2.1 below), which in turn will crucially result in a fully discrete EDI inequality. Passing next to the limit in an appropriate sense, we will establish full convergence of the scheme towards a dissipative EDI solution. In the time-continuous setting, this idea of proving convergence by passing to the limit in a semi-discrete EDI was already implemented in [23] for a similar finite volume discretization of the Fokker-Planck. This essentially boils down to proving asymptotic lower bounds on the two dissipation functionals involved in a continuous-in-time discrete-in-space EDI, which in our setting will be given by Propositions 4.5 and 4.6 below. Nonetheless, the possibility of vacuum (vanishing of the densities $\rho, \theta$) makes the analysis more delicate in our setting compared to [23], and the methods of proof are different.

**Remark 1.1.** Exploiting the estimates and the resulting compactness properties established later on, one could directly prove that the solutions produced in the limit by our numerical scheme are distributional solutions of the Fokker-Planck equation, following for instance the methodology proposed in [9, 11] and even for non-convex domains $\Omega$. At the price of the convexity assumption on $\Omega$, the recovery of (1.1) in the distributional sense is for free (cf. Proposition A.2 in the Appendix), together with uniqueness of the EDI solutions [30].

Finally, let us stress that our time extrapolation is local and only requires pointwise evaluation ($\rho_K^{n+1} = \Xi(\rho_K^n, \theta_K^{n+1/2})$) for each cell $K$ in the finite volume discretization), in contrast to the nonlocal one introduced in [29]. Similarly, the mid-point rule proposed in [39] relies on the computationally expensive evaluation of intermediate Wasserstein geodesics, which we completely dispense from. Our approach also shares some features with [43], but with different relations between $\theta^{n+1/2}$, $\rho^n$ and $\rho^{n+1}$. The specific choices we make for $\Theta$ and $\Xi$ in this paper allow us to rigorously prove the convergence or our scheme, beyond the partial consistency and stability results provided in [43].

The paper is organized as follows: The scheme is introduced in Section 2. After introducing usual concepts related to TPFA finite volumes in Section 2.1, the space and time discretization are presented in Section 2.2, where the extrapolation is constructed. Then elements of numerical analysis at fixed grid are presented in Section 3. This encompasses the well-posed character of the scheme in Section 3.1 as well as the fully discrete EDI in Section 3.2. The latter plays a crucial role in Section 4, where the convergence of the scheme towards an EDI solution is established under some restriction on the mesh detailed in Section 4.1. Compactness properties on the approximate reconstructions are then derived in Section 4.2 and refined in 4.3 thanks to some discrete Aubin-Lions-Simon argument. We pass to the limit and establish two separate Gamma-liminf's for the dissipation functionals in Sections 4.4 and 4.5, and the full convergence is then detailed in Section 4.6. Numerical results are then presented in Section 5, showing that our scheme is second order accurate in time and space. Finally, we defer two technical parts to the appendix: Appendix A recalls
basic properties of EDI solutions, while Appendix B contains an extension adapted to our needs of an Aubin-Lions-Simon lemma by Moussa [50].

## 2. A SPACE-TIME DISCRETIZATION FOR THE FOKKER-PLANCK EQUATION

### 2.1. Finite volume discretization.
The space discretization of our scheme falls into the framework of TPFA finite volumes. It requires the definition of an admissible mesh of the domain $\Omega \subset \mathbb{R}^d$, which is assumed to be polyhedral with Lebesgue measure $m_\Omega > 0$.

An admissible mesh of $\Omega$ is a triplet $(\mathcal{T}, \bar{\Sigma}, (x_K)_{K \in \mathcal{T}})$, consisting in a set of cells $K \in \mathcal{T}$, facets $\sigma \in \bar{\Sigma}$, and cell centers $x_K \in \Omega$, satisfying in addition the conditions in [20, Definition 9.1]. Specifically, we require the following:

1) The cells $K \in \mathcal{T}$ are open disjoint polyhedra with positive $d$-dimensional Lebesgue measure $m_K > 0$. They form a tessellation of $\Omega$, i.e.

$$ \bigcup_{K \in \mathcal{T}} \bar{K} = \bar{\Omega} \quad \text{and} \quad \sum_{K \in \mathcal{T}} m_K = m_\Omega. $$

2) The facets $\sigma \in \bar{\Sigma}$ are closed subsets of $\bar{\Omega}$ contained in an hyperplane of $\mathbb{R}^d$, and with strictly positive $(d-1)$-dimensional Hausdorff (or Lebesgue) measure denoted by $m_\sigma > 0$. Every facet $\sigma \in \bar{\Sigma}$ satisfies either $\sigma = K|L := \partial K \cup \partial L$ or $\sigma = \partial K \cup \partial \Omega$, for some $K, L \in \mathcal{T}$ with $K \neq L$. The subset of interior facets $\Sigma \subset \bar{\Sigma}$ is the set of facets $\sigma$ for which there exists $K, L \in \mathcal{T}$ such that $\sigma = K|L$.

3) For any cell $K \in \mathcal{T}$, there exists a subset $\bar{\Sigma}_K \subset \bar{\Sigma}$ such that

$$ \partial K = \bigcup_{\sigma \in \bar{\Sigma}_K} \sigma \quad \text{and} \quad \bar{\Sigma} = \bigcup_{K \in \mathcal{T}} \bar{\Sigma}_K. $$

We denote the interior facets associated with a cell $K$ by $\Sigma_K = \bar{\Sigma}_K \cap \Sigma$.

4) Two cell centers $x_K$ and $x_L$ coincide if and only if $K = L$. Moreover, if $\sigma = K|L$ then $x_K - x_L$ is orthogonal to $\sigma$, and denoting $d_\sigma := |x_K - x_L|$, the outward normal to the cell $K$ on the facet $\sigma \in \Sigma_K$ is given by

$$ (2.1) \qquad n_{K\sigma} = \frac{x_L - x_K}{d_\sigma}. $$

Discrete densities are represented by collections of degrees of freedom $\rho = (\rho_K)_{K \in \mathcal{T}} \in \mathbb{R}_+^T$, where $\rho_K$ is the degree of freedom associated with the cell $K \in \mathcal{T}$. Similarly, fluxes are represented by the collection of outward fluxes through the inner facets, and denoted as follows: $F = ((F_{K\sigma}, F_{L\sigma}))_{\sigma=K|L \in \Sigma} \subset \mathbb{R}^{2\Sigma}$. We also define the space of conservative fluxes as follows

$$ \mathbb{F}_{\Sigma} := \{ \boldsymbol{F} \in \mathbb{R}^{2\Sigma} ; F_{K\sigma} + F_{L\sigma} = 0 \quad \forall \sigma = K|L \in \Sigma \}. $$

For any $\boldsymbol{F} \in \mathbb{F}_\Sigma$ we denote $F_\sigma := |F_{K\sigma}| = |F_{L\sigma}|$.

We discretize a fixed time interval $[0, T]$ in $N \in \mathbb{N}^*$ time steps of size $\tau = T/N$. A discrete time-dependent density is described by a collection $(\rho^n)_{n=0}^N$, where $\rho^n$ is the discrete density associated with the time $t^n = n\tau$. Discrete time-dependent fluxes are staggered in time with respect to the densities and they are therefore described by $(F^{n+1/2})_{n=0}^{N-1}$, with $F^{n+1/2}$ representing now the discrete fluxes at time $t^{n+1/2} = (n + 1/2)\tau$.
2.2. **Numerical scheme.** As already mentioned, the dissipation properties of our scheme will only be guaranteed by the correct choice of an extrapolation $\rho^{n+1} = \Xi(\rho^n, \theta^{n+1/2})$. In order to construct $\Xi$ we first define a specific nonlinear mean $\Theta : \mathbb{R}_+ \times \mathbb{R}_+ \to \mathbb{R}_+$

$$ (2.2) \quad \Theta(a, b) := H^{*'} \left( \frac{H(b) - H(a)}{b - a} \right) = \exp \left( \frac{b \log b - a \log a}{b - a} - 1 \right), \quad \forall a, b > 0, a \neq b, $$

where the entropy function

$$ H(a) := \begin{cases} a \log a - a + 1 & \text{if } a > 0 \\ 1 & \text{if } a = 0, \\ +\infty & \text{if } a < 0 \end{cases}, $$

has explicit Legendre-Fenchel transform $H^*(p) = \exp(p) - 1$. We naturally extend by continuity

$$ (2.3) \quad \Theta(a, a) := a, \quad \Theta(a, 0) = \Theta(0, a) := e^{-1}a, \quad \forall a \ge 0. $$

Note that by usual properties of convex duality there holds $(H^*)'(p) = (H')^{-1}(p)$ for all $p \in \mathbb{R}$. This fact together with (2.2)-(2.3) implies that

$$ (2.4) \quad (b - a)H'(\Theta(a, b)) = H(b) - H(a), \quad \forall a, b \ge 0, (a, b) \ne (0, 0), $$

which will precisely entail the discrete chain rule.

At least formally, our extrapolation is simply given by inverting the mean, i.e. $\Xi(a, \cdot) = \Theta(a, \cdot)^{-1}$. However, as is clear from Figure 1, $\Theta(a, 0) > 0$ prevents any global invertibility and some extra care is needed in order to obtain a well-posed scheme. To this end, one can check that $\Theta \in C(\mathbb{R}_+ \times \mathbb{R}_+; \mathbb{R}_+)$ defined by (2.2) and (2.3) is jointly concave in its arguments and 1-homogeneous. In particular, defining the concave, increasing function $f \in C(\mathbb{R}_+; [e^{-1}, \infty))$ as

$$ f(r) := \Theta(1, r) = \begin{cases} \exp \left( \frac{r \log(r) - r + 1}{r - 1} \right) & \text{if } r > 0, \\ e^{-1} & \text{if } r = 0 \end{cases}, $$

we have that for any $a, b > 0$

$$ \Theta(a, b) = af\left(\frac{b}{a}\right) = bf\left(\frac{a}{b}\right) = \Theta(b, a). $$

Since $f$ is concave and increasing (see Figure 1), its inverse $r = f^{-1}(s)$ is unambiguously defined at least on $[e^{-1}, +\infty)$. Extending this inverse to the whole real line $s \in \mathbb{R}$, our extrapolation $\Xi$ is finally defined as

$$ (2.5) \quad g(s) := \begin{cases} f^{-1}(s) & \text{if } s > e^{-1} \\ 0 & \text{otherwise} \end{cases} \quad \text{and} \quad \Xi(a, c) := \begin{cases} ag\left(\frac{c}{a}\right) & \text{if } a > 0, \\ ec & \text{if } a = 0 \end{cases}, \quad c \in \mathbb{R}. $$

Note that $f$ has vertical tangent at $r = 0^+$, which implies that $g(\cdot)$ and $\Xi(a, \cdot)$ are $C^1$, convex functions for any fixed $a \ge 0$ as depicted in Figure 1. For any fixed $a \ge 0$, $\Theta(a, \cdot)$ is an invertible map from $\mathbb{R}^+$ to $[ae^{-1}, \infty)$ and $\Xi(a, \cdot)$ coincides with its inverse when restricted on $[ae^{-1}, \infty)$. In other words, $c = \Theta(a, b)$ is a mean between $a$ and $b$, whereas $b = \Xi(a, c)$ is the corresponding extrapolation. We stress again that, for any $a \ge 0$, $\Xi(a, \cdot)$ is a well-defined $C^1$, convex, non-decreasing function on the whole $\mathbb{R}$. We have moreover

$$ \Xi(a, \Theta(a, b)) = b, \quad \Theta(a, \Xi(a, c)) = c \quad \forall a, b \ge 0, c \ge e^{-1}a, $$

but this invertibility relation may fail if $c < e^{-1}a$. This turns out to be quite delicate because our numerical scheme primarily solves for $\theta^{n+1/2}$, and then extrapolates $\rho^{n+1} =$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. A graphical representation of the functions f, g, Θ(a, ·) and Ξ(a, ·).</figcaption></figure>

Ξ(ρⁿ, θⁿ¹/²): If for some reason 0 ≤ θⁿ¹/² < e⁻¹ρⁿ, which does happen at least from our numerical experiments, then the invertibility relation fails and the chain rule (2.4) does not hold as such. Fortunately, and this is the whole cornerstone of our subsequent analysis, one still has an upper chain-rule (2.6). For convenience we collect here useful properties of Θ, Ξ.

**Lemma 2.1.** *There holds*

$$
(2.6) \quad \forall a \ge 0, c > 0 : \quad b = \Xi(a, c) \Rightarrow (b-a)H'(c) \ge H(b) - H(a)
$$

with equality if $c \ge e^{-1}a$, and moreover

$$
(2.7) \quad \forall a, b \ge 0, c \in \mathbb{R} : \quad \Theta(a, b) \ge \min(a, b) \quad \text{and} \quad \Xi(a, c) \ge 2c - a.
$$

The possible failure of equality in (2.6) is precisely what makes our scheme non-variational
in the sense that ρⁿ+¹ cannot be characterized as the minimizer of some functional (cf.
Remark 3.3). Whenever the scheme produces a value θₖⁿ+¹/² < e⁻¹ρₖⁿ (which eventually
happens at least in our simulations) an entropy release

$$
r_K^{n+1/2} = H(\rho_K^{n+1}) - H(\rho_K^n) - H'(\theta_K^{n+1/2})(\rho_K^{n+1} - \rho_K^n) < 0
$$

occurs in (2.6) for ρ<sub>K</sub><sup>n+1</sup> = Ξ(ρ<sub>K</sub><sup>n</sup>, θ<sub>K</sub><sup>n+1/2</sup>), compared to the expected variational equality.
Note however that our scheme keeps some variational character as it amounts to a mini-
mization problem in θ<sup>n+1/2</sup>/π, cf. the proof of Proposition 3.1 below.

Proof. Let us begin with (2.6) and fix $a \ge 0$. Since $c > 0$ one always has $b = \Xi(a, c) \ge 0$, including if $a = 0$ (in which case $b = \Xi(0, c) = ec$). If $b > 0$ then by definition of $\Xi$ we have $c > e^{-1}a$, thus one can legitimately write $\Theta(a, b) = \Xi(a, \cdot)^{-1}(b) \Rightarrow c = \Theta(a, b)$ and from (2.4) we see that equality holds in (2.6). If $b = 0$ and $a = 0$, (2.6) is trivially satisfied. If now $b = 0$ and $a > 0$ then, again by definition of $\Xi$, we see that necessarily $c \le \Theta(a, 0)$ and therefore by convexity $H'(c) \le H'(\Theta(a, 0))$. Whence

$$
(b - a)H'(c) = (0 - a)H'(c) \geq (0 - a)H'(\Theta(a, 0)) = H(0) - H(a)
$$

as desired, where the last equality follows again from (2.4).

As for (2.7), consider first the case $a \le b$. Then $\Theta(a,b) = af(b/a) \ge af(1) = a$, and thus
by symmetry $\Theta(a,b) \ge \min(a,b)$. The second inequality in (2.7) follows by convexity: for
$a \ge 0$ simply write $\Xi(a,c) = ag(c/a) \ge a[g(1)+g'(1)(c/a-1)] = a[1+2(c/a-1)] = 2c-a$,
and the proof is complete.
□
We are now in position of defining the scheme. Let $V \in C^2(\bar{\Omega})$ be a given potential and $\rho^0 \in L^1(\Omega; \mathbb{R}_+)$ a nonnegative density with finite entropy and positive total mass

$$ \mathcal{H}(\rho^0) = \int_{\Omega} H(\rho^0) < \infty, \quad \rho^0[\Omega] := \int_{\Omega} \rho^0 > 0, $$

where as before $H(a) = a \log a - a + 1$. Denote by $V \in \mathbb{R}^T$, $\pi \in \mathbb{R}^T$ and $\rho^0 \in \mathbb{R}^T$ the discrete functions defined by

$$ (2.8) \qquad V_K := V(x_K), \quad \pi_K := \exp(-V_K), \quad \rho_K^0 := \frac{1}{m_K} \int_K \rho^0, \quad \forall K \in \mathcal{T}. $$

A discrete solution is a pair of discrete curves $(\rho^n)_{n=0}^N$ and $(\theta^{n+1/2})_{n=0}^{N-1}$ satisfying for all $n = 0, \dots, N-1$,

$$ (2.9) \qquad m_K \frac{\rho_K^{n+1} - \rho_K^n}{\tau} + \sum_{\sigma \in \Sigma_K} m_\sigma F_{K\sigma}^{n+1/2} = 0, \quad K \in \mathcal{T}, $$

where $F^{n+1/2} \in \mathbb{F}_\Sigma$ is the square-root approximation (SQRA) finite volume flux [41, 32]

$$ (2.10) \qquad F_{K\sigma}^{n+1/2} = \frac{1}{d_\sigma} \pi_\sigma \left( \frac{\theta_K^{n+1/2}}{\pi_K} - \frac{\theta_L^{n+1/2}}{\pi_L} \right) \quad \text{with} \quad \pi_\sigma = \sqrt{\pi_K \pi_L}, \quad K \in \mathcal{T}, \sigma \in \Sigma_K, $$

constructed on the intermediate densities $\theta^{n+1/2} = (\theta_K^{n+1/2})_{K \in \mathcal{T}}$ at time $t^{n+1/2} = t^n + \tau/2$. To complete the scheme, the discrete density $\rho^{n+1}$ at time $t^{n+1}$ is defined from $\rho^n$ and $\theta^{n+1/2}$ by extrapolation:

$$ (2.11) \qquad \rho_K^{n+1} = \Xi(\rho_K^n, \theta_K^{n+1/2}), \quad K \in \mathcal{T}. $$

We stress again that this can be considered as a problem in the single primary variable $\theta^{n+1/2}$, from which $F^{n+1/2}$, $\rho^{n+1}$ can be explicitly obtained whenever needed.

**Remark 2.2.** Our scheme can be thought of as an extension of the usual Crank-Nicolson scheme, which corresponds to the linear time-extrapolation

$$ \rho_K^{n+1} = \Xi_{CN}(\rho_K^n, \theta_K^{n+1/2}) = 2\theta_K^{n+1/2} - \rho_K^n. $$

This scheme is known to be second-order accurate in time and energy stable for quadratic energies. However, it is neither positivity preserving nor entropy-stable for Boltzmann type energies, and its extension to our entropic framework thus requires the introduction of the nonlinear extrapolation (2.5).

Let us also mention that our approach shares similarities with the so-called discrete variational derivative method [26], at least when the relation $\theta^{n+1/2} = \Theta(\rho^n, \rho^{n+1})$ holds true (i.e. when equality holds in (2.6)). However our choice to use $\theta^{n+1/2}$ as an unknown and then to extrapolate to reconstruct $\rho^{n+1}$ allows us to deal with the degenerate geometry stemming from optimal transportation and to incorporate the positivity constraint in the scheme, while allowing the entropy release leading to the inequality in (2.6). This is a cornerstone to carry out the rigorous convergence analysis presented in this paper. The choice of keeping $\theta^{n+1/2}$ as the main unknown is also key in the implementation strategy, which shows great robustness despite the singularly nonlinear character of the scheme.

### 3. DISCRETE WELL-POSEDNESS AND DISSIPATIVE STRUCTURE

In this section we prove the main properties of the scheme (2.9)-(2.11). We establish existence and uniqueness of solutions, as well as a discrete version of the energy dissipation inequality which will be crucial for the convergence analysis carried out in Section 4.
**3.1. Existence and uniqueness of solutions.** We first establish one-step well-posedness of the scheme, and therefore global existence and uniqueness of the whole discrete curve by immediate recursion.

**Proposition 3.1.** For any $\rho^n \ge 0$ with $\sum_K m_K \rho_K^n > 0$ there exists a unique $\theta^{n+1/2}$ and $\rho^{n+1}$ verifying (2.9)-(2.11). Moreover there holds

$$ \theta_K^{n+1/2} > 0, \quad \rho_K^{n+1} \ge 0, \quad \sum_{K \in \mathcal{T}} m_K \rho_K^{n+1} = \sum_{K \in \mathcal{T}} m_K \rho_K^n, $$

and

$$ (3.1) \qquad \frac{\rho_K^{n+1} + \rho_K^n}{2} \ge \theta_K^{n+1/2} \ge \min(\rho_K^{n+1}, \rho_K^n) $$

for all $K \in \mathcal{T}$.

Note in particular that our scheme is positivity and mass preserving.

*Proof.* Recall that on can view (2.9)-(2.11) as a single equation for $\theta^{n+1/2}$. Changing variables $s_K = \theta_K^{n+1/2}/\pi_K$ for all $K \in \mathcal{T}$, it is easy to see that the former problem is equivalent to finding a critical point of

$$ (3.2) \qquad \mathcal{J}(\mathbf{s}) := \frac{1}{\tau} \sum_{K \in \mathcal{T}} m_K J_K(s_K) + \sum_{\sigma \in \Sigma_{\text{int}}} \frac{m_\sigma \pi_\sigma}{2d_\sigma} |s_K - s_L|^2, \quad \mathbf{s} \in \mathbb{R}^T, $$

where $J_K(\cdot)$ is any primitive of the function $s \mapsto \Xi(\rho_K^n, s\pi_K) - \rho_K^n$. Note that $\mathcal{J}$ is $C^1$ and convex. Hence critical points are necessarily global minima, and by compactness at least one minimum exists. (By definition of $\Xi$ it is not difficult to check that all the $J_K$'s are coercive, regardless of the particular value of $\rho_K^n$)

Let $\mathbf{s}$ be any minimizer and write $\theta^{n+1/2}$, $F^{n+1/2}$, $\rho^{n+1}$ for the corresponding auxiliary variables. Summing (2.9) over $K \in \mathcal{T}$ immediately guarantees mass conservation

$$ (3.3) \qquad \sum_{K \in \mathcal{T}} m_K \rho_K^{n+1} = \sum_{K \in \mathcal{T}} m_K \rho_K^n > 0. $$

This implies that, for any minimizer $\mathbf{s}$, there exists a $K \in \mathcal{T}$ such that $\theta_K^{n+1/2} = s_K \pi_K > e^{-1} \rho_K^n \ge 0$. For if not, then $\rho_K^{n+1} \le 0$ for all $K \in \mathcal{T}$ by definition (2.5) of $\Xi$, which in turn would contradict (3.3). As a consequence for any minimizer $\mathbf{s}$ at least one of the $J_K$'s is strictly convex in a neighborhood of $s_K$. This improved convexity in at least one direction suffices to compensate for the lack of strict convexity of the discrete Dirichlet energy in (3.2), and $\mathcal{J}$ is thus strictly convex in the neighborhood of $\mathbf{s}$. Since $\mathcal{J}$ is also globally convex, this proves existence and uniqueness of $\theta^{n+1/2}$ as claimed.

In order to show that $\theta^{n+1/2} \ge 0$, set $K^* = \operatorname*{argmin}_K \theta_K^{n+1/2}/\pi_K$. Then by (2.9)-(2.10)

$$ (3.4) \qquad \Xi(\rho_{K^*}^n, \theta_{K^*}^{n+1/2}) = \rho_{K^*}^{n+1} = \rho_{K^*}^n + \frac{\tau}{m_{K^*}} \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{d_\sigma} \pi_\sigma \left( \frac{\theta_L^{n+1/2}}{\pi_L} - \frac{\theta_{K^*}^{n+1/2}}{\pi_{K^*}} \right) \ge \rho_{K^*}^n $$

From this we see that if $\rho_{K^*}^n = 0$ then by (2.5) $e\theta_{K^*}^{n+1/2} = \Xi(0, \theta_{K^*}^{n+1/2}) \ge 0$. If now $\rho_{K^*}^n > 0$ and $\theta_{K^*}^{n+1/2} \le 0$ then, again by definition of $\Xi$, we would have that $\Xi(\rho_{K^*}^n, \theta_{K^*}^{n+1/2}) = 0$ (3.4), and this would contradict (3.4) since $\rho_{K^*}^n > 0$. Whence $\theta_{K^*}^{n+1/2} \ge 0$, and therefore $\theta_K^{n+1/2} \ge 0$ for all $K \in \mathcal{T}$.
Let us now implement a strong maximum principle-type argument in order to improve this nonnegativity to strict positivity. Assuming by contradiction that $\theta_{K^*}^{n+1/2} = 0$, we see that $\rho_{K^*}^{n+1} = \Xi(\rho_{K^*}^n, 0) = 0$, and evaluating (2.9) for $K = K^*$ yields

$$0 \geq m_K \frac{0 - \rho_{K^*}^n}{\tau} = \sum_{\sigma \in \Sigma_{K^*}} \frac{m_\sigma \pi_\sigma}{d_\sigma} \left( \frac{\theta_L^{n+1/2}}{\pi_L} - 0 \right) \geq 0.$$

This would imply $\theta_L^{n+1/2} = 0$ for all $L \in \mathcal{T}$ sharing a facet with $K$, thus also $\rho_L^{n+1} = 0$ since $\Xi(\rho_L^n, 0) = 0$ always. Propagating from neighboring cell to neighboring cell we would conclude that $\rho^{n+1} \equiv 0$, which would in turn contradict the mass conservation (3.3). Hence, we have shown that $\theta_K^{n+1/2} > 0$ for all $K \in \mathcal{T}$, and as a consequence $\rho_K^{n+1} = H(\rho_K^n, \theta_K^{n+1}) \geq 0$, yet again by definition (2.5) of $\Xi$.

Let us finally establish the bounds (3.1) for $\theta^{n+1/2}$. Since $\theta_K^{n+1/2} > 0$, clearly the lower bound $\theta_K^{n+1/2} \geq \min(\rho_K^{n+1}, \rho_K^n)$ only needs to be checked when both $\rho_K^{n+1}, \rho_K^n > 0$. However, in this case we are necessarily in the “invertibility regime” $\theta_K^{n+1/2} = \Theta(\rho_K^{n+1}, \rho_K^n)$, and the claim immediately follows from the first bound in (2.7). The second bound in (2.7) also gives $\frac{1}{2}[\rho_K^n + \rho_K^{n+1}] = \frac{1}{2}[\rho_K^n + \Xi(\rho_K^n, \theta_K^{n+1/2})] \geq \frac{1}{2}[\rho_K^n + (2\theta_K^{n+1/2} - \rho_K^n)] = \theta_K^{n+1/2}$ and the proof is complete. $\square$

**3.2. Discrete energy dissipation equality.** For any nonnegative discrete density $\rho \geq 0$ we define the discrete total energy of the system

$$\mathcal{E}_{\mathcal{T}}(\rho) := \mathcal{H}_{\mathcal{T}}(\rho) + \sum_{K \in \mathcal{T}} m_K V_K \rho_K, \quad \text{where} \quad \mathcal{H}_{\mathcal{T}}(\rho) := \sum_K m_K H(\rho_K).$$

In this section we show that the solutions of our scheme satisfy a fully discrete energy dissipation inequality with respect to the discrete energy $\mathcal{H}$. We will strongly rely on the following convex real-valued conjugate functions

$$\psi(z) = 2z \operatorname{arcsinh}(z/2) - 2\sqrt{4 + z^2} + 4, \quad \psi^*(\xi) = 4(\cosh(\xi/2) - 1),$$

which emerge naturally in (electro-)chemistry [25, 48, 8], large deviations of jump processes [49], multi-scale limits of diffusion processes [42, 24], and more [55]. Note in particular that for any $a, b > 0$ we have identity

$$\sqrt{ab}(\psi^*)'(\log a - \log b) = a - b.$$

This allows to recast the SQRA fluxes $F^{n+1/2}$ in (2.10) as

$$(3.5) \qquad F_{K\sigma}^{n+1/2} = \frac{\theta_{\sigma}^{n+1/2}}{d_{\sigma}} (\psi^*)' \left( \log \left( \frac{\theta_{K}^{n+1/2}}{\pi_{K}} \right) - \log \left( \frac{\theta_{L}^{n+1/2}}{\pi_{L}} \right) \right)$$

with

$$\theta_{\sigma}^{n+1/2} := \sqrt{\theta_{K}^{n+1/2} \theta_{L}^{n+1/2}}.$$

Next, observe from the critical upper-chain rule (2.6) and $\rho_K^{n+1} = \Xi(\rho_K^n, \theta_K^{n+1/2})$ that

$$(3.6) \qquad H(\rho_K^{n+1}) - H(\rho_K^n) \leq \log(\theta_K^{n+1/2})(\rho_K^{n+1} - \rho_K^n), \quad \forall K \in \mathcal{T}, n \geq 0.$$

Adding $V_K(\rho_K^{n+1} - \rho_K^n) = -\log \pi_K(\rho_K^{n+1} - \rho_K^n)$ on both sides, multiplying by $m_K$, and denoting for convenience

$$\phi_K^{n+1/2} := \log \left( \frac{\theta_K^{n+1/2}}{\pi_K} \right),$$
we find

$$
\begin{align*}
& (3.7) \qquad \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^{n+1}) - \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^n) = \sum_K m_K [H(\rho_K^{n+1}) - H(\rho_K^n)] - \sum_K m_K \log \pi_K [\rho_K^{n+1} - \rho_K^n] \\
& \leq \sum_K m_K \log \left( \frac{\theta_K^{n+1/2}}{\pi_K} \right) (\rho_K^{n+1} - \rho_K^n) \\
& \stackrel{(3.6)}{=} -\tau \sum_K \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{d_\sigma} \phi_K^{n+1/2} \theta_\sigma^{n+1/2} (\psi^*)' (\phi_K^{n+1/2} - \phi_L^{n+1/2}) \\
& = -\tau \sum_{\sigma \in \Sigma} \frac{m_\sigma \theta_\sigma^{n+1/2}}{d_\sigma} (\phi_K^{n+1/2} - \phi_L^{n+1/2}) (\psi^*)' (\phi_K^{n+1/2} - \phi_L^{n+1/2}),
\end{align*}
$$

where we used the fact that the function $(\psi^*)'$ is odd in the last equality. Let us define, for all $\boldsymbol{\rho} \in \mathbb{R}_+^T$ and $\mathbf{F}, \boldsymbol{\xi} \in \mathbb{F}_{\mathcal{T}}$,

$$
\mathcal{D}_{\psi}(\boldsymbol{\rho}, \mathbf{F}) := \sum_{\sigma \in \Sigma} \frac{m_{\sigma} \rho_{\sigma}}{d_{\sigma}} \psi \left( \frac{d_{\sigma} F_{\sigma}}{\rho_{\sigma}} \right), \quad \mathcal{D}_{\psi}^{*}(\boldsymbol{\rho}, \boldsymbol{\xi}) := \sum_{\sigma \in \Sigma} \frac{m_{\sigma} \rho_{\sigma}}{d_{\sigma}} \psi^{*} (d_{\sigma} \xi_{\sigma}),
$$

where as before $\rho_\sigma = \sqrt{\rho_K \rho_L}$. By definition $\mathcal{D}_\psi^*(\rho, \cdot)$ is nothing but the Legendre transform of $\mathcal{D}_\psi(\rho, \cdot)$ with respect to the pairing

$$
\langle \boldsymbol{\xi}, \mathbf{F} \rangle_{\Sigma} = \sum_{\sigma \in \Sigma} m_{\sigma} d_{\sigma} \xi_{K\sigma} F_{K\sigma}.
$$

We also define

$$
(3.8) \qquad \mathcal{R}_{\psi}(\boldsymbol{\rho}) := \mathcal{D}_{\psi}^{*}(\boldsymbol{\rho}, -\nabla_{\Sigma} \phi),
$$

with $\phi \in \mathbb{R}^T$ and $\nabla_{\Sigma}\phi \in \mathbb{F}_{\Sigma}$ given by

$$
\phi_K := \log \left( \frac{\rho_K}{\pi_K} \right), \quad (\nabla_{\Sigma} \phi)_{K\sigma} := \frac{\phi_L - \phi_K}{d_\sigma}.
$$

With these definitions, the calculations above imply altogether:

**Proposition 3.2.** Any discrete solution satisfies the one-step discrete EDI

$$
(3.9) \qquad \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^{n+1}) + \tau \mathcal{D}_{\psi}(\boldsymbol{\theta}^{n+1/2}, \mathbf{F}^{n+1/2}) + \tau \mathcal{R}_{\psi}(\boldsymbol{\theta}^{n+1/2}) \le \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^n),
$$

and equality holds if $\theta_K^{n+1/2} \ge e^{-1} \rho_K^n$ for all $K \in \mathcal{T}$.

*Proof.* Leveraging the expression (3.5) for the fluxes $F_{K\sigma}$, we obtain from (3.7)

$$
\begin{align*}
\mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^{n+1}) - \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^n) &\leq -\tau \sum_{\sigma \in \Sigma} m_{\sigma} d_{\sigma} \frac{\phi_K^{n+1/2} - \phi_L^{n+1/2}}{d_{\sigma}} F_{K\sigma}^{n+1/2} \\
&= -\tau \langle -\nabla_{\Sigma} \phi^{n+1/2}, \mathbf{F}^{n+1/2} \rangle \\
&= -\tau \mathcal{D}_{\psi} \left( \boldsymbol{\theta}^{n+1/2}, \mathbf{F}^{n+1/2} \right) - \tau \mathcal{D}_{\psi}^* \left( \boldsymbol{\theta}^{n+1/2}, -\nabla_{\Sigma} \phi^{n+1/2} \right).
\end{align*}
$$

In the last equality we simply used the equality case in the $\mathcal{D}_\psi, \mathcal{D}_\psi^*$ Fenchel duality, which stands owing to $F_{K\sigma}^{n+1/2} = \frac{\theta_\sigma^{n+1/2}}{d_\sigma} (\psi^*)' - (\nabla_\Sigma \phi)_{K\sigma}$ in (3.5). $\square$
**Remark 3.3.** By analogy with the continuous setting, and similarly to [35, 58], an alternative scheme could consist in defining recursively $\tilde{\rho}^{n+1}$ as a solution to the following variational problem:

$$ (3.10) \qquad \tilde{\rho}^{n+1} \in \underset{\rho \ge 0}{\arg\min} \inf_{F} \{ \mathcal{E}_{\mathcal{T}}(\rho) + \tau \mathcal{D}_{\psi}(\theta, F) + \tau \mathcal{R}_{\psi}(\theta) \}, $$

in which the continuity equation $m_K \frac{\rho_K - \tilde{\rho}_K^n}{\tau} + \sum_{\sigma \in \Sigma_K} m_\sigma F_{K\sigma} = 0$ is imposed as a constraint and $\theta = \Theta(\tilde{\rho}^n, \rho)$. Note that this problem admits indeed minimizers since, by the same calculations as above, we always have that the function minimized in (3.10) is bounded from below by $\mathcal{E}_{\mathcal{T}}(\tilde{\rho}^n)$, and the set of admissible discrete densities is compact. Note also that, discarding $\mathcal{R}_{\psi}(\theta)$, one is left with a discretized version of the classical JKO scheme.

In general, the solution obtained via (3.10) is different from the solution $(\rho^n)_n$ obtained using our scheme. In fact, due to (2.11), we may have $\theta^{n+1/2} \neq \Theta(\rho^n, \rho^{n+1})$ if $\rho^{n+1}$ is not strictly positive. On the other hand, if and whenever our scheme outputs $\theta^{n+1/2} \ge e^{-1}\rho^n$, then the invertibility $\rho^{n+1} = \Xi(\rho^n, \theta^{n+1/2}) \Leftrightarrow \theta^{n+1/2} = \Theta(\rho^n, \rho^{n+1})$ holds and the equality holds in (3.9). As a consequence $\rho^{n+1}$ solves (3.10) with $\tilde{\rho}^n = \rho^n$, since it realizes the lower bound $\mathcal{E}_{\mathcal{T}}(\rho^n)$. Our scheme is somehow “almost variational”, in the sense that it is locally variational except in those situations when entropy releases occur due to equality failure in (2.6). The significant advantage of using our scheme is that the optimality conditions for (3.10) are much harder to manage than the system (2.9)-(2.11), and both the theoretical analysis and numerical implementation for (3.10) become more intricate.

Starting from the expression of $\psi^*$, easy algebra allows to recast the discrete Fisher functional (3.8) as

$$ (3.11) \qquad \mathcal{R}_{\psi}(\rho) = 2 \sum_{\sigma \in \Sigma} \frac{m_{\sigma} \pi_{\sigma}}{d_{\sigma}} \left| \sqrt{\frac{\rho_K}{\pi_K}} - \sqrt{\frac{\rho_L}{\pi_L}} \right|^2. $$

Clearly this is a consistent approximation of the dissipation rate

$$ (3.12) \qquad 2 \int_{\Omega} \pi \left| \nabla \sqrt{\frac{\rho}{\pi}} \right|^2 = \frac{1}{2} \int_{\Omega} \rho \left| \nabla \log \frac{\rho}{\pi} \right|^2 $$

appearing in (1.2).

In order to gain compactness in the next section we exploit Proposition 3.2 to retrieve uniform bounds for the discrete curves $(\rho^n)_n$, $(F^{n+1/2})_n$, $(\theta^{n+1/2})_n$.

**Lemma 3.4.** There exists a constant $C > 0$ only depending on $m_\Omega$, $\mathcal{H}(\rho^0)$, $(\max(V) - \min(V))$, and the total mass $\rho^0[\Omega]$, such that

$$ \sum_{n=0}^{N-1} \tau \left[ \mathcal{D}_{\psi}(\theta^{n+1/2}, F^{n+1/2}) + \mathcal{R}_{\psi}(\theta^{n+1/2}) \right] \le C $$

and

$$ \sup_{0 \le n \le N-1} \left( \mathcal{H}_{\mathcal{T}}(\rho^{n+1}) + \mathcal{H}_{\mathcal{T}}(\theta^{n+1/2}) \right) \le C. $$
Proof. Summing Proposition 3.2 over $n$ we get

$$ 
\begin{aligned}
& \sum_{n=0}^{N-1} \tau \left[ \mathcal{D}_{\psi}(\boldsymbol{\theta}^{n+1/2}, \boldsymbol{F}^{n+1/2}) + \mathcal{R}_{\psi}(\boldsymbol{\theta}^{n+1/2}) \right]^{(3.9)} \le \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^0) - \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^N) \\
& = \left[ \mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^0) + \sum_K m_K V_K \rho_K^0 \right] - \left[ \mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^N) + \sum_K m_K V_K \rho_K^N \right] \\
& \le \mathcal{H}(\rho^0) + (\max(V) - \min(V)) \int_{\Omega} \rho^0,
\end{aligned}
$$

where in the last inequality we used successively Jensen's inequality to bound $\mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^0) \le \mathcal{H}(\rho^0)$, $\mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^N) \ge 0$, and the mass conservation $\sum_K m_K \rho_K^N = \sum_K m_K \rho_K^0 = \int_{\Omega} \rho^0$. For the bound on $\mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^{n+1})$, note first that $\mathcal{D}_{\psi}, \mathcal{R}_{\psi} \ge 0$ in (3.9) and therefore

$$ \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^{n+1}) \le \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^n) \le \dots \le \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^0). $$

This gives similarly

$$ 
\begin{aligned}
\mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^{n+1}) &= \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^{n+1}) - \sum m_K \rho_K^{n+1} V_K \\
&\le \mathcal{E}_{\mathcal{T}}(\boldsymbol{\rho}^0) - \sum m_K \rho_K^{n+1} V_K \\
&= \mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^0) + \sum m_K \rho_K^0 V_K - \sum m_K \rho_K^{n+1} V_K \\
&\le \mathcal{H}(\rho^0) + (\max V - \min V) \rho^0[\Omega].
\end{aligned}
$$

As for the bound on $\mathcal{H}_{\mathcal{T}}(\boldsymbol{\theta}^{n+1/2})$, let us first recall the elementary but useful property of the entropy function

$$ (3.13) \qquad H(c) \le 1 + \frac{1}{2}[H(a) + H(b)] \quad \text{for } a, b, c \ge 0, \quad c \le \frac{a+b}{2} $$

(For $c \le 1$ one has trivially $H(c) \le H(0) = 1$, while for $c \ge 1$ one can simply use the monotonicity $H(c) \le H(a + b/2)$ and conclude by convexity.) Owing to $\theta_K^{n+1/2} \le \frac{1}{2}(\rho_K^n + \rho_K^{n+1})$ from Proposition 3.1, we get

$$ 
\begin{aligned}
\mathcal{H}_{\mathcal{T}}(\boldsymbol{\theta}^{n+1/2}) &= \sum_K m_K H(\theta_K^{n+1/2}) \\
&\le \sum_K m_K \left[ 1 + \frac{1}{2} \left( H(\rho_K^n) + H(\rho_K^{n+1}) \right) \right] \\
&= m_{\Omega} + \frac{1}{2} \left[ \mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^n) + \mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^{n+1}) \right],
\end{aligned}
$$

and the previous uniform bound on $\mathcal{H}_{\mathcal{T}}(\boldsymbol{\rho}^n)$ concludes the proof. $\square$

## 4. CONVERGENCE VIA THE ENERGY DISSIPATION EQUALITY

In this section we establish the convergence of the discrete solutions associated with a sequence of meshes $\mathcal{T}$ and time steps $\tau$, in the limit $\tau, \text{size}(\mathcal{T}) \to 0$, where

$$ \text{size}(\mathcal{T}) := \max\{\text{diam}(K); K \in \mathcal{T}\}. $$

For technical reasons that will appear later in the proofs, we require the sequence of meshes to satisfy some asymptotic isotropy condition inspired from [31], up to some subset of vanishing $d$-dimensional Lebesgue measure as in [18]. We further have to assume some
CFL-type condition, cf. (4.5) in what follows, which for quasi-uniform meshes would simply write $\tau = o(\text{size}(\mathcal{T}))$.

Throughout the section we will denote by $Q_T := [0, T] \times \Omega$ the space-time domain.

4.1. **Assumptions on the sequence of meshes.** Our convergence result relies on the following assumptions on the sequence of meshes:

1) mesh regularity: there exists a constant $\zeta > 0$ uniform w.r.t. $\mathcal{T}$ such that

$$ (4.1) \quad \sum_{\sigma \in \Sigma_K} \frac{m_\sigma d_\sigma}{2d} \le \zeta m_K \quad \text{and} \quad \zeta^{-1} \text{dist}(x_K, K) \le \text{diam}(K) \le \zeta \min_{\sigma \in \Sigma_K} d_\sigma, \quad \forall K \in \mathcal{T} $$

whereas

$$ (4.2) \qquad d_\sigma \le \zeta \text{ size}(\mathcal{T}), \quad \forall \sigma \in \bar{\Sigma}. $$

2) asymptotic isotropy: there exists a subset $\mathcal{T}_{\text{iso}} \subset \mathcal{T}$ and a nonnegative $\varepsilon_{\mathcal{T}} \to 0$ as size($\mathcal{T}$) $\to 0$, such that

$$ (4.3) \quad (1 - \varepsilon_{\mathcal{T}}) |v|^2 \le \frac{1}{2m_K} \sum_{\sigma \in \bar{\Sigma}_K} m_\sigma d_\sigma (v \cdot n_{K\sigma})^2 \le (1 + \varepsilon_{\mathcal{T}}) |v|^2, \quad \forall K \in \mathcal{T}_{\text{iso}}, \forall v \in \mathbb{R}^d; $$

moreover, denoting $\Omega_{\text{iso}} := \cup\{K; K \in \mathcal{T}_{\text{iso}}\}$ we have

$$ (4.4) \qquad \lim_{\text{size}(\mathcal{T}) \to 0} \text{meas}(\Omega \setminus \Omega_{\text{iso}}) = 0; $$

3) CFL-type condition: denoting by $d_{\Sigma}^{\min} = \min_{\sigma \in \Sigma} d_\sigma$, then we assume that

$$ (4.5) \qquad \frac{\tau}{d_{\Sigma}^{\min}} \to 0 \quad \text{as } \text{size}(\mathcal{T}) \to 0. $$

Conditions (4.1) and (4.2) are satisfied by usual discretizations based on Delaunay triangulations (or dual Voronoi diagrams) under mild regularity assumptions. Condition (4.3) is much more restrictive. A weighted version of condition (4.3) was introduced in [31] under the name of asymptotic isotropy to study convergence of discrete optimal transport models to their continuous counterparts. In order to ensure convergence, such weights need to be chosen consistently with the reconstruction operator mapping densities from cells to edges. In our case, the reconstruction is defined by the map $(\theta_K, \theta_L) \mapsto \theta_\sigma = \sqrt{\theta_K \theta_L}$, and for this specific choice the isotropy assumption in [31] takes precisely the form (4.3). This condition also imposes a strong regularity requirement on the meshes. In particular, taking $v = e_i$ with $\{e_i\}_{i=1}^d$ an orthonormal basis and summing over all $i = 1, \dots, d$, this implies

$$ \frac{1}{2m_K} \sum_{\sigma \in \bar{\Sigma}_K} m_\sigma d_\sigma \le d(1 + \varepsilon_{\mathcal{T}}), \quad \forall K \in \mathcal{T}_{\text{iso}}. $$

This is verified if, at least in the limit size($\mathcal{T}$) $\to 0$, each edge $\sigma$ divides the corresponding diamond subcell in two parts of equal area $m_\sigma d_\sigma/2d$. However, in contrast with [31], we allow the isotropy condition to fail in an asymptotically negligible volume $\Omega_{\text{iso}}$, which is precisely the meaning of (4.4). This improved flexibility allows us to consider a practical refinement strategy and generate a sequence of meshes for which the assumption is verified; see Remark 4.1 and [18]. As noted already in [31], condition (4.3) can be obtained by requiring a stronger condition, which is usually referred to as superadmissibility [21] or center of mass condition. Specifically, denoting by $x_\sigma$ the barycenter of the facet $\sigma$, suppose that

$$ (4.6) \qquad x_\sigma = \frac{x_K + x_L}{2}, \quad \forall \sigma = (K|L) \in \Sigma_K $$
Then applying Gauss's theorem to the vector fields $\langle x - x_K, e_i \rangle e_j$ for $i, j = 1, \dots, d$, we recover

$$ \frac{2}{m_K} \sum_{\sigma \in \bar{\Sigma}_K} m_\sigma d_\sigma(x_\sigma - x_K) \otimes (x_\sigma - x_K) = \mathrm{Id}, $$

which directly implies (4.3) on the cell $K$. This suggests that the classical refinement procedure by subsequent subdivisions, described below in Remark 4.1, generates a sequence of meshes for which the assumption holds.

Finally, condition (4.5) is introduced for purely technical reasons in order to guarantee that the reconstructions based on $(\rho^{n+1})_{n \ge 0}$ and $(\theta^{n+1/2})_{n \ge 0}$, defined in (4.7) below, share their cluster points as size($T$) $\to 0$.

**Remark 4.1 (Refinement by subdivision).** Given a bounded polygonal set $\Omega \subset \mathbb{R}^2$, consider an admissible mesh $T^0$ made of acute triangles, and subdivide each control volume by partitioning its edges using a fixed number of points and joining the corresponding points on all edges. Choosing as cell centers the triangles' circumcenters, the superadmissibility condition (4.6) holds for all triangles not sharing an edge with the initial partition $T^0$. Consequently, increasing the number of subdivisions yields a sequence of admissible meshes verifying the asymptotic isotropy assumption above.

**4.2. Compactness and limit densities.** Let us define a reconstruction for the discrete densities and fluxes. For the densities we define, for $0 \le n \le N - 1$,

$$ \begin{align} \rho_{T,\tau}(t,x) &= \rho_K^{n+1} \quad \text{for a.e. } x \in K, t \in (t^n, t^{n+1}], \tag{4.7} \\ \theta_{T,\tau}(t,x) &= \theta_K^{n+1/2} \quad \text{for a.e. } x \in K, t \in (t^n, t^{n+1}], \nonumber \\ \theta_{\Sigma,\tau}(t,x) &= \theta_\sigma^{n+1/2} \quad \text{for a.e. } x \in \Delta_\sigma, t \in (t^n, t^{n+1}], \nonumber \end{align} $$

where

$$ \theta_{\sigma}^{n+1/2} = \begin{cases} \sqrt{\theta_{K}^{n+1/2} \theta_{L}^{n+1/2}} & \text{if } \sigma = K|L \\ \theta_{K}^{n+1/2} & \text{if } \sigma = K \cap \partial\Omega \end{cases}. $$

The diamond cell $\Delta_\sigma$ corresponding to the edge $\sigma$ is a polytope included in $\Omega$, the vertices of which being $x_K$ and those of $\sigma$ if $\sigma \subset \partial\Omega$, and additionally $x_L$ if $\sigma = K|L \in \Sigma$. Note that we do not require $\Delta_\sigma$ to be convex as $x_K$ can lie outside of $K$.

For the initial density profile $\rho^0$, which has been discretized into $\rho^0$ by (2.8), we build the approximation $\rho_T^0$ defined by

$$ \rho_T^0(x) = \rho_K^0 \quad \text{for a.e. } x \in K. $$

Then one readily checks that $\rho_T^0$ converges (strongly) in $L^1(\Omega)$ towards $\rho^0$. We will also need a reconstruction for the terminal discrete density at $t = T$, which will be given by

$$ (4.8) \qquad \rho_{T,\tau}^T(x) = \rho_K^N \quad \text{for a.e. } x \in K. $$

Finally, for the fluxes we use the following reconstruction:

$$ (4.9) \qquad F_{\Sigma,\tau}(t,x) = dF_{K\sigma}^{n+1/2} n_{K\sigma} \quad \text{for a.e. } x \in \Delta_\sigma, t \in (t^n, t^{n+1}], $$

where $F_{K\sigma}^{n+1/2}$ is defined in (3.5) for $\sigma = K|L$, and $F_{K\sigma}^{n+1/2} = 0$ on the boundary $\sigma \subset \partial\Omega$. Note that this is well-defined since $F_{K\sigma}^{n+1/2} = -F_{K\sigma}^{n+1/2}$ and $n_{K\sigma} = -n_{L\sigma}$ for $\sigma = K|L$.
By Lemma 3.4, the total space-time entropies of $\rho_{\mathcal{T},\tau}$ and $\theta_{\mathcal{T},\tau}$ are uniformly bounded, i.e. there exists a constant $C > 0$ independent of $\mathcal{T}, \tau$ such that

$$ (4.10) \qquad \int_{Q_T} H(\rho_{\mathcal{T},\tau}) \le CT \quad \text{and} \quad \int_{Q_T} H(\theta_{\mathcal{T},\tau}) \le CT. $$

Therefore, given any family of admissible meshes $\mathcal{T}^k$ and time steps $\tau^k$ with size($\mathcal{T}^k$) $\to 0$ and $\tau^k \to 0$ as $k \to \infty$, there exists $\rho, \theta \in L^1(Q_T)$ such that, up to extraction of a subsequence if needed and as $k \to \infty$,

$$ (4.11) \qquad \theta_{\mathcal{T}^k, \tau^k} \to \theta, \quad \rho_{\mathcal{T}^k, \tau^k} \to \rho, \quad \text{weakly in } L^1(Q_T). $$

Similarly, since the entropy of $\rho_{\mathcal{T}^k, \tau^k}^T$ is uniformly bounded, we have that there exists $\rho^T \in L^1(\Omega)$ such that, up to a further extraction as $k \to \infty$,

$$ (4.12) \qquad \rho_{\mathcal{T}^k, \tau^k}^T \to \rho^T, \quad \text{weakly in } L^1(\Omega). $$

We claim now that the entropy of $\theta_{\Sigma, \tau}$ is also uniformly bounded. Indeed, as

$$ \theta_{\sigma}^{n+1/2} = \sqrt{\theta_{K}^{n+1/2} \theta_{L}^{n+1/2}} \le \frac{1}{2} (\theta_{K}^{n+1/2} + \theta_{L}^{n+1/2}), $$

it follows from our previous entropy bound (4.10) that

$$ \begin{align} \int_{Q_T} H(\theta_{\Sigma, \tau}) &= \sum_{n=0}^{N-1} \tau \sum_{\sigma \in \Sigma} \mathrm{meas}(\Delta_{\sigma}) H(\theta_{\sigma}^{n+1/2}) \nonumber \\ &= \sum_{n=0}^{N-1} \tau \sum_{\sigma \in \Sigma} \frac{m_{\sigma} d_{\sigma}}{d} H(\theta_{\sigma}^{n+1/2}) \nonumber \\ &\le \sum_{n=0}^{N-1} \tau \sum_{\sigma \in \Sigma} \frac{m_{\sigma} d_{\sigma}}{d} \left[ 1 + \frac{1}{2} \left( H(\theta_{K}^{n+1/2}) + H(\theta_{L}^{n+1/2}) \right) \right] \tag{4.13} \\ &= \sum_{n=0}^{N-1} \tau \sum_{K} \sum_{\sigma \in \Sigma_K} \frac{m_{\sigma} d_{\sigma}}{2d} \left[ 1 + H(\theta_{K}^{n+1/2}) \right] \nonumber \\ &\le \zeta \sum_{n=0}^{N-1} \tau \sum_{K} m_K \left[ 1 + H(\theta_{K}^{n+1/2}) \right] \nonumber \\ &= \zeta \left( m_{\Omega} T + \int_{Q_T} H(\theta_{\mathcal{T}, \tau}) \right) \le CT. \nonumber \end{align} $$

This gives equiintegrability of $\{\theta_{\Sigma, \tau}\}_{\mathcal{T}, \tau}$ for any family of admissible meshes and time steps. We use this to show that the fluxes $\{F_{\Sigma, \tau}\}_{\mathcal{T}, \tau}$ are also equiintegrable. To this end, denote by $d_{\Sigma} \in L^{\infty}(\Omega)$ the piecewise constant function equal to $d_{\sigma}$ in each diamond subcell, and define

$$ (4.14) \qquad D_{\psi} := \sum_{n=0}^{N-1} \tau \mathcal{D}_{\psi}(\theta^{n+1/2}, \mathbf{F}^{n+1/2}) = \int_{Q_T} \frac{\theta_{\Sigma, \tau}}{d_{\Sigma}^2} \psi \left( \frac{|F_{\Sigma, \tau}| d_{\Sigma}}{d \theta_{\Sigma, \tau}} \right). $$

Observe that $D_{\psi}$ is uniformly bounded due to Lemma 3.4. Let now $A \subset Q_T$ be an arbitrary measurable subset, and for any $\varepsilon > 0$ write

$$ \|F_{\Sigma, \tau}\|_{L^1(A)} = \int_A |F_{\Sigma, \tau}| = \frac{d\varepsilon}{2D_{\psi}} \int_A \frac{\theta_{\Sigma, \tau}}{d_{\Sigma}^2} \left( \frac{|F_{\Sigma, \tau}| d_{\Sigma}}{d \theta_{\Sigma, \tau}} \right) \left( \frac{2d_{\Sigma} D_{\psi}}{\varepsilon} \right). $$
By $\psi, \psi^*$-Young's inequality and the expression for $D_\psi$ in (4.14), we obtain

$$ \begin{aligned} \|F_{\Sigma, \tau}\|_{L^1(A)} &\le \frac{d\epsilon}{2D_{\psi}} \int_A \frac{\theta_{\Sigma, \tau}}{d_{\Sigma}^2} \left[ \psi \left( \frac{|F_{\Sigma, \tau}|d_{\Sigma}}{d\theta_{\Sigma, \tau}} \right) + \psi^* \left( \frac{2d_{\Sigma}D_{\psi}}{\epsilon} \right) \right] \\ &\le \frac{\epsilon}{2} + \frac{d\epsilon}{2D_{\psi}} \int_A \frac{\theta_{\Sigma, \tau}}{d_{\Sigma}^2} \psi^* \left( \frac{2d_{\Sigma}D_{\psi}}{\epsilon} \right) \\ &\le \frac{\epsilon}{2} + \frac{d\epsilon}{2D_{\psi}} \max_{\sigma} \left( \frac{1}{d_{\sigma}^2} \psi^* \left( \frac{2D_{\psi}d_{\sigma}}{\epsilon} \right) \right) \int_A \theta_{\Sigma, \tau} \end{aligned} $$

by definition of $d_\Sigma$. For all $\epsilon \ge 0$, there exists a constant $C_\epsilon$ such that $\psi^*(\xi/\epsilon) = 4(\cosh(\xi/2\epsilon) - 1) \le C_\epsilon \xi^2/2$ if $\xi \le 1$. Since $D_\psi$ is bounded we have that eventually $\xi_\sigma = 2D_\psi d_\sigma \le 1$ in the $\max_\sigma(\dots)$ term if $\text{size}(\mathcal{T})$ is sufficiently small (using (4.2)), hence

$$ (4.15) \qquad \|F_{\Sigma, \tau}\|_{L^1(A)} \le \frac{\epsilon}{2} + d\epsilon C_\epsilon D_\psi \int_A \theta_{\Sigma, \tau}. $$

Since $\{\theta_{\Sigma, \tau}\}_{\mathcal{T}, \tau}$ is equiintegrable one can pick $\delta = \delta(\epsilon) > 0$ such that $\text{meas}(A) \le \delta \implies \|\theta_{\Sigma, \tau}\|_{L^1(A)} \le \frac{1}{2dC_\epsilon D_\psi}$, and therefore

$$ \text{meas}(A) \le \delta \implies \|F_{\Sigma, \tau}\|_{L^1(A)} \le \frac{\epsilon}{2} + d\epsilon C_\epsilon D_\psi \|\theta_{\Sigma, \tau}\|_{L^1(A)} \le \epsilon. $$

This means precisely that $\{F_{\Sigma, \tau}\}_{\mathcal{T}, \tau}$ is equiintegrable, and as a consequence we can assume up to extraction of a further subsequence that

$$ (4.16) \qquad F_{\Sigma^k, \tau^k} \rightharpoonup F \quad \text{weakly in } L^1(Q_T; \mathbb{R}^d) $$

for some vector field $F \in L^1(Q_T; \mathbb{R}^d)$.

The next lemma shows that the previous weak limits $\rho, \theta$ from (4.11) coincide, and as of now one should keep in mind $\theta = \rho$. Note carefully that this requires a condition $\tau = o(d_{\Sigma}^{\min})$ on the mesh.

**Lemma 4.2.** Consider a sequence of solutions associated with $(\mathcal{T}^k, \tau^k)_k$ satisfying size$(\mathcal{T}^k) \to 0$ and (4.5), i.e. $\tau^k = o(d_{\Sigma}^{\min})$ as $k \to \infty$. Then

$$ \lim_{k \to \infty} \|\theta_{\mathcal{T}^k, \tau^k} - \rho_{\mathcal{T}^k, \tau^k}\|_{L^1(Q_T)} = 0. $$

*Proof.* Let us consider the solution $(\rho^n)_{n=0}^N, (\theta^{n+1/2})_{n=0}^{N-1}$ obtained for fixed $\mathcal{T}, \tau$. By Proposition 3.1 and (2.9)-(2.10) we control first by $\psi, \psi^*$-Young inequality

$$ \begin{aligned} \sum_K m_K |\theta_K^{n+1/2} - \rho_K^{n+1}|^{(3.1)} &\le \sum_{K \in \mathcal{T}} m_K |\rho_K^{n+1} - \rho_K^n| \\ &\le 2\tau \sum_{\sigma \in \Sigma_{\text{int}}} \frac{m_\sigma}{d_\sigma} \pi_\sigma \left| \frac{\theta_K^{n+1/2}}{\pi_K} - \frac{\theta_L^{n+1/2}}{\pi_L} \right| \\ &\le 2\tau \sum_{\sigma \in \Sigma_{\text{int}}} \frac{m_\sigma}{d_\sigma^2} \theta_\sigma^{n+1/2} \left[ \psi \left( \frac{\pi_\sigma}{\theta_\sigma^{n+1/2}} \left| \frac{\theta_K^{n+1/2}}{\pi_K} - \frac{\theta_L^{n+1/2}}{\pi_L} \right| \right) + \psi^*(d_\sigma) \right] \\ &\le \frac{2\tau}{d_{\Sigma}^{\min}} \sum_{\sigma \in \Sigma} \frac{m_\sigma}{d_\sigma} \theta_\sigma^{n+1/2} \left[ \psi \left( \frac{\pi_\sigma}{\theta_\sigma^{n+1/2}} \left| \frac{\theta_K^{n+1/2}}{\pi_K} - \frac{\theta_L^{n+1/2}}{\pi_L} \right| \right) + \psi^*(d_\sigma) \right]. \end{aligned} $$
Pick now $C > 0$ such that $\psi^*(\xi) = 4(\cosh(\xi/2) - 1) \le C\xi^2$ for $\xi \le 1$. Summing the above estimate over $n$, and recalling the definition (4.14) of $D_\psi$, we obtain

$$ \begin{aligned} \|\theta_{\mathcal{T},\tau} - \rho_{\mathcal{T},\tau}\|_{L^1(Q_T)} &= \sum_{n=0}^{N-1} \tau \sum_{K \in \mathcal{T}} m_K |\theta_K^{n+1/2} - \rho_K^{n+1}| \\ &\le \frac{2\tau}{d_{\Sigma}^{\min}} \left( D_\psi + C \sum_n \tau \sum_{\sigma \in \Sigma} m_\sigma d_\sigma \theta_\sigma^{n+1/2} \right) = \frac{2\tau}{d_{\Sigma}^{\min}} (D_\psi + Cd \|\theta_{\Sigma,\tau}\|_{L^1(Q_T)}). \end{aligned} $$

Recall now that $D_\psi$ is bounded (Lemma 3.4), and observe that since $\{\theta_{\Sigma,\tau}\}$ is equiintegrable (owing to the entropy bound (4.13)) it has bounded $L^1(Q_T)$ norm. Due to our standing assumption (4.5) the last term is $o(1)$ as $k \to \infty$ and the proof is complete. $\square$

**4.3. Strong convergence of the approximate densities.** The main goal of this section is to establish improved compactness and therefore strong convergence of the reconstructions $\theta_{\mathcal{T},\tau}$ and $\rho_{\mathcal{T},\tau}$. Since the Fokker-Planck equation is linear this is not strictly required in order to prove the convergence of the scheme in Section 4.6 (see also Remark 4.8), and this should rather be read as a separate result of independent interest.

As a preliminary, we derive a uniform estimate on the discrete $L_t^2 \dot{H}_x^1$ semi-norm of $\sqrt{\theta_{\mathcal{T},\tau}}$, where the discrete TPFA $\dot{H}^1$ semi-norm is classically defined as

$$ |\boldsymbol{u}|_{1,\mathcal{T}}^2 = \sum_{\sigma \in \Sigma} \frac{m_\sigma}{d_\sigma} (u_K - u_L)^2, \quad \boldsymbol{u} \in \mathbb{R}^\mathcal{T}. $$

Starting from (3.11), we first rearrange the Fisher information as the sum of a linear part plus the $\dot{H}^1$ semi-norm

$$ (4.17) \qquad \frac{1}{2} \mathcal{R}_\psi(\boldsymbol{\rho}) = \sum_{\sigma \in \Sigma} \frac{m_\sigma \pi_\sigma}{d_\sigma} \left| \sqrt{\frac{\rho_K}{\pi_K}} - \sqrt{\frac{\rho_L}{\pi_L}} \right|^2 = \mathcal{I}(\boldsymbol{\rho}) + |\sqrt{\boldsymbol{\rho}}|_{1,\mathcal{T}}^2, $$

where

$$ (4.18) \qquad \mathcal{I}(\boldsymbol{\rho}) := \sum_{K \in \mathcal{T}} \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{d_\sigma} \rho_K \left( \sqrt{\frac{\pi_L}{\pi_K}} - 1 \right). $$

**Lemma 4.3.** There exists $C$ depending on $V$ and $\zeta$ (but neither on $\tau$ nor on size($\mathcal{T}$)) such that

$$ (4.19) \qquad |\mathcal{I}(\boldsymbol{\rho})| \le C \boldsymbol{\rho}[\Omega] + \frac{1}{2} |\sqrt{\boldsymbol{\rho}}|_{1,\mathcal{T}}^2, \quad \forall \boldsymbol{\rho} \in \mathbb{R}_+^\mathcal{T}. $$

As a consequence, there exists $C_T$ uniform with respect to $\mathcal{T}$, $\tau$ such that the solution of our scheme satisfies

$$ (4.20) \qquad \sum_{n=0}^{N-1} \tau \left| \sqrt{\boldsymbol{\theta}^{n+1/2}} \right|_{1,\mathcal{T}}^2 \le C_T. $$

*Proof.* Let us first focus on (4.19). Bearing in mind that $\pi_K = e^{-V_K}$, we write first $\sqrt{\frac{\pi_L}{\pi_K}} - 1 = e^{\frac{V_K-V_L}{2}} - 1$ in (4.18). Applying the mean value theorem $e^x - 1 = x + \frac{x^2}{2}e^y$ for some $y$ between 0 and $x$, we further split, for any $\boldsymbol{\rho} \in \mathbb{R}_+^\mathcal{T}$

$$ \mathcal{I}(\boldsymbol{\rho}) = \mathcal{I}_1(\boldsymbol{\rho}) + \mathcal{I}_2(\boldsymbol{\rho}) $$
with

$$ \mathcal{I}_1(\boldsymbol{\rho}) := \frac{1}{2} \sum_{\sigma \in \Sigma} \frac{m_{\sigma}}{d_{\sigma}} (\rho_K - \rho_L)(V_K - V_L) $$

$$ \mathcal{I}_2(\boldsymbol{\rho}) := \frac{1}{8} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \Sigma_K} \frac{m_{\sigma}}{d_{\sigma}} \rho_K (V_K - V_L)^2 e^{y_{K,\sigma}^{n+1/2}} $$

for some $y_{K,\sigma}^{n+1/2}$ between 0 and $(V_K - V_L)/2$. As $\rho_K \ge 0$, and using the regularity of $V$, we get that

$$ \mathcal{I}_2(\boldsymbol{\rho}) \le e^{\frac{\max V - \min V}{2}} \|\nabla V\|_{\infty}^2 \frac{1}{8} \sum_{K \in \mathcal{T}} \rho_K \sum_{\sigma \in \Sigma_K} m_{\sigma} d_{\sigma} \stackrel{(4.1)}{\le} C \sum_{K \in \mathcal{T}} \rho_K m_K = C \boldsymbol{\rho}[\Omega] $$

with $C$ depending only on $V$ and $\zeta$. For the $\mathcal{I}_1$ term we use next the identities $a - b = (\sqrt{a} - \sqrt{b})(\sqrt{a} + \sqrt{b})$ and $ab \le (a^2 + b^2)/4$ to get

$$ \mathcal{I}_1(\boldsymbol{\rho}) \le \frac{1}{2} |\sqrt{\boldsymbol{\rho}}|_{1,\mathcal{T}}^2 + \frac{1}{8} \sum_{\sigma \in \Sigma} \frac{m_{\sigma}}{d_{\sigma}} (\sqrt{\rho_K} + \sqrt{\rho_L})^2 (V_K - V_L)^2. $$

The Lipschitz continuity of $V$ and the elementary inequality $(\sqrt{a} + \sqrt{b})^2 \le 2(a + b)$ then yield

$$ \begin{aligned} \mathcal{I}_1(\boldsymbol{\rho}) &\le \frac{1}{2} |\sqrt{\boldsymbol{\rho}}|_{1,\mathcal{T}}^2 + \frac{\|\nabla V\|_{\infty}^2}{4} \sum_{K \in \mathcal{T}} \rho_K \sum_{\sigma \in \Sigma_K} m_{\sigma} d_{\sigma} \\ &\stackrel{(4.1)}{\le} \frac{1}{2} |\sqrt{\boldsymbol{\rho}}|_{1,\mathcal{T}}^2 + \frac{\|\nabla V\|_{\infty}^2}{4} \sum_{K \in \mathcal{T}} \rho_K 2d\zeta m_K \\ &= \frac{1}{2} |\sqrt{\boldsymbol{\rho}}|_{1,\mathcal{T}}^2 + C \boldsymbol{\rho}[\Omega], \end{aligned} $$

for $C$ depending again only on $V$, $\zeta$. Combining the above elements provides (4.19).

Turning now to (4.20), observe from (4.17)-(4.19) that

$$ \frac{1}{2} |\sqrt{\boldsymbol{\rho}}|_{1,\mathcal{T}}^2 \le \frac{1}{2} \mathcal{R}_{\psi}(\boldsymbol{\rho}) + C \boldsymbol{\rho}[\Omega], \quad \forall \boldsymbol{\rho} \in \mathbb{R}_+^{\mathcal{T}}. $$

Summing over $n$ gives, for any discrete curve $(\boldsymbol{\theta}^{n+1/2})_{n=0}^{N-1}$,

$$ (4.21) \qquad \sum_{n=0}^{N-1} \tau \left| \sqrt{\boldsymbol{\theta}^{n+1/2}} \right|_{1,\mathcal{T}}^2 \le \sum_{n=0}^{N-1} \tau \mathcal{R}_{\psi}(\boldsymbol{\theta}^{n+1/2}) + 2C \|\theta_{\mathcal{T},\tau}\|_{L^1(Q_T)} $$

When evaluated for our particular solution of the discrete scheme, the first term in the right-hand side is bounded by Lemma 3.4. Recalling $\theta_K^{n+1/2} \le \frac{\rho_K^{n+1} + \rho_K^n}{2}$ from (3.1) and the mass conservation $\rho^{n+1}[\Omega] = \rho^n[\Omega] = \rho^0[\Omega]$ from Proposition 3.1, we see that $\|\theta_{\mathcal{T},\tau}\|_{L^1(Q_T)} \le T\rho^0[\Omega]$ and the proof is complete. $\square$

We can now upgrade the previous weak $L^1(Q_T)$ convergence of the approximate densities $\rho_{\mathcal{T},\tau}$, $\theta_{\mathcal{T},\tau}$ into strong convergence to a common limit.

**Proposition 4.4.** Let $\rho, \theta$ be as in (4.11) and assume (4.5) as in Lemma 4.2. Then $\rho = \theta$ and, up to extraction of a subsequence,

$$ (4.22) \qquad \theta_{\mathcal{T}^k, \tau^k} \underset{k \to +\infty}{\longrightarrow} \rho \qquad \text{strongly in } L^1(Q_T), $$

$$ (4.23) \qquad \rho_{\mathcal{T}^k, \tau^k} \underset{k \to +\infty}{\longrightarrow} \rho \qquad \text{strongly in } L^1(Q_T). $$
*Proof.* Observe first from Lemma 4.2 that the weak limits must coincide $\rho = \theta$, so it suffices to prove that $\theta_{\mathcal{T}^{k,\tau^k}} \to \rho$ strongly in $L^1(Q_T)$. Our proof relies on a combination of an Aubin-Lions-Simon concentration-compactness argument and a monotone Minty's trick, already proposed in [3]. This will however need some adaptation of results from [50] to our specific setup, which we defer to Proposition B.1 in the appendix.

Let $f$ be an increasing and bounded function from $\mathbb{R}_+$ to $\mathbb{R}$ such that $f(0) = 0$ and such that $z \mapsto f(z^2)$ is 1-Lipschitz continuous (typically $f(z) = \tanh\sqrt{z}$). Define next the piecewise constant and discrete functions

$$f_{\mathcal{T},\tau} := f(\theta_{\mathcal{T},\tau}) \quad \text{and} \quad \mathbf{f}^{n+1/2} := f(\boldsymbol{\theta}^{n+1/2}).$$

Recall from (4.11) that

$$\rho_{\mathcal{T}^k, \tau^k} \xrightarrow{k \to +\infty} \rho \quad \text{weakly in } L^1(Q_T),$$

and observe that, since $f$ is bounded, $\{f_{\mathcal{T},\tau}\}$ is bounded in $L^\infty(Q_T)$ and therefore

$$f_{\mathcal{T}^k, \tau^k} \xrightarrow{k \to +\infty} \mathbf{f} \quad \text{weakly-} * \text{ in } L^\infty(Q_T)$$

for some $\mathbf{f} \in L^\infty(Q_T)$ and possibly up to extraction of a subsequence. We aim to use Proposition B.1 from the Appendix to guarantee that, with suitable time-compactness on $\{\rho_{\mathcal{T},\tau}\}$ and space compactness on $\{f_{\mathcal{T},\tau}\}$, we can pass to the limit in the product $\rho_{\mathcal{T},\tau} f_{\mathcal{T},\tau} \to \rho \mathbf{f}$ in the sense of measures.

We first focus on the space compactness. Since $z \mapsto f(z^2)$ is 1-Lipschitz, we have that

$$|\mathbf{f}^{n+1/2}|_{1, \mathcal{T}^k} \leq |\sqrt{\boldsymbol{\theta}^{n+1/2}}|_{1, \mathcal{T}^k}, \quad n \geq 0,$$

hence we deduce from Lemma 4.3 that

$$(4.24) \qquad \sum_{n=0}^{N-1} \tau |\mathbf{f}^{n+1/2}|_{1, \mathcal{T}^k}^2 \leq C_T.$$

A slight adaptation of [20, Lemma 9.3] shows first of all that the limit $\mathbf{f} \in L^2(0, T; H^1(\Omega))$, and moreover controls the $L^2$ space difference quotients by the discrete $H^1$ norm in the following quantitative sense: there exists a constant $C > 0$ only depending on $\Omega$ such that, for any compact subset $\omega \subset \subset \Omega$ and $h \in \mathbb{R}^d$ such that $|h| < \text{dist}(\omega, \partial\Omega)$,

$$\int_0^T \int_\omega |f_{\mathcal{T},\tau}(t, x+h) - f_{\mathcal{T},\tau}(t, x)|^2 \leq |h| [|h| + C \text{ size}(\mathcal{T})] \sum_n \tau |\sqrt{\mathbf{f}^{n+1/2}}|_{1, \mathcal{T}}^2 \leq C_T |h|.$$

This gives in turn the (suboptimal) $L^1$ difference quotient estimate

$$(4.25) \qquad \int_0^T \int_\omega |f_{\mathcal{T}^k, \tau^k}(t, x+h) - f_{\mathcal{T}^k, \tau^k}(t, x)| \leq C \sqrt{|h|}, \quad \text{uniformly in } k.$$

Turning now to the compactness in time for $\rho$, take any arbitrary test function $\varphi \in C_c^\infty(Q_T)$, define $(\varphi^n)_n \in \mathbb{R}^T$ by

$$\varphi_K^n := \frac{1}{m_K} \int_K \varphi(t^n, x)dx, \quad K \in \mathcal{T}, n = 0, \dots, N-1,$$
and compute

$$
\begin{aligned}
\left| \sum_{n=0}^{N-1} \sum_{K \in \mathcal{T}} m_K (\rho_K^{n+1} - \rho_K^n) \varphi_K^n \right| & \stackrel{(2.9)}{=} \left| \sum_{n=0}^{N-1} \tau \sum_{\sigma \in \Sigma} m_\sigma F_{K\sigma}^{n+1/2} (\varphi_L^n - \varphi_K^n) \right| \\
& \le \max_{n, \sigma} \left| \frac{\varphi_L^n - \varphi_K^n}{d_\sigma} \right| \left| \sum_{n=0}^{N-1} \tau \sum_{\sigma \in \Sigma} m_\sigma d_\sigma \left| F_{K\sigma}^{n+1/2} \right| \right| \\
& = \max_{n, \sigma} \left| \frac{\varphi_L^n - \varphi_K^n}{d_\sigma} \right| \| F_{\Sigma, \tau} \|_{L^1(Q_T)}.
\end{aligned}
$$

Taking $\varepsilon = 1$ and $A = Q_T$ in (4.15) gives that $\|F_{\Sigma,\tau}\|_{L^1(Q_T)} \le C$ is bounded uniformly in $\mathcal{T}, \tau$. Moreover, it is shown in [3, §4.4] that

$$
\left| \frac{\varphi_L^n - \varphi_K^n}{d_\sigma} \right| \le (1 + 2\zeta) \|\nabla \varphi\|_{\infty},
$$

so that altogether we get the $\mathcal{M}(0, T; (W^{1,\infty}(\Omega')))$ estimate

$$
(4.26) \quad \left| \langle \partial_t \rho_{\mathcal{T},\tau}, \varphi \rangle \right| = \left| \sum_{n=0}^{N-1} \tau \sum_{K \in \mathcal{T}} m_K \frac{\rho_K^{n+1} - \rho_K^n}{\tau} \varphi_K^n \right| \le C \|\nabla \varphi\|_{\infty}.
$$

We are now in position of rigorously applying our Aubin-Lions compactness from Proposition B.1: (i) $\rho_{\mathcal{T}^k, \tau^k}$ is bounded in $L^\infty(0, T; L^1(\Omega))$ (conservation of mass), hence in $L^1(Q_T)$, and $f_{\mathcal{T}^k, \tau^k}$ is bounded in $L^\infty(Q_T)$, the equiintegrability (ii) follows from the entropy bound (4.10), the space compactness (iii) is exactly given by (4.25), and the time compactness (iv) is just (4.26). We conclude that, up to extraction of a subsequence if need be,

$$
(4.27) \qquad \rho_{\mathcal{T}^k, \tau^k} f_{\mathcal{T}^k, \tau^k} \xrightarrow{k \to +\infty} \rho f \quad \text{in } \mathcal{M}(Q_T),
$$

and from Lemma 4.2 also

$$
(4.28) \qquad \theta_{\mathcal{T}^k, \tau^k} f_{\mathcal{T}^k, \tau^k} \xrightarrow{k \to +\infty} \rho f \quad \text{in } \mathcal{M}(Q_T).
$$

Fix now any $z \in \mathbb{R}_+$ and $\varphi \in C(\bar{Q}_T)$ with $\varphi \ge 0$. Recalling that $f_{\mathcal{T},\tau} = f(\theta_{\mathcal{T},\tau})$, the monotonicity of $f$ gives

$$
\int_{Q_T} (\theta_{\mathcal{T}^k, \tau^k} - z)(f_{\mathcal{T}^k, \tau^k} - f(z)) \varphi \ge 0
$$

and therefore owing to (4.28)

$$
\int_{Q_T} (\rho - z)(f - f(z)) \varphi = \lim_{k \to +\infty} \int_{Q_T} (\theta_{\mathcal{T}^k, \tau^k} - z)(f_{\mathcal{T}^k, \tau^k} - f(z)) \varphi \ge 0.
$$

Since $\varphi \ge 0$ is arbitrary we see that

$$
\forall z \ge 0, \quad (\rho - z)(f - f(z)) \ge 0 \quad \text{a.e. in } Q_T,
$$

which implies

$$
(4.29) \qquad f = f(\rho) \quad \text{a.e. in } Q_T.
$$

Finally, we claim that the non-negative sequence $(R_k)_k \subset L^1(Q_T)$ defined by

$$
R_k := (\theta_{\mathcal{T}^k, \tau^k} - \rho)(f(\theta_{\mathcal{T}^k, \tau^k}) - f(\rho)) \ge 0
$$
converges to 0 strongly in $L^1(Q_T)$. Indeed, thanks to (4.28) and (4.29):

$$ 
\begin{aligned}
0 \le \|R_k\|_{L^1(Q_T)} &= \int_{Q_T} R_k = \int_{Q_T} (\theta_{\mathcal{T}^k, \tau^k} - \rho)(f(\theta_{\mathcal{T}^k, \tau^k}) - f(\rho)) \\
&= \int_{Q_T} \theta_{\mathcal{T}^k, \tau^k} f_{\mathcal{T}^k, \tau^k} - \int_{Q_T} \rho f_{\mathcal{T}^k, \tau^k} - \int_{Q_T} \theta_{\mathcal{T}^k, \tau^k} f(\rho) + \int_{Q_T} \rho f(\rho) \\
&\xrightarrow{k \to \infty} \int_{Q_T} \rho f - \int_{Q_T} \rho f - \int_{Q_T} \theta f(\rho) + \int_{Q_T} \rho f(\rho) = 0
\end{aligned}
$$

since we already proved that $\rho = \theta$ and $f = f(\rho)$. This strong convergence implies almost everywhere convergence $R_k(t, x) \to 0$ in $Q_T$, up to a subsequence if needed. As $f$ is increasing, this also implies almost everywhere convergence of $\theta_{\mathcal{T}^k, \tau^k}$ towards $\rho$. Vitali's convergence theorem then provides (4.22), and (4.23) finally follows from Lemma 4.2. $\square$

### 4.4. Asymptotic lower bound for the discrete Fisher information.
In this section we show that the $\Gamma$-lim inf of the total dissipation functional

$$ \sum_{n=0}^{N-1} \tau \mathcal{R}_{\psi}(\boldsymbol{\theta}^{n+1/2}) $$

with respect to the weak $L^1$ convergence is bounded from below by the total dissipation of the continuous system, i.e. the dissipation rate (3.12) integrated in time. As this is a statement on the functional itself, throughout this section we will consider arbitrary discrete curves $(\boldsymbol{\theta}^{n+1/2})_n$ that possibly do not solve (2.9)–(2.10)–(2.11).

Let us first consider the easy case of a trivial background potential $V \equiv 0$. In that case, according to (3.11), the dissipation is exactly given by the semi-norm

$$ \frac{1}{2} \mathcal{D}_{\psi^*}(\boldsymbol{\rho}) = |\sqrt{\boldsymbol{\rho}}|_{1, \tau}^2 = \sum_{\sigma \in \Sigma} \frac{m_{\sigma}}{d_{\sigma}} |\sqrt{\rho_K} - \sqrt{\rho_L}|^2. $$

Using the previous adaptation of [20, Lemma 9.3] to control $L^2$ space difference quotients by the discrete seminorms, there exists a constant $C > 0$ only depending on $\Omega$ such that, for any compact subset $\omega \subset \Omega$ and $h \in \mathbb{R}^d$ such that $|h| < \text{dist}(\omega, \partial\Omega)$,

$$ (4.30) \quad \int_0^T \int_\omega \left| \sqrt{\theta_{\mathcal{T},\tau}(t, x+h)} - \sqrt{\theta_{\mathcal{T},\tau}(t, x)} \right|^2 \le \sum_n \tau |\sqrt{\theta^{n+1/2}}|_{1,\tau}^2 |h| [|h| + C \, \text{size}(\mathcal{T})]. $$

In particular, consider any sequence of meshes and time steps such that $\text{size}(\mathcal{T}^k) \to 0$ and $\tau^k \to 0$, and let $(\theta_k^{n+1/2})_{n=0}^{N-1}$ be any associated discrete curve (again, not necessarily solution to our discrete scheme). Assume that the reconstruction $\theta_{\mathcal{T}^k, \tau^k} \in L^1(Q_T)$ from (4.7) converges as

$$ \theta_{\mathcal{T}^k, \tau^k} \to \theta \quad \text{weakly in } L^1(Q_T). $$

Because $f(a, b) = |\sqrt{a} - \sqrt{b}|^2$ is convex and continuous, the left-hand side of (4.30) is lower-semicontinuous for the weak $L^1$ convergence, hence

$$ \int_0^T \left\| \sqrt{\theta(t, \cdot + h)} - \sqrt{\theta(t, \cdot)} \right\|_{L^2(\omega)}^2 \mathrm{d}t \le |h|^2 \liminf_{k \to \infty} \sum_n \tau^k |\sqrt{\theta_k^{n+1/2}}|_{1, \mathcal{T}^k}^2 $$
and therefore, by classical characterization of $H^1(\Omega)$ by difference quotients,

$$ (4.31) \qquad \int_0^T \|\nabla \sqrt{\theta(t, \cdot)}\|_{L^2(\Omega)}^2 dt \le \liminf_{k \to \infty} \sum_n \tau^k |\sqrt{\theta_k^{n+1/2}}|_{1, \mathcal{T}^k}^2. $$

This settles the case $V \equiv 0$.

In order to prove the analogue result in the presence of a non-zero potential $V$, we will rewrite below the dissipation functional as the sum of the above $H^1$ seminorm plus a linear term, which can be dealt with easily at least when $\nabla V \cdot n_{\partial\Omega} = 0$ on the boundary. More precisely, at the continuous level, if $\nabla V \cdot n_{\partial\Omega} = 0$ on the boundary we have the identity

$$ (4.32) \qquad \int \pi \left| \nabla \sqrt{\rho} \right|^2 = \frac{1}{2} \int \rho \left( \frac{|\nabla V|^2}{2} - \Delta V \right) + \int |\nabla \sqrt{\rho}|^2. $$

This formula can be directly related to the expression for the discrete Fisher information (3.11), decomposed into $\frac{1}{2}\mathcal{R}_\psi(\rho) = \mathcal{I}(\rho) + |\sqrt{\rho}|_{1, \mathcal{T}}^2$ as in (4.17)-(4.18). The case $\nabla V \cdot n_{\partial\Omega} \ne 0$ on the boundary will be handled via an approximation argument from [17].

**Proposition 4.5.** Let $(\theta_k^{n+1/2})_{n=0}^{N-1} \in (\mathbb{R}_+^T)^N$ be a given discrete curve associated with $\mathcal{T}^k, \tau^k$ with reconstruction $\theta_{\mathcal{T}^k, \tau^k} \in L^1(Q_T)$ as in (4.7), and suppose that

$$ \theta_{\mathcal{T}^k, \tau^k} \rightharpoonup \theta, \quad \text{weakly in } L^1(Q_T) $$

for some $\theta \in L^1(Q_T)$. Then

$$ \liminf_{k \to \infty} \sum_n \tau^k \mathcal{R}_\psi(\theta_k^{n+1/2}) \ge 2 \int_{Q_T} \pi \left| \nabla \sqrt{\frac{\theta}{\pi}} \right|^2. $$

*Proof.* Consider first $\nabla V \cdot n_{\partial\Omega} = 0$ on the boundary, and assume that the $\liminf$ in the statement is finite (otherwise the statement is vacuous). By (4.21) we readily obtain that

$$ \liminf_{k \to \infty} \sum_n \tau^k \left| \sqrt{\theta_k^{n+1/2}} \right|_{1, \mathcal{T}}^2 \le 2C \limsup_{k \to \infty} \|\theta_{\mathcal{T}^k, \tau^k}\|_{L^1(Q_T)} + \liminf_{k \to \infty} \sum_n \tau^k \mathcal{R}_\psi(\theta_k^{n+1/2}) < +\infty $$

is finite. By the previous considerations for $V \equiv 0$ we see that (4.31) holds, hence comparing (4.17)-(4.18) on the one hand and (4.32) on the other hand, clearly it suffices to show:

$$ (4.33) \qquad \liminf_{k \to \infty} \sum_n \tau^k \mathcal{I}(\theta_k^{n+1/2}) \ge \frac{1}{2} \int \theta \left( \frac{|\nabla V|^2}{2} - \Delta V \right). $$

To this end we first observe that

$$ \sqrt{\frac{\pi_L}{\pi_K}} - 1 = \exp \left( \frac{V_K - V_L}{2} \right) - 1 \ge \frac{V_K - V_L}{2} + \frac{(V_K - V_L)^2}{8}. $$

Therefore we can write $\mathcal{I}(\rho) \ge \mathcal{I}_1(\rho) + \mathcal{I}_2(\rho)$ where

$$ \mathcal{I}_1(\rho) := \sum_{K \in \mathcal{T}} \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{d_\sigma} \rho_K \frac{(V_K - V_L)^2}{8}, $$

$$ \mathcal{I}_2(\rho) := \sum_{K \in \mathcal{T}} \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{d_\sigma} \rho_K \frac{V_K - V_L}{2}. $$

For fixed $K \in \mathcal{T}$ we have

$$ \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{d_\sigma} \frac{V_K - V_L}{2} \ge -\frac{1}{2} \int_K \Delta V - C \, \text{size}(\mathcal{T}^k), $$
where C only depends on V.

Since V belongs to $C^2(\bar{\Omega})$, and because of the orthogonality condition (2.1), there holds

$$ \frac{V_K - V_L}{d_\sigma} = -\nabla V(x_K) \cdot n_{K\sigma} + \mathcal{O}(d_\sigma) = -\nabla V(x) \cdot n_{K\sigma} + \mathcal{O}(d_\sigma + |x - x_K|) $$

for any $x \in K$. Therefore, we deduce from (4.1) and (4.2) that

$$ \frac{V_K - V_L}{d_\sigma} = -\nabla V(x) \cdot n_{K\sigma} + \mathcal{O}(\mathrm{size}(\mathcal{T}^k)). $$

Using the isotropy condition (4.3) and then integrating w.r.t. $x \in K$ yields

$$ (4.34) \quad \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{d_\sigma} \frac{|V_K - V_L|^2}{8} \ge (1 - \varepsilon_{\mathcal{T}^k}) \left( \int_K \frac{|\nabla V|^2}{4} - C \, \mathrm{size}(\mathcal{T}^k) \right) = \int_K \frac{|\nabla V|^2}{4} + o(1), $$

where the $o(1)$ remainder is uniform both in $K \in \mathcal{T}^k$ and $n = 1, \dots, N-1$ as $k \to \infty$. Hence we obtain that

$$ \sum_n \tau^k \mathcal{I}_1(\boldsymbol{\theta}_k^{n+1/2}) \ge \int_{Q_T^k} \theta_{\mathcal{T}^k, \tau^k} \frac{|\nabla V|^2}{4} + o(1) $$

and

$$ \sum_n \tau^k \mathcal{I}_2(\boldsymbol{\theta}_k^{n+1/2}) \ge - \int_{Q_T} \theta_{\mathcal{T}^k, \tau^k} \Delta V + o(1), $$

where we have set

$$ Q_T^k := (0, T) \times (\Omega \setminus \Omega_{\mathrm{iso}}^k). $$

Since we assume $\theta_{\mathcal{T}^k, \tau^k} \to \theta$ weakly in $L^1(Q_T)$, the inequality involving $\mathcal{I}_2$ immediately passes to the limit. For the $\mathcal{I}_1$ lower bound, the weak convergence implies that $\{\theta_{\mathcal{T}^k, \tau^k}\}_k$ is equi-integrable. Owing to our standing assumption (4.4) we see that $\mathrm{meas}(Q \setminus Q_T^k) \to 0$ hence the $\mathcal{I}_1$ inequality also passes to the limit and our claim (4.33) follows.

Let us finally settle the case $\nabla V \cdot n_{\partial\Omega} \ne 0$ on the boundary. Fix $\varepsilon > 0$ and take an approximation $V^\varepsilon \in C_c^\infty(\mathbb{R}^d)$ of V satisfying

$$ \nabla V^\varepsilon \cdot n_{\partial\Omega} = 0 \quad \text{on } \partial\Omega, \quad \|V - V^\varepsilon\|_{W^{1,p}(\Omega)} \le \varepsilon, $$

for some $p \in [1, \infty)$ to be chosen later. The existence of such a function for arbitrarily large but finite $p \ge 1$ is due to Droniou [17]. Define $\mathcal{I}_2^\varepsilon$ in the obvious way, simply substituting $V^\varepsilon$ for V in the previous definition of $\mathcal{I}_2$, and for any discrete function $\rho \in \mathbb{R}_+^T$ decompose now $\frac{1}{2}\mathcal{R}_\psi(\rho) = \mathcal{I}(\rho) + |\sqrt{\rho}|_{1,\tau}^2 \ge \mathcal{I}_1(\rho) + \mathcal{I}_2(\rho) + |\sqrt{\rho}|_{1,\tau}^2$ as

$$ (4.35) \quad \frac{1}{2} \mathcal{R}_\psi(\rho) \ge [\mathcal{I}_1(\rho) + \mathcal{I}_2^\varepsilon(\rho) + |\sqrt{\rho}|_{1,\tau}^2] + [\mathcal{I}_2(\rho) - \mathcal{I}_2^\varepsilon(\rho)]. $$
Let us first estimate the difference $\mathcal{I}_2 - \mathcal{I}_2^\varepsilon$. To this end, pick an arbitrary $\eta > 0$ small, and for any $\theta \in \mathbb{R}_+^T$ with spatial reconstruction $\theta_T \in L^1(\Omega)$ write

$$
\begin{aligned}
|\mathcal{I}_2(\boldsymbol{\theta}) - \mathcal{I}_2^\varepsilon(\boldsymbol{\theta})| & \le \sum_{\sigma \in \Sigma} \frac{m_\sigma}{2d_\sigma} |\theta_K - \theta_L| |V_K - V_K^\varepsilon - V_L + V_L^\varepsilon| \\
& = \sum_{\sigma \in \Sigma} \frac{m_\sigma}{2d_\sigma} |\sqrt{\theta_K} - \sqrt{\theta_L}| \cdot (\sqrt{\theta_K} + \sqrt{\theta_L}) |V_K - V_K^\varepsilon - V_L + V_L^\varepsilon| \\
& \le \left( \eta |\sqrt{\boldsymbol{\theta}}|_{1,T}^2 + \frac{1}{4\eta} \sum_{\sigma \in \Sigma} \frac{m_\sigma}{d_\sigma} |\sqrt{\theta_K} + \sqrt{\theta_L}|^2 |V_K - V_K^\varepsilon - V_L + V_L^\varepsilon|^2 \right) \\
& \le \left( \eta |\sqrt{\boldsymbol{\theta}}|_{1,T}^2 + \frac{1}{2\eta} \sum_{\sigma \in \Sigma} \frac{m_\sigma}{d_\sigma} (\theta_K + \theta_L) |V_K - V_K^\varepsilon - V_L + V_L^\varepsilon|^2 \right).
\end{aligned}
$$

Similar arguments as those employed to establish (4.34) show that

$$ \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{2d_\sigma} |V_K - V_K^\varepsilon - V_L + V_L^\varepsilon|^2 \le \frac{1}{m_K} \int_K |\nabla(V - V^\varepsilon)|^2 + C_\varepsilon (\varepsilon_T + \text{size}(T)), \quad \forall K \in \mathcal{T}_{\text{iso}} $$

where the constant $C_\varepsilon$ depends on $V - V^\varepsilon$ but not on the mesh. On the other hand for $K \notin \mathcal{T}_{\text{iso}}^k$ we simply write

$$
\begin{aligned}
\forall K \notin \mathcal{T}_{\text{iso}}^k, \quad & \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{2d_\sigma} |V_K - V_K^\varepsilon - V_L + V_L^\varepsilon|^2 \le \sum_{\sigma \in \Sigma_K} \frac{m_\sigma d_\sigma}{2} (\|\nabla V\|_\infty + \|\nabla V^\varepsilon\|_\infty)^2 \\
& \stackrel{(4.1)}{\le} d \zeta m_K (\|\nabla V\|_\infty + \|\nabla V^\varepsilon\|_\infty)^2 \\
& \le C_\varepsilon m_K,
\end{aligned}
$$

where $C_\varepsilon$ again depends on $\varepsilon$ but not on $\mathcal{T}^k$. Hence we find

$$
|\mathcal{I}_2(\boldsymbol{\theta}) - \mathcal{I}_2^\varepsilon(\boldsymbol{\theta})| \le \eta |\sqrt{\boldsymbol{\theta}}|_{1,T}^2 + \frac{1}{\eta} \int_{\Omega} \theta_T |\nabla(V - V^\varepsilon)|^2 \\
+ C_\varepsilon \left( (\varepsilon_T + \text{size}(T)) \|\theta_T\|_{L^1(\Omega)} + \|\theta_T\|_{L^1(\Omega \setminus \Omega_{\text{iso}})} \right),
$$

which inserted into (4.35) yields

$$
\begin{aligned}
\frac{1}{2} \mathcal{R}_\psi(\boldsymbol{\theta}) \ge & \left[ \mathcal{I}_1(\boldsymbol{\theta}) + \mathcal{I}_2^\varepsilon(\boldsymbol{\theta}) + (1-\eta)|\sqrt{\boldsymbol{\theta}}|_{1,T}^2 \right] \\
& - \frac{1}{\eta} \int_{\Omega} \theta_T |\nabla(V - V^\varepsilon)|^2 - C_\varepsilon \left( (\varepsilon_T + \text{size}(T)) \|\theta_T\|_{L^1(\Omega)} + \|\theta_T\|_{L^1(\Omega \setminus \Omega_{\text{iso}})} \right).
\end{aligned}
$$

Evaluating for a discrete curve $(\theta_k^{n+1/2})$ and summing over $n$, the first three terms in the right-hand side pass to the $\liminf$ as soon as $\theta_{\mathcal{T}^k, \tau^k} \to \theta$ weakly in $L^1(Q_T)$, exactly as in the previous case where $\nabla V \cdot n_{\partial\Omega} = 0$ on the boundary. Moreover the non-isotropic term vanishes as $k \to \infty$ – as before, $\theta_{\mathcal{T}^k, \tau^k}$ is equiintegrable and $Q_T^k = (0, T) \times (\Omega \setminus \Omega_{\text{iso}})$ has vanishing measure – hence for fixed $\eta > 0$ and $\varepsilon > 0$ we obtain

$$
\begin{aligned}
\liminf_{k \to \infty} \sum_n \tau^k \mathcal{R}_\psi(\theta_k^{n+1/2}) & \\
\ge & \int_{Q_T} \theta \left( \frac{|\nabla V|^2}{2} - \Delta V^\varepsilon \right) + 2(1-\eta) \int_{Q_T} |\nabla \sqrt{\theta}|^2 - \frac{2}{\eta} \int_{Q_T} \theta |\nabla(V - V^\varepsilon)|^2
\end{aligned}
$$
This shows as a byproduct that, whenever the left-hand side is finite, $||\nabla\sqrt{\theta}||_{L^2(Q_T)}$ is finite too and $\nabla\theta = 2\sqrt{\theta}\nabla\sqrt{\theta} \in L_t^2 L_x^1$ due to $\sqrt{\theta} \in L_t^\infty L_x^2$. At this stage one would wish to substitute $-\int \theta \Delta V^\epsilon$ by the desired $+\int \nabla\theta \cdot \nabla V$. To this end we use the exact same strategy as before, but this time at the continuous level: Since $\nabla V^\epsilon \cdot n_{\partial\Omega}$ on the boundary and $\theta \in L_t^2 W_x^{1,1}$ we can legitimately integrate by parts

$$ 
\begin{aligned}
\left| \int_{Q_T} \theta \Delta V^\epsilon + \int_{Q_T} \nabla\theta \cdot \nabla V \right| &= \left| \int_{Q_T} \nabla\theta \cdot (\nabla V^\epsilon - \nabla V) \right| \\
&= 2 \left| \int_{Q_T} \nabla\sqrt{\theta} \cdot \sqrt{\theta}(\nabla V^\epsilon - \nabla V) \right| \\
&\le 2\eta \int_{Q_T} \left|\nabla\sqrt{\theta}\right|^2 + \frac{2}{\eta} \int_{Q_T} \theta |\nabla V^\epsilon - \nabla V|^2
\end{aligned}
$$

and thus

$$ 
\begin{aligned}
\liminf_{k \to \infty} \sum_n \tau^k \mathcal{R}(\theta_k^{n+1/2}) & \\
& \ge \int_{Q_T} \left( \theta \frac{|\nabla V|^2}{2} + \nabla\theta \cdot \nabla V \right) + (2-4\eta) \int_{Q_T} |\nabla\sqrt{\theta}|^2 - \frac{4}{\eta} \int_{Q_T} \theta |\nabla(V - V^\epsilon)|^2.
\end{aligned}
$$

Fix now $q > 1$ such that, for any $f \in H^1(\Omega)$,

$$ \|f\|_{L^{2q}(\Omega)} \le C_q \|f\|_{H^1(\Omega)} $$

for some constant $C_q > 0$ depending only on $\Omega$ and $q$. Choosing $p$ as $1/p + 1/q = 1$ we estimate in the last term

$$ 
\begin{aligned}
\int_{Q_T} \theta |\nabla(V - V^\epsilon)|^2 &\le \|\nabla(V - V^\epsilon)\|_{L^{2p}(\Omega)}^2 \int_0^T \|\theta(t, \cdot)\|_{L^q(\Omega)} \mathrm{d}t \\
&= \|\nabla(V - V^\epsilon)\|_{L^{2p}(\Omega)}^2 \int_0^T \|\sqrt{\theta}(t, \cdot)\|_{L^{2q}(\Omega)}^2 \mathrm{d}t \\
&\le C_q^2 \|\nabla(V - V^\epsilon)\|_{L^{2p}(\Omega)}^2 \|\sqrt{\theta}\|_{L^2(0,T; H^1(\Omega))}^2
\end{aligned}
$$

Choosing $\epsilon > 0$ sufficiently small so that

$$ C_q^2 \|\nabla(V - V^\epsilon)\|_{L^{2p}(\Omega)}^2 \le \eta^2 $$

gives

$$ 
\begin{aligned}
\liminf_{k \to \infty} \sum_n \tau^k \mathcal{R}(\theta_k^{n+1/2}) & \\
& \ge \int_{Q_T} \left( \theta \frac{|\nabla V|^2}{2} + \nabla\theta \cdot \nabla V + 2 |\nabla\sqrt{\theta}|^2 \right) - 4\eta \int_{Q_T} |\nabla\sqrt{\theta}|^2 - 4\eta, \\
& = 2\pi |\nabla\sqrt{\frac{\theta}{\pi}}|^2
\end{aligned}
$$

and since $\eta > 0$ was arbitrary the proof is complete. $\square$

<figure><img src="image_3.png" /></figure>
4.5. **Asymptotic lower bound for the discrete Benamou-Brenier functional.** In this section we establish a lower bound for the kinetic part of the dissipation, analogous to the previous section but now for

$$ \sum_{n=0}^{N-1} \tau \mathcal{D}_{\psi}(\boldsymbol{\theta}^{n+1/2}, \boldsymbol{F}^{n+1/2}) = \sum_{n=0}^{N-1} \tau \sum_{\sigma \in \Sigma} \frac{m_{\sigma}}{d_{\sigma}} \theta_{\sigma}^{n+1/2} \psi \left( \frac{d_{\sigma} F_{K\sigma}^{n+1/2}}{\theta_{\sigma}^{n+1/2}} \right). $$

The lower bound will be given by the Benamou-Brenier functional (see [57, Section 5.3.1]), which is the map given by $(\rho, F) \in L^1(Q_T; \mathbb{R}^{d+1}) \mapsto \int_{Q_T} B(\rho, F)$ with the function $B$ defined as in (1.5).

Just as in section 4.4 for the Fisher information, this statements is about how the geometry discretization allows a consistent interplay between the discrete $\mathcal{D}_{\psi}$ dissipation and the continuous Benamou-Brenier functional, and does not deal with particular solutions of our scheme properly speaking. Accordingly, consider $\theta \in L^1(Q_T)$ and $F \in L^1(Q_T; \mathbb{R}^d)$ that are obtained as weak $L^1$ limits of arbitrary sequences of reconstructed densities $\{\theta_{\tau^k, \tau^k}\}_k$ and fluxes $\{F_{\Sigma^k, \Sigma^k}\}_k$. Our strategy below is similar to [37], where unconditional convergence of discrete to continuous optimal transport models is proved. Recall that the Benamou-Brenier functional can be classically written (e.g. [57, Proposition 5.18]) as

$$ (4.36) \qquad \int_{Q_T} B(\theta, F) = \sup_{b \in C(\bar{Q}_T; \mathbb{R}^d)} \left\{ \langle b, F \rangle - \int_{Q_T} \theta \frac{|b|^2}{2} \right\}, $$

where the duality pairing is $\langle b, F \rangle = \int_{Q_T} b \cdot F$. Consider first the case when (4.36) is finite. Then, for any arbitrary small $\eta > 0$ we can find $b \in C(\bar{Q}_T; \mathbb{R}^d)$ such that

$$ (4.37) \qquad \int_{Q_T} B(\theta, F) \le \langle b, F \rangle - \int_{Q_T} \theta \frac{|b|^2}{2} + \eta \le \lim_{k \to \infty} \langle b, F_{\Sigma^k, \tau^k} \rangle - \int_{Q_T} \theta_{\tau^k, \tau^k} \frac{|b|^2}{2} + \eta, $$

and by density we can actually assume that $b \in C^1(\bar{Q}_T; \mathbb{R}^d)$. At the discrete level, the analogue of (4.36) is (by definition of the $\psi, \psi^*$ convex duality)

$$ \sum_n \tau \sum_\sigma \frac{m_\sigma \theta_\sigma^{n+1/2}}{d_\sigma} \psi \left( \frac{F_\sigma^{n+1/2} d_\sigma}{\theta_\sigma^{n+1/2}} \right) = \sup_{b \in \mathbb{F}_\tau} \left\{ \sum_{n, \sigma} \tau m_\sigma d_\sigma b_{K\sigma}^{n+1/2} F_{K\sigma}^{n+1/2} - \sum_n \tau \sum_\sigma \frac{m_\sigma \theta_\sigma^{n+1/2}}{d_\sigma} \psi^* \left( b_{K\sigma}^{n+1/2} d_\sigma \right) \right\}. $$

Let us take $b$ as in (4.37), and define

$$ b_{K\sigma}^{n+1/2} := \frac{1}{\tau \operatorname{meas}(\Delta_{\sigma})} \int_{t^n}^{t^{n+1}} \int_{\Delta_{\sigma}} b \cdot n_{K\sigma} $$

so that

$$ \langle b, F_{\Sigma, \tau} \rangle = \sum_{n, \sigma} \tau m_{\sigma} d_{\sigma} b_{K\sigma}^{n+1/2} F_{K\sigma}^{n+1/2}. $$

Moreover owing to $\theta_{\sigma}^{n+1/2} = \sqrt{\theta_{K}^{n+1/2} \theta_{L}^{n+1/2}} \le \frac{1}{2}(\theta_{K}^{n+1/2} + \theta_{L}^{n+1/2})$ there holds

$$ \sum_n \tau \sum_\sigma \frac{m_\sigma \theta_\sigma^{n+1/2}}{d_\sigma} \psi^* \left( b_{K\sigma}^{n+1/2} d_\sigma \right) \le \sum_n \tau \sum_K \theta_K^{n+1/2} \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{2d_\sigma} \psi^* \left( b_{K\sigma}^{n+1/2} d_\sigma \right). $$
Note that $|b_{K\sigma}^{n+1/2}| \le \|b\|_{\infty}$ and, owing to our assumption (4.2) on the meshes, $d_{\sigma} \le \zeta$ size($\mathcal{T}$) is small. Since $\psi^*(\xi) = 4(\cosh(\xi/2) - 1) = |\xi|^2/2 + \mathcal{O}(|\xi|^4) \le \frac{1}{2}|\xi|^2(1 + |\xi|)$ for small $\xi$ we can bound

$$ \psi^*(b_{K\sigma}^{n+1/2} d_{\sigma}) \le \frac{d_{\sigma}^2}{2} |b_{K\sigma}^{n+1/2}|^2 [1 + \|b\|_{\infty} d_{\sigma}]. $$

Since $b \in C^1(\bar{Q}_T)$ is Lipschitz and $d_{\sigma} \le \zeta$ size($\mathcal{T}$) we obtain

$$ \psi^*(b_{K\sigma}^{n+1/2} d_{\sigma}) \le \frac{d_{\sigma}^2}{2} \left[ (b(t,x) \cdot n_{K\sigma})^2 + C(\tau + \text{size}(\mathcal{T})) \right], \quad (t,x) \in [t^n, t^{n+1}] \times K $$

for some constant $C > 0$ depending on $b, \zeta$, but not on $\mathcal{T}, \tau$. Leveraging one last time the Lipschitz regularity of $b$ we see that

$$ \psi^*(b_{K\sigma}^{n+1/2} d_{\sigma}) \le \frac{d_{\sigma}^2}{2} \left[ \frac{1}{\tau m_K} \int_{t^n}^{t^{n+1}} \int_K (b(t,x) \cdot n_{K\sigma})^2 + C(\tau + \text{size}(\mathcal{T})) \right]. $$

Hence, using yet again the mesh isotropy, we get

$$ \begin{aligned} & \sum_{n,K} \tau \theta_K^{n+1/2} \sum_{\sigma \in \Sigma_K} \frac{m_{\sigma}}{2d_{\sigma}} \psi^*(b_{K\sigma}^n d_{\sigma}) \\ & \le \sum_{n,K} \theta_K^{n+1/2} \sum_{\sigma \in \Sigma_K} \left[ \int_{t^n}^{t^{n+1}} \int_K \frac{m_{\sigma} d_{\sigma}}{4m_K} (b(t,x) \cdot n_{K\sigma})^2 + \tau \frac{m_{\sigma} d_{\sigma}}{2} C(\tau + \text{size}(\mathcal{T})) \right] \\ & \stackrel{(4.3)-(4.1)}{\le} (1 + \varepsilon \tau) \int_{Q_T} \theta_{\tau, \tau} \frac{|b|^2}{2} + \frac{\zeta \|b\|_{\infty}^2}{2} \int_0^T \int_{\Omega \setminus \Omega_{\text{iso}}} \theta_{\tau, \tau} + \frac{Cd\zeta}{2} (\tau + \text{size}(\mathcal{T})) \|\theta_{\tau, \tau}\|_{L^1(Q_T)} \\ & \le \int_{Q_T} \theta_{\tau, \tau} \frac{|b|^2}{2} + r_{\tau, \tau} \end{aligned} $$

where the remainder

$$ r_{\tau, \tau} := \frac{\zeta \|b\|_{\infty}^2}{2} \int_0^T \int_{\Omega \setminus \Omega_{\text{iso}}} \theta_{\tau, \tau} + \left[ \frac{Cd\zeta}{2}(\tau + \text{size}(\mathcal{T})) + \frac{\|b\|_{\infty}^2}{2}\varepsilon\tau \right] \|\theta_{\tau, \tau}\|_{L^1(Q_T)}. $$

Therefore, combining the previous estimates with (4.37) we deduce that

$$ \sum_{n=0}^{N-1} \tau \mathcal{D}_{\psi}(\boldsymbol{\theta}_k^{n+1/2}) \geq \langle b, F_{\Sigma^k, \tau^k} \rangle - \int_{Q_T} \theta_{\tau^k, \tau^k} \frac{|b|^2}{2} - r_{\tau^k, \tau^k}. $$

Since $(\theta_{\tau^k, \tau^k})_k$ is converging weakly in $L^1$ it is also equiintegrable and $L^1$-bounded. Due to our assumption (4.4) that meas((0, T) × (Ω \ Ω$_{iso}$)) → 0 we see that $r_{\tau^k, \tau^k} \to 0$ as $k \to \infty$, and as a consequence

$$ \liminf_{k \to \infty} \sum_{n=0}^{N-1} \tau \mathcal{D}_{\psi}(\boldsymbol{\theta}_k^{n+1/2}) \geq \langle b, F \rangle - \int_{Q_T} \theta \frac{|b|^2}{2} \geq \int_{Q_T} B(\theta, F) - \eta. $$

Since $\eta > 0$ was arbitrary the claim follows.

If now (4.36) is infinite we can proceed in a similar fashion. For any fixed $M > 0$ large there is $b \in C^1(\bar{Q}_T; \mathbb{R}^d)$ such that

$$ \langle b, F \rangle - \int_{Q_T} \theta \frac{|b|^2}{2} \geq M. $$
Following the same line of thought as above we find

$$ \liminf_{k \to \infty} \sum_{n=0}^{N-1} \tau \mathcal{D}_\psi(\boldsymbol{\theta}_k^{n+1/2}) \geq M. $$

Since $M$ is arbitrary, this $\liminf$ is also infinite. We have just proven the following:

**Proposition 4.6.** Let $(\boldsymbol{\theta}_k^{n+1/2})_{n=0}^{N-1} \in (\mathbb{R}_+^T)^N$ and $(\boldsymbol{F}_k^{n+1/2})_{n=0}^{N-1} \in (\mathbb{F}_T)^N$ be given discrete curves associated with the mesh $\mathcal{T}^k$ and time step $\tau^k$. Suppose that their reconstruction $\theta_{\mathcal{T}^k, \tau^k} \in L^1(Q_T)$ and $F_{\Sigma^k, \tau^k} \in L^1(Q_T; \mathbb{R}^d)$ from (4.7) and (4.9) converge as

$$ \theta_{\mathcal{T}^k, \tau^k} \to \theta \quad \text{weakly in } L^1(Q_T), \quad F_{\Sigma^k, \tau^k} \to F \quad \text{weakly in } L^1(Q_T; \mathbb{R}^d). $$

Then,

$$ \liminf_{k \to \infty} \sum_{n=0}^{N-1} \tau \mathcal{D}_\psi(\boldsymbol{\theta}_k^{n+1/2}, \boldsymbol{F}_k^{n+1/2}) \geq \int_{Q_T} B(\theta, F). $$

**4.6. Convergence of the scheme.** It remains now to prove that the curve $\rho = \theta$, constructed in the previous sections as the limit of solutions to our discrete scheme, is actually an EDI solution.

**Theorem 4.7.** Let $(\rho_k^n)_{n=0}^N \in (\mathbb{R}_+^T)^N$ and $(\theta_k^{n+1/2})_{n=0}^{N-1} \in (\mathbb{R}_+^T)^N$ be the densities obtained as the unique solution of the scheme (2.9)-(2.10)-(2.11), associated with a mesh $\mathcal{T}^k$ and time step $\tau^k$ satisfying $\tau^k = o(d_{\Sigma^k}^{\min})$. Let $\rho_{\mathcal{T}^k, \tau^k} \in L^1(Q_T)$, $\theta_{\mathcal{T}^k, \tau^k} \in L^1(Q_T)$ and $\rho_{\mathcal{T}^k, \tau^k}^T \in L^1(\Omega)$ be the reconstructions defined in (4.7) and (4.8). Then, there exists $\rho \in L^1(Q_T) \cap C([0, T]; L_w^1(\Omega))$ such that

$$ \rho_{\mathcal{T}^k, \tau^k}, \theta_{\mathcal{T}^k, \tau^k} \to \rho \text{ in } L^1(Q_T) \quad \text{and} \quad \rho_{\mathcal{T}^k, \tau^k}^T \to \rho^T = \rho(T) \text{ weakly in } L^1(\Omega) $$

as $k \to +\infty$. Moreover, the limiting density $\rho$ is the unique EDI solution with initial datum $\rho^0$ in the sense of Definition 1.

**Proof.** First, notice that no subsequence is involved in the above statement, as a by-product of the uniqueness of EDI solutions, cf. Proposition A.2 in the Appendix. Showing compactness for $(\rho_{\mathcal{T}^k, \tau^k})_k$ and $(\theta_{\mathcal{T}^k, \tau^k})_k$, and the fact that any cluster point $\rho$ is an EDI solution then automatically gives the convergence of the whole sequence. So in what follows, we will not indicate when convergences hold up to a subsequence.

First of all, the convergence of the reconstructions is a consequence of Proposition 4.4. Summing Proposition 3.2 in time we get the discrete EDI estimate

$$ (4.38) \qquad \mathcal{E}_T(\rho_k^N) + \sum_n \tau \left\{ \mathcal{D}_\psi(\boldsymbol{\theta}^{n+1/2}, \boldsymbol{F}^{n+1/2}) + \mathcal{R}_\psi(\boldsymbol{\theta}^{n+1/2}) \right\} \leq \mathcal{E}_T(\rho_k^0). $$

Jensen's inequality and the definition (2.8) of $\rho_k^0$ gives

$$ \mathcal{H}(\rho_{\mathcal{T}^k}^0) = \sum_{K \in \mathcal{T}^k} m_K H(\rho_K^0) \leq \mathcal{H}(\rho^0), $$

whereas

$$ \left| \sum_{K \in \mathcal{T}^k} m_K \rho_K^0 V_K - \int_{\Omega} \rho^0 V \right| \leq \sum_{K \in \mathcal{T}^k} \int_K \rho^0(x) |V(x) - V(x_K)| \, d\,x \leq C \, \text{size}(\mathcal{T}^k) $$

thanks to the regularity of $V$ and (4.1). Therefore,

$$ \mathcal{E}_T(\rho_k^0) \leq \mathcal{E}(\rho^0) + \mathcal{O}(\text{size}(\mathcal{T}^k)). $$
We deduce from similar arguments that

$$ \left| \mathcal{E}_T(\rho_k^N) - \mathcal{E}(\rho_{T^k, \tau^k}^T) \right| = \left| \sum_{K \in \mathcal{T}^k} m_K \rho_K^N V_K - \int_{\Omega} \rho_{T^k, \tau^k}^T V \right| \le C \, \text{size}(\mathcal{T}^k). $$

Taking the above estimates into account in (4.38) leads to

$$ (4.39) \quad \mathcal{E}(\rho_{T^k, \tau^k}^T) + \sum_n \tau \left\{ \mathcal{D}_\psi(\boldsymbol{\theta}^{n+1/2}, \mathbf{F}^{n+1/2}) + \mathcal{R}_\psi(\boldsymbol{\theta}^{n+1/2}) \right\} \le \mathcal{E}(\rho^0) + C \, \text{size}(\mathcal{T}^k). $$

The first term on the left immediately passes to the limit by standard weak-L¹ lower semi-continuity of the convex functional $\mathcal{E}$ together with the weak convergence (4.12). Moreover, let us recall from (4.16) that $F_{\Sigma^k, \tau^k} \to F$ weakly in $L^1(Q_T)$. Propositions 4.5 and 4.6 thus allow to take the limit in the dissipation terms and conclude that

$$ \mathcal{E}(\rho^T) + 2 \int_{Q_T} \pi \left| \nabla \sqrt{\frac{\rho}{\pi}} \right|^2 + \int_{Q_T} B(\rho, F) \le \mathcal{E}(\rho^0). $$

It only remains to show that the pair $(\rho, F)$ solves the continuity equation with initial and terminal data $\rho^0$ and $\rho^T$ taken in the $C([0, T]; L^1(\Omega))$ sense. For this, take any test-function $\varphi \in C^2(\bar{Q}_T)$ and denote

$$ \varphi_K^n = \frac{1}{m_K} \int_K \varphi(t^n, \cdot). $$

For fixed $\mathcal{T}, \tau$, observe on the one hand that

$$ \begin{aligned} \langle \rho_{\mathcal{T}, \tau}, \partial_t \varphi \rangle &= \sum_n \sum_K m_K \rho_K^{n+1} (\varphi_K^{n+1} - \varphi^n) \\ &= \sum_n \sum_K m_K (\rho_K^n - \rho_K^{n+1}) \varphi_K^n + \sum_K m_K (\rho_K^N \varphi_K^N - \rho_K^0 \varphi_K^0) \\ &= \sum_n \sum_K m_K (\rho_K^n - \rho_K^{n+1}) \varphi_K^n + \int_{\Omega} \varphi(T, \cdot) \rho_{\mathcal{T}, \tau}^T - \int_{\Omega} \varphi(0, \cdot) \rho_{\mathcal{T}}^0, \end{aligned} $$

and on the other hand that

$$ \begin{aligned} \langle F_{\Sigma, \tau}, \nabla \varphi \rangle &= \frac{1}{2} \sum_{n,K} \sum_{\sigma \in \Sigma_K} m_{\sigma} d_{\sigma} F_{K\sigma}^{n+1/2} \frac{1}{\Delta_{\sigma}} \int_{t^n}^{t^{n+1}} \int_{\Delta_{\sigma}} \nabla \varphi \cdot n_{K\sigma} \\ &= \frac{1}{2} \sum_{n,K} \sum_{\sigma \in \Sigma_K} \tau m_{\sigma} d_{\sigma} F_{K\sigma}^{n+1/2} \frac{\varphi_L^n - \varphi_K^n}{d_{\sigma}} + \mathcal{O} \left( C_{\varphi}(\text{size}(\mathcal{T}) + \tau) \|F_{\Sigma, \tau}\|_{L^1(Q_T)} \right) \\ &= - \sum_{n,K} \sum_{\sigma \in \Sigma_K} \tau m_{\sigma} F_{K\sigma}^{n+1/2} \varphi_K^n + \mathcal{O} \left( C_{\varphi}(\text{size}(\mathcal{T}) + \tau) \|F_{\Sigma, \tau}\|_{L^1(Q_T)} \right), \end{aligned} $$

where $C_{\varphi}$ is a constant depending on $\|D_{t,x}^2 \varphi\|_{\infty}$ only. Hence by (2.9)

$$ \left| \langle \rho_{\mathcal{T}, \tau}, \partial_t \varphi \rangle + \langle F_{\Sigma, \tau}, \nabla \varphi \rangle - \int_{\Omega} \varphi(T, \cdot) \rho_{\mathcal{T}, \tau}^T + \int_{\Omega} \varphi(0, \cdot) \rho_{\mathcal{T}}^0 \right| = \mathcal{O} \left( C_{\varphi}(\text{size}(\mathcal{T}) + \tau) \|F_{\Sigma, \tau}\|_{L^1(Q_T)} \right). $$

Since $F_{\Sigma^k, \tau^k} \to F$ weakly in $L^1(Q_T; \mathbb{R}^d)$, since $\rho_{T^k, \tau^k} \to \rho$ in $L^1(Q_T)$, since $\rho_{T^k}^0 \to \rho^0$ in $L^1(\Omega)$ and since $\rho_{T^k, \tau^k}^T \to \rho^T$ weakly in $L^1(\Omega)$ we can take the limit to retrieve the weak formulation (1.7) of the continuity equation for $C^2(\bar{Q}_T)$ test functions, and therefore for all $\varphi \in C^1(\bar{Q}_T)$ by density. Finally, it is well-known [1] that any pair $(\rho, F)$ solving the
continuity equation with finite kinetic energy $\int_{Q_T} B(\rho, F) < +\infty$ is $L_w^1(\Omega)$-continuous in time with initial/terminal data $\rho^0, \rho^T$, and the proof is complete. $\square$

**Remark 4.8.** Note that, since Proposition 4.5 and 4.6 were established using weak convergence only, the strong convergence of the reconstructions in Proposition 4.4 is itself not strictly needed to prove that the limit is an EDI solution.

## 5. NUMERICAL IMPLEMENTATION AND RESULTS

In this section we describe the implementation of the scheme defined by (2.9)–(2.11). We also present some numerical tests confirming the second order accuracy, both in time and space.

**5.1. Nested Newton method.** Given $\rho^n$, computing $\rho^{n+1}$ requires solving the nonlinear system (2.9)–(2.11) with in particular $\rho_K^{n+1} = \Xi(\rho_K^n, \theta_K^{n+1/2})$. For practical numerical purposes, solving for $\theta^{n+1/2}$ as the primary variable requires solving a nonlinear scalar system in each cell in order to evaluate the function $\Xi(\rho_K^n, \cdot)$ itself and its derivatives. Since $\Xi$ is $C^1$, convex, and $|\Xi'(a, \cdot)| \le e$ this can be achieved efficiently with a Newton-Rhapson method. In order to limit the number of linear systems to be solved in practice, however, we follow here an alternative reparametrization strategy inspired from [7]. Recall from Section 2.2 that our interpolation $\Xi(a, \cdot)$ is defined in terms of the convex function $g = f^{-1}$ in (2.5), whose graph $\{y = g(x)\} \subset \mathbb{R} \times \mathbb{R}_+$ we choose to reparametrize as

$$ \mathbb{R} \ni s \mapsto (x(s), y(s)) := \begin{cases} (s + \xi, g(s + \xi)) & \text{if } s \le 0, \\ (f(\lambda s + g(\xi)), \lambda s + g(\xi)) & \text{else.} \end{cases} $$

Here $\xi > e^{-1}$ is an arbitrary cutoff threshold: for $s \le 0$ one runs the graph $y = g(x)$ at unit speed, while for $s > \xi$ one rather chooses to run the inverse graph $x = f(y)$ at speed $\lambda$. We impose $\lambda = g'(\xi)$, so that $x(\cdot)$ and $y(\cdot)$ are $C^1$ across $s = 0$. Note that, by construction, $\rho_K^{n+1} = \Xi(\rho_K^n, \theta_K^{n+1/2})$ if and only if there exists (a unique) $s_K \in \mathbb{R}$ such that

$$ (5.1) \qquad \rho_K^{n+1} = X(\rho_K^n; s_K) := \begin{cases} s_K & \text{if } \rho_K^n = 0 \\ \rho_K^n x(s_K) & \text{otherwise} \end{cases} $$

and

$$ \theta_K^{n+1/2} = Y(\rho_K^n; s_K) := \begin{cases} e^{-1}s_K & \text{if } \rho_K^n = 0 \\ \rho_K^n y(s_K) & \text{otherwise} \end{cases}. $$

At each time step, we first look for $s \in \mathbb{R}^T$ solving

$$ (5.2) \qquad m_K \frac{X(\rho_K^n; s_K) - \rho_K^n}{\tau} + \sum_{\sigma \in \Sigma_K} \frac{m_\sigma}{d_\sigma} \pi_\sigma \left( \frac{Y(\rho_K^n; s_K)}{\pi_K} - \frac{Y(\rho_L^n; s_L)}{\pi_L} \right) = 0, \quad \forall K \in \mathcal{T}, $$

by the Newton-Raphson method, and then update $\rho^{n+1}$ according to (5.1). Note that this reparametrization does not change the exact solution, which by Proposition 3.1 should satisfy $\rho_K^{n+1} \ge 0$ and $\theta_K^{n+1/2} > 0$ for all $K \in \mathcal{T}$.

**Remark 5.1.** The evaluation of $X(\rho_K^n; \cdot)$, $Y(\rho_K^n; \cdot)$ and their derivatives requires the solution of an inner nonlinear system only when $\rho_K^n > 0$ and $s_K \le 0$, which at convergence corresponds to $\rho_K^{n+1} \le \xi \rho_K^n$. In practice we verified numerically that, even choosing $\xi \approx e^{-1}$, the outer Newton method for (5.2) only require very few iterations. In that case, the inner Newton method is merely required in order to guarantee robustness of the scheme when dealing with solutions with steep gradients or vanishing densities.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Illustration of the refinement patterns used in the numerical tests; 2a is the base mesh used for both refinement patterns, while 2b and 2c are the second meshes for the two different refinement patterns used in the tests.</figcaption></figure>

5.2. Test-cases and numerical results. We set $\Omega = [0, 1]^2$ and consider the two refinement patterns illustrated in Figure 2. Note that only the subdivision refinement satisfies all the requirements on the mesh geometry from Section 4.1.

5.2.1. Convergence test. In order to investigate convergence and accuracy of our scheme we consider a test-case with gravitational potential $V(x) = -gx_1$ for $g \in \mathbb{R}_+$ and $x = (x_1, x_2) \in \Omega$, in which case the continuous model reduces to the following linear Fokker-Planck equation:

$$ \partial_t \rho - \Delta \rho - g \partial_{x_1} \rho = 0 \quad \text{in } Q_T. $$

An exact solution given by

$$ \rho(t, x) = \exp\left(-\alpha(t + \delta) + \frac{g}{2}x_1\right) \left(\pi \cos(\pi x_1) + \frac{g}{2} \sin(\pi x_1)\right) + \pi \exp\left(g\left(x_1 - \frac{1}{2}\right)\right) $$

with $\alpha = \pi^2 + g^2/4$ and $\delta > 0$. Note that the initial condition satisfies

$$ \rho^0(1, x_2) = \pi \exp\left(\frac{g}{2}\right) \left(1 - \exp(-\alpha\delta)\right) $$

for all $x_2 \in [0, 1]$, and in particular $\rho^0(1, x_2) = 0$ if $\delta = 0$, and $\rho(t, x_1, x_2) > 0$ for all $t > 0$ and $(x_1, x_2) \in \Omega$. In the following we fix $g = 1$. Tables 1–4 show the $L^1(Q_T)$, $L^2(Q_T)$, and $L^\infty(Q_T)$ errors computed for the reconstructed density $\rho_{\tau, \mathcal{T}}$ for the two values $\delta = 0.0001$, $\delta = 0$ and both refinement methods. We observe that the aimed second order in time and space convergence is practically achieved for the $L^1(Q_T)$ norm whatever the initialization and the mesh refinement strategy. In the presence of vacuum, when the initial profile $\rho^0$ partially vanishes along $x_1 = 1$ for $\delta = 0$, second order accuracy is lost for the $L^\infty(Q_T)$ norm, and to a lesser extent for the $L^2(Q_T)$ norm.

In figure 3 we plot the error in the dissipation balance

$$ \Delta(t^n) := \frac{1}{\tau} [\mathcal{E}_T(\rho^n) - \mathcal{E}_T(\rho^{n+1})) - \mathcal{D}_\psi(\theta^{n+1/2}) - \mathcal{D}_{\psi^*}(\theta^{n+1/2})] $$

as a function of time. Recall that $\Delta(t) \equiv 0$ is expected in the limit for continuous EDI solutions of the Fokker-Planck equation, so the smaller the numerical $\Delta$ the more accurately dissipation is captured by the scheme.
<table><caption>TABLE 1. $\delta = 0.001$, refinement by subdivision.</caption><thead><tr><th>τ</th><th>size(T)</th><th>L<sup>1</sup></th><th>rate</th><th>L<sup>2</sup></th><th>rate</th><th>L<sup>∞</sup></th><th>rate</th><th>ρ<sub>min</sub></th></tr></thead><tbody><tr><td>2.50e-02</td><td>3.06e-01</td><td>4.48e-03</td><td></td><td>1.04e-02</td><td></td><td>5.96e-02</td><td></td><td>6.68e-02</td></tr><tr><td>1.25e-02</td><td>1.53e-01</td><td>9.84e-04</td><td>2.19</td><td>2.44e-03</td><td>2.09</td><td>3.65e-02</td><td>7.08e-01</td><td>5.55e-02</td></tr><tr><td>6.25e-03</td><td>7.65e-02</td><td>2.39e-04</td><td>2.04</td><td>6.01e-04</td><td>2.03</td><td>1.46e-02</td><td>1.32</td><td>5.28e-02</td></tr><tr><td>3.13e-03</td><td>3.82e-02</td><td>5.95e-05</td><td>2.01</td><td>1.48e-04</td><td>2.02</td><td>4.94e-03</td><td>1.56</td><td>5.22e-02</td></tr><tr><td>1.56e-03</td><td>1.91e-02</td><td>1.48e-05</td><td>2.01</td><td>3.65e-05</td><td>2.02</td><td>1.42e-03</td><td>1.79</td><td>5.21e-02</td></tr><tr><td>7.81e-04</td><td>9.56e-03</td><td>3.67e-06</td><td>2.01</td><td>9.05e-06</td><td>2.01</td><td>3.38e-04</td><td>2.07</td><td>5.21e-02</td></tr></tbody></table>

<table><caption>TABLE 2. $\delta = 0.001$, refinement by repetition.</caption><thead><tr><th>τ</th><th>size(T)</th><th>L<sup>1</sup></th><th>rate</th><th>L<sup>2</sup></th><th>rate</th><th>L<sup>∞</sup></th><th>rate</th><th>ρ<sub>min</sub></th></tr></thead><tbody><tr><td>2.50e-02</td><td>3.06e-01</td><td>2.75e-03</td><td></td><td>7.64e-03</td><td></td><td>7.52e-02</td><td></td><td>6.68e-02</td></tr><tr><td>1.25e-02</td><td>1.53e-01</td><td>6.10e-04</td><td>2.17</td><td>1.87e-03</td><td>2.03</td><td>2.81e-02</td><td>1.42</td><td>5.55e-02</td></tr><tr><td>6.25e-03</td><td>7.65e-02</td><td>1.50e-04</td><td>2.03</td><td>4.85e-04</td><td>1.95</td><td>1.43e-02</td><td>9.76e-01</td><td>5.28e-02</td></tr><tr><td>3.13e-03</td><td>3.82e-02</td><td>3.75e-05</td><td>2.00</td><td>1.24e-04</td><td>1.97</td><td>5.29e-03</td><td>1.43</td><td>5.22e-02</td></tr><tr><td>1.56e-03</td><td>1.91e-02</td><td>9.36e-06</td><td>2.00</td><td>3.10e-05</td><td>2.00</td><td>1.64e-03</td><td>1.69</td><td>5.21e-02</td></tr><tr><td>7.81e-04</td><td>9.56e-03</td><td>2.34e-06</td><td>2.00</td><td>7.73e-06</td><td>2.00</td><td>4.31e-04</td><td>1.93</td><td>5.21e-02</td></tr></tbody></table>

<table><caption>TABLE 3. $\delta = 0$, refinement by subdivision.</caption><thead><tr><th>τ</th><th>size(T)</th><th>L<sup>1</sup></th><th>rate</th><th>L<sup>2</sup></th><th>rate</th><th>L<sup>∞</sup></th><th>rate</th><th>ρ<sub>min</sub></th></tr></thead><tbody><tr><td>2.50e-02</td><td>3.06e-01</td><td>4.60e-03</td><td></td><td>1.09e-02</td><td></td><td>8.49e-02</td><td></td><td>1.61e-02</td></tr><tr><td>1.25e-02</td><td>1.53e-01</td><td>1.05e-03</td><td>2.13</td><td>2.85e-03</td><td>1.93</td><td>6.43e-02</td><td>4.00e-01</td><td>4.06e-03</td></tr><tr><td>6.25e-03</td><td>7.65e-02</td><td>2.64e-04</td><td>1.99</td><td>8.26e-04</td><td>1.79</td><td>3.73e-02</td><td>7.85e-01</td><td>1.02e-03</td></tr><tr><td>3.13e-03</td><td>3.82e-02</td><td>6.68e-05</td><td>1.98</td><td>2.45e-04</td><td>1.75</td><td>2.00e-02</td><td>8.99e-01</td><td>2.55e-04</td></tr><tr><td>1.56e-03</td><td>1.91e-02</td><td>1.66e-05</td><td>2.01</td><td>7.31e-05</td><td>1.75</td><td>1.04e-02</td><td>9.48e-01</td><td>6.39e-05</td></tr><tr><td>7.81e-04</td><td>9.56e-03</td><td>4.13e-06</td><td>2.01</td><td>2.18e-05</td><td>1.74</td><td>5.29e-03</td><td>9.72e-01</td><td>1.60e-05</td></tr></tbody></table>

<table><caption>TABLE 4. $\delta = 0$, refinement by repetition.</caption><thead><tr><th>τ</th><th>size(T)</th><th>L<sup>1</sup></th><th>rate</th><th>L<sup>2</sup></th><th>rate</th><th>L<sup>∞</sup></th><th>rate</th><th>ρ<sub>min</sub></th></tr></thead><tbody><tr><td>2.50e-02</td><td>3.06e-01</td><td>2.84e-03</td><td></td><td>8.16e-03</td><td></td><td>8.17e-02</td><td></td><td>1.61e-02</td></tr><tr><td>1.25e-02</td><td>1.53e-01</td><td>6.74e-04</td><td>2.07</td><td>2.33e-03</td><td>1.81</td><td>5.55e-02</td><td>5.57e-01</td><td>4.06e-03</td></tr><tr><td>6.25e-03</td><td>7.65e-02</td><td>1.74e-04</td><td>1.96</td><td>7.23e-04</td><td>1.69</td><td>3.43e-02</td><td>6.94e-01</td><td>1.02e-03</td></tr><tr><td>3.13e-03</td><td>3.82e-02</td><td>4.43e-05</td><td>1.97</td><td>2.25e-04</td><td>1.69</td><td>1.91e-02</td><td>8.43e-01</td><td>2.55e-04</td></tr><tr><td>1.56e-03</td><td>1.91e-02</td><td>1.12e-05</td><td>1.99</td><td>6.92e-05</td><td>1.70</td><td>1.01e-02</td><td>9.15e-01</td><td>6.39e-05</td></tr><tr><td>7.81e-04</td><td>9.56e-03</td><td>2.82e-06</td><td>1.99</td><td>2.11e-05</td><td>1.71</td><td>5.23e-03</td><td>9.55e-01</td><td>1.60e-05</td></tr></tbody></table>

5.2.2. *A test case with a steep potential.* We define now in polar coordinates

$$V(x) = 5(1 - \exp(-r^2/\sigma^2))(1 - \cos^6(20r - \theta)), \qquad \rho^0(x) = \frac{1}{\sqrt{2\pi}\sigma} \exp(-r^2/\sigma^2),$$

where the origin $r = 0$ is set at the center of the domain $(1/2, 1/2)$ and $\sigma = 10^{-2}$. The solution for these data concentrates on a small area of the domain centered around the curve $20r = \theta$. The density profile at different times is displayed in figure 5. In figure 6 we show the time evolution of the minimum of the density and the number of outer Newton iterations.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Error in the energy balance $\Delta$, for $\delta = 0$. The labels refer to different meshes and time steps corresponding to the row numbers of Table 3 (left, refinement by subdivision) and 4 (right, refinement by repetition)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Dissipation rates $D_\psi(\theta^{n+1/2})$ (a), and $D_{\psi^*}(\theta^{n+1/2})$ (b), for $\delta = 0$ and using the finest mesh and time step, corresponding to the last row of Table 3 (left, refinement by subdivision) and 4 (right, refinement by repetition)</figcaption></figure>

## 6. CONCLUSION AND PROSPECTS

The numerical strategy we propose in this paper shows great promises, as its structure allows to carry out the rigorous numerical analysis. The convergence can even be established under reasonable assumptions on the mesh that might be relaxed even further. The scheme preserves positivity and is compatible with thermodynamics in the sense that free energy decays in time at a rate which accurately approaches the exact one. Moreover, our method is computationally efficient thanks to the local-in space character of the time extrapolation. In particular, the resolution of the nonlinear systems arising from our scheme does not seem to be computationally demanding in practice, at least in the test cases we have run.

We only validated our approach so far on the simple linear Fokker-Planck equation. The extension to nonlinear (possibly degenerate) parabolic equations (or even systems) with gradient structure remains to be done. We also plan to extend our strategy to the case of Poisson-Nernst-Planck systems, where the extrapolation step is no longer local but still computable at reasonable cost.
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Density evolution for the second test case (note that the color scale is renormalized in each picture).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Time evolution of density minimum $\min_{\tau, \tau'} \rho_{\tau, \tau'}(t, \cdot)$ and number of outer Newton iterations to solve the nonlinear system at each time step.</figcaption></figure>

## APPENDIX A. PROPERTIES OF EDI SOLUTIONS

We discuss here basic properties of EDI solutions needed for our purpose. The starting point is the following chain rule:
**Proposition A.1 (Chain rule).** Let $\Omega \subset \mathbb{R}^d$ be convex and $V \in C^2(\bar{\Omega})$. Take any $(\rho, F)$ satisfying the continuity equation (1.7) in time $[0, T]$, with moreover finite kinetic energy and Fisher information

$$ \int_{0}^{T} \mathcal{R}(\rho) + \int_{Q_T} B(\rho, F) < +\infty. $$

Then $t \mapsto \mathcal{E}(\rho^t)$ is absolutely continuous with distributional derivative

$$ (A.1) \qquad \frac{d}{dt} \mathcal{E}(\rho^t) = \int_{\Omega} F^t \cdot \nabla \log \left( \frac{\rho^t}{\pi} \right) \quad \in L^1(0, T). $$

As $\Omega$ is convex, one can directly apply results from [1] (see in particular §10.1.2.E, Proposition 10.3.18, and thm. 10.4.9 therein). Roughly speaking, the convexity of $\Omega$ combined with the regularity of $V$ guarantee that the relative entropy $\rho \mapsto \mathcal{E}(\rho) = \mathcal{H}(\rho|\pi)$ is $\lambda$-displacement convex for some $\lambda \in \mathbb{R}$, which then opens the way to the subdifferential calculus developed in [1]. The extension to non-smooth and non-convex Lipschitz domains of the above chain rule is an open problem up to our knowledge.

The following properties of EDI solutions is then an easy corollary:

**Proposition A.2.** Under the same assumptions, take in addition $\rho^0$ with finite energy $\mathcal{E}(\rho^0) < +\infty$. Then EDI solutions with initial datum $\rho^0$ are unique, solve the Fokker-Planck equation (1.1) at least in the distributional sense, and satisfy in fact Energy Dissipation Equality in the sense that $t \mapsto \mathcal{E}(\rho^t)$ is absolutely continuous with

$$ \frac{d}{dt} \mathcal{E}(\rho^t) = -\mathcal{R}(\rho^t) - \int_{\Omega} B(\rho^t, F^t) \quad \in L^1(0, T). $$

*Proof.* Let us first show that any EDI solution is a distributional solution, which amounts to proving that the flux driving the continuity equation (1.7) is $F = -\nabla\rho - \rho\nabla V = -\rho\nabla\log(\frac{\rho}{\pi})$. To this end we first note from (A.1) and Young's inequality that

$$ |\mathcal{E}(\rho^{t_1}) - \mathcal{E}(\rho^{t_0})| \le \frac{1}{2} \int_{t_0}^{t_1} \int_{\Omega} \left\{ \frac{|F|^2}{\rho} + \rho \left| \nabla \log \left( \frac{\rho}{\pi} \right) \right|^2 \right\} = \int_{t_0}^{t_1} \int_{\Omega} B(\rho, F) + \int_{t_0}^{t_1} \mathcal{R}(\rho) $$

for arbitrary subinterval $[t_0, t_1] \subset [0, T]$. As a consequence

$$ J(\rho, F; t_0, t_1) := \mathcal{E}(\rho^{t_1}) + \int_{t_0}^{t_1} \mathcal{R}(\rho^t) + \int_{t_0}^{t_1} \int_{\Omega} B(\rho^t, F^t) - \mathcal{E}(\rho^{t_0}) \ge 0 $$

is nonnegative. By additivity, if $(\rho, F)$ is an EDI solution in time $[0, T]$ - meaning $J(\rho, F; 0, T) \le 0$ - we have that

$$ 0 \le J(\rho, F; t_0, t_1) = J(\rho, F; 0, T) - J(\rho, F; 0, t_0) - J(\rho, F; t_1, T) \le 0 $$

and $(\rho, F)$ is therefore an EDI solution in any subinterval. With the absolute continuity from (A.1) this gives

$$ \begin{aligned} & -\frac{1}{2} \int_{\Omega} \frac{|F^t|^2}{\rho^t} - \frac{1}{2} \int_{\Omega} \rho^t \left| \nabla \log \left( \frac{\rho^t}{\pi} \right) \right| \\ & \le \int_{\Omega} F^t \cdot \nabla \log \left( \frac{\rho^t}{\pi} \right) = \frac{d}{dt} \mathcal{E}(\rho^t) \stackrel{\text{EDI}}{\le} - \int_{\Omega} B(\rho^t, F^t) - \mathcal{R}(\rho^t) \\ & = -\frac{1}{2} \int_{\Omega} \frac{|F^t|^2}{\rho^t} - \frac{1}{2} \int_{\Omega} \rho^t \left| \nabla \log \left( \frac{\rho^t}{\pi} \right) \right|. \end{aligned} $$
This forces equality in Young's inequality, hence $F^t = -\rho^t \nabla \log(\frac{\rho^t}{\pi})$ and $(\rho, F)$ is indeed a distributional solution. This also gives equality in EDI, and $(\rho, F)$ is in fact an EDE solution as in our statement.

For the uniqueness we implement the so-called *Gigli's trick* [30]. Let $(\rho_1, F_1), (\rho_2, F_2)$ be two EDI solutions with common initial datum $\rho^0$ and set

$$ \tilde{\rho} := \frac{1}{2}(\rho_1 + \rho_2) \quad \text{and} \quad \tilde{F} := \frac{1}{2}(F_1 + F_2). $$

By linearity $(\tilde{\rho}, \tilde{F})$ still solves the continuity equation with initial datum $\tilde{\rho}^0 = \rho^0$. Fix any $\tau \in (0, T]$ and recall from the first part of the proof that any EDI solution in $[0, T]$ is also an EDI solution in $[0, \tau]$. Summing the two EDI inequalities for $i = 1, 2$ and leveraging the joint convexity of the Fisher information and Benamou-Brenier functionals, we see that

$$ \begin{aligned} & \frac{1}{2}(\mathcal{E}(\rho_1^\tau) + \mathcal{E}(\rho_2^\tau)) + \int_0^\tau \mathcal{R}(\tilde{\rho}) + \int_{Q_\tau} B(\tilde{\rho}, \tilde{F}) \\ & \le \frac{1}{2}(\mathcal{E}(\rho_1^\tau) + \mathcal{E}(\rho_2^\tau)) + \int_0^\tau \frac{1}{2}[\mathcal{R}(\rho_1) + \mathcal{R}(\rho_2)] + \int_{Q_\tau} \frac{1}{2}[B(\rho_1, F_1) + B(\rho_2, F_2)] \\ & \le \frac{1}{2}(\mathcal{E}(\rho_1^0) + \mathcal{E}(\rho_2^0)) = \mathcal{E}(\rho^0) = \mathcal{E}(\tilde{\rho}^0). \end{aligned} $$

This shows in particular that $(\tilde{\rho}, \tilde{F})$ has finite kinetic energy and Fisher information in $[0, \tau]$. Owing to (A.1) and Young's inequality we see that $\mathcal{E}(\tilde{\rho}^0) \le \mathcal{E}(\tilde{\rho}^\tau) + \int_0^\tau \mathcal{R}(\tilde{\rho}) + \int_{Q_\tau} B(\tilde{\rho}, \tilde{F})$, which substituted into the above right-hand side yields

$$ \mathcal{E}\left(\frac{\rho_1^\tau + \rho_2^\tau}{2}\right) = \mathcal{E}(\tilde{\rho}^\tau) \ge \frac{1}{2}(\mathcal{E}(\rho_1^\tau) + \mathcal{E}(\rho_2^\tau)). $$

By strict convexity of $\mathcal{E}$ this implies $\rho_1^\tau = \rho_2^\tau$, and since $\tau \in (0, T)$ was arbitrary the proof is complete. $\square$

## APPENDIX B. AN AUBIN-LIONS-MOUSSA COMPENSATION-COMPACTNESS RESULT

The technical statement below is a kind of compensation-compactness argument, allowing to pass to the limit for the product of two weakly converging sequences.

**Proposition B.1.** Let $\Omega \subset \mathbb{R}^d$ be Lipschitz and bounded, and consider two sequences $\rho_k, f_k$ such that

(i) $\rho_k$ is bounded in $L^1(Q_T)$ and $f_k$ is bounded in $L^\infty(Q_T)$

(ii) $\rho_k$ is $L^1(Q_T)$-equiintegrable

(iii) for some fixed modulus of continuity $\eta(\cdot)$ and all $\omega \subset \subset \Omega$, there holds

$$ (B.1) \quad \sup_k \int_0^T \int_\omega |f_k(t, x) - f_k(t, x + h)| \le \eta(|h|) \quad \text{for } |h| \le \operatorname{dist}(\omega, \partial\Omega) $$

(iv) $\partial_t \rho_k$ is bounded in $M(0, T; (W^{m,q}(\Omega))')$ for some fixed $m \in \mathbb{N}, q \ge 1$

Then, up to extraction of a subsequence, $\rho_k \rightharpoonup \rho$ weakly in $L^1(Q_T)$, $f_k \stackrel{*}{\rightharpoonup} f$ weakly-* in $L^\infty(Q_T)$, and $\rho_k f_k \rightharpoonup \rho f$ in $M(Q_T)$ in the sense that

$$ \int_{Q_T} \rho_k f_k \varphi \to \int_{Q_T} \rho f \varphi, \quad \forall \varphi \in C(\bar{Q}_T). $$
We stress that this is merely a variant on [50, Proposition 1], with however the subtle difference that the space compactness is obtained therein via suitable $W^{1,p}(\Omega)$ bounds, whereas we use here the more versatile difference quotient estimate (iii). This has two main advantages: first, this covers the case where space compactness is obtained via discrete difference quotients, as is typical for finite volume schemes (see our practical application for the proof of Proposition 4.4). Second, and more importantly, this allows more general scenarios since in [50] some restriction is imposed between the (spatial) Sobolev exponent $p \in [1, d]$ on $f$ and the limitation $\alpha \in [1, p^*)$ on some dual $L^{\alpha'}(\Omega)$ estimates on $\rho_k$. This will be circumvented in our particular setting via the equiintegrability assumption (ii) on $\rho_k$, but at the expense of our $L^\infty$ bounds on $f_k$.

*Proof.* Assume first that $\Omega = \mathbb{R}^d$ and that our space compactness (B.1) holds with $\omega = \mathbb{R}^d$. We reproduce below the proof of [50, Proposition 1] almost verbatim, and only point out the main differences. Pick a mollifying sequence $\zeta_n(x) = n^d\zeta(nx)$ (acting in space only and supported in $B_{1/n}$), and write for any test-function $\varphi$

$$ \begin{aligned} \int_0^T \int_{\mathbb{R}^d} (\rho f - \rho_k f_k)\varphi &= \int_0^T \int_{\mathbb{R}^d} [\rho f - (\rho * \zeta_n)f]\varphi \\ &\quad + \int_0^T \int_{\mathbb{R}^d} [(\rho * \zeta_n)f - (\rho_k * \zeta_n)f_k]\varphi \\ &\quad + \int_0^T \int_{\mathbb{R}^d} [(\rho_k * \zeta_n)f_k - (\rho_k f_k) * \zeta_n]\varphi \\ &\quad + \int_0^T \int_{\mathbb{R}^d} [(\rho_k f_k) * \zeta_n - \rho_k f_k]\varphi \\ &=: I_1(n) + I_2(k, n) + I_3(k, n) + I_4(k, n). \end{aligned} $$

Arguing as in [50] one shows without too much trouble that $I_1(n) \to 0$ as $n \to \infty$, that $I_2(k, n) \to 0$ for fixed $n$ as $k \to \infty$, and that $I_4(k, n) \to 0$ uniformly in $k$ as $n \to \infty$. The main difference lies here in the Friedrich commutator $I_3$, which we handle now with care. To this end we will show that

$$ S_{k,n}(t,x) := [(\rho_k * \zeta_n)f_k - (\rho_k f_k) * \zeta_n](t,x) $$

converges strongly to zero in $L^1((0,T) \times \mathbb{R}^d)$ as $n \to \infty$, uniformly in $k$. We first observe that, by Fubini's theorem,

$$ \|S_{k,n}\|_{L^1((0,T) \times \mathbb{R}^d)} \le \int_{B_{1/n}} \int_0^T \int_{\mathbb{R}^d} |\rho_k(t, x-y)| |f_k(t, x) - f_k(t, x-y)| \zeta_n(y) dx dt dy. $$

Take $r(n) \to +\infty$ as $n \to \infty$, and let $E_{k,n}^y = \{(t, x) \in (0, T) \times \mathbb{R}^d : |\rho_k(t, x-y)| > r(n)\}$. Owing to our $L^1$ bound (i) we have that meas$(E_{k,n}^y) \le \frac{1}{r(n)} \int_{Q_T} \rho_k \le \frac{C}{r(n)} \to 0$ uniformly in $k$ as $n \to \infty$. Whence by the equiintegrability assumption (ii)

$$ \begin{aligned} & \int_{E_{k,n}^y} |\rho_k(t, x-y)| |f_k(t, x) - f_k(t, x-y)| dx dt \\ & \le 2 \|f_k\|_{L^\infty((0,T) \times \mathbb{R}^d)} \int_{E_{k,n}^y} |\rho_k(t, z)| dz dt \underset{n \to \infty}{\longrightarrow} 0, \quad \text{uniformly in } y, k, \end{aligned} $$
and integrating in $y \in B_{1/n}$

$$ (B.2) \quad \int_{B_{1/n}} \int_{E_{k,n}^y} |\rho_k(t, x-y)| |f_k(t, x) - f_k(t, x-y)| \zeta_n(y) dx dt dy \xrightarrow{n \to \infty} 0, \quad \text{uniformly in } k. $$

In $(E_{k,n}^y)^c$ we have by definition $|\rho_k(t, x-y)| \le r(n)$ and we use instead the space compactness (B.1) (recalling also that we can take $\omega = \mathbb{R}^d$ at this stage) to estimate

$$ \begin{aligned} & \int_{(E_{k,n}^y)^c} |\rho_k(t, x-y)| |f_k(t, x) - f_k(t, x-y)| dx dt \\ & \le \int_{(E_{k,n}^y)^c} r(n) |f_k(t, x) - f_k(t, x-y)| dx dt \\ & \le r(n) \sup_k \int_{(0,T) \times \mathbb{R}^d} |f_k(t, x) - f_k(t, x-y)| dx dt \\ & \le r(n) \eta(|y|). \end{aligned} $$

Choosing the speed $r(n) = o(\frac{1}{\eta(1/n)})$ as $n \to \infty$ and integrating in $y$, we end up with

$$ (B.3) \quad \int_{B_{1/n}} \int_{(E_{k,n}^y)^c} |\rho_k(t, x-y)| |f_k(t, x) - f_k(t, x-y)| \zeta_n(y) dx dt dy \\ \le \int_{B_{1/n}} r(n) \eta(|y|) \zeta_n(y) dy \le r(n) \eta(1/n) \xrightarrow{n \to \infty} 0, \quad \text{uniformly in } k. $$

Gathering (B.2)–(B.3) gives $\|S_{k,n}\|_{L^1} \to 0$ and therefore

$$ |I_3(k, n)| \le \|S_{k,n}\|_{L^1} \|\varphi\|_{L^\infty} \xrightarrow{n \to \infty} 0, \quad \text{uniformly in } k $$

at least in the whole space $\Omega = \mathbb{R}^d$. The rest of the proof is then identical to [50, Proposition 1] and we omit the details for the sake of brevity.

Coming back now to the case of bounded domains $\Omega \subset \mathbb{R}^d$, we can apply the same localization argument from [50, Proposition 3]: Take an exhausting sequence of compact sets $K_l \subset \subset \Omega$ and a sequence of bump functions $0 \le \chi_l(x) \le 1$ such that $\chi \equiv 1$ on $K_l$, with meas($\Omega \setminus K_l$) $\le 1/l$ as $l \to +\infty$. Extending $\rho_k, f_k$ by zero outside of $\Omega$, it is easy to check that the sequences $\{\chi_l \rho_k\}_k, \{\chi_l f_k\}_k$ satisfy the assumptions in the previous step for fixed $l$, hence $\chi_l^2 \rho_k f_k \to \chi_l^2 \rho f$ in the sense of measures as $k \to \infty$. Writing $E_l = \Omega \setminus K_l$, we get

$$ \begin{aligned} \left| \int_{Q_T} (\rho_k f_k - \rho f) \varphi \right| &= \left| \int_{Q_T} (\chi_l^2 \rho_k f_k - \chi_l^2 \rho f) \varphi \right| + \int_{Q_T} (1 - \chi_l^2) \rho_k f_k \varphi - \int_{Q_T} (1 - \chi_l^2) \rho f \varphi \right| \\ &\le \left| \int_{Q_T} (\chi_l^2 \rho_k f_k - \chi_l^2 \rho f) \varphi \right| + \int_{Q_T} (1 - \chi_l^2) |\rho_k f_k \varphi| + \int_{Q_T} (1 - \chi_l^2) |\rho f \varphi| \\ &\le \left| \int_{Q_T} (\chi_l^2 \rho_k f_k - \chi_l^2 \rho f) \varphi \right| + \|\varphi\|_{L^\infty(Q_T)} \int_0^T \int_{E_l} |\rho_k f_k| + |\rho f| \cdot \end{aligned} $$

Pick an arbitrary $\varepsilon > 0$. With our assumptions (i)–(ii) it is easy to check that

$$ \sup_k \int_0^T \int_E |\rho_k f_k| \to 0 \quad \text{as meas}(E) \to 0, \quad E \subset \Omega. $$
Owing to meas$(E_l) \le 1/l \to 0$ and $\rho f \in L^1(Q_T)$, we can first choose $l = l_0$ large enough so that the second term in the r.h.s. is less than $\varepsilon/2$. For this fixed $l = l_0$ the first term can also be made smaller than $\varepsilon/2$ if $k \ge k_0$ is large enough, and the claim finally follows. □

**Acknowledgements** This work was supported in part by the Labex CEMPI (ANR-11-LABX-0007-01) and by the project MATHSOUT of the PEPR Mathematics in Interaction (ANR-23-EXMA-0010) funded by the French National Research Agency. LM was supported by Fundação para a Ciência e a Tecnologia through grants 10.54499/UIDB/00208/2020 and 10.54499/2020.00162.CEECIND/CP1595/CT0008.

## REFERENCES

[1] L. Ambrosio, N. Gigli, and G. Savaré. Gradient flows in metric spaces and in the space of probability measures. Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, second edition, 2008.

[2] L. Ambrosio and S. Serfaty. A gradient flow approach to an evolution problem arising in superconductivity. *Comm. Pure Appl. Math.*, 61(11):1495-1539, 2008.

[3] B. Andreianov, C. Cancès, and A. Moussa. A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs. *J. Funct. Anal.*, 273(12):3633-3670, 2017.

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

[5] J.-D. Benamou, G. Carlier, and M. Laborde. An augmented Lagrangian approach to Wasserstein gradient flows and applications. In *Gradient flows: from theory to application*, volume 54 of *ESAIM Proc. Surveys*, pages 1-17. EDP Sci., Les Ulis, 2016.

[6] A. Blanchet. A gradient flow approach to the Keller-Segel systems. RIMS Kokyuroku's lecture notes, vol. 1837, pp. 52-73, June 2013.

[7] K. Brenner and C. Cancès. Improving Newton's method performance by parametrization: The case of the Richards equation. *SIAM J. Numer. Anal.*, 55(4):1760-1785, 2017.

[8] C. Cancès, C. Chainais-Hillairet, B. Merlet, F. Raimondi, and J. Venel. Mathematical analysis of a thermodynamically consistent reduced model for iron corrosion. *Z. Angew. Math. Phys.*, 74(96), 2023.

[9] C. Cancès and C. Guichard. Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure. *Found. Comput. Math.*, 17(6):1525-1584, 2017.

[10] C. Cancès, D. Matthes, and F. Nabet. A two-phase two-fluxes degenerate Cahn-Hilliard model as constrained Wasserstein gradient flow. *Arch. Ration. Mech. Anal.*, 233(2):837-866, 2019.

[11] C. Cancès and J. Venel. On the square-root approximation finite volume scheme for nonlinear drift-diffusion equations. *Comptes Rendus. Mathématique*, 361:525-558, 2023.

[12] C. Cancès, T. O. Gallouët, and G. Todeschi. A variational finite volume scheme for Wasserstein gradient flows. *Numer. Math.*, 146(3):437-480, 2020.

[13] C. Cancès, D. Matthes, F. Nabet, and E.-M. Rott. Finite elements for Wasserstein $W_p$ gradient flows. HAL: hal-03719189, 2024.

[14] J. A. Carrillo, K. Craig, L. Wang, and C. Wei. Primal dual methods for Wasserstein gradient flows. *Found. Comput. Math.*, 2021. Online first.

[15] J. A. Carrillo, B. Düring, D. Matthes, and M. S. McCormick. A Lagrangian scheme for the solution of nonlinear diffusion equations using moving simplex meshes. *J. Sci. Comput.*, 73(3):1463-1499, 2018.

[16] J.-B. Casteras and L. Monsaingeon. Hidden dissipation and convexity for Kimura equations. *SIAM J. Math. Anal.*, 55(6):7361-7398, 2023.

[17] J. Droniou. A density result in Sobolev spaces. *J. Math. Pures Appl.*, 81(7):697-714, 2002.

[18] J. Droniou and N. Nataraj. Improved $L^2$ estimate for gradient schemes and super-convergence of the TPFA finite volume scheme. *IMA J. Numer. Anal.*, 38(3):1254-1293, 2018.

[19] A. Esposito, F. S. Patacchini, A. Schlichting, and D. Slepčev. Nonlocal-interaction equation on graphs: Gradient flow structure and continuum limit. *Arch. Ration. Mech. Anal.*, 240:699-760, 2021.

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

[21] R. Eymard, T. Gallouët, and R. Herbin. Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces. *IMA J. Numer. Anal.*, 30(4):1009-1043, 2010.
[22] C. Falcó, R.E. Baker, and J.A. Carrillo. A local continuum model of cell-cell adhesion. *SIAM J. Appl. Math*, 84(3):S17-S42, 2023.

[23] D. Forkert, J. Maas, and L. Portinale. Evolutionary Γ-convergence of entropic gradient flow structures for Fokker-Planck equations in multiple dimensions. *SIAM J. Math. Anal.*, 54(4):4297-4333, 2022.

[24] T. Frenzel and M. Liero. Effective diffusion in thin structures via generalized gradient systems and EDP-convergence. *Discr. Cont. Dyn. Syst. S*, 14(1):395-425, 2021.

[25] A. T. Jr. Fromhold and E. L. Cook. Diffusion currents in large electric fields for discrete lattices. *J. Appl. Phys.*, 38(4):1546-1553, 1967.

[26] D. Furihata and T. Matsuo. *Discrete variational derivative method. A structure-preserving numerical method for partial differential equations*. Boca Raton, FL: CRC Press, 2011.

[27] A. Galichon. *Optimal transport methods in economics*. Princeton University Press, 2018.

[28] T. O. Gallouët, Q. Mérigot, and A. Natale. Convergence of a Lagrangian discretization for barotropic fluids and porous media flow. *SIAM J. Math. Anal.*, 54(3):2990-3018, 2022.

[29] T. O. Gallouët, A. Natale, and G. Todeschi. From geodesic extrapolation to a variational BDF2 scheme for Wasserstein gradient flows. *Math. Comp.*, 2024.

[30] N. Gigli. On the heat flow on metric measure spaces: existence, uniqueness and stability. *Calc. Var. Partial Differential Equations*, 39(1):101-120, 2010.

[31] P. Gladbach, E. Kopfer, and J. Maas. Scaling limits of discrete optimal transport. *SIAM J. Math. Anal.*, 52(3):2759-2802, 2020.

[32] M. Heida. Convergences of the squareroot approximation scheme to the Fokker-Planck operator. *Math. Models Methods Appl. Sci.*, 28(13):2599-2635, 2018.

[33] A. Hraivoronska, A. Schlichting, and O. Tse. Variational convergence of the Scharfetter-Gummel scheme to the aggregation-diffusion equation and vanishing diffusion limit. *arXiv:2306.02226*, 2023.

[34] R. Jordan, D. Kinderlehrer, and F. Otto. The variational formulation of the Fokker-Planck equation. *SIAM J. Math. Anal.*, 29(1):1-17, 1998.

[35] A. Jüngel, U. Stefanelli, and L. Trussardi. Two structure-preserving time discretizations for gradient flows. *Appl. Math. Opt.*, 80:733-764, 2019.

[36] D. Kinderlehrer, L. Monsaingeon, and X. Xu. A Wasserstein gradient flow approach to Poisson-Nernst-Planck equations. *ESAIM Control Optim. Calc. Var.*, 23(1):137-164, 2017.

[37] H. Lavenant. Unconditional convergence for discretizations of dynamical optimal transport. *Math. Comp.*, 90(328):739-786, 2021.

[38] H. Leclerc, Q. Mérigot, F. Santambrogio, and F. Stra. Lagrangian discretization of crowd motion and linear diffusion. *SIAM J. Numer. Anal.*, 58(4):2093-2118, 2020.

[39] G. Legendre and G. Turinici. Second-order in time schemes for gradient flows in Wasserstein and geodesic metric spaces. *Comptes Rendus. Mathematique*, 355(3):345-353, 2017.

[40] W. Li, J. Lu, and L. Wang. Fisher information regularization schemes for Wasserstein gradient flows. *J. Comput. Phys.*, 416:109449, 2020.

[41] H. C. Lie, K. Fackeldey, and M. Weber. A square root approximation of transition rates for a Markov state model. *SIAM J. Matrix Anal. Appl.*, 34(2):738-756, 2013.

[42] M. Liero, A. Mielke, M. A. Peletier, and D. R. M. Renger. On microscopic origins of generalized gradient structures. *Discr. Cont. Dyn. Syst. S*, 10(1):1-35, 2017.

[43] C. Liu, C. Wang, S. M. Wise, X. Yue, and Zhou S. A second order accurate numerical method for the Poisson-Nernst-Planck system in the energetic variational formulation. *arXiv:2208.06123*, 2022.

[44] J. Maas. Gradient flows of the entropy for finite Markov chains. *J. Funct. Anal.*, 261(8):2250-2292, 2011.

[45] D. Matthes and H. Osberger. Convergence of a variational Lagrangian scheme for a nonlinear drift diffusion equation. *ESAIM Math. Model. Numer. Anal.*, 48(3):697-726, 2014.

[46] D. Matthes, E.-M. Rott, G. Savaré, and A. Schlichting. A structure preserving discretization for the Derrida-Lebowitz-Speer-Spohn equation based on diffusive transport. *arXiv:2312.13284*, 2023.

[47] B. Maury, A. Roudneff-Chupin, and F. Santambrogio. A macroscopic crowd motion model of gradient flow type. *Math. Models Methods Appl. Sci.*, 20(10):1787-1821, 2010.

[48] A. Mielke, R. I. A. Patterson, and M. A. Peletier. Non-equilibrium thermodynamical principles for chemical reactions with mass-action kinetics. *SIAM J. Appl. Math.*, 77(4):1562-1585, 2017.

[49] A. Mielke, M. A. Peletier, and D. R. M. Renger. On the relation between gradient flows and the large-deviation principle, with applications to Markov chains and diffusion. *Potential Anal.*, 41(4):1293-1327, 2014.

[50] A. Moussa. Some variants of the classical Aubin-Lions lemma. *J. Evol. Equ.*, 16:65-93, 2016.
[51] A. Natale. Gradient flows of interacting Laguerre cells as discrete porous media flows. *arXiv:2304.05069*, 2023.

[52] A. Natale and G. Todeschi. TPFA finite volume approximation of Wasserstein gradient flows. In R. Klöfkorn, E. Keilegavlen, F. A. Radu, and J. Fuhrmann, editors, *Finite Volumes for Complex Applications IX – Methods, Theoretical Aspects, Examples*, volume 323 of *Springer Proceedings in Mathematics & Statistics*, pages 193–202, 2020.

[53] F. Otto. The geometry of dissipative evolution equations: the porous medium equation. *Comm. Partial Differential Equations*, 26(1-2):101–174, 2001.

[54] M. A. Peletier, R. Rossi, G. Savaré, and O. Tse. Jump processes as generalized gradient flows. *Calc. Var. Partial Differential Equations*, 61:33, 2022.

[55] M. A. Peletier and A. Schlichting. Cosh gradient systems and tilting. *Nonlinear Anal.*, 231:113094, 2023.

[56] J. W. Portegies and M. A. Peletier. Well-posedness of a parabolic moving-boundary problem in the setting of Wasserstein gradient flows. *Interfaces Free Bound.*, 12(2):121–150, 2010.

[57] F. Santambrogio. *Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling*. Progress in Nonlinear Differential Equations and Their Applications 87. Birkhäuser Basel, 1 edition, 2015.

[58] U. Stefanelli. A new minimizing-movements scheme for curves of maximal slope. *ESAIM Control Optim. Calc. Var.*, 28:59, 2022.

INRIA, UNIV. LILLE, CNRS, UMR 8524 - LABORATOIRE PAUL PAINLEVÉ, F-59000 LILLE, FRANCE.
*Email address: clement.cances@inria.fr, andrea.natale@inria.fr*

GROUP OF MATHEMATICAL PHYSICS, DEPARTAMENTO DE MATEMÁTICA, INSTITUTO SUPERIOR TÉCNICO,
AV. ROVISCO PAIS, 1049-001 LISBOA, PORTUGAL

IECL, FACULTÉ DES SCIENCES ET TECHNOLOGIE, CAMPUS, BOULEVARD DES AIGUILLETTES, BP 70239,
54506 VANDŒUVRE-LÈS-NANCY, FRANCE.

*Email address: leonard.monsaingeon@ciencias.ulisboa.pt*