# Using Sinkhorn in the JKO scheme adds linear diffusion

Aymeric Baradat, Anastasiia Hraivoronska, Filippo Santambrogio

► To cite this version:

Aymeric Baradat, Anastasiia Hraivoronska, Filippo Santambrogio. Using Sinkhorn in the JKO scheme adds linear diffusion. 2025. hal-04948845

HAL Id: hal-04948845

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

Preprint submitted on 14 Feb 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Using Sinkhorn in the JKO scheme adds linear diffusion

Aymeric BARADAT*,

Anastasiia HRAIVORONSKA*,

Filippo SANTAMBROGIO*

February 14, 2025

Abstract

The JKO scheme is a time-discrete scheme of implicit Euler type that allows to construct weak solutions of evolution PDEs which have a Wasserstein gradient structure. The purpose of this work is to study the effect of replacing the classical quadratic optimal transport problem by the Schrödinger problem (a.k.a. the entropic regularization of optimal transport, efficiently computed by the Sinkhorn algorithm) at each step of this scheme. We find that if $\epsilon$ is the regularization parameter of the Schrödinger problem, and $\tau$ is the time step parameter, considering the limit $\tau, \epsilon \to 0$ with $\epsilon/\tau \to \alpha \in \mathbb{R}_+$ results in adding the term $\alpha/2 \Delta \rho$ on the right-hand side of the limiting PDE. In the case $\alpha = 0$ we improve a previous result by Carlier, Duval, Peyré and Schmitzer (2017).

# 1 Introduction

**The JKO scheme.** The Jordan-Kinderlehrer-Otto scheme (a.k.a. JKO scheme), introduced for the first time in [17] for linear diffusion equations (Fokker-Planck), then applied in [21] to the case of non-linear diffusions, has been studied in a wide variety of contexts since then. We refer to [2] and [27] for a more detailed presentation. We can say that the JKO scheme is one of the most beautiful applications of optimal transport, as it allows us to solve Cauchy problems – and hence to predict later states of a system – using as a brick the Monge-Kantorovich problem (which involves the optimal interpolation between the configurations at two fixed times).

The idea is that if $\Omega$ is a bounded domain in the Euclidean space or a Riemannian manifold (say for simplicity compact, without boundaries), optimal transport is a way to endow the space of probability measures $\mathcal{P}(\Omega)$ with a formal Riemannian structure. Therefore, given an energy functional $\mathcal{F}: \mathcal{P}(\Omega) \to \mathbb{R}$, it is possible to wonder what is the formal gradient flow of $\mathcal{F}$ in the space of probability measures. Formally, this gradient flow, that we could write:

$$ \dot{\rho} = -\nabla_{\mathcal{P}(\Omega)} \mathcal{F}(\rho), $$

takes the form of the PDE

$$ \partial_t \rho - \mathrm{div} \left( \rho \nabla \left( \frac{\delta \mathcal{F}}{\delta \rho} [\rho] \right) \right) = 0, \quad (1) $$

where $\frac{\delta \mathcal{F}}{\delta \rho}[\rho]$ is the first variation of $\mathcal{F}$ that is the function implicitly defined by the formal condition

$$ \left. \frac{d}{ds} \mathcal{F}(\rho + s\varphi) \right|_{s=0} = \int_{\Omega} \frac{\delta \mathcal{F}}{\delta \rho}[\rho](x) \varphi(x) \, dx, \quad \forall \varphi \in \mathcal{C}(\Omega) \text{ with } \int \varphi(x) \, dx = 0. $$

<footnote>*Universite Claude Bernard Lyon 1, CNRS, Centrale Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, 43 bd du 11 Noembre 1918, 69622 Villeurbanne, France, {baradat,hraivoronska,santambrogio}@math.univ-lyon1.fr</footnote>
As a consequence, many well-known PDEs, such as for instance the heat equation, are formally gradient flows of well-chosen energy functionals (the Boltzmann entropy in the latter case). Therefore, a very natural idea consists in trying to solve these so-called Wasserstein gradient flows using an implicit Euler scheme.

In this context, once $\Omega$ and $\mathcal{F}$ are chosen, given a time step parameter $\tau > 0$ and an initial condition $\rho_0 \in \mathcal{P}(\Omega)$, the JKO scheme is:

$$ \left\{ \begin{array}{l} \rho_0^\tau := \rho_0, \\ \forall n \in \mathbb{N}, \rho_{n+1}^\tau \in \arg\min \left\{ \frac{D(\rho_n^\tau, \rho)^2}{2\tau} + \mathcal{F}(\rho) \middle| \rho \in \mathcal{P}(\Omega) \right\}, \end{array} \right. \qquad (2) $$

where $D$ stands for the 2-Wasserstein distance, that we will define below at Definition 2.1. The goal of the theory is hence to prove that under some assumptions on $\mathcal{F}$, if $(\rho_n^\tau)_{n \in \mathbb{N}}$ satisfies the induction relation (2), the family of piecewise constant curves defined by

$$ \bar{\rho}^{\tau}(t) := \rho_{[t/\tau]}^{\tau} $$

is compact in the set of curves valued in $\mathcal{P}(\Omega)$ as $\tau \to 0$ (i.e. there exists $\tau_k \to 0$ such that the corresponding curves uniformly converge - for the distance $D$ - to a limit curve), and that the corresponding limit points are weak solution of the PDE (1). Once again, we refer to [27] for a review of the available assumptions and techniques allowing to prove such results.

**Our contribution.** In this work, we study the effect of replacing the usual Wasserstein distance in this scheme by its so-called entropic regularization, a.k.a. Schrödinger cost studied for instance in [19, 20]. This quantity is labeled by a positive parameter $\varepsilon$, and we denote it by $D_\varepsilon$. We will define it properly below at Definition 2.3, but let us already mention that it is obtained from the Wasserstein distance either by adding to a transport cost a penalization of the form $\varepsilon$ times the Boltzmann entropy of the transport plan (see formula (7) for more precisions), or equivalently by adding a diffusive term of the form $\frac{\varepsilon}{2}\Delta\rho$ in the right-hand side of the continuity equation in the Benamou-Brenier formulation of optimal transport (which is the formulation that we will use). As we will explain in Subsection 2.2, it is a good approximation of the Wasserstein distance since $D_\varepsilon \to D$ as $\varepsilon \to 0$, and it is very useful for numerics since it can be computed using the so-called Sinkhorn algorithm.

Our main result, precisely stated in Theorem 2.9, can be described as follows. For given $\tau, \varepsilon > 0$, let $(\rho_n^{\tau, \varepsilon})_{n \in \mathbb{N}}$ satisfy the following iterative relations (the entropic JKO scheme):

$$ \left\{ \begin{array}{l} \rho_0^{\tau, \varepsilon} := \rho_0, \\ \forall n \in \mathbb{N}, \rho_{n+1}^{\tau, \varepsilon} \in \arg\min \left\{ \frac{D_\varepsilon(\rho_n^{\tau, \varepsilon}, \rho)^2}{2\tau} + \mathcal{F}(\rho) \middle| \rho \in \mathcal{P}(\Omega) \right\}. \end{array} \right. \qquad (3) $$

Then, let us choose $(\tau_k)$ and $(\varepsilon_k)$ two sequences of positive numbers converging to 0, and such that the ratio $\varepsilon_k/\tau_k$ converges to some $\alpha \in \mathbb{R}_+$ as $k \to +\infty$ (note that we allow $\alpha$ to vanish). Then, under rather usual assumptions for $\mathcal{F}$ that will be detailed in Subsection 2.3, the family of piecewise constant curves

$$ \bar{\rho}^k(t) := \rho_{[t/\tau_k]}^{\tau_k, \varepsilon_k} $$
is compact in the set of continuous curves valued in $\mathcal{P}(\Omega)$ as $k \to +\infty$, and the corresponding limit points are weak solutions of the PDE

$$ \partial_t \rho - \operatorname{div} \left( \rho \nabla \left( \frac{\delta \mathcal{F}}{\delta \rho} [\rho] \right) \right) = \frac{\alpha}{2} \Delta \rho. \qquad (4) $$

In words, we identify the correct regime between $\varepsilon$ and $\tau$ to see the effect of the regularization at the limit $\tau \to 0$, and we show that this effect is the addition of linear diffusion in the limiting PDE.

Using $D_\varepsilon$ instead of $D$ as a way of approximating the JKO scheme is not new. It was, for instance, done in [22], with convergence results proven in [7]. But in all the references we know so far, the entropic regularization was just considered as a technical, computational tool to improve the numerical feasibility, and the idea was that $\varepsilon$ had to be very small. This also impacted the quality of the numerical results, as the Sinkhorn algorithm is much less stable and converges slowly for small $\varepsilon$. We remark that even though the case $\varepsilon/\tau \to 0$ is not the main focus of this paper, our result covers this case and holds without the technical assumption $\varepsilon|\log \varepsilon| \le C\tau^2$ that was present in [7].

What we prove here is that the entropic approximation in the term $D_\varepsilon$ provides linear diffusion in the limit. An early result studying this connection appeared in [1] in the quite restrictive setting corresponding in our language to $\mathcal{F} = 0$, in one dimension. The validity of the same result in more general setting was somehow conjectured because adding to the functional of the gradient flow a Boltzmann entropy means adding linear diffusion to the PDE, but the delicate point is that the Kantorovich formulation of $D_\varepsilon$ involves the entropy of the transport plan, and not of its marginals, see formula (7).

We also note that some recent results exist about the use of $D_\varepsilon$ in gradient flows when $\varepsilon$ is not small at all: this is what is done in [18] and [16], using the so-called *debiasied Sinkhorn divergence*, which consists in defining a distance using

$$ (\mu, \nu) \in \mathcal{P}(\Omega) \times \mathcal{P}(\Omega) \mapsto D_{\varepsilon}^{2}(\mu, \nu) - \frac{1}{2} D_{\varepsilon}^{2}(\mu, \mu) - \frac{1}{2} D_{\varepsilon}^{2}(\nu, \nu), $$

a quantity which contrarily to $D_\varepsilon$ vanishes if $\mu = \nu$ and which is a convex, smooth, positive definite loss function that metrizes the weak convergence of probability measures on a compact set, see [23, 13, 11]. However, the spirit of [18] and [16] is quite different from what we do in this paper, as the goal in [16] is to obtain a new PDE, corresponding to the differential structure of the space of probability measures endowed with the geometry of debiased Sinkhorn divergence studied in [18], for fixed $\varepsilon > 0$.

**Context and notations.** In this article, we will always work in the flat torus, that is, $\Omega = \mathbb{T}^d$. Here are a few notations that we will constantly use.

* We denote by Leb the Lebesgue measure on $\mathbb{T}^d$, normalized so that $\mathrm{Leb}(\mathbb{T}^d) = 1$.

* We call $(\sigma_t)_t \ge 0$ the heat kernel on the torus, that is, the solution of

$$ \left\{ \begin{array}{l} \partial_t \sigma_t = \frac{1}{2} \Delta \sigma_t, \\ \sigma_0 = \delta_0. \end{array} \right. $$
• Whenever we consider a function $f = f(t, x)$ depending on time and space, $f(t)$ stands for the map $x \mapsto f(t, x)$ depending on space only.

**Outline of the paper.** In Section 2, we define properly the Wasserstein distance and the Schrödinger cost, explain the link between the two and the role of the Sinkhorn algorithm to compute the second, present the assumptions we use and state the main result of this work. In Section 3, we provide the properties of the Schrödinger cost that we need in the proof of Theorem 2.9. Most of the results included in this section are classical, and we will often only provide references to the relevant literature for proofs, but we also include some formal proofs to ensure the text remains self-contained and accessible. In Section 4, we detail the proof of our main result, Theorem 2.9.

# 2 Key notions and main result

## 2.1 The Wasserstein distance and its entropic regularization

We have not yet described what is the Wasserstein distance we use. There are several ways to define it, but in the present article, we will mainly work in so-called *Benamou-Brenier formulations* [5]. According to this choice, the Wasserstein distance is defined as follows. See [12, Section 4.1] for a more detailed presentation.

**Definition 2.1.** Let $\mu, \nu \in \mathcal{P}(\mathbb{T}^d)$. The Wasserstein distance between $\mu$ and $\nu$ is defined by

$$ \frac{D(\mu, \nu)^2}{2} := \inf \left\{ \int_0^1 \int \frac{|v|^2}{2} \rho \, dt \left| \begin{array}{l} \partial_t \rho + \operatorname{div}(\rho v) = 0, \\ \rho_0 = \mu, \rho_1 = \nu \end{array} \right. \right\}. $$

*Remark 2.2.* This quantity is indeed a distance, and it will actually be our reference distance for $\mathcal{P}(\mathbb{T}^d)$. Namely, when dealing with $\mathcal{P}(\mathbb{T}^d)$, we always consider it as a metric space endowed with $D$. Note that the topology induced by $D$ is the one of weak convergence (that is, in duality with the set of continuous functions).

A rescaling in time in the continuity equation ensures the following equivalent formula which depends on a parameter $\tau > 0$ and which is often useful when studying the JKO scheme:

$$ \frac{D(\mu, \nu)^2}{2\tau} = \inf \left\{ \int_0^\tau \int \frac{|v|^2}{2} \rho \, dt \left| \begin{array}{l} \partial_t \rho + \operatorname{div}(\rho v) = 0, \\ \rho_0 = \mu, \rho_\tau = \nu \end{array} \right. \right\}. $$

In this work, our goal is to study the JKO scheme when we replace this Wasserstein distance by an entropic regularized version of it, implying the Schrödinger problem instead of the quadratic optimal transport problem. In Benamou-Brenier formulation, this entropic regularization is defined as follows.

**Definition 2.3.** Let $\mu, \nu \in \mathcal{P}(\mathbb{T}^d)$ and $\varepsilon > 0$. The Schrödinger cost of parameter $\varepsilon$ between $\mu$ and $\nu$ is

$$ \frac{D_{\varepsilon}(\mu, \nu)^2}{2} := \inf \left\{ \varepsilon H(\mu) + \int_0^1 \int \frac{|v|^2}{2} \rho \, dt \left| \begin{array}{l} \partial_t \rho + \operatorname{div}(\rho v) = \frac{\varepsilon}{2} \Delta \rho, \\ \rho_0 = \mu, \rho_1 = \nu \end{array} \right. \right\}, \quad (5) $$
where $H: \mathcal{P}(\mathbb{T}^d) \to [0, +\infty]$ is the Boltzmann entropy, defined via

$$H(\mu) := \begin{cases} \int_{\mathbb{T}^d} \log(\mu(x)) \mu(x) \, dx & \text{if } \mu \ll \text{Leb}, \\ + \infty & \text{otherwise,} \end{cases}$$

where we identified measures with their densities with respect to Leb. As $\text{Leb}(\mathbb{T}^d) = 1$, this quantity is always nonnegative.

**Remark 2.4.** We call the quantity $D_\varepsilon$ a cost and not a distance because it is not a distance. For instance, $D_\varepsilon(\mu, \mu) \neq 0$, unless $\mu = \text{Leb}$. On the other hand, the Schrödinger cost is symmetric as we will recall in Corollary 3.3. Obtaining this convenient symmetry is the reason for the term $\varepsilon H(\mu)$ in (5).

There are several ways to define the Schrödinger cost between $\mu$ and $\nu$, even once chosen to work in Benamou-Brenier formulations. In our work, the formulation above is the most useful one, at least for the formal computations. Yet, we will recall in Proposition 3.2 below other equivalent formulations.

Once again, a time rescaling provides the following equivalent formulation for the Schrödinger cost:

$$\frac{D_\varepsilon(\mu, \nu)^2}{2\tau} = \inf \left\{ \frac{\varepsilon}{\tau} H(\mu) + \int_0^\tau \int \frac{|v|^2}{2} \rho \, dt \left| \begin{array}{l} \partial_t \rho + \text{div}(\rho v) = \frac{\varepsilon}{2\tau} \Delta \rho, \\ \rho_0 = \mu, \rho_\tau = \nu \end{array} \right. \right\}. \quad (6)$$

We can already see why formulation (5) is so practical for proving our result: the linear diffusion added in the right hand side of equation (4) is precisely the one added to the continuity equation in (6).

## 2.2 Static and dynamic entropic regularization of optimal transport problems

The goal of this section is to explain the interest in using the Schrödinger cost in place of the Wasserstein distance. The main interest is indeed computational.

First of all, as observed in [14] (see also [15] for a more abstract result), the dynamical formulation of the Schrödinger problem and the associated cost $D_\varepsilon$ is the Benamou-Brenier version of the problem of minimizing the relative entropy of the law of a stochastic process with given initial and final marginals $\mu$ and $\nu$ with respect to the law of the reversible Brownian motion, that is, whose initial law is the Lebesgue measure. In other terms, following the celebrated Sanov theorem [25], the Schrödinger problem can be interpreted as a large deviation problem [20].

Also, a short reasoning involving the additivity of the logarithm implies that the latter problem can be reformulated in terms of a “static” entropy minimization. Let us consider all probability measures $\gamma \in \mathcal{P}(\Omega \times \Omega)$ whose marginal measures are given by $(\pi_x)_\#\gamma = \mu$ and $(\pi_y)_\#\gamma = \nu$. These probability measures on the product space are called transport plans in the Monge-Kantorovich theory and the set of transport plans is often denoted by $\Pi(\mu, \nu)$. The reformulation in [14] writes

$$\frac{1}{2} D_{\varepsilon}^{2}(\mu, \nu) = \min\{\varepsilon H(\gamma|R_{\varepsilon}) : \gamma \in \Pi(\mu, \nu)\},$$

where $H(\cdot|\cdot)$ denotes the relative entropy and $R_\varepsilon$ is the joint law of a random variable of the form $(B_0, B_\varepsilon)$, $(B_t)_{t \ge 0}$ being a reversible Brownian motion. Knowing the transition kernel of the
Brownian motion, let us define the cost $c_{\varepsilon}(x, y)$ through the formula:

$$ \frac{dR_{\varepsilon}}{dx \otimes dy}(x, y) = \sigma_{\varepsilon}(y - x) := \frac{1}{\sqrt{2\pi\varepsilon^d}} \exp\left(-\frac{c_{\varepsilon}(x, y)}{\varepsilon}\right). $$

Then, it is possible to see that our problem can be reformulated as follows:

$$ \frac{1}{2} D_{\varepsilon}^{2}(\mu, \nu) = \min \left\{ \int c_{\varepsilon} d\gamma + \varepsilon H(\gamma) + \frac{d\varepsilon}{2} \log(2\pi\varepsilon) : \gamma \in \Pi(\mu, \nu) \right\}, \quad (7) $$

where the entropy is now a Boltzmann entropy for measures on the product space. This is the so-called entropic regularization of the optimal transport problem with cost $c_{\varepsilon}$, as forgetting about the constant term which has an effect on the optimal value but not on the optimizers, it consists in adding a penalization of the form of an entropy to an optimal transport problem.

Solving problems of this type can be attacked numerically via the Sinkhorn algorithm, also called iterative proportional fitting procedure (IPFP), which can be seen as an alternate maximization on the dual problem or as alternate projections for the Kullback-Leibler divergence in terms of the primal one. This was introduced long ago in [28, 29], and then popularized in optimal transport in [8] and [6]. This iterative algorithm is incredibly efficient, especially when $\varepsilon$ is not too small.

It is interesting to observe that, if we replace the torus $\mathbb{T}^d$ with the Euclidean space $\mathbb{R}^d$ (which leads to problems since the Lebesgue measures is not a probability measure on $\mathbb{R}^d$, but let us not comment on this technical issue), then we obtain the following explicit expression for $c_{\varepsilon}$, independent of $\varepsilon$:

$$ c_{\varepsilon}(x, y) = \frac{|x - y|^2}{2}. \quad (8) $$

This makes it clear that, in the limit $\varepsilon \to 0$, the computed value should converge to the squared Wasserstein distance (in its Kantorovich formulation, and not in the dynamic formulation that we presented here). This is rigorously established in [7] under the necessary and sufficient condition that $H(\mu)$ and $H(\nu)$ are finite (otherwise, $D_{\varepsilon}$ is infinite whatever $\varepsilon > 0$).

In the case of the torus, identity (8) does not hold anymore, but is replaced by

$$ c_{\varepsilon}(x, y) = -\varepsilon \log \left( \sum_{k \in \mathbb{Z}^d} \exp \left( -\frac{|Y + k - X|^2}{2\varepsilon} \right) \right), \quad (9) $$

where $X$ and $Y$ are any lifts of $x$ and $y$ from $\mathbb{T}^d$ to $\mathbb{R}^d$. Still:

* The convergence of $D_{\varepsilon}$ towards $D$ at the level of the Kantorovich formulation can be proved using the ideas of [7] in virtue of the uniform convergence

$$ c_{\varepsilon}(x, y) \xrightarrow{\varepsilon \to 0} \frac{d(x, y)^2}{2}, $$

where $d$ is the canonical distance on $\mathbb{T}^d$.

* The Sinkhorn algorithm can be used efficiently to compute $D_{\varepsilon}$, as this efficiency is not affected by replacing the squared distance with $c_{\varepsilon}$. The only difference with the standard case is that one needs to compute $c_{\varepsilon}$ once for all, either by using formula (9), or by solving the heat equation.
* We believe that our result remains the same replacing $D_\epsilon$ by the quantity (7) in which $c_\epsilon$ is replaced by $d(x, y)^2/2$, but we did not try to prove this fact. Once our result is known, such a proof would reduce to establishing some quantitative stability in the entropic regularization of optimal transport with respect to the cost.

The same could also be done in bounded domains using reflected Brownian motions and the heat kernel with Neumann boundary conditions, but we decided not to insist in this paper on the boundary issues (in general, we should impose everywhere no-flux boundary conditions).

We close this section by noticing that the proof of convergence $D_\epsilon \to D$ under the condition that $H(\mu)$ and $H(\nu)$ are finite is also possible to perform at the level of the dynamical formulations [19, 4].

With all these observations, we believe that it is natural to study the effect of replacing $D$ by its regularized version in the JKO scheme, when letting both $\epsilon$ and $\tau$ go to 0 in an appropriate regime.

## 2.3 Precise statement of the main result

In this article, we consider energy functionals $\mathcal{F}$ of the form

$$ \mathcal{F}(\rho) := \begin{cases} \int \left\{ V(x)\rho(x) + \frac{1}{2}W * \rho(x)\rho(x) + f(\rho(x)) \right\} dx, & \text{if } \rho \ll \text{Leb}, \\ +\infty, & \text{otherwise,} \end{cases} \qquad (10) $$

where $V, W: \mathbb{T}^d \to \mathbb{R}$ and $f: \mathbb{R}_+ \to \mathbb{R}$ satisfy assumptions to be specified in a few lines. At the formal level, it can be checked that for a given $\rho \in \mathcal{P}(\mathbb{T}^d)$ with $\rho \ll \text{Leb}$,

$$ \frac{\delta \mathcal{F}}{\delta \rho}[\rho] = V + W * \rho + f'(\rho), $$

so that (1) writes

$$ \partial_t \rho - \operatorname{div}(\rho(\nabla V + \nabla W * \rho)) = \Delta g(\rho), $$

where $g: \mathbb{R}_+ \to \mathbb{R}$ is the function defined for all $s \in \mathbb{R}_+$ by

$$ g(s) = sf'(s) - f(s). \qquad (11) $$

Of course, equation (4) then becomes

$$ \partial_t \rho - \operatorname{div}(\rho(\nabla V + \nabla W * \rho)) = \Delta g(\rho) + \frac{\alpha}{2} \Delta \rho. \qquad (12) $$

### 2.3.1 Assumptions on the internal energy potential

To state our main result, we will need some assumptions on the function $f$ appearing in the definition of $\mathcal{F}$ in (10). Our goal is to consider rather general functions $f$, imposing conditions on its smoothness and growth, but not relying on some specific kind of nonlinearity. The general assumptions we choose are clearly not optimal (see Remark 2.8) but lead to rather direct proofs. Two different assumptions will appear at different stages of the proof of Theorem 2.9.

The first one is necessary for proving that for positive $\tau, \epsilon$, equation (12) is satisfied up to a term converging to 0 as $\tau, \epsilon \to 0$. The second one concerns strong compactness and is therefore necessary to prove that the limiting curve is a distributional solution of the (nonlinear) limiting PDE (12).
**Assumption 2.5.** The function $f: \mathbb{R}_+ \to \mathbb{R} \cup \{+\infty\}$ is convex, and the domain where it is finite $\text{Dom}(f) := \{f < +\infty\}$ is a nonempty and nontrivial interval of $\mathbb{R}_+$. In other words, calling $d_-, d_+ \in \mathbb{R}_+ \cup \{+\infty\}$ the lower and upper bound of $\text{Dom}(f)$, we have $d_- < d_+$. Moreover, we assume that $f$ is continuous on $\text{Dom}(f)$, of class $C^3$ on its interior $(d_-, d_+)$, and that

$$ \text{if } d_- > 0, \text{ then } f'(s) \underset{s \to d_-}{\longrightarrow} -\infty \quad \text{and} \quad \text{if } d_+ < +\infty, \text{ then } f'(s) \underset{s \to d_+}{\longrightarrow} +\infty. \quad (13) $$

Finally, we assume that $d_- < 1 < d_+$.

**Remark 2.6.** Under Assumption 2.5, the function $g$ defined by (11) is well defined on $(d_-, d_+)$. As, it is increasing, we can extend it with potentially infinite values at $d_-$ and $d_+$.

The Assumption $d_- < 1 < d_+$ is very natural. Indeed, in this paper, we consider measures $\rho \in \mathcal{P}(\mathbb{T}^d)$ such that $\mathcal{F}(\rho) < +\infty$, that is, measures with densities satisfying

$$ \int f(\rho) \, dx < +\infty. $$

But by Jensen's inequality, if $f$ is convex,

$$ \int f(\rho) \geq f(1). $$

So in particular, for such a $\rho$ to exist, $1 \in \text{Dom}(f) \subset [d_-, d_+]$. In addition, under condition (13), it is easy to see that if $1 \in \{d_-, d_+\}$, then the Lebesgue measure is the only measure for which $\mathcal{F}$ is finite. We do not want to work in this situation, and therefore we assume $d_- < 1 < d_+$.

Our second assumption is the following.

**Assumption 2.7.** Either $\alpha > 0$ and we have nothing more to assume. Or $\alpha = 0$, and then, we need two assumptions. First, we need the following growth condition: either the upper bound $d_+$ is finite, or the function $h$ defined for all $s \in (d_-, d_+)$ by

$$ h(s) := \int_{1}^{s} \sqrt{r} f''(r) \, dr, $$

satisfies

$$ \frac{h(s)}{s^{\frac{d-2}{2d}}} \underset{s \to +\infty}{\longrightarrow} +\infty. \quad (14) $$

Second, we assume that at least one of the following conditions holds

(i) $f$ is superlinear,

(ii) $d_- = 0$,

(iii) $f(d_-) = +\infty$.

**Remark 2.8.** The fact that we do not need any further assumption in the case when $\alpha > 0$ can be understood as follows. In that case, the strong compactness needed to pass to the limit directly comes from the regularity properties of the Schrödinger problem. On the other hand, when $\alpha = 0$, the strong compactness has to be provided by the limit PDE, and we have to make assumptions similar to the ones of the non-entropic JKO scheme. We chose conditions that provide a rather
simple proof in the classical case, but which are not optimal. We believe that any assumption which would make the classical JKO scheme converge would be enough to prove our convergence result with $\alpha = 0$, but we wanted to stay in a rather easy setting in order to emphasize the arguments specific to the entropic case.

A quick computation shows that if $f: s \mapsto s \log s$, then $h(s) = \mathcal{O}_{s \to +\infty}(s^{1/2})$, so the condition (14) is satisfied. In the case where $f: s \mapsto s^m/(m-1)$, $m \neq 1$, that is the porous medium or fast diffusion case, then $h(s) = \mathcal{O}_{s \to +\infty}(s^{m-1/2})$. The assumption becomes $m > 1 - 1/d$, which is essentially the same as [2, Condition (10.4.21)], used in [2, Example 11.2.4]. Here, we will provide a straightforward proof.

### 2.3.2 Precise statement of the main result

We are now in a position to state our main result.

**Theorem 2.9.** Let us assume that $V$ and $W$ are of class $C^2$, that $f$ satisfies Assumption 2.5, and let us define $\mathcal{F}$ as in (10) and $g$ as in (11). Let us consider two sequences of positive numbers $(\tau_k)_{k \in \mathbb{N}}$ and $(\varepsilon_k)_{k \in \mathbb{N}}$ such that

$$ \tau_k \underset{k \to +\infty}{\longrightarrow} 0 \quad \text{and} \quad \alpha_k := \frac{\varepsilon_k}{\tau_k} \underset{k \to +\infty}{\longrightarrow} \alpha \in [0, +\infty). $$

Let $\rho_0 \in \mathcal{P}(\mathbb{T}^d)$ satisfy $\mathcal{F}(\rho_0) < +\infty$ and $H(\rho_0) < +\infty$. Let $(\rho_n^k)_{n \in \mathbb{N}}$ be a sequence satisfying the inductive relation (3) with $\varepsilon = \varepsilon_k$ and $\tau = \tau_k$. Finally, let us define for all $t \ge 0$

$$ \bar{\rho}^k(t) := \rho_{[t/\tau_k]}^k. $$

Then, for any $T > 0$, the following assertions hold:

* The sequence $(\bar{\rho}^k)$ is compact for the uniform convergence of functions on $[0, T]$ valued into $\mathcal{P}(\mathbb{T}^d)$, when this set of probabilities is endowed with the distance $D$. Moreover, its limit points belong to $C([0, T]; \mathcal{P}(\mathbb{T}^d))$.

* For all $k \in \mathbb{N}$ and $t \in [0, T]$, $\bar{\rho}^k(t)$ has a density, $g(\bar{\rho}_k)$ belongs to $L^1([0, T] \times \mathbb{T}^d)$ and the following convergence holds in $D'((0, T) \times \mathbb{T}^d)$:

$$ \partial_t \bar{\rho}^k - \operatorname{div}(\bar{\rho}^k(\nabla V + \nabla W * \bar{\rho}^k)) - \Delta g(\bar{\rho}^k) - \frac{\alpha_k}{2} \Delta \bar{\rho}^k \underset{k \to +\infty}{\longrightarrow} 0. \quad (15) $$

* If in addition, $f$ satisfies Assumption 2.7, any limit point $\rho$ of $(\bar{\rho}^k)$ given by the first point has a density, $g(\rho)$ belongs to $L^1([0, T] \times \mathbb{T}^d)$, and $\rho$ is a distributional solution of

$$ \left\{ \begin{array}{l} \partial_t \rho - \operatorname{div}(\rho(\nabla V + \nabla W * \rho)) = \Delta g(\rho) + \frac{\alpha}{2} \Delta \rho, \\ \\ \rho|_{t=0} = \rho_0. \end{array} \right. \quad (16) $$

**Remark 2.10.** In the case when $\alpha = 0$, that is, when the regularization parameter converges to 0 faster than the time step, we improve a previous result from [7] where the authors had to assume the much stronger condition $\varepsilon_k |\log \varepsilon_k| = \mathcal{O}_{k \to +\infty}(\tau_k^2)$.

Our result shows that to get a linear diffusion term $\Delta\rho$ in the limiting PDE, we can either act at the level of the energy functional by including some Boltzmann entropy (if $f(s) = s \log s$, then
$g(s) = s$), act at the level of the "distance" part of the scheme by keeping $\alpha > 0$, or mix both possibilities. As it is easier to compute Schrödinger costs than Wasserstein distances, especially for large values of $\varepsilon$, it seems relevant to choose to act at the level of the distance only. Yet, a rigorous comparison of the two approaches would also need to know the convergence rates for both schemes: it is useful to compute each time step more quickly only provided the convergence rate is not dramatically worsened.

The picky reader could observe that the case $\alpha = 0$ and $f = 0$ is not covered in this theorem. Of course, in that case, the convergence is even easier as strong compactness for $\rho$ is not needed to justify the limiting PDE.

### 3 Preliminaries on the Schrödinger problem

In this section, we gather a few standard results about the Schrödinger problem that will be useful in the proof of Theorem 2.9. In that proof, there will always be a rescaling by $\tau$ as in (6), so the only somehow unusual aspect of our presentation is that we will carry this rescaling all along to be able to use the results directly later. Of course, to get more standard formulations, it suffices to take $\tau = 1$. We give sketches of proofs when we think it can be useful, but mainly give references to where the results can be found in the literature. Only Proposition 3.9 is mainly useful in the context of JKO, and hence somehow new, although expected. Therefore, for this last proposition, we provide a complete proof.

First, the Schrödinger cost is convex and lower semi-continuous in the topology of narrow convergence.

**Proposition 3.1.** Let $\varepsilon > 0$. The quantities $(\mu, \nu) \mapsto D_{\varepsilon}^2(\mu, \nu)$ and $(\mu, \nu) \mapsto D_{\varepsilon}^2(\mu, \nu)/2 - \varepsilon H(\mu)$ are convex and lower semi-continuous in the topology of narrow convergence.

*Formal proof.* This is a direct consequence of the convexity and lower semi-continuity of the cost and linearity of the constraints in (5), when considered in the variables $(\rho, \rho\nu)$. See [4, Section 2.1] for more details. □

Actually, the Schrödinger problem has several formulations, even in the Benamou-Brenier style. The one we presented in Definition 2.3 is forward in time in the sense that the density follows a forward parabolic equation, but there is another formulation that is backward in time, and a third one that is symmetric in time. We refer to [20, 14, 15] for the rigorous proofs.

**Proposition 3.2.** Let $\mu, \nu \in \mathcal{P}(\mathbb{T}^d)$ and $\varepsilon, \tau > 0$. We have

$$ \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} = \inf \left\{ \frac{\varepsilon}{\tau} H(\nu) + \int_0^\tau \int \frac{|w|^2}{2} \rho \, dt \left| \begin{array}{l} \partial_t \rho + \operatorname{div}(\rho w) = -\frac{\varepsilon}{2\tau} \Delta \rho, \\ \rho_0 = \mu, \rho_\tau = \nu \end{array} \right. \right\} \quad (17) $$

$$ = \inf \left\{ \frac{\varepsilon H(\mu) + H(\nu)}{2} + \frac{1}{2} \int_0^\tau \int \left[ |c|^2 + \left| \frac{\varepsilon}{2\tau} \nabla \log \rho \right|^2 \right] \rho \, dt \left| \begin{array}{l} \partial_t \rho + \operatorname{div}(\rho c) = 0, \\ \rho_0 = \mu, \rho_\tau = \nu \end{array} \right. \right\}. \quad (18) $$

*Formal proof.* Let us explain the main idea without checking all the regularity needed to justify the computations. Let us consider $(\rho, \nu)$ a competitor for (6) and call

$$ c := v - \frac{\varepsilon}{2\tau} \nabla \log \rho, \quad w := v - \frac{\varepsilon}{\tau} \nabla \log \rho. $$
The first observation is that

$$ \partial_t \rho + \operatorname{div}(\rho c) = 0 \quad \text{and} \quad \partial_t \rho + \operatorname{div}(\rho w) = -\frac{\varepsilon}{2\tau} \Delta \rho. \qquad (19) $$

In addition, we have

$$ \frac{1}{2} \int_{0}^{\tau} \int |v|^2 \rho \, dt = \frac{1}{2} \int_{0}^{\tau} \int |c|^2 \rho \, dt + \frac{1}{2} \int_{0}^{\tau} \int \left| \frac{\varepsilon}{2\tau} \nabla \log \rho \right|^2 \rho \, dt + \frac{\varepsilon}{2\tau} \int_{0}^{\tau} \int \{c \cdot \nabla \log \rho\} \rho \, dt. $$

On the other hand, equation (19) implies

$$ \frac{d}{dt} \int \rho \log \rho = \int \{c \cdot \nabla \log \rho\} \rho. $$

Therefore,

$$ \frac{\varepsilon}{2\tau} \int_{0}^{\tau} \int \{c \cdot \nabla \log \rho\} \rho \, dt = \frac{\varepsilon}{2\tau} (H(\nu) - H(\mu)). $$

Identity (18) follows, and (17) is proved in the same way. □

As a corollary of this result, we obtain the symmetry of $D_\varepsilon$.

**Corollary 3.3.** For any $\mu, \nu \in \mathcal{P}(\mathbb{T}^d)$ and $\varepsilon > 0$, $D_\varepsilon(\mu, \nu) = D_\varepsilon(\nu, \mu)$.

*Proof.* Formulation (18) is clearly invariant through time reversal. □

The next result is concerned with the existence of optimizers for the minimization problems in play. It can be deduced from the existence results in [20] or [4].

**Proposition 3.4.** Let $\mu, \nu \in \mathcal{P}(\mathbb{T}^d)$ and $\varepsilon, \tau > 0$. The Schrödinger cost $D_\varepsilon(\mu, \nu)$ as defined in Definition 2.3 is finite if and only if $H(\mu)$ and $H(\nu)$ are finite. In that case, there exist unique minimizers $(\rho, v)$, $(\rho, c)$ and $(\rho, w)$ achieving the infimum in (5), (17) and (18) respectively, of common density. In particular

$$ \frac{\varepsilon}{\tau} H(\mu) + \int_{0}^{\tau} \int \frac{|v|^2}{2} \rho \, dt = \frac{\varepsilon}{\tau} \frac{H(\mu) + H(\nu)}{2} + \frac{1}{2} \int_{0}^{\tau} \int \left\{ |c|^2 + \left| \frac{\varepsilon}{2\tau} \nabla \log \rho \right|^2 \right\} \rho \, dt = \frac{\varepsilon}{\tau} H(\nu) + \int_{0}^{\tau} \int \frac{|w|^2}{2} \rho \, dt. \qquad (20) $$

In addition, $\rho$ is in $AC^2([0, \tau]; \mathcal{P}(\mathbb{T}^d))$, with metric derivative

$$ |\dot{\rho}(t)|^2 = \int |c(t)|^2 \rho(t), \quad \text{for a.e. } t \in [0, \tau]. $$

Finally, we have almost everywhere

$$ v - \frac{\varepsilon}{2\tau} \nabla \log \rho = c = w + \frac{\varepsilon}{2\tau} \nabla \log \rho. \qquad (21) $$

The next proposition is a duality result, which is a direct application of the results in [9] once the equivalence between static and Benamou-Brenier formulations of the Schrödinger problem – as presented in [14] – is understood. Here, we will propose a very different formal proof based on duality arguments starting from the Benamou-Brenier formulation.
**Proposition 3.5.** Let $\mu, \nu \in \mathcal{P}(\mathbb{T}^d)$ and $\varepsilon, \tau > 0$. Then

$$ \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} = \frac{\varepsilon}{\tau} H(\mu) + \sup_{\varphi \in L^1(\nu)} \langle \varphi, \nu \rangle - \left\langle \frac{\varepsilon}{\tau} \log \left( \left\{ \exp \left( \frac{\tau\varphi}{\varepsilon} \right) \right\} * \sigma_{\varepsilon} \right), \mu \right\rangle. \quad (22) $$

Moreover, as soon as $H(\mu)$ and $H(\nu)$ are finite, the maximization problem on the right-hand side admits a maximizer $\varphi$, and if we call

$$ \psi := \frac{\varepsilon}{\tau} \log \mu - \frac{\varepsilon}{\tau} \log \left( \left\{ \exp \left( \frac{\tau \varphi}{\varepsilon} \right) \right\} * \sigma_{\varepsilon} \right), \quad (23) $$

we have

$$ \varphi = \frac{\varepsilon}{\tau} \log \nu - \frac{\varepsilon}{\tau} \log \left( \left\{ \exp \left( \frac{\tau \psi}{\varepsilon} \right) \right\} * \sigma_{\varepsilon} \right). \quad (24) $$

Stated differently, $\varphi$ satisfies the fixed point condition:

$$ \varphi = \frac{\varepsilon}{\tau} \log \nu - \frac{\varepsilon}{\tau} \log \left( \left\{ \frac{\mu}{\exp\left(\frac{\tau}{\varepsilon}\varphi\right) * \sigma_{\varepsilon}} \right\} * \sigma_{\varepsilon} \right). \quad (25) $$

In particular, $\varphi - \frac{\varepsilon}{\tau} \log \nu$ is of class $C^{\infty}$.

Finally, there is a unique $\varphi$ satisfying (25) up to an additive constant and hence a unique maximizer in (22) up to an additive constant.

**Remark 3.6.** In (22), it looks like we broke the symmetry of the problem. This is because as we will see in the sketch of proof, these formulas are obtained by looking at the dual of the Schrödinger problem formulated as in (6). If we were starting from (17), we would obtain a symmetric dual problem.

**Formal proof.** Let us write problem (6) formally as an inf-sup problem, in the spirit of [3, Section 1.4]:

$$ \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} = \frac{\varepsilon}{\tau} H(\mu) + \inf_{(\rho, v)} \sup_{\varphi} \int_0^{\tau} \int \frac{|v|^2}{2} \rho \, dt \\ + \langle \varphi(\tau), \nu \rangle - \langle \varphi(0), \mu \rangle - \int_0^{\tau} \langle \partial_t \varphi + v \cdot \nabla \varphi + \frac{\varepsilon}{2\tau} \Delta \varphi, \rho \rangle \, dt, $$

where forgetting any regularity issue, the inf is taken over functions $\rho : [0, \tau] \times \mathbb{T}^d \to \mathbb{R}_+$ and $v : [0, \tau] \times \mathbb{T}^d \to \mathbb{R}^d$, and the sup is taken over $\varphi : [0, \tau] \times \mathbb{T}^d \to \mathbb{R}$. As the functional inside the inf-sup is convex in $(\rho, \rho v)$ and concave (actually, affine) in $\varphi$, it seems legitimate to swap the inf and the sup, ending up with:

$$ \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} = \frac{\varepsilon}{\tau} H(\mu) + \sup_{\varphi} \langle \varphi(\tau), \nu \rangle - \langle \varphi(0), \mu \rangle \\ + \inf_{(\rho, v)} \int_0^{\tau} \int \frac{|v|^2}{2} \rho \, dt - \int_0^{\tau} \langle \partial_t \varphi + v \cdot \nabla \varphi + \frac{\varepsilon}{2\tau} \Delta \varphi, \rho \rangle \, dt. $$
Minimizing in $v$ once $\varphi$ is fixed, we find $v = \nabla\varphi$, and then minimizing in $\rho$, we find that the minimizer is $-\infty$ as soon as the inequality $\partial_t\varphi + \frac{1}{2}|\nabla\varphi|^2 + \frac{\varepsilon}{2\tau}\Delta\varphi \le 0$ does not hold everywhere, and cancels in the opposite case. We conclude that we have

$$ \frac{D_{\varepsilon}(\mu, v)^2}{2\tau} = \frac{\varepsilon}{\tau}H(\mu) + \sup_{\substack{\varphi \text{ s.t.} \\ \partial_t\varphi + \frac{1}{2}|\nabla\varphi|^2 + \frac{\varepsilon}{2\tau}\Delta\varphi \le 0}} \langle\varphi(\tau), v\rangle - \langle\varphi(0), \mu\rangle. \quad (26) $$

Now, given a competitor $\varphi$ satisfying the constraint, let us call

$$ \theta := \exp\left(\frac{\tau\varphi}{\varepsilon}\right) \quad \text{so that} \quad \partial_t\theta + \frac{\varepsilon}{2\tau}\Delta\theta \le 0, $$

and the problem becomes

$$ \frac{D_{\varepsilon}(\mu, v)^2}{2\tau} = \frac{\varepsilon}{\tau}H(\mu) + \sup_{\substack{\theta \text{ s.t.} \\ \partial_t\theta + \frac{\varepsilon}{2\tau}\Delta\theta \le 0}} \left\langle \frac{\varepsilon}{\tau} \log \theta(\tau), v \right\rangle - \left\langle \frac{\varepsilon}{\tau} \log \theta(0), \mu \right\rangle. $$

Due to the sign in the constraint, given a final condition $\theta(\tau)$, we increase the value of the right-hand side by replacing $\theta(0)$ by $\theta(\tau) * \sigma_{\varepsilon}$. But this change corresponds to considering only the equality case in the constraint, which is still admissible. This means that we do not change the value of the problem by replacing the inequality sign with an equality sign in the constraint, and we deduce

$$ \begin{aligned} \frac{D_{\varepsilon}(\mu, v)^2}{2\tau} &= \frac{\varepsilon}{\tau}H(\mu) + \sup_{\theta:\mathbb{T}^d\to\mathbb{R}} \left\langle \frac{\varepsilon}{\tau}\log\theta, v \right\rangle - \left\langle \frac{\varepsilon}{\tau}\log(\theta * \sigma_{\varepsilon}), \mu \right\rangle \\ &= \frac{\varepsilon}{\tau}H(\mu) + \sup_{\varphi:\mathbb{T}^d\to\mathbb{R}} \left\langle \varphi, v \right\rangle - \left\langle \frac{\varepsilon}{\tau}\log\left(\exp\left(\frac{\tau\varphi}{\varepsilon}\right) * \sigma_{\varepsilon}\right), \mu \right\rangle, \end{aligned} $$

where in the second line, we changed the variable according to $\varphi = \frac{\varepsilon}{\tau}\log\theta$. This leads to (22).

The final step consists in deriving the optimality condition for $\varphi$ in this last formula. A direct first-order computation gives

$$ \nu = \exp\left(\frac{\tau}{\varepsilon}\varphi\right) \left( \frac{\mu}{\exp\left(\frac{\tau}{\varepsilon}\varphi\right) * \sigma_{\varepsilon}} \right) * \sigma_{\varepsilon}, $$

which implies (25) by rearranging the terms. In particular, defining $\psi$ as in (23), we get (24). The regularity result follows directly from these formulas.

Uniqueness is obtained by showing that (25) admits a unique solution up to an additive constant, and we will not comment on that issue further. $\square$

In the next proposition, we explicit the link between the so-called Schrödinger potentials $\varphi$ and $\psi$ obtained in Proposition 3.5, and the optimal pairs $(\rho, v)$ and $(\rho, w)$ from Proposition 3.4. These properties can be found in [20].

**Proposition 3.7.** Let $\mu, v \in \mathcal{P}(\mathbb{T}^d)$ with $H(\mu)$ and $H(v)$ finite, and $\varepsilon, \tau > 0$. Let $(\rho, v)$ and $(\rho, w)$ be the minimizers from (6) and (17) given by Proposition 3.4. Let $\bar{\varphi}$ be a maximizer of (22) given by Proposition 3.5 and $\bar{\psi}$ defined as in (23). For any $t \in [0, \tau]$, call

$$ \varphi(t) := \frac{\varepsilon}{\tau} \log \left( \exp \left( \frac{\tau \bar{\varphi}}{\varepsilon} \right) * \sigma_{\frac{\varepsilon}{\tau}(\tau-t)} \right) \quad \text{and} \quad \psi(t) := \frac{\varepsilon}{\tau} \log \left( \exp \left( \frac{\tau \bar{\psi}}{\varepsilon} \right) * \sigma_{\frac{\varepsilon}{\tau}t} \right). \quad (27) $$
Then, we have

$$ \text{for almost all } t \in [0, \tau], \qquad v(t) = \nabla \varphi(t) \quad \text{and} \quad w(t) = \nabla \psi(t) \qquad (28) $$

$$ \text{for all } t \in [0, \tau], \qquad \rho(t) = \exp\left(\frac{\tau}{\varepsilon}(\varphi(t) - \psi(t))\right). \qquad (29) $$

In particular, up to modifying $v$ and $w$ on a negligible set of times, $\rho$, $v$ and $w$ are of class $C^\infty$ on $(0, \tau) \times \mathbb{T}^d$.

Finally,

$$ \varphi(\tau) = \bar{\varphi}, \qquad \psi(0) = \bar{\psi}, \qquad (30) $$

and $\varphi$ and $\psi$ satisfy the following identities in the classical sense on $(0, \tau) \times \mathbb{T}^d$

$$ \partial_t \varphi + \frac{1}{2} |\nabla \varphi|^2 + \frac{\varepsilon}{2\tau} \Delta \varphi = 0 \quad \text{and} \quad \partial_t \psi + \frac{1}{2} |\nabla \psi|^2 - \frac{\varepsilon}{2\tau} \Delta \psi = 0. \qquad (31) $$

Formal proof. The arguments of the previous proof show that if $(\rho, v)$ is optimal in (6), $\varphi$ is an optimizer in (26) and $\bar{\varphi}$ is the corresponding optimizer of (22), then:

* $v = \nabla \varphi$,

* $\varphi$ is a solution of the first equation in (31),

* the auxiliary function $\theta := \exp(\frac{\tau}{\varepsilon}\varphi)$ solves $\partial_t \theta + \frac{\varepsilon}{2\tau}\Delta\theta = 0$,

* $\bar{\varphi} = \varphi(\tau)$.

Therefore, the first equations in (27), (28), (30) and (31) are proved. The second equations – concerning $\psi$ – are proved in the same way inverting the direction of time. Finally, to prove (29), observe that because of (21) and (28), we have

$$ \partial_t \rho + \text{div} \left( \rho \nabla \left( \frac{\varphi + \psi}{2} \right) \right) = 0 \quad \text{and hence} \quad \partial_t \log \rho + \nabla \log \rho \cdot \nabla \left( \frac{\varphi + \psi}{2} \right) + \Delta \left( \frac{\varphi + \psi}{2} \right) = 0. $$

Therefore, using this PDE and (31), a quick computation shows

$$ \begin{aligned} \partial_t \left( \varphi - \psi - \frac{\varepsilon}{\tau} \log \rho \right) &= -\frac{1}{2} \left\{ \left|\nabla \varphi\right|^2 - \left|\nabla \psi\right|^2 \right\} + \frac{\varepsilon}{\tau} \nabla \log \rho \cdot \nabla \left( \frac{\varphi + \psi}{2} \right) \\ &= -\nabla \left( \frac{\varphi + \psi}{2} \right) \cdot \nabla \left( \varphi - \psi - \frac{\varepsilon}{\tau} \log \rho \right) = 0, \end{aligned} $$

where we used that, by (21) and (28), we have $\nabla(\varphi - \psi - \frac{\varepsilon}{\tau} \log \rho) = 0$. We conclude by observing that at time $t = \tau$, $\varphi(\tau) - \psi(\tau) - \frac{\varepsilon}{\tau} \log \rho(\tau) \equiv 0$ as a consequence of (27), (30) and (24). $\square$

As a corollary of this proposition, several evolution PDEs can be derived for relevant quantities of the problem. Here is a list of some that will be useful for us.

**Corollary 3.8.** In the situation of Proposition 3.7, calling $\vec{m} := \rho\nabla\varphi$ and $\vec{m} := \rho\nabla\psi$ the so called forward and backward momentums, the following identities hold classicaly in $(0, \tau) \times \mathbb{T}^d$:

$$ \partial_t \vec{m} + \mathbf{div}(\vec{m} \otimes \nabla \psi) + \frac{\varepsilon}{2\tau} \Delta \vec{m} = 0, \qquad \partial_t \vec{m} + \mathbf{div}(\vec{m} \otimes \nabla \varphi) - \frac{\varepsilon}{2\tau} \Delta \vec{m} = 0 \qquad (32) $$
(here the notation **div** stand for the divergence of a matrix field, which provides a vector field as a result). Moreover, denoting by $\vec{e} = \rho|\nabla\varphi|^2/2$ and $\vec{e} = \rho|\nabla\psi|^2/2$ the so-called forward and backward local energies, we have classically in $(0, \tau) \times \mathbb{T}^d$

$$ \partial_t \vec{e} + \mathrm{div}(\vec{e}\nabla\psi) + \frac{\varepsilon}{2\tau}\Delta\vec{e} = \frac{1}{2}\rho|D^2\varphi|^2, \quad \partial_t \vec{e} + \mathrm{div}(\vec{e}\nabla\varphi) - \frac{\varepsilon}{2\tau}\Delta\vec{e} = -\frac{1}{2}\rho|D^2\psi|^2. \quad (33) $$

*Proof.* This proof is only algebraic. Recall that all the following PDEs hold in a classical sense on $(0, \tau) \times \mathbb{T}^d$, and that all the quantities involved are of class $C^\infty$ because of (27) and (29):

$$ \left\{ \begin{array}{l} \partial_t \rho + \mathrm{div}(\rho \nabla \varphi) - \frac{\varepsilon}{2\tau} \Delta \rho = 0, \\ \partial_t \varphi + \frac{1}{2} |\nabla \varphi|^2 + \frac{\varepsilon}{2\tau} \Delta \varphi = 0, \end{array} \right. \qquad \left\{ \begin{array}{l} \partial_t \rho + \mathrm{div}(\rho \nabla \psi) + \frac{\varepsilon}{2\tau} \Delta \rho = 0, \\ \partial_t \psi + \frac{1}{2} |\nabla \psi|^2 - \frac{\varepsilon}{2\tau} \Delta \psi = 0. \end{array} \right. $$

Therefore, we have

$$ \begin{aligned} \partial_t \vec{m} + \mathbf{div}(\vec{m} \otimes \nabla \psi) &= \partial_t(\rho \nabla \varphi) + \mathbf{div}(\nabla \varphi \otimes (\rho \nabla \psi)) \\ &= (\partial_t \rho + \mathrm{div}(\rho \nabla \psi)) \nabla \varphi + \rho (\nabla(\partial_t \varphi) + D^2 \varphi \cdot \nabla \psi) \\ &= -\frac{\varepsilon}{2\tau} \Delta \rho \nabla \varphi + \rho (\nabla(\partial_t \varphi) + D^2 \varphi \cdot \nabla \psi). \end{aligned} $$

Recalling that

$$ \nabla \psi = \nabla \varphi - \frac{\varepsilon}{\tau} \nabla \log \rho, $$

we deduce that

$$ \begin{aligned} \partial_t \vec{m} + \mathbf{div}(\vec{m} \otimes \nabla \psi) &= -\frac{\varepsilon}{2\tau} \Delta \rho \nabla \varphi + \rho \left( \nabla \left( \partial_t \varphi + \frac{1}{2} |\nabla \varphi|^2 \right) - \frac{\varepsilon}{\tau} D^2 \varphi \cdot \nabla \log \rho \right) \\ &= -\frac{\varepsilon}{2\tau} \left( \Delta \rho \nabla \varphi + 2 D^2 \varphi \cdot \nabla \rho + \rho \nabla \Delta \varphi \right) = -\frac{\varepsilon}{2\tau} \Delta \vec{m}, \end{aligned} $$

which is the first equation in (32). The second one is proved in the same way.

The same computations let us compute:

$$ \begin{aligned} \partial_t \vec{e} + \mathrm{div}(\vec{e}\nabla\psi) &= \frac{1}{2} \left( \partial_t(\vec{m} \cdot \nabla\varphi) + \mathrm{div}((\nabla\psi \otimes \vec{m}) \cdot \nabla\varphi) \right) \\ &= \frac{1}{2} \left( \partial_t \vec{m} + \mathbf{div}(\vec{m} \otimes \nabla\psi) \right) \cdot \nabla\varphi + \vec{m} \cdot \left( \nabla(\partial_t\varphi) + D^2\varphi \cdot \nabla\psi \right) \\ &= -\frac{\varepsilon}{4\tau} \Delta \vec{m} \cdot \nabla\varphi + \vec{m} \cdot \left( \nabla \left( \partial_t\varphi + \frac{1}{2}|\nabla\varphi|^2 \right) - \frac{\varepsilon}{\tau} D^2\varphi \cdot \nabla \log \rho \right) \\ &= -\frac{\varepsilon}{4\tau} \left( \Delta \vec{m} \cdot \nabla\varphi + 2\vec{m} \cdot (D^2\varphi \cdot \nabla \log \rho) + \vec{m} \cdot \nabla \Delta\varphi \right). \end{aligned} $$

On the other hand,

$$ \Delta \vec{e} = \frac{1}{2} \Delta \vec{m} \cdot \nabla \varphi + \mathrm{tr}(D^2 \varphi \cdot D \vec{m}) + \frac{1}{2} \vec{m} \cdot \nabla \Delta \varphi. $$

Therefore, also using that $\vec{m} = \rho \nabla \varphi$ and hence $D\vec{m} = \rho D^2 \varphi + \nabla \varphi \otimes \nabla \rho$, we find

$$ \begin{aligned} \partial_t \vec{e} + \mathrm{div}(\vec{e}\nabla\psi) + \frac{\varepsilon}{2\tau}\Delta\vec{e} &= \frac{\varepsilon}{2\tau}\left(\mathrm{tr}(D^2\varphi \cdot D\vec{m}) - \vec{m} \cdot (D^2\varphi \cdot \nabla \log \rho)\right) \\ &= \frac{\varepsilon}{2\tau}\mathrm{tr}\left(D^2\varphi \cdot \{D\vec{m} - \vec{m} \otimes \nabla \log \rho\}\right) \\ &= \frac{\varepsilon}{2\tau}\mathrm{tr}\left(D^2\varphi \cdot \{D\vec{m} - \nabla\varphi \otimes \nabla \rho\}\right) = \frac{\varepsilon}{2\tau}\rho|D^2\varphi|^2. \end{aligned} $$

This is the first equation in (33), and the second one follows the same lines.

<figure><img src="image_2.png" /></figure>
The final result of this section is a keystone when deriving the optimality conditions of one step of the regularized JKO scheme. It essentially asserts that the subdifferential of the convex functional $\nu \mapsto D_{\varepsilon}(\mu, \nu)^2$ is the set of maximizers in (22). We provide a detailed proof of this result.

**Proposition 3.9.** Let $\mu, \nu \in \mathcal{P}(\mathbb{T}^d)$ with $H(\mu)$ and $H(\nu)$ finite, and let $\varepsilon, \tau > 0$. Finally, let $\varphi$ be a maximizer of (22).

* If $\nu' \in \mathcal{P}(\mathbb{T}^d)$ is not absolutely continuous with respect to $\nu$ and $H(\nu') < +\infty$, then calling $\nu_u := (1-u)\nu + u\nu'$ for $u \in [0, 1]$, we have:

$$ \left. \frac{d}{du} \frac{D_{\varepsilon}(\mu, \nu_u)^2}{2\tau} \right|_{u=0} = -\infty. $$

* If $\zeta \in L^{\infty}(\nu)$ and $\int \zeta d\nu = 0$, then calling $\nu_u := (1+u\zeta)\nu$ for $u \in \mathbb{R}$ with $|u|$ small enough, we have:

$$ \left. \frac{d}{du} \frac{D_{\varepsilon}(\mu, \nu_u)^2}{2\tau} \right|_{u=0} = \int \varphi \zeta \nu. $$

*Proof.* First point. Let us call $(\rho, w)$ and $(\rho', w')$ the optimizers for problem (17) between $\mu$ and $\nu$, and $\mu$ and $\nu'$ respectively as given by Proposition 3.4. By Proposition 3.1 and Corollary 3.3, we easily see that for all $u \in [0, 1]$,

$$ \frac{D_{\varepsilon}(\mu, \nu_u)^2}{2\tau} \le \frac{\varepsilon}{\tau} H(\nu_u) + (1-u) \left\{ \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} - \frac{\varepsilon}{\tau} H(\nu) \right\} + u \left\{ \frac{D_{\varepsilon}(\mu, \nu')^2}{2\tau} - \frac{\varepsilon}{\tau} H(\nu') \right\}. $$

Rearranging the terms and neglecting terms with the right sign, we find for all $u \in (0, 1]$:

$$ \frac{D_{\varepsilon}(\mu, \nu_u)^2 - D_{\varepsilon}(\mu, \nu)^2}{2\tau u} \le \frac{\varepsilon}{\tau} \frac{H(\nu_u) - H(\nu)}{u} + \frac{\varepsilon}{\tau} H(\nu) + \frac{D_{\varepsilon}(\mu, \nu')^2}{2\tau}. $$

Therefore, we only need to show that

$$ \left. \frac{d}{du} H(\nu_u) \right|_{u=0} = -\infty. $$

This is an easy consequence of the monotone convergence theorem since

$$ u \mapsto \frac{\nu_u \log \nu_u - \nu \log \nu}{u} $$

is nonincreasing and converges to $-\infty$ on the set $\{\nu' > 0\} \cap \{\nu = 0\}$, which is of positive Lebesgue measure.

**Second point.** Take $u \in \mathbb{R}$ sufficiently small so that we have $(1+u\zeta)\nu \in \mathcal{P}(\mathbb{T}^d)$ (i.e. $1+u\zeta \ge 0$). By Proposition 3.5, we have

$$ \begin{aligned} \frac{D_{\varepsilon}(\mu, \nu_u)^2}{2\tau} &\ge \frac{\varepsilon}{\tau} H(\mu) + \langle \varphi, \nu_u \rangle - \left\langle \frac{\varepsilon}{\tau} \log \left( \left\{ \exp \left( \frac{\tau\varphi}{\varepsilon} \right) \right\} * \sigma_{\varepsilon} \right), \mu \right\rangle \\ &= \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} + \langle \varphi, \nu_u - \nu \rangle = \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} + u \langle \varphi, \zeta \nu \rangle. \end{aligned} $$
On the other hand, let us call $\varphi_u$ the optimizer in (22) with $\nu_u$ in place of $\nu$, chosen so that $\int \exp(\frac{\tau}{\varepsilon}\varphi_u) = 1$ (recall that it is unique up to an additive constant). First, let us prove that $\varphi_u$ converges to $\varphi$ as $u \to 0$ up to an additive constant. Indeed, the optimality condition (25) for $\varphi_u$ writes:

$$ \varphi_u = \frac{\varepsilon}{\tau} \log \nu_u - \frac{\varepsilon}{\tau} \log \left( \left\{ \frac{\mu}{\exp(\frac{\tau}{\varepsilon}\varphi_u) * \sigma_{\varepsilon}} \right\} * \sigma_{\varepsilon} \right). \quad (34) $$

Yet, $(\exp(\frac{\tau}{\varepsilon}\varphi_u))_u$ is precompact in the narrow topology of $\mathcal{P}(\mathbb{T}^d)$, so $(\exp(\frac{\tau}{\varepsilon}\varphi_u) * \sigma_{\varepsilon})_u$ is precompact in $C(\mathbb{T}^d)$ with a uniform below bound, and so the second term in the right-hand side of (34) is precompact in $C(\mathbb{T}^d)$. Also, the first term in the right hand side of (34) clearly converges to $\frac{\varepsilon}{\tau} \log \nu$ in $L^1(\mathbb{T}^d)$. In turn, $(\varphi_u)_u$ is precompact in $L^1(\mathbb{T}^d)$. Let $\bar{\varphi}$ be a corresponding limit point. Passing to the limit in (34), we deduce that $\bar{\varphi}$ satisfies the optimality condition (25) associated with the marginals $(\mu, \nu)$ and by uniqueness, $\varphi - \bar{\varphi}$ is constant.

Now, the same argument as before provides

$$ \frac{D_{\varepsilon}(\mu, \nu_u)^2}{2\tau} \le \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} + u\langle \varphi_u, \zeta\nu \rangle. $$

Using (34) makes it easy to show

$$ \langle \varphi_u, \zeta\nu \rangle \underset{u \to 0}{\longrightarrow} \langle \bar{\varphi}, \zeta\nu \rangle = \langle \varphi, \zeta\nu \rangle, $$

where the last equality follows from $\int \zeta\nu = 0$ and the fact that $\varphi - \bar{\varphi}$ is constant.

All in all,

$$ \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} + u\langle \varphi, \zeta\nu \rangle \le \frac{D_{\varepsilon}(\mu, \nu_u)^2}{2\tau} \le \frac{D_{\varepsilon}(\mu, \nu)^2}{2\tau} + u\langle \varphi, \zeta\nu \rangle + o(u), $$

and the result follows.

□

# 4 Proof of Theorem 2.9

Let $(\bar{\rho}^k)$ be a piecewise constant curve as in the statement of the theorem, and $T > 0$ is a finite time horizon given once for all. For the sake of simplicity, we will assume that for all $k \in \mathbb{N}$, $T$ is an integer multiple of $\tau_k$, but the proof is easy to adapt otherwise.

It will often be convenient to work with another curve than $\rho^k$ valued in $\mathcal{P}(\mathbb{T}^d)$, interpolating between the sequence defined through the JKO scheme. We define this interpolation as follows. Let $k, n \in \mathbb{N}$ and $(\mu, \xi) = (\mu(s), \xi(s))$, $s \in [0, \tau_k]$ be the optimal pair in (6) as given by Proposition 3.4 with $\varepsilon = \varepsilon_k$ and $\tau = \tau_k$, between $\rho_n^k$ and $\rho_{n+1}^k$. For all $s \in (0, \tau_k]$, we define

$$ \rho^k(n\tau_k + s) := \mu(s) \quad \text{and} \quad v^k(n\tau_k + s) := \xi(s). $$

By definition, it is clear that for every $k \in \mathbb{N}$ we have:

* By Proposition 3.4, the curve $\rho^k$ belongs to $AC^2([0, T]; \mathcal{P}(\mathbb{T}^d))$ and so in particular is continuous.
• Distributionally, on $(0, T) \times \mathbb{T}^d$, recalling that $\alpha_k = \varepsilon_k/\tau_k$,

$$ \partial_t \rho^k + \operatorname{div}(\rho^k v^k) = \frac{\alpha_k}{2} \Delta \rho^k. \qquad (35) $$

(This equation is obviously satisfied on sets of type $(n\tau_k, (n+1)\tau_k)$, $n \in \mathbb{N}$, and we deduce that it is satisfied in the whole $(0, T) \times \mathbb{T}^d$ using continuity of $\rho^k$.)

• For all $n \in \mathbb{N}$,

$$ \rho^k(n\tau_k) = \bar{\rho}^k(n\tau_k) = \rho_n^k. $$

From now on, the proof is divided into several steps following quite closely the proof for the classical JKO scheme (see [26, Section 8.3] or [27, Section 4.4]), with two really new ideas. In Steps 1 and 2, we derive and exploit estimates showing that $(\rho^k)_{k \in \mathbb{N}}$ is precompact in the uniform topology of $C([0, T]; \mathcal{P}(\mathbb{T}^d))$. This is made by comparing $\rho^k$ between two time steps with solutions of the heat equation, which is the first novelty of our proof. Then, in Steps 3 to 6, we prove the distributional convergence (15). The main arguments to this convergence are the derivation of the optimality conditions (46) for each step of the problem and a quantitative study of the continuity of the momentums $\rho^k v^k$ in some adapted topology. Finally, Steps 7 to 9 are devoted to proving that every term in (15) converges as $k \to 0$. To this aim, the only difficulty is to show the convergence

$$ g(\bar{\rho}^k) \xrightarrow{k \to +\infty} g(\rho). $$

This is done thanks to the estimate (51) that we will introduce later, which ensures strong compactness in space for $\rho^k$. This estimate is a consequence of some monotonicity property of the kinetic energy along Schrödinger bridges, the second main new ingredient of the proof. We conclude with a classical argument of the Aubin-Lions type that we decided to deduce directly from [24].

**Step 1.** A bound for $v^k$ in $L^2(\rho^k)$.

Let us prove that $(\rho^k)$ is precompact in the uniform topology of $C([0, T]; \mathcal{P}(\mathbb{T}^d))$. As in the non-regularized case, this compactness relies on the fact that the kinetic action along one time step is controlled by the dissipation of the energy functional during this time step. The purpose of this step is to derive this control. We will deduce compactness at Step 2 below.

Our bound is obtained by comparing $(\rho^k, v^k)$ to another competitor at each step, and hence it does not use the optimality conditions for the Schrödinger problem. Therefore, at this step, we will only use the Propositions 3.2 and 3.4.

In the classical case, compactness is obtained by considering stationary curves as competitors. Here, this is done by taking solutions of the heat equation. For given $k, n \in \mathbb{N}$ and $s \in (0, \tau_k]$, let us call $\mu(s) := \rho_n^k * \sigma_{\alpha_k s}$ and $\xi(s) := 0$. Comparing $\rho_n^k * \sigma_{\tau_k}$ to $\rho_{n+1}^k$ in the optimization problem (3) and then using $(\mu, \xi)$ as a competitor in (6) for computing $D_{\varepsilon_k}(\rho_n^k, \rho_n^k * \sigma_{\tau_k})$, we find

$$ \int_{n\tau_k}^{(n+1)\tau_k} \int \frac{|v^k|^2}{2} \rho^k ds + \mathcal{F}(\rho_{n+1}^k) \le \mathcal{F}(\rho_n^k * \sigma_{\tau_k}). \qquad (36) $$

(Observe that the initial entropic term from (6) is the same on both sides of this inequality, and hence simplifies.)
On the other hand, for all $\rho \in \mathcal{P}(\mathbb{T}^d)$ with $\mathcal{F}(\rho) < +\infty$ and $t \ge 0$, we have

$$ \begin{aligned} \mathcal{F}(\rho * \sigma_t) - \mathcal{F}(\rho) &\le \int \{V * \sigma_t - V\} \rho + \int \{(W * \sigma_t - W) * \rho\} \rho * \sigma_t + \int \{f(\rho * \sigma_t) - f(\rho)\} \mathrm{d}x \\ &\le \sup\{V * \sigma_t - V\} + \sup\{(W * \sigma_t - W)\} \end{aligned} $$

(where we used Jensen's inequality to bound the last term of the first line.) But now, it is clear that for all $C^2$ function $A$ on $\mathbb{T}^d$ (actually $(\Delta A)_+ \in L^\infty$ is enough), there exists $C > 0$ such that

$$ \sup\{A * \sigma_t - A\} \le C_t, \quad \text{for all } t \ge 0. $$

Plugging this observation into (36), we find $C$ only depending on $V$ and $W$ such that

$$ \int_{n\tau_k}^{(n+1)\tau_k} \int \frac{|v^k|^2}{2} \rho^k \mathrm{d}s \le \mathcal{F}(\rho_n^k) - \mathcal{F}(\rho_{n+1}^k) + C\tau_k. \quad (37) $$

This is the main estimate needed the proof of compactness of $(\rho^k)$.

Before giving this proof of compactness, let us provide a few direct consequences of (37) that will be useful in the sequel. Following Proposition 3.4, let us introduce

$$ c^k := v^k - \frac{\alpha_k}{2} \nabla \log \rho^k \quad \text{and} \quad w^k := v^k - \alpha_k \nabla \log \rho^k. $$

Identity (20) applies in this context, and plugging it in (37) for each of these vector fields, we get

$$ \begin{aligned} \frac{1}{2} \int_{n\tau_k}^{(n+1)\tau_k} \int \left\{ |c^k|^2 + \left| \frac{\alpha_k}{2} \nabla \log \rho^k \right|^2 \right\} \rho^k \mathrm{d}s \le \mathcal{F}(\rho_n^k) - \mathcal{F}(\rho_{n+1}^k) + \frac{\alpha_k}{2} \left( H(\rho_n^k) - H(\rho_{n+1}^k) \right) + C\tau_k, \\ \int_{n\tau_k}^{(n+1)\tau_k} \int \frac{|w^k(s)|^2}{2} \rho^k(s) \mathrm{d}s \le \mathcal{F}(\rho_n^k) - \mathcal{F}(\rho_{n+1}^k) + \alpha_k \left( H(\rho_n^k) - H(\rho_{n+1}^k) \right) + C\tau_k. \end{aligned} $$

Summing such inequalities over values of $n$, we get

$$ \mathcal{F}(\rho^k(T)) + \int_0^T \int \frac{|v^k|^2}{2} \rho^k \mathrm{d}t \le \mathcal{F}(\rho_0) + CT, \quad (38) $$

$$ \mathcal{F}(\rho^k(T)) + \frac{\alpha_k}{2} H(\rho^k(T)) + \frac{1}{2} \int_0^T \int \left[ |c^k|^2 + \left| \frac{\alpha_k}{2} \nabla \log \rho^k \right|^2 \right] \rho^k \mathrm{d}t \le \mathcal{F}(\rho_0) + \frac{\alpha_k}{2} H(\rho_0) + CT, \quad (39) $$

$$ \mathcal{F}(\rho^k(T)) + \alpha_k H(\rho^k(T)) + \int_0^T \int \frac{|w^k|^2}{2} \rho^k \mathrm{d}t \le \mathcal{F}(\rho_0) + \alpha_k H(\rho_0) + CT. \quad (40) $$

**Step 2.** $(\rho^k)$ is compact in the uniform topology of $C([0, T]; \mathcal{P}(\mathbb{T}^d))$.

Recalling the link between the metric derivative of $\rho^k$ and the vector field $c^k$ given by Proposition 3.4, estimate (39) provides:

$$ \frac{1}{2} \int_0^T |\dot{\rho}^k(s)|^2 \mathrm{d}s \le \mathcal{F}(\rho_0) + \frac{\alpha_k}{2} H(\rho_0) + CT. \quad (41) $$

We conclude using the usual arguments that:
* The sequence $(\rho^k)$ is uniformly bounded in $AC^2([0, T]; \mathcal{P}(\mathbb{T}^d))$, and hence relatively compact in $C([0, T]; \mathcal{P}(\mathbb{T}^d))$.

* The distance $D(\bar{\rho}^k(t), \rho^k(t)) = D(\rho^k(\lceil t/\tau_k \rceil \tau_k), \rho^k(t))$ vanishes as $k \to +\infty$, uniformly in $[0, T]$.

* Therefore, $(\bar{\rho}^k)$ is compact in $L^\infty([0, T]; \mathcal{P}(\mathbb{T}^d))$ and has the same limit points as $(\rho^k)$.

From now on, up to extraction, we can assume that there is a limiting curve in $C(\mathbb{R}_+; \mathcal{P}(\mathbb{T}^d))$ such that both $\bar{\rho}^k$ and $\rho^k$ converge to $\rho$ in $\mathcal{P}(\mathbb{T}^d)$, locally uniformly in time.

**Step 3.** Distributional convergence, first observations.

Here, we prove that $g(\bar{\rho}^k)$ belongs to $L^1([0, T] \times \mathbb{T}^d)$ and that the convergence (15) holds.

At Step 4 and 5, we will find that for all $k \in \mathbb{N}$, $n \in \mathbb{N}^*$, $\rho_n^k$ is of class $C^2$ and with values in $(d_-, d_+)$, the interior of $\text{Dom}(f)$. Therefore, we will get for free that $g(\bar{\rho}^k) \in L^1([0, T] \times \mathbb{T}^d)$, and even that for all $t > 0$, in the classical sense, using (11):

$$ \nabla g(\bar{\rho}^k(t)) = g'(\bar{\rho}^k(t)) \nabla \bar{\rho}^k(t) = \bar{\rho}^k(t) f''(\bar{\rho}^k(t)) \nabla \bar{\rho}^k(t) = \bar{\rho}^k(t) \nabla f'(\bar{\rho}^k(t)). \quad (42) $$

So let us temporary admit this regularity result and gather some information necessary to prove (15).

The very first observation is that $\rho^k$ solves (35), and as $\bar{\rho}^k - \rho^k$ converge to 0 in $L^\infty([0, T]; \mathcal{P}(\mathbb{T}^d))$, this convergence also holds distributionally. Therefore, we can replace $\rho^k$ by $\bar{\rho}^k$ in the $\partial_t$ and $\Delta$ terms of (35) and deduce

$$ \partial_t \bar{\rho}^k + \text{div}(\rho^k v^k) - \frac{\alpha_k}{2} \Delta \bar{\rho}^k \underset{k \to +\infty}{\longrightarrow} 0. $$

Therefore, using (42), we conclude that proving (15) reduces to prove that distributionally,

$$ \rho^k v^k + \bar{\rho}^k \nabla (V + W * \bar{\rho}^k + f'(\bar{\rho}^k)) \underset{k \to +\infty}{\longrightarrow} 0. \quad (43) $$

This is where the optimality conditions of our optimization problem comes into play. At Step 4 and 5, we will see that this condition writes as follows. Up to considering a left continuous version of $v^k$, for all $n \in \mathbb{N}^*$, we have at time $n\tau_k$:

$$ v^k = -\nabla (V + W * \bar{\rho}^k + f'(\bar{\rho}^k)). \quad (44) $$

Therefore, at that times, convergence (43) is actually an equality.

Moreover, let us call

$$ \vec{m}^k := \rho^k v^k. \quad (45) $$

At Step 6, we will show that for a given $n \in \mathbb{N}$ and $t \in (n\tau_k, (n+1)\tau_k]$, the quantity in (43) – which as a consequence of (44) is nothing but $\vec{m}^k(t) - \vec{m}^k((n+1)\tau_k)$ – can be controlled using (32).

**Step 4.** Optimality conditions.

The purpose of this step is to show that for all $k, n \in \mathbb{N}$:

1. The measure $\rho_{n+1}^k$ has Leb-almost everywhere values in $(d_-, d_+)$, the interior of $\text{Dom}(f)$.

2. If $\varphi$ is a maximizer for the dual problem (22) with $\mu = \rho_n^k$, $\nu = \rho_{n+1}^k$, $\varepsilon = \varepsilon_k$ and $\tau = \tau_k$, the function

$$ \zeta := \varphi + V + W * \rho_{n+1}^k + f'(\rho_{n+1}^k), \quad (46) $$

which is well defined almost everywhere thanks to the first point, is constant.
The consequences of these two points will be given in the next step. Our proof relies on the inequality

$$ \frac{D_{\varepsilon_k}(\rho_n^k, \rho_{n+1}^k)^2}{2\tau_k} + \mathcal{F}(\rho_{n+1}^k) \le \frac{D_{\varepsilon_k}(\rho_n^k, \rho)^2}{2\tau_k} + \mathcal{F}(\rho) \quad (47) $$

written for well chosen competitors $\rho \in \mathcal{P}(\mathbb{T}^d)$, and computing derivatives as in Proposition 3.9.

Let us first prove that Leb-almost everywhere, $\rho_{n+1}^k \in (d_-, d_+)$. We know that

$$ \int f(\rho_{n+1}^k) < +\infty. $$

Therefore, almost everywhere, $\rho_{n+1}^k \in \text{Dom}(f)$. We now need to show that under Assumption 2.5, $\rho_{n+1}^k$ cannot touch the boundaries of $\text{Dom}(f)$ on a set of positive Lebesgue measure.

To do so, for $u \in [0, 1]$, we define $\rho_u := (1-u)\rho_{n+1}^k + u \text{Leb}$. Applying (47) to $\rho_u$, and using the convexity of $\rho \mapsto \int f(\rho)$, we get for all $u \in [0, 1]$

$$ \frac{D_{\varepsilon_k}(\rho_n^k, \rho_{n+1}^k)^2}{2\tau_k} + \int \left(V + \frac{1}{2}W * \rho_{n+1}^k\right) \rho_{n+1}^k \le \frac{D_{\varepsilon_k}(\rho_n^k, \rho_u)^2}{2\tau_k} + \int \left(V + \frac{1}{2}W * \rho_u\right) \rho_u + u \int \left\{f(1) - f(\rho_{n+1}^k)\right\}. $$

Rearranging the terms, we find for all $u \in (0, 1]$:

$$ \int \{f(\rho_{n+1}^k) - f(1)\} + \int \left(V + W * \rho_{n+1}^k\right) (\rho_{k+1}^n - 1) + \mathcal{o}(1) \le \frac{1}{u} \left( \frac{D_{\varepsilon_k}(\rho_n^k, \rho_u)^2}{2\tau_k} - \frac{D_{\varepsilon_k}(\rho_n^k, \rho_{n+1}^k)^2}{2\tau_k} \right). $$

So the right hand side cannot converge to $-\infty$ as $u \to 0$, and because of Proposition 3.9, we need to have $\text{Leb} \ll \rho_{n+1}^k$, and hence $\rho_{n+1}^k > 0$ almost everywhere. We obtain $\rho_{n+1}^k > d_-$ almost everywhere in the particular case where $d_- = 0$.

Now, we apply once again (47) to $\rho_u$, but this time using the convexity of $\rho \mapsto D_{\varepsilon_k}(\rho_n^k, \rho)^2$. Rearranging the terms in a similar way as what we just did, we find for all $u \in (0, 1]$:

$$ \frac{D_{\varepsilon_k}(\rho_n^k, \rho_{n+1}^k)^2 - D_{\varepsilon_k}(\rho_n^k, \text{Leb})^2}{2\tau_k} + \int \left(V + W * \rho_{n+1}^k\right) (\rho_{k+1}^n - 1) + \mathcal{o}(1) \le \int \frac{f((1-u)\rho_{n+1}^k + u) - f(\rho_{n+1}^k)}{u}. $$

We conclude by monotone convergence that

$$ \int f'(\rho_{n+1}^k)(1 - \rho_{n+1}^k) > -\infty, $$

and with (13), we conclude that almost everywhere, $\rho_{n+1}^k > d_-$ also when $d_- > 0$, and $\rho_{n+1}^k < d_+$, even when $d_+ < +\infty$.

Let us prove the second point of our claim. For all $\eta > 0$ sufficiently small, let us call

$$ A_{\eta} := \{ x \in \mathbb{T}^d \text{ such that } \rho_{n+1}^k(x) \in (d_- + \eta, \min(d_+ - \eta, 1/\eta)] \}, $$

and let us call $A$ the increasing limit of $A_\eta$, which is of full Lebesgue measure. On $A$, formula (46) defines $\zeta$ properly. For a given $\eta > 0$ sufficiently small to have $\rho_{n+1}^k(A_\eta) > 0$, let us call

$$ \zeta_{\eta} := \left( \zeta - \frac{1}{\rho_{n+1}^k(A_{\eta})} \int_{A_{\eta}} \zeta \rho_{n+1}^k \right) 1_{A_{\eta}}. $$
It is clearly a bounded function (recall that $\varphi - \alpha_k \log \rho_{n+1}^k$ is smooth in virtue of Proposition 3.5) with $\int \zeta_\eta \rho_{n+1}^k = 0$. Hence, for $u \in \mathbb{R}$ with $|u|$ sufficiently small, $\rho_u := (1 + u\zeta_\eta)\rho_{n+1}^k \in \mathcal{P}(\mathbb{T}^d)$, so (47) applies. Moreover, by Proposition 3.9, we have

$$ \left. \frac{d}{du} \frac{D_\varepsilon(\rho_n^k, \rho_u)^2}{2\tau} \right|_{u=0} = \int \varphi \zeta_\eta \rho_{n+1}^k. $$

On the other hand, the following derivative is straightforward to compute using the dominated convergence theorem:

$$ \left. \frac{d}{du} \mathcal{F}(\rho_u) \right|_{u=0} = \int \left( V + W * \rho_{n+1}^k + f'(\rho_{n+1}^k) \right) \zeta_\eta \rho_{n+1}^k. $$

Therefore, for all $\eta > 0$ small enough, we have

$$ \int \left( \varphi + V + W * \rho_{n+1}^k + f'(\rho_{n+1}^k) \right) \zeta_\eta \rho_{n+1}^k = 0. $$

This can be re-written as $\int \zeta_\eta \rho_{n+1}^k = 0$ and, using the fact that $\zeta_\eta$ has zero mean, we can subtract a constant from $\zeta$. Yet, $\zeta_\eta$ vanishes outside $A_\eta$, and on this set it differs from $\zeta$ by a constant: this means that we have $\int \zeta_\eta^2 \rho_{n+1}^k = 0$, i.e. $\zeta_\eta = 0$. Hence, for all $\eta > 0$ sufficiently small, $\zeta$ is constant on $A_\eta$. As $A$ is of full Lebesgue measure on the torus because of the first point of this step, we get the result. From now on, we set

$$ \varphi_{n+1}^k := - \left( V + W * \rho_{n+1}^k + f'(\rho_{n+1}^k) \right). \quad (48) $$

Such a function differs by a constant from any maximizer $\varphi$ of (22), therefore, by Proposition 3.5, choosing $\varphi$ such that $\varphi - \alpha_k \log \rho_{n+1}^k$ is of class $C^\infty$, the same holds for $\varphi_{n+1}^k$.

**Step 5.** Consequences of (46): regularity of $\rho_n^k$, $n \in \mathbb{N}^*$ and identity (44).

A first consequence of the optimality condition (46) is that for all $k \in \mathbb{N}$, $n \in \mathbb{N}^*$, $\rho_n^k$ and $\varphi_n^k$ are of class $C^2$ and $\rho_n^k$ takes values in a compact subset of $(d_-, d_+)$. Indeed, we already saw that $\varphi_n^k - \alpha_k \log \rho_n^k$ is of class $C^\infty$, and by assumptions, both $V$ and $W * \rho_{n+1}^k$ are of class $C^2$. Therefore,

$$ f'(\rho_{n+1}^k) + \alpha_k \log \rho_{n+1}^k $$

is of class $C^2$. In particular it is bounded, and using (13) we obtain that $\rho_{n+1}^k$ stays away from both $d_-$ and $d_+$. Moreover, since $f'$ is $C^2$ and nondecreasing on $(d_-, d_+)$, the function $f' + \alpha_k \log$ is a $C^2$ diffeomorphism of $(d_-, d_+)$ into its image, which allows to invert it. Our claim follows.

Following Proposition 3.7, let us define for all $n \in \mathbb{N}$ and $s \in (0, \tau]$:

$$ \varphi^k(n\tau + s) := \alpha_k \log \left( \exp \left( \frac{\varphi_{n+1}^k}{\alpha_k} \right) * \sigma_{\alpha_k(\tau-s)} \right). \quad (49) $$

As $\varphi_{n+1}^k$ is of class $C^2$, $\varphi^k$ is $C^1$ with respect to time and $C^2$ with respect to space on all intervals of type $(n\tau_k, (n+1)\tau_k]$. Moreover, up to changing $v^k$ on a negligible set of $\mathbb{R}_+ \times \mathbb{T}^d$, we have everywhere for all $k, n \in \mathbb{N}$ and $t \in (n\tau_k, (n+1)\tau_k]$:

$$ v^k(t) = \nabla \varphi^k(t). $$
At time $t = (n + 1)\tau_k$, we find (44).

**Step 6.** Proof of (43).

The last step in order to prove the distributional convergence (15) is to show (43). Let us define

$$ \vec{m}^k := \rho^k \nabla \varphi^k. $$

This is the same definition as in (45). By Step 5, it is easy to see that $\vec{m}^k$ is continuous in time and space on all intervals of type $(n\tau_k, (n + 1)\tau_k]$, and by Corollary 3.8, calling

$$ \psi^k := \varphi^k - \alpha_k \log \rho^k, \qquad (50) $$

we have in the classical sense on all intervals of type $(n\tau_k, (n + 1)\tau_k)$:

$$ \partial_t \vec{m}^k + \mathbf{div}(\vec{m}^k \otimes \nabla \psi^k) + \frac{\alpha_k}{2} \Delta \vec{m}^k = 0. $$

Therefore, for all $t > 0$, at least in a distributional sense,

$$ \begin{aligned} \rho^k(t)v^k(t) + \bar{\rho}^k(t)\nabla(V + W * \bar{\rho}^k(t) + f'(\bar{\rho}^k(t))) &= \vec{m}^k(t) - \vec{m}^k(\lceil t/\tau_k \rceil \tau_k) = -\int_t^{\lceil t/\tau_k \rceil} \partial_t \vec{m}^k \\ &= \mathbf{div} \left( \int_t^{\lceil t/\tau_k \rceil \tau_k} \rho^k(s) \nabla \varphi^k(s) \otimes \nabla \psi^k(s) \, ds \right) + \frac{\alpha_k}{2} \Delta \left( \int_t^{\lceil t/\tau_k \rceil \tau_k} \rho^k(s) \nabla \varphi^k(s) \, ds \right). \end{aligned} $$

Therefore, to get the convergence we want, it is enough to prove that both

$$ A(t, x) := \int_t^{\lceil t/\tau_k \rceil \tau_k} \rho^k(s, x) \nabla \varphi^k(s, x) \otimes \nabla \psi^k(s, x) \, ds \quad \text{and} \quad B(t, x) := \int_t^{\lceil t/\tau_k \rceil \tau_k} \rho^k(s, x) \nabla \varphi^k(s, x) \, ds $$

converge to 0 as $k \to +\infty$ in $L^1([0, T] \times \mathbb{T}^d)$. We only treat the case of $A$, because the case of $B$ is similar but easier. We have

$$ \begin{aligned} \|A\|_{L^1([0,T] \times \mathbb{T}^d)} &\le \int_0^T \int_{\lceil t/\tau_k \rceil \tau_k}^{\lceil t/\tau_k \rceil \tau_k} \int |\nabla \varphi^k(s)| |\nabla \psi^k(s)| \rho^k(s) \, dx \, ds \, dt \\ &\le \tau_k \int_0^T \int |\nabla \varphi^k(t)| |\nabla \psi^k(t)| \rho^k(t) \, dx \, dt \\ &\le \tau_k \sqrt{\int_0^T \int |\nabla \varphi^k(t)|^2 \rho^k(t) \, dx \, dt} \sqrt{\int_0^T \int |\nabla \psi^k(t)|^2 \rho^k(t) \, dx \, dt} \\ &\le \tau_k \sqrt{\mathcal{F}(\rho_0) + CT} \sqrt{\mathcal{F}(\rho_0) + \alpha_k H(\rho_0) + CT} \underset{k \to +\infty}{\longrightarrow} 0. \end{aligned} $$

In the last line, we used the bounds (38) and (40), as well as $w^k = \nabla \psi^k$ which is a consequence of Proposition 3.7.

**Step 7.** Strong compactness: first observation.

By Step 2, $\bar{\rho}^k \to \rho$ in the uniform topology of $C([0, T]; \mathcal{P}(\mathbb{T}^d))$. This implies the following convergence in the sense of distributions, easy to justify term by term:

$$ \partial_t \bar{\rho}^k - \mathrm{div}(\bar{\rho}^k(\nabla V + \nabla W * \bar{\rho}^k)) - \frac{\alpha_k}{2} \Delta \bar{\rho}^k \underset{k \to +\infty}{\longrightarrow} \partial_t \rho - \mathrm{div}(\rho(\nabla V + \nabla W * \rho)) - \frac{\alpha}{2} \Delta \rho. $$
Note that because of (38), for all $t \in \mathbb{R}_+$, $\mathcal{F}(\bar{\rho}^k(t))$ is bounded uniformly in $k$. Using the lower semi-continuity of $\rho \mapsto \int f(\rho) dx$ and the bounds on $V$ and $W$ and letting $k \to +\infty$, we deduce that for all $t \in \mathbb{R}_+$, $\rho(t)$ is absolutely continuous with respect to the Lebesgue measure, has values in $[d_-, d_+]$, and hence $g(\rho)$ is well defined almost everywhere, with possibly infinite values, see Remark 2.6.

Given (15), in order to justify that $\rho$ solves (16), we need to prove:

1. The function $g(\rho)$ belongs to $L^1([0, T] \times \mathbb{T}^d)$.

2. The convergence

$$g(\bar{\rho}^k) \underset{k \to +\infty}{\longrightarrow} g(\rho)$$

holds in the sense of distributions (we will actually show strong convergence in $L^1([0, T] \times \mathbb{T}^d)$).

**Step 8.** Strong compactness: main estimate.

The main tool towards this goal is the following estimate: There exists $C > 0$ only depending on $V$ and $W$ such that for all $k \in \mathbb{N}$,

$$\int_0^T \int \{|\nabla h(\bar{\rho}^k)|^2 + \alpha_k^2 |\nabla \sqrt{\bar{\rho}^k}|^2\} dx dt \le C (\mathcal{F}(\rho_0) + \alpha_k H(\rho_0) + T), \quad (51)$$

where for all $s \in (d_-, d_+)$,

$$h(s) := \int_1^s \sqrt{r} f''(r) \, dr. \quad (52)$$

(Recall that by Assumption 2.5, we have $d_- < 1 < d_+$, and that by Step 4 and the regularity deduced at Step 5, $\bar{\rho}^k$ only takes values in $(d_-, d_+)$.)

To prove this estimate, the starting point is (40): what we will do is to bound from below the left-hand side of (40), using the decrease of the backward kinetic energy deduced from Corollary 3.8. Indeed, thanks to (33) and the regularity proved at Step 5, calling $\tilde{e}^k := \frac{1}{2}\rho^k |\nabla \psi^k|^2 = \frac{1}{2}\rho^k |w^k|^2$, we see that for all $n \in \mathbb{N}$,

$$\tau_k \int \frac{|\nabla \psi^k((n+1)\tau_k)|^2}{2} \rho^k((n+1)\tau_k) \le \int_{n\tau_k}^{(n+1)\tau_k} \int \frac{|w^k|^2}{2} \rho^k dt.$$

On the other hand, we know by (48), (49) and (50) that

$$\psi^k((n+1)\tau_k) = -(V + W * \rho_{n+1}^k + f'(\rho_{n+1}^k) + \alpha_k \log \rho_{n+1}^k) + C.$$

Plugging this identity in the previous estimate, we find

$$\tau_k \int \frac{|\nabla (V + W * \rho_{n+1}^k + f'(\rho_{n+1}^k) + \alpha_k \log \rho_{n+1}^k)|^2}{2} \rho_{n+1}^k \le \int_{n\tau_k}^{(n+1)\tau_k} \int \frac{|w^k|^2}{2} \rho^k dt,$$

and finally, because $\bar{\rho}^k$ is constant on the time step intervals

$$\int_{n\tau_k}^{(n+1)\tau_k} \int \frac{|\nabla (V + W * \bar{\rho}^k + f'(\bar{\rho}^k) + \alpha_k \log \bar{\rho}^k)|^2}{2} \bar{\rho}^k dt \le \int_{n\tau_k}^{(n+1)\tau_k} \int \frac{|w^k|^2}{2} \rho^k dt,$$
Summing this inequality over $n$ and using (40), we get:

$$ \int_0^T \int \frac{|\nabla (V + W * \bar{\rho}^k + f'(\bar{\rho}^k) + \alpha_k \log \bar{\rho}^k)|^2}{2} \bar{\rho}^k dt \le \mathcal{F}(\rho_0) + \alpha_k H(\rho_0) + CT. $$

Now, using the inequality $|a|^2 \le 2|a + b + c|^2 + 4(|b|^2 + |c|^2)$, we find

$$ \int_0^T \int |\nabla (f'(\bar{\rho}^k) + \alpha_k \log \bar{\rho}^k)|^2 \bar{\rho}^k dt \\ \le 4(\mathcal{F}(\rho_0) + \alpha_k H(\rho_0) + CT) + 4 \int_0^T \int \{|\nabla V|^2 + |\nabla W * \bar{\rho}^k|^2\} \bar{\rho}^k dt. $$

To treat the left-hand side, observe that by the regularity obtained at Step 5 and the definition (52) of $h$, we have

$$ \begin{aligned} |\nabla(f'(\bar{\rho}^k) + \alpha_k \log \bar{\rho}^k)|^2 \bar{\rho}^k &= \left| \left( \sqrt{\bar{\rho}^k} f''(\bar{\rho}^k) + \frac{\alpha_k}{\sqrt{\bar{\rho}^k}} \right) \nabla \bar{\rho}^k \right|^2 \\ &= \left( \sqrt{\bar{\rho}^k} f''(\bar{\rho}^k) + \frac{\alpha_k}{\sqrt{\bar{\rho}^k}} \right)^2 |\nabla \bar{\rho}^k|^2 \\ &\ge \left( \sqrt{\bar{\rho}^k} f''(\bar{\rho}^k) \right)^2 |\nabla \bar{\rho}^k|^2 + \left( \frac{\alpha_k}{\sqrt{\bar{\rho}^k}} \right)^2 |\nabla \bar{\rho}^k|^2 \\ &= |h'(\bar{\rho}^k) \nabla \bar{\rho}^k|^2 + \left| \frac{\alpha_k \nabla \bar{\rho}^k}{\sqrt{\bar{\rho}^k}} \right|^2 \\ &= |\nabla h(\bar{\rho}^k)|^2 + 2\alpha_k^2 |\nabla \sqrt{\bar{\rho}^k}|^2, \end{aligned} $$

where the inequality follows from the convexity of $f$. For the right-hand side, we get our estimate by using uniform bounds for $\nabla V$ and $\nabla W$.

**Step 9.** A technical lemma of the Poincaré type.

Estimate (51) gives a control of $\nabla h(\bar{\rho}^k)$ in $L^2([0, T] \times \mathbb{T}^d)$. A crucial point in the coming steps is to deduce a control for $h(\bar{\rho}^k)$ in some Lebesgue space. We will use the following lemma.

**Lemma 4.1.** 1. Let $\rho \in \mathcal{P}(\mathbb{T}^d) \cap L^1(\mathbb{T}^d)$ with $h(\rho) \in H^1(\mathbb{T}^d)$. Let $p := 2d/(d-2)$. There exists $C$ only depending on $d$ and $h$ such that:

$$ \|h(\rho)\|_{L^p(\mathbb{T}^d)} \le C(1 + \|\nabla h(\rho)\|_{L^2(\mathbb{T}^d)}). $$

2. Let $A \subset L^1(\mathbb{T}^d) \cap \mathcal{P}(\mathbb{T}^d)$ be a family of functions such that there exists $s_- \in (d_-, 1)$ and $\delta_- \in (0, 1)$ such that for all $\rho \in A$

$$ \mathrm{Leb}(\rho \ge s_-) \ge \delta_- . $$

Then, there is $C$ only depending on $d$, $h$, $s_-$ and $\delta_-$ such that for all $\rho \in A$ with $h(\rho) \in H^1(\mathbb{T}^d)$:

$$ \|h(\rho)\|_{L^p(\mathbb{T}^d)} \le C(1 + \|\nabla h(\rho)\|_{L^2(\mathbb{T}^d)}). $$
3. Finally, any uniformly equi-integrable family $\mathcal{A} \subset L^1(\mathbb{T}^d) \cap \mathcal{P}(\mathbb{T}^d)$ satisfies the condition of the second point.

Proof of Lemma 4.1. Of course, by Sobolev embeddings, it is enough to estimate the $L^2$ norms instead of the $L^p$ norms of $(h(\rho))_{\pm}$. The main tool is a small adaptation of [10, Exercise 15, Section 5.10]: for all $\delta \in (0, 1)$, there exists a constant $C = C(\delta) > 0$ such that for all $f \in H^1(\mathbb{T}^d)$, if $\text{Leb}(f = 0) \ge \delta$, we have

$$ \|f\|_{L^2(\mathbb{T}^d)} \le C(\delta) \|\nabla f\|_{L^2(\mathbb{T}^d)}. \qquad (53) $$

An important observation for the proof of this lemma is that $s \mapsto h(s)$ is nonnegative when $s \ge 1$ and nonpositive when $s \le 1$.

Let us consider $\rho \in \mathcal{P}(\mathbb{T}^d) \cap L^1(\mathbb{T}^d)$ with $h(\rho) \in H^1(\mathbb{T}^d)$ and prove the first point of the lemma. Let $s_+ \in (1, d_+)$ and $\delta_+ := 1 - 1/s_+$. We have

$$ \text{Leb}(\rho \ge s_+) \le \frac{1}{s_+}, $$

so that

$$ \text{Leb}(h(\rho) \le h(s_+)) \ge \delta_+. $$

We conclude with (53) that

$$ \|(h(\rho) - h(s_+))_+\|_{L^2(\mathbb{T}^d)} \le C(\delta_+) \|\nabla h(\rho)\|_{L^2(\mathbb{T}^d)}, $$

and hence that

$$ \|h(\rho)_+\|_{L^2} \le h(s_+) + C(\delta_+) \|\nabla h(\rho)\|_{L^2(\mathbb{T}^d)}. $$

The second point follows the same lines. Let us assume in addition that $\rho \in \mathcal{A}$ and let us choose $s_- \in (d_-, 1)$ and $\delta_- \in (0, 1)$ as in the statement of the lemma. We find that $\text{Leb}(\rho \ge s_-) \ge \delta_-$. We conclude as before that

$$ \|h(\rho)\|_{L^2} \le |h(s_-)| + C(\delta_-) \|\nabla h(\rho)\|_{L^2(\mathbb{T}^d)}, $$

and the conclusion follows.

Let us finally prove the third point. Let $\mathcal{A} \subset L^1(\mathbb{T}^d) \cap \mathcal{P}(\mathbb{T}^d)$ be a uniformly equi-integrable set. Let $s_- \in (d_-, 1)$. Let $\delta_- = \delta_-(\mathcal{A}) > 0$ be such that for all $\rho \in \mathcal{A}$, for all Borelian set $A$, if $\text{Leb}(A) < \delta_-$ then,

$$ \int_A \rho < 1 - s_- $$

For all $\rho \in \mathcal{A}$,

$$ 1 = \int \rho \le s_- + \int_{\{\rho \ge s_\-\}} \rho. $$

So

$$ \int_{\{\rho \ge s_\-\}} \rho \ge 1 - s_\-, $$

and hence $\text{Leb}(\rho \ge s_-) \ge \delta_-$ by definition of $\delta_-$. $\square$

<figure><img src="image_3.png" /></figure>
Step 10. Strong compactness of $(\bar{\rho}^k)$ in $L^1([0, T] \times \mathbb{T}^d)$.

The purpose of this step is to use the estimate (51) to show that:

$$ \bar{\rho}^k \underset{k \to +\infty}{\longrightarrow} \rho \quad \text{in the strong topology of } L^1([0, T] \times \mathbb{T}^d). \qquad (54) $$

We will use extensively the results of [24]. The idea is to apply [24, Theorem 2] to show that $(\bar{\rho}^k)$ is precompact in the topology of convergence in measure of $M(0, T; L^1(\mathbb{T}^d))$ (see [24, Section 1.1]). As $(\bar{\rho}^k)$ is clearly uniformly equi-integrable in $L^1(0, T; L^1(\mathbb{T}^d))$ (simply because for all $k \in \mathbb{N}$ and all $t \in \mathbb{R}_+$, $\|\bar{\rho}^k(t)\|_{L^1(\mathbb{T}^d)} = 1$), by [24, Proposition 1.10], our claim (54) will follow.

We define a functional $\mathfrak{F}$ which can take two different forms according to $\alpha > 0$ or $\alpha = 0$. In the case $\alpha > 0$, for all $\mu \in L^1(\mathbb{T}^d)$, let us set

$$ \mathfrak{F}(\mu) := \begin{cases} \int |\nabla \sqrt{\mu}|^2 dx & \text{if } \mu \in \mathcal{P}(\mathbb{T}^d), \text{ and } \sqrt{\mu} \in H^1(\mathbb{T}^d), \\ + \infty & \text{otherwise,} \end{cases} $$

while in the case $\alpha = 0$ we set

$$ \mathfrak{F}(\mu) := \begin{cases} \int |\nabla h(\mu)|^2 dx & \text{if } \mu \in \mathcal{P}(\mathbb{T}^d), \text{ and } h(\mu) \in H^1(\mathbb{T}^d), \\ + \infty & \text{otherwise.} \end{cases} $$

We will apply [24, Theorem 2] with this $\mathfrak{F}$ and $g = D^2$. Let us check the assumptions of the theorem with these choices. For what concerns the estimates, we have on the one hand, in virtue of (51):

$$ \sup_k \int_0^T \mathfrak{F}(\bar{\rho}_t^k) dt < +\infty. $$

We look at bounds in time. The estimate (41) implies that all the curves $t \mapsto \rho_k(t)$ are Hölder continuous with a constant independent of $k$, more precisely $D^2(\rho_k(t), \rho_k(s)) \le K|t-s|$. Moreover, the curves $t \mapsto \bar{\rho}_k(t)$ are piecewise constants on intervals of length $\tau_k$ and coincide with the previous curves at the end of these intervals. Hence we obtain

$$ D^2(\bar{\rho}_k(s), \bar{\rho}_k(t)) = D^2(\rho_k(\lceil s/\tau_k \rceil \tau_k), \rho_k(\lceil t/\tau_k \rceil \tau_k)) \le K(\lceil t/\tau_k \rceil - \lceil s/\tau_k \rceil)\tau_k. $$

This allows to obtain easily the condition

$$ \lim_{\tau \to 0} \sup_k \int_0^{T-\tau} D^2(\bar{\rho}_k(t), \bar{\rho}_k(t+\tau)) dt = 0, \qquad (55) $$

since we have

$$ \int_0^{T-\tau} D^2(\bar{\rho}_k(t), \bar{\rho}_k(t+\tau)) dt \le K \tau_k \int_0^T (\lceil (t+\tau)/\tau_k \rceil - \lceil t/\tau_k \rceil) dt $$

and the integral in the right hand side equals $\tau T/\tau_k$ if $\tau < \tau_k$ and is bounded by $T[\tau/\tau_k]$ otherwise. In both case, the right hand side in the equation above is bounded by $2K(T + \sup_k \tau_k)\tau$ and (55) follows.
Therefore, in order to meet the assumptions of [24, Theorem 2] and hence to conclude (54), we just need to show that $\mathfrak{F}$ has compact sublevels in $L^1(\mathbb{T}^d)$ and that $D^2$ is lower semi-continuous in $L^1(\mathbb{T}^d) \times L^1(\mathbb{T}^d)$ (the compatibility condition (1.21c) from [24] is obvious here).

Concerning $D^2$, this is obvious, since it is well known that the Wasserstein distance metrizes the weak convergence of measure, which is weaker than the strong $L^1$ topology.

Concerning $\mathfrak{F}$, we show that for a given $M > 0$, the set

$$ \mathcal{K}_M := \{ \mu \in L^1(\mathbb{T}^d) \text{ such that } \mathfrak{F}(\mu) \le M \} $$

is a compact set of $L^1(\mathbb{T}^d)$.

When $\alpha > 0$, let $(\mu_n)_{n \in \mathbb{N}}$ be a sequence in $\mathcal{K}_M$ and for all $n$, $a_n := \sqrt{\mu_n}$. We have $\int a_n^2 dx = 1$ and $\int |\nabla a_n|^2 dx \le M$. Therefore, $(a_n)_{n \in \mathbb{N}}$ is bounded in $H^1(\mathbb{T}^d)$, hence precompact in $L^2$, and therefore $(\mu_n)_{n \in \mathbb{N}}$ is precompact in $L^1(\mathbb{T}^d)$. Using classical semi-continuity arguments, we see that it limit points are in $\mathcal{K}_M$ and the conclusion follows.

When $\alpha = 0$, this is slightly more difficult and we have to use condition (14) as well as Lemma 4.1. We consider once again $(\mu_n)_{n \in \mathbb{N}}$ a sequence in $\mathcal{K}_M$ and we show that it has a limit point in $\mathcal{K}_M$ in the strong topology of $L^1$. For all $n$, let $a_n := h(\rho_n)$.

We first prove that $(a_n)$ is bounded in $L^p(\mathbb{T}^d)$ with $p = 2d/(d - 2)$. By the first point of Lemma 4.1, $((a_n)_+)$ is bounded in $L^p(\mathbb{T}^d)$. Then, we deduce that $\mathcal{K}_M$ is uniformly equi-integrable. Indeed, whenever $K > 1$, using (14),

$$ \int_{\mu_n \ge K} \mu_n \, dx = \int_{\mu_n \ge K} \frac{\mu_n}{h(\mu_n)^p} h(\mu_n)^p \, dx \le \|(a_n)_+\|_{L^p(\mathbb{T}^d)}^p \sup_{s \ge K} \frac{s}{h(s)^p} \xrightarrow{s \to +\infty} 0. $$

Therefore, the second and third point of the lemma apply with $A = \mathcal{K}_M$, so we find that $(a_n)$ is bounded in $L^p(\mathbb{T}^d)$ and hence in $H^1(\mathbb{T}^d)$.

Thus, up to considering a subsequence, we have $h(\mu_n) \to h(\mu)$ in $L^2(\mathbb{T}^d)$, and so this convergence also holds in measure. But $h^{-1}: \mathbb{R} \to \mathbb{R}_+$ is continuous, so $\mu_n \to \mu$ in measure as well. As in addition, as we already saw, $\mathcal{K}_M$ is uniformly equi-integrable, $\mu_n \to \mu$ in $L^1(\mathbb{T}^d)$. Once again, the fact that the limit points belong to $\mathcal{K}_M$ follow from standard semi-continuity arguments.

**Step 11.** Strong compactness of $g(\bar{\rho}^k)$ in $L^1([0, T] \times \mathbb{T}^d)$.

We can finally conclude the proof by showing that

$$ g(\bar{\rho}^k) \xrightarrow{k \to +\infty} g(\rho) \quad \text{in} \quad L^1([0, T] \times \mathbb{T}^d). $$

By the previous step, $\bar{\rho}^k \to \rho$ in the strong topology of $L^1([0, T] \times \mathbb{T}^d)$. On the other hand, $g$ is continuous. Therefore, if we can show that $(g(\bar{\rho}^k))$ is uniformly equi-integrable on $[0, T] \times \mathbb{T}^d$, then we can conclude that $g(\bar{\rho}^k) \to g(\rho)$ in $L^1([0, T] \times \mathbb{T}^d)$.

This uniform equi-integrability is obtained showing that both the positive and negative parts of $(g(\bar{\rho}^k) - g(1))$ are uniformly integrable. Observe that $s \mapsto g(s) - g(1)$ is nonnegative when $s \ge 1$ and nonpositive when $s \le 1$. First, we have for all $s \ge 1$:

$$ g(s) = g(1) + \int_{1}^{s} \sqrt{r} h'(r) \, dr \le g(1) + \sqrt{s} h(s). $$
We deduce that

$$ (g(\bar{\rho}^k) - g(1))_+ \le |g(1)| + \sqrt{\bar{\rho}^k} h(\bar{\rho}^k)_+. $$

But $(\sqrt{\bar{\rho}^k})$ is bounded in $L^\infty(\mathbb{R}_+; L^2(\mathbb{T}^d))$ and by (51) and the first point of Lemma 4.1, $(h(\bar{\rho}^k)_+)$ is bounded in $L^2([0, T]; L^p(\mathbb{R}^d))$ for $p = 2d/(d-2)$. So $(g(\bar{\rho}^k)_+)$ is bounded in $L^2([0, T]; L^p(\mathbb{T}^d))$ for some $p' > 1$ and uniform equi-integrability for the positive part of $(g(\bar{\rho}^k))$ on $[0, T] \times \mathbb{T}^d$ follows.

It remains to prove uniform equi-integrability of the negative part of $(g(\bar{\rho}^k))$. This is where we need to use the second point of Assumption 2.7. The easiest case is when $d_- = 0$ and $f(0) < +\infty$. Indeed, in this case, by convexity of $f$, for all $s \in [0, 1]$,

$$ -g(s) = f(s) - s f'(s) \le f(0), $$

so that for all $k$

$$ (g(\bar{\rho}^k) - g(1))_- \le |g(1)| + |f(0)|. $$

Therefore, $((g(\bar{\rho}^k) - g(1))_-)$ is bounded in $L^\infty(\mathbb{R}_+ \times \mathbb{T}^d)$ and uniform equi-integrability follows.

In all the other cases, we will use the following inequality: if $s < 1$,

$$ g(1) - g(s) = \int_s^1 \sqrt{r} h'(r) dr \le h(s)_{-}. $$

So we need to prove uniform equi-integrability for $(h(\bar{\rho}^k)_-)$.

We actually prove that $(h(\bar{\rho}^k)_-) $ is bounded in $L^2([0, T]; L^2(\mathbb{T}^d))$ using (39) and the second and third points of Lemma 4.1 with a well chosen $\mathcal{A}$. First, let us fix $K$ a large constant bounding from above the right-hand side of (39) uniformly in $k$. We call

$$ \mathcal{A} := \begin{cases} \{ \rho \in L^1(\mathbb{T}^d) \cap \mathcal{P}(\mathbb{T}^d) \text{ such that } H(\rho) \le 2K/(\inf_k \alpha_k) \}, & \text{if } \alpha > 0, \\ \{ \rho \in L^1(\mathbb{T}^d) \cap \mathcal{P}(\mathbb{T}^d) \text{ such that } \int f(\rho) \le K + \|V\|_{\infty} + \|W\|_{\infty} \}, & \text{if } \alpha = 0. \end{cases} $$

If $\alpha > 0$ or if $\alpha = 0$ but $f$ is superlinear, $\mathcal{A}$ is uniformly equi-integrable and hence, according to the third point of Lemma 4.1, $\mathcal{A}$ satisfies the second point of Lemma 4.1. If $\alpha = 0$ but $f(d_-) = +\infty$, let us call $\tilde{f}: s \mapsto f(s) - f(1) - (s-1)f'(1) \in \mathbb{R}_+$. Then for all $s_- \in (d_-, 1)$ and all $\rho \in \mathcal{A}$,

$$ \text{Leb}(\rho < s_{-}) \tilde{f}(s_{-}) \le \int \tilde{f}(\rho) = \int f(\rho) - f(1) \le K + \|V\|_{\infty} + \|W\|_{\infty} + |f(1)|. $$

Therefore, for $s_-$ sufficiently close to $d_-$, $\text{Leb}(\rho < s_-) \le 1/2$ uniformly in $\rho \in \mathcal{A}$. In that case as well, $\mathcal{A}$ satisfies the condition of the second point of Lemma 4.1.

We conclude that in all these cases, there exists $C > 0$ such that for all $\rho \in \mathcal{A}$ such that $h(\rho) \in H^1(\mathbb{T}^d)$,

$$ \|h(\rho)\|_{L^2(\mathbb{T}^d)} \le C(1 + \|\nabla h(\rho)\|_{L^2(\mathbb{T}^d)}). $$

Using (39) and (51), we find that $(\nabla h(\bar{\rho}^k))$ is bounded in $L^2([0, T]; L^2(\mathbb{T}^d))$, and that for all $k$ and almost all $t \in [0, T]$, $\bar{\rho}^k(t) \in \mathcal{A}$. Therefore, by the previous inequality, $(h(\bar{\rho}^k)_-)$ is bounded in $L^2([0, T]; L^2(\mathbb{T}^d))$. This last observation concludes the proof of Theorem 2.9. $\square$

**Acknowledgments.** The authors acknowledge the support of European union via the ERC AdG 101054420 EYAWKAJKOS.
## References

[1] S. Adams, N. Dirr, M. A Peletier, and J. Zimmer. From a large-deviations principle to the Wasserstein gradient flow: a new micro-macro passage. *Communications in Mathematical Physics*, 307:791–815, 2011.

[2] L. Ambrosio, N. Gigli, and G. Savaré. *Gradient flows in metric spaces and in the space of probability measures*. Springer, 2006.

[3] A. Baradat and H. Lavenant. Regularized unbalanced optimal transport as entropy minimization with respect to branching brownian motion. To appear in Astérisque, preprint https://arxiv.org/abs/2111.01666.

[4] A. Baradat and L. Monsaingeon. Small noise limit and convexity for generalized incompressible flows, Schrödinger problems, and optimal transport. *Arch. Ration. Mech. Anal.*, 235(2):1357–1403, 2020.

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

[6] J-D. Benamou, G. Carlier, M. Cuturi, L. Nenna, and G. Peyré. Iterative Bregman projections for regularized transportation problems. *SIAM Journal on Scientific Computing*, 37(2):A1111–A1138, 2015.

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

[8] M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In *Advances in neural information processing systems*, pages 2292–2300, 2013.

[9] S. Di Marino and A. Gerolin. An optimal transport approach for the Schrödinger bridge problem and convergence of Sinkhorn algorithm. *J. Sci. Comput.*, 85(2):Paper No. 27, 28, 2020.

[10] L. C. Evans. *Partial differential equations*, volume 19 of *Graduate Studies in Mathematics*. American Mathematical Society, Providence, RI, second edition, 2010.

[11] J. Feydy, T. Séjourné, F-X. Vialard, S. Amari, A. Trouvé, and G. Peyré. Interpolating between Optimal Transport and MMD using Sinkhorn Divergences. In *International Conference on Artificial Intelligence and Statistics*, 2018.

[12] A. Figalli and F. Glaudo. *An invitation to optimal transport, Wasserstein distances, and gradient flows*. EMS Textbooks in Mathematics. EMS Press, Berlin, 2021.

[13] A. Genevay, G. Peyré, and M. Cuturi. Learning generative models with Sinkhorn divergences. In Amos Storkey and Fernando Perez-Cruz, editors, *Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics*, volume 84 of *Proceedings of Machine Learning Research*, pages 1608–1617. PMLR, 09–11 Apr 2018.
[14] I. Gentil, C. Léonard, and L. Ripani. About the analogy between optimal transport and minimal entropy. *Annales de la Faculté des Sciences de Toulouse. Mathématiques. Série 6*, 3:569-600, 2017.

[15] N. Gigli and L. Tamanini. Benamou-Brenier and duality formulas for the entropic cost on RCD*(K, N) spaces. *Probability Theory and Related Fields*, pages 1-34, 2019.

[16] M. Hardion. Gradient flows in the geometry of the Sinkhorn divergences, 2025. master thesis.

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

[18] H. Lavenant, J. Luckhardt, G. Mordant, B. Schmitzer, and L. Tamanini. The Riemannian geometry of Sinkhorn divergences. Preprint, https://arxiv.org/abs/2405.04987.

[19] C. Léonard. From the Schrödinger problem to the Monge-Kantorovich problem. *Journal of Functional Analysis*, 262(4):1879-1920, 2012.

[20] C. Léonard. A survey of the Schrödinger problem and some of its connections with optimal transport. *Discrete & Continuous Dynamical Systems-A*, 34(4):1533-1574, 2014.

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

[22] G. Peyré. Entropic approximation of Wasserstein gradient flows. *SIAM J. Imaging Sci.*, 8(4):2323-2351, 2015.

[23] A. Ramdas, N. G. Trillos, and M. Cuturi. On Wasserstein two-sample testing and related families of nonparametric tests. *Entropy*, 19(2), 2017.

[24] R. Rossi and G. Savaré. Tightness, integral equicontinuity and compactness for evolution problems in Banach spaces. *Annali della Scuola Normale Superiore di Pisa - Classe di Scienze*, Ser. 5, 2(2):395-431, 2003.

[25] I. N. Sanov. On the probability of large deviations of random variables. Technical report, North Carolina State University. Dept. of Statistics, 1958.

[26] F. Santambrogio. Optimal transport for applied mathematicians. *Birkhäuser, NY*, 2015.

[27] F. Santambrogio. {Euclidean, metric, and Wasserstein} gradient flows: an overview. *Bull. Math. Sci.*, 7(1):87-154, 2017.

[28] R. Sinkhorn. A relationship between arbitrary positive matrices and doubly stochastic matrices. *The annals of mathematical statistics*, 35(2):876-879, 1964.

[29] R. Sinkhorn. Diagonal equivalence to matrices with prescribed row and column sums. *The American Mathematical Monthly*, 74(4):402-405, 1967.