---
abstract: |
    The dynamical formulation of optimal transport, also known as Benamou-Brenier formulation or Computational Fluid Dynamic formulation, amounts to write the optimal transport problem as the optimization of a convex functional under a PDE constraint, and can handle *a priori* a vast class of cost functions and geometries. Several disretizations of this problem have been proposed, leading to computations on flat spaces as well as Riemannian manifolds, with extensions to mean field games and gradient flows in the Wasserstein space.

    In this article, we provide a framework which guarantees convergence under mesh refinement of the solutions of the space-time discretized problems to the one of the infinite-dimensional one for quadratic optimal transport. The convergence holds without condition on the ratio between spatial and temporal step sizes, and can handle arbitrary positive measures as input, while the underlying space can be a Riemannian manifold. Both the finite volume discretization proposed by Gladbach, Kopfer and Maas, as well as the discretization over triangulations of surfaces studied by the present author in collaboration with Claici, Chien and Solomon fit in this framework.
address: 'Department of Mathematics, University of British Columbia, Vancouver BC Canada V6T 1Z2'
author:
- Hugo Lavenant
bibliography:
- 'bibliography.bib'
title: Unconditional convergence for discretizations of dynamical optimal transport
---

Introduction
============

In this article, we study a problem of calculus of variations of the form $$\label{equation_problem_unformal_v}
\min_{\rho, {\mathbf{v}}} \left\{ \iint \frac{1}{2} \rho |{\mathbf{v}}|^2 \, {\mathrm{d}}t \, {\mathrm{d}}x  \ \text{ such that } \ {\partial}_t \rho + \nabla \cdot ( \rho {\mathbf{v}}) = 0 \right\},$$ where the unknowns are $\rho$ a space-time dependent density and ${\mathbf{v}}$ a space-time dependent velocity field. These two unknowns are linked by the continuity equation ${\partial}_t \rho + \nabla \cdot ( \rho {\mathbf{v}}) = 0$, supplemented with no-flux boundary conditions, which says that the mass with density $\rho$ is advected by the velocity field ${\mathbf{v}}$. Provided the velocity field is smooth enough (which is not always the case for the optimizer), $\rho$ is uniquely determined by its temporal boundary values and ${\mathbf{v}}$, hence ${\mathbf{v}}$ could be though as a control. The functional which is optimized is the *action*, that is the space-time integral of the density of kinetic energy. The initial and final temporal values of $\rho$ can be either fixed, or a penalization depending on them can be added in the objective functional.

Instances of this problem appear in at least two situations. If the initial and final values $\rho_0$ and $\rho_1$ of $\rho$ are fixed, the problem amounts to compute a geodesic in the Wasserstein space (see [@Benamou2000] or [@SantambrogioOTAM Chapter 5]), also known as the displacement interpolation (or McCann’s interpolation [@Mccann1997]) between the initial and the final value of $\rho$. In short, the optimal $\rho$ is for some geometry the shortest curve joining its endpoints $\rho_0$ and $\rho_1$, while the value of the problem is the squared Wasserstein distance between $\rho_0$ and $\rho_1$. If the final value is penalized rather than fixed, reads as one step of the minimizing movement scheme (also known as the JKO scheme after [@Jordan1998]) for gradient flows in the Wasserstein space (see for instance [@AGS] or [@SantambrogioOTAM Chapter 8]). On the other hand, one could also add a running cost depending on the density to the action, and then the problem would read as the variational formulation of an instance of Mean Field Games [@Benamou2016]: roughly, $\rho$ describes the evolution in time of a density of a crowd of agents trying to reach a final destination with minimum effort while avoiding other agents.

In this article, we will mainly focus on the case where the initial and final values of $\rho$ are fixed (that is the computation of a geodesic in the Wasserstein space) because this case concentrates most of the difficulties. In the last section, we will also deal with a final penalization of the density which does not contain additional technicalities, and we will comment on the issues when on tries to add a running cost depending on the density.

Since the seminal work of Benamou and Brenier [@Benamou2000], it is understood that can be recast as a convex problem: the idea is to introduce a new unknown ${\mathbf{m}}= \rho {\mathbf{v}}$ as the continuity equation becomes the linear constraint ${\partial}_t \rho + \nabla \cdot {\mathbf{m}}= 0$, while the density of kinetic energy $$\frac{1}{2} \rho |{\mathbf{v}}|^2 = \frac{|{\mathbf{m}}|^2}{2 \rho}$$ is a jointly convex function of $\rho$ and ${\mathbf{m}}$. Hence, we will rather consider the problem $$\label{equation_problem_unformal}
\min_{\rho, {\mathbf{m}}} \left\{ \iint \frac{|{\mathbf{m}}|^2}{2 \rho} \, {\mathrm{d}}t \, {\mathrm{d}}x  \ \text{ such that } \ {\partial}_t \rho + \nabla \cdot {\mathbf{m}}= 0 \right\},$$ whose rigorous formulation is given below in Section \[section\_framework\]. With this convex reformulation, one can hope to compute numerically solutions of the problem : it has been done first in [@Benamou2000] and in many subsequent works detailed below. In all of these works, problem is discretized in space and time, yielding a finite dimensional convex optimization problem. Thanks to the general theory of convex optimization, there are strong guarantees that iterative schemes will indeed provide a very good approximation of the finite dimensional problems resulting from the discretization of . However, to the best of the author’s knowledge, there are few proofs that solutions of the finite dimensional problems will indeed converge to solutions of as the spatial and temporal grids are refined.

This problem is delicate and does not fit in classical theories for two (linked) reasons. First, the natural framework for the continuous problem is the set of positive measures: indeed, admits a solution if $\rho_0$ and $\rho_1$ the initial and final values of the density are fixed and equal to arbitrary positive measures sharing the same total mass. Moreover, there is no spatial regularizing effect: if $\rho_0$ is a Dirac mass in $x$ while $\rho_1$ is a Dirac mass is $y$, then one solution is $t \mapsto \rho_t = \delta_{\gamma(t)}$, where $\gamma$ is a constant-speed geodesic joining $x$ to $y$. Thus a satisfying convergence result should work even if the datum are spatially very singular which makes all the theory of finite elements phrased in Hilbert spaces is not available. In addition, the function $(\rho, {\mathbf{m}}) \mapsto |{\mathbf{m}}|^2/(2 \rho)$ is not continuous at the point $(0,0)$, meaning that one has to take care of what happens when there is void, i.e. when the density vanishes. As there is no spatial regularization effect, one really have to face this issue. In the infinite-dimensional theory this issue is now well-understood (see for instance [@AGS Chapter 8]), but it is not immediate to see how such theory will apply for discretizations of .

Related work
------------

(1,0) – (19,0) ; (10,-0.5) node\[below\][$\tau$ temporal step size]{} ;

(19,0) circle (0.2) ; (19,0) node\[below\][$\tau = 0$]{} ;

(0,1) – (0,11) ; (-0.5,6) node\[above, rotate = 90\][ $\sigma$ spatial step size]{} ;

(0,11) circle (0.2) ; (0,11) node\[above\][$\sigma = 0$]{} ;

(19,11) node\[draw\] ;

(6,3) node\[draw\] ;

(19,3) node\[draw\] ;

(10,3) – (15,3) ; (12.5,3) node\[below\][ ]{} ;

(19,5) – (19,9) ; (19,7) node\[right\][ ]{} ;

(6,5) – (15,10) ; (12,6.5) node\[right\][ ]{} ;

(6,5.5) – (15,10.5) ; (10,8) node\[left\][ ]{} ;

Since the seminal article [@Benamou2000], formulation and its variations have been used and solved numerically after discretization in a variety of context: congested dynamic [@Buttazzo2009], Wasserstein gradient flows [@Benamou2016gradient; @Carrillo2019], variational mean field games [@Benamou2015; @Benamou2016], unbalanced optimal transport [@Chizat2018], Wasserstein geodesics over surfaces [@Lavenant2018] to give a few examples.

Since the augmented Lagrangian method proposed in [@Benamou2000] to solve the (finite-dimensional) problem obtained after discretization, variants have been proposed in the more general framework of *proximal splitting* algorithms [@Papadakis2014; @Carrillo2019], and recently with the use of a clever Helmholtz-Hodge decomposition to handle the divergence constraint [@Henry2019]. Proofs of convergence of these convex optimization algorithms in the infinite dimensional case were provided by [@Guittet2003; @Hug2016] even though they are phrased in the framework of Hilbert spaces (that is the authors work with densities being $L^2$ w.r.t. space and time) and not in the one of positive measures.

Starting with considerations far from numerics, Maas [@Maas2011] defined a notion of optimal transport on graphs, which can be read as a semi-discretization of : the time variable is kept continuous, but the spatial variable is replaced by a discrete one, namely the objects are defined on either the vertices or edges of a graph. In this framework, it was later proved in [@Gigli2013] and [@Gladbach2018] that, if the discretization of space is refined (if one interprets the graphs as a finite volume discretization of the space), then on gets a finer and finer approximation of : actually it is phrased as a Gromov-Hausdorff convergence of metric spaces (namely, convergence of the spaces of probability measures endowed with the Wasserstein distance). More is said in Subsection \[subsection\_finite\_volumes\], as the discretization proposed by Maas and collaborators fits into our framework. We mention that [@Erbar2017] offers a temporal discretization of Problem if the space is already discretized as a graph, and proves convergence of solutions when one refines the temporal discretization, see Figure \[figure\_previous\_results\] for a schematic view of these different results. However, if one tries to combine convergence under refinement of temporal discretizations [@Erbar2017] and spatial discretization [@Gladbach2018], it leads to a restriction on the ratio between the temporal and spatial step sizes. On the contrary, we will show in the present work that convergence holds without any such restriction.

Eventually, in [@Carrillo2019], the authors prove convergence of a fully discretized version of to the original problem, though this is not the main point of their article. However, they work only with a finite difference discretization and make strong regularity assumptions on the solution of the problem (density bounded from below by a strictly positive constant, and smooth density and momentum) which are not satisfied if one works with arbitrary positive measures as temporal boundary conditions. As the authors point out [@Carrillo2019 Remark 8], they observe numerically that their method seems to work even if their strong regularity assumptions are not satisfied.

Contribution and organization of this article
---------------------------------------------

The goal of this article is to provide a framework which automatically guarantees converge of solutions of a fully discretized version of to the ones of . We will allow initial and final values of the density to be arbitrary positive measures (sharing the same total mass), hence do not require any regularity of the solution of . Moreover, we will give generic conditions on the spatial discretization for this convergence to hold: these conditions will be generic enough to encompass the finite volume discretization of [@Gladbach2018] as well as the discretization on triangulations of surfaces proposed by the present author in [@Lavenant2018]. The underlying spatial space will be allowed to be either a convex domain in a Euclidean space or a Riemannian manifold as in [@Lavenant2018].

However our result is not quantitative at all: a study of the speed of convergence would likely depend on the regularity of the solution of the infinite dimensional problem but we prefer to focus on the most generic case for the original problem . Moreover, the functionals that we are minimizing are $1$-homogeneous, hence not strictly convex: even a speed of convergence for the value of the problems would not automatically lead to a speed of convergence for the solutions. To the best of our knowledge, this question is completely open.

In the rest of this article, we first state precisely the problem we aim to discretize, and sufficient conditions on the spatial discretization to ensure convergence of solutions of fully space-time discretized problem: this is the object of Section \[section\_framework\]. The proof of our convergence result is provided in Section \[section\_proof\]. Then in Section \[section\_examples\] we show that our framework applies to spatial discretizations already present in the litterature, namely the one of Gladbach, Kopfer and Maas [@Gladbach2018], as well as the one we proposed with Claici, Chien and Solomon [@Lavenant2018]. Eventually, we show in Section \[section\_extension\] that we few additional work, one can add a penalization of the final density in the functional to be minimized (retrieving one step of the JKO scheme for Wasserstein gradient flows) while still having guarantees of convergence, but that the addition of a running cost depending on the density yields more involved issues that we do not cover.

We want to emphasize that the methods of proof and general ideas of this article are very much inspired from the aforementioned works [@Gladbach2018; @Erbar2017; @Carrillo2019]. Our techniques are closely related to theirs. However, we think that the framework that we propose shades a new light on their ideas and clarifies between what is really necessary to get this kind of convergence and what was contingent to the specific choice of discretization.

Framework and statement of the result
=====================================

If $A$ is a subset of $B$, the function ${\mathbbm{1}}_A$, defined on $B$, takes the value $1$ on $A$ and $0$ elsewhere. In the rest of this article, we will use $C$ to denote a constant independent on some parameters (specified in the context) whose valued may change from line to line. Similarly, we will use $(\varepsilon_\sigma)_\sigma$ to denote a generic function which tends to $0$ as $\sigma \to 0$, but which may change from line to line.

\[section\_framework\]

Infinite dimensional problem
----------------------------

As far as space is concerned, we work with $(X,g)$ a smooth compact Riemannian manifold possibly with a non-empty boundary. The boundary of $X$, denoted by ${\partial}X$ is also assumed to be infinitely smooth and, importantly, *convex*. For the definition of a convex boundary, we refer to [@Wang2009], see also [@Bartolo2002] for a geometric point of view. As an example, the reader can have in mind convex bounded domains in the Euclidean space with smooth boundary or smooth compact Riemannian manifold without boundary.

The smoothness and convexity assumptions on the boundary will be used only in Proposition \[prop\_regularization\] to regularize curves of measures. However, experiments [@Lavenant2018 Fig. 10] suggest that they probably are not optimal.

For each point $x \in X$, we denote the tangent space by $T_x X$, this space is equipped with a scalar product $g_x$ and a norm $| \ |_x$. The volume measure is denoted by ${\mathrm{d}}x$, if not specified integration is always performed w.r.t. this measure. The Riemannian distance is $d_g$. The tangent bundle is denoted by $T X$ and $\pi : TX \to X$ is the canonical projection. The space $C(X)$ is the one of continuous functions over $X$, and by a slight abuse of notations we denote by $C(X, TX)$ the space of continuous vector fields over $X$, i.e. of mappings ${\mathbf{v}}: X \to TX$ such that $\pi({\mathbf{v}}(x)) =x$ for all $x \in X$. We will use bold letters to denote vector fields. These spaces are endowed with the supremum norm, generating the topology of uniform convergence.

We denote by respectively ${\mathcal{M}}(X)$ and ${\mathcal{M}}(TX)$ the (topological) duals of $C(X)$ and $C(X, TX)$ and the duality products will be denoted $\langle \cdot, \cdot \rangle$. We use the same letters for duality products in $C(X)$ and $C(X,TX)$ as it will be clear from context which one is used. The space ${\mathcal{M}}(X)$, which is nothing else than the space of Borel measures on $X$, contains the convex subspace ${\mathcal{M}}_+(X)$ made of positive measures. The operator norm of an element $\rho \in {\mathcal{M}}(X)$ is nothing else than the total variation norm of the measure and is denoted by $\|\rho \|$, and analogously for ${\mathcal{M}}(TX)$. In particular, if $\rho \in {\mathcal{M}}_+(X)$, $\|\rho \| = \langle \rho, 1 \rangle$ and such a quantity is called the *mass* of $\rho$. We will often identify a measure with its density w.r.t. the volume measure.

An important particular case (where our results are already new and interesting) is the flat one, namely when $X = \Omega \subset {\mathbb{R}}^d$ is a closed convex subset of ${\mathbb{R}}^d$ with smooth boundary. In this case, the volume measure is the Lebesgue measure, $C(X) = C(\Omega)$ and $C(X, TX) = C(\Omega, {\mathbb{R}}^d)$ is the set of continuous functions valued into ${\mathbb{R}}^d$. Moreover, ${\mathcal{M}}(TX) = {\mathcal{M}}(\Omega)^d$ is nothing else than the set of measures valued in ${\mathbb{R}}^d$ and defined over $\Omega$.

If $\rho \in {\mathcal{M}}(X)$ and ${\mathbf{m}}\in {\mathcal{M}}(TX)$, we define $$\label{equation_Benamou_Brenier_x}
A(\rho, {\mathbf{m}}) := \sup_{a,{\mathbf{b}}} \left\{ \langle  \rho, a \rangle + \langle  {\mathbf{m}}, {\mathbf{b}}\rangle \ : \ a \in C(X), \ {\mathbf{b}}\in C(X, TX) \text{ and } \forall x \in X, \  a(x) + \frac{|{\mathbf{b}}(x)|^2_x}{2} \leqslant 0 \right\}.$$ This is the so-called Benamou-Brenier formula: it is clearly a convex lower semi-continuous functional, and (see [@SantambrogioOTAM Proposition 5.18]) it is finite if and only if $\rho \in {\mathcal{M}}_+(X)$ and ${\mathbf{m}}$ has a density ${\mathbf{v}}: X \to TX$ w.r.t. $\rho$, and in such a case $$A(\rho, {\mathbf{m}}) = A(\rho, \rho {\mathbf{v}}) = \int_X \frac{1}{2} |{\mathbf{v}}(x)|^2_x ~ \rho({\mathrm{d}}x).$$ In other words one can write $A(\rho, {\mathbf{m}}) = \int_X |{\mathbf{m}}|^2/(2 \rho)$ with the convention that $$\label{equation_convention_b22a}
\frac{|b|^2}{2a} := \begin{cases}
+ \infty & \text{if } a < 0, \\
+ \infty & \text{if } a = 0 \text{ and } |b| > 0, \\
0 & \text{if } a=|b| = 0, \\
{\displaystyle{\frac{|b|^2}{2a}}} & \text{otherwise}.
\end{cases}$$ We will keep such a convention in the rest of the article. Moreover, let us remark that if $\rho \in {\mathcal{M}}_+(X)$, then clearly one can take $a = - |{\mathbf{b}}|^2/2$ in formulation . Hence, if one defines $A^\star$ on ${\mathcal{M}}(X) \times C(X,TX)$ by $$A^\star(\rho, {\mathbf{b}}) := \frac{1}{2} \int_X |{\mathbf{b}}(x)|_x^2 ~ \rho({\mathrm{d}}x) = \frac{1}{2} \langle \rho, |{\mathbf{b}}|^2 \rangle,$$ then one can simply write $$A(\rho, {\mathbf{m}}) = \sup_{{\mathbf{b}}\in C(X,TX)} \left( \langle {\mathbf{m}}, {\mathbf{b}}\rangle - A^\star(\rho, {\mathbf{b}}) \right)$$ at least if $\rho \in {\mathcal{M}}_+(X)$, that is if we assume that $\rho$ is positive. The functional $A$ is nothing else than the Legendre transform of $A^\star$ w.r.t. its second variable.

In the article [@Gladbach2018], the notations for $A$ and $A^\star$ are swapped compared to here. This is because we try to work as much a possible with $(\rho, {\mathbf{m}})$ as unknowns which are thought as “primal” variables, even though rigorously they should be considered as “dual” variables.

As far as time is concerned, we assume that the initial time is $t=0$ and the final time $t=1$. Hence the temporal domain is $[0,1]$. This space is endowed with the Lebesgue measure ${\mathrm{d}}t$, and we denote by ${\mathrm{d}}t \otimes {\mathrm{d}}x$ the measure on $[0,1] \times X$ which is the tensorial product of the Lebesgue measure on $[0,1]$ and the volume measure on $X$. We naturally define the spaces $C([0,1] \times X)$ and $C([0,1] \times X, TX)$ of respectively space-time dependent scalar functions and space-time dependent velocity fields. They are endowed with the supremum norm and their topological duals are denoted by ${\mathcal{M}}([0,1] \times X)$ and ${\mathcal{M}}([0,1] \times TX)$. Similarly to the previous case, ${\mathcal{M}}_+([0,1] \times X)$ denotes the subset of ${\mathcal{M}}([0,1] \times X)$ made of positive measures. The norm (operator norm which is the total variation norm) of $\rho \in {\mathcal{M}}([0,1] \times X)$ is denoted by $\|\rho \|$, analogously for ${\mathcal{M}}([0,1] \times TX)$. To make the distinction apparent, we will use double brackets $\llangle \cdot, \cdot \rrangle$ for the duality products of space-time dependent objects. The Benamou-Brenier formula can be extended to space-time dependent objects: if $(\rho, {\mathbf{m}}) \in {\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$, $$\begin{gathered}
{\mathcal{A}}(\rho, {\mathbf{m}}) := \sup_{a,{\mathbf{b}}} \Bigg\{  \llangle \rho, a \rrangle + \llangle {\mathbf{m}}, {\mathbf{b}}\rrangle \ : \ a \in C([0,1] \times X), \ {\mathbf{b}}\in C([0,1] \times X, TX) \\
\text{ and } \forall (t,x) \in [0,1] \times X, \  a(t,x) + \frac{|{\mathbf{b}}(t,x)|^2_x}{2} \leqslant 0 \Bigg\}.\end{gathered}$$ Similarly, we define ${\mathcal{A}}^\star$ on ${\mathcal{M}}([0,1] \times X) \times C([0,1] \times X, X)$ by $${\mathcal{A}}^\star(\rho, {\mathbf{b}}) := \frac{1}{2} \iint_{[0,1] \times X} |{\mathbf{b}}(t,x)|_x^2 ~ \rho({\mathrm{d}}t, {\mathrm{d}}x) = \frac{1}{2} \llangle \rho, |{\mathbf{b}}|^2 \rrangle,$$ and one can simply write, if $\rho \in {\mathcal{M}}_+([0,1] \times X)$ $${\mathcal{A}}(\rho, {\mathbf{m}}) = \sup_{{\mathbf{b}}\in C([0,1] \times X,TX)} \left( \llangle {\mathbf{m}}, {\mathbf{b}}\rrangle - {\mathcal{A}}^\star(\rho, {\mathbf{b}}) \right).$$ Actually, thanks to a smoothing argument we can even take the supremum over vector fields ${\mathbf{b}}$ in $C^1([0,1] \times X, TX)$ and the identity stays valid.

The gradient $\nabla$ maps $C^1(X)$ into $C(X,TX)$ while the divergence $\nabla \cdot$ is the adjoint of the gradient w.r.t. the $L^2$ scalar product (weighted by the volume measure). Namely, if $f \in C^1(X)$ and $\mathbf{g} \in C^1(X, TX)$ vanish on a neighborhood of ${\partial}X$, $$\int_X g_x(\nabla f(x), \mathbf{g}(x)) ~ {\mathrm{d}}x = - \int_X f(x) ( \nabla \cdot \mathbf{g} )(x) ~ {\mathrm{d}}x.$$

If $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$, we define ${\mathcal{CE}}(\rho_0, \rho_1)$ the set of pairs $(\rho, {\mathbf{m}})$ which satisfies the continuity equation with initial and final values given by respectively $\rho_0$ and $\rho_1$ as follows: $$\begin{aligned}
{\mathcal{CE}}(\rho_0, \rho_1) := \Bigg\{ (\rho, {\mathbf{m}}) \ : \ \forall \phi \in C^1([0,1] \times X), \ 
\llangle \rho, {\partial}_t \phi \rrangle + \llangle {\mathbf{m}}, \nabla \phi \rrangle = \langle \phi(1, \cdot) , \rho_1 \rangle - \langle \phi(0, \cdot),  \rho_0  \rangle \Bigg\} \\ \subset{\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX). \end{aligned}$$ Indeed, it is nothing else than a weak formulation of the continuity equation ${\partial}_t \rho + \nabla \cdot {\mathbf{m}}= 0$ with no-flux boundary conditions on ${\partial}X$. In particular, if $\rho_0$ and $\rho_1$ do not have the same total mass, it is easy to see that this set is empty by testing with a constant function.

\[rmk\_rho\_curve\] Let $(\rho, {\mathbf{m}})$ such that $(\rho, {\mathbf{m}}) \in {\mathcal{CE}}(\bar{\rho}_0, \bar{\rho}_1)$ for some $\bar{\rho}_0, \bar{\rho}_1 \in {\mathcal{M}}_+(X)$ and ${\mathcal{A}}(\rho, {\mathbf{m}}) < + \infty$. Thanks to the continuity equation, it is not difficult to see that the temporal marginal of $\rho$ is proportional to the Lebesgue measure on $[0,1]$. In particular, we can define $(\rho_t)_{t \in [0,1]}$ a curve valued in ${\mathcal{M}}_+(X)$ as the disintegration of $\rho$ w.r.t. its temporal marginal [@AGS Theorem 5.3.1]. Thanks to the continuity equation and the estimate ${\mathcal{A}}(\rho, {\mathbf{m}}) < + \infty$, in fact the map $t \mapsto \rho_t \in {\mathcal{M}}_+(X)$ is continuous for the topology of weak convergence (actually even $1/2$-Hölder w.r.t. the Wasserstein distance) and $\bar{\rho}_0$ (resp. $\bar{\rho}_1$) is the limit as $t \to 0$ (resp. $t \to 1$) of $\rho_t$. In short, $\rho \in {\mathcal{M}}_+([0,1 \times X])$ can also be seen as a continuous curve valued in ${\mathcal{M}}_+(X)$, and $\bar{\rho}_0, \bar{\rho}_1$ are the values of the curve at the temporal boundaries.

Let $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ sharing the same total mass. We define the functional ${\mathcal{J}}_{\rho_0, \rho_1}$ on ${\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ by $${\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) := \begin{cases}
{\mathcal{A}}(\rho, {\mathbf{m}}) & \text{if } (\rho, {\mathbf{m}}) \in {\mathcal{CE}}(\rho_0, \rho_1), \\
+ \infty & \text{otherwise}.
\end{cases}$$ Given what is said above, minimizing ${\mathcal{J}}_{\rho_0, \rho_1}$ amounts to solve problem described in the introduction. This is a convex optimization problem defined on the space ${\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ endowed with the topology of weak convergence.

We will denote by $W_2$ the quadratic Wasserstein distance, see Appendix \[section\_wasserstein\]. The functional ${\mathcal{J}}_{\rho_0, \rho_1}$ is closely related to this distance as mentioned in the introduction and the next theorem summarizes such a link. However, we will not rely on it in the sequel: the reader unfamiliar with the theory of optimal transport can skip it for a first reading.

\[theorem\_J\_W2\] Under the assumption that $\rho_0$ and $\rho_1$ share the same total mass, there exists a minimizer to the functional ${\mathcal{J}}_{\rho_0, \rho_1}$ and $$\label{equation_J_Wasserstein}
\min_{{\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)} \  {\mathcal{J}}_{\rho_0, \rho_1} = \frac{1}{2} W_2^2(\rho_0, \rho_1).$$ Moreover, at least in the case where $X$ is a flat convex domain or a Riemannian manifold without boundary and if $\rho_0$ or $\rho_1$ is absolutely continuous w.r.t. the volume measure then the minimizer is unique.

Existence is really easy from the direct method of calculus if variations. The only non trivial thing to check is the existence of at least one competitor, but this it is a byproduct of the proof of .

For the proof of , one can refer to [@SantambrogioOTAM Theorem 5.28] for $X$ being a flat convex domain and [@Erbar2010 Proposition 2.5] when $X$ is a Riemannian manifold without boundary. We have not found the case of a Riemannian manifold *with* boundary written explicitly in the literature, but following [@Villani2003 Remark 8.3] the only requirement is to prove that a manifold with boundary can be isometrically embedded in a Euclidean space. Such a property is true, as a Riemannian manifold with boundary can always be seen as a subset of a Riemannian manifold without boundary [@Pigola2016 Corollary B], and the latter can be embedded isometrically in a Euclidean space thanks to Nash’s theorem.

Eventually, uniqueness of the minimizer comes from the proof of uniqueness for the optimal transport problem: see [@SantambrogioOTAM Proposition 5.32 and Theorem 1.22] for the flat case and [@Villani2008 Corollary 7.23] combined with [@Mccann2001 Theorems 8 and 9] on a manifold without boundary.

The goal of this article is to explain how one can approximate ${\mathcal{J}}_{\rho_0, \rho_1}$ with functionals defined on finite dimensional spaces.

Discretization
--------------

The discretizations already proposed in the literature share some similar structure as they try to mimick the one of the infinite dimensional problem. They split the temporal and spatial variables, and although the temporal discretization is rather straightforward, more diverse propositions have been made for the spatial one. As far as the latter is concerned, what is needed is at least an equivalent of the action, and a divergence operator. We assume that we have some finite dimensional approximations of such objects. Namely, we will have a family of models indexed by a parameter $\sigma$ thought as a spatial step size. The parameter $\sigma$ belongs to a subset $\Sigma$ of $(0, + \infty)$ which contains $0$ as an accumulation point and in the sequel “for all $\sigma > 0$” means “for all $\sigma$ in the the set $\Sigma$”.

\[definition\_approximation\] A family of finite dimensional models of dynamical optimal transport is the datum, for any $\sigma > 0$ of $$( \underbrace{{\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma}_{\text{``spaces''}}, 
\underbrace{A_\sigma}_{\text{``action''}},
\underbrace{{\mathrm{Div}}_\sigma}_{\text{``derivation''}} )$$ where:

1.  ${\mathcal{X}}_\sigma$ and ${\mathcal{Y}}_\sigma$ are finite dimensional vector spaces.

2.  $A_\sigma : {\mathcal{X}}_\sigma \times {\mathcal{Y}}_\sigma \to [0, + \infty]$ is a proper convex functional $(-1)$-positively homogeneous in its first variable and $2$-homogeneous in its second variable. We require that $A_\sigma$ is non increasing on a cone ${\mathcal{X}}_{\sigma,+} \subset {\mathcal{X}}_\sigma$ in the sense that for every $P_1, P_2 \in {\mathcal{X}}_{\sigma,+}$ and any ${\mathbf{M}}\in {\mathcal{Y}}_\sigma$, $$A_\sigma(P_1+ P_2, {\mathbf{M}}) \leqslant A_\sigma(P_1, {\mathbf{M}}).$$

3.  ${\mathrm{Div}}_\sigma : {\mathcal{Y}}_\sigma \to {\mathcal{X}}_\sigma$ is a linear operator.

We will think of ${\mathcal{X}}_\sigma$ and ${\mathcal{Y}}_\sigma$ as approximations of the space ${\mathcal{M}}(X)$ and ${\mathcal{M}}(TX)$ respectively, while $A_\sigma$ mimics the action $A$ and ${\mathrm{Div}}_\sigma$ corresponds the divergence operator. The space ${\mathcal{X}}_{\sigma,+}$ is thought as ${\mathcal{M}}_+(X)$. We emphasize that we do not require any duality structure on ${\mathcal{X}}_\sigma$ and ${\mathcal{Y}}_\sigma$, nor any norm on them. For the moment, we have made no explicit link between these objects and the original ones living on the manifold $X$. As we expect the approximation to be finer and finer as $\sigma \to 0$, the dimension of ${\mathcal{X}}_\sigma$ and ${\mathcal{Y}}_\sigma$ are thought as increasing as $\sigma$ decreases.

By convention, we will use capital letters to denote the counterpart in the spaces ${\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma$ of elements in ${\mathcal{M}}(X)$ and ${\mathcal{M}}(TX)$. Hence a generic element of ${\mathcal{X}}_\sigma$ is $P$ (a capitalized $\rho$), and a generic element in ${\mathcal{Y}}_\sigma$ is ${\mathbf{M}}$.

Next, we take $N+1 \geqslant 2$ the number of spatial time steps, and we use $\tau = 1/N$ to denote the temporal step size.

Let $N+1 \geqslant 2$ be given and let us define $\tau=1/N$. Moreover, let $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_{\sigma}$ a family of finite dimensional models of dynamical optimal transport.

If $\sigma > 0$ is given, and $\bar{P}_0$ and $\bar{P}_1$ are elements in ${\mathcal{X}}_{\sigma,+}$, we define the functional ${\mathcal{J}}^{N, \sigma}_{\bar{P}_0, \bar{P}_1}$ on $({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^{N}$ by $$\label{equation_discrete_cost}
{\mathcal{J}}^{N, \sigma}_{\bar{P}_0, \bar{P}_1}( (P_k)_{0 \leqslant k \leqslant N}, ({\mathbf{M}}_k)_{1 \leqslant k \leqslant N}  ) = \tau \sum_{k=1}^{N} A_\sigma \left( \frac{P_{k-1} + P_{k}}{2}, {\mathbf{M}}_k \right)$$ if $(P_k)_{0 \leqslant k \leqslant N} \in ({\mathcal{X}}_{\sigma,+})^{N+1}$ and the discrete continuity equation is satisfied, that is $$\begin{cases}
\tau^{-1} (P_{k} - P_{k-1})  + {\mathrm{Div}}_\sigma ({\mathbf{M}}_k) = 0, & \forall k \in \{1,2, \ldots, N \} \\
P_{0} = \bar{P}_0, & \\
P_{N} = \bar{P}_1, & \\
\end{cases}$$ and $+ \infty$ otherwise.

The function ${\mathcal{J}}^{N, \sigma}_{\bar{P}_0, \bar{P}_1}$ is mimicking the continuous one ${\mathcal{J}}_{\rho_0, \rho_1}$. Actually what we are doing, given the continuity equation, is piecewise affine interpolation in time for the densities, and piecewise constant interpolation in time for the momentum. On the other hand, for the cost functional , we transform the piecewise affine interpolation of the density into a piecewise constant one by averaging in time before computing the action. This averaging is necessary because we want to avoid at any cost putting momentum where there is no mass. To understand the issue, let us go back for a moment to the continuous case and let us focus on a single time step: assume that $(\rho, {\mathbf{m}}) \in {\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ are respectively affine and constant w.r.t. the time variable. That is, $\rho(t,\cdot) = (1-t) \rho(0,\cdot) + t \rho(1, \cdot)$ and ${\mathbf{m}}(t, \cdot) = {\mathbf{m}}(\cdot)$. To simplify the exposition, assume that they have densities w.r.t. the volume measure. Then $${\mathcal{A}}(\rho, {\mathbf{m}}) = \iint_{[0,1] \times X} \frac{|{\mathbf{m}}(x)|^2}{2 \rho(t,x)} \, {\mathrm{d}}t \, {\mathrm{d}}x =  \int_X \frac{|{\mathbf{m}}(x)|^2}{2} \underbrace{\left(\int_0^1 \frac{{\mathrm{d}}t}{ (1-t) \rho(0,x) + t \rho(1,x) } \right)}_{= + \infty \text{ if  } \rho(0,x) \text{ or }\rho(1,x) = 0} \, {\mathrm{d}}x.$$ Hence ${\mathcal{A}}(\rho, {\mathbf{m}}) = + \infty$ if $|{\mathbf{m}}(x)|^2 > 0$ and *either* $\rho(0,x)$ or $\rho(1,x)$ vanishes for $x$ belonging to a set of positive measure. That is, ${\mathbf{m}}$ must have a density w.r.t. both $\rho(0, \cdot)$ and $\rho(1, \cdot)$ for the action to be finite. Such a condition is very restrictive (think that $\rho(0, \cdot)$ could be a Dirac mass), and things are much better with our choice of discretization where ${\mathbf{m}}$ must have a density w.r.t. $\rho(0, \cdot) + \rho(1, \cdot)$. We refer the reader to the proof of Proposition \[prop\_regularization\] where the particular choice of averaging in time for the densities is crucial.

Minimizing the function ${\mathcal{J}}^{N, \sigma}_{\bar{P}_0, \bar{P}_1}$ is a finite dimensional convex optimization problem. As ${\mathcal{J}}^{N, \sigma}_{\bar{P}_0, \bar{P}_1} \geqslant 0$, it always has a solution if admissible competitors exist. A consequent effort has been devoted to solve it efficiently as recalled earlier. We emphasize that for the actual computation of a minimizer of ${\mathcal{J}}^{N, \sigma}_{\bar{P}_0, \bar{P}_1}$, one can (and very often does) introduce an additional duality structure on ${\mathcal{X}}_\sigma$ and ${\mathcal{Y}}_\sigma$.

Before stating the assumptions sufficient to make the link between these finite dimensional problems and the infinite dimensional one, we need to introduce some additional vocabulary.

First, in the sequel we will require some estimates to hold uniformly if some test functions are uniformly regular. Hence we need to define some way to measure regularity of test functions on the manifold $X$. To that extent, we use the following definition.

Let $(U^\alpha, \varphi^\alpha)_{\alpha}$ be a *finite* atlas of $X$, where each open set $U^\alpha \subset X$ is mapped by $\varphi_\alpha$ on a open set of either a Euclidean space or a half Euclidean space.

For $q \geqslant 1$ integer, we say that ${\mathcal{B}}\subset C(X)$ is a *bounded* set of $C^q(X)$ if $$\sup_\alpha \sup_{f \in {\mathcal{B}}} \| (f \circ (\varphi^\alpha)^{-1}) \|_{C^q( \varphi^\alpha(U^\alpha)  )} < + \infty$$ Similarly, we say that ${\mathcal{B}}\subset C(X, TX)$ is a *bounded* set of $C^q(X, TX)$ if $$\sup_\alpha \sup_{f \in {\mathcal{B}}} \| D \varphi^\alpha (f \circ (\varphi^\alpha)^{-1}) \|_{C^q( \varphi^\alpha (U^\alpha)  )} < + \infty$$

In the definition above, the $C^q$ norm of a function defined on a Euclidean space is just the supremum norm of the function and all its partial derivatives up to order $q$. As $X$ is compact and smooth, these definitions do not depend on the choice of the atlas.

W will also need to introduce a duality structure on ${\mathcal{Y}}_\sigma$. Indeed, at least when studying the discretization of [@Gladbach2018], we found this duality structure necessary. We emphasize, however, that to *define* the finite dimensional functional ${\mathcal{J}}^{N, \sigma}$ this duality structure was not needed, nor in the *consequences* of Theorem \[theo\_main\] below.

Specifically, let ${\mathcal{Y}}'_\sigma$ be the dual of ${\mathcal{Y}}_\sigma$, that is the set of linear form on ${\mathcal{Y}}_\sigma$. If $L : {\mathcal{Y}}_\sigma \to {\mathcal{M}}(TX)$, we can define its *adjoint* operator $L^\top : C(X,TX) \to {\mathcal{Y}}'_\sigma$ by the following identity: for any ${\mathbf{b}}\in C(X,TX)$ and ${\mathbf{M}}\in {\mathcal{Y}}_\sigma$, $$\langle {\mathbf{M}}, L^\top({\mathbf{b}})  \rangle_{{\mathcal{Y}}_\sigma \times {\mathcal{Y}}'_\sigma} = \langle L({\mathbf{M}}), {\mathbf{b}}\rangle_{ {\mathcal{M}}(TX) \times C(X,TX) }.$$ As ${\mathcal{Y}}_\sigma$ is finite dimensional, $L^\top({\mathbf{b}})$ is well-defined for any ${\mathbf{b}}\in C(X, TX)$. Moreover, for a fixed $\sigma > 0$ we define $A^\star_\sigma : {\mathcal{X}}_\sigma \times {\mathcal{Y}}'_\sigma \to [0, + \infty)$ to be the Legendre transform of $A_\sigma$ w.r.t. its second variable: if ${\mathbf{B}}\in {\mathcal{Y}}'_\sigma$ and $P \in {\mathcal{X}}_\sigma$ then $$A^\star_\sigma(P, {\mathbf{B}}) = \sup_{{\mathbf{M}}\in {\mathcal{Y}}_\sigma} \left( \langle {\mathbf{M}}, {\mathbf{B}}\rangle_{{\mathcal{Y}}_\sigma \times {\mathcal{Y}}'_\sigma} - A_\sigma(P, {\mathbf{M}}) \right).$$ As ${\mathcal{Y}}_\sigma$ is a finite-dimensional and $A_\sigma$ is proper, the following identity holds [@Rockafellar1970 Theorem 12.2]: for any $P \in {\mathcal{X}}_{\sigma,+}$ and ${\mathbf{M}}\in {\mathcal{Y}}_\sigma$, $$A_\sigma(P, {\mathbf{M}}) = \sup_{{\mathbf{B}}\in {\mathcal{Y}}'_\sigma} \left( \langle {\mathbf{M}}, {\mathbf{B}}\rangle_{{\mathcal{Y}}_\sigma \times {\mathcal{Y}}'_\sigma} - A^\star_\sigma(P, {\mathbf{B}}) \right).$$ Moreover, one can easily check that $A^\star_\sigma$ is also $2$-homogeneous in its second variable. In particular, a scaling argument leads to, for any $P \in {\mathcal{X}}_\sigma$, ${\mathbf{M}}\in {\mathcal{Y}}_\sigma$ and ${\mathbf{B}}\in {\mathcal{Y}}'_\sigma$, $$\label{equation_control_A_Astar}
\langle {\mathbf{M}}, {\mathbf{B}}\rangle_{{\mathcal{Y}}_\sigma \times {\mathcal{Y}}'_\sigma} \leqslant \sqrt{A_\sigma(P, {\mathbf{M}})} \sqrt{A^\star_\sigma(P, {\mathbf{B}})}.$$

Now we can state the assumption necessary to make the link between this finite dimensional problems and the infinite dimensional one. It amounts to give ourselves reconstruction and sampling operators sending ${\mathcal{X}}_\sigma$ onto ${\mathcal{M}}(X)$, ${\mathcal{Y}}_\sigma$ onto ${\mathcal{M}}(TX)$ and *vice versa*, and assume some sort of *commutation* relations between them.

\[definition\_adapated\] Let $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_{\sigma}$ be a family of finite dimensional models of dynamical optimal transport. Such a family is said *adapted* to the manifold $X$ if there exists, for any $\sigma > 0$, “reconstruction operators” $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}, R_{{\mathcal{Y}}_\sigma}$ and “sampling operators” $S_{{\mathcal{X}}_\sigma}, S_{{\mathcal{Y}}_\sigma}$ where:

1.  $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}  : {\mathcal{X}}_\sigma \to {\mathcal{M}}(X)$ and $R_{{\mathcal{Y}}_\sigma} : {\mathcal{Y}}_\sigma \to {\mathcal{M}}(TX)$ are linear operators. They are defined everywhere on ${\mathcal{X}}_\sigma$ and ${\mathcal{Y}}_\sigma$. Moreover, both $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ and $R^{A}_{{\mathcal{X}}_\sigma}$ map ${\mathcal{X}}_{\sigma,+}$ into ${\mathcal{M}}_+(X)$.

2.  $S_{{\mathcal{X}}_\sigma} : {\mathcal{M}}(X) \to {\mathcal{X}}_\sigma$ and $S_{{\mathcal{Y}}_\sigma} : {\mathcal{M}}(TX) \to {\mathcal{Y}}_\sigma$ are linear operators. We require $S_{{\mathcal{X}}_\sigma}$ to be defined everywhere on ${\mathcal{M}}(X)$ and to map ${\mathcal{M}}_+(X)$ into ${\mathcal{X}}_{\sigma,+}$, while $S_{{\mathcal{Y}}_\sigma}$ can have a domain ${\mathcal{D}}(S_{{\mathcal{Y}}_\sigma})$ different from ${\mathcal{M}}(TX)$ which must at least contains measures which have a continuous densities (hence is dense in ${\mathcal{M}}(TX)$).

Moreover, these operators must satisfy the following properties.

1.  \[asmp\_interp\_sampling\] *(Reconstruction is asymptotically the inverse of sampling)* For any $\rho \in {\mathcal{M}}_+(X)$, there holds $$\lim_{\sigma \to 0} (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ S_{{\mathcal{X}}_\sigma})(\rho) = \rho$$ weakly in ${\mathcal{M}}(X)$.

2.  \[asmp\_interps\] *(Reconstructions of the density are asymptotically equivalent)* Let ${\mathcal{B}}$ be a bounded set of $C^1(X)$. There exists $(\varepsilon_\sigma)_{\sigma > 0}$ tending to $0$ as $\sigma \to 0$ such that, for any $\sigma > 0$, any $\phi \in {\mathcal{B}}$ and any $P \in {\mathcal{X}}_\sigma$, $$\langle ( R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} - R^A_{{\mathcal{X}}_\sigma} )(P), \phi  \rangle \leqslant  \varepsilon_\sigma  \|R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} (P) \|.$$

3.  \[asmp\_interp\_der\] *(Asymptotic commutation of reconstruction and derivation)* Let ${\mathcal{B}}$ be a bounded set of $C^2(X)$. There exists $(\varepsilon_\sigma)_{\sigma > 0}$ tending to $0$ as $\sigma \to 0$ such that, for any ${\mathbf{M}}\in {\mathcal{Y}}_\sigma$ and any $\phi \in {\mathcal{B}}$, $$\left| \langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ {\mathrm{Div}}_\sigma) ({\mathbf{M}}_\sigma), \phi \rangle + \langle R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}}_\sigma), \nabla \phi \rangle \right| \leqslant \varepsilon_\sigma   \| R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}})\|.$$

4.  \[asmp\_samp\_der\] *(Exact commutation of sampling and derivation)* For any ${\mathbf{m}}\in {\mathcal{M}}(TX)$ which has a density (still denoted by ${\mathbf{m}}$) w.r.t. volume measure in $C^1(X,TX)$ satisfying ${\mathbf{m}}(x) \cdot {\mathbf{n}}_X(x) = 0$ for all $x \in {\partial}X$, then for any $\sigma > 0$ $$S_{{\mathcal{X}}_\sigma} (\nabla \cdot {\mathbf{m}}) =  ({\mathrm{Div}}_\sigma \circ S_{{\mathcal{Y}}_\sigma})({\mathbf{m}}).$$

5.  \[asmp\_interp\_action\] *(One-sided estimates for the Legendre transform of the action)* Let ${\mathcal{B}}$ be a bounded set of $C^1(X,TX)$. There exists $(\varepsilon_\sigma)_{\sigma > 0}$ tending to $0$ as $\sigma \to 0$ such that the following holds. If ${\mathbf{b}}\in {\mathcal{B}}$ then for any $\sigma > 0$ and any $P \in {\mathcal{X}}_\sigma$, $$A_\sigma^\star( P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}) ) \leqslant  A^\star( R^A_{{\mathcal{X}}_\sigma} (P), {\mathbf{b}}) + \varepsilon_\sigma \| R^A_{{\mathcal{X}}_\sigma}(P) \|.$$ Moreover, there exists $C \geqslant 1$ such that for any ${\mathbf{b}}\in C(X,TX)$, there holds $$A_\sigma^\star( P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}) ) \leqslant \frac{C}{2} \| R^A_{{\mathcal{X}}_\sigma} (P)   \| \left( \sup_{x \in X} |{\mathbf{b}}(x)|_x \right)^2.$$

6.  \[asmp\_samp\_actio\] *(Asymptotic one-sided estimate for the action on smooth densities)* Let ${\mathcal{B}}$ a bounded set of $C^1(X)$ such that all functions in ${\mathcal{B}}$ are uniformly bounded from below by a positive constant. Let ${\mathcal{B}}'$ a bounded subset of $C^1(X,TX)$. Then there exists $(\varepsilon_\sigma)_{\sigma > 0}$ tending to $0$ as $\sigma \to 0$ such that, if $\rho \in {\mathcal{M}}(X)$ has a density in ${\mathcal{B}}$ and ${\mathbf{m}}\in {\mathcal{M}}(TX)$ has a density in ${\mathcal{B}}'$ then for any $\sigma > 0$ $$A_\sigma( S_{{\mathcal{X}}_\sigma} (\rho), S_{{\mathcal{Y}}_\sigma}({\mathbf{m}}) )  \leqslant A(\rho, {\mathbf{m}}) + \varepsilon_\sigma.$$

7.  \[asmp\_controllability\] *(Controllability)* There exists $(\varepsilon_\sigma)_{\sigma > 0}$ which tends to $0$ as $\sigma \to 0$ and $\omega$ is a continuous functions satisfying $\omega(0) = 0$, such that the following holds. If $x, y$ are points in $X$, then for any $\sigma > 0$ there exists $\hat{P} \in {\mathcal{X}}^+_\sigma$ and two elements $\hat{{\mathbf{M}}}_1, \hat{{\mathbf{M}}}_2 \in {\mathcal{Y}}_\sigma$ such that $$\label{equation_controllability}
    \begin{cases}
    {\mathrm{Div}}_\sigma (\hat{{\mathbf{M}}}_1) = \hat{P} - S_{{\mathcal{X}}_\sigma} (\delta_x) \\
    {\mathrm{Div}}_\sigma ({\mathbf{M}}_2) = \hat{P} - S_{{\mathcal{X}}_\sigma} (\delta_y)
    \end{cases} \ \ \text{ and } \ \forall i \in \{ 1,2 \}, \ \ A_\sigma(\hat{P}, \hat{{\mathbf{M}}}_i) \leqslant \omega(d_g(x,y)) + \varepsilon_\sigma.$$

Notice that we give ourselves $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ and $R^A_{{\mathcal{X}}_\sigma}$ two different reconstructions operators for the density. Indeed, even if they are assumed to be asymptotically equivalent, each of these operators will be used to pass to the limit different terms in the objective functional (namely the continuity equation with $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ and the action with $R^A_{{\mathcal{X}}_\sigma}$).

This list of conditions can seem quite long, the reader can take a look at the examples in Section \[section\_examples\] to see how one can check them in practice. Definition \[definition\_adapated\] is only concerned with the spatial discretization, but with its help we will be able to say something about the limit of the functional ${\mathcal{J}}^{N, \sigma}$ defined for space-time dependent objects.

Except for \[asmp\_controllability\] which is discussed in Remark \[rmk\_controllability\] below, all requirements are some sort of *commutation* relations. Except for \[asmp\_samp\_der\], all allow for some leeway, that is the presence of an error term. Moreover, although some differential operators are involved in \[asmp\_interp\_der\] and \[asmp\_samp\_der\], the remaining properties only deal with $0$-th order quantities: importantly the assumptions \[asmp\_interp\_action\] and \[asmp\_samp\_actio\] involving the action are only one-sided estimates. Eventually, notice that the condition of no-flux boundary condition appears only in \[asmp\_samp\_der\], and not in the way the action is discretized: in fact with these assumptions ${\mathcal{Y}}_\sigma$ is supposed to be a discretization of the space of momenta with no-flux boundary conditions.

In practice, one chooses the sampling operators in order for \[asmp\_samp\_der\] to be satisfied, namely exact commutation between sampling and derivation. For the reconstruction operators, $R^{A}_{{\mathcal{X}}_\sigma}$ and $R_{{\mathcal{Y}}_\sigma}$ are chosen such that the one-sided estimation on the action \[asmp\_interp\_action\] holds. Once this is done, one can choose $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ to get \[asmp\_interp\_der\] asymptotic commutation between reconstruction and derivation. Then, properties \[asmp\_interp\_sampling\], \[asmp\_interps\] and \[asmp\_samp\_actio\] are simple *consistency* properties to check, involving no derivatives.

(-4,5) node\[left\][Density]{} ; (0,-7) node\[left\][Momentum]{} ;

(0,0) rectangle (10, 10) ; (12,0) rectangle (22, 10) ; (24,0) rectangle (34, 10) ;

(-5,-1) – (39,-1) ;

(6,-2) rectangle (17, -12) ; (18,-2) rectangle (28, -12) ;

(2,2) circle (0.3) ; (5,10.5) node\[above\][$\delta_x$]{} ;

(32,8) circle (0.3) ; (29,10.5) node\[above\][$\delta_y$]{} ;

(14,2) – (20,8) ; (17,10.5) node\[above\][$\hat{\rho}$]{} ;

(7.85,-9.9) – (14.85,-3.9) – (15,-4.1) – (8,-10.1) – cycle ; (11,-12.5) node\[below\][$|\hat{{\mathbf{m}}}_1|$]{} ;

(19.85,-9.9) – (25.85,-3.9) – (26,-4.1) – (20,-10.1) – cycle ; (23,-12.5) node\[below\][$|\hat{{\mathbf{m}}}_2|$]{} ;

\[rmk\_controllability\] The less standard assumption seems to be \[asmp\_controllability\]. It is used to show how one can control the minimal value of ${\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}$ in terms of the Wasserstein distance between $\rho_0$ and $\rho_1$ even if $N$ is very small, see Proposition \[prop\_controllability\]. Checking such a condition relies on the particular discretization that is chosen, usually it consists in mimicking the proof of \[asmp\_controllability\] if there is no discretization, that is if ${\mathcal{X}}_\sigma = {\mathcal{M}}(X)$ and ${\mathcal{Y}}_\sigma = {\mathcal{M}}(TX)$. Indeed, if $x, y \in X$ and $\gamma : [0,1] \to X$ is a constant-speed geodesic joining these two points, let us define $\hat{\rho} \in {\mathcal{M}}_+(X)$ and $\hat{{\mathbf{m}}}_1, \hat{{\mathbf{m}}}_2 \in {\mathcal{M}}(TX)$ by, for $\phi \in C(X)$ and ${\mathbf{b}}\in C(X,TX)$, $$\begin{gathered}
\label{equation_controllability_continuous}
\langle \hat{\rho}, \phi \rangle = \int_0^1 \phi(\gamma(t)) ~ {\mathrm{d}}t, \ \ 
\langle \hat{{\mathbf{m}}}_1, {\mathbf{b}}\rangle =  - \int_0^1 (1-t) ~ g_{\gamma(t)}({\mathbf{b}}(\gamma(t)), \dot{\gamma}(t)) ~ {\mathrm{d}}t \\
\text{ and } 
\langle \hat{{\mathbf{m}}}_2, {\mathbf{b}}\rangle = \int_0^1 t ~ g_{\gamma(t)} ({\mathbf{b}}(\gamma(t)) , \dot{\gamma}(t)) ~ {\mathrm{d}}t.\end{gathered}$$ That is, $\hat{\rho}$ is the $\mathcal{H}^1$ Haussdorff measure restricted to $\gamma([0,1])$ (with correct scaling factor to make it a probability measure), while $\hat{{\mathbf{m}}}_1$ and $\hat{{\mathbf{m}}}_2$ have a $L^\infty$ density w.r.t. $\hat{\rho}$ (see Figure \[figure\_controllablity\] for a representation of this construction). One can check that $\nabla \cdot \hat{{\mathbf{m}}}_1 = \hat{\rho} - \delta_x$ while $\nabla \cdot \hat{{\mathbf{m}}}_2 = \hat{\rho} - \delta_y$. Moreover, the density of $\hat{{\mathbf{m}}}_1$ w.r.t. $\hat{\rho}$ at the point $\gamma(t)$ is given by $(1-t) \dot{\gamma}(t)$, whose norm is bounded by $|\dot{\gamma}(t)|_{\gamma(t)} = d_g(x,y)$ the Riemannian distance between $x$ and $y$. Hence we get $A(\hat{\rho}, \hat{{\mathbf{m}}}_1) \leqslant d_g(x,y)^2/2$, and similarly for $A( \hat{\rho}, \hat{{\mathbf{m}}}_2)$. In practice to prove \[asmp\_controllability\] one chooses $\hat{P}$ and $\hat{{\mathbf{M}}}_i$ that are mimicking $\hat{\rho}$ and $\hat{{\mathbf{m}}}_i$. We cannot use directly the sampling operators $S_{{\mathcal{X}}_\sigma}$ and $S_{{\mathcal{Y}}_\sigma}$ as the control that we get on the action thanks \[asmp\_samp\_actio\] is only valid for measures with smooth densities, while $\hat{\rho}$ and $\hat{{\mathbf{m}}}_i$ are far from such regularity.

As said earlier, the introduction of a duality structure in \[asmp\_interp\_action\] is not entirely natural. Let us state an other assumption which implies \[asmp\_interp\_action\] and which involves only the action and not its Legendre transform.

\[lemma\_interp\_action\] Assume that the following property holds:

1.  \[asmp\_interp\_action\_stronger\] *(One-sided estimate for the action)* There exists $(\varepsilon_\sigma)_{\sigma}$ tending to $0$ as $\sigma \to 0$ such that, for any $\sigma > 0$ and any $(P, {\mathbf{M}}) \in {\mathcal{X}}_\sigma \times {\mathcal{Y}}_\sigma$, $$A(R^A_{{\mathcal{X}}_\sigma} (P), R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}})) \leqslant (1 + \varepsilon_\sigma) A_\sigma(P, {\mathbf{M}}).$$

Then condition \[asmp\_interp\_action\] is satisfied.

The reason we introduced \[asmp\_interp\_action\] instead of \[asmp\_interp\_action\_stronger\] is that, for the discretization proposed in [@Gladbach2018], only \[asmp\_interp\_action\] seems to be provable.

Let ${\mathbf{b}}\in C(X,TX)$ and $P \in {\mathcal{X}}_\sigma$. We know that there exists ${\mathbf{M}}\in {\mathcal{Y}}_\sigma$ such that $A^\star_\sigma(P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}})) = \langle R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}), {\mathbf{M}}\rangle - A_\sigma(P, {\mathbf{M}})$. Using the assumption \[asmp\_interp\_action\_stronger\], we deduce that $$\begin{gathered}
A^\star_\sigma(P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}})) = \langle R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}), {\mathbf{M}}\rangle - A_\sigma(P, {\mathbf{M}}) 
\leqslant \langle R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}}) , {\mathbf{b}}\rangle - \frac{1}{1+ \varepsilon_\sigma} A( R^A_{{\mathcal{X}}_\sigma} (P), R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}})  )  \\
= (1+ \varepsilon_\sigma) \left( \langle \tilde{{\mathbf{m}}}, {\mathbf{b}}\rangle -  A( R^A_{{\mathcal{X}}_\sigma} (P),\tilde{{\mathbf{m}}} ) \right)\end{gathered}$$ with $\tilde{{\mathbf{m}}} = (1+ \varepsilon_\sigma)^{-1} R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}}) $. Using the very definition of the Legendre transform, we see that $$A^\star_\sigma(P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}})) \leqslant (1+ \varepsilon_\sigma) A^\star ( R^A_{{\mathcal{X}}_\sigma} (P), {\mathbf{b}}).$$ Now if ${\mathbf{b}}$ belongs to ${\mathcal{B}}$ a bounded set of $C^1(X,TX)$, we can easily bound the error term $\varepsilon_\sigma A^\star ( R^A_{{\mathcal{X}}_\sigma} (P), {\mathbf{b}})$ by $C \varepsilon_\sigma \| R^A_{{\mathcal{X}}_\sigma}(P) \|$ where $C$ depends only on ${\mathcal{B}}$. Getting the second estimate is easy as in the formula above ${\mathbf{b}}$ can be any continuous function. We then simply use the estimate for the continuous action: $$A^\star ( R^A_{{\mathcal{X}}_\sigma} (P), {\mathbf{b}}) = \frac{1}{2} \langle R^A_{{\mathcal{X}}_\sigma} (P), |{\mathbf{b}}|^2  \rangle \leqslant \frac{C}{2} \| R^A_{{\mathcal{X}}_\sigma} (P)  \|  \left( \sup_{x \in X} |{\mathbf{b}}(x)|_x \right)^2. \qedhere$$

Before stating the convergence result, we need to explain how to build space-time reconstruction operators from the spatial ones. These operators will be denoted ${\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}, {\mathcal{R}}^{A}_{N, {\mathcal{X}}_\sigma} : ({\mathcal{X}}_\sigma)^{N+1} \to {\mathcal{M}}([0,1] \times X)$ and ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma} : ({\mathcal{Y}}_\sigma)^N \to {\mathcal{M}}([0,1] \times TX)$. Let $(P,{\mathbf{M}}) = ((P_k)_{0 \leqslant k \leqslant N}, ({\mathbf{M}}_k)_{1 \leqslant k \leqslant N})$ a point in $({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^{N}$. For ${\mathcal{R}}^{{\mathcal{CE}}}_{N,{\mathcal{X}}_\sigma}$ we use picewise linear interpolations in time on intervals $[(k-1) \tau, k \tau]$, while for ${\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma}$ as well as ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}$, we use piecewise constant interpolation on these same intervals. Specifically for any test functions $\phi$ and ${\mathbf{b}}$ in respectively $C([0,1] \times X)$ and $C([0,1] \times X, TX)$, $$\begin{aligned}
\llangle {\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}(P), \phi \rrangle & = \sum_{k=1}^N \int_{(k-1) \tau}^{k \tau} \left\langle \frac{k \tau - t}{\tau} R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P_{k-1}) + \frac{t - (k-1) \tau}{\tau} R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P_k), \phi(t, \cdot)  \right\rangle ~ {\mathrm{d}}t, \\
\llangle {\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma}(P), \phi \rrangle & = \sum_{k=1}^N \left\langle R^A_{{\mathcal{X}}_\sigma} \left( \frac{P_{k-1} + P_k}{2} \right), \int_{(k-1)\tau}^{k \tau} \phi(t, \cdot) ~ {\mathrm{d}}t \right\rangle, \\ 
\llangle {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma }({\mathbf{M}}), {\mathbf{b}}\rrangle & = \sum_{k=1}^N \left\langle R_{{\mathcal{Y}}_\sigma}({\mathbf{M}}_k), \int_{(k-1) \tau}^{k \tau} {\mathbf{b}}(t, \cdot) ~ {\mathrm{d}}t \right\rangle.    \end{aligned}$$

Statement of the convergence result
-----------------------------------

The main result of the present work is that, properly reconstructed, the minimizers of ${\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}$ will converge to the ones of ${\mathcal{J}}_{\rho_0, \rho_1}$ in the limit $N \to + \infty, \sigma \to 0$. In the sequel, by the limit $N \to + \infty, \sigma \to 0$, it means that we can take $(N_n)_{n \in {\mathbb{N}}}$ and $(\sigma_n)_{n \in {\mathbb{N}}}$ two sequences, tending respectively to $+ \infty$ and $0$ with *no* specification of the speed at which these convergences occur, and the result holds when we study the limit $n \to + \infty$ of minimizers of ${\mathcal{J}}^{N_n, \sigma_n}_{S_{{\mathcal{X}}_{\sigma_n}}(\rho_0), S_{{\mathcal{X}}_{\sigma_n}}(\rho_1)}$.

\[theo\_main\] Let $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_\sigma$ a family of finite dimensional models of dynamical optimal transport which is adapted to the Riemannian manifold $(X,g)$ in the sense of Definition \[definition\_adapated\], with reconstruction and sampling operators $(R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}, R_{{\mathcal{Y}}_\sigma}, S_{{\mathcal{X}}_\sigma}, S_{{\mathcal{Y}}_\sigma})_\sigma$.

Let $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ sharing the same total mass be given. For any $N$ and $\sigma > 0$, let $(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) \in ({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^N$ be minimizers of ${\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}$.

Then, in the limit $N \to + \infty$ and $\sigma \to 0$, up to the extraction of a subsequence, ${\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$ and ${\mathcal{R}}^{A}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$ will converge weakly to the same limit $\rho \in {\mathcal{M}}([0,1] \times X)$, and ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma})$ will converge weakly to ${\mathbf{m}}\in {\mathcal{M}}([0,1] \times TX)$, in addition the pair $(\rho, {\mathbf{m}})$ is a minimizer of ${\mathcal{J}}_{\rho_0, \rho_1}$ and $$\lim_{N \to + \infty, \sigma \to 0} {\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)} ( P^{N, \sigma}, {\mathbf{M}}^{N, \sigma} ) = {\mathcal{J}}_{\rho_0, \rho_1} (\rho, {\mathbf{m}}).$$

In particular, if the minimizer of ${\mathcal{J}}_{\rho_0, \rho_1}$ is unique (which can be guaranteed in some situation, see Theorem \[theorem\_J\_W2\]), we can remove the “up to the extraction of a subsequence” in the statement.

We emphasize that the reconstructed densities ${\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$, ${\mathcal{R}}^{A}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$ and the reconstructed momentum ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma})$ do not necessarily satisfy a continuity equation, nor have a finite action. Only their weak limit do.

Such a result is in fact obtained by combining two sharper results, which look like a $\Gamma$-convergence result.

\[theo\_GammaLiminf\] Let $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_\sigma$ and $(R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}, R_{{\mathcal{Y}}_\sigma}, S_{{\mathcal{X}}_\sigma}, S_{{\mathcal{Y}}_\sigma})_\sigma$ be as in Theorem \[theo\_main\].

Let $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ sharing the same total mass be given. For any $N$ and $\sigma > 0$, let $(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) \in ({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^N$ such that $$\sup_{N, \sigma} {\mathcal{J}}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}^{N, \sigma}(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) < + \infty.$$ Then, in the limit $N \to + \infty$ and $\sigma \to 0$, up to the extraction of a subsequence, ${\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$ and ${\mathcal{R}}^{A}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$ will converge weakly to the same limit $\rho \in {\mathcal{M}}([0,1] \times X)$, and ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma})$ will converge weakly to ${\mathbf{m}}\in {\mathcal{M}}([0,1] \times TX)$, in addition $${\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) \leqslant \liminf_{N \to + \infty, \sigma \to 0} {\mathcal{J}}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}^{N, \sigma}(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}).$$

\[theo\_GammaLimsup\] Let $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_\sigma$ and $(R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}, R_{{\mathcal{Y}}_\sigma}, S_{{\mathcal{X}}_\sigma}, S_{{\mathcal{Y}}_\sigma})_\sigma$ be as in Theorem \[theo\_main\].

Let $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ sharing the same total mass be given, and let $(\rho, {\mathbf{m}}) \in {\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ be such that ${\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) < + \infty$. For any $N$ and $\sigma > 0$, we can build $(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) \in ({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^N$ such that $$\limsup_{N \to + \infty, \sigma \to 0}  {\mathcal{J}}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}^{N, \sigma}(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) \leqslant {\mathcal{J}}_{\rho_0, \rho_1} (\rho, {\mathbf{m}}).$$

It is of course very simple to obtain Theorem \[theo\_main\] from Theorems \[theo\_GammaLiminf\] and \[theo\_GammaLimsup\], the delicate part being the proof of the latter.

It would be tempting to write this result as $\Gamma$-convergence one. However, in Theorem \[theo\_GammaLimsup\], there is no guarantee that the sampled momentum ${\mathbf{M}}^{N, \sigma}$ converges to ${\mathbf{m}}$, especially if the density of ${\mathbf{m}}$ w.r.t. $\rho$ does not belong to the closure of vector fields which are gradient, see Remark \[rmk\_m\_not\_optimal\]. Moreover, we think that a lesson of our approach is that there is no canonical way to inject ${\mathcal{X}}_\sigma$ and ${\mathcal{Y}}_\sigma$ into ${\mathcal{M}}(X)$ and ${\mathcal{M}}(TX)$, only different ways which are all useful and asymptotically equivalent.

In the sequel, we first prove Theorems \[theo\_GammaLiminf\] and \[theo\_GammaLimsup\], and then show how previously proposed discretizations satisfy the assumptions of Definition \[definition\_adapated\]: this is the object of Section \[section\_proof\] and Section \[section\_examples\] respectively. Given the way our statements are phrased, these two sections are entirely independent.

Proof of the convergence result {#section_proof}
===============================

The goal of this section is to prove Theorems \[theo\_GammaLiminf\] as well as \[theo\_GammaLimsup\]. As explained in the introduction, the main issue is the discontinuity of the function $(\rho, {\mathbf{m}}) \mapsto |{\mathbf{m}}|^2 /(2 \rho)$ at the point $(0,0)$. As this function is nevertheless lower semi-continuous, Theorem \[theo\_GammaLiminf\] is still quite straightforward. However, for the proof Theorem \[theo\_GammaLimsup\] we rely on careful approximation arguments.

Proof of Theorem \[theo\_GammaLiminf\]
--------------------------------------

Recall that we assume that $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ sharing the same total mass are given. For any $N, \sigma$, we have $(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) \in ({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^N$ such that $$\sup_{N, \sigma} {\mathcal{J}}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}^{N, \sigma}(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) < + \infty.$$

*First step: uniform estimates on the mass of the measures*. First, we need to prove that, as measures, ${\mathcal{R}}^{{\mathcal{CE}}}_{N, \sigma}(P^{N, \sigma})$, ${\mathcal{R}}^{A}_{N, \sigma}(P^{N, \sigma})$ and ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma})$ are uniformly bounded. To that extent, we rely both on the continuity equation and the one-sided estimate on the action \[asmp\_interp\_action\]. Given a particular discretization (like the ones in Section \[section\_examples\]), a more direct proof may exist, but we want to stress out that it can already be proved only with the assumptions in Definition \[definition\_adapated\] (even though it leads to a more technical estimates).

Let us use the continuity equation to get an estimate on $\| {\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma}) \|$ and $\|{\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma}) \|$. To that extent, let us take the discrete continuity equation, first apply $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ and then take the duality product with $\phi$ the constant function equal to $1$. We get, for any $k \in \{ 1,2, \ldots, N \}$, $$\langle R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{k}), \phi \rangle = \langle R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{k-1}), \phi \rangle + \tau \langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ {\mathrm{Div}}_\sigma) ({\mathbf{M}}^{N, \sigma}_k), \phi \rangle.$$ As $P^{N, \sigma}_k \in {\mathcal{X}}_{\sigma,+}$ we know that $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{k}) \in {\mathcal{M}}_+(X)$ hence taking the duality product against $\phi$ the constant function equal to $1$ yields the mass of this measure. On the other hand, because of the assumption \[asmp\_interp\_der\] and as $\nabla \phi = 0$, up to an error $\varepsilon_\sigma$ which goes to $0$ as $\sigma \to 0$, $$\left| \|R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{k}) \| - \|R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{k-1})\| \right| \leqslant \varepsilon_\sigma \tau \|R_{{\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}_k) \|.$$ Moreover, we know that $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{0}) = (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ S_{{\mathcal{X}}_\sigma})(\rho_0)$ converges to $\rho_0$, which has a given mass. Summing the estimates above, we see that for any $k \in \{ 0,1, \ldots, N \}$, $$\left| \|R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{k}) \| - \| \rho_0 \| \right| \leqslant \varepsilon_\sigma + \varepsilon_\sigma \tau \sum_{l=1}^N \|R_{{\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}_k) \| =  \varepsilon_\sigma +  \varepsilon_\sigma \| {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}) \|.$$ Using \[asmp\_interps\] tested against the constant function, we deduce that $\| R^{A}_{X_\sigma} (P^{N, \sigma}_k) \| \leqslant \| R^{{\mathcal{CE}}}_{X_\sigma} (P^{N, \sigma}_k) \| + \varepsilon_\sigma$ for all $k \in \{ 0,1, \ldots N \}$. Summing these estimates over $k$, we get that $$\label{zz_aux_2}
\|{\mathcal{R}}^{A}_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma}) \| + \| {\mathcal{R}}^{{\mathcal{CE}}}_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma}) \|  \leqslant 2 \| \rho_0 \| + \varepsilon_\sigma + \varepsilon_\sigma \| {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}) \|,$$ where $\varepsilon_\sigma$ may have changed from one line to another but keeps the property of tending to $0$ as $\sigma \to 0$.

Now we have to control the norm of the reconstructed momentum $\| {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}) \|$. Let ${\mathbf{b}}\in C([0,1] \times X,TX)$ such that $|{\mathbf{b}}(t,x)| \leqslant 1$ for all $(t,x) \in [0,1] \times X$. Not to overburden the equations, let us introduce for $k \in \{ 1,2, \ldots, N \}$ the functions ${\mathbf{b}}_k \in C(X,TX)$ defined by $${\mathbf{b}}_k(x) = \frac{1}{\tau} \int_{(k-1) \tau}^{k \tau} {\mathbf{b}}(t,x) ~ {\mathrm{d}}t,$$ they also satisfy the pointwise estimate $|{\mathbf{b}}_k(x)|_x \leqslant 1$ for all $k$ and $x$. Using estimate followed by Cauchy-Schwarz, we get: $$\begin{aligned}
\llangle  {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}), {\mathbf{b}}\rrangle 
&= \sum_{k=1}^N \tau \langle R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}}^{N, \sigma}_k), {\mathbf{b}}_k \rangle \\
& \leqslant \sum_{k=1}^N \tau \sqrt{ A_\sigma( P^{N, \sigma}_k, {\mathbf{M}}^{N, \sigma}_k  )} \sqrt{A^\star_\sigma(  P^{N, \sigma}_k, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}_k)   )   }  \\
& \leqslant \sqrt{ \sum_{k=1}^N \tau A_\sigma( P^{N, \sigma}_k, {\mathbf{M}}^{N, \sigma}_k  )  } \sqrt{ \sum_{k=1}^N \tau A^\star_\sigma(  P^{N, \sigma}_k, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}_k)   )  }.\end{aligned}$$ The first term in the product is controlled by ${\mathcal{J}}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}^{N, \sigma}(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma})$ which is assumed to be uniformly bounded in $N, \sigma$, while for the second we use the assumption \[asmp\_interp\_action\] and the pointwise estimates on ${\mathbf{b}}_k$. Denoting by $C > 0$ a constant which may change from one equation to another, we end up with $$\llangle {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}), {\mathbf{b}}\rrangle \leqslant \sqrt{C} \sqrt{ \frac{C}{2} \sum_{k=1}^N \tau \| R^A_{{\mathcal{X}}_\sigma} (P^{N, \sigma}_k) \|   } \leqslant C \sqrt{ \| {\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma} (P^{N, \sigma}) \|  }.$$ Now we can take the supremum over all ${\mathbf{b}}$ and use that we got from the discrete continuity equation to end up with $$\| {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}) \| \leqslant C \sqrt{C + C  \| {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}) \| },$$ with $C > 0$ a constant independent on $N$ and $\sigma$.

It allows us to conclude that $\| {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}) \| $ is bounded uniformly in $N$ and $\sigma$. Using , we see that $\| {\mathcal{R}}^{A}_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma}) \|$ and $ \| {\mathcal{R}}^{{\mathcal{CE}}}_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma}) \|$ are also bounded uniformly in $N$ and $\sigma$. Hence, up to extraction, ${\mathcal{R}}^{{\mathcal{CE}}}_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma})$, ${\mathcal{R}}^A_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma})$ and ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma})$ converge weakly in respectively ${\mathcal{M}}([0,1] \times X), {\mathcal{M}}([0,1] \times X)$ and ${\mathcal{M}}([0,1] \times TX)$ to some limits that we call respectively $\rho, \tilde{\rho}$ and ${\mathbf{m}}$.

*Second step: the limits $\rho$ and $\tilde{\rho}$ are the same*. Thanks to the assumption \[asmp\_interps\], this is easy to see. Let us take $\phi \in C^1([0,1] \times X)$ a smooth function. In particular, $\{ \phi(t, \cdot) \ : \ t \in [0,1] \}$ is a bounded set of $C^1(X,TX)$. Moreover, up to an error controlled by $C \tau$ in $L^\infty$ norm, $\phi(t, \cdot)$ can be replaced by $\phi(k \tau, \cdot)$ for all $t \in [(k-1) \tau, k \tau]$. Given the uniform bound on the mass of ${\mathcal{R}}^{{\mathcal{CE}}}_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma})$ and ${\mathcal{R}}^A_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma})$, we conclude that $$\left| \left\llangle ( {\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma} - {\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma} )(P^{N, \sigma}), \phi \right\rrangle \right| \leqslant \sum_{k=1}^N \tau \left| \left\langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} - R^{A}_{{\mathcal{X}}_\sigma}) \left( \frac{P_{k-1} + P_k}{2} \right), \phi(k \tau, \cdot) \right\rangle \right| + C \tau.$$ Now, given \[asmp\_interps\] and the uniform bound on $\| {\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma}) \| $ we see easily that, for some error $\varepsilon_\sigma$ which goes to $0$ as $\sigma \to 0$, $$\left| \left\llangle ( {\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma} - {\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma} )(P^{N, \sigma}), \phi \right\rrangle \right| \leqslant \varepsilon_\sigma + C \tau.$$ Sending $N \to + \infty$ (hence $\tau \to 0$) and $\sigma \to 0$, we end up with $\langle \rho - \tilde{\rho}, \phi \rangle$. As $\phi$ was an arbitrary smooth function, it gives us the equality between $\rho$ and $\tilde{\rho}$.

*Third step: passing to the limit in the continuity equation*. Let $\phi \in C^2([0,1] \times X)$ a smooth test function. We will test the continuity equation against $({\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma}), {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}))$ and pass to the limit. Indeed, an integration by parts followed by the use of the discrete continuity equation leads to $$\begin{aligned}
& \llangle  {\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma} (P^{N, \sigma}), {\partial}_t \phi \rrangle  
= \sum_{k=1}^N \int_{(k-1) \tau}^{k \tau} \left\langle \frac{k \tau - t}{\tau} R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{k-1}) + \frac{t - (k-1) \tau}{\tau} R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_k), {\partial}_t \phi(t, \cdot)  \right\rangle \, {\mathrm{d}}t \\
&= \langle R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_N), \phi(1, \cdot) \rangle - \langle R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_0), \phi(0, \cdot) \rangle - \sum_{k=1}^N \left\langle R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \left( \frac{P^{N, \sigma}_k - P^{N, \sigma}_{k-1}}{\tau} \right), \int_{(k-1) \tau}^{k \tau} \phi(t, \cdot) \,{\mathrm{d}}t  \right\rangle \\
&=\langle ( R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ S_{{\mathcal{X}}_\sigma})(\rho_1), \phi(1, \cdot) \rangle - \langle ( R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ S_{{\mathcal{X}}_\sigma})(\rho_0), \phi(0, \cdot) \rangle + \sum_{k=1}^N \left\langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ {\mathrm{Div}}_\sigma)( {\mathbf{M}}^{N, \sigma}_k ), \int_{(k-1) \tau}^{k \tau} \phi(t, \cdot) \, {\mathrm{d}}t  \right\rangle.\end{aligned}$$ As $\{ \phi(t, \cdot) \ : \ t \in [0,1] \}$ is a bounded set of $C^2(X)$, thanks to \[asmp\_interp\_der\], with an error $\varepsilon_\sigma$ tending to $0$ as $\sigma \to 0$, we can write that $$\begin{gathered}
\left| \underbrace{\sum_{k=1}^N \left\langle R_{{\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}_k), \int_{(k-1) \tau}^{k \tau} \nabla \phi(t, \cdot) \, {\mathrm{d}}t \right\rangle}_{ = \llangle {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}), \nabla \phi \rrangle }  + \sum_{k=1}^N \left\langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ {\mathrm{Div}}_\sigma)( {\mathbf{M}}^{N, \sigma}_k ), \int_{(k-1) \tau}^{k \tau} \phi(t, \cdot) \, {\mathrm{d}}t  \right\rangle \right| \\
\leqslant \varepsilon_\sigma \| {\mathcal{R}}_{N,{\mathcal{Y}}_\sigma} ({\mathbf{M}}^{N, \sigma}) \|.\end{gathered}$$ Plugging back this information, we see that $$\begin{gathered}
\left| \llangle {\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma} (P^{N, \sigma}), {\partial}_t \phi \rrangle + \llangle {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}), \nabla \phi \rrangle - \langle ( R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ S_{{\mathcal{X}}_\sigma})(\rho_1), \phi(1, \cdot) \rangle + \langle ( R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ S_{{\mathcal{X}}_\sigma})(\rho_0), \phi(0, \cdot) \rangle   \right| \\
\leqslant  \varepsilon_\sigma \| {\mathcal{R}}_{N,{\mathcal{Y}}_\sigma} ({\mathbf{M}}^{N, \sigma}) \|.\end{gathered}$$ As $\| {\mathcal{R}}_{N,{\mathcal{Y}}_\sigma} ({\mathbf{M}}^{N, \sigma}) \|$ is uniformly bounded in $N$ and $\sigma$, sending $N \to + \infty$ and $\sigma \to 0$, with the help of \[asmp\_interp\_sampling\] to handle the boundary terms, we end up with $$\llangle \rho, {\partial}_t \phi \rrangle + \llangle {\mathbf{m}}, \nabla \phi \rrangle - \langle \rho_1, \phi(1, \cdot) \rangle + \langle \rho_0, \phi(0, \cdot) \rangle = 0.$$ As $\phi$ is an arbitrary smooth function, it means that $(\rho, {\mathbf{m}}) \in {\mathcal{CE}}(\rho_0, \rho_1)$, that is it satisfies the continuity equation with boundary conditions $(\rho_0, \rho_1)$.

*Fourth step: passing to the limit the action*. In this abstract setting, this is simple. Even though \[asmp\_interp\_action\_stronger\] implies \[asmp\_interp\_action\], let us first show the result if \[asmp\_interp\_action\_stronger\] holds as it is a one-line estimate. Indeed, given \[asmp\_interp\_action\] and the way ${\mathcal{R}}^{A}_{N, {\mathcal{X}}_\sigma}, {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}$ are defined, $$\begin{gathered}
{\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}( P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}  ) = \tau \sum_{k=1}^{N} A_\sigma \left( \frac{P^{N, \sigma}_{k-1} + P^{N, \sigma}_{k}}{2}, {\mathbf{M}}^{N, \sigma}_k \right) \\ 
\geqslant \frac{1}{1+ \varepsilon_\sigma} \tau \sum_{k=1}^N A \left( R^{A}_{\sigma} \left( \frac{P^{N, \sigma}_{k-1} + P^{N, \sigma}_{k}}{2} \right), R_{{\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}_k)  \right) 
= \frac{1}{1+ \varepsilon_\sigma} {\mathcal{A}}( {\mathcal{R}}^{A}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma}), {\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma}) ).\end{gathered}$$ As the action ${\mathcal{A}}$ is lower semi-continuous on ${\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$, we can easily pass to the limit.

On the other hand, let us just assume \[asmp\_interp\_action\]. Let us fix $\eta > 0$, we take ${\mathbf{b}}\in C(^1[0,1] \times X, TX)$ smooth such that $${\mathcal{A}}(\tilde{\rho}, {\mathbf{m}}) \leqslant \llangle {\mathbf{m}}, {\mathbf{b}}\rrangle - {\mathcal{A}}^\star( \tilde{\rho}, {\mathbf{b}}) + \eta.$$ Now, by continuity of ${\mathcal{A}}^\star(\cdot, {\mathbf{b}})$ (w.r.t. weak convergence) and the weak convergence of ${\mathcal{R}}_{N,{\mathcal{Y}}_\sigma}( {\mathbf{M}}^{N, \sigma} )$ and ${\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$, one sees that $$\begin{gathered}
\llangle {\mathbf{m}}, {\mathbf{b}}\rrangle - {\mathcal{A}}^\star( \tilde{\rho}, {\mathbf{b}})
= \lim_{N \to + \infty, \sigma \to 0} \left(  \llangle {\mathcal{R}}_{N,{\mathcal{Y}}_\sigma}( {\mathbf{M}}^{N, \sigma} ), {\mathbf{b}}\rrangle - {\mathcal{A}}^\star( {\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma}), {\mathbf{b}}) \right) \\
=  \lim_{N \to + \infty, \sigma \to 0}  \sum_{k=1}^N \int_{(k-1) \tau}^{k \tau} \left(  \langle {\mathbf{M}}^{N, \sigma}_k, R_{{\mathcal{Y}}_\sigma}^\top ( {\mathbf{b}}(t, \cdot) ) \rangle - A^\star \left( R^A_{{\mathcal{X}}_\sigma} \left( \frac{P_{k-1}^{N, \sigma} + P_{k}^{N, \sigma}}{2} \right)  , {\mathbf{b}}(t, \cdot)  \right) \right) \,{\mathrm{d}}t. \end{gathered}$$ As $\{ {\mathbf{b}}(t, \cdot), \ t \in [0,1] \}$ is a bounded set of $C^1(X,TX)$, we can use \[asmp\_interp\_action\] to write that $$\begin{gathered}
\llangle {\mathbf{m}},  {\mathbf{b}}\rrangle - {\mathcal{A}}^\star( \tilde{\rho}, {\mathbf{b}})  
\leqslant  \liminf_{N \to + \infty, \sigma \to 0}  \sum_{k=1}^N \int_{(k-1) \tau}^{k \tau} \Bigg( \langle {\mathbf{M}}^{N, \sigma}_k, R_{{\mathcal{Y}}_\sigma}^\top ( {\mathbf{b}}(t, \cdot) ) \rangle   -  A^\star_\sigma \Bigg(  \frac{P_{k-1}^{N, \sigma} + P_{k}^{N, \sigma}}{2}  ,   R_{{\mathcal{Y}}_\sigma}^\top ( {\mathbf{b}}(t, \cdot))  \Bigg) \\  
+ \frac{\varepsilon_\sigma}{2} \| R^A_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_{k-1}) + R^A_{{\mathcal{X}}_\sigma}(P^{N, \sigma}_k) \| \Bigg) \, {\mathrm{d}}t  \end{gathered}$$ Then, using the definition of the Legendre transform, we conclude that $$\begin{gathered}
{\mathcal{A}}(\tilde{\rho} , {\mathbf{m}}) - \eta \leqslant \liminf_{N \to + \infty, \sigma \to 0}  \left[ \left\{ \tau \sum_{k=1}^N  A_\sigma \left( \frac{P_{k-1}^{N, \sigma} + P_{k}^{N, \sigma}}{2}, {\mathbf{M}}_k^{N, \sigma}  \right) \right\} + \varepsilon_\sigma \| {\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma}(P^{N,\sigma}) \| \right]    \\ 
\leqslant \liminf_{N \to + \infty, \sigma \to 0}  {\mathcal{J}}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}^{N, \sigma}(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}). \end{gathered}$$ Notice that the error term has disappeared as $\| {\mathcal{R}}^A_{N, {\mathcal{X}}_\sigma}(P^{N,\sigma}) \|$ is uniformly bounded in $N$ and $\sigma$. As $\eta$ can be taken arbitrary small, we end up with $$\label{equation_GammaLiminf_A}
\liminf_{N \to + \infty, \sigma \to 0} {\mathcal{J}}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}^{N, \sigma}(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) \geqslant {\mathcal{A}}(\tilde{\rho}, {\mathbf{m}}).$$ We have seen $\tilde{\rho} = \rho$ and $(\rho, {\mathbf{m}}) \in {\mathcal{CE}}(\rho_0, \rho_1)$: it makes the right hand side equal to ${\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}})$.

As one can see in the proof of the first and third steps, in fact $(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma})$ do not need to satisfy exactly the discrete continuity equation, one could allow for some leeway. As pointed out in [@Carrillo2019 Remark 2], allowing for some leeway can help to speed up the numerical computation of a minimum of ${\mathcal{J}}^{N, \sigma}$.

Proof of Theorem \[theo\_GammaLimsup\]
--------------------------------------

Such a proof is more involved as it relies on a careful regularization procedure, as well as the use of the controllability assumption \[asmp\_controllability\] to handle what is happening near the temporal boundaries. In the sequel, by an abuse of notations, we will identify a measure with its density w.r.t. the volume measure. Let us take $\rho_0, \rho_1$ sharing the same total mass and $(\rho, {\mathbf{m}}) \in {\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ such that ${\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) < + \infty$. To prove our result, we will first regularize $(\rho, {\mathbf{m}})$ and then sample via $S_{{\mathcal{X}}_\sigma}$ and $S_{{\mathcal{Y}}_\sigma}$. Because of the regularization we lose the temporal boundary conditions: this is to remedy to this problem that we use the assumption \[asmp\_controllability\], see Proposition \[prop\_controllability\] below.

We recall that $W_2$ denotes the quadratic Wasserstein distance, see Appendix \[section\_wasserstein\].

*First tool: regularization of $(\rho, {\mathbf{m}})$*. The first tool is to regularize the continuous pair $(\rho, {\mathbf{m}})$. This is object of the following proposition.

\[prop\_regularization\] Let $(\rho, {\mathbf{m}}) \in {\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ and $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ such that ${\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) < + \infty$. For any $\eta > 0$, there exists $(\tilde{\rho}, \tilde{{\mathbf{m}}}) \in {\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ such that:

1.  The densities of $\tilde{\rho}$ and $\tilde{{\mathbf{m}}}$ w.r.t. ${\mathrm{d}}t \otimes {\mathrm{d}}x$ (still denoted by $\tilde{\rho}$ and $\tilde{{\mathbf{m}}}$) are smooth in the sense that $\tilde{\rho} \in C^1([0,1] \times X)$ and is bounded from below uniformly on $[0,1] \times X$ by a strictly positive constant; while $\tilde{{\mathbf{m}}} \in C^1([0,1] \times X, TX)$.

2.  The pair $(\tilde{\rho}, \tilde{{\mathbf{m}}})$ belongs to ${\mathcal{CE}}(\tilde{\rho}(0, \cdot) , \tilde{\rho}(1, \cdot))$ and $W_2(\tilde{\rho}(0, \cdot) , \rho_0) \leqslant \eta$, as well as $W_2(\tilde{\rho}(1, \cdot), \rho_1) \leqslant \eta$.

3.  The following estimate on the action holds: $${\mathcal{A}}(\tilde{\rho}, \tilde{{\mathbf{m}}}) \leqslant (1+ \eta) {\mathcal{A}}(\rho, {\mathbf{m}}).$$

Though we have not found this result phrased like this in the literature, closely related ones are available. As it can be considered as standard and it is not the core of the argument, we delay the proof of such a result until Appendix \[section\_regularization\_curves\]. Notice that we have not mentioned in the Proposition whether $\tilde{\rho}$ and $\tilde{{\mathbf{m}}}$ are closed to $\rho$ and ${\mathbf{m}}$ as it is irrelevant for the rest of the analysis, more is said in Remark \[rmk\_m\_not\_optimal\].

*Second tool: controllability*. The issue with the regularization procedure is that it does not keep the boundary values. Assumption \[asmp\_controllability\] will help us to do a little surgery near the temporal boundaries.

\[prop\_controllability\] There exists a continuous function $\hat{\omega} : [0, + \infty) \to [0, + \infty)$ independent on $N$ and $\sigma$, such that $\hat{\omega}(0) = 0$, and an error $\varepsilon_{N, \sigma}$ depending only on $N, \sigma$ and going to $0$ in the limit $N \to + \infty, \sigma \to 0$ such that for any $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ sharing the same total mass, there holds $$\min_{({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^N} {\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)} \ \leqslant  \hat{\omega}(W_2(\rho_0, \rho_1)) + \varepsilon_{N, \sigma} \| \rho_0 \|.$$

Here the key point is that $\rho_0, \rho_1$ are *arbitrary* positive measures and there is no condition on the ratio between the temporal and spatial step sizes.

For simplicity (the only effect would be the appearance of errors of order $\tau$ which can be absorbed in $\varepsilon_{N, \sigma}$), we assume that $N$ is even. Let us also fix $\sigma > 0$.

We first consider the case where $\rho_0 = \delta_x$ and $\rho_1 = \delta_y$ with $x,y \in X$. Given assumption \[asmp\_controllability\], there exists $\hat{{\mathbf{M}}}_1, \hat{{\mathbf{M}}}_2 \in {\mathcal{Y}}_\sigma$ and $\hat{P} \in {\mathcal{X}}_{\sigma,+}$ satisfying . Now let us take $\chi : [0,1] \to [0,1]$ the continuous function defined by $$\chi(t) = \begin{cases}
4 t^2 & \text{if } t \leqslant 1/2, \\ 
4 (1-t)^2 & \text{if } t \geqslant 1/2.
\end{cases}$$ In particular, and it was chosen for that, $\int_0^1 \dot{\chi}^2/\chi < + \infty$ while $\chi(0) = \chi(1) = 0$ and $\chi(1/2) = 1$. We define $((P_k)_{0 \leqslant k \leqslant N}, ({\mathbf{M}}_k)_{1 \leqslant k \leqslant N}) \in {\mathcal{X}}_\sigma^{N+1} \times {\mathcal{Y}}_\sigma^N$ by $$\begin{gathered}
P_k := \begin{cases} 
(1 - \chi(k \tau)) S_{{\mathcal{X}}_\sigma}(\delta_x) + \chi(k \tau) \hat{P} & \text{if } k \leqslant N/2, \\
(1-\chi(k \tau)) S_{{\mathcal{X}}_\sigma}(\delta_y) + \chi(k \tau) \hat{P} & \text{if } k \geqslant N/2, \\
\end{cases} \\
\text{ and } {\mathbf{M}}_k = 
\begin{cases}
{\displaystyle{ \frac{ ( \chi(k \tau) - \chi((k-1) \tau))  }{\tau}}} \hat{{\mathbf{M}}}_1  & \text{if } k \leqslant N/2, \\
{\displaystyle{ \frac{ ( \chi(k \tau) - \chi((k-1) \tau))  }{\tau}}} \hat{{\mathbf{M}}}_2 & \text{if } k > N/2.
\end{cases}\end{gathered}$$ Notice that, for $k \leqslant N/2$, $P_k$ is just a convex combination of $S_{{\mathcal{X}}_\sigma}(\delta_x)$ and $\hat{P}$, while all the ${\mathbf{M}}_k$ are proportional to $\hat{{\mathbf{M}}}_1$; a similar situation occurs for $k \geqslant N/2$. By construction and in particular thanks to , it is quite straightforward to see that $((P_k)_{0 \leqslant k \leqslant N}, ({\mathbf{M}}_k)_{1 \leqslant k \leqslant N})$ satisfies the discrete continuity equation. We need to estimate the action: for the sake of the exposition, let us just estimate on the part $k \leqslant N/2$, the other part being completely symmetric. We use first that $S_{{\mathcal{X}}_\sigma}(\delta_x) \in {\mathcal{X}}_{\sigma,+}$ and the action only decreases if we add an element of ${\mathcal{X}}_{\sigma,+}$ to its first variable, and then the $(-1,2)$ homogeneity of the action $A_\sigma$: $$\begin{aligned}
& {\mathcal{A}}_{N, \sigma} ((P_k)_{0 \leqslant k \leqslant N/2}, ({\mathbf{M}}_k)_{1 \leqslant k \leqslant N/2}) 
\\ 
& = \tau \sum_{k=1}^{N/2} A_\sigma \left( \left[ 1 - \frac{\chi((k-1) \tau) + \chi(k \tau)}{2} \right] S_{{\mathcal{X}}_\sigma}(\delta_x) + \frac{\chi((k-1) \tau) + \chi(k \tau)}{2} \hat{P}, \frac{ ( \chi(k \tau) - \chi((k-1) \tau))  }{\tau} \hat{{\mathbf{M}}}_1  \right) \\
& \leqslant  \tau \sum_{k=1}^{N/2} A_\sigma \left( \frac{\chi((k-1) \tau) + \chi(k \tau)}{2} \hat{P}, \frac{ ( \chi(k \tau) - \chi((k-1) \tau))  }{\tau} \hat{{\mathbf{M}}}_1  \right) \\
& =  A_\sigma(\hat{P}, \hat{{\mathbf{M}}}_1) \sum_{k=1}^{N/2} \frac{(\chi(k \tau) - \chi((k-1) \tau))^2}{\tau (\chi(k \tau) + \chi((k-1) \tau))}.\end{aligned}$$ By assumption we can control $A_\sigma(\hat{P}, {\mathbf{M}}_1)$ by $\omega(d_g(x,y)) + \varepsilon_\sigma$. About the remaining sum, notice that it is nothing else than a discretization of $\int_0^1 \dot{\chi}^2/\chi$ which is finite. More precisely given the explicit expression of $\chi$, $$\frac{(\chi(k \tau) - \chi((k-1) \tau))^2}{\tau (\chi(k \tau) + \chi((k-1) \tau))} = 4 \frac{\tau (2k-1)^2}{k^2 + (k-1)^2} \leqslant 16 \tau$$ hence the sum will be bounded independently on $N$. A similar computation can be performed exactly in the same way for $k \geqslant N/2$. In short the sum is bounded independently on $N$, which translates in $${\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\delta_x), S_{{\mathcal{X}}_\sigma}(\delta_y)}(P, {\mathbf{M}}) \leqslant C \omega(d_g(x,y)) + \varepsilon_{N,\sigma}.$$ and it is exactly what we wanted to prove.

Next, let us consider the general case: we take $\rho_0, \rho_1 \in {\mathcal{M}}(X)$ with the same total mass. Let $\pi \in {\mathcal{M}}(X \times X)$ an optimal transport plan between them (see Appendix \[section\_wasserstein\]). For each $(x,y) \in X \times X$, we consider the pair $(P^{xy}, {\mathbf{M}}^{xy})$ built as above which has an energy ${\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\delta_x), S_{{\mathcal{X}}_\sigma}(\delta_y)}(P^{xy}, {\mathbf{M}}^{xy})$ bounded by $C \omega(d_g(x,y)) + \varepsilon_{N,\sigma}$. Then we simply set $$(P, {\mathbf{M}}) := \iint_{X \times X} (P^{xy}, {\mathbf{M}}^{xy}) \, \pi({\mathrm{d}}x, {\mathrm{d}}y) \in ({\mathcal{X}}^+_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^N.$$ By linearity of the discrete continuity equation, the latter is still satisfied with boundary conditions $(S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1))$ as $\pi$ has appropriate marginals. Moreover, by convexity of the action $A_\sigma$ and its $1$-homogeneity, $$\begin{gathered}
{\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}(P, {\mathbf{M}}) \leqslant \iint_{X \times X} {\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\delta_x), S_{{\mathcal{X}}_\sigma}(\delta_y)}(P^{xy}, {\mathbf{M}}^{xy}) \, \pi({\mathrm{d}}x, {\mathrm{d}}y) \\ 
\leqslant \iint_{X \times X} ( \omega( d_g(x,y) ) + \varepsilon_{N, \sigma}) \, \pi({\mathrm{d}}x, {\mathrm{d}}y).  \end{gathered}$$ Now, as $\omega$ is continuous and $\omega(0) = 0$, it is clear that if $\iint d_g(x,y)^2 \, \pi({\mathrm{d}}x, {\mathrm{d}}y)$ goes to $0$, so does $\iint_{X \times X} \omega( d_g(x,y) ) \, \pi({\mathrm{d}}x, {\mathrm{d}}y)$, hence the right hand side can be written $\hat{\omega}(W_2(\rho_0, \rho_1))$ with $\hat{\omega}$ continuous and $\hat{\omega}(0) = 0$, up to an error $\varepsilon_{N, \sigma} \| \pi \| = \varepsilon_{N, \sigma} \| \rho_0 \|$.

(0,0) – (14,0) ; (0,0) circle (0.1) ; (14,0) circle (0.1) ;

(3,0) node\[above\][$(\rho, {\mathbf{m}})$]{} ; (0,0) node\[left\][$\rho_0$]{} ; (14,0) node\[right\][$\rho_1$]{} ;

(4, -0.3) – (4,-1.7) ; (4.3,-1) node\[right\] ; (3,-2) node\[above\][($\tilde{\rho}, \tilde{{\mathbf{m}}})$]{} ;

(0,-2) – (14,-2) ; (0.1,-1.9) rectangle (-0.1, -2.1) ; (14.1,-1.9) rectangle (13.9, -2.1) ;

(0,-2) node\[left\][$\tilde{\rho}_0$]{} ; (14,-2) node\[right\][$\tilde{\rho}_1$]{} ;

(4, -2.3) – (4,-3.7) ; (4.3,-3) node\[right\] ; (2,-4) – (12,-4) ; (2.1,-3.9) rectangle (1.9, -4.1) ; (12.1,-3.9) rectangle (11.9, -4.1) ;

(2,-4) node\[above\][$\tilde{\rho}_0$]{} ; (12,-4) node\[above\][$\tilde{\rho}_1$]{} ;

(0,-4) – (2,-4) ; (11,-4) – (14,-4) ;

(0,-4) circle (0.1) ; (14,-4) circle (0.1) ;

(0,-3.9) node\[left\][$\rho_0$]{} ; (14,-3.9) node\[right\][$\rho_1$]{} ;

(0,-6) circle (0.1) ; (14,-6) circle (0.1) ; in [1,2,...,27]{} [([0.5]{},-6) circle (0.05) ;]{}

(0,-6) node\[left\][$S_{{\mathcal{X}}_\sigma}(\rho_0)$]{} ; (14,-6) node\[right\][$S_{{\mathcal{X}}_\sigma}(\rho_1)$]{} ;

(2.1,-5.9) rectangle (1.9, -6.1) ; (12.1,-5.9) rectangle (11.9, -6.1) ;

(2,-6.1) node\[below\][$S_{{\mathcal{X}}_\sigma}(\tilde{\rho}_0)$]{} ; (12,-6.1) node\[below\][$S_{{\mathcal{X}}_\sigma}(\tilde{\rho}_1)$]{} ;

(0.3, -4.3) – (0.3,-5.7) ; (0.4,-5) node\[right\] ;

(12.3, -4.3) – (12.3,-5.7) ; (12.4,-5) node\[right\] ;

(5, -4.3) – (5,-5.7) ; (5.1,-5) node\[right\] ;

*Proof of Theorem \[theo\_GammaLimsup\]*. We have now all the tools to prove the desired result. The idea is simple: we regularize the curve, use controllability to adjust the temporal endpoints, while in the interior we can sample quite easily as we have something regular. Figure \[figure\_sketch\_proof\] represents graphically how we proceed. More specifically, we take $(\rho, {\mathbf{m}}) \in {\mathcal{M}}([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ a given pair such that ${\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) < + \infty$. Let us fix $T > 0$ and $\eta > 0$. By simplicity, we will assume that $N$ is always chosen such that $T$ is a multiple of $\tau = 1/N$.

Let $(\tilde{\rho}, \tilde{{\mathbf{m}}})$ the regularized curved given by Proposition \[prop\_regularization\]. On the other hand, let $(\hat{P}^{N, \sigma}, \hat{{\mathbf{M}}}^{N, \sigma}) \in {\mathcal{X}}_\sigma^{T/\tau +1} \times {\mathcal{Y}}_\sigma^{T/\tau}$ the curve joining $S_{{\mathcal{X}}_\sigma}(\rho_0)$ onto $S_{{\mathcal{X}}_\sigma}(\tilde{\rho}_0)$ in $T/ \tau$ time steps with a controlled cost given by Proposition \[prop\_controllability\]. Similarly with $(\bar{P}^{N, \sigma}, \bar{{\mathbf{M}}}^{N, \sigma}) \in {\mathcal{X}}_\sigma^{T/\tau +1} \times {\mathcal{Y}}_\sigma^{T/\tau}$ joining $S_{{\mathcal{X}}_\sigma}(\tilde{\rho}_1)$ onto $S_{{\mathcal{X}}_\sigma}(\rho_1)$ in $T/ \tau$ time steps with a controlled cost. Eventually, we set, for $N, \sigma$ and $k \in \{ 0,1, \ldots, N \}$ given $$P^{N, \sigma}_k := \begin{cases}
\hat{P}^{N, \sigma}_k & \text{if } k \tau \in [0,T] \\
{\displaystyle{S_{{\mathcal{X}}_\sigma} \left( \tilde{\rho}\left( T + \frac{k \tau - T}{1 - 2T}, \cdot \right)  \right)}} & \text{if } k \tau \in [T, 1-T] \\
\bar{P}^{N, \sigma}_{k - (1-T)/\tau} & \text{if } k \tau \in [1-T,T].
\end{cases}$$ We chose the momentum accordingly: if $k \in \{ 1,2, \ldots, N \}$, $${\mathbf{M}}^{N, \sigma}_k := \begin{cases}
T^{-1}\hat{{\mathbf{M}}}^{N, \sigma}_k & \text{if } (k+1/2) \tau \in [0,T] \\
(1-2T)^{-1} {\displaystyle{S_{{\mathcal{Y}}_\sigma} \left( \frac{1}{\tau} \int_{(k-1) \tau}^{k \tau}  \tilde{{\mathbf{m}}}\left( T + \frac{t - T}{1-2T}  , \cdot \right) \, {\mathrm{d}}t  \right)}} & \text{if } (k + 1/2) \tau \in [T, 1-T] \\
T^{-1} \bar{{\mathbf{M}}}^{N, \sigma}_{k  - (1-T)/\tau}  & \text{if } (k+1/2) \tau \in [1-T,T].
\end{cases}$$ Notice that we have inserted a factor $1/T$ in the momentum $\hat{{\mathbf{M}}}^{N, \sigma}$ and $\bar{{\mathbf{M}}}^{N, \sigma}$ because of the difference of temporal scaling. With our choice, $(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma})$ satisfies the discrete continuity equation. Indeed, it is clear if $k \tau \notin [T, 1-T]$ as $\hat{{\mathbf{M}}}^{N, \sigma}$, $\bar{{\mathbf{M}}}^{N, \sigma}$ are chosen for that. On the other hand, for $k$ such that $k \tau \in [T, 1-T]$, denoting by $t_k = T + (k \tau - T)/(1-2T)$, $$\begin{gathered}
P^{N, \sigma}_{k+1} - P^{N, \sigma}_k = S_{{\mathcal{X}}_\sigma} \left( \int_{t_{k-1}}^{t_k} {\partial}_t \tilde{\rho}(t, \cdot) \, {\mathrm{d}}t \right) = S_{{\mathcal{X}}_\sigma} \left(-  \int_{t_{k-1}}^{t_k} \nabla \cdot \tilde{{\mathbf{m}}}(t, \cdot) \, {\mathrm{d}}t \right) \\
= - ({\mathrm{Div}}_\sigma \circ S_{{\mathcal{Y}}_\sigma}) \left( \int_{t_{k-1}}^{t_k} \tilde{{\mathbf{m}}}((t, \cdot) \, {\mathrm{d}}t \right) =  - \tau {\mathrm{Div}}_\sigma ({\mathbf{M}}^{N, \sigma}_k),\end{gathered}$$ where the third identity comes from \[asmp\_samp\_der\] and the last one is a change of variables.

It remains to estimate the action of $(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma})$. Given the way we have built it (see Proposition \[prop\_controllability\]), and taking in account the change in the temporal scaling, the contribution of the action for $k \tau \notin [T, 1-T]$ does not exceed $T^{-1}(\hat{\omega}(\eta) + \varepsilon_{N, \sigma} \| \rho_0 \|)$ (recall that $W_2(\rho_0, \tilde{\rho}(0, \cdot)) \leqslant \eta$ and similarly for the final value). On the other hand, for the part on $[T, 1-T]$ we will use the consistency property \[asmp\_samp\_actio\]. As $\tilde{\rho}$ and $\tilde{{\mathbf{m}}}$ are smooth (they are $C^1$ and $\tilde{\rho}$ is uniformly bounded from below) we can write, up to an error $\varepsilon_\sigma$ which depends only on $\sigma$ and which tends to $0$ as $\sigma \to 0$, that for all $N, \sigma, k$, $$A_\sigma \left( \frac{P^{N, \sigma}_{k-1} + P^{N, \sigma}_{k}}{2}, {\mathbf{M}}^{N, \sigma}_k \right) \leqslant A \left( \frac{\tilde{\rho}(t_{k-1}, \cdot) + \tilde{\rho}(t_k, \cdot)}{2}, \frac{1}{\tau} \int_{t_{k-1}}^{t_k} \tilde{{\mathbf{m}}}(t, \cdot) \, {\mathrm{d}}t \right)  + \varepsilon_\sigma$$ On the other hand, as the action $A$ depends smoothly of its input when restricted to the set of measures which have a smooth density, and as $\tilde{\rho}(t, \cdot)$ and $\tilde{{\mathbf{m}}}(t, \cdot)$ depend smoothly on the temporal variable, $$\left| A \left( \frac{\tilde{\rho}(t_{k-1}, \cdot) + \tilde{\rho}(t_k, \cdot)}{2}, \frac{1}{\tau} \int_{t_{k-1}}^{t_k} \tilde{{\mathbf{m}}}(t, \cdot) \, {\mathrm{d}}t \right)  - \frac{1}{\tau(1-2T)} \int_{t_{k-1}}^{t_k} A(\tilde{\rho}(t, \cdot), \tilde{{\mathbf{m}}}(t, \cdot)) \, {\mathrm{d}}t  \right| \leqslant \frac{C}{N^2},$$ where $C$ depends on $\tilde{\rho}$ and $\tilde{{\mathbf{m}}}$, but not on $N$. Hence, summing over $k$, one can see that $$\sum_{k \ : \ k \tau \in [T, 1-2T]} \tau A_\sigma \left( \frac{P^{N, \sigma}_{k-1} + P^{N, \sigma}_{k}}{2}, {\mathbf{M}}^{N, \sigma}_k \right) \leqslant \frac{1}{1-2T} {\mathcal{A}}(\tilde{\rho}, \tilde{{\mathbf{m}}}) + \varepsilon_\sigma + \frac{C}{N}.$$ Now, given that ${\mathcal{A}}(\tilde{\rho}, \tilde{{\mathbf{m}}}) \leqslant (1 + \eta){\mathcal{A}}(\rho, {\mathbf{m}})$ by construction, putting all these information together, $$\limsup_{N \to + \infty, \sigma \to 0} \ {\mathcal{J}}^{N, \sigma}_{S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1)}(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) \leqslant \frac{2}{T} \hat{\omega}(\eta) + \frac{1+ \eta}{1- 2T} {\mathcal{A}}(\rho, {\mathbf{m}}).$$ We recall that $\eta$ and $T$ are arbitrary. If we choose first $T$ very small, and then $\eta$ small enough, we can see that we reach the desired conclusion.

On the relation with already proposed discretizations {#section_examples}
=====================================================

Now that we have proved Theorems \[theo\_GammaLiminf\] and \[theo\_GammaLimsup\], to justify the interest of our result, we will show how previous works can be embedded in this framework. This is the case for the the discretization on triangulations of surfaces used by the present author in [@Lavenant2018] with Claici, Chien and Solomon, and our general framework was mainly designed for this case. Moreover, we show that it can also be used to analyze the finite volume discretization proposed by Gladbach, Kopfer and Maas [@Gladbach2018].

Triangulations of surfaces {#subsection_finite_elements}
--------------------------

In this subsection, we want to show how our method applies to the discretization that we proposed in [@Lavenant2018]. We stronly encourgage the reader to take a look at the figures of this article, as we have numerically solved the fully discretized problem. The underlying space $X$ is now a smooth surface in ${\mathbb{R}}^3$ and we assume that we have a triangulation that approximates it in the $C^1$-sense, see [@Hildebrandt2006] and some explanations below. We use a discretization that is reminiscent of finite elements as velocity fields are defined over the triangles. The triangulation will be assumed to be *regular* but not necessarily *uniformly regular* (see below). In particular, it could allow one to try a multiscale discretization, with a refinement of the mesh only where needed. Moreover, compared to the finite volume discretization described in Subsection \[subsection\_finite\_volumes\] below, no isotropy condition is required.

Before diving into the details of the spatial discretization, let us do a remark (that can be skipped at first reading) about the temporal discretization used in [@Lavenant2018] which differs from the one of the present article. Indeed, in [@Lavenant2018] we rather focused on discretizing the dual problem, hence the resulting functional (see [@Lavenant2018 Equation (16)]) is slightly different, but can be expressed as follows by a duality argument. With the notation of the present article, if $\bar{P}_0, \bar{P}_1 \in {\mathcal{X}}_{\sigma,+}$ are given then we define $\tilde{{\mathcal{J}}}^{N, \sigma}_{\bar{P}_0, \bar{P}_1}$ on $({\mathcal{X}}_\sigma)^N \times ({\mathcal{Y}}_\sigma)^{N+1}$ by $$\tilde{{\mathcal{J}}}^{N, \sigma}_{\bar{P}_0, \bar{P}_1}( (P_k)_{1 \leqslant k \leqslant N}, ({\mathbf{M}}_k)_{0 \leqslant k \leqslant N}  ) = \frac{\tau}{2} A_\sigma( P_1, {\mathbf{M}}_0  ) + \frac{\tau}{2} A_\sigma(P_{N}, {\mathbf{M}}_N) + \tau \sum_{k=1}^{N-1} A_\sigma \left( \frac{P_{k} + P_{k+1}}{2}, {\mathbf{M}}_k \right)$$ if $(P_k)_{1 \leqslant k \leqslant N} \in ({\mathcal{X}}_{\sigma,+})^N$ and the discrete continuity equation is satisfied, that is $$\begin{cases}
\tau^{-1} (P_{k+1} - P_{k})  + {\mathrm{Div}}_\sigma ({\mathbf{M}}_{k}) = 0, & \forall k \in \{1,2, \ldots, N-1 \} \\
(\tau/2)^{-1} (P_{1} - \bar{P}_0) + {\mathrm{Div}}_\sigma( {\mathbf{M}}_0 ) = 0, & \\
(\tau/2)^{-1} (\bar{P}_1 - P_{N}) + {\mathrm{Div}}_\sigma( {\mathbf{M}}_N ) = 0 , & \\
\end{cases}$$ and $+ \infty$ otherwise. As the reader can see, the difference appears near the temporal boundary conditions (in some sense the first and last temporal step sizes are $\tau/2$ and not $\tau$), but the techniques of this article still work if one replaces ${\mathcal{J}}^{N, \sigma}$ by $\tilde{{\mathcal{J}}}^{N, \sigma}$.

Let us go through the spatial discretization. We take $X \subset {\mathbb{R}}^3$ a smooth compact $2$-dimensional submanifold of ${\mathbb{R}}^3$ without boundary. We approximate $X$ with a triangulation. Namely, we assume that for each $\sigma > 0$ small enough we have a pair $({\mathcal{T}}^\sigma, {\mathcal{V}}^\sigma)$ where ${\mathcal{V}}^\sigma \subset {\mathbb{R}}^3$ is a finite subset of ${\mathbb{R}}^3$ (the vertices), and ${\mathcal{T}}^\sigma$ is a set of triangles in ${\mathbb{R}}^3$.

Each $K \in {\mathcal{T}}^\sigma$ is a triangle in ${\mathbb{R}}^3$ whose vertices belong to ${\mathcal{V}}^\sigma$. We denote $X^\sigma = \bigcup_{K \in {\mathcal{T}}^\sigma} K$ the polyhedral surface that it generates. Notice that both $X^\sigma$ and $X$ are subsets of ${\mathbb{R}}^3$, but their intersection may be empty. We assume that $\sigma$ is the maximum diameter of elements in ${\mathcal{T}}^\sigma$. Following [@Brenner2007 Definition (4.4.13)], we assume that the family of triangulations is regular: it means that there exists a constant $c > 0$ (independent on $\sigma$) such that for every $\sigma > 0$ and every $K \in {\mathcal{T}}^\sigma$, one can fit a disk of radius $c ~ \mathrm{diam}(K)$ in $K$. It implies that all angles in the triangles are uniformly bounded from below independently on $\sigma$; that for every $\sigma > 0$ and every $K \in {\mathcal{T}}^\sigma$, the number of triangles sharing a vertex with $K$ is bounded independently on $\sigma$, and all these neighboring triangles have a diameter comparable to the one of $K$. However, contrary to the next section we do not assume that the triangulation is *uniformly regular* (that is, quasi regular in the sense of [@Brenner2007 Definition (4.4.13)]) which means that the minimal diameter of triangles in ${\mathcal{T}}^\sigma$ can be much smaller than $\sigma$.

Following [@Hildebrandt2006], we assume that the family of triangulations approximate the manifold $X$ in the $C^1$ sense. Specifically $X^\sigma$ is assumed to be a graph over $X$: it means that $X^\sigma$ always lie in the set of points in ${\mathbb{R}}^3$ for which the projection onto $X$ is smooth and uniquely defined, and that the projection onto $X$ is one to one once restricted to $X^\sigma$. In particular we can define $\Psi_\sigma : X \to X_\sigma$ as the inverse of the projection map $\Psi^{-1}_\sigma : X^\sigma \to X$, and both $\Psi_\sigma$ and $\Psi_\sigma^{-1}$ are continuous and one-to-one. The differential of $\Psi_\sigma$ at the point $x \in X$ is denoted by $D \Psi_\sigma(x) : T_x X \to {\mathbf{n}}_K^\perp \subset {\mathbb{R}}^3$, where $K$ is the triangle in which $\Psi_\sigma(x)$ lies and ${\mathbf{n}}_K$ is the normal to $K$. Notice that $D \Psi_\sigma$ is well defined up to an negligible set on $X$ (the image by $\Psi_\sigma^{-1}$ of the edges of the triangles of ${\mathcal{T}}^\sigma$). We assume that $X^\sigma$ converges in Hausdorff distance to $X$, that is $0$-th order convergence. We also assume that the convergence holds at first order in the following sense: if we denote by ${\mathbf{n}}_X : X \to {\mathbb{R}}^3$ the normal mapping and ${\mathbf{n}}_K$ the normal to the triangle $K$ then $$\lim_{\sigma \to 0} \left( \sup_{K \in {\mathcal{T}}^\sigma} \sup_{x \in K} | \mathbf{n}_K - \mathbf{n}_X( \Psi^{-1}_\sigma(x)  )  | \right) = 0.$$ We refer to [@Hildebrandt2006 Theorem 2] for alternative (but equivalent) formulations of such a convergence.

Our result will be already interesting if $X$ were flat: the reader not so familiar with differential geometry can think as $\Psi_\sigma$ being the identity mapping to translate the results in the flat case.

Let us set some additional notations first, see Figure \[figure\_finite\_elements\] for an explanation about some of them. The area (i.e. the $2$-dimensional Haussdorff measure) of a triangle $K \in {\mathcal{T}}^\sigma$ is denoted by $|K|$. The set of vertices of a triangle $K$ is ${\mathcal{V}}^\sigma_K \subset {\mathcal{V}}^\sigma$, while the set of triangles to which a vertex $v$ belongs is ${\mathcal{T}}^\sigma_v \subset {\mathcal{T}}^\sigma$. If $v \in {\mathcal{V}}^\sigma$, we denote by $\hat{\phi}_v$ the function defined on $X_\sigma$, piecewise linear on each triangle, such that $\hat{\phi}_v(v) = 1$ and $\hat{\phi}_v$ is $0$ everywhere else on ${\mathcal{V}}^\sigma$. In particular, $\hat{\phi}_v$ has a gradient constant over each triangle, and non zero only on ${\mathcal{T}}^\sigma_v$. We define $|{\mathcal{T}}^\sigma_v| = \frac{1}{3} \sum_{K \in {\mathcal{T}}^\sigma_v} |K|$ the “area” of the vertex $v$. It can also be expressed by $|{\mathcal{T}}^\sigma_v| = \int_{X^\sigma} \hat{\phi}_v$, where integration is performed w.r.t. the $2$-dimensional Haussdorff measure.

We define our finite dimensional model of dynamical optimal transport as follows.

-   ${\mathcal{X}}_\sigma = {\mathbb{R}}^{{\mathcal{V}}^\sigma}$ and ${\mathcal{X}}_{\sigma,+} = ({\mathbb{R}}_+)^{{\mathcal{V}}^\sigma}$ while $${\mathcal{Y}}_\sigma = \left\{ {\mathbf{M}}= ({\mathbf{M}}_K)_{K \in {\mathcal{T}}^\sigma} \in ({\mathbb{R}}^3)^{{\mathcal{T}}^\sigma} \ : \ \forall K \in {\mathcal{T}}^\sigma, \ {\mathbf{M}}_K \cdot \mathbf{n}_K = 0  \right\}.$$ That is, density is defined by one scalar per vertex, while momentum is one vector by triangle, lying in the vector space parallel to the triangle.

-   The divergence operator, for ${\mathbf{M}}\in {\mathcal{Y}}_\sigma$ and $v \in {\mathcal{V}}^\sigma$, is given by $$({\mathrm{Div}}_\sigma {\mathbf{M}})_v = - \frac{1}{|{\mathcal{T}}^\sigma_v|}\sum_{K \in {\mathcal{T}}^\sigma_v} |K| \left( \left. \nabla \hat{\phi}_v \right|_{K} \cdot {\mathbf{M}}_K \right).$$ This divergence is, for proper scalar products, the dual operator to the gradient restricted to the set of piecewise affine functions.

-   The action $A_\sigma$ is defined for $(P, {\mathbf{M}}) \in {\mathcal{X}}_\sigma \times {\mathcal{Y}}_\sigma$ as $$A_\sigma(P, {\mathbf{M}}) = \sum_{K \in {\mathcal{T}}^\sigma} \frac{1}{2} \frac{ |{\mathbf{M}}_K|^2  }{ (\sum_{v \in {\mathcal{V}}^\sigma_K} P_v)/3} |K|.$$ The discrete action $A_\sigma$ easily satisfies the homogeneity and monotonicity requirements of Definition \[definition\_approximation\]. In the formula above, we follow the convention of to define the quotient: such a choice makes $A_\sigma$ convex and lower semi-continuous.

Up to now, we have copied the definitions of [@Lavenant2018]. Both the values of the density and the momentum on vertices or triangles are thought as *intensive*, that is they are normalized by the volume of the vertex or the triangle they live on. Of course, all of these definitions do not require the knowledge of $X$, only the one of the triangulation. Let us now switch to reconstruction and sampling.

-   To sample, we assume that ${\mathbf{m}}\in {\mathcal{M}}(TX)$ is absolutely continuous w.r.t. volume measure, and $\rho \in {\mathcal{M}}(X)$ is arbitrary. We define $$(S_{{\mathcal{X}}_\sigma} (\rho))_v = \frac{1}{|{\mathcal{T}}^\sigma_v|} \langle  \rho, \hat{\phi}_v \circ \Psi_\sigma \rangle  \ \text{ and } \ (S_{{\mathcal{Y}}_\sigma} ({\mathbf{m}}))_K = \frac{1}{|K|} \int_{ \Psi_\sigma^{-1}(K)} D \Psi_\sigma(x) {\mathbf{m}}(x) \, {\mathrm{d}}x .$$ Notice that $S_{{\mathcal{X}}_\sigma} (\rho)$ is some sort of projection onto the set of piecewise affine functions on $X_\sigma$ of the image measure of $\rho$ by $\Psi_\sigma$.

-   There are two ways to reconstruct the density, and one way for the momentum. Roughly, for $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ we put a Dirac mass at every vertex, while for $R^A_{{\mathcal{X}}_\sigma}$ we take a density constant over each triangle, equal to the mean of the densities of the vertices of the triangle. On the other hand, for $R_{{\mathcal{Y}}_\sigma}$ we take a momentum constant over each triangle. Then, we pull everything back onto $X$ via the mapping $\Psi_\sigma$. Specifically, if $\phi \in C(X)$ and ${\mathbf{b}}\in C(X,TX)$ are test functions defined over $X$, $$\begin{gathered}
    \langle R^A_{{\mathcal{X}}_\sigma} (P), \phi \rangle = \sum_{K \in {\mathcal{T}}^\sigma}  \left( \frac{1}{3} \sum_{v \in {\mathcal{V}}^\sigma_K} P_v \right) \int_{\Psi_\sigma^{-1}(K)} \phi(x) \,{\mathrm{d}}x,
     \ \  \ \langle \ R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} (P), \phi \rangle = \sum_{v \in {\mathcal{V}}^\sigma} P_v |{\mathcal{T}}^\sigma_v| \phi( \Psi_\sigma^{-1}(v) ), \\
    \text{ and } \ \langle R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}}), {\mathbf{b}}\rangle =  \sum_{K \in {\mathcal{T}}^\sigma} {\mathbf{M}}_K \cdot \left( \int_{\Psi_\sigma^{-1}(K)}  D \Psi_\sigma(x)  {\mathbf{b}}(x) \, {\mathrm{d}}x \right).  \end{gathered}$$

It is really straightforward to check that $R^A_{{\mathcal{X}}_\sigma}$ and $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ map ${\mathcal{X}}_{\sigma,+}$ into ${\mathcal{M}}_+(X)$ while $S_{{\mathcal{X}}_\sigma}({\mathcal{M}}_+(X)) \subset {\mathcal{X}}_{\sigma,+}$.

\[proposition\_finite\_elements\_adapated\] Let $X \subset {\mathbb{R}}^3$ be a smooth compact $2$-dimensional submanifold of ${\mathbb{R}}^3$ without boundary. For each $\sigma > 0$, assume that we have $({\mathcal{T}}^\sigma, {\mathcal{V}}^\sigma)_\sigma$ a triangulation in ${\mathbb{R}}^3$, and that this family of triangulations is regular and converges as $\sigma \to 0$ to $X$ in the $C^1$ sense defined in [@Hildebrandt2006].

Then, provided that $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_\sigma$ are defined as above, we get a family of finite dimensional models of dynamical optimal transport which is adapted to $X$ (in the sense of Definition \[definition\_adapated\]), with reconstruction and sampling operators $(R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}, R_{{\mathcal{Y}}_\sigma}, S_{{\mathcal{X}}_\sigma}, S_{{\mathcal{Y}}_\sigma})_\sigma$ defined as above.

The reconstruction operators still require the knowledge of $X$, which could be disappointing. However, we can define reconstruction operators valued in the set of measures over ${\mathbb{R}}^3$ (namely, concentrated on the polyhedral surface $X^\sigma$): densities and momenta are taken for instance taken constant on each triangle. Using the notations of Theorem \[theo\_main\], in the limit $N \to + \infty, \sigma \to 0$ it is easy to see (because of the $C^1$ convergence of the surface) that the reconstructed densities and momenta on the surfaces $X_\sigma$ will converge to the same limit as ${\mathbb{R}}^A_{N,{\mathcal{X}}_\sigma}(P^{N, \sigma})$ and ${\mathbb{R}}_{N,{\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma})$, provided the latter exist. As Theorem \[theo\_main\] combined with Proposition \[proposition\_finite\_elements\_adapated\] justify that these objects have a limit, so do the reconstructions which are not pulled back on $X$.

Before doing computations, let us define some distortion coefficients for areas: if $\sigma > 0$ and $v \in {\mathcal{V}}^\sigma$ while $K \in {\mathcal{T}}^\sigma$, $$\alpha^\sigma_v := \frac{\int_{X} \hat{\phi}_v (\Psi_\sigma(x)) \, {\mathrm{d}}x  }{ |{\mathcal{T}}^\sigma_v|} - 1 \ \text{ and } \beta^\sigma_K = \frac{\int_{\Psi_\sigma^{-1}(K)} \, {\mathrm{d}}x}{|K|} - 1.$$ Thanks to the $C^1$ convergence of the surface which implies convergence of areas [@Hildebrandt2006 Theorem 2], $\sup_{v \in {\mathcal{V}}^\sigma} |\alpha^\sigma_v|$ tends to $0$ as $\sigma \to 0$, as well as $\sup_{K \in {\mathcal{T}}^\sigma} |\beta^\sigma_K|$. Analogously, we define the distortion coefficients for the metric by, for $\sigma > 0$ and $K \in {\mathcal{T}}^\sigma$, $$\theta^\sigma_K = \sup_{ x \in \Psi^{-1}_\sigma(K) }\sup_{w \in T_x X, w \neq 0} \left|  \frac{|D\Psi_\sigma(x) w|}{|w|} - 1 \right|,$$ and again still by the $C^1$ convergence of the surface $\sup_K \theta^\sigma_K$ tends to $0$ as $\sigma \to 0$.

We first check \[asmp\_samp\_der\], that is the exact commutation between derivation and sampling. If ${\mathbf{m}}\in {\mathcal{M}}(TX)$ has a smooth density then for $v \in {\mathcal{V}}^\sigma$, $$\begin{aligned}
(S_{{\mathcal{X}}_\sigma} (\nabla \cdot {\mathbf{m}}) )_v & = \frac{1}{|{\mathcal{T}}^\sigma_v|} \int_{X} \hat{\phi}_v(\Psi_\sigma(x)) (\nabla \cdot {\mathbf{m}})(x) \, {\mathrm{d}}x \\
&= - \frac{1}{|{\mathcal{T}}^\sigma_v|}  \int_X \nabla \left( \hat{\phi}_v \circ \Psi_\sigma \right) \cdot {\mathbf{m}}(x) \, {\mathrm{d}}x \\
& = - \sum_{K \in {\mathcal{T}}^\sigma_v} \frac{1}{|{\mathcal{T}}^\sigma_v|} \int_{\Psi_\sigma^{-1}(K)} \left( (D \Psi_\sigma)^\top(x) \left. \nabla \hat{\phi}_v \right|_K \right) \cdot {\mathbf{m}}(x) \,  {\mathrm{d}}x \\ 
& = - \frac{1}{|{\mathcal{T}}^\sigma_v|} \sum_{K \in {\mathcal{T}}^\sigma_v} |K| \left( \left. \nabla \hat{\phi}_v \right|_K \right) \cdot \frac{1}{|K|} \int_{\Psi_\sigma^{-1}(K)}  \left( D \Psi_\sigma (x) {\mathbf{m}}(x) \right) \, {\mathrm{d}}x 
= (( {\mathrm{Div}}_\sigma \circ S_{{\mathcal{Y}}_\sigma} ) ({\mathbf{m}}))_v. \end{aligned}$$ Next we check \[asmp\_samp\_actio\], which amounts to say that we use a legit quadrature formula for the action. Let ${\mathcal{B}}$ and ${\mathcal{B}}'$ be bounded sets of $C^1(X)$ and $C^1(X,TX)$ respectively, with functions in ${\mathcal{B}}$ uniformly bounded from below by a positive constant. Let $\rho \in {\mathcal{M}}(X)$ and ${\mathbf{m}}\in {\mathcal{M}}(TX)$ whose densities w.r.t. volume measure (still denoted by $\rho$ and ${\mathbf{m}}$) belong to ${\mathcal{B}}$ and ${\mathcal{B}}'$ respectively. Thanks to the boundedness of ${\mathcal{B}}$, and given the definition of the coefficients $\alpha^\sigma_v$, we see that $|(S_{{\mathcal{X}}_\sigma} (\rho))_v - \rho(\Psi_\sigma^{-1}(v))| \leqslant C ( \sigma + |\alpha^\sigma_v|)$ where $C$ depends only on ${\mathcal{B}}$. For sampling of the momentum, we have to use also the distortions coefficients $\theta^\sigma_K$ for the metric tensor to find $\left| |(S_{{\mathcal{Y}}_\sigma} ({\mathbf{m}}))_K| - |{\mathbf{m}}( \Psi_\sigma^{-1}(z) )| \right| \leqslant C ( \sigma + |\beta^\sigma_K| + \theta^\sigma_K  )$ provided $z$ is any point of $K$. Putting these estimates together and using that $\rho$ is bounded from below, $$\left| \frac{ |{\mathbf{M}}_K|^2  }{2 (\sum_{x \in {\mathcal{V}}^\sigma_K} P_x)/3} |K| - \frac{|{\mathbf{m}}(z)|^2}{2\rho(z)} |K| \right| \leqslant C |K| \varepsilon_\sigma$$ where $z$ is any point of $\Psi_\sigma^{-1}(K)$ and $\varepsilon_\sigma$ depends on $\sigma, \alpha^{\sigma}_v$, $\beta^\sigma_K$ and $\theta^\sigma_K$ and tends to $0$ uniformly in $K$ as $\sigma \to 0$. On the other hand, as $\rho$ and ${\mathbf{m}}$ have uniformly smooth densities, if $z \in \Psi_\sigma^{-1}(K)$ $$\left| \frac{|{\mathbf{m}}(z)|^2}{2\rho(z)} |K| - \int_{\Psi_\sigma^{-1}(K)} \frac{|{\mathbf{m}}(x)|^2}{2 \rho(x)} \,{\mathrm{d}}x  \right| \leqslant C |K| (\sigma + |\alpha^\sigma_K| ).$$ Putting these two estimates together and summing over $K$ we get \[asmp\_samp\_actio\].

Now we turn to reconstruction. The one-sided estimate for the action is clear. Actually, thanks to Lemma \[lemma\_interp\_action\], we will check \[asmp\_interp\_action\_stronger\] rather than \[asmp\_interp\_action\]: the reconstruction operator $R^A_{{\mathcal{X}}_\sigma}$ is built for that. Let us take $P \in {\mathcal{X}}_\sigma$ and ${\mathbf{M}}\in {\mathcal{Y}}_\sigma$. For every $K \in {\mathcal{T}}^\sigma$, the measure $R^A_{{\mathcal{X}}_\sigma} (P)$ has a density constant over $\Psi_\sigma^{-1}(K)$ and equal to $\frac{1}{3} \sum_{v \in {\mathcal{V}}^\sigma_K} P_v$, while the measure $R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}})$ has a density given by $(D \Psi_\sigma)^\top {\mathbf{M}}_K$. As the linear operator $D \Psi_\sigma$ is almost an isometry (up to an error $\theta^\sigma_K$), $|(D \Psi_\sigma)^\top {\mathbf{M}}_K|^2 \leqslant (1 + \theta^\sigma_K)^2 |{\mathbf{M}}_K|^2$. Hence, $$\begin{aligned}
A(R^A_{{\mathcal{X}}_\sigma} (P),R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}}) ) & = \sum_{K \in {\mathcal{T}}^\sigma} \int_{\Psi_\sigma^{-1}(K)} \frac{|(R_{{\mathcal{Y}}_\sigma} {\mathbf{M}}) (x)|^2}{ 2 (R^A_{{\mathcal{X}}_\sigma} P)(x) } \,{\mathrm{d}}x \\
& \leqslant \sum_{K \in {\mathcal{T}}^\sigma} (1+ \theta^\sigma_K)^2   \frac{ |{\mathbf{M}}_K|^2  }{2 (\sum_{v \in {\mathcal{V}}^\sigma_K} P_v)/3}  \int_{\Psi_\sigma^{-1}(K)} \,{\mathrm{d}}x \\ 
& \leqslant  \sum_{K \in {\mathcal{T}}^\sigma} (1+ \theta^\sigma_K)^2 (1+ |\beta^\sigma_K|)  \frac{ |{\mathbf{M}}_K|^2  }{2 (\sum_{x \in {\mathcal{V}}^\sigma_K} P_x)/3} |K| \\
&\leqslant \left( \sup_{K \in {\mathcal{T}}^\sigma}  (1+ \theta^\sigma_K)^2 (1+ |\beta^\sigma_K|) \right) A_\sigma(P, {\mathbf{M}}),\end{aligned}$$ and the factor in front of $ A_\sigma(P, {\mathbf{M}})$ in the last line tends to $1$ as $\sigma \to 0$. We emphasize that we have used the convention in the computation above. On the other hand, to check \[asmp\_interp\_der\] commutation between reconstruction and derivation, let us take ${\mathcal{B}}$ a bounded set of $C^2(X)$ and $\phi \in {\mathcal{B}}$ and compute: $$\begin{gathered}
\langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ {\mathrm{Div}}_\sigma)(P), \phi \rangle 
= - \sum_{v \in {\mathcal{V}}^\sigma} \left[ \sum_{K \in {\mathcal{T}}^\sigma_v} |K| \left( \left. \nabla \hat{\phi}_v \right|_{K} \cdot {\mathbf{M}}_K \right) \right] \phi( \Psi_\sigma^{-1}(v) ) \\
= - \sum_{K \in {\mathcal{T}}^\sigma} |K| {\mathbf{M}}_K \cdot \left( \sum_{v \in {\mathcal{V}}^\sigma_K} \phi(\Psi^{-1}_\sigma(v)) \left. \nabla \hat{\phi}_v \right|_K  \right).\end{gathered}$$ Now we notice that $\sum_{v \in {\mathcal{V}}^\sigma_K} \phi(\Psi^{-1}_\sigma(v)) \left. \nabla \hat{\phi}_v \right|_K$ is the gradient restricted to $K$ of the function $\hat{\phi} : X^\sigma \to {\mathbb{R}}$ which is continuous, piecewise affine on each triangle, and whose value on the vertex $v \in {\mathcal{V}}^\sigma$ coincide with $\phi( \Psi_\sigma^{-1}(v) )$. Thus, $$\langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ {\mathrm{Div}}_\sigma)(P), \phi \rangle  = - \sum_{K \in {\mathcal{T}}^\sigma} |K| {\mathbf{M}}_K \cdot \left. \nabla \hat{\phi} \right|_K.$$ As $\hat{\phi}$ coincides with $\phi \circ \Psi^{-1}_\sigma : X_\sigma \to {\mathbb{R}}$ on ${\mathcal{V}}^\sigma$, thanks to the regularity of the mesh we know [@Brenner2007 Corollary (4.4.24)] that at least on a triangle $K$, $$\sup_{y \in K} \left| \left. \nabla \hat{\phi} \right|_K - \nabla \left( \phi \circ \Psi_\sigma^{-1} \right)(y) \right| \leqslant C \sigma,$$ where the constant $C$ depends on the second derivatives of $\phi \circ \Psi_\sigma^{-1}$ over $K$. As $K$ is flat, the second derivatives of $\Psi_\sigma^{-1}$ restricted to $K$ depend only on the derivatives of the metric tensor of the underlying space $X$. The latter is assumed to be smooth, while the derivatives of $\phi$ are bounded by a constant that depends only on ${\mathcal{B}}$. Hence the constant $C$ depends only on $X$ and ${\mathcal{B}}$, but does not depend on $\sigma$. On the other hand, by definition of $R_{{\mathcal{Y}}_\sigma}$, there holds $$\langle R_{{\mathcal{Y}}_\sigma} ({\mathbf{M}}), \nabla \phi \rangle =  \sum_{K \in {\mathcal{T}}^\sigma} {\mathbf{M}}_K \cdot \int_{\Psi_\sigma^{-1}(K)}  D \Psi_\sigma(x) \nabla \phi(x) \,{\mathrm{d}}x.$$ Putting these two information together, $$\begin{gathered}
\left| \langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ {\mathrm{Div}}_\sigma)(P), \phi \rangle + \langle (R_{{\mathcal{Y}}_\sigma} {\mathbf{M}}), \nabla \phi \rangle  \right| \\
\leqslant \left| \sum_{K \in {\mathcal{T}}^\sigma} {\mathbf{M}}_K \cdot \left( - \int_K \nabla \left( \phi \circ \Psi_\sigma^{-1} \right)(y) \,{\mathrm{d}}y + \int_{\Psi_\sigma^{-1}(K)}  D \Psi_\sigma(x) \nabla \phi(x) \,{\mathrm{d}}x  \right) \right| 
+ C \sigma \left( \sum_{K \in {\mathcal{T}}^\sigma}  |K| |{\mathbf{M}}_K| \right)  . \end{gathered}$$ By $C^1$ convergence of the triangulations, we can deduce that each term in front of ${\mathbf{M}}_K$ is bounded by an error $|K| C \varepsilon_\sigma$, where $\varepsilon_\sigma$ depends on $\Psi_\sigma$ but tends to $0$ as $\sigma$ tends to $0$, while $C$ depends on ${\mathcal{B}}$. This is enough to yield \[asmp\_interp\_der\].

The first point \[asmp\_interp\_sampling\] is very straightforward and left to the reader. On the other hand, for the second one \[asmp\_interps\], that is to check that $R^A_{{\mathcal{X}}_\sigma}$ and $R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ are asymptotically equivalent, let us take ${\mathcal{B}}$ a bounded set of $C^1(X)$ and $\phi \in {\mathcal{B}}$, then $$\langle ( R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} - R^A_{{\mathcal{X}}_\sigma} )(P), \phi  \rangle = \sum_{v \in {\mathcal{V}}^\sigma} P_v \left( |{\mathcal{T}}^\sigma_v| \phi(\Psi_\sigma^{-1}(v)) - \sum_{K \in {\mathcal{T}}^\sigma_v} \frac{1}{3} \int_{\Psi_\sigma^{-1}(K)} \phi  \right).$$ As $\phi \in {\mathcal{B}}$, for each $v \in {\mathcal{V}}^\sigma$ the discrepancy between $\fint_{\Psi_\sigma^{-1}(K)} \phi$ and $\phi(\Psi_\sigma^{-1}(v))$ is smaller than $C \sigma$. Moreover, using the fact that the distortions coefficients $\alpha^\sigma_v$ are tending to $0$ when $\sigma \to 0$, we get the result.

Eventually we need to prove controllability \[asmp\_controllability\]. The idea is to start with controllability at the continuous level: if $x,y \in X$, there exists $\gamma : [0,1] \to X$ a constant-speed joining them. Up to a small modification of the curve, we can assume that $\Psi_\sigma \circ \gamma$ intersects edges of the triangulation ${\mathcal{T}}^\sigma$ in a finite number of points while $|\dot{\gamma}(t)| \leqslant 2 d_g(x,y)$ for all $t \in [0,1]$. Let $\hat{\rho} \in {\mathcal{M}}([0,1] \times X)$ and $\hat{{\mathbf{m}}}_1, \hat{{\mathbf{m}}}_2 \in {\mathcal{M}}(TX)$ the continuous objects built according to . We set $\hat{P} = S_{{\mathcal{X}}_\sigma}(\hat{\rho})$ and $\hat{{\mathbf{M}}}_i = S_{{\mathcal{Y}}_\sigma}(\hat{{\mathbf{m}}}_i)$ for $i \in \{1,2\}$. The reader can check, as $\Psi_\sigma \circ \gamma$ intersects finitely many edges, that $S_{{\mathcal{Y}}_\sigma}(\hat{{\mathbf{m}}}_i)$ is well defined and, by \[asmp\_samp\_der\], that ${\mathrm{Div}}_\sigma( \hat{{\mathbf{M}}}_1 ) = \hat{P} - S_{{\mathcal{X}}_\sigma}(\delta_x)$, and ${\mathrm{Div}}_\sigma( \hat{{\mathbf{M}}}_2 ) = \hat{P} - S_{{\mathcal{X}}_\sigma}(\delta_y)$. We need to estimate the action (and we cannot use \[asmp\_samp\_actio\] because of a lack of smoothness), let us do it for $A_\sigma(\hat{P}, \hat{{\mathbf{M}}}_1)$. Let $K$ be one of the triangle that $\Psi_\sigma \circ \gamma$ crosses and $\hat{\rho}_K = \hat{\rho}( \Psi_\sigma^{-1}(K))$. As the mesh is regular, all the triangles sharing one vertex with $K$ have an area comparable to the one of $K$, hence $|{\mathcal{T}}^\sigma_v|$ is also comparable to $|K|$ for $v$ any vertex of $K$. Thus we can estimate, up to a constant $C$ independent on $\sigma$ and changing from one equation to an other, $$\frac{1}{3} \sum_{v \in {\mathcal{V}}^\sigma_K} \hat{P}_v \geqslant \frac{1}{C |K|} \sum_{v \in {\mathcal{V}}^\sigma_K} \langle \hat{\rho}, \hat{\phi}_v \circ \Psi_\sigma  \rangle 
\geqslant  \frac{1}{C |K|} \sum_{v \in {\mathcal{V}}^\sigma_K} \langle \hat{\rho}, {\mathbbm{1}}_{K} \circ \Psi_\sigma \rangle \geqslant \frac{\hat{\rho}_K}{C|K|}$$ where we have first used that $\sum_{v \in {\mathcal{V}}^\sigma_K} \hat{\phi}_v \geqslant {\mathbbm{1}}_K$ and then the definition of $\hat{\rho}_K$ as well as the low distortion of areas induced by $\Psi_\sigma$. On the other hand, recall that $\hat{{\mathbf{m}}}_1$ has a density bounded by $2 d_g(x,y)$ w.r.t. $\hat{\rho}$. Hence, still on this triangle $K$ and starting from the definition of $\hat{{\mathbf{m}}}_1$ and $S_{{\mathcal{Y}}_\sigma}$ one can easily estimate $|\hat{{\mathbf{M}}}_1| \leqslant C \hat{\rho}_K d_g(x,y) / |K|$, where $C$ is independent on $\sigma$ and depends on the operator norm of $D \Psi_\sigma$. Putting these estimates together, $$A_\sigma(\hat{P}, \hat{{\mathbf{M}}}_1) \leqslant \sum_{K \in {\mathcal{V}}^\sigma} \frac{1}{2 ( \hat{\rho}_{K} / C|K| )  } \left( C \frac{\hat{\rho}_{K} d_g(x,y)}{|K|} \right)^2 |K| \leqslant C d_g(x,y)^2 \sum_{K \in {\mathcal{V}}^\sigma}  \hat{\rho}_{K} = C d_g(x,y)^2,$$ where the last identity comes from the definition of the $\hat{\rho}_K$ and the fact that $\hat{\rho}$ has unit mass.

Finite volumes {#subsection_finite_volumes}
--------------

This part focuses on the discretization proposed by Gladbach, Kopfer and Maas in [@Gladbach2018]. We strongly advise the reader to read the article and reference therein to get a better understanding of this discretization and the results known about it. We take $X = \Omega \subset {\mathbb{R}}^d$ a convex compact subset of ${\mathbb{R}}^d$ with smooth boundary. In particular ${\mathcal{M}}(TX) = {\mathcal{M}}(\Omega)^d$ is the space of $d$-dimensional vectorial measures.

For each $\sigma > 0$, we assume that we have a mesh $({\mathcal{T}}^\sigma, (x_K)_{K \in {\mathcal{T}}^\sigma})$ that is: ${\mathcal{T}}^\sigma$ is a partition of $\Omega$ in convex sets (called cells) with non empty interior; and for each $K \in {\mathcal{T}}^\sigma$, we have a point $x_K$ located in the interior of the cell $K$. For the sake of clarity we will work in the simplest framework proposed by the aforementioned article: we will only consider symmetric means. As identified by [@Gladbach2018], we will need to make an *isotropy* condition for the convergence to hold, see below.

We refer to Figure \[figure\_finite\_volumes\] to visualize some notations. The volume of each cell $K$ is denoted by $|K|$. For each pair $(K,L) \in ({\mathcal{T}}^\sigma)^2$, we assume that $\bar{K} \cap \bar{L} =: (K|L)$ is a flat polytope with $d-1$-dimensional Haussdorff measure $|(K|L)|$. Although cells included in the interior of $\Omega$ are polytopes, we allow for the boundary of a cell to coincide with ${\partial}\Omega$, hence be curved. If $|(K|L)| > 0$, we write $K \sim L$ and denote by ${\mathcal{E}}^\sigma = \{ (K, L) \ : \ K \sim L \} \subset {\mathcal{T}}^\sigma \times {\mathcal{T}}^\sigma$ the set of “edges”. For a cell $K$, we denote by ${\mathcal{T}}^\sigma_K \subset {\mathcal{T}}^\sigma$ the set of neighboring cells, i.e. ${\mathcal{T}}^\sigma_K = \{ L \in {\mathcal{T}}^\sigma \ : \ K \sim L \}$. We assume that if $K \sim L$ then $x_K - x_L$ is orthogonal to $(K|L)$. We denote by $d_{KL}$ the distance between $x_K$ and $x_L$ and ${\mathbf{n}}_{KL} = (x_L -x_K)/d_{KL}$ the normal to $(K|L)$.

We assume that the mesh is *uniformly regular* in the following sense. There exists a constant $c > 0$ independent on $\sigma$ such that, for any $\sigma > 0$ the following holds: the diameter of each cell is bounded by $c^{-1} \sigma$; additionaly for any $K \in {\mathcal{T}}^\sigma$ the ball of center $x_K$ and radius $c \sigma$ is contained in $K$; and if $(K,L) \in {\mathcal{E}}^\sigma$ then $|(K|L)| \geqslant c \sigma^{d-1}$. The terminology slightly differs from [@Gladbach2018]: we have added the word “uniformly” to make the distinction with the regularity assumption of the discretization presented in the previous section.

Eventually, following [@Gladbach2018 Definition 1.3] we make an isotropy assumption. Specifically there exists $(\varepsilon_\sigma)_\sigma$, tending to $0$ as $\sigma \to 0$ such that for every $\sigma > 0$, every cell $K \in {\mathcal{T}}^\sigma$ and every vector $v \in {\mathbb{R}}^d$, there holds $$\label{equation_isotropy}
\frac{1}{2 |K|} \sum_{L \in {\mathcal{T}}^\sigma_K} |(K|L)| d_{KL} (v \cdot {\mathbf{n}}_{KL})^2 \leqslant (1+ \varepsilon_\sigma) |v|^2.$$ We refer to [@Gladbach2018] to get a better understanding of this condition, which the authors have identified as being necessary for convergence to hold. As they indicate, if the mesh satisfies the “center of mass” condition, that is if for $(K,L) \in {\mathcal{E}}^\sigma$ the mean between $x_K$ and $x_L$ is also the center of mass of $(K|L)$, then the isotropy condition holds.

We give also ourselves a *symmetric* mean $\theta : {\mathbb{R}}_+ \times {\mathbb{R}}_+  \to {\mathbb{R}}_+$ admissible in the sense of [@Gladbach2018 Definition 2.3], that is: a smooth, positively $1$-homogeneous, symmetric, jointly concave function such that $\theta(1,1) = 1$. It implies in particular that $\theta(a,b) \leqslant (a+b)/2$.

Then, for $\sigma$ given, that is for a given mesh, we define the approximation of dynamical optimal transport as follows.

-   ${\mathcal{X}}_\sigma = {\mathbb{R}}^{{\mathcal{T}}^\sigma}$ while ${\mathcal{X}}_{\sigma,+} = ({\mathbb{R}}_+)^{{\mathcal{T}}^\sigma}$ and ${\mathcal{Y}}_\sigma \subset {\mathbb{R}}^{{\mathcal{E}}^\sigma}$ is the set of $(M_{KL})_{(K,L) \in {\mathcal{E}}^\sigma}$ such that $M_{KL} = - M_{LK}$ for any $(K,L) \in {\mathcal{E}}^\sigma$. As elements of ${\mathcal{Y}}_\sigma$ are scalar fields rather than vectorial ones, a generic element in ${\mathcal{Y}}_\sigma$ will be denoted by $M$ and not ${\mathbf{M}}$.

-   The divergence operator is defined by, if $M \in {\mathcal{Y}}_\sigma$ and $K \in {\mathcal{T}}^\sigma$, $$({\mathrm{Div}}_\sigma M)_K := \sum_{L \in {\mathcal{T}}^\sigma_K} \frac{|(K|L)|}{|K|} M_{KL}.$$

-   The discrete action is defined on ${\mathcal{X}}_\sigma \times {\mathcal{Y}}_\sigma$ as $$A_\sigma(P, M) := \sum_{(K,L) \in {\mathcal{E}}^\sigma} \frac{M_{KL}^2}{ 4 \theta(P_K, P_L)} |(K|L)| d_{KL} 
    $$ We have a supplemental $1/2$ factor as each edge appears twice in ${\mathcal{E}}^\sigma$. The discrete action $A_\sigma$ easily satisfies the homogeneity and monotonicity requirements of Definition \[definition\_approximation\]. In the formula above we follow the convention of to define the quotient: such a choice makes $A_\sigma$ convex and lower semi-continuous.

In other words, the density is defined on the cells and the momentum on the edges, i.e. at interfaces between cells. Both the values of the density and the momentum on cells or edges are thought as *intensive*, that is they are normalized by the volume or area of the cell or edge they live on. Up to now, we have just copied the definitions of [@Gladbach2018].

As we will have to consider assumption \[asmp\_interp\_action\], we need to specify a representation of ${\mathcal{Y}}'_\sigma$ the dual of ${\mathcal{Y}}_\sigma$. To that extent we consider ${\mathcal{Y}}_\sigma$ as a (finite dimensional) Hilbert space with scalar product define by, for $B,M \in {\mathcal{Y}}_\sigma$, $$\langle M, B \rangle = \sum_{(K,L) \in {\mathcal{E}}^\sigma} \frac{M_{KL} B_{KL}}{2} d_{KL} |(K|L)|.$$ With that choice which allows us to identify ${\mathcal{Y}}'_\sigma$ with ${\mathcal{Y}}_\sigma$, the Legendre transform of the action is simply given by, for $P \in {\mathcal{X}}_{\sigma, +}$ and $B \in {\mathcal{Y}}_\sigma$, $$A^\star_\sigma (P, B) = \sum_{(K,L) \in {\mathcal{E}}^\sigma} \frac{\theta(P_K, P_L)}{ 4 } (B_{KL})^2 d_{KL} |(K|L)|.$$

Now we need to define reconstruction and sampling operators.

-   To sample, we assume that ${\mathbf{m}}$ has a continuous density w.r.t. the $d$-dimensional Lebesgue measure, still denoted by ${\mathbf{m}}$, while $\rho \in {\mathcal{M}}(X)$ is arbitrary. We define $$(S_{{\mathcal{X}}_\sigma} (\rho))_K = \frac{\rho(K)}{|K|} = \fint_K {\mathrm{d}}\rho \ \text{ and } \ (S_{{\mathcal{Y}}_\sigma} ({\mathbf{m}}))_{KL} = \frac{1}{|(K |L)|} \int_{(K|L)} {\mathbf{m}}\cdot \mathbf{n}_{KL} \, {\mathrm{d}}\mathcal{H}^{d-1} = \fint_{(K|L)} {\mathbf{m}}\cdot \mathbf{n}_{KL},$$ where $\mathcal{H}^{d-1}$ is the $(d-1)$-dimensional Hausdorff measure.

-   On the other hand, reconstruction operators are defined as follows. Let $P \in {\mathcal{X}}_\sigma$ and $M \in {\mathcal{Y}}_\sigma$. We define $$\begin{gathered}
    R^A_{{\mathcal{X}}_\sigma} (P) = \sum_{K \in {\mathcal{T}}^\sigma} P_K {\mathbbm{1}}_K,
     \ \ R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} (P) = \sum_{K \in {\mathcal{T}}^\sigma} P_K |K| \delta_{x_K}, \\
    \text{ and } \ R_{{\mathcal{Y}}_\sigma} (M) = \frac{1}{2} \sum_{(K,L) \in {\mathcal{E}}^\sigma}  M_{KL} d_{KL} {\mathbf{n}}_{KL} {\mathbbm{1}}_{(K|L)}.  \end{gathered}$$ More specifically, the measure $R_{{\mathcal{Y}}_\sigma} (M)$ is a $(d-1)$-dimensional measure concentrated on the interfaces between cells: if ${\mathbf{b}}\in C(\Omega, {\mathbb{R}}^d)$ then $$\langle R_{{\mathcal{Y}}_\sigma} (M), {\mathbf{b}}\rangle = \frac{1}{2} \sum_{(K,L) \in {\mathcal{E}}^\sigma} M_{KL} d_{KL} \int_{(K|L)} {\mathbf{b}}(x) \cdot {\mathbf{n}}_{KL}  \, \mathcal{H}^{d-1}({\mathrm{d}}x).$$ Again notice the factor $1/2$ as each edge is counted twice.

It is really straightforward to check that $R^A_{{\mathcal{X}}_\sigma}, R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}$ map ${\mathcal{X}}_{\sigma,+}$ into ${\mathcal{M}}_+(X)$ while $S_{{\mathcal{X}}_\sigma}({\mathcal{M}}_+(X)) \subset {\mathcal{X}}_{\sigma,+}$. It is also very easy to check that in this case, given the scalar product that we put on ${\mathcal{Y}}_\sigma$, that there holds $R_{{\mathcal{Y}}_\sigma}^\top = S_{{\mathcal{Y}}_\sigma}$ (provided one identifies measures with their density w.r.t. Lebesgue measure).

Let $\Omega \subset {\mathbb{R}}^d$ compact convex with a smooth boundary. For each $\sigma > 0$, assume that we have $({\mathcal{T}}^\sigma, (x_K)_{K \in {\mathcal{T}}^\sigma})$ an admissible uniformly regular mesh of diameter $\sigma$ and which satisfies the isotropy condition. Let $\theta$ be an admissible symmetric mean in the sense of of [@Gladbach2018].

Then, provided that $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_\sigma$ are defined as above, we get a family of finite dimensional models of dynamical optimal transport which is adapted to $\Omega$ (in the sense of Definition \[definition\_adapated\]), with reconstruction and sampling operators $(R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}, R_{{\mathcal{Y}}_\sigma}, S_{{\mathcal{X}}_\sigma}, S_{{\mathcal{Y}}_\sigma})_\sigma$ defined as above.

Let us show that what we have defined satisfies all the requirements of Definition \[definition\_adapated\]. We will not check the different requirements in the order they are stated, but rather highlight the way the objects are chosen.

We first check commutation between derivation and sampling \[asmp\_samp\_der\]. This is just the application of the divergence theorem. Indeed, if ${\mathbf{m}}\in {\mathcal{M}}(\Omega)^d $ has a smooth density and satisfies ${\mathbf{m}}\cdot {\mathbf{n}}_\Omega = 0$ on ${\partial}\Omega$ then for $\sigma > 0$ and $K \in {\mathcal{T}}^\sigma$, $$(S_{{\mathcal{X}}\sigma} (\nabla \cdot {\mathbf{m}}))_K = \fint_K (\nabla \cdot {\mathbf{m}}) = \sum_{L \sim K} \frac{1}{|K|} \int_{(K|L)} {\mathbf{m}}\cdot \mathbf{n}_{KL} = (({\mathrm{Div}}_\sigma \circ S_{{\mathcal{Y}}_\sigma})({\mathbf{m}}))_K.$$ The asymptotic exactitude for the action $A_\sigma$ when evaluated on samples of smooth densities and momenta \[asmp\_samp\_actio\] amounts to check that the quadrature formula is correct. Let ${\mathcal{B}}$ and ${\mathcal{B}}'$ be bounded sets of $C^1(\Omega)$ and $C^1(\Omega,{\mathbb{R}}^d)$ respectively, with functions in ${\mathcal{B}}$ uniformly bounded from below by a positive constant. Let $\rho \in {\mathcal{M}}(\Omega)$ and ${\mathbf{m}}\in {\mathcal{M}}(\Omega)^d$ whose density w.r.t. Lebesgue measure (still denoted by $\rho$ and ${\mathbf{m}}$) belong to ${\mathcal{B}}$ and ${\mathcal{B}}'$ respectively. In particular, up to an error controlled by $C \sigma$, we can replace $(S_{{\mathcal{Y}}_\sigma}({\mathbf{m}}))_{KL}$ by ${\mathbf{m}}(x_K) \cdot {\mathbf{n}}_{KL}$ and $(S_{{\mathcal{X}}_\sigma}(\rho))_{K}$ by $\rho(x_K)$ hence $$\begin{aligned}
A_\sigma(S_{{\mathcal{X}}_\sigma}(P), S_{{\mathcal{Y}}_\sigma}({\mathbf{m}})) & = \sum_{(K,L) \in {\mathcal{E}}^\sigma} \frac{ |(S_{{\mathcal{Y}}_\sigma}({\mathbf{m}}))_{KL}|^2  }{ 4 \theta( (S_{{\mathcal{X}}_\sigma}(\rho))_{K} , (S_{{\mathcal{X}}_\sigma}(\rho))_{L} )} d_{KL} |(K|L)| \\
& \leqslant \sum_{K \in {\mathcal{T}}^\sigma} \frac{1}{4 \rho(x_K)} \sum_{L \in {\mathcal{T}}^\sigma_K} |{\mathbf{m}}(x_k) \cdot \mathbf{n}_{KL}|^2 d_{KL} |(K|L)| + C \sigma \\
& \leqslant (1+ \varepsilon_\sigma) \sum_{K \in {\mathcal{T}}^\sigma} |K| \frac{|{\mathbf{m}}(x_k)|^2}{2 \rho(x_k)}  + C \sigma 
= (1+\varepsilon_\sigma)(A(\rho, {\mathbf{m}}) + C \sigma) + C \sigma ,  \end{aligned}$$ where the third inequality comes from and the last one is just a quadrature formula. As $A(\rho, {\mathbf{m}}) $ is bounded by a constant which depends only on ${\mathcal{B}}$ and ${\mathcal{B}}'$, we get the result. Although we have used it was not necessary and we could have relied only on [@Gladbach2018 Lemma 5.4] at the price of slightly more involved estimates.

Now we turn to reconstruction. For the one-sided estimation on the action, we use \[asmp\_interp\_action\] rather than \[asmp\_interp\_action\_stronger\], as we have not been able to build reconstruction operators such that the latter holds. We rely crucially on the isotropy assumption for this step. Let ${\mathcal{B}}$ a bounded set of $C^1(\Omega, {\mathbb{R}}^d)$ and ${\mathbf{b}}\in {\mathcal{B}}$, as well as $P \in {\mathcal{X}}_\sigma$. Using the property $\theta(a,b) \leqslant (a+b)/2$, we get easily that $$\begin{gathered}
A_\sigma^\star( P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}) ) 
= \sum_{(K,L) \in {\mathcal{E}}^\sigma} \frac{\theta(P_K, P_L)}{ 4 } ((R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}))_{KL})^2 d_{KL} |(K|L)| \\
\leqslant  \frac{1}{4} \sum_{K \in {\mathcal{T}}^\sigma} P_K \sum_{L \in {\mathcal{T}}^\sigma_K} ((R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}))_{KL})^2 d_{KL} |(K|L)|.  \end{gathered}$$ Given the regularity of ${\mathbf{b}}$, it is clear that $(R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}))_{KL}$ can be replaced by ${\mathbf{b}}(x_K) \cdot {\mathbf{n}}_{KL}$ up to an error controlled by $\sigma$. Hence we can write, using the isotropy condition , that $$\begin{gathered}
A_\sigma^\star( P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}) ) 
\leqslant \frac{1}{4} \sum_{K \in {\mathcal{T}}^\sigma} P_K \sum_{L \in {\mathcal{T}}^\sigma_K} ( {\mathbf{b}}(x_K) \cdot {\mathbf{n}}_{KL} )^2 d_{KL} |(K|L)| + C \sigma \| R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} (P) \|  \\ 
\leqslant  (1+ \varepsilon_\sigma) \frac{1}{2} \sum_{K \in {\mathcal{T}}^\sigma} P_K |K| \left| {\mathbf{b}}(x_K) \right|^2 + C \sigma \| R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} (P) \| . \end{gathered}$$ On the other hand, by smoothness of ${\mathbf{b}}$, it is clear that $|K| |{\mathbf{b}}(x_K)|^2 \leqslant \int_{K} |{\mathbf{b}}|^2 + C |K| \sigma$, and that ${\mathbf{b}}$ is uniformly bounded by a constant that depends only on ${\mathcal{B}}$. Hence we conclude with the definition of $R^A_{{\mathcal{X}}_\sigma}$ that $$A_\sigma^\star( P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}) )  \leqslant \frac{1}{2} \sum_{K \in {\mathcal{T}}^\sigma} P_K \int_K |{\mathbf{b}}|^2 + \varepsilon_\sigma \| R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P) \| = A^\star( R^A_{{\mathcal{X}}_\sigma}(P), {\mathbf{b}}) + \varepsilon_\sigma \| R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}(P) \|.$$ To check the second estimate in \[asmp\_interp\_action\], let us take ${\mathbf{b}}\in C(\Omega,{\mathbb{R}}^d)$, by homogeneity of the formula w.r.t. ${\mathbf{b}}$ we assume that $|{\mathbf{b}}(x)| \leqslant 1$ for all $x \in \Omega$. In particular, $((R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}))_{KL})^2 \leqslant 1$ for all $(K,L) \in {\mathcal{E}}^\sigma$, hence $$A_\sigma^\star( P, R_{{\mathcal{Y}}_\sigma}^\top ({\mathbf{b}}) ) 
\leqslant \frac{1}{2} \sum_{K \in {\mathcal{T}}^\sigma} P_K \sum_{L \in {\mathcal{T}}^\sigma_K} d_{KL} |(K|L)|.$$ By uniform regularity of the mesh, up to a constant independent on $\sigma$ the quantity $\sum_{L \in {\mathcal{T}}_K^\sigma} d_{KL} |(K|L)|$ can be controlled by $|K|$, hence we get the second estimate in \[asmp\_interp\_action\].

On the other hand, to check \[asmp\_interp\_der\] the commutation between reconstruction and derivation, let us take ${\mathcal{B}}$ a bounded set of $C^2(\Omega)$ and $\phi \in {\mathcal{B}}$, as well as $M \in {\mathcal{Y}}_\sigma$. Then $$\begin{aligned}
\langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} & \circ {\mathrm{Div}}_\sigma)  (M), \phi \rangle + \langle R_{{\mathcal{Y}}_\sigma} (M), \nabla \phi \rangle \\ 
& = \sum_{K \in {\mathcal{T}}^\sigma} \sum_{L \in {\mathcal{T}}^\sigma_K} \phi(x_K) |(K|L)| M_{KL} + \frac{1}{2} \sum_{(K,L) \in {\mathcal{E}}^\sigma} \left(  M_{KL} d_{KL} \fint_{(K|L)} \nabla \phi \cdot {\mathbf{n}}_{KL}  \right) \\
& = \sum_{(K,L) \in {\mathcal{E}}^\sigma} |(K|L)| \frac{M_{KL}}{2} \left(  ( \phi(x_K) - \phi(x_L) ) + d_{KL}  \fint_{(K|L)} \nabla \phi \cdot {\mathbf{n}}_{KL} \right).\end{aligned}$$ By smoothness of $\phi$, one can say that $$\left| ( \phi(x_K) - \phi(x_L) ) + d_{KL}  \fint_{(K|L)} \nabla \phi \cdot {\mathbf{n}}_{KL} \right| \leqslant C \sigma^2,$$ where $C$ depends only on ${\mathcal{B}}$, hence we deduce that the total error can be bounded by $$\left| \langle (R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}  \circ {\mathrm{Div}}_\sigma)  M, \phi \rangle + \langle R_{{\mathcal{Y}}_\sigma} (M), \nabla \phi \rangle \right| \leqslant C \sigma \sum_{(K,L) \in {\mathcal{E}}^\sigma} |M_{KL}| |(K | L)| d_{KL} \leqslant C \sigma \| R_{{\mathcal{Y}}_\sigma} (M) \|,$$ which is the desired result.

The first and second points \[asmp\_interp\_sampling\] and \[asmp\_interps\] are very straightforward to check, this is left to the reader.

Lastly, we need to check \[asmp\_controllability\], that is controllability. As explained in Remark \[rmk\_controllability\], we mimic the construction in the continuous case. Let $x, y \in \Omega$ and $K_x, K_y \in {\mathcal{T}}^\sigma$ the cells to which they respectively belong. By regularity of the mesh, there exists a path $K_x = K_1, K_2,  \ldots , K_Q = K_y$ of length $Q$ such that $K_{i-1} \sim K_{i}$ for every $i$ and $C^{-1} (Q-1) \sigma \leqslant \sum_i |x_{K_{i-1}} - x_{K_i}| \leqslant C (|x-y| + \sigma)$. The density $\hat{P} \in {\mathcal{X}}_\sigma$ is chosen as $$\hat{P}_{K_i} = \frac{1}{Q |K_i|},$$ and $\hat{P}_K = 0$ if $K$ is not in the chain $\{ K_1, \ldots, K_Q \}$. On the other hand, the momentum is chosen to be $0$ on every edge not connecting a two consecutive cells in $\{ K_1, \ldots, K_R \}$ and for every $i \in \{ 2, \ldots, Q \}$ we set $$(\hat{M}_1)_{K_{i-1} K_i} = - \frac{Q+1-i}{Q |(K_{i-1}|K_i)|} \  \text{ and } \ (\hat{M}_2)_{K_{i-1} K_i} = \frac{i-1}{Q |(K_{i-1}|K_i)|},$$ and of course the opposite on the edge $(K_i,K_{i-1})$. It is just a matter of a straightforward computation to check that ${\mathrm{Div}}_\sigma ({\mathbf{M}}_1) = \hat{P} - S_{{\mathcal{X}}_\sigma} (\delta_x)$ and ${\mathrm{Div}}_\sigma ({\mathbf{M}}_2) = \hat{P} - S_{{\mathcal{X}}_\sigma} (\delta_y)$. On the other hand, as the mesh is $\sigma$-regular, we know that $|K| \sim C \sigma^d$ and $|(K|L)| \sim C \sigma^{d-1}$ for any $(K,L) \in {\mathcal{E}}^\sigma$. We can thus estimate $$\begin{gathered}
A_\sigma(\hat{P}, \hat{M}_1) = \sum_{i=2}^{Q} \frac{\left( \frac{Q+1-i}{Q |(K_{i-1}|K_i)|} \right)^2}{2 \theta \left( \frac{1}{Q|K_{i-1}|}, \frac{1}{Q|K_i|} \right) } |(K_{i-1}|K_i)| d_{K_{i-1} K_i} \\ \leqslant \frac{C}{Q} \left( \sum_{i=2}^Q (Q+1-i)^2 \right) \frac{\sigma}{\sigma^{d-1} \sigma^{-d}} \leqslant C ((Q-1) \sigma)^2.\end{gathered}$$ Given the estimate $(Q-1) \sigma \leqslant C( |x-y| + \sigma )$, we indeed get the announced result. The estimate for $A_\sigma(\hat{P}, \hat{M}_2)$ is entirely analogous.

Extensions {#section_extension}
==========

In this section we investigate two natural extensions: instead of fixing the final value of the density, we add a cost depending on it in the objective functional; alternatively we can also add a running cost depending on the density. Although the first extension is rather easy, the second one is more involved and we only explain what the additional difficulties are without fully addressing them.

One alternative that we not explore is the use of other Lagrangians, that is to replace $|{\mathbf{m}}|^2 / (2 \rho)$ by $\rho L(\cdot, {\mathbf{m}}/\rho)$ where the Lagrangian $L : TX \to [0, + \infty]$ is convex in its second variable but not necessarily quadratic. In such a situation, applying the heat flow does not necessarily decrease the action hence the techniques of proof for Proposition \[prop\_regularization\] would break down. Nevertheless, if $X$ is the torus, one could hope to use convolution instead of applying the heat flow and rely on the convexity of the Lagrangian. We have preferred to explore complicated geometries (that is Riemannian manifolds like in [@Lavenant2018]) restricted to quadratic Lagrangians rather than simple geometries (taking $X$ to be a torus) with more general Lagrangians: we find the results more interesting in the first case rather than in the second one, but this is mainly a matter of personal taste.

Final penalization of the density
---------------------------------

Given the *minimizing movement scheme* for gradient flows [@Jordan1998] in the Wasserstein, it is natural to take $G : {\mathcal{M}}_+(X) \to [0, + \infty]$ and to consider the following problem, whose unknowns are $\rho \in {\mathcal{M}}([0,1] \times X)$, ${\mathbf{m}}\in {\mathcal{M}}([0,1] \times TX)$ and $\rho_1 \in {\mathcal{M}}_+(X)$ (but $\rho_0 \in {\mathcal{M}}_+(X)$ is given): $$\label{equation_problem_JKO}
\min_{\rho, {\mathbf{m}}, \rho_1} {\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) + G(\rho_1).$$ In other words, the final value of the density is no longer fixed but penalized with the help of $G$. Of course to approximate the problem the idea is to use $G_\sigma$ an approximation of $G$ which is defined over ${\mathcal{X}}_\sigma$. The natural requirements for it to be a good approximation of $G$ are the followings.

\[definition\_adapted\_final\] Let $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_\sigma$ be a finite dimensional model of dynamical optimal transport with reconstruction and sampling operators $(R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}, R_{{\mathcal{Y}}_\sigma}, S_{{\mathcal{X}}_\sigma}, S_{{\mathcal{Y}}_\sigma})_\sigma$.

Assume that $G : {\mathcal{M}}_+(X) \to [0, + \infty]$ is a given functional over the space of measures. For every $\sigma > 0$ we give ourselves $G_\sigma : {\mathcal{X}}_\sigma \to [0, + \infty]$. The family of functions $(G_\sigma)_\sigma$ is an adapted discretization of $G$ if:

1.  \[asmp\_G\_GammaLiminf\] There exists $(\varepsilon_\sigma)_\sigma$ which tends to $0$ when $\sigma \to 0$ such that, for any $\sigma > 0$ and any $P \in {\mathcal{X}}_\sigma$, $$G( R^A_{{\mathcal{X}}_\sigma} (P)  ) \leqslant (1+ \varepsilon_\sigma) G_\sigma( P).$$

2.  \[asmp\_G\_GammaLimsup\] For any $\rho \in {\mathcal{M}}_+(X)$, there holds $$\limsup_{\sigma \to 0} G_\sigma( S_{{\mathcal{X}}_\sigma} (\rho) ) \leqslant G( \rho ).$$

This definition was chosen for the following result to hold.

Let $({\mathcal{X}}_\sigma, {\mathcal{Y}}_\sigma, A_\sigma, {\mathrm{Div}}_\sigma)_\sigma$ a family of finite dimensional models of dynamical optimal transport which is adapted to $X$ in the sense of Definition \[definition\_adapated\], with reconstruction and sampling operators $(R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma}, R^A_{{\mathcal{X}}_\sigma}, R_{{\mathcal{Y}}_\sigma}, S_{{\mathcal{X}}_\sigma}, S_{{\mathcal{Y}}_\sigma})_\sigma$. Assume also that $(G_\sigma)_\sigma$ is an adapted discretization of a lower semi-continuous functional $G : {\mathcal{M}}_+(X) \to [0, + \infty]$ in the sense of Definition \[definition\_adapted\_final\].

Let $\rho_0 \in {\mathcal{M}}_+(X)$ be given.

For any $N$ and any $\sigma > 0$, let $(P^{N, \sigma}, {\mathbf{M}}^{N, \sigma})$ a minimizer of the functional $$(P, {\mathbf{M}}) \in ({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^N \mapsto {\mathcal{J}}^{N, \sigma}_{ S_{{\mathcal{X}}_\sigma}(\rho_0), P_N } (P, {\mathbf{M}}) + G_\sigma(P_N).$$ Then, in the limit $N \to + \infty$ and $\sigma \to 0$, up to the extraction of a subsequence, ${\mathcal{R}}^{{\mathcal{CE}}}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$ and ${\mathcal{R}}^{A}_{N, {\mathcal{X}}_\sigma}(P^{N, \sigma})$ will converge weakly to the same limit $\rho \in {\mathcal{M}}([0,1] \times X)$, the measure $R^{A}_{{\mathcal{X}}_\sigma}(P_N^{N, \sigma})$ will converge weakly to $\rho_1 \in {\mathcal{M}}_+(X)$ and ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma})$ will converge weakly to ${\mathbf{m}}\in {\mathcal{M}}([0,1] \times TX)$, in addition the triplet $(\rho, {\mathbf{m}}, \rho_1)$ is a solution of the problem and $$\lim_{N \to + \infty, \sigma \to 0} \left( {\mathcal{J}}^{N, \sigma}_{ S_{{\mathcal{X}}_\sigma}(\rho_0), P^{N, \sigma}_N } (P^{N, \sigma}, {\mathbf{M}}^{N, \sigma}) + G_\sigma(P^{N, \sigma}_N) \right) =  {\mathcal{J}}_{\rho_0, \rho_1} ( \rho, {\mathbf{m}}) + G(\rho_1).$$

As the proof closely follows the proofs of Section \[section\_proof\], we will only sketch it.

We can copy almost line by line the proof of Theorem \[theo\_GammaLiminf\]. Indeed, the uniform control on the masses of the measures ${\mathcal{R}}^{{\mathcal{CE}}}_{N, \sigma}(P^{N, \sigma})$, ${\mathcal{R}}^{A}_{N, \sigma}(P^{N, \sigma})$ and ${\mathcal{R}}_{N, {\mathcal{Y}}_\sigma}({\mathbf{M}}^{N, \sigma})$ is obtained only through the assumption that $(R^{{\mathcal{CE}}}_{{\mathcal{X}}_\sigma} \circ S_{{\mathcal{X}}_\sigma} )(\rho_0)$ is bounded uniformly in $\sigma$, which is still the case. In particular, we deduce that $R^A_{{\mathcal{X}}_\sigma}( P^{N, \sigma}_N )$ is bounded in ${\mathcal{M}}_+(X)$, hence converges to a limit $\rho_1$. Passing to the limit in the continuity equation is done exactly in the same way. Then to pass to the limit in the objective functional, we use the estimate \[asmp\_G\_GammaLiminf\] as well as the lower semi-continuity of $G$ to see that $$\liminf_{N \to + \infty, \sigma \to 0}  G_\sigma(P^{N, \sigma}_N) \geqslant G(\rho).$$ Combining with , we get the analogous of Theorem \[theo\_GammaLiminf\].

To go the other way around, that is to prove the analogous of Theorem \[theo\_GammaLimsup\], if $(\rho, {\mathbf{m}}, \rho_1)$ is a solution of the problem , we can use exactly the same discrete competitors as those built in the proof of Theorem \[theo\_GammaLimsup\]. The only new term to handle is the boundary term, but this is easily done with \[asmp\_G\_GammaLimsup\] if we take $P^{N, \sigma}_N = S_{{\mathcal{X}}_\sigma}(\rho_1)$.

Comments on how to add a running cost depending on the density
--------------------------------------------------------------

Starting from considerations in mean field games [@Benamou2016], it is also natural to take $F : {\mathcal{M}}_+(X) \to [0, + \infty]$ convex lower semi-continuous and to consider problems of the form $$\label{equation_problem_MFG}
\min_{\rho, {\mathbf{m}}} {\mathcal{J}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) + \int_0^1 F(\rho_t) \,{\mathrm{d}}t,$$ where $\rho_0$ and $\rho_1$ are fixed or penalized. For the notation $\rho_t$, we refer to Remark \[rmk\_rho\_curve\]. For instance, one can take $F$ to be a quadratic penalization of the density, or (minus) the Boltzmann entropy. The adaptation of the strategy of this article to this case is more involved, and would imply to dive deeper in the details of the discretization. For this reason, and not to overburden the article we will not give precise statements but only give hints about what the additional difficulties are. Let us point out that one can see numerical results in [@Lavenant2018 Section 5.4] when $F$ is a quadratic penalization.

Of course, the approach would be to choose $F_\sigma$ an adapted discretization of $F$ in the sense of Definition \[definition\_adapted\_final\], and then to approximate Problem with the minimization of $$(P, {\mathbf{M}}) \in ({\mathcal{X}}_\sigma)^{N+1} \times ({\mathcal{Y}}_\sigma)^N \mapsto {\mathcal{J}}^{N, \sigma}_{ S_{{\mathcal{X}}_\sigma}(\rho_0), S_{{\mathcal{X}}_\sigma}(\rho_1) }(P,{\mathbf{M}}) + \tau \sum_{k=1}^N F_\sigma \left( \frac{P_{k-1} + P_k}{2} \right).$$

The analogue of Theorem \[theo\_GammaLiminf\] is very easy to prove once we assume that $F$ is lower semi-continuous. However, in order to prove an analogue of Theorem \[theo\_GammaLimsup\], one faces two additional difficulties.

-   Now, in the analogue of the controllability property (Proposition \[prop\_controllability\]), one also needs to control (using the notations of \[asmp\_controllability\]) $F_\sigma(\hat{P})$ with $F(\rho_0)$ and $F(\rho_1)$. If we use the construction of the discretization on triangulations of surfaces (see Section \[subsection\_finite\_elements\]), such control would rely on *geodesic convexity* of the functional $F$ (i.e. convexity along geodesics in the Wasserstein space, which is linked to the Ricci curvature of the underlying space), as well as the convexity of $F$, and a finer analysis of the reconstruction and sampling operators.

-   Moreover, one would want that the regularization procedure (Proposition \[prop\_regularization\]) does not increase too much the running cost $\int F(\rho_t) \, {\mathrm{d}}t$. Given the proof of Proposition \[prop\_regularization\], this is easy if $F$ is either continuous for the topology of weak convergence, or decreases when composed with the heat flow. For instance, $F$ being a potential energy or an internal energy (provided the function giving the density of energy in terms of the density is convex) would easily fit in the framework.

Acknowledgments {#acknowledgments .unnumbered}
===============

The present article was mostly written while the author was a PhD student at *Université Paris-Sud*, and completed while he was a postdoctoral fellow of the *Pacific Institute for the Mathematical Sciences*. The author acknowledges the support of ANR MAGA (ANR-16-CE40-0014). He wants to thank Quentin Mérigot for fruitful discussions.

Wasserstein distance {#section_wasserstein}
====================

We briefly recall the definition of the quadratic Wasserstein distance, see [@AGS; @Villani2008; @SantambrogioOTAM] to learn more. We recall that $d_g$ is the geodesic distance on the Riemannian manifold $(X,g)$. The quadratic Wasserstein distance is denoted by $W_2$: if $\mu, \nu \in {\mathcal{M}}_+(X)$ are measures on $X$ sharing the same total mass then $$W_2^2(\mu, \nu) := \min_{\pi} \iint_{X \times X} d_g(x,y)^2 \,\pi( {\mathrm{d}}x, {\mathrm{d}}y)$$ where $\pi$ is a measure on the product space $X \times X$ whose first marginal is $\mu$ and second marginal is $\nu$. There always exists an optimal $\pi$, and it is called an optimal transport plan between $\mu$ and $\nu$.

For $a > 0$, let ${\mathcal{P}}_a(X) \subset {\mathcal{M}}_+(X)$ the set of positive measures over $X$ whose total mass is $a$. Then $W_2$ metrizes the weak convergence of measures on ${\mathcal{P}}_a(X)$.

See for instance [@Villani2008 Theorem 6.9], and we emphasize that the compactness of $X$ is necessary for this statement to hold.

Although, as highlighted in Theorem \[theorem\_J\_W2\], the Wasserstein distance is closely related to the problem we are studying, we have only used the definition above as a convenient tool to quantify convergence in measure, see Propositions \[prop\_regularization\] and \[prop\_controllability\].

Regularizing curves of measures {#section_regularization_curves}
===============================

Regularizing curves valued in the space of probability measures has been attacked in a variety of setting. In Proposition \[prop\_regularization\], as we want both a regularization of the density and the momentum, and we are in a Riemannian manifold possibly with a smooth and convex boundary, we do not think that a statement for this precise situation already exists. We refer the reader to [@AGS Chapter 8] or [@SantambrogioOTAM Chapter 5] for regularization arguments on flat spaces, and to [@Otto2005 Proposition 4.3] or [@Erbar2010] for regularization on Riemannian manifolds without boundaries.

The idea is to use the heat flow to regularize because it has the effect of decreasing the action. However, we apply the heat flow only to the density, and we choose the momentum *a posteriori*, once the density is regularized. To keep it consistent with the point of view chosen in this article, we prefer to use the Bakry-Émery estimate to prove that the action is decreases along the heat flow.

For $\rho \in {\mathcal{M}}_+([0,1] \times X)$ and $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ we define $${\mathcal{A}}_{\rho_0, \rho_1}(\rho) = \min_{{\mathbf{m}}} \left\{ {\mathcal{A}}(\rho, {\mathbf{m}}) \ : \ {\mathbf{m}}\in {\mathcal{M}}([0,1] \times TX) \text{ and } (\rho, {\mathbf{m}}) \in {\mathcal{CE}}(\rho_0, \rho_1) \right\}.$$

In fact $\rho_0$ and $\rho_1$ are superfluous because if ${\mathbf{m}}\in {\mathcal{M}}([0,1] \times X)$ is such that $(\rho, {\mathbf{m}})$ satisfies the continuity equation (let’s say tested again functions compactly supported in $(0,1) \times X$) and ${\mathcal{A}}_{\rho_0, \rho_1}(\rho, {\mathbf{m}}) < + \infty$, then one can recover $\rho_0$ and $\rho_1$ from $\rho$, see Remark \[rmk\_rho\_curve\]. However, as we want to avoid to deal with disintegration, we prefer to stick to this redundant formulation.

The quantity ${\mathcal{A}}_{\rho_0, \rho_1}(\rho)$ admits a “dual” formulation which will be helpful in the subsequent estimates.

\[lemma\_dual\_action\] Let $\rho \in {\mathcal{M}}_+([0,1] \times X)$ and $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$. Then there holds $${\mathcal{A}}_{\rho_0, \rho_1}(\rho) = \sup_{\phi} \left\{ \langle \rho_1, \phi(1, \cdot) \rangle - \langle \rho_0, \phi(0, \cdot) \rangle - \left\llangle \rho, {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2  \right\rrangle  \ : \ \phi \in C^1([0,1] \times X)  \right\}.$$

We introduce the Hilbert space $H = L^2_\rho([0,1] \times X, TX)$: it is the space of measurable vector fields ${\mathbf{v}}: [0,1] \times X \to TX$ endowed with scalar product $$\llangle {\mathbf{v}}, {\mathbf{w}}\rrangle_H = \iint_{[0,1] \times X} g_x({\mathbf{v}}(t,x), {\mathbf{w}}(t,x)) \, \rho({\mathrm{d}}t, {\mathrm{d}}x)$$ and with a norm $\| {\mathbf{v}}\|_H^2 = \llangle {\mathbf{v}}, {\mathbf{v}}\rrangle_H$, where vector fields which coincide $\rho$-a.e. are considered equal. We know that ${\mathcal{A}}(\rho, {\mathbf{m}})$ is finite if and only if ${\mathbf{m}}$ has a density w.r.t. $\rho$, hence we can switch back again to take ${\mathbf{w}}$ the density of ${\mathbf{m}}$ w.r.t. $\rho$ as an unknown. Given the definition of ${\mathcal{CE}}(\rho_0, \rho_1)$ we can write $${\mathcal{A}}_{\rho_0, \rho_1}(\rho) = \min_{{\mathbf{w}}\in H} \left\{ \frac{1}{2} \| {\mathbf{w}}\|^2_H \ : \  \forall \phi \in C^1([0,1] \times X), \ \llangle {\mathbf{w}}, \nabla \phi \rrangle_H = \langle \rho_1, \phi(1, \cdot) \rangle - \langle \rho_0, \phi(0, \cdot) \rangle -\llangle \rho, {\partial}_t \phi \rrangle \right\}.$$ Let us call ${\mathbf{v}}$ the optimal vector field in the problem above. We see that ${\mathbf{v}}$ is the projection, in $H$, of the vector $0$ onto the affine space $$Y = \left\{ {\mathbf{w}}\in H \ : \ \forall \phi \in C^1([0,1] \times X), \ \llangle {\mathbf{w}}, \nabla \phi \rrangle_H = \langle \rho_1, \phi(1, \cdot) \rangle - \langle \rho_0, \phi(0, \cdot) \rangle -\llangle \rho, {\partial}_t \phi \rrangle  \right\}.$$ In particular ${\mathbf{v}}$ is orthogonal to the linear space $Y^\perp$, hence it belongs to the closure in $H$ of $\{ \nabla \phi \ : \ \phi \in C^1([0,1] \times X) \}$.

To get the conclusion, we simply write that $\frac{1}{2} \| {\mathbf{v}}\|^2_H = \sup_{{\mathbf{w}}\in H} \llangle {\mathbf{v}}, {\mathbf{w}}\rrangle_H - \frac{1}{2} \| {\mathbf{w}}\|^2_H$. Hence, for all $\phi \in C^1([0,1] \times X)$, taking ${\mathbf{w}}= \nabla \phi$, $${\mathcal{A}}_{\rho_0, \rho_1}(\rho) = \frac{1}{2} \| {\mathbf{v}}\|^2_H \geqslant \llangle {\mathbf{v}}, \nabla \phi \rrangle_H - \frac{1}{2} \| \nabla \phi \|^2_H = \langle \rho_1, \phi(1, \cdot) \rangle - \langle \rho_0, \phi(0, \cdot) \rangle - \llangle \rho, {\partial}_t \phi \rrangle - \frac{1}{2} \llangle \rho, |\nabla \phi|^2  \rrangle.$$ Moreover, taking a sequence $(\phi_n)_{n \in {\mathbb{N}}}$ in $C^1([0,1] \times X)$ such that $\nabla \phi_n$ converges to ${\mathbf{v}}$ in $H$ as $n \to + \infty$, we get the result.

Now we can introduce the main tool to regularize on manifold, that is the heat flow. Let $\Phi : [0, \infty) \times C(X) \to C(X)$ denote the heat flow on $X$ with Neumann boundary conditions. Specifically, if $a \in C(X)$, let us call $u : (0, + \infty) \times X \to X$ the unique solution to the heat equation with initial condition $a$, that is the Cauchy problem $$\begin{cases}
{\partial}_s u = \Delta u & \text{ on } (0, + \infty) \times \ring{X}, \\
u \cdot {\mathbf{n}}_{{\partial}X} = 0 & \text{ on } (0, + \infty) \times {\partial}X, \\
{\displaystyle{ \lim_{s \to 0} u(s, \cdot)  }} = a & \text{ in } C(X), 
\end{cases}$$ where ${\mathbf{n}}_X : {\partial}X \to TX$ denotes the outward normal and $\Delta$ is the Laplace-Beltrami operator on $X$. Then for $s \geqslant 0$ we define $\Phi_s a \in C(X)$ as $u(s, \cdot)$.

\[theorem\_bakry\_emery\] Let $(X,g)$ be a smooth compact Riemannian manifold with a (possibly empty) smooth and convex boundary. Then there exists $C \in {\mathbb{R}}$ such that for all $\phi \in C^1(X)$ and any $s \geqslant 0$, $$\left|\nabla \left( \Phi_s \phi \right) \right|^2 \leqslant e^{C s} \Phi_s \left( |\nabla \phi|^2 \right).$$

As $X$ is smooth and compact, its Ricci curvature is bounded from below by $-C/2 \in {\mathbb{R}}$. To conclude, we refer the reader to [@Wang2011 Theorem 1.1] for the case where $X$ may have a convex boundary.

We have all the tools at our disposal to prove Proposition \[prop\_regularization\]. In the proof, if $a$ is a function defined on ${\mathbb{R}}$ and $b$ a function defined on $X$, then $a \otimes b$ is the function defined on ${\mathbb{R}}\times X$ by $(a \otimes b)(t,x) = a(t) b(x)$. By a slight abuse of notations, if $\rho \in {\mathcal{M}}([0,1] \times X)$ and $\phi \in C({\mathbb{R}}\times X)$ then $\llangle \rho, \phi \rrangle$ denotes the duality product between $\rho$ and the restriction of $\phi$ to $[0,1] \times X$.

Let $(\rho, {\mathbf{m}}) \in {\mathcal{M}}_+([0,1] \times X) \times {\mathcal{M}}([0,1] \times TX)$ be given, and $\rho_0, \rho_1 \in {\mathcal{M}}_+(X)$ such that $(\rho, {\mathbf{m}}) \in {\mathcal{CE}}(\rho_0, \rho_1)$ and ${\mathcal{A}}(\rho, {\mathbf{m}}) < + \infty$.

To regularize, we extend the curve $\rho$ by setting it equal to $\rho_0$ for $t \leqslant 0$ and $\rho_1$ for $t \geqslant 1$, and then we take a convolution in time and apply the heat flow in space. To be rigorous, we do that by duality. Let $(\chi_n)_{n \in {\mathbb{N}}}$ be a smooth even approximation of unity on ${\mathbb{R}}$ such that $\chi_n$ is compactly supported in $[-1/n, 1/n]$ for any $n$. For any $\phi \in C([0,1] \times X)$, we define $\Theta_n(\phi) \in C({\mathbb{R}}\times X)$ by, for any $t \in {\mathbb{R}}$, $$\Theta_n(\phi)(t, \cdot) = \int_0^1 \chi_n(t-s) \Phi_{1/n}( \phi(s, \cdot) ) \, {\mathrm{d}}s.$$ Then, we define $\tilde{\rho}_n \in {\mathcal{M}}([0,1] \times X)$ by $$\llangle \tilde{\rho}_n, \phi \rrangle := \left\langle \rho_0, \int_{- \infty}^0 \Theta_n(\phi)(t, \cdot) \, {\mathrm{d}}t \right\rangle + \llangle \rho, \Theta_n(\phi) \rrangle + \left\langle \rho_1, \int_{1}^{+ \infty} \Theta_n(\phi)(t, \cdot) \, {\mathrm{d}}t \right\rangle.$$ By parabolic regularity and regularity of $\chi_n$, for a fixed $n$ we know that $\Theta_n(\phi)$ belongs to $C^\infty([0,1] \times X)$, with derivatives controlled by $\| \phi \|_{L^1([0,1] \times X)}$. It allows easily, by duality, to say that $\tilde{\rho}_n$ has a density w.r.t. ${\mathrm{d}}x \otimes {\mathrm{d}}t$, and that this density (still denoted by $\tilde{\rho}_n$) belongs to $C^\infty([0,1] \times X)$. Moreover, as $\Theta_n(\phi)(t, \cdot)$ is uniformly bounded from below by a strictly positive constant depending only on $\| \phi(t, \cdot) \|_{L^1(X)}$ if $\phi$ is nonnegative, we know that $\tilde{\rho}_n$ is uniformly bounded from below by a strictly positive constant.

Now we have to estimate the action, we will do it with the help of Lemma \[lemma\_dual\_action\]. First we identify the boundary conditions: let us define $\tilde{\rho}_{n,0}$ and $\tilde{\rho}_{n,1}$ by, for any $\phi \in C(X)$, $$\begin{gathered}
\langle \tilde{\rho}_{n,0}, \phi \rangle := \frac{1}{2} \langle \rho_0, \Phi_{1/n} \phi \rangle + \left\llangle \rho, \chi_n  \otimes ( \Phi_{1/n} \phi) \right\rrangle, \\ 
\text{and } \langle \tilde{\rho}_{n,1}, \phi \rangle := \frac{1}{2} \langle \rho_1, \Phi_{1/n} \phi \rangle + \left\llangle \rho, (\chi_n(1 - \cdot)) \otimes ( \Phi_{1/n} \phi) \right\rrangle. \end{gathered}$$ They could be guessed by the way we have heuristically chosen $\tilde{\rho}_n$, or simply retro-engineered from the computations below.

Next, we need to understand how $\Theta_n$ commutes with derivatives. For the temporal derivative, a integration by parts leads to, for $\phi \in C^1([0,1] \times X)$, $${\partial}_t [\Theta_n(\phi)](t, \cdot) = \int_0^1 \chi_n'(t-s) \Phi_{1/n}( \phi(s, \cdot) ) \, {\mathrm{d}}s 
= - \chi_n(t-1) \Phi_{1/n} \phi(1, \cdot) + \chi_n(t) \Phi_{1/n} \phi(0, \cdot) + \Theta_n({\partial}_t \phi)(t, \cdot).$$ On the other hand, using Theorem \[theorem\_bakry\_emery\] as well as the convexity of the squared norm, $$|\nabla \Theta_n(\phi)|^2 \leqslant e^{C/n} \Theta_n(|\nabla \phi|^2).$$

Now let us take $\phi \in C^1([0,1] \times X)$ as a test function to evaluate ${\mathcal{A}}_{\tilde{\rho}_{n,0}, \tilde{\rho}_{n,1}}(\tilde{\rho}_n)$. We compute: $$\begin{aligned}
\langle \tilde{\rho}_{n,1} & , \phi(1, \cdot) \rangle - \langle \tilde{\rho}_{n,0}, \phi(0, \cdot) \rangle - \left\llangle \tilde{\rho}_n, {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2  \right\rrangle \\
& =  \frac{1}{2} \langle \rho_1, \Phi_{1/n} \phi(1, \cdot) \rangle + \left\llangle \rho, (\chi_n(1 - \cdot) ) \otimes ( \Phi_{1/n} \phi(1, \cdot)) \right\rrangle \\
& - \frac{1}{2} \langle \rho_0, \Phi_{1/n} \phi(0, \cdot) \rangle - \left\llangle \rho, \chi_n  \otimes ( \Phi_{1/n} \phi(0, \cdot)) \right\rrangle \\
& -  \left\langle \rho_0, \int_{- \infty}^0 \Theta_n \left( {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2 \right)  (t, \cdot) \, {\mathrm{d}}t \right\rangle  - \left\langle \rho_1, \int_{0}^{+ \infty} \Theta_n \left( {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2 \right)(t, \cdot) \, {\mathrm{d}}t \right\rangle \\
& - \left\llangle \rho, \Theta_n \left( {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2 \right)  \right\rrangle.\end{aligned}$$ We do some cleaning: in the fourth line, we discard the terms in $|\nabla \phi|^2$ as they have a sign. In the last line, we use the estimate coming from Theorem \[theorem\_bakry\_emery\]. Also, in the fourth and the last lines, we use the explicit expression that we have for $\Theta_n({\partial}_t \phi)$. Specifically, for the fourth line: $$\begin{gathered}
\left\langle  \rho_0  , \int_{- \infty}^0 \Theta_n \left( {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2 \right)  (t, \cdot) \, {\mathrm{d}}t \right\rangle  \\
\leqslant \left\langle \rho_0, \int_{- \infty}^0 \left( {\partial}_t [\Theta_n(\phi)](t, \cdot) + \chi_n(t-1) \Phi_{1/n} \phi(1, \cdot) - \chi_n(t) \Phi_{1/n} \phi(0, \cdot)  \right) \, {\mathrm{d}}t \right\rangle \\ 
= \langle \rho_0, \Theta_n(\phi)(0, \cdot) \rangle - \frac{1}{2} \langle \rho_0, \Phi_{1/n} \phi(0, \cdot) \rangle.\end{gathered}$$ Of course, a very similar computation holds for the term involving $\rho_1$. On the other hand, we can also estimate in the last line $$\begin{gathered}
\left\llangle \rho, \Theta_n \left( {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2 \right)  \right\rrangle \\
 \geqslant \left\llangle \rho, \left( {\partial}_t [\Theta_n(\phi)] + (\chi_n(\cdot-1)) \otimes \Phi_{1/n} \phi(1, \cdot) - \chi_n \otimes \Phi_{1/n} \phi(0, \cdot)  \right)  \right\rrangle + e^{-C/n} \left\llangle \rho, \frac{1}{2} \left| \nabla \Theta_n(\phi) \right|^2 \right\rrangle.\end{gathered}$$ Putting all these computations together, we find many cancellations and end up with $$\begin{gathered}
\langle \tilde{\rho}_{n,1}  , \phi(1, \cdot) \rangle - \langle \tilde{\rho}_{n,0}, \phi(0, \cdot) \rangle - \left\llangle \tilde{\rho}_n, {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2  \right\rrangle \\
\leqslant \langle \rho_1, \Theta_n(\phi)(1, \cdot) \rangle - \langle \rho_0, \Theta_n(\phi)(0, \cdot) \rangle - \llangle \rho,  {\partial}_t [\Theta_n(\phi)] \rrangle - e^{-C/n} \left\llangle \rho, \frac{1}{2} \left| \nabla \Theta_n(\phi) \right|^2 \right\rrangle.\end{gathered}$$ The right hand side can be bounded using Lemma \[lemma\_dual\_action\] for ${\mathcal{A}}_{\rho_0, \rho_1}(\rho)$ with test function $e^{-C/n} \Theta_n(\phi)$. Hence we get $$\langle \tilde{\rho}_{n,1}  , \phi(1, \cdot) \rangle - \langle \tilde{\rho}_{n,0}, \phi(0, \cdot) \rangle - \left\llangle \tilde{\rho}_n, {\partial}_t \phi + \frac{1}{2} |\nabla \phi|^2  \right\rrangle \leqslant e^{C/n} {\mathcal{A}}_{\rho_0, \rho_1}(\rho).$$ Taking the supremum over $\phi$ and given Lemma \[lemma\_dual\_action\], we see that $${\mathcal{A}}_{\tilde{\rho}_{n,0}, \tilde{\rho}_{n,1}}(\tilde{\rho}_n) \leqslant e^{C/n} {\mathcal{A}}_{\rho_0, \rho_1}(\rho) \leqslant e^{C/n} {\mathcal{A}}(\rho, {\mathbf{m}}).$$

As we have seen above, the measure $\tilde{\rho}_n$ admits a smooth density bounded from below, and thanks to the equation above its action can be controlled. They are two things left to do: check that the boundary values of $\tilde{\rho}_n$ are almost $(\rho_0, \rho_1)$, and choose the optimal momentum.

For the boundary values, we need to check that $\tilde{\rho}_{n,0}$ converges weakly to $\rho_0$ as $n \to + \infty$ (the computation is almost identical for $\tilde{\rho}_{n,1}$). Let $\phi \in C^1(X)$ a test function. Let $\Gamma_n : [0, + \infty) \to {\mathbb{R}}$ a primitive of $-\chi_n$ whose value is $1/2$ for $t=0$ (hence $\Gamma_n(t) = 0$ for $t \geqslant 1/n$). As a test function in the continuity equation, we use $\Gamma_n \otimes (\Phi_{1/n} \phi)$ to write $$\frac{1}{2} \langle \rho_0, \Phi_{1/n} \phi \rangle = \llangle \rho, \chi_n  \otimes ( \Phi_{1/n} \phi ) \rrangle - \llangle {\mathbf{m}}, \Gamma_n \otimes (\nabla (\Phi_{1/n} \phi)) \rrangle.$$ Notice that $\Gamma_n \leqslant {\mathbbm{1}}_{[0,1/n]}$ while $\nabla (\Phi_{1/n} \phi)$ is bounded uniformly, for instance thanks to Theorem \[theorem\_bakry\_emery\]. As ${\mathbf{m}}$ is a finite measure, by dominated convergence $$\lim_{n \to + \infty} \left( \llangle {\mathbf{m}}, \Gamma_n \otimes (\nabla (\Phi_{1/n} \phi) \rrangle \right) = 0.$$ With the help of this information, we can conclude that $$\lim_{n \to + \infty} \langle \tilde{\rho}_{n,0}, \phi \rangle = \langle \rho_0, \phi \rangle.$$ As the function $\phi$ is smooth but arbitrary, we conclude to the weak convergence of $\tilde{\rho}_{n,0}$ to $\rho_0$, hence to convergence in Wasserstein distance as the mass of $\tilde{\rho}_{n,0}$ coincides with the one of $\rho_0$.

Eventually, we need to choose the momentum: we will choose the optimal one given $\tilde{\rho}_n$. Recall that $\tilde{\rho}_n \in C^\infty([0,1] \times X)$ denotes, by abuse of notations, the density of $\tilde{\rho}_n$ w.r.t. ${\mathrm{d}}t \otimes {\mathrm{d}}x$. For each $t \in [0,1]$ we choose $\psi_n(t, \cdot)$ the unique solution with $0$-mean of the elliptic equation $$\label{equation_continuity_equation_smooth}
\nabla \cdot (  \tilde{\rho}_n(t,\cdot) \nabla \psi_n(t, \cdot) ) = -{\partial}_t \tilde{\rho}_n(t, \cdot)$$ with Neumann boundary conditions on ${\partial}X$. As $\tilde{\rho}_n(t, \cdot)$ is smooth and bounded from below, this equation always admits a solution. Such a solution is smooth by elliptic regularity, up to the boundary as ${\partial}X$ is smooth. By setting $\tilde{{\mathbf{m}}}_n = \nabla \psi_n \tilde{\rho}_n$, we get that $(\tilde{\rho}_n, \tilde{{\mathbf{m}}}_n)$ satisfies the continuity equation ($\psi_n$ was built for that in ) and $${\mathcal{A}}(\tilde{\rho}_n, \tilde{{\mathbf{m}}}_n) = {\mathcal{A}}_{\tilde{\rho}_{n,0}, \tilde{\rho}_{n,1}}(\tilde{\rho}_n).$$ (For the latter, one can use for instance Lemma \[lemma\_dual\_action\] with $\psi_n$ as a test function). As a consequence, ${\mathcal{A}}(\tilde{\rho}_n, \tilde{{\mathbf{m}}}_n) \leqslant e^{C/n} {\mathcal{A}}(\rho, {\mathbf{m}})$, and we set $(\tilde{\rho}, \tilde{{\mathbf{m}}}) = (\tilde{\rho}_n, \tilde{{\mathbf{m}}}_n)$ for $n$ large enough to reach the conclusion of the proposition (large enough so that the temporal boundary conditions almost match and the increase in the action is small enough).

\[rmk\_m\_not\_optimal\] Using the notations of the proof of Proposition \[prop\_regularization\] above, it is easy to check that $\tilde{\rho}_n$ converges to $\rho$ weakly in ${\mathcal{M}}_+([0,1] \times X)$. However, and in particular if $${\mathcal{A}}_{\rho_0, \rho_1}(\rho) < {\mathcal{A}}(\rho, {\mathbf{m}})$$ then $\tilde{{\mathbf{m}}}_n$ will *not* converge to ${\mathbf{m}}$ as $n \to + \infty$. In fact, ${\mathbf{m}}$ was almost superfluous in our analysis as we have not tried to regularize it and rather choose the momentum $\tilde{{\mathbf{m}}}_n$ thanks to .
