# Wasserstein Gradient Flows of semi-discret energies: evolution of urban areas and uniform quantization

Joao Miguel Machado

► To cite this version:

Joao Miguel Machado. Wasserstein Gradient Flows of semi-discret energies: evolution of urban areas and uniform quantization. 2026. hal-05535338

HAL Id: hal-05535338

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

Preprint submitted on 4 Mar 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# WASSERSTEIN GRADIENT FLOWS OF SEMI-DISCRETE ENERGIES: EVOLUTION OF URBAN AREAS AND UNIFORM QUANTIZATION

JOÃO MIGUEL MACHADO

**ABSTRACT.** We study the Wasserstein gradient flow of semi-discrete energies in the space of probability measures, that is functionals depending on two measures– one being an absolutely continuous density and the other an atomic measure. These energies appear naturally in the field of urban planning. This is done via the celebrated JKO scheme, for which we prove convergence to a limiting system composed of a parabolic PDE with singular advection coupled with an ODE, also presenting singular dynamics. This is first done under more general assumptions using classical tools, and in a second moment convergence is proven to hold in $L_x^2 H_x^1$ for the cases of linear and Porous-Medium type diffusions. We then pass to the study of some qualitative properties of this system, such as the convergence of the atoms towards the baricenters of their corresponding Laguerre cells. We finish this work with extensive numerical simulations that aid in formulating conjectures for the qualitative behavior of this system; in the case of linear diffusion, for instance, we observe a dynamic crystallization phenomenon.

**Keywords.** Optimal Transport, Gradient Flows, Urban Planning, Optimal Quantization

**2020 Mathematics Subject Classification.** 49Q22, 35A15, 91B52

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>1.1. Contributions</td><td>3</td></tr><tr><td>Acknowledgments</td><td>4</td></tr><tr><td>2. Optimal Transport and Minimizing Movement Schemes</td><td>4</td></tr><tr><td>2.1. Optimal Transport, Wasserstein distances and the semi-discrete problem</td><td>4</td></tr><tr><td>2.2. Minimizing movement schemes</td><td>6</td></tr><tr><td>2.3. Two interpolations</td><td>7</td></tr><tr><td>3. Convergence of MMS</td><td>8</td></tr><tr><td>3.1. Optimality conditions</td><td>9</td></tr><tr><td>3.2. Strong convergence of the pressure</td><td>14</td></tr><tr><td>3.3. The limit PDE-ODE system</td><td>15</td></tr><tr><td>3.4. Strong L<sup>2</sup>H<sup>1</sup> convergence in Porous Medium case</td><td>18</td></tr><tr><td>4. Qualitative Properties</td><td>23</td></tr><tr><td>4.1. Invariant properties of the atomic dynamics</td><td>24</td></tr><tr><td>4.2. On the dynamics of the optimal quantization equation</td><td>26</td></tr><tr><td>5. Numerical Simulations and Conjectures</td><td>29</td></tr><tr><td>A crystallization phenomenon</td><td>30</td></tr><tr><td>Discussion on the long-time behavior of the densities</td><td>31</td></tr><tr><td>References</td><td>32</td></tr></table>

## 1. INTRODUCTION

In the present work we study a class of coupled PDE-ODE system that is relevant in the mathematical modeling for the evolution of an urban area and for the optimal quantization of probability measures.

A central question in the mathematical modeling of urban systems is how to describe the interplay between a population distributed over a territory and the location of a finite number of working or service sites. A variational approach to this problem was proposed in [12] by Buttazzo and Santambrogio, who introduced an energy functional describing the competition of three major
effects: congestion of the population through the penalization of its density, the cost of operating working sites, and the global cost of transportation from residences to workplaces.

In this scenario, the population density is described by an absolutely continuous probability measure $\varrho$, while the distribution of working or service sites is given by an atomic measure $\mu = \sum_{i=1}^{N} a_i \delta_{x_i}$, where $x_i$ corresponds to the position of the i-th center and $a_i$ the population percentage of population attending it. The energy described above is given by

$$ (1.1) \qquad \mathcal{E}(\varrho, \mu) \stackrel{\text{def.}}{=} \mathcal{F}(\varrho) + \mathcal{G}(\mu) + W_2^2(\varrho, \mu), $$

where

$$ (1.2) \qquad \mathcal{F}(\varrho) \stackrel{\text{def.}}{=} \begin{cases} \int_{\Omega} F(\varrho(x))dx, & \text{if } \varrho \ll \mathcal{L}^d \llcorner \Omega, \\ +\infty, & \text{otherwise,} \end{cases} $$

represents a congestion term for the population density and

$$ (1.3) \qquad \mathcal{G}(\mu) \stackrel{\text{def.}}{=} \begin{cases} \sum_{i=1}^{N} g(a_i), & \text{if } \mu = \sum_{i=1}^{N} a_i \delta_{x_i} \text{ with } (x_i)_{i=1}^{N} \subset \bar{\Omega}, \\ +\infty, & \text{otherwise,} \end{cases} $$

represents the cost of operation of a given atomic distribution of working sites. The coupling term given by the squared-Wasserstein distance represents the global transportation cost of the population to their working sites. See [40, 43] for a definition and properties, or Section 3 for a concise presentation. Besides the work of Butazzo and Santambrogio, optimal transport techniques have been vastly used to model the coupling of an absolutely continuous probability density and an atomic measure [7–10], following therefore in the category of semi-discrete transport, see [28, 35] for more details on the semi-discrete setting.

Starting from the variational principle described above, our first goal is to derive a system of evolution equations for the population density, workplaces and the corresponding proportions of the population that work on each site, which is done through the gradient flow of the energy 1.1. The equation obtained with this methodology is the following

$$ (1.4) \qquad \left\{ \begin{aligned} \partial_t \varrho_t &= \Delta P(\varrho_t) + \text{div} \left( \varrho_t \left( \sum_{i=1}^N (x - x_i(t)) \mathbb{1}_{\Omega_i(t)} \right) \right) \\ \dot{x}_i(t) &= -a_i(t)x_i(t) + \int_{\Omega_i(t)} x d\varrho_t + N_{\Omega}(x_i), \quad \text{for } i = 1, \dots, N \\ \dot{a}_i(t) &= (-g'(a_i(t)) + \psi_i(t)) \mathbb{1}_{\{a_i > 0\}}, \quad \text{for } i = 1, \dots, N \\ \Omega_i(t) &= \text{Lag}_i(\psi_i, \mathbf{x}_i), \quad \psi_i(t) \text{ is a potential for } W_2^2(\varrho_t, \mu_t), \end{aligned} \right. $$

where $N_{\Omega}(x)$ the normal cone of $\Omega$ at $x$, and $P(\rho) \stackrel{\text{def.}}{=} \rho F'(\rho) - F(\rho)$ denotes the pressure. Notice the indicator function multiplying the dynamics of $a_i$, by this we mean that if at a given time $\bar{t}$ the atom $i$ reaches null mass, then it vanishes from the dynamics from this point on.

This will be done with a minimizing movement scheme (MMS) in the product topology induced by the $\ell^2$ Wasserstein distance for the densities $\varrho$ and the $\ell^2$ euclidean norm for the sets of atoms and weights in $\Omega^N$ and the $(N-1)$-dimensional simplex. Although it can be defined in any metric space [3], in the context of optimal transport the MMS is better known as the JKO scheme, first introduced in [27, 37] to interpret the Fokker-Planck and Porous Medium equations as gradient flows in this topology. Since [3], there is a standard theory for gradient flows in the Wasserstein topology with respect to geodesically convex energies in the sense of McCann [32]. Although many interesting properties can be proven for solutions of well-established equations having this variational structure, it has also been used as a tool to prove existence of solutions to more complicated models, see [17, 31]. In these references, as well as in the present work, although the energy we study is not geodesically convex, we can still define the minimizing movement scheme and study its properties.

System 1.4 can be interpreted as follows: The first equation on $\varrho_t$ describes a competition between the tendency of diffusion of the pressure $P(\varrho_t)$ and a concentration of its mass inside each Laguerre cell $\Omega_i$ towards the atom $x_i$. The evolution describing the atoms' positions is proportional to their distance to the barycenter of there corresponding Laguerre cell. On the other hand, the
evolution of their corresponding weights is balanced by the fluctuations between the resistance to growth at the current size $g'(a_i)$ and the Kantorovitch potential $\psi_t$.

The hole system is interconnected through the optimality of the optimal transportation problem which is encoded by $\psi(t)$ and the definition of Laguerre cells

$$ \text{Lag}_i(\psi_t, \mathbf{x}_t) \stackrel{\text{def.}}{=} \left\{ x \in \Omega : \frac{1}{2}|x - x_i(t)|^2 - \psi_i(t) \le \frac{1}{2}|x - x_j(t)|^2 - \psi_j(t) \quad j = 1, \dots, N \right\}. $$

Such optimality conditions have a clear economic equilibrium interpretation; indeed $\psi_i$ can be interpreted as the average salary an individual should expect working at $x_i$. Therefore, an optimal $\psi$ achieves the equilibrim between compensation and geographic displacement under a distribution

of employment demand $\mu_t = \sum_{i=1}^N a_i(t)\delta_{x_i(t)}$ and a geographic distribution of workforce offer given by

the population $\rho_t$. In addition, it is known from the theory of optimal transport that Kantorovitch potentials are unique up to an additive constant; in our case $\psi$ is uniquely determined by the fact

$$ \text{that } \sum_{i=1}^{N} \dot{a}_i(t) = 0. $$

*Hypothesis.* In order for this problem to be meaningful from a modelling perspective and mathematically challenging, we make the following hypothesis.

(H1) $\Omega$ is a convex and bounded subset of $\mathbb{R}^d$;

(H2) $F \in \mathcal{C}^2(\mathbb{R}_+)$, is convex, with superlinear growth

$$ \lim_{t \to +\infty} \frac{F(t)}{t} = +\infty $$

and satisfies McCann’s condition of displacement convexity, namely

$$ t \mapsto t^d F(t^{-d}) \text{ is convex and nonincreasing.} $$

(H3) $g \in \mathcal{C}_{\text{loc}}^1((0, 1])$, satisfies $g(0) = 0$ and has a cusp at 0, in the sense that

$$ \lim_{t \to 0^+} \frac{g(t)}{t} = +\infty. $$

By fixing the weights to be constant equal to $a_i \stackrel{\text{def.}}{=} 1/N$ and the density penalization to be Boltzmann entropy $\mathcal{F}(\rho) = \mathcal{H}(\rho) \stackrel{\text{def.}}{=} \int_{\Omega} \rho \log \rho dx$, we obtain the following *dynamic quantization equation*

$$ (1.5) \qquad \begin{cases} \partial_t \rho = \Delta \rho + \text{div} \left( \rho_t \sum_{i=1}^N (x - x_i(t)) 1_{\Omega_i(t)} \right), \\ \dot{x}_i(t) = b_i(t) - x_i(t), \text{ for } i = 1, \dots, N, \end{cases} $$

where $b_i(t) \stackrel{\text{def.}}{=} N \int_{\Omega_i(t)} x d\rho_t$ corresponds to the barycenter of the i-th optimal Laguerre cells.

**1.1. Contributions.** As mentioned above, one of the major contributions of the present manuscript is the derivation, and proof of existence of weak solutions to the equation (1.4) via the JKO scheme. The coupling with the discrete system introduces many non-trivial steps. The first difficulty is as dealing with boundary effects, in principle the atoms can evolve on the boundary of $\Omega$ making the normal cone constraints active in their dynamics. What prevents this from happening is the fact that the atoms are pulled towards the barycenter of their Laguerre cell and the uniform in time integrability of $\rho_t$ conferred by the internal energy functional $\mathcal{F}$. The fact that $g$ has a cusp at 0 also becomes problematic to the dynamics of $a_i$. Our approach is to show that any limiting curve of the discrete scheme has the property that, whenever an atom reaches 0, it remains 0 for the rest of the evolution. Therefore, it suffices to characterize their dynamics in any open interval until the first time they vanish. In any such interval, $g'(a_i(t))$ never blows-up.

Next, following newer developments in the theory of the JKO scheme [41], we show that the JKO scheme associated to the energy (1.1) converges strongly in $L^2([0, T]; H^1(\Omega))$, instead of the classical convergence in the strong topology of $L^1([0, T] \times \Omega)$, obtained with more standard tools. This is done when the internal energy $\mathcal{F}$ is either Boltzman’s entropy or of porous medium type,
that is a diffusion of the form $\Delta\rho^m$ for some $m > 1$. Even with smoother advections, than the one coming from the semi-discrete coupled dynamics, is an improvement from [41], since up until now these types of strong convergence results have only been obtained in the linear diffusion regime.

This through convergence analysis allows us to study the qualitative properties of the equation (1.4), as well as some simplified versions of it, where the boundary effects do not take place. This is done in Section 4, where we show that atoms never touch the boundary, unless at the exact moment that their mass reaches 0. We also show that if an atom is initialized in the boundary, then it is immediately pushed inside the domain, in other words if $x_i(0) \in \partial\Omega$, then $x_i(t) \in \text{int}\Omega$ for $t$ in a sufficiently small neighborhood of 0. These properties immediately imply that the variant equation (1.5) is globally well-posed over $\mathbb{R}_+$. In this case, we show in Theorem 4.5 that the distance between atoms and their respective barycenters converge to 0 as $t \to +\infty$. This is done by combining the dissipation of energy of the gradient flow with a careful analysis of the regularity of the evolution of potentials $t \mapsto \psi_i(t)$, which implies global absolute continuity of the evolution of barycenters. This done in Lemma 2.1.

We finish the present work with extensive numerical simulations in Section 5 that not only corroborate the theoretical results we have proven, but also allow us to formulate many conjectures on the long time behavior of the system, for both linear and porous medium type diffusions.

## ACKNOWLEDGMENTS

The author wishes to thank Guillaume Carlier, Quentin Mérigot and Filippo Santambrogio for suggesting this problem and their numerous remarks that enhanced this work. He also warmly thanks the support of the Lagrange Mathematical and Computational Research Center.

## 2. OPTIMAL TRANSPORT AND MINIMIZING MOVEMENT SCHEMES

In this section we review some well known results in the literature of Optimal Transport, Wasserstein gradient flows and establish some notation used in the sequel.

### 2.1. Optimal Transport, Wasserstein distances and the semi-discrete problem.
Let $\Omega \subset \mathbb{R}^d$ be a compact, convex set, we let $\mathcal{M}(\Omega)$ denote the space of finite (scalar) Radon measures on $\Omega$ and $\mathcal{M}^d(\Omega) \stackrel{\text{def.}}{=} (\mathcal{M}(\Omega))^d$ the space of finite vector-valued Radon measures. The set of (Borel) probability measures on $\Omega$ is written as $\mathcal{P}(\Omega) \subset \mathcal{M}(\Omega)$. Given $\mu, \nu \in \mathcal{P}(\Omega)$. The 2-Wasserstein distance $W_2^2$ is defined via the value function of the quadratic optimal transport problem with the cost $c(x, y) = |x - y|^2$ and admits three equivalent formulations:

$$ (2.1) \qquad W_2^2(\mu, \nu) \stackrel{\text{def.}}{=} \min_{\gamma \in \Pi(\mu, \nu)} \int_{\Omega \times \Omega} |x - y|^2 d\gamma(x, y) = \sup_{\substack{\varphi, \psi \in \mathcal{C}_b(\Omega) \\ \varphi \oplus \psi \le c}} \int_{\Omega} \varphi d\mu + \int_{\Omega} \psi d\nu, $$

where $\Pi(\mu, \nu)$ is the set of transport plans with marginals $\mu, \nu$. It is well known that this quantity defines a distance in $\mathcal{P}(\Omega)$ which metrizes the narrow convergence of probability measures, i.e. convergence in duality with the continuous and bounded functions $\mathcal{C}_b(\Omega)$. The reader is referred to the numerous monographs on the topic for further information [3, 40, 43] for more information. In the sequel we shall give a brief discussion of the relevant properties for the subsequent work.

For the quadratic case we can always obtain existence of optimal Kantorovich potentials $(\varphi, \psi)$ for the dual problem from (2.1), see for instance [40, Prop. 1.11] for a proof in a compact domain. Optimizers for the primal and dual problems are related as

$$ \varphi(x) + \psi(y) = |x - y|^2 \quad \gamma \text{-a.e.} $$

and whenever $\mu \ll \mathcal{L}^d$, Brenier's theorem says that $\gamma = (\text{id}, T)_*\mu$ is concentrated on the graph of a map $T$ which is the gradient of a convex function

$$ T = \nabla u, \text{ and } u \text{ is related to the optimal potential } u = \frac{1}{2} |\cdot|^2 - \varphi. $$

In this case we also know that the optimal map $T$ is unique and the Kantorovich potentials $(\varphi, \psi)$ unique up to constants.
The semi-discrete variant of the optimal transportation problem consists of assuming that the first marginal is given by an absolutely continuous measure $\rho \ll \mathcal{L}^d$ and a discrete measure. Given $\mathbf{x} = (x_i)_{i=1}^N \in \Omega^N$ and $\mathbf{a} = (a_i)_{i=1}^N \in \Delta_{N-1}$ the simplex, we let

$$ \mu_{\mathbf{x},\mathbf{a}} \stackrel{\text{def.}}{=} \sum_{i=1}^{N} a_i \delta_{x_i}, $$

denote the atomic measure with atoms concentrated over $\mathbf{x}$ and weights $\mathbf{a}$. In this case, dual formulation of the transportation cost can be written as

$$ (2.2) \qquad \frac{1}{2} W_2^2(\rho, \mu_{\mathbf{x},\mathbf{a}}) = \sup_{\psi \in \mathbb{R}^N} \sum_{i=1}^N \int_{\text{Lag}_i(\mathbf{x},\psi)} \left( \frac{1}{2} |x - x_i|^2 - \psi_i \right) d\rho + \sum_{i=1}^N a_i \psi_i, $$

where $\mathbf{x} = (x_i)_{i=1}^N$ and $\text{Lag}_i(\mathbf{x}, \psi)$ denotes the $i$-th Laguerre cell induced by $\psi$ and $\mathbf{x}$, defined as

$$ (2.3) \qquad \text{Lag}_i(\mathbf{x}, \psi) \stackrel{\text{def.}}{=} \left\{ x \in \Omega : \frac{1}{2} |x - x_i|^2 - \psi_i \le \frac{1}{2} |x - x_j|^2 - \psi_j, \forall j = 1, \dots, N \right\}. $$

The maximizer is the unique $\psi \in \mathbb{R}^N$, up to adding a constant vector, such that $\rho(\text{Lag}_i(\mathbf{x}, \psi)) = a_i$ for all $i = 1, \dots, N$. These facts can be obtained via Brenier's Theorem applied to the semi-discrete case. Throughout the text, we shall fix the notation

$$ \Omega_i = \text{Lag}_i(\mathbf{x}, \psi), \quad \Sigma_{ij} = \partial \text{Lag}_i(\mathbf{x}, \psi) \cap \partial \text{Lag}_j(\mathbf{x}, \psi) \text{ for } i \ne j, $$

whenever the context is clear for instance the Laguerre tesselation is obtained from $W_2^2(\rho, \sum a_i \delta_{x_i})$ and $\psi$ is optimal, and define the corresponding optimal Kantorovich potential induced by this tesselation

$$ (2.4) \qquad \Phi[\mathbf{x}, \psi](x) \stackrel{\text{def.}}{=} \sum_{i=1}^{N} \left[ \frac{1}{2} |x - x_i|^2 - \psi_i \right] \mathbf{1}_{\Omega_i}(x). $$

Since Laguerre cells are convex polygons whose faces are determined by the points $\mathbf{x}$ and the potential $\psi$, they are stable with respect to these quantities, as stated in the following.

**Lemma 2.1.** Given a sequence $(\mathbf{x}_n, \psi_n)_{n \in \mathbb{N}}$ converging to $(\mathbf{x}, \psi)$, it holds that

$$ \text{Lag}_i(\mathbf{x}_n, \psi_n) \xrightarrow{n \to \infty} \text{Lag}_i(\mathbf{x}, \psi). $$

It will be particularly important for us to perform variations of the optimal transport cost with respect to a measure, or in the semi-discrete case compute gradients with respect to the atoms. The former can be found for instance in [40, Prop. 7.17]. To the derivative of the semi-discrete cost can be expressed in terms of the barycenter of the optimal Laguerre tesselation. Therefore, we introduce the following notation

$$ (2.5) \qquad \mathbf{b}[\mathbf{x}, \psi] = (b_i[\mathbf{x}, \psi])_{i=1}^N \quad \text{where } b_i[\mathbf{x}, \psi] \stackrel{\text{def.}}{=} \int_{\text{Lag}_i(\mathbf{x}, \psi)} x d\rho(x). $$

Differentiating the semi-discrete transport cost is a very delicate matter, the major issue arises when two points coincide, with the derivatives being more and more singular as the become too close. For this it is useful to introduce the following notation for the generalized diagonal

$$ (2.6) \qquad \mathbb{D}_{N,\varepsilon} \stackrel{\text{def.}}{=} \{|x_i - x_j| \le \varepsilon \text{ for some pair } i \ne j\}, $$

In the sequel we summarize the results which are relevant to us.

**Lemma 2.2.** Let $\rho \in \mathcal{P}_{\text{ac}}(\Omega)$, then the following hold:

(i) The first variation of $W_2^2(\rho, \nu)$ is given by the unique, up to constants, Kantorovitch potential $\varphi$ w.r.t. $\rho$, i.e. for any $\bar{\rho} \in \mathcal{P}_{\text{ac}}(\Omega)$ we have

$$ \left. \frac{d}{d\varepsilon} \right|_{\varepsilon=0^+} W_2^2(\rho + \varepsilon(\bar{\rho} - \rho), \nu) = \int_{\Omega} \varphi d(\bar{\rho} - \rho). $$
(ii) In the semi-discrete transport case, for $\rho \in \mathcal{P}_{ac}(\Omega)$ fixed, set

$$F(\mathbf{x}, \mathbf{a}) \stackrel{\text{def.}}{=} \frac{1}{2} W_2^2(\rho, \mu_{\mathbf{x},\mathbf{a}}).$$

If $\mathbf{x} \in \Omega^N \setminus \mathbb{D}_{N,\varepsilon}$ and $a_i > 0$ for all $i = 1, \dots, N$ then $\mathbf{x} \mapsto F(\mathbf{x}, \mathbf{a})$ is twice differentiable, with gradients given by

$$\nabla_{x_i} F(\mathbf{x}, \mathbf{a}) = \int_{\Omega_i} (x_i - x) d\rho = a_i (x_i - b_i[\mathbf{x}, \psi]),$$

where $\psi$ is the unique Kantorovitch potential associated with $\mu_{\mathbf{x},\mathbf{a}}$. The second derivatives are given by: for $i \neq j$

$$\nabla_{x_i, x_j} F(\mathbf{x}, \mathbf{a}) = \int_{\Sigma_{ij}} (x - x_i) \otimes (x - x_j) \frac{\rho(x)}{|x_i - x_j|} d\mathcal{H}^{d-1}(x),$$

and in the diagonal for all $i = 1, \dots, N$

$$\nabla_{x_i, x_i} F(\mathbf{x}, \mathbf{a}) = a_i I_d - \sum_{j \neq i} \nabla_{x_i, x_j} F(\mathbf{x}, \mathbf{a}).$$

In addition, if $\rho$ is continuous, then $\mathbf{x} \mapsto F(\mathbf{x}, \mathbf{a})$ is $\mathcal{C}^2$.

We do not give a proof for this result, for item (i) the reader is referred to [40, Prop.7.17], while for item (ii) the reader is referred to [19]. In the latter, there is a simple argument using the envelope theorem: letting $\psi$ denote an optimal Kantorovich potential we have

$$\begin{aligned} \nabla_{x_i} F(\mathbf{x}, \mathbf{a}) &= \nabla_{x_i} \left\{ \int_{\Omega} \min_{j=1,\dots,N} \left( \frac{1}{2} |x-x_j|^2 - \psi_j \right) d\rho + \sum_{j=1}^N \psi_j a_j \right\} \\ &= \int_{\text{Lag}_i(\mathbf{x}, \psi)} (x_i - x) d\rho = a_i(x_i - b_i[\mathbf{x}, \psi]). \end{aligned}$$

With this characterization of the gradient there is a very simple argument to see that $\mathbf{x} \mapsto F(\mathbf{x}, \mathbf{a})$ is semi-concave, see also [33]. For simplicity, let us take $a_i = 1/N$ for all $i$, let $F(\mathbf{x}), F(\mathbf{y})$ denote the minimal cost of transporting $\rho$ to the empirical measures over $\mathbf{x}$ and $\mathbf{y}$, respectively. Let $\Omega_i$ denote the optimal Laguerre tesselation associated with $F(\mathbf{x})$, then setting $\bar{F}(\mathbf{x}) = F(\mathbf{x}) - \frac{1}{2N} \|\mathbf{x}\|^2$ we have

$$\begin{aligned} F(\mathbf{y}) &\le \sum_{i=1}^N \int_{\Omega_i} \frac{1}{2} |x - y_i|^2 d\rho = F(\mathbf{x}) + \sum_{i=1}^N \int_{\Omega_i} \left( \frac{1}{2} |x - y_i|^2 - \frac{1}{2} |x - x_i|^2 \right) d\rho \\ &\le F(\mathbf{x}) - \frac{1}{N} \sum_{i=1}^N b_i \cdot (y_i - x_i) + \frac{1}{2N} \sum_{i=1}^N |y_i|^2 - |x_i|^2. \end{aligned}$$

As a result

$$\bar{F}(\mathbf{y}) \le \bar{F}(\mathbf{x}) + \nabla_{\mathbf{x}} \bar{F}(\mathbf{x}) \cdot (\mathbf{y} - \mathbf{x}),$$

which shows that $\mathbf{x} \mapsto F(\mathbf{x})$ is $1/(2N)$ semi-concave.

## 2.2. Minimizing movement schemes.
Our approach to show existence of the PDE-ODE system, we use the so called minimizing movement scheme. It consists of an implicit Euler scheme that can be easily adapted to metric spaces, the reader is referred to [3] for more details on this general formulation.

For simplicity, we introduce the notation $\mathbf{z} = (\rho, \mathbf{x}, \mathbf{a}) \in \mathcal{P}_{ac}(\Omega) \times \Omega^{\otimes N} \times \Delta_{N-1}$, where $\mathbf{x} = (x_i)_{i=1}^N$ and $\mathbf{a} = (a_i)_{i=1}^N$. The set $\Delta_{N-1}$ denotes the $(N-1)$-dimensional simplex. Consider the distance between $\mathbf{z} = (\rho, \mathbf{x}, \mathbf{a})$ and $\mathbf{z}' = (\rho', \mathbf{x}', \mathbf{a}')$ in this product space defined as

$$(2.7) \qquad d_{W_2, \ell_2}^2(\mathbf{z}, \mathbf{z}') \stackrel{\text{def.}}{=} W_2^2(\rho, \bar{\rho}) + |\mathbf{x} - \mathbf{x}'|^2 + |\mathbf{a} - \mathbf{a}'|^2,$$

where $|\cdot|$ denotes the euclidean distance. Given some $\mathbf{z} = (\rho, \mathbf{x}, \mathbf{a})$, and we make the slight abuse of writing $\mathcal{E}(\mathbf{z}) = \mathcal{E}(\rho, \mu_{\mathbf{x},\mathbf{a}})$.
A minimizing movement scheme in this topology consists of a sequence $(\varrho_k^\tau, \mathbf{x}_k^\tau, \mathbf{a}_k^\tau)_{k \in \mathbb{N}}$ such that $(\varrho_0^\tau, \mathbf{x}_0^\tau, \mathbf{a}_0^\tau) = (\varrho_0, \mathbf{x}_0, \mathbf{a}_0)$ is given and for $k \ge 0$ we have

$$ (2.8) \qquad (\varrho_{k+1}^\tau, \mathbf{x}_{k+1}^\tau, \mathbf{a}_{k+1}^\tau)_{k \in \mathbb{N}} \in \underset{\tau \ge 0}{\operatorname{argmin}} \mathcal{F}(\varrho) + \mathcal{G}(\mu_{\mathbf{x},\mathbf{a}}) + W_2^2(\varrho, \mu_{\mathbf{x},\mathbf{a}}) \\ + \frac{1}{2\tau} \left( W_2^2(\varrho_k^\tau, \varrho) + \|\mathbf{x}_k^\tau - \mathbf{x}\|^2 + \|\mathbf{a}_k^\tau - \mathbf{a}\|^2 \right). $$

It is not difficult to show existence of minimizers for each step of this scheme, however we cannot in general expect uniqueness due to the fact that this energy is semi-concave in $\mathbf{x}$, as discussed in the previous section. Whenever the context is clear we will omit the dependence on $\tau$, and we use also the notation $\mu_k = \mu_{\mathbf{x}_k, \mathbf{a}_k}$ to designate the atomic measure induced by an element of the minimizing movement scheme.

The following result comes directly from the definition and is standard.

**Lemma 2.3.** *Assuming that $\mathcal{E}(\varrho_0, \mu_0) < +\infty$, there is a constant $C > 0$ depending only on $\Omega$ and $T$ such that*

$$ (2.9) \qquad \mathcal{E}(\varrho_{k+1}^\tau, \mu_{k+1}^\tau) \le \mathcal{E}(\varrho_k^\tau, \mu_k^\tau) \le \dots \le \mathcal{E}(\varrho_0, \mu_0) \le C, \quad \text{for all } k \in \mathbb{N}, \tau > 0 $$

$$ (2.10) \qquad \sum_{k=0}^{N_\tau-1} d_{W_2, \ell_2}^2 (\mathbf{z}_k^\tau, \mathbf{z}_{k+1}^\tau) \le C\tau. $$

*Proof.* To prove (2.9), use $(\varrho_k^\tau, \mu_k^\tau)$ as a competitor for the problem defining $(\varrho_{k+1}^\tau, \mu_{k+1}^\tau)$ to obtain for all $k$ that

$$ \frac{1}{2\tau} d_{W_2, \ell_2}^2 (\mathbf{z}_k^\tau, \mathbf{z}_{k+1}^\tau) \le \mathcal{E}(\varrho_k^\tau, \mu_k^\tau) - \mathcal{E}(\varrho_{k+1}^\tau, \mu_{k+1}^\tau). $$

Since the distance $d_{W_2, \ell_2}$ is non-negative (2.9) follows from a simple induction.

Conversely, sum the above inequality in $k$, since the right-hand side telescopes and $\mathcal{E} \ge 0$, we get that

$$ \sum_{k=0}^{N_\tau-1} d_{W_2, \ell_2}^2 (\mathbf{z}_k^\tau, \mathbf{z}_{k+1}^\tau) \le 2\mathcal{E}(\varrho_0, \mu_0)\tau, $$

and the result follows.

□

**2.3. Two interpolations.** Given a time interval $[0, T]$ and a sequence $(\varrho_k^\tau, \mathbf{x}_k^\tau, \mathbf{a}_k^\tau)_{k \in \mathbb{N}}$ obtained via the minimizing movement scheme, we consider two time interpolations namely

*   **Staircase interpolation:**

$$ (2.11) \qquad (\varrho_t^\tau, \mathbf{x}_t^\tau, \mathbf{a}_t^\tau) \stackrel{\text{def.}}{=} (\varrho_k^\tau, \mathbf{x}_k^\tau, \mathbf{a}_k^\tau), \quad \text{if } t \in (k\tau, (k+1)\tau] $$

*   **Geodesic interpolation:** The second family of curves $(\bar{\varrho}^\tau, \bar{\mathbf{x}}^\tau, \bar{\mathbf{a}}^\tau)_{\tau>0}$ now contained in $\mathcal{C}([0, T]; \mathcal{P}_{\text{ac}}(\Omega))$ is defined as the geodesic between $(\varrho_k^\tau, \mathbf{x}_k^\tau, \mathbf{a}_k^\tau)$ and $(\varrho_{k+1}^\tau, \mathbf{x}_{k+1}^\tau, \mathbf{a}_{k+1}^\tau)$ at times $k\tau$ and $(k+1)\tau$ in the topology induced by $d_{W_2, \ell_2}$. Hence, setting

$$ (2.12) \qquad s_t \stackrel{\text{def.}}{=} \frac{1}{\tau}((k+1)\tau - t), \quad T_{k+1}^\tau \text{ is the OT map from } \varrho_{k+1}^\tau \text{ to } \varrho_k^\tau, $$

we define the geodesic interpolation as

$$ (2.13) \qquad \bar{\varrho}_t^\tau \stackrel{\text{def.}}{=} (T_{k+1}^\tau(t))_\# \varrho_{k+1}^\tau, \quad \text{with } T_{k+1}^\tau(t) \stackrel{\text{def.}}{=} ((1-s_t)\text{id} + s_t T_{k+1}^\tau) $$

$$ (2.14) \qquad \bar{\mathbf{x}}_t^\tau \stackrel{\text{def.}}{=} (1 - s_t)\mathbf{x}_{k+1}^\tau + s_t \mathbf{x}_k^\tau $$

$$ (2.15) \qquad \bar{\mathbf{a}}_t^\tau \stackrel{\text{def.}}{=} (1 - s_t)\mathbf{a}_{k+1}^\tau + s_t \mathbf{a}_k^\tau. $$

While the staircase interpolation is easier to relate to the optimality conditions, see Prop. 3.1, the geodesic interpolation solves by construction the continuity equation with a suitable velocity, as this is a general property of geodesics in the Wasserstein space, see [40, Thm. 5.14] and [2, Chap. 8].

Hence we define also the set of velocities, starting with a staircase velocity field

$$ (2.16) \qquad \mathbf{v}_t^\tau = \mathbf{v}_{(k+1)}^\tau = \frac{\text{id} - T_{k+1}^\tau}{\tau}, \quad \text{for } t \in (k\tau, (k+1)\tau]. $$
The velocity of the geodesic interpolation can be computed as follows: for an arbitrary $\phi \in \mathcal{C}_c^\infty(\Omega)$ and $t \in (k\tau, (k+1)\tau]$ we have

$$ \begin{aligned} \frac{d}{dt} \int_{\Omega} \phi d\bar{\varrho}_{t}^{\tau} &= \frac{d}{dt} \int_{\Omega} \phi \circ T_{k+1}^{\tau}(t) d\varrho_{k+1}^{\tau} = \int_{\Omega} \nabla \phi \circ T_{k+1}^{\tau}(t) \cdot \left( \frac{\mathrm{id} - T_{k+1}^{\tau}}{\tau} \right) d\varrho_{k+1}^{\tau} \\ &= \int_{\Omega} \nabla \phi \cdot \underbrace{\left( \frac{\mathrm{id} - T_{k+1}^{\tau}}{\tau} \right) \circ (T_{k+1}^{\tau}(t))^{-1}}_{=v_t^\tau} d\bar{\varrho}_t^{\tau}. \end{aligned} $$

Therefore, the geodesic velocity interpolation can be equivalently rewritten as

$$ (2.17) \qquad \mathbf{v}_t^{\tau} \stackrel{\text{def.}}{=} \mathbf{v}_{(k+1)}^{\tau} \circ \left( \mathrm{id} + (t - (k+1)\tau)\mathbf{v}_{(k+1)}^{\tau} \right)^{-1}, \quad \text{for } t \in (k\tau, (k+1)\tau], $$

and $(\varrho^{\tau}, \mathbf{v}^{\tau})$ solve the continuity equation.

This feature can also be expressed in terms of the momentum variables, that is the vector measures in $\mathcal{M}^d([0, T] \times \Omega)$, defined as

$$ (2.18) \qquad E^{\tau} \stackrel{\text{def.}}{=} \mathbf{v}^{\tau} \varrho^{\tau} \quad \text{and} \quad \bar{E}^{\tau} \stackrel{\text{def.}}{=} \bar{\mathbf{v}}^{\tau} \bar{\varrho}^{\tau}, $$

it holds trivially that

$$ \partial_t \bar{\varrho}_t^{\tau} + \operatorname{div} \bar{E}^{\tau} = 0, \quad \text{for all } \tau > 0, $$

in the sense of distributions.

Using the a priori estimates from Lemma 2.3 we have the following result which says that both families of interpolations $(\varrho^{\tau}, \mathbf{x}^{\tau}, \mathbf{a}^{\tau}, E^{\tau})_{\tau>0}$ and $(\bar{\varrho}^{\tau}, \bar{\mathbf{x}}^{\tau}, \bar{\mathbf{a}}^{\tau}, \bar{E}^{\tau})_{\tau>0}$ enjoy compactness properties and their limits coincide.

**Proposition 2.4 (Chap.8 of [40]).** *There is a curve*

$$ (\varrho_t, \mathbf{x}_t, \mathbf{a}_t, E_t)_{t \in [0, T]} \in \mathcal{C}^{0,1/2}([0, T]; \mathcal{P}_{\text{ac}}(\Omega) \times \bar{\Omega} \times \Delta_{N-1}) \times \mathcal{M}^d([0, T] \times \Omega) $$

which is the limit, up to subsequences that are not relabelled, of both families of interpolations $(\varrho^{\tau}, \mathbf{x}^{\tau}, \mathbf{a}^{\tau}, E^{\tau})_{\tau>0}$ and $(\bar{\varrho}^{\tau}, \bar{\mathbf{x}}^{\tau}, \bar{\mathbf{a}}^{\tau}, \bar{E}^{\tau})_{\tau>0}$. More precisely, it holds that

$$ \varrho^{\tau}, \bar{\varrho}^{\tau} \xrightarrow{\tau \to 0} \varrho \text{ uniformly in } [0, T], $$

$$ (\mathbf{x}^{\tau}, \mathbf{a}^{\tau}), (\bar{\mathbf{x}}^{\tau}, \bar{\mathbf{a}}^{\tau}) \xrightarrow{\tau \to 0} (\mathbf{x}, \mathbf{a}) \text{ uniformly in } [0, T], $$

$$ E^{\tau}, \bar{E}^{\tau} \xrightarrow{\tau \to 0} \text{ in the narrow topology of } \mathcal{M}^d([0, T] \times \Omega). $$

In addition, $\partial_t \varrho_t + \operatorname{div} E_t = 0$, in the sense of distributions.

### 3. CONVERGENCE OF MMS

Given the results of Section 2, in other to prove the convergence of the minimizing movement scheme to the coupled system (1.4), we only need to characterize the limit momentum variable $E$ and the dynamics of the variables $\mathbf{x}, \mathbf{a}$. This will be done in the sequel with the Euler-Lagrange equations characterizing the optimality of the minimizing movement.

However, there are still multiple difficulties in characterizing the limit curve. To deal with the non-linear term $\Delta P(\varrho)$, we will require a stronger convergence of $\varrho^{\tau}$ to $\varrho$, which can be done with well-established results in the literature.

The evolution of the discrete measure is more subtle. Regarding the evolution of $\mathbf{x}$, we need to be careful with the boundary effects, since the points are restricted to $\Omega$, the optimality conditions must push points to the interior of the domain if the minimizer of any step is on the boundary. These boundary effects however cannot be passed onto the limit as $\tau \to 0$, unless we assume that $\partial\Omega$ is smooth. Instead, we show in Prop. 3.2 that if at an iteration of the minimizing movement scheme a point is in the interior of $\Omega$, then it remains in the interior as long as its mass remains positive.

This leads us to the heart of the matter: what about the masses $\mathbf{a}$? As the dynamics around $\mathbf{a}$ contain the singular term $g'$, it is not clear what happens to a point when its mass reaches 0. We show in Prop. 3.3 that if $(\varrho_t, \mathbf{x}_t, \mathbf{a}_t)_{t \ge 0}$ is a curve obtained from the minimizing movement scheme, then if at some time $a_i(t) = 0$, then it remains null of any $s > t$. This means that we only need to characterize the dynamics of $a_i$ until the first time that it vanishes.
**3.1. Optimality conditions.** In order to derive optimality conditions, we will need the notion of *Bouligand's tangent cone*

$$ (3.1) \qquad T_{\Omega}(z) \stackrel{\text{def.}}{=} \{v = \lim_{n \to \infty} v_n : z + \varepsilon_n v_n \in \Omega, \varepsilon_n \to 0\} $$

see [26, Def. 5.1.1], which gives the admissible directions to perform variations in $\Omega$. Likewise, the normal cone is then given by the polar cone of $T_{\Omega}(z)$, see for instance [26, Part I, Prop. 5.2.4]. In other words, it can be written as

$$ (3.2) \qquad N_{\Omega}(z) \stackrel{\text{def.}}{=} \{n \in \mathbb{R}^d : \langle n, v \rangle \le 0 \text{ for all } v \in T_{\Omega}(z)\}. $$

Whenever $\Omega$ has a smooth boundary $N_{\Omega}(z)$ is generated by the outwards normal vector.

**Proposition 3.1.** Let $(\varrho_{k+1}^{\tau}, \mathbf{x}_{k+1}^{\tau}, \mathbf{a}_{k+1}^{\tau})$ be a minimizer of (2.8), and define the set

$$ (3.3) \qquad I_{k+1} \stackrel{\text{def.}}{=} \{i = 1, \dots, N : 0 < a_{k+1,i}^{\tau} < 1\}, $$

corresponding to the indexes for which the constraints on the weights are not active. Then it holds that

$$ \left\{ \begin{array}{ll} -\mathbf{v}_{(k+1)}^{\tau} = \nabla F'(\varrho_{k+1}^{\tau}) + \sum_{i=1}^{N} (x - x_{k+1,i}^{\tau}) \mathbb{1}_{\Omega_{k+1,i}^{\tau}} & \varrho_{k+1}^{\tau} - \text{a.e.} \\ \\ 0 \in \left( \frac{x_{k+1,i}^{\tau} - x_{k,i}^{\tau}}{\tau} + a_{k+1,i}^{\tau} \cdot x_{k+1,i}^{\tau} - \int_{\Omega_{k+1,i}^{\tau}} x d\varrho_{k+1}^{\tau} \right) + N_{\Omega}(x_{k+1,i}), & i = 1, \dots, N \\ \\ -\frac{a_{k+1,i}^{\tau} - a_{k,i}^{\tau}}{\tau} = g'(a_{k+1,i}^{\tau}) + \psi_{k+1,i}^{\tau}, & \text{for all } i \in I_{k+1} \\ \\ \varrho_{k+1,i}^{\tau} = \text{Lag}_i (\psi_{k+1,i}^{\tau}, \mathbf{x}_{k+1,i}^{\tau}), & i = 1, \dots, N. \end{array} \right. $$

where $\psi_{k+1}^{\tau}$ is a Kantorovitch potential associated with $\mu_{k+1}^{\tau}$ in $W_2^2(\varrho_{k+1}^{\tau}, \mu_{k+1}^{\tau})$.

In the proof of this result, since $\tau$ is fixed we omit shall it, in order to simplify notation.

**Proof.** The optimality conditions for $\varrho_{k+1}$ are very similar to the classical theory in bounded domains, see [40, Chap. 7,8], hence we focus mostly on the analysis of $\mathbf{x}_{k+1}$, and $\mathbf{a}_{k+1}$.

**Optimality conditions for $\varrho$:** Since for any $k \in \mathbb{N}$, the internal energy $\mathcal{F}$ enforces that $\varrho_{k+1} \ll \mathcal{L}^d$, we can use the first variation formulas from Lemma 2.2. Let us first derive the optimality conditions for $\varrho_{k+1}$, so given $\rho \in \mathcal{P}_{\text{ac}}(\Omega)$ consider a variation $\varrho_{\varepsilon} \stackrel{\text{def.}}{=} \varrho_{k+1} + \varepsilon(\rho - \varrho_{k+1})$. Optimality gives that

$$ -\frac{1}{2\tau} \left. \frac{\mathrm{d}}{\mathrm{d}\varepsilon} \right|_{\varepsilon=0^+} W_2^2(\varrho_k, \varrho_\varepsilon) \le \left. \frac{\mathrm{d}}{\mathrm{d}\varepsilon} \right|_{\varepsilon=0^+} \mathcal{F}(\varrho_\varepsilon) + \left. \frac{\mathrm{d}}{\mathrm{d}\varepsilon} \right|_{\varepsilon=0^+} \frac{1}{2} W_2^2(\varrho_\varepsilon, \mu_{k+1}). $$

Using Lemma 2.2, we obtain that

$$ -\frac{\psi_{k+1}^{\varrho}}{\tau} = F'(\varrho_{k+1}) + \varphi_{k+1}^{\varrho, \mu}, $$

where $\psi_{k+1}^{\varrho}$ and $\varphi_{k+1}^{\varrho, \mu}$ are the Kantorovitch potentials associated with $\varrho_{k+1}$ in $W_2^2(\varrho_k, \varrho_{k+1})$ and $W_2^2(\varrho_{k+1}, \mu_{k+1})$, respectively.

On the other hand, from the optimality conditions of the optimal transportation problem, it holds that

$$ \varphi_{k+1}^{\varrho, \mu}(x) + \psi_{k+1,i} = \frac{1}{2} \|x - x_{k+1,i}\|^2, \quad \text{for } \varrho_{k+1} - \text{a.e. } x \in \Omega_{k+1,i}. $$

Since $\Omega_{k+1,i}$ is a convex polyhedra, its boundary is $\varrho_{k+1}$ negligible, and it follows that

$$ \nabla \varphi_{k+1}^{\varrho, \mu}(x) = \sum_{i=1}^{N} (x - x_{k+1,i}) \mathbb{1}_{\Omega_{k+1,i}}. $$
In addition, from Brenier's Theorem we know that the optimal map from $\varrho_{k+1}$ to $\varrho_k$ is given by $T_{k+1} = \text{id} - \nabla \psi_{k+1}^\varrho$. Combining these results we obtain that

$$ \frac{T_{k+1} - \text{id}}{\tau} = \nabla F'(\varrho_{k+1}) + \sum_{i=1}^{N} (x - x_{k+1,i}) \mathbb{1}_{\Omega_{k+1,i}}. $$

**Optimality conditions for x:** To derive the optimality conditions w.r.t. the atoms $\mathbf{x}_{k+1} = (x_{k+1,i})_{i=1}^N$, fix some $i$ and consider a direction in the *Bouligand's tangent cone*, $v \in T_\Omega(x_{k+1,i})$, where

$$ T_\Omega(z) \stackrel{\text{def.}}{=} \{v = \lim_{n \to \infty} v_n : z + \varepsilon_n v_n \in \Omega, \varepsilon_n \to 0\}, $$

see [26, Def. 5.1.1]. Hence, take $\varepsilon_n$ and $v_n$ converging to $v$ as in the definition above and set

$$ \bar{x}_{j, \varepsilon_n} \stackrel{\text{def.}}{=} \begin{cases} x_{k+1,j}, & \text{if } j \neq i, \\ x_{k+1,i} + \varepsilon_n v_n, & \text{j = i,} \end{cases} $$

and the define a variation of the atomic measures as

$$ \mu_{\varepsilon_n} \stackrel{\text{def.}}{=} \sum_{j=1}^{N} a_{k+1,j} \delta_{\bar{x}_{j, \varepsilon_n}}. $$

Comparing the energies of $\mu_{k+1}$ and $\mu_{\varepsilon_n}$ with $\varrho_{k+1}$ fixed, we obtain from Lemma 2.2 that

$$ \begin{aligned} 0 &\le \lim_{n \to \infty} \left( \frac{1}{2\tau} \frac{\|\bar{x}_{i, \varepsilon_n} - x_{k,i}\|^2}{\varepsilon_n} - \|x_{k+1,i} - x_{k,i}\|^2 + \frac{W_2^2(\varrho_{k+1}, \mu_{\varepsilon_n}) - W_2^2(\varrho_{k+1}, \mu_{k+1})}{\varepsilon_n} \right) \\ &= \lim_{n \to \infty} \left( \left\langle \frac{x_{k+1,i} - x_{k,i}}{\tau}, v_n \right\rangle + \partial_{x_i} \frac{1}{2} W_2^2(\varrho_{k+1}, \mu_{k+1}) \cdot v_n + \frac{o(\varepsilon_n)}{\varepsilon_n} \right) \\ &= \left\langle \frac{x_{k+1,i} - x_{k,i}}{\tau} + \int_{\text{Lag}_i(\psi_{k+1}; \mathbf{x}_{k+1})} (x_{k+1,i} - x) d\varrho_{k+1}, v \right\rangle. \end{aligned} $$

Recalling, for instance from [26, Part I, Prop. 5.2.4], that the normal space $N_\Omega(x_{k+1,i})$ is given by the polar cone to $T_\Omega(x_{k+1,i})$, that $\Omega_{k+1,i} = \text{Lag}_i(\psi_{k+1}, X_{k+1})$ and $\varrho_{k+1}(\Omega_{k+1,i}) = a_{k+1,i}$ we obtain that for all $i = 1, \dots, N$

$$ - \left( \frac{x_{k+1,i} - x_{k,i}}{\tau} + a_{k+1,i} \cdot x_{k+1,i} - \int_{\Omega_{k+1,i}} x d\varrho_{k+1} \right) = n_{k+1,i} \in N_\Omega(x_{k+1,i}). $$

**Optimality conditions for a:** In the sequel, we perform variations for the mass variables $(a_{k+1,i})_{i=1}^N$. Consider two indexes $i, j \in I_{k+1}$ and construct the variation $\mathbf{a}_\varepsilon = (a_{\varepsilon,h})_{h=1}^N$ given by

$$ a_{\varepsilon,h} \stackrel{\text{def.}}{=} \begin{cases} a_{k+1,h}, & h \neq i, j, \\ a_{k+1,j} - \varepsilon, & h = j, \\ a_{k+1,i} + \varepsilon, & h = i. \end{cases} $$

And define the new measure $\mu_\varepsilon \stackrel{\text{def.}}{=} \mu_{\mathbf{x}_{k+1}, \mathbf{a}_\varepsilon}$. For $\varepsilon$ small enough this variation belongs to the simplex $\Delta_{N-1}$ and hence is admissible. As a result, comparing the energies of $(\varrho_{k+1}, \mu_{k+1})$ and $(\varrho_{k+1}, \mu_\varepsilon)$ we obtain that there is a Kantorovitch potential $\bar{\psi}_{k+1}$ such that

$$ \frac{a_{k+1,i} - a_{k,j}}{\tau} + g'(a_{k+1,i}) + \bar{\psi}_{k+1,i} \ge \frac{a_{k+1,j} - a_{k,j}}{\tau} + g'(a_{k+1,j}) + \bar{\psi}_{k+1,j}. $$

In addition, we can assume that $\bar{\psi}_{k+1}$ is bounded in $\ell^\infty$-norm by a constant $C_\Omega$, since by Brenier's Theorem it the evaluation of a Lipschitz function over the points $\{x_{k+1,i}\}$ and uniquely defined up to constants, so that we can assume that $\bar{\psi}_{k+1,1} = 0$ and the bound follows from the Lipschitz continuity and the fact that $\Omega$ has finite diameter.

Changing the role of $i$ and $j$, it follows that we have the equality

$$ \begin{aligned} & \frac{a_{k+1,i} - a_{k,j}}{\tau} + g'(a_{k+1,i}) + \bar{\psi}_{k+1,i} \\ & = \frac{a_{k+1,j} - a_{k,j}}{\tau} + g'(a_{k+1,j}) + \bar{\psi}_{k+1,j} = c_{k+1}, \text{ for all } i, j \in I_{k+1}. \end{aligned} $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. If $x_{k,i}$ belongs in the interior of $\Omega$ and $x_{k+1,i}$ is in the boundary $\partial\Omega$, we can produce a better competitor by projecting $x_{k+1,i}$ to the convex set $\Omega_\delta$ of points whose distance to the boundary is at least $\delta$.</figcaption></figure>

The constant $c_{k+1}$ above can be obtained by averaging all these quantities, namely

$$c_{k+1} = \frac{1}{|I_{k+1}|} \left( \sum_{j \in I_{k+1}} \bar{\psi}_{k+1,j} + g'(a_{k+1,j}) + \sum_{j \in I_{k+1}} \frac{a_{k+1,j} - a_{k,j}}{\tau} \right) \\ = \frac{1}{|I_{k+1}|} \left( \sum_{j \in I_{k+1}} \bar{\psi}_{k+1,j} + g'(a_{k+1,j}) + \sum_{j \in I_k \setminus I_{k+1}} \frac{a_{k,j}}{\tau} \right),$$

where the second equality comes from the fact that both $a_{k+1,i}$ and $a_{k,i}$ sum to 1.

As a result, by adding the constant $c_{k+1}$ to the Kantorovitch potentials, i.e. $\psi_{k+1,i} \stackrel{\text{def.}}{=} \bar{\psi}_{k+1,i} + c_{k+1}$, they remain optimal we have the equality

$$-\frac{a_{k+1,i} - a_{k,i}}{\tau} = g'(a_{k+1,i}) + \psi_{k+1,i} \text{ for all } i \in I_{k+1},$$

which concludes the proof. $\square$

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

In the previous proposition, one would like to obtain more information on the minimizers to remove the dependence on the normal cones of $\Omega$ and $\Delta_{N-1}$. The next result gives a natural criterion for a point $x_{k+1,i}$ to be in the interior of $\Omega$, so that the normal cone $N_\Omega(x_{k+1,i})$ is null.

**Proposition 3.2.** Let $\Omega \subset \mathbb{R}^d$ be compact and convex, and let $(\varrho_{k+1}^\tau, \mathbf{x}_{k+1}^\tau, \mathbf{a}_{k+1}^\tau)$ be a minimizer of (2.8). Then for all $i = 1, \dots, N$ if $a_{k+1,i} > 0$ and $x_{k,i} \in \text{int}\Omega$, it also holds that $x_{k+1,i} \in \text{int}\Omega$.

Once again we omit the dependence on $\tau$ in this proof since it is fixed.

*Proof.* First notice that from the à priori estimates it holds that $\mathcal{F}(\varrho_{k+1}) < +\infty$ and hence, since $F$ is superlinear, $\varrho_{k+1}$ is uniformly integrable, so that for any $\varepsilon > 0$ there is $\delta' > 0$ such that if

$$E \text{ is a Borelian with } |E| < \delta', \text{ then } \varrho_{k+1}(E) < \varepsilon.$$

Assume by contradiction that $x_{k+1,i} \in \partial\Omega$. In addition, since $x_{k,i} \in \text{int}\Omega$, it holds that $\text{dist}(x_{k,i}, \partial\Omega) > 0$. Defining the family of closed and convex sets

$$\Omega_\delta \stackrel{\text{def.}}{=} \{x \in \Omega : \text{dist}(x, \partial\Omega) \ge \delta\},$$
it holds that for $\delta$ small enough $x_{k,i} \in \Omega_\delta$, and as the Minkowski content of $\partial\Omega$ coincides with its perimeter, we have that $|\Omega \setminus \Omega_\delta| \xrightarrow{\delta \to 0} 0$. So that we can choose $\delta$ small enough so that $\varrho_{k+1}(\Omega \setminus \Omega_\delta) < \varepsilon$.

For such a choice of $\delta$, we set $x_\delta \stackrel{\text{def.}}{=} \text{proj}_{\Omega_\delta}(x_{k+1,i})$, the unique orthogonal projection onto the convex set $\Omega_\delta$, and define a new competitor $\mathbf{x}_\delta$ by replacing the entry $x_{k+1,i}$ by $x_\delta$, see Figure 3.1. As a result, since $x_{k,i} \in \Omega_\delta$, it follows from convexity that

$$ \|x_{\delta} - x_{k,i}\|^2 \le \|x_{k+1,i} - x_{k,i}\|^2. $$

Hence, comparing the energies of this competitor and the minimizer, it follows that

$$ \begin{aligned} & 0 \le W_2^2(\varrho_{k+1}, \mu_{\mathbf{x}_\delta, \mathbf{a}_{k+1}}) - W_2^2(\varrho_{k+1}, \mu_{\mathbf{x}_{k+1}, \mathbf{a}_{k+1}}) \\ & \le \int_{\Omega_{k+1,i} \cap \Omega_\delta} (|x - x_\delta|^2 - |x - x_{k+1,i}|^2) \, d\varrho_{k+1} \\ & + \int_{\Omega_{k+1,i} \setminus \Omega_\delta} (|x - x_\delta|^2 - |x - x_{k+1,i}|^2) \, d\varrho_{k+1}. \end{aligned} $$

To estimate the first term, notice that for any $x \in \Omega_{k+1,i} \cap \Omega_\delta$ we have that

$$ \begin{aligned} |x - x_\delta|^2 - |x - x_{k+1,i}|^2 &= |x_\delta|^2 - |x_{k+1,i}|^2 - 2 \langle x, x_\delta - x_{k+1,i} \rangle \pm 2 \|x_\delta - x_{k+1,i}\|^2 \\ &= |x_\delta|^2 - |x_{k+1,i}|^2 - 2 \langle x_\delta, x_\delta - x_{k+1,i} \rangle \\ & \quad - 2|x_\delta - x_{k+1,i}|^2 - 2 \langle x - x_\delta, x_\delta - x_{k+1,i} \rangle \\ &= - |x_\delta - x_{k+1,i}|^2 - 2 \langle x - x_\delta, x_\delta - x_{k+1,i} \rangle \le -\delta^2, \end{aligned} $$

where in the last inequality we have used that $x_{k+1,i} \in \partial\Omega$ and classical properties of the projection onto closed convex sets.

On the other hand, to estimate the second term, for any $x \in \Omega_{k+1,i} \setminus \Omega_\delta$ it holds that $|x - x_\delta|^2 - |x - x_{k+1,i}|^2 \le \delta^2$. Indeed, for any such $x$, we have that

$$ \begin{aligned} |x - x_\delta|^2 &= |x - x_{k+1,i}|^2 + |x_{k+1,i} - x_\delta|^2 + 2 \langle x - x_{k+1,i}, x_{k+1,i} - x_\delta \rangle \\ &\le |x - x_{k+1,i}|^2 + |x_{k+1,i} - x_\delta|^2 \\ &\le |x - x_{k+1,i}|^2 + \delta^2, \end{aligned} $$

once again the last inequality coming from the properties of projection onto convex sets.

Coming back to the estimation of the energies, we obtain that

$$ \begin{aligned} & 0 \le W_2^2(\varrho_{k+1}, \mu_{\mathbf{x}_\delta, \mathbf{a}_{k+1}}) - W_2^2(\varrho_{k+1}, \mu_{\mathbf{x}_{k+1}, \mathbf{a}_{k+1}}) \\ & \le \delta^2 (-\varrho_{k+1}(\Omega_{k+1,i} \cap \Omega_\delta) + \varrho_{k+1}(\Omega_{k+1,i} \setminus \Omega_\delta)) \\ & \le \delta^2 (-a_{k+1,i} + 2\varrho_{k+1}(\Omega_{k+1,i} \setminus \Omega_\delta),). \end{aligned} $$

Choosing $\delta$ small enough so that $2\varrho_{k+1}(\Omega_{k+1,i} \setminus \Omega_\delta) \le 2\varepsilon < a_{k+1,i}$, the previous construction contradicts the minimality, so that it must hold that $x_{k+1,i} \in \text{int}\Omega$. $\square$

Concerning the behavior of the weights $a_i^\tau$, we use their optimality conditions to show that if $a_i$ is a limit curve as $\tau \to 0$, then it has the property that if at some time $a_i$ reaches 0, then it remains null for any subsequent time instant.

**Proposition 3.3.** Let $\Omega \subset \mathbb{R}^d$ be compact and convex, and let $(\varrho_{k+1}^\tau, \mathbf{x}_{k+1}^\tau, \mathbf{a}_{k+1}^\tau)$ be a minimizer of (2.8). For all $i = 1, \dots, N$ let $t \mapsto a_i(t)$ be a limit trajectory of $a_i^\tau(\cdot)$, then if $a_i(t) = 0$ we have $a_i(s) = 0$ for all $s \ge t$.

*Proof.* Let us first recall the inactive set $I_{k+1} = \{h: 0 < a_{k+1,h} < 1\}$ defined in Proposition 3.1. Given a limit curve $a_i(\cdot)_{i=1,...,N}$, define

$$ t_1 \stackrel{\text{def.}}{=} \min \{t \ge 0 : a_i(t) = 0 \text{ for } i \in \{1, \dots, N\}\} \text{ and } \Lambda_1 \stackrel{\text{def.}}{=} \{i : a_i(t_1) = 0\}, $$

and by induction we define for $n > 1$

$$ t_n \stackrel{\text{def.}}{=} \min \{t \ge 0 : a_i(t) = 0 \text{ for } i \in \{1, \dots, N\} \setminus \Lambda_{n-1}\} \text{ and } \Lambda_n \stackrel{\text{def.}}{=} \{i : a_i(t_n) = 0\}. $$

In particular, we expect that $\Lambda_n \subset \Lambda_{n+1}$.
Let us start by proving that $a_i(t) = 0$ for all times $t \in (t_1, t_2)$ and all indexes $i \in \Lambda_1$. Fix some time $t \in (t_1, t_2)$, hence by definition there is some $\delta$ such that for all $j \notin \Lambda_1$ it holds that $a_j(\cdot) > \delta$ over $(t_1, t)$. In particular, let $k \in \mathbb{N}$ be such that $k\tau < t_1 \le (k+1)\tau$, so that for $\tau$ small enough if holds, from the uniform convergence of $a_i^\tau$ to $a_i$, that

$$a_{k,i}^{\tau} \le C\sqrt{\tau}, \text{ for all } i \in \Lambda_1 \text{ and } a_{n,j}^{\tau} \ge \delta/2, \text{ for all } j \notin \Lambda_1, t_1 \le n\tau \le t.$$

Recall the definition of the set $I_{k+1} \stackrel{\text{def.}}{=} \{i: 0 < a_{k+1,i}^{\tau} < 1\}$, and notice that either $I_{k+1}$ is empty or it has at least two elements since if it where to have only one element, there would either be a single positive $a_{k+1,i}^{\tau}$ but smaller than 1, or there would be an index with mass 1 and another in $I_{k+1}$, either way contradicting the fact that the total mass sums to 1. We will now prove that

$$\sum_{i \in \Lambda_1 \cap I_{k+1}} a_{k+1,i}^{\tau} \le \sum_{i \in \Lambda_1 \cap I_k} a_{k,i}^{\tau} \le CN\sqrt{\tau},$$

with an argument that can be carried by induction for all $n \ge k$ such that $n\tau < t$. The case that $I_{k+1}$ is empty is trivially true, hence we suppose that $|I_{k+1}| \ge 2$.

From the optimality conditions of $a_{k+1,i}^{\tau}$, we have that

$$\sum_{i \in \Lambda_1 \cap I_{k+1}} a_{k+1,i}^{\tau} - a_{k,i}^{\tau} = \tau \left( \sum_{i \in \Lambda_1 \cap I_{k+1}} -g'(a_{k+1,i}^{\tau}) - \psi_{k+1,i}^{\tau} \right).$$

And now we can use the characterization of the Kantorovitch potential $\psi_{k+1}$ in terms of the potential $\bar{\psi}_{k+1}$, which is bounded by a constant depending on $\Omega$. Indeed, from the proof of Proposition 3.1 we have that

$$\sum_{i \in \Lambda_1 \cap I_{k+1}} a_{k+1,i}^{\tau} - a_{k,i}^{\tau} = \tau \left( \sum_{i \in \Lambda_1 \cap I_{k+1}} -g'(a_{k+1,i}^{\tau}) - \bar{\psi}_{k+1,i}^{\tau} \right. \\ \left. + \frac{1}{|I_{k+1}|} \left( \sum_{j \in I_{k+1}} \bar{\psi}_{k+1,j} + g'(a_{k+1,j}^{\tau}) + \sum_{j \in I_k \setminus I_{k+1}} \frac{a_{k,j}^{\tau}}{\tau} \right) \right) \\ \le \tau \left( C_{\Omega} - \left(1 - \frac{1}{|I_{k+1}|}\right) \sum_{i \in \Lambda_1 \cap I_{k+1}} g'(a_{k+1,i}^{\tau}) + \frac{1}{|I_{k+1}|} \sum_{j \in I_{k+1} \setminus \Lambda_1} g'(a_{k+1,j}^{\tau}) \right) \\ + \frac{1}{|I_{k+1}|} \sum_{j \in I_k \setminus I_{k+1}} a_{k,j}^{\tau}.$$

But notice that, for all indexes $j \in \Lambda_1 \setminus I_{k+1}$, the mass $a_{k+1,j}^{\tau}$ is bounded from below away from zero, hence $g'(a_{k+1,j}^{\tau})$ is bounded from above by a finite positive constant. On the other hand, as $\tau \to 0$ we have that $g'(a_{k+1,i}^{\tau}) \to \infty$ for $i \in \Lambda_1 \cap I_{k+1}$. As a result the term being multiplied by $\tau$ in the above inequality is negative for $\tau$ small enough. It then follows that

$$\sum_{i \in \Lambda_1 \cap I_{k+1}} a_{k+1,i}^{\tau} \le \sum_{i \in \Lambda_1 \cap I_{k+1}} a_{k,i}^{\tau} + \sum_{j \in I_k \setminus I_{k+1}} a_{k,j}^{\tau} \\ = \sum_{i \in \Lambda_1 \cap I_{k+1}} a_{k,i}^{\tau} + \sum_{j \in I_k \setminus I_{k+1}} a_{k,j}^{\tau} = \sum_{i \in \Lambda_1 \cap I_k} a_{k,i}^{\tau},$$

where the last equality comes from the fact that $I_k \setminus I_{k+1} \subset \Lambda_1$. Indeed, in our current case $|I_{k+1}| \ge 2$, if $i \in I_k \setminus I_{k+1}$ it must hold that $a_{k+1,i}^{\tau} = 0$ and $0 < a_{k,i}^{\tau} \le C\sqrt{\tau}$, which cannot be the case if $i \in \Lambda_1^c$, since in this complement $a_i^\tau(t)$ is uniformly bounded away from 0 over $(t_1, t)$.

We can iterate this argument, obtaining that

$$\sum_{i \in \Lambda_1 \cap I_{n+1}} a_{n+1,i}^{\tau} < \sum_{i \in \Lambda_1 \cap I_n} a_{n,i}^{\tau}, \text{ while } t_1 < (n+1)\tau < t_2.$$

This implies that for all $t_1 < t < t_2$ we have

$$\sum_{i \in \Lambda_1 \cap I_{n+1}} a_i^{\tau}(t) < \sum_{i \in \Lambda_1 \cap I_n} a_i^{\tau}(t_1) \le |\Lambda_1| C \sqrt{\tau}.$$
From uniform convergence we conclude that $a_i(\cdot) = 0$ for any $i \in \Lambda_1$ over $(t_1, t)$.

Since the curves $t \mapsto a_i(t)$ are continuous, it must hold that $\Lambda_1 \subset \Lambda_2$, and we can repeat the same argument for bigger set of indexes $\Lambda_2$ on the time interval $(t_2, t_3)$. It then follows by the same argument and an induction principle that $a_i(t) = 0$ for any $i \in \Lambda_n$ and $t \ge t_n$. The result follows. $\square$

**3.2. Strong convergence of the pressure.** Recalling that the pressure variable is defined as $P(\varrho) = \varrho F'(\varrho) - F(\varrho)$, and $P(\cdot)$ is a continuous function, we notice that its gradient is given by

$$\nabla P(\varrho) = \varrho \nabla F'(\varrho).$$

As a result, the optimality conditions from Prop. 3.1 yield a uniform estimate for $(P(\varrho^\tau))_{\tau>0}$.

**Lemma 3.4.** There exists a constant $C$, depending only on $\Omega$ and $T$ such that

$$\|P(\varrho^\tau)\|_{L^1([0,T];W^{1,1}(\Omega))} \le C,$$

$$\|\nabla P(\varrho^\tau)\|_{L^1([0,T] \times \Omega)} \le C,$$

for all $\tau > 0$.

*Proof.* The bound on $\|P(\varrho^\tau)\|_{L^1([0,T];L^1(\Omega))}$ is easily obtainable with the a priori energy estimates (2.9). From the optimality conditions on $\varrho_{k+1}^\tau$, we obtain for the stair case interpolation that

$$|\nabla P(\varrho_t^\tau)| \le (|\mathbf{v}_t^\tau| + \operatorname{diam} \Omega)\varrho_t^\tau, \text{ a.e. in } [0, T] \times \Omega,$$

so that

$$\begin{align*}
\int_0^T \|\nabla P(\varrho_t^\tau)\|_{L^1(\Omega)} dt &\le \sum_{k=0}^{N_\tau-1} \int_{k\tau}^{(k+1)\tau} \left( \int_\Omega |\nabla P(\varrho_{k+1}^\tau)| dx \right) dt \\
&\le \sum_{k=0}^{N_\tau-1} \tau \left( \int_\Omega (|\mathbf{v}_{(k+1)}^\tau| + \operatorname{diam} \Omega) \varrho_{k+1}^\tau dx \right) dt \\
&\le \tau N_\tau + \frac{1}{\tau} \sum_{k=0}^{N_\tau-1} \int_\Omega |\mathrm{id} - T_{k+1}^\tau|^2 d\varrho_{k+1}^\tau \\
&\le \tau N_\tau + \frac{1}{\tau} \sum_{k=0}^{N_\tau-1} W_2^2(\varrho_k^\tau, \varrho_{k+1}^\tau) \\
&\le 2T + \mathcal{E}(\varrho_0, \mu_0),
\end{align*}$$

and the result follows. $\square$

In the sequel, we are ready to strengthen the convergence of $\rho^\tau$ to strong convergence in $L^1$ and as a result, prove that $P(\varrho^\tau)$ and $\nabla P(\varrho^\tau)$ also pass to the limit. This will be done, as in [16] and [21], resorting to the following result by Savaré and Rossi [39, Thm. 2].

**Theorem 3.5.** On a Banach space $\mathcal{X}$, let be given

* a normal coercive integrand, $\Phi : \mathcal{X} \to [0, +\infty]$, i.e. lower semi-continuous with relatively compact sublevels in $\mathcal{X}$;

* a pseudo-distance $d : \mathcal{X} \times \mathcal{X} \to [0, +\infty]$, i.e., $d$ is lower semi-continuous, and if $d(\rho, \rho') = 0$ with $\Phi(\rho), \Phi(\rho') < \infty$, then $\rho = \rho'$.

Let $\mathcal{U}$ be a set of measurable functions $u : (0, T) \to \mathcal{X}$. Assuming that

$$\sup_{u \in \mathcal{U}} \int_0^T \Phi(u(t)) dt < \infty \quad \text{and} \quad \limsup_{h \to 0} \sup_{u \in \mathcal{U}} \int_0^{T-h} d(u(t+h), u(t)) dt = 0,$$

then $\mathcal{U}$ contains a subsequence sequence that converges in measure (with respect to $t \in (0, T)$) to a limit $u_\star$.
**Proposition 3.6.** Up to subsequences, the staircase approximation $(\varrho^\tau)_{\tau>0}$ satisfies

$$
\begin{align*}
\varrho^\tau \underset{\tau \to 0}{\longrightarrow} \varrho \text{ in } L^1([0, T] \times \Omega) \text{ and pointwise,} \\
P(\varrho^\tau) \underset{\tau \to 0}{\longrightarrow} P(\varrho) \text{ in } L^1([0, T] \times \Omega) \text{ and pointwise,} \\
\nabla P(\varrho^\tau) \underset{\tau \to 0}{\stackrel{*}{\longrightarrow}} \nabla P(\varrho) \text{ in } \mathcal{M}^d([0, T] \times \Omega).
\end{align*}
$$

*Proof.* Throughout this proof, many different convergence results are proven up to subsequences, so at each step we assume without saying explicitly that we have selected a subsequence for which all the previous convergences hold.

First, we claim that the conditions of Thm. 3.5 hold for the family $(\varrho^\tau, \mathbf{x}^\tau, \mathbf{a}^\tau)_{\tau>0}$ with the choice

$$
d(\xi, \xi') \stackrel{\text{def.}}{=} \begin{cases} d_{W_2, \ell_2}(\xi, \xi'), & \text{if } \xi, \xi' \in \mathcal{P}_{\text{ac}}(\Omega) \times \mathbb{R}^{dN} \times \Delta_{N-1} \\ +\infty, & \text{otherwise} \end{cases}
$$

and

$$
\Phi(\xi) \stackrel{\text{def.}}{=} \begin{cases} \mathcal{F}(\varrho) + \|P(\varrho)\|_{\text{BV}(\Omega)} + \|\mathbf{x}\|^2 + \|\mathbf{a}\|^2, & \text{if } \xi \in \mathcal{P}_{\text{ac}}(\Omega) \times \mathbb{R}^{dN} \times \Delta_{N-1} \\ +\infty, & \text{otherwise,} \end{cases}
$$

as shown in [16, 21]. As a result, there is a subsequence such that $\varrho^\tau \xrightarrow{\tau \to 0} \varrho$ in measure in $[0, T]$, therefore up to a further subsequence, convergence also holds in $L^1([0, T]; L^1(\Omega))$. From Lebesgue's dominated convergence theorem, it also holds that $\varrho^\tau$ converges to $\varrho$ in $L^1([0, T] \times \Omega)$ and pointwise for a.e. $(t, x) \in [0, T] \times \Omega$.

Moving on to the convergence of $P(\varrho^\tau)$, first notice that since $\rho \mapsto P(\rho)$ is continuous, we have a.e. convergence of $P(\varrho^\tau)$ to $P(\varrho)$. In addition, notice that Lemma 2.3 gives a uniform bound

$$
\|P(\varrho^\tau)\|_{L^1([0,T];BV(\Omega))} \leq C,
$$

and thanks a Sobolev embedding, a uniform bound on $\|P(\varrho^\tau)\|_{L^1([0,T];L^{d/(d-1)})}$. A suitable interpolation inequality, proved for instance in [16, Lemma 5.3], gives that

$$
\|P(\varrho^\tau)\|_{L^{(d+1)/d}([0,T]\times\Omega)} \le C,
$$

so that by De la Vallée Poussin's Theorem [24, Thm. 2.29] $(P(\varrho^\tau))_{\tau>0}$ is uniformly integrable in $[0, T] \times \Omega$, since

$$
\sup_{\tau>0} \int_{[0,T]\times\Omega} \gamma(|P(\varrho^\tau(t,x))|) dt dx < +\infty
$$

with the superlinear function $\rho \mapsto \gamma(\rho) = |\rho|^{\frac{d+1}{d}}$. If follows from Vitali's convergence theorem [24, Thm. 2.2.4] that $P(\varrho^\tau)$ converges to $P(\varrho)$ in $L^1([0, T] \times \Omega)$.

In addition, the uniform bound in $L^1([0, T] \times \text{BV}(\Omega))$ along with the strong convergence in $L^1$, allows us to use weak compactness in $\text{BV} [1, \text{Prop. 3.13}]$ to conclude that

$$
\nabla P(\varrho^\tau) \xrightarrow{\tau \to 0} \nabla P(\varrho) \text{ in } \mathcal{M}^d([0, T] \times \Omega),
$$

the result follows.

□

**3.3. The limit PDE-ODE system.** Now we can capitalize on the previous results to conclude with the convergence of the proposed minimizing movement scheme and consequently to the existence of a weak solution to the coupled system (1.4).

**Theorem 3.7.** Assume that $\mathcal{E}(\varrho_0, \mu_0) < +\infty$ and either

• the initial conditions $\mathbf{x}_0 \in (\text{int} \, \Omega)^{\otimes N}$;

• $\Omega$ has a $\mathcal{C}^1$ boundary.

Under these assumptions, the minimizing movement scheme defined in (2.8) admits subsequences converging to weak solutions of (1.4).

*Proof.* Recall that from Proposition 2.4, the pair $(\varrho, E)$ solves the continuity equation. Hence, to obtain the desired limit it suffices to characterize the momentum measure $E$.
To this end, recall that $E^\tau \xrightarrow[\tau \to 0]{} E$ in $\mathcal{M}^d([0, T] \times \Omega)$, where $E^\tau$ is the family of staircase momenta, characterized from the Euler-Lagrange equations (Prop. 3.1) as

$$E^\tau = \mathbf{v}^\tau \varrho^\tau = -\nabla P(\varrho^\tau) - \varrho^\tau \sum_{i=1}^{N} (x - x_{t,i}^\tau) \mathbf{1}_{\Omega_{t,i}^\tau}.$$

In the sequel, consider $\varphi \in \mathcal{C}_c^\infty((0, T) \times \Omega; \mathbb{R}^d)$, we have that

$$\int_{[0,T]\times\Omega} \varphi \cdot dE = \lim_{\tau \to 0} \int_{[0,T]\times\Omega} \varphi \cdot dE^\tau \\ = \lim_{\tau \to 0} \left( \int_{[0,T]\times\Omega} \mathrm{div}\varphi P(\varrho^\tau) dx - \sum_{i=1}^{N} \int_{[0,T]\times\Omega} \varphi \cdot (x - x_{t,i}^\tau) \mathbf{1}_{\Omega_{t,i}^\tau} d\varrho^\tau \right).$$

From the strong convergence of $P(\varrho^\tau)$ to $P(\varrho)$ in $L^1([0, T] \times \Omega)$, the first term has the desired limit. As for the second, recall that $x_{t,i}^\tau \xrightarrow[\tau \to 0]{} x_{t,i}$ uniformly in $[0, T]$. As for the Kantorovitch potentials, we recall as in the proof of Proposition 3.1 that the potentials are unique up to constants and that the potentials used in the optimality conditions for $\mathbf{a}_{k+1}^\tau$ are for the form

$$\psi_{k+1,i}^\tau = \bar{\psi}_{k+1,i}^\tau + c_{k+1}^\tau,$$

where the first term corresponds to another potential bounded by a constant depending on the domain, $|\bar{\psi}_{k+1,i}^\tau| \le C_\Omega$, and $c_{k+1}^\tau$ is a constant that might explode as $\tau \to 0$. We can therefore define the curves

$$\psi_t^\tau \stackrel{\mathrm{def.}}{=} \psi_{k+1}^\tau \quad \text{and} \quad \bar{\psi}_t^\tau \stackrel{\mathrm{def.}}{=} \bar{\psi}_{k+1}^\tau, \quad \text{if } t \in (k\tau, (k+1)\tau].$$

So that, $\bar{\psi}_t^\tau$ is uniformly bounded and converges point-wise in $[0, T]$ to $\bar{\psi}_t$, a bounded Kantorovitch potential associated with $W_2^2(\varrho_t, \mu_t)$. But recall from the Lemma 2.1 that

$$(\psi, \mathbf{x}) \mapsto \mathrm{Lag}_i(\psi, \mathbf{x})$$

is continuous from $\mathbb{R}^d \times \mathbb{R}^d$ to $L^1$, and since Laguerre cells are invariant by the addition of constants, it follows that

$$\Omega_{t,i}^\tau = \mathrm{Lag}_i(\psi_t^\tau, \mathbf{x}_t^\tau) \xrightarrow{\tau \to 0} \mathrm{Lag}_i(\psi_t, \mathbf{x}_t) = \Omega_{t,i} \text{ point-wise in } [0, T].$$

From Egorov's Theorem, for each $\varepsilon > 0$, there a set $I_\varepsilon \subset [0, T]$, with $|[0, T] \setminus I_\varepsilon| \le \varepsilon$, where this point-wise convergence can be strengthened to uniform converge. It then follows that for every $\varepsilon > 0$, we have

$$\limsup_{\tau \to 0} \left| \int_{[0,T]\times\Omega} \varphi \cdot (x - x_{t,i}^\tau) \mathbf{1}_{\Omega_{t,i}^\tau} d\varrho^\tau - \int_{[0,T]\times\Omega} \varphi \cdot (x - x_{t,i}^\tau) \mathbf{1}_{\Omega_{t,i}^\tau} d\varrho^\tau \right| \\ \le \limsup_{\tau \to 0} \left| \int_{I_\varepsilon \times \Omega} \varphi \cdot (x - x_{t,i}^\tau) \mathbf{1}_{\Omega_{t,i}^\tau} d\varrho^\tau - \int_{I_\varepsilon \times \Omega} \varphi \cdot (x - x_{t,i}^\tau) \mathbf{1}_{\Omega_{t,i}^\tau} d\varrho^\tau \right| + C_{\varphi, \Omega, \varrho\varepsilon}.$$

From the strong convergence of $\varrho^\tau$ to $\varrho$ in $L^1([0, T] \times \Omega)$, Prop. 3.6, the uniform convergence of $x_{t,i}^\tau$ and of $\mathbf{1}_{\Omega_{t,i}^\tau}$, the second $\limsup$ above is null. As a consequence, we conclude that

$$E = -\nabla P(\varrho) - \varrho \sum_{i=1}^{N} (x - x_{t,i}^\tau) \mathbf{1}_{\Omega_{t,i}^\tau}.$$

Moving on to the dynamics of the atoms, recall that from Proposition 2.4 it holds that $\bar{\mathbf{x}}^\tau$ converges to $\mathbf{x}$ weakly in $H^1([0, T])$ both $\bar{\mathbf{x}}^\tau, \mathbf{x}^\tau$ converge uniformly (and hence also strongly in $L^2([0, T])$) to $\mathbf{x}$. To characterize the limit dynamics of the $\mathbf{x}$ variables, consider $\varphi \in \mathcal{C}_c^\infty((0, T))$ and we treat two cases: either

* the initial conditions $\mathbf{x}_0 \in (\mathrm{int}\,\Omega)^{\otimes N}$;

* $\Omega$ has a $\mathcal{C}^1$ boundary.
For simplicity, let us discuss the first one. In this case, Prop. 3.2 ensures that the entire sequence $(x_k)_k$ is contained in the interior of $\Omega$ and therefore the boundary effects of the normal cone does not intervene in the Euler-Lagrange equations from Prop. 3.1, so that we have

$$ 
\begin{align*}
-\int_0^T \varphi'(t)x_{t,i}dt &= \lim_{\tau \to 0} -\int_0^T \varphi'(t)\bar{x}_{t,i}^\tau dt = \lim_{\tau \to 0} \int_0^T \varphi(t)\dot{\bar{x}}_{t,i}^\tau dt \\
&= \lim_{\tau \to 0} -\int_0^T \varphi(t) \left( a_{k+1,i}^\tau \cdot x_{k+1,i}^\tau - \int_{\Omega_{k+1,i}^\tau} x d\varrho_{k+1}^\tau \right) dt \\
&= -\int_0^T \varphi(t) \left( a_{t,i} \cdot x_{t,i} - \int_{\Omega_{t,i}} x d\varrho_t \right) dt,
\end{align*}
$$

where the last limit hold, as before, because of the uniform convergence in $t$ of the Laguerre cells and the convergence of $\varrho^\tau$ to $\varrho_t$. This gives the desired characterization of the flow of $x_t$. The calculations are analogous in the second case, with the addition of the normal vector $n_\Omega(\bar{x}_{k+1,i}^\tau)$, which converges to a vector in $N_\Omega(x_{t,i})$ since the boundary is smooth.

The dynamics of the weights can be deduced in a similar way to the convergence of the atoms. Also from Proposition 2.4, we know that $\bar{a}^\tau$ converges to $\mathbf{a}$ weakly in $H^1([0, T])$ and both $\bar{a}^\tau, \mathbf{a}^\tau$ converge uniformly to $\mathbf{a}$. On the other hand, we must take into account the singularity at zero presented by the optimality conditions of the $a_i$ variables given by the term $g'(a_i)$, which is compensated with the Kantorovitch potentials. The issue as we have seen is that the Kantorovitch potentials consists of a bounded potential added by a constant that might diverge to compensate for the singularity in any single $g'(a_i(t))$ that goes to infinity.

As in the proof of Proposition 3.3, define the quantities

$$ t_i \stackrel{\text{def.}}{=} \inf \{t \ge 0 : a_i(t) = 0\} \quad \text{for } i = 1, \dots, N, $$

and for a fixed $i$, take some other index $j$ such that $t_j \ge t_i$. If such $j$ does not exist it means that $a_i \equiv 1$ and $t_1 = +\infty$. Hence, to characterize the dynamics of $a_i$ it suffices to consider test functions in $\varphi \in C_c^\infty((0, t_i))$, since after this time we know that $a_i$ vanishes. As a result, for such test function we obtain that $\inf_{\text{supp } \varphi} a_i$, $\inf_{\text{supp } \varphi} a_j > 0$ and from the uniform convergence the same holds for $a_i^\tau$ and $a_j^\tau$. As a result we obtain that

$$ 
\begin{align*}
\int_0^T (a_i(t) - a_j(t)) \dot{\varphi}(t) dt &= \lim_{\tau \to 0} \int_0^T (a_i^\tau(t) - a_j^\tau(t)) \dot{\varphi}(t) dt \\
&= \lim_{\tau \to 0} - \int_0^T (\dot{a}_i^\tau(t) - \dot{a}_j^\tau(t)) \varphi(t) dt \\
&= \lim_{\tau \to 0} \int_0^T \varphi(t) (g'(a_i^\tau(t)) - \psi_{t,i}^\tau - g'(a_j^\tau(t)) + \psi_{t,j}^\tau) dt \\
&= \lim_{\tau \to 0} \int_0^T \varphi(t) (g'(a_i^\tau(t)) - \bar{\psi}_{t,i}^\tau - g'(a_j^\tau(t)) + \bar{\psi}_{t,j}^\tau) dt \\
&= \int_0^T \varphi(t) (g'(a_i(t)) - \bar{\psi}_{t,i} - g'(a_j(t)) + \bar{\psi}_{t,j}) dt,
\end{align*}
$$

since $g \in C_{\text{loc}}^1((0, 1])$, where $a_i^\tau, a_j^\tau$ remains bounded from below in the support of $\varphi$, and the Kantorovitch potentials $\bar{\psi}_t$ being uniformly bounded and Lipschitz, converge to the Kantorovitch potential $\bar{\psi}_t$.

We then conclude that for any pair $i, j$ such that $t_j \ge t_i$ their corresponding dynamics satisfy

$$ \dot{a}_i(t) + g'(a_i(t)) - \bar{\psi}_{t,i} = \dot{a}_j(t) + g'(a_j(t)) - \bar{\psi}_{t,j} = c_t, $$

where $c_t \in L^2((0, t))$ for $t < t_i$. We can now do as in Proposition 3.1 and define a new Kantorovitch potential $\psi_t = \bar{\psi}_t + c_t$, in such a way that

$$ \dot{a}_i(t) = -g'(a_i(t)) + \psi_{t,i} \text{ in } L^2((0, t_i)). $$
In addition, since $a_h(t) = 0$ for $t \ge t_h$, we obtain that

$$0 = \frac{d}{dt} \sum_{i \sim a_i(t) > 0} a_i(t) = \sum_{i \sim a_i(t) > 0} -g'(a_i(t)) + \psi_{t,i},$$

which finishes the desired characterization of the dynamics and the result follows. □

**3.4. Strong $L^2H^1$ convergence in Porous Medium case.** In the previous section, as in seminal papers that introduced the JKO scheme [27, 37], the convergence of the minimizing movement scheme was proven in the $L^1$ topology. This is natural when working with probability measures, but certainly the limit PDE (e.g. heat, Fokker-Plack, Porous Medium equation) enjoys much more regularity, which suggests that convergence in stronger topologies also holds. In recent papers [22, 41] the authors show that the JKO scheme of the Boltzman entropy functional

$$(3.4) \qquad \mathcal{F}(\varrho) = \mathcal{H}(\varrho) \stackrel{\text{def.}}{=} \begin{cases} \int_{\Omega} \varrho \log \varrho dx, & \text{if } \varrho \ll \mathcal{L}^d \llcorner \Omega, \\ +\infty, & \text{otherwise,} \end{cases}$$

(resp. with a potential energy term), converges to the unique solution of the heat (resp. Fokker-Plank) equation in the strong $L^2H^2$ topology.

In the case of the PDE from (1.4), the advection term coming from the semi-discrete transport is not sufficiently smooth to adapt the arguments from [22, 41], so that we can expect $L^2H^1$ convergence instead. In fact, we shall prove convergence of the staircase interpolation of the pressure variable $P(\varrho^\tau)$ to $P(\varrho)$ in the strong topology of $L^2H^1$ when the internal energy is given by a Porous Medium term of the form

$$(3.5) \qquad \mathcal{F}_m(\varrho) \stackrel{\text{def.}}{=} \begin{cases} \frac{1}{m-1} \int_{\Omega} \varrho^m dx, & \text{if } \varrho \ll \mathcal{L}^d \llcorner \Omega, \\ +\infty, & \text{otherwise,} \end{cases}$$

obtained for $F(\rho) = F_m(\rho) = \frac{1}{m-1}\rho^m$, and yielding a pressure $P_m(\rho) = \rho^m$.

Following [22, 41], the strategy of the proof will consist on first proving $L^p$ estimates for iterates of the JKO scheme, by bounding their $L^p$ norms by the one of the previous estimate multiplied by a factor that remains summable as $\tau \to 0$. This allows to pass to the limit as $p \to \infty$ to obtain also $L^\infty$ bounds whenever the initial condition is bounded. Our argument consists of exploiting the displacement convexity of $\rho \to \|\rho\|_p^p$ as in [22] in order to deal with the Porous Medium term, which also gives the desired bounds for the gradient. The strong convergence then will follow with an adaptation of the arguments in [41]. The same convergence can be obtained for linear diffusion, but we have chosen to focus on the Porous Medium case for a more concise presentation.

**Theorem 3.8.** Assume that $\Omega$ be a convex and bounded domain of $\mathbb{R}^d$ with Lipschitz boundary. Let $(\varrho_k^\tau, \mathbf{x}_k^\tau, \mathbf{a}_k^\tau)_{k \in \mathbb{N}}$ be a sequence obtained via the JKO scheme, with the internal energy given by the Boltzmann entropy (3.4) or the Porous Medium term $\mathcal{F}_m$ from (3.5). Then the following estimates hold

(1) for all $1 < p < +\infty$, if $\varrho_0 \in L^p(\Omega)$ it holds that for all $k \ge 1$ that

$$\|\varrho_k^\tau\|_{L^p(\Omega)} \le (1 - \tau d(p-1))^{-1/p} \|\varrho_{k-1}^\tau\|_{L^p(\Omega)}.$$

If $\varrho_0 \in L^\infty(\Omega)$, taking the limit as $p \to \infty$ we obtain

$$\|\varrho_k^\tau\|_{L^\infty(\Omega)} \le e^{\tau d} \|\varrho_{k-1}^\tau\|_{L^\infty(\Omega)}.$$

(2) For $\varrho_0 \in L^{m+1}(\Omega)$, there exists a constant $C$ depending on $m, T, \Omega$ and $\varrho_0$ such that

$$\int_0^T \int_{\Omega} |\nabla P_m(\varrho^\tau)|^2 dx dt \le C,$$

for all $\tau$ small enough.

In addition, if $\varrho_0 \in L^{m+1}(\Omega)$, up to subsequences the pressures $P_m(\varrho^\tau)$ converge in the strong topology of $L^2H^1$ to $P(\varrho)$, where $(\varrho, \mathbf{x}, \mathbf{a})$ is a solution of (1.4).
These estimates can be proved with the same techniques for the classical JKO scheme for the Porous Medium, without the semi-discrete transport term and with a Lipschitz potential independent of the $\rho$. In this case, such estimates imply that any solution of the PME obtained with the JKO scheme is an energy solution, see [42, Section 5.3.2] and [4, Def. 7.1], defined as follows

**Definition 3.9.** Let $m > 1$, and $V \in W^{1,\infty}(\Omega)$ be a vector field, which is Lipschitz in space for a.e. $t$ and for all $\rho \in \mathcal{P}(\Omega)$. Then $\rho$ is an energy solution of the Porous Medium equation with advection

$$ (PME_m) \qquad \partial_t \rho = \Delta \rho^m + \operatorname{div}(\rho \nabla V), $$

if it is a solution in the sense of distributions such that

$$ (1) \quad \rho \in \mathcal{C}([0, T]; L^1(\Omega)) \cap L^\infty([0, T] \times \Omega); $$

$$ (2) \quad P_m(\rho) \in L^2([0, T]; H^1(\Omega)). $$

Although the PME might admit more than one weak solution, there is a unique energy solution, see [42, Section 5.3.2] and [4]. This discussion gives the following

**Corollary 3.10.** Let $\Omega$ be a convex and bounded domain of $\mathbb{R}^d$ with Lipschitz boundary. Given a potential $V \in W^{1,\infty}(\Omega)$ and an initial condition $\rho_0 \in L^{m+1}(\Omega) \cap \mathcal{P}(\Omega)$. Let $(\rho^\tau)_{\tau>0}$ be the staircase interpolation of JKO scheme referent to the Porous Medium energy with advection

$$ \rho \mapsto \frac{1}{m-1} \int_{\Omega} \rho^m dx + \int_{\Omega} V \rho dx, $$

then the pressures $(P_m(\rho^\tau))_{\tau>0}$ converge in the strong topology of $L^2([0, T]; H^1(\Omega))$ to the unique energy solution of the PME with advection introduced by the potential $V$.

The $L^2 H^2$ convergence for the Porous Medium case remains open, even with a smooth advection term. The main ingredient to prove such strong convergence in the Fokker-Planck case is the 5-gradients inequality introduced in [20]. In [41, Lemma 5.1], the inequality is derived with a remainder term which can be controlled and then exploited to obtain the stronger convergence.

The major difficulty of the case with a semi-discrete transportation term is that it introduces to the dynamics a vector field with a singular Laplacian, which appear in the computations even for the discrete-time scheme. Hence, our first step in the following Lemma is to compute this Laplacian explicitly.

**Lemma 3.11.** Given $(\rho, \mathbf{x}, \mathbf{a})$, define the potential

$$ (3.6) \qquad \Phi[\rho, \mathbf{x}, \mathbf{a}](x) \stackrel{\text{def.}}{=} \sum_{i=1}^{N} \left( \frac{1}{2} |x - x_i|^2 - \psi_i \right) \mathbf{1}_{\Omega_i}, $$

where $(\Omega_i)_{i=1}^N$ corresponds to the optimal Laguerre tessellation associated to the transportation of $\rho$ to $\mu_{\mathbf{x},\mathbf{a}}$ and $(\psi_i)_{i=1}^N$ is the corresponding Kantorovitch potential.

The gradient vector field of $\Phi[\rho, \mathbf{x}, \mathbf{a}]$ is given by

$$ (3.7) \qquad \nabla \Phi[\rho, \mathbf{x}, \mathbf{a}](x) = \sum_{i=1}^{N} (x - x_i) \mathbf{1}_{\Omega_i}, $$

and its Laplacian is the following measure

$$ (3.8) \qquad \Delta\Phi[\rho, \mathbf{x}, \mathbf{a}] = d + \sum_{j<i} (x_i - x_j) \cdot n_{\Omega_i} \mathcal{H}^{d-1} \llcorner \partial\Omega_j \cap \partial\Omega_i. $$

In particular, it holds that $\Delta\Phi[\rho, \mathbf{x}, \mathbf{a}] \le d$.

*Proof.* For simplicity of notation we shall refer to $\Phi[\rho, \mathbf{x}, \mathbf{a}]$ as $\Phi$ omitting the arguments. Notice that from the definition of Laguerre tessellations, the potential defined in (3.6) can be rewritten as

$$ \Phi[\rho, \mathbf{x}, \mathbf{a}](x) = \min_{i=1,\dots,N} \frac{1}{2} |x - x_i|^2 - \psi_i. $$

As an optimal Kantorovitch potential, it is Lipschitz continuous, and therefore it is a.e. differentiable. As a result, by the envelope theorem its gradient must coincide with the expression in (3.7).
We must only compute its Laplacian in the sense of distributions, for this consider some $\varphi \in C^1(\Omega)$ and compute

$$ \int_{\Omega} \varphi d(\Delta \Phi) dx = - \int_{\Omega} \nabla \varphi \cdot \nabla \Phi dx + \int_{\partial \Omega} \varphi \nabla \Phi \cdot n_{\Omega} d\mathcal{H}^{d-1} \\ = - \sum_{i=1}^{N} \int_{\Omega_i} \nabla \varphi \cdot (x - x_i) dx + \sum_{i=1}^{N} \int_{\partial \Omega_i \cap \partial \Omega} \varphi(x - x_i) \cdot n_{\Omega} d\mathcal{H}^{d-1}. $$

Each integral in the first sum above can be developed as

$$ \int_{\Omega_i} \nabla \varphi \cdot (x - x_i) dx = - d \int_{\Omega_i} \varphi dx + \int_{\partial \Omega_i \cap \partial \Omega} \varphi(x - x_i) \cdot n_{\Omega} d\mathcal{H}^{d-1} \\ + \int_{\partial \Omega_i \setminus \partial \Omega} \varphi(x - x_i) \cdot n_{\Omega_i} d\mathcal{H}^{d-1}, $$

where we have used the fact that $n_{\Omega_i}$ coincides with $n_{\Omega}$ at $\partial \Omega$. Define also $\Sigma_{ij} \stackrel{\text{def.}}{=} \partial \Omega_i \cap \partial \Omega_j$ for $i \neq j$, and set the convention that $n_{\Sigma_{ij}} = n_{\Omega_i} = -n_{\Omega_j}$. Since the intersections $\partial \Omega_i \cap \partial \Omega_j \cap \partial \Omega$ are $\mathcal{H}^{d-1}$ negligible, summing all these contributions we obtain the following expression for the Laplacian

$$ \Delta \Phi = d + \sum_{j<i} (x_i - x_j) \cdot n_{\Omega_i} \mathcal{H}^{d-1} \llcorner \partial \Omega_j \cap \partial \Omega_i. $$

However, notice that for any point $x \in \partial \Omega_j \cap \partial \Omega_i$, we have from convexity that

$$ \begin{aligned} (x_i - x_j) \cdot n_{\Omega_i}(x) &= (x_i - x) \cdot n_{\Omega_i}(x) - (x_j - x) \cdot n_{\Omega_i}(x) \\ &= (x_i - x) \cdot n_{\Omega_i}(x) + (x_j - x) \cdot n_{\Omega_j}(x) \le 0, \end{aligned} $$

and we conclude that $\Delta \Phi \le d$. $\square$

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

Now we pass to the question of $L^p, L^\infty$ estimates for iterates of the JKO scheme and to $H^1$ estimates for the pressure gradient. The arguments are an adaptation of the ones in [22,41], relying on the flow interchange technique.

*Proof of items (1) and (2) of Theorem 3.8:* To simplify notation we drop the dependence on $\tau$ from the iterates of the JKO scheme. As in [22], the strategy of the proof consists of combining Lemma 3.11 above with the geodesic convexity of the functional

$$ \varrho \mapsto \mathcal{F}_p(\varrho) = \frac{1}{p-1} \int_{\Omega} \varrho^p dx. $$

Indeed, if $P_p(\varrho), P_p(\eta) \in W^{1,1}(\Omega)$ and $\varrho_t \stackrel{\text{def.}}{=} ((1-t)\text{id} + tT)_\sharp \varrho$ is the geodesic between them, it holds, for instance from [2, Lemma 10.4.4.], that

$$ \left. \frac{d}{dt} \right|_{t=0^+} \mathcal{F}_p(\varrho_t) = \int_{\Omega} \nabla P_p(\varrho) \cdot (T - \text{id}) dx - \int_{\partial \Omega} P_p(\varrho) (T - \text{id}) \cdot n_{\Omega} d\mathcal{H}^{d-1} \\ \ge \int_{\Omega} \nabla P_p(\varrho) \cdot (T - \text{id}) dx. $$

As in the notation of Lemma 3.11, we set $\Phi_{k+1} \stackrel{\text{def.}}{=} \Phi[\varrho_{k+1}, \mathbf{x}_{k+1}, \mathbf{a}_{k+1}]$ and let $T_{k+1} = \text{id} - \nabla \varphi_{k+1}$ denote the optimal transportation map from $\varrho_{k+1}$ to $\varrho_k$. The optimality conditions for the sequence obtained from the JKO scheme tells us that

$$ m \varrho_{k+1}^{m-1} = -\Phi_{k+1} - \frac{\varphi_{k+1}}{\tau} \quad \text{in} \ \{\varrho_{k+1} > 0\}. $$

As a result, we conclude that $\varrho_{k+1}^{m-1}$ is Lipschitz continuous and hence $\varrho_{k+1}$ is bounded. Hence for all $p > 1$ the pressure $P_p(\varrho_{k+1}) \in L^1(\Omega)$ and

$$ \begin{aligned} \nabla P_p(\varrho_{k+1}) &= p \frac{m-1}{m} \varrho_{k+1}^{p-m+1} \nabla F'_m(\varrho_{k+1}) \\ &= -p \frac{m-1}{m} \varrho_{k+1}^{p-m+1} \left( \nabla \Phi_{k+1} + \frac{\nabla \varphi_{k+1}}{\tau} \right) \in L^\infty(\Omega). \end{aligned} $$
It follows that $P_p(\Omega) \in W^{1,\infty}(\Omega)$ for all $p > 1$, and the geodesic convexity of $\mathcal{F}_p$ gives us that

$$ 
\begin{aligned}
\frac{1}{p-1} \int_{\Omega} \varrho_k^p dx - \frac{1}{p-1} \int_{\Omega} \varrho_{k+1}^p dx &\geq \int_{\Omega} \nabla P_p(\varrho_{k+1}) \cdot (T_{k+1} - \mathrm{id})dx \\
&= p \int_{\Omega} \varrho_{k+1}^{p-1} \nabla \varrho_{k+1} \cdot (-\nabla \varphi_{k+1})dx \\
&= \tau p \int_{\Omega} \varrho_{k+1}^{p-2} \nabla \varrho_{k+1} \cdot (\nabla P_m(\varrho_{k+1}) + \varrho_{k+1} \nabla \Phi_{k+1}) dx
\end{aligned}
$$

And we obtain the key estimate

$$ (3.9) \quad \int_{\Omega} \varrho_k^p dx - \int_{\Omega} \varrho_{k+1}^p dx \geq \tau(p-1)p \left( m \int_{\Omega} \varrho_{k+1}^{p+m-3} |\nabla \varrho_{k+1}|^2 dx + \int_{\Omega} \nabla (\varrho_{k+1})^p \cdot \nabla \Phi_{k+1} dx \right) $$

To obtain the $L^p$ estimates, we bound from below the first term in the RHS of (3.9) by 0 and use integration by parts in the second term in order to obtain

$$ 
\begin{aligned}
\int_{\Omega} \varrho_k^p dx - \int_{\Omega} \varrho_{k+1}^p dx &\geq \tau(p-1) \left( \int_{\Omega} \varrho_{k+1}^p (-\Delta\Phi_{k+1}) dx + \int_{\partial\Omega} \varrho_{k+1}^p \nabla\Phi_{k+1} \cdot n_{\Omega} d\mathcal{H}^{d-1} \right) \\
&\geq -\tau(p-1)d \int_{\Omega} \varrho_{k+1}^p dx,
\end{aligned}
$$

where we have used the fact that $-\Delta\Phi_{k+1} \geq -d$ from Lemma 3.11 and the boundary terms are non-negative from the convexity of the domain as done many times in the proof of Lemma 3.11. This gives the desired estimate from item (1), and passing to the limit as $p \to \infty$, we obtain the $L^\infty$ bounds from item (2).

Moving on to the pressure gradient estimates, assuming that $\varrho_0 \in L^{m+1}(\Omega)$ and considering $p = m + 1$ in (3.9) we obtain

$$ 
\begin{aligned}
\int_{\Omega} \varrho_k^{m+1} dx - \int_{\Omega} \varrho_{k+1}^{m+1} dx &\geq \tau m \left( m \int_{\Omega} \varrho_{k+1}^{2m-2} |\nabla \varrho_{k+1}|^2 dx + \int_{\Omega} \nabla (\varrho_{k+1}^{m+1}) \cdot \nabla \Phi_{k+1} dx \right) \\
&= \tau \left( m^2 \int_{\Omega} |\nabla P_m(\varrho_{k+1})|^2 dx + (m+1)m \int_{\Omega} \nabla P_m(\varrho_{k+1}) \cdot \nabla \Phi_{k+1} \varrho_{k+1} dx \right) \\
&\geq \tau \left( \frac{m^2}{2} \int_{\Omega} |\nabla P_m(\varrho_{k+1})|^2 dx - \frac{m^2}{2} \int_{\Omega} |\nabla \Phi_{k+1}|^2 \varrho_{k+1}^2 dx \right),
\end{aligned}
$$

where the last step comes from the $\varepsilon$-Young's inequality, with $\varepsilon = (m+1)/m$. Finally, summing both sides of the above inequality over $k$, the LHS telescopes and and recalling that $|\nabla\Phi_{k+1}| \leq \text{diam}(\Omega)$ and that $\varrho^\tau = \varrho_{k+1}$ over the interval $(k\tau, (k+1)\tau)$ we obtain

$$ 
\begin{aligned}
\int_0^T \int_{\Omega} |\nabla P(\varrho^\tau)|^2 dx dt &= \sum_k \tau \int_{\Omega} |\nabla P(\varrho_{k+1})|^2 dx \\
&\leq \sum_k \tau \int_{\Omega} |\nabla \Phi_{k+1}|^2 \varrho_{k+1}^2 dx + \frac{2}{m} \int_{\Omega} \varrho_0^{m+1} dx \leq C_{m, \Omega, \varrho_0}.
\end{aligned}
$$

Finally, to get an $L^2 H^1$ bound on $P_m(\varrho_{k+1})$, notice that since $\varrho_0 \in L^m(\Omega)$, from item (1), this remains true for all subsequent $\varrho_k$, for $\tau$ small enough. So the mean of $P_m(\varrho_{k+1})$ remains uniformly bounded and Poincaré-Wirtinger inequality gives

$$ 
\begin{aligned}
\| P_m(\varrho_k) \|_{L^2(\Omega)} &\leq \left\| P_m(\varrho_k) - \int_{\Omega} P_m(\varrho_k) \right\|_{L^2(\Omega)} + \int_{\Omega} P_m(\varrho_k) \\
&\leq C \| \nabla P_m(\varrho_k) \|_{L^2(\Omega)} + \int_{\Omega} P_m(\varrho_k).
\end{aligned}
$$

It follows that $P_m(\varrho^\tau)$ remains uniformly bounded in $\tau$ inside $L^2([0, T]; H^1(\Omega))$. $\square$

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

Consider now a subsequence, not relabelled, of $(\varrho^\tau, \mathbf{x}^\tau, \mathbf{a}^\tau)_{\tau>0}$ for which the convergence to the limiting equation (1.4) holds from Thm. 3.7. From the previous reasoning $(\varrho^\tau)_{\tau>0}$ is bounded in $L^2([0, T]; H^1(\Omega))$, so we conclude that it converges in the weak topology of this Hilbert space to $\varrho$. Indeed, since Thm. 3.7 implies $L^1([0, T]; \Omega)$ convergence of this fixed subsequence of $(\varrho^\tau)_{\tau>0}$ to some solution, we can extract a further convergent subsequence in $L^2([0, T]; H^1(\Omega))$. This limit
must then be the same solution, otherwise we would have a contradiction with the $L^1$ convergence. From the Urysohn property of the weak $L^2([0, T]; H^1(\Omega))$-convergence, it follows that the original sequence must also converge weakly in $L^2([0, T]; H^1(\Omega))$.

Since it is a Hilbert space, to obtain strong convergence it suffices to prove convergence of the norms. This is done with the interpolation introduced in [41, Section 4] given by

$$ (3.10) \qquad \varrho_t^{\tau, \varepsilon} \stackrel{\text{def.}}{=} \begin{cases} \varrho_{k+1}^{\tau}, & t \in (k\tau, k\tau + (1-\varepsilon)\tau], \\ \varrho_{k+1}^{\tau} \frac{(k+1)\tau-t}{\varepsilon\tau} + \varrho_{k+2}^{\tau} \frac{t-(k+1)\tau+\varepsilon\tau}{\varepsilon\tau}, & t \in (k\tau + (1-\varepsilon)\tau, (k+1)\tau], \\ \varrho_{\lceil \frac{T}{\tau} \rceil}^{\tau}, & t \in (\lceil \frac{T}{\tau} \rceil, T]. \end{cases} $$

The advantage of this new family is that it is well adapted to apply a suitable version of the Aubin-Lions-Simon compactness theorem in $L^p([0, T]; B)$, where $B$ is a Banach space, while still being uniformly close to $\varrho^\tau$ in the $L^2([0, T], H^1(\Omega))$ topology.

Since the application of this compactness theorem depends only on the construction of the interpolation, we summarize in the following Lemma the convergence properties obtained in [41, Prop. 4.2 and Cor. 4.4].

**Lemma 3.12.** Let $(\varrho^\tau)_{\tau>0}$ be a subsequence of the staircase interpolation of the JKO scheme which is strongly convergent in $L^1([0, T] \times \Omega)$, as in Thm. 3.7, to $\varrho$. Then both $(\varrho^\tau)_{\tau>0}$ and $(P_m(\varrho^\tau))_{\tau>0}$ converge to $\varrho$ and $P_m(\varrho)$, respectively, strongly in $L^2([0, T]; L^2(\Omega))$.

With these elements we can finish the proof the strong $L^2([0, T]; H^1(\Omega))$ convergence, whose proof structure is strongly inspired by [41, Thm. 4.5].

*Proof of Thm. 3.8.* Let $(\varrho, \mathbf{x}, \mathbf{a})_{t \in [0, T]}$ be a solution of the coupled system obtained as the limit of a staircase interpolation $(\varrho^\tau, \mathbf{x}^\tau, \mathbf{a}^\tau)_{\tau>0}$. To facilitate notation, recall the definition of the potential $\Phi[\varrho, \mathbf{x}, \mathbf{a}]$ introduced in (3.6) and define

$$ \begin{aligned} \Phi[\varrho] &\stackrel{\text{def.}}{=} \Phi[\varrho, \mathbf{x}, \mathbf{a}], \\ \Phi^\tau &\stackrel{\text{def.}}{=} \Phi[\varrho^\tau, \mathbf{x}^\tau, \mathbf{a}^\tau], \\ \Phi_{k+1} &\stackrel{\text{def.}}{=} \Phi[\varrho_{k+1}^\tau, \mathbf{x}_{k+1}^\tau, \mathbf{a}_{k+1}^\tau]. \end{aligned} $$

In Thm 3.7, we show that $\varrho$ is a solution in a weak sense to the equation

$$ (3.11) \qquad \partial_t \varrho = \nabla P_m(\varrho) + \operatorname{div}(\varrho \nabla \Phi[\varrho]), $$

with non-flux boundary conditions and such that $P_m(\varrho) \in L^1([0, T]; \text{BV}(\Omega))$. In the first part of the proof of the current theorem, we actually show that $P_m(\varrho) \in L^2([0, T]; H^1(\Omega))$. As a result, we can use $P_m(\varrho)$ as a test function for (3.11) giving that

$$ \begin{aligned} & - \int_0^T \int_{\Omega} (|\nabla P_m(\varrho)|^2 + \varrho \nabla P_m(\varrho) \cdot \nabla \Phi[\varrho]) \,dx\,dt = \int_0^T \int_{\Omega} P_m(\varrho) \partial_t \varrho \,dx\,dt \\ & = \int_0^T \frac{d}{dt} \left( \frac{1}{m+1} \int_{\Omega} \varrho_t^{m+1} \,dx \right) \,dt = \frac{1}{m+1} \int_{\Omega} \varrho_T^{m+1} \,dx - \frac{1}{m+1} \int_{\Omega} \varrho_0^{m+1} \,dx. \end{aligned} $$

In other words we have the energy relation

$$ (3.12) \qquad \int_{\Omega} \varrho_T^{m+1} \,dx - \int_{\Omega} \varrho_0^{m+1} \,dx = -(m+1) \int_0^T \int_{\Omega} (|\nabla P_m(\varrho)|^2 + \varrho \nabla P_m(\varrho) \cdot \nabla \Phi[\varrho]) \,dx\,dt $$

On the other hand, using the case $p = m+1$ in the proof of the first part of the present theorem, for instance in equation (3.9), it holds for the discrete iterates of the JKO scheme that

$$ \begin{aligned} & \int_{\Omega} (\varrho_k^{\tau})^{m+1} \,dx - \int_{\Omega} (\varrho_{k+1}^{\tau})^{m+1} \,dx \\ & \geq \tau \left( \int_{\Omega} |\nabla P_m(\varrho_{k+1}^{\tau})|^2 \,dx + \int_{\Omega} \varrho_{k+1}^{\tau} \nabla P_m(\varrho_{k+1}^{\tau}) \cdot \nabla \Phi_{k+1} \,dx \right). \end{aligned} $$

Hence, summing over $k$ the LHS telescopes, and we obtain that

$$ (3.13) \qquad \int_{\Omega} (\varrho_0^\tau)^{m+1} \,dx - \int_{\Omega} (\varrho_T^\tau)^{m+1} \,dx \geq \int_0^T \int_{\Omega} (|\nabla P_m(\varrho^\tau)|^2 + \varrho^\tau \nabla P_m(\varrho^\tau) \cdot \nabla \Phi^\tau) \,dx\,dt. $$
From the pointwise convergence of $\varrho_t^\tau$ to $\varrho_t$ in Wasserstein topology, for every time $t$, and the weak convergence of $P_m(\varrho^\tau)$ we obtain that

$$ (3.14) \qquad \int_{\Omega} (\varrho_T)^{m+1} dx \le \liminf_{\tau \to 0^+} \int_{\Omega} (\varrho_T^\tau)^{m+1} dx, \\ \int_0^T \int_{\Omega} |\nabla P_m(\varrho)|^2 dx dt \le \liminf_{\tau \to 0^+} \int_0^T \int_{\Omega} |\nabla P_m(\varrho^\tau)|^2 dx dt $$

We can rewrite the following term as

$$ \begin{aligned} \lim_{\tau \to 0} \int_0^T \int_{\Omega} \varrho^\tau \nabla P_m(\varrho^\tau) \cdot \nabla \Phi^\tau dx dt &= \sum_{i=1}^N \lim_{\tau \to 0} \int_0^T \int_{\Omega} \mathbb{1}_{\Omega_i^\tau(t)} \varrho^\tau \nabla P_m(\varrho^\tau) \cdot (x - x_i^\tau) dx dt \\ &= \sum_{i=1}^N \int_0^T \int_{\Omega} \mathbb{1}_{\Omega_i(t)} \varrho \nabla P_m(\varrho) \cdot (x - x_i) dx dt \\ &= \int_0^T \int_{\Omega} \varrho \nabla P_m(\varrho) \cdot \nabla \Phi[\varrho] dx dt, \end{aligned} $$

where the limit is computed with weak-strong convergence, since $\mathbb{1}_{\Omega_i^\tau} \xrightarrow{L^1([0,T] \times \Omega)} \mathbb{1}_{\Omega_i}$ for all $i = 1, \dots, N$, so that $\varrho^\tau \mathbb{1}_{\Omega_i^\tau}$ converges strongly in $L^2([0,T]; L^2(\Omega))$ and $\nabla P_m(\varrho^\tau) \cdot (x - x_i^\tau)$ converges weakly in the same space.

As a result, going back to (3.13), using the previous convergence, (3.14) and the energy relation (3.12), we obtain

$$ \begin{aligned} \int_0^T \int_{\Omega} |\nabla P_m(\varrho)|^2 dx dt &\le \liminf_{\tau \to 0^+} \int_0^T \int_{\Omega} |\nabla P_m(\varrho^\tau)|^2 dx dt \le \limsup_{\tau \to 0^+} \int_0^T \int_{\Omega} |\nabla P_m(\varrho^\tau)|^2 dx dt \\ &\le \limsup_{\tau \to 0^+} \left( \int_{\Omega} (\varrho_0^\tau)^{m+1} dx - \int_{\Omega} (\varrho_T^\tau)^{m+1} dx - \int_0^T \int_{\Omega} \varrho^\tau \nabla P_m(\varrho^\tau) \cdot \nabla \Phi^\tau dx dt \right) \\ &\le \int_{\Omega} (\varrho_0)^{m+1} dx - \int_{\Omega} (\varrho_T)^{m+1} dx - \int_0^T \int_{\Omega} \varrho \nabla P_m(\varrho) \cdot \nabla \Phi[\varrho] dx dt \\ &= \int_0^T \int_{\Omega} |\nabla P_m(\varrho)|^2 dx dt. \end{aligned} $$

From Lemma 3.12 we already had the strong convergence of $P_m(\varrho^\tau)$ to $P_m(\varrho)$ in $L^2([0, T]; L^2(\Omega))$, therefore the previous computation implies convergence of the norms of $P_m(\varrho^\tau)$ in $L^2([0, T]; H^1(\Omega))$, which combined with the weak convergence gives the desired result. $\square$

## 4. QUALITATIVE PROPERTIES

In this section we investigate the qualitative behavior of a simplified PDE obtained as the limit equation of the minimizing movement scheme introduced earlier, but without the boundary effects encoded by the normal cone. All assumptions on the domain and data are those stated in the introduction. The resulting coupled PDE-ODE system reads

$$ (4.1) \qquad \left\{ \begin{array}{l} \partial_t \varrho_t = \Delta P(\varrho_t) + \operatorname{div} \left( \varrho_t \left( \sum_{i=1}^N (x - x_i(t)) \mathbb{1}_{\Omega_i(t)} \right) \right) \\ \dot{x}_i(t) = -a_i(t)x_i(t) + \int_{\Omega_i(t)} x d\varrho_t \\ \dot{a}_i(t) = (-g'(a_i(t)) + \psi_i(t)) \mathbb{1}_{\{a_i > 0\}} \\ \Omega_i(t) = \operatorname{Lag}_i(\psi_t, \mathbf{x}_t), \quad \psi(t) \text{ is a potential for } W_2^2(\varrho_t, \mu_t). \end{array} \right. $$

We first prove that for any solution of (4.1), not only those obtained from the minimizing movement scheme, the atoms remain strictly inside the domain. Moreover, when the boundary $\partial\Omega$ is smooth, any atom initially on the boundary is instantaneously pushed into the interior (see Theorem 4.1). We then focus on the uniform quantization case $a_i \equiv 1/N$, where we show that the atoms stay uniformly separated and that the distance between any point and its Laguerre cell's barycenter converges to zero.
**4.1. Invariant properties of the atomic dynamics.** In this subsection we establish several qualitative properties of the coupled system (1.4). The main result is the following theorem.

**Theorem 4.1.** Let $\Omega$ be an open, bounded, convex set, and let $(\varrho_t, \mathbf{x}_t, \mathbf{a}_t)_{t \in [0,T]}$ be a solution of (1.4). Then the following properties hold:

(1) The family $(\varrho_t)_{t \ge 0}$ is uniformly integrable.

(2) If for some $\bar{t} \ge 0$ we have $a_i(\bar{t}) = 0$, then $a_i(t) = 0$ for all $t \ge \bar{t}$.

(3) If $\partial\Omega$ is of class $\mathcal{C}^1$ and $x_i(0) \in \partial\Omega$, then $x_i(t)$ is instantaneously pushed towards the interior, that is $x_i(t) \in \text{int}\Omega$ for all $t$ in a sufficiently small neighborhood of 0.

(4) Suppose that either $\partial\Omega$ is of class $\mathcal{C}^1$, or $(x_i(0))_{i=1}^N \subset \text{int}\Omega$. Then, for any $t > 0$, if $x_i(t) \in \partial\Omega_i(t)$ it follows that $a_i(t) = 0$.

<figure><img src="image_6.png" /><figcaption>FIGURE 2. Construction of the invariant region $\mathcal{D}$. The uniform integrability of $\varrho_t$ ensures that the barycenter $b_i(t, y)$ remains inside $\Omega_\delta$, allowing the construction of a smooth and convex invariant set $\mathcal{D}$ between the boundaries of $\Omega_\delta$ and $\Omega_{\delta/2}$.</figcaption></figure>

The proof relies on a dynamical systems approach: we construct a positively invariant domain for the dynamics of the atomic positions $(x_i)_{i=1}^N$, while keeping $(\varrho_t, a_i(t))$ fixed. In this context, a positively invariant domain is a region that cannot be left by the trajectories generated by the system.

**Definition 4.2.** Let $V: \mathbb{R} + \mathbb{R}^d \to \mathbb{R}^d$ be a vector field, measurable in time and continuous in space. A set $A \subset \mathbb{R}^d$ is said to be positively invariant for the dynamics induced by $V$ if every solution $(x(t))_{t \ge 0}$ of the ODE $\dot{x}_t = V(t, x_t)$ with $x(0) \in A$ satisfies $x(t) \in A$ for all $t \ge 0$.

The theory of positively invariant sets originates from the work of Nagumo [36] (see also [5, Chap. 4]). Here we exploit the convexity of $\Omega$ together with the uniform integrability of $(\varrho_t)_{t \ge 0}$, inherited from the energy functional $\mathcal{F}$, to apply a simpler classical criterion. In particular, if $A$ has a smooth boundary, a sufficient condition for $A$ to be invariant under the flow induced by $V$ is that $\langle V(t, y), n_A(y) \rangle < 0$ for all $y \in \partial A$, where $n_A(y)$ denotes the outwards unity normal of $\partial A$ evaluated at $y$, see for instance [23, Lemma 1].

**Proof. Item (1):** The uniform integrability of $(\varrho_t)_{t \in [0,T]}$ follows from De la Vallée Poussin's Theorem and the gradient flow structure of the equation on $\varrho_t$. Indeed, using the gradient flow structure of the equation satisfied by $\varrho_t$ we have

$$ \frac{d}{dt} \mathcal{E}(\varrho_t, \mu_t) \le 0, \text{ so that } \mathcal{F}(\varrho_t) \le \mathcal{E}(\varrho_t, \mu_t) \le \mathcal{E}(\varrho_0, \mu_0). $$
**Item (3):** In the case that $\partial\Omega$ is $\mathcal{C}^1$, the normal cone at each point of the boundary is a subspace of dimension 1, that is given by $N_\Omega(x) = n_\Omega(x)\mathbb{R}$, whereas the tangent space to $\partial\Omega$ at $x$, is given by the orthogonal subspace of this vector $T_x\partial\Omega = \{n_\Omega(x)\}^\perp$, but should not be mistaken with Bouligand's tangent cone (3.1) $T_\Omega(x)$. To perform a proof by contradiction, assume that there is an interval $(0, \tau)$ for which $x_i(t) \in \partial\Omega$, then it follows that $\dot{x}_i(0) \in T_{x_i(0)}\partial\Omega$, and we have that $\dot{x}_i(0) \cdot n_\Omega(x_i(0)) = 0$.

Define the family of convex sets

$$ \Omega^\delta \stackrel{\text{def.}}{=} \{x \in \Omega : \text{dist}(x, \partial\Omega) \geq \delta\}, $$

so that $|\Omega \setminus \Omega^\delta| \leq C\delta$. We conclude from the uniform integrability of $(\varrho_t)_{t \in [0,T]}$ that, the baricenter

$$ b_i(t) \stackrel{\text{def.}}{=} \frac{1}{a_i(t)} \int_{\Omega_i(t)} x d\varrho_t \in \text{int} \, \Omega_\delta, $$

for all $\delta$ sufficiently small.

To finish our construction that will lead to a contradiction, notice that $\dot{x}_i(0)$ can be written as

$$ \dot{x}_i(0) = a_i(0)(b_i - x_i(0)) + n_i, \text{ where } n_i \text{ is parallel to } n_\Omega(x_i(0)). $$

Let $x_\delta$ denote the projection of $x_i(0)$ onto the convex set $\Omega^\delta$, so that $x_i(0) - x_\delta = \delta \frac{n_i}{\|n_i\|}$. As a result, the fact that $\dot{x}_i(0) \cdot n_i = 0$ implies, for all $\delta > 0$ small enough so that $b_i \in \Omega^\delta$, that

$$ \begin{aligned} 0 &= \|n_i\|^2 + a_i(0) \langle b_i - x_\delta, n_i \rangle + a_i(0) \langle x_\delta - x_i(0), n_i \rangle \\ &= \|n_i\|^2 + a_i(0) \langle b_i - x_\delta, n_i \rangle - a_i(0) \|n_i\| \delta, \\ &= \|n_i\|^2 + \frac{a_i(0)}{2} \left( \delta^2 + \|b_i - x_\delta\|^2 - \|x_i(0) - b_i\|^2 \right) - a_i(0) \|n_i\| \delta, \end{aligned} $$

where the last equality was obtained via the cosine law for the triangle formed by the points $b_i, x_i(0), x_\delta$. But notice that this equality cannot be true for infinitely many values of $\delta$, which leads to a contradiction with the hypothesis that $x_i(t) \in \partial\Omega$ for all $t \in (0, \tau)$, and we conclude that $x_i$ enters the interior instantaneously.

**Item (4):** Assume by contradiction that $\bar{t} > 0$ is the first time that $x_i$ reaches the boundary $\partial\Omega$ and is such that $a_i(\bar{t}) > 0$. Our approach is to construct an invariant region for the dynamics satisfied by $x_i$. Fixing the trajectories of $(\varrho_t, a_t)_{t \in [0,T]}$, notice that for each $i = 1, \dots, N$ the ODE describing the evolution of $t \mapsto x_i(t)$ is given by

$$ \dot{x}_i(t) = V_i(t, x_i(t)), \text{ where } V_i(t, y) \stackrel{\text{def.}}{=} -\nabla_y W_2^2 \left( \varrho_t, a_i(t)\delta_y + \sum_{j \neq i} a_j(t)\delta_{x_j(t)} \right). $$

Therefore, we have that

$$ V_i(t, y) = a_i(t)(\beta_i(t, y) - y) \text{ where } \beta_i(t, y) \stackrel{\text{def.}}{=} \frac{1}{a_i(t)} \int_{\Gamma_i(t, y)} x d\varrho_t $$

and $\Gamma_i(t, y)$ corresponds to the $i$-th Laguerre cell of the tessellation induced by the semi-discrete transport of $\varrho_t$ to $a_i(t)\delta_y + \sum_{j \neq i} a_j(t)\delta_{x_j(t)}$. As a result, $V_i$ is clearly measurable in time and the $\mathcal{C}^0$ regularity in space comes from the differentiability properties of the semi-discrete transport term [19, Thm 1 and Prop. 2].

To construct an invariant region, set $\delta \stackrel{\text{def.}}{=} \min_{i=1, \dots, N} \text{dist}(x_i(0), \partial\Omega) > 0$ and define

$$ \Omega_s \stackrel{\text{def.}}{=} \{x \in \Omega : \text{dist}(x, \partial\Omega) \geq s\}. $$

Hence, from the uniform integrability and the convexity of $\Omega$, for any $\varepsilon > 0$ there is some $\delta < \bar{\delta}$ and a compact, convex set with smooth boundary $\mathcal{D}_{\bar{\delta}}$ such that

$$ \Omega_{\bar{\delta}/2} \subseteq \mathcal{D}_{\bar{\delta}} \subseteq \Omega, \text{ with } \varrho_t(\Omega \setminus \Omega_{\bar{\delta}}) \leq \varepsilon \text{ for all } t \in [0, T]. $$

Set $\bar{t} \stackrel{\text{def.}}{=} \min\{t \geq 0 : x_i(t) \in \partial\Omega\}$ and suppose by contradiction that $\bar{a}_i \stackrel{\text{def.}}{=} a_i(\bar{t}) > 0$. Performing the previous construction with $\varepsilon \ll \bar{a}_i/4$ it follows that

$$ \beta_i(t, y) \in \Omega_{\bar{\delta}} \text{ for all } t \in [0, T] \text{ and all } y \in \Omega. $$
As a result, since $\mathcal{D}_{\bar{\delta}}$ has smooth boundary and is convex, for any $y \in \partial\mathcal{D}_{\bar{\delta}}$ we get that

$$n_{\mathcal{D}_{\bar{\delta}}}(y) \cdot a_i(t) (\beta_i(t) - y) < 0,$$

meaning that $\mathcal{D}_{\bar{\delta}}$ is invariant for the flow associated with the dynamics of $t \mapsto x_i(t)$. This contradicts the definition of $\bar{t}$, so that $x_i(\bar{t}) \in \partial\Omega$ can only happen if $a_i(\bar{t}) = 0$. $\square$

**4.2. On the dynamics of the optimal quantization equation.** In this paragraph we focus on a simplified form of the coupled system by fixing the weights $a_i = \frac{1}{N}$ for all $i = 1, \dots, N$. In this case, equation (4.1) becomes the gradient flow of the energy

$$\mathcal{E}(\varrho, \mathbf{x}) \stackrel{\text{def.}}{=} \mathcal{H}(\varrho) + F_N(\varrho, \mathbf{x}), \text{ where } F_N(\varrho, \mathbf{x}) \stackrel{\text{def.}}{=} \frac{1}{2} W_2^2 \left( \varrho, \frac{1}{N} \sum_{i=1}^N \delta_{x_i} \right),$$

which corresponds to the uniform quantization problem. This way the equation (4.1) assumes the form

$$(4.2) \qquad \begin{cases} \partial_t \varrho = \Delta \varrho + \operatorname{div}(\varrho \nabla \Phi_t[\varrho_N, \mathbf{x}_t]), & \text{in } (0, \infty) \times \Omega, \\ -\dot{\mathbf{x}}_t = \nabla_{\mathbf{x}} F_N(\varrho_t, \mathbf{x}_t) = \frac{1}{N}(\mathbf{x}_t - \mathbf{b}_t). \end{cases}$$

Existence of solutions is guaranteed from the JKO scheme, provided that the atoms $(x_{i,t})_{i=1}^N$ remain apart and away from the boundary $\partial\Omega$ for all $t > 0$. Using item (4) of Theorem 4.1 with $a_i(t) \equiv 1/N$, this holds if the initial atoms $\mathbf{x}_0$ are distinct and away from the boundary. In addition, since the equation satisfied by $\mathbf{x}_t$ is a gradient flow of a semi-concave functional, atoms cannot collide in finite time. Therefore, we have existence of solutions to (4.2) for all times $t \ge 0$. Similar statements to this were already remarked for instance in [40, Section 6.4.2], [30], but we resume this discussion in further details in Lemma 4.3 below.

A particularly relevant case is when the density is frozen $\varrho_t \equiv \varrho$, so the dynamics of $\mathbf{x}_t$ become a continuous-time analog of Loyd's algorithm for uniform barycenter quantization. For the discrete Loyd algorithm, it is known (see [10, 34, 38]) that the iterates converge toward configurations where each $x_i$ coincides with its barycenter. Our goal here is to establish the analogous continuous result,

$$(4.3) \qquad \|\mathbf{x}_t - \mathbf{b}_t\| \xrightarrow{t \to \infty} 0.$$

and to describe the regularity and stability properties of the flow leading to this convergence.

The evolution of $\varrho_t$ will only affect the analysis through mild regularity assumptions: as the solution of a Fokker-Planck equation with bounded drift, it is Hölder continuous in time and has uniformly bounded densities. Hence several intermediate results, such as Lemma 4.3, remain valid for any continuous curve $\varrho_t$ with these properties, including the stationary case $\varrho_t \equiv \varrho$.

Finally, the gradient flow structure yields finite kinetic energy for $\mathbf{x}_t$

$$(4.4) \qquad \int_0^\infty \|\dot{\mathbf{x}}_t\|^2 dt = \int_0^\infty \|\mathbf{x}_t - \mathbf{b}_t\|^2 dt < +\infty,$$

but this alone does not imply (4.3). Uniform continuity of $t \mapsto \mathbf{x}_t - \mathbf{b}_t$ is also needed. While $\cdot \mathbf{x}_t$ it is clearly Lipschitz with constant diam $\Omega$, the evolution of $\mathbf{b}_t$ can degenerate close to the generalized diagonal

$$(4.5) \qquad \mathbb{D}_N \stackrel{\text{def.}}{=} \{\mathbf{x} \in \Omega^N : x_i = x_j \text{ for some } i \neq j\},$$

see Figure 4.2. Indeed, the dynamics of the atoms become singular when they approach each other. A refined analysis of neighboring cells $i \sim j$, that is cells such that $\mathcal{H}^{d-1}(\Sigma_{ij}) > 0$ where $\Sigma_{ij} \stackrel{\text{def.}}{=} \partial\Omega_i \cap \partial\Omega_j$, leads to quantitative lower bounds on pairwise distances and to global well-posedness, which we establish in Lemma 4.3 below.

**Lemma 4.3.** Let $\Omega$ be a compact and convex subset of $\mathbb{R}^d$, and fix a continuous curve $(\varrho_t)_{t \ge 0}$ in $\mathcal{P}_{ac}(\Omega)$ with uniformly bounded densities. Then the following hold:

(1) The maximal interval of existence for the ODE

$$-\dot{\mathbf{x}}_t = \nabla_{\mathbf{x}} F_N(\varrho_t, \mathbf{x}_t)$$

is given by $\mathbb{R}_+$, provided that the initial atoms $\mathbf{x}_0$ are distinct and either $\partial\Omega$ is $\mathcal{C}^1$ or the initial conditions are away from the boundary $\partial\Omega$.
<figure><img src="image_7.png" /><figcaption>FIGURE 3. A small variation of the atoms does not impact the barycenter too much, unless they are too close. In this situation, where the points are very close, swapping their labels implies a small variation in their positions, but exchanges the barycenters that were very far apart.</figcaption></figure>

(2) For all times $t \ge 0$ the barycenters are away from each other, i.e. there exists $C_\Omega$ such that
$$|b_i(t) - b_j(t)| \ge \frac{C_{\Omega,\varrho}}{N^{d-1}}.$$

(3) If $(\varrho_t, \mathbf{x}_t)_{t \ge 0}$ has finite kinectic energy (4.4), then there exists a constant $c$ for which
$$|x_i(t) - x_j(t)| \ge c,$$
for all $i \sim j$ and for all $t \ge 0$.

*Proof.* To prove the first item, we recall that from [19, Thm. 1], if $\varrho \in \mathcal{P}_{\text{ac}}(\Omega)$ then $\mathbf{x} \mapsto F_N(\varrho, \mathbf{x})$ is has up to 2 derivatives in $\Omega^N \setminus \mathbb{D}_N$, where $\mathbb{D}_N$ is the generalized diagonal defined in (4.5), the set of colliding atoms, and the Hessian is given by

$$\nabla_{x_i, x_j} F_N(\varrho, \mathbf{x}) = \int_{\Sigma_{ij}} (x - x_i) \otimes (x - x_j) \frac{\varrho(x)}{|x_i - x_j|} d\mathcal{H}^{d-1}(x), \text{ for } i \neq j,$$

and

$$\nabla_{x_i, x_i} F_N(\varrho, \mathbf{x}) = \frac{1}{N} I - \sum_{j \neq i} \nabla_{x_i, x_j} F_N(\varrho, \mathbf{x}), \text{ for } i = 1, \dots, N.$$

If, in addition, $\varrho$ is $\mathcal{C}^0(\Omega)$ then $F_N(\varrho, \cdot)$ is $\mathcal{C}^2$ in $\Omega^N \setminus \mathbb{D}_N$. On the other hand, $t \mapsto \nabla F_N(\varrho_t, \mathbf{x})$ is continuous, given $\mathbf{x} \in \Omega^N \setminus \mathbb{D}_N$. Either way, one can check that the Hessian is bounded in compact subsets of $\Omega^N \setminus \mathbb{D}_N$, so that $\mathbf{x} \mapsto \nabla_{x F_N(\varrho, \mathbf{x})}$ is locally Lipschitz, so by the Cauchy-Lipschitz theorem the flow is locally well-posed in $\text{int} \, \Omega^N \setminus \mathbb{D}_N$.

The maximal interval of existence theorem, see for instance [11, Thm. 2.3], says that if $T$ is the supremum of the times for which a solution $\mathbf{x}_t$ to the ODE exists, then either $T = +\infty$ or

$$\lim_{t \to T^{-}} \left( |\mathbf{x}_t| + \frac{1}{\text{dist}(\mathbf{x}_t, \partial(\text{int} \, \Omega^N \setminus \mathbb{D}_N))} \right) = +\infty.$$

For all $t < T$, the atoms $\mathbf{x}_t$ remain from a positive distance to the boundary $\partial\Omega$ due to item (4) of theorem 4.1, provided that $\partial\Omega$ is $\mathcal{C}^1$ or the initial atoms are away from the boundary. In addition, to control the distance of $\mathbf{x}_t$ away from $\mathbb{D}_N$, we proceed as in [30, 34]: since $b_{i,t} \in \Omega_{i,t}$ we have that for $i \neq j$ that

$$\frac{\mathrm{d}}{\mathrm{d}t} \frac{1}{2} |x_{i,t} - x_{j,t}|^2 = - \langle x_{i,t} - x_{j,t}, b_{i,t} - x_{i,t} - (b_{j,t} - x_{j,t}) \rangle \\ = |x_{i,t} - x_{j,t}|^2 + \underbrace{\langle x_{i,t} - x_{j,t}, b_{j,t} - b_{i,t} \rangle}_{\ge 0} \ge |x_{i,t} - x_{j,t}|^2.$$

From the reverse Grönwall inequality it follows that $|x_{i,t} - x_{j,t}|^2 \ge e^{-t}|x_{i,0} - x_{j,0}|^2$ for all $t < T$. We conclude that $\text{dist}(\mathbf{x}_t, \partial(\text{int} \, \Omega^N \setminus \mathbb{D}_N)) \ge c_T > 0$ for all $t < T$ and $T > 0$. Therefore, the only possibility is that $T = +\infty$.

To prove item (2), we use the upper bound on the density $\|\varrho_t\|_{\infty} \le C_{\varrho}$. Given $i \sim j$, assume w.l.o.g. that $\Sigma_{ij} \subset \{x_d = 0\}$. Set $d_{ij}$ such that $\varrho_t(\{x_d > d_{ij}0\}) = 1/2N$. Then it must hold that
$b_i \in \{x_d = d_{ij}\}$ and $\text{dist}(b_i, \partial\Sigma_{ij}) = d_{ij}$. As a result, we have that

$$ \frac{1}{2N} = \varrho_t(\{x_d > d_{ij}0\}) \leq C_{\varrho}|\{x_d > d_{ij}\} \cap \Omega| \leq C_{\varrho} \text{diam} \Omega d_{ij}^{d-1}, $$

and the result follows by rearranging the above inequality and bounding $|b_i - b_j| \geq \text{dist}(b_i, \Sigma_{ij}) + \text{dist}(b_j, \Sigma_{ij})$.

Moving on to item (3), suppose by contradiction that for every $\varepsilon > 0$ there exists a sequence $t_k \xrightarrow{k \to \infty} +\infty$ for which there is a pair $i_k \sim j_k$ such that

$$ |x_{i_k}(t_k) - x_{j_k}(t_k)| \leq \frac{\varepsilon}{3}. $$

Since $t \mapsto \mathbf{x}_t$ is Lipschitz continuous with constant $\text{diam} \Omega$, we see that there is an interval $I_k$ centered at $t_k$ of size $|I_k| = 2\varepsilon/3 \text{diam} \Omega$ such that

$$ |x_{i_k}(t) - x_{j_k}(t)| \leq |x_{i_k}(t_k) - x_{i_k}(t)| + |x_{i_k}(t_k) - x_{j_k}(t_k)| + |x_{j_k}(t_k) - x_{j_k}(t)| \leq \varepsilon $$

for all $t \in I_k$. In particular, taking $\varepsilon$ much smaller than $\frac{C_{\Omega,\varrho}}{N^{d-1}} \leq |b_i - b_j|$, so that $\frac{C_{\Omega,\varrho}}{N^{d-1}} \leq \|\mathbf{x}_t - \mathbf{b}_t\|$ inside $I_k$ for all $k \in \mathbb{N}$. In this case we would have

$$ \int_{0}^{\infty} \|\mathbf{x}_t - \mathbf{b}_t\|^2 dt \geq \sum_{k \in \mathbb{N}} \int_{I_k} \|\mathbf{x}_t - \mathbf{b}_t\|^2 dt = +\infty. $$

This gives a contraction and hence $|x_i - x_j|$ must be uniformly bounded from below. $\square$

The above Lemma gives a uniform lower bound on the distance of $\mathbf{x}_t$ to the generalized diagonal, meaning that the dynamics does not degenerate asymptotically. However, the evolution of the $\mathbf{b}_t$ depends on the smoothness of the evolution of the free boundaries $\Omega_i(t) = \text{Lag}_i(\mathbf{x}_t, \psi_t)$, that is on the regularity of $\psi_t$. Since $t \mapsto \mathbf{x}_t$ is uniformly Lipschitz, we shall see below that the smoothness of $t \mapsto \psi_t$ depends essentially on the evolution of $\varrho_t$.

**Lemma 4.4.** Let $\Omega$ be a compact and convex subset of $\mathbb{R}^d$, and fix a continuous curve $(\varrho_t)_{t \geq 0}$ in $\mathcal{P}_{\text{ac}}(\Omega)$ with uniformly bounded densities, and such that $t \mapsto \varrho_t$ is uniformly continuous in time, with modulus of continuity $\omega_{\varrho}$. Then $t \mapsto \psi_t$, and consequently $t \mapsto \mathbf{b}_t$, are uniformly continuous with modulus of continuity $\omega_{\varrho}$.

*Proof.* The regularity of evolution of barycenters is dictated by the regularity of the evolution of the curve $\varrho_t$ and of the optimal Kantorovich potentials $\psi_t$. To see this, define the function $G_t: \Omega^N \times \mathbb{R}^N \to \mathbb{R}^{N+1}$ as

$$ G_t: \Omega^N \times \mathbb{R}^N \to E \stackrel{\text{def.}}{=} \{\sum \psi_i = 0\} \quad \text{where } G_{i,t}(\mathbf{x}, \psi) \stackrel{\text{def.}}{=} \varrho_t(\text{Lag}_i(\mathbf{x}, \psi)) - \frac{1}{N}. \quad (\mathbf{x}, \psi) \mapsto (G_{i,t}(\mathbf{x}, \psi))_{i=1}^N, $$

It then follows that for all $t \geq 0$, at the evolution of the gradient flow we have $G_t(\mathbf{x}_t, \psi_t) \equiv 0$. In addition, the Jacobian of $G_t$ can be computed explicitly as

$$ \frac{\partial G_{t,i}}{\partial \psi_j} = - \int_{\Sigma_{ij}} \frac{\varrho_t(x)}{|x_j - x_i|} d\mathcal{H}^{d-1}(x) \text{ for } j \neq i, \quad \frac{\partial G_{t,i}}{\partial \psi_i} = - \sum_{j \neq i} \frac{\partial G_{t,i}}{\partial \psi_j}, $$

$$ \frac{\partial G_{t,i}}{\partial x_j} = - \int_{\Sigma_{ij}} \frac{x - x_i}{|x_j - x_i|} \varrho_t(x) d\mathcal{H}^{d-1}(x) \text{ for } j \neq i, \quad \frac{\partial G_{t,i}}{\partial x_i} = - \sum_{j \neq i} \frac{\partial G_{t,i}}{\partial x_j}. $$

Therefore, for each $t \geq 0$ the Jacobian $D_{\psi}G_t$ has rank $N-1$ and we can apply the implicit function theorem, because $E$ is a $N-1$-dimensional hyperplane, to write $\psi_t$ as a $\mathcal{C}^1$ function of $\mathbf{x}_t$, $\psi_t = g_t(\mathbf{x}_t)$. In addition, from Lemma 4.3, since the points $x_i$ are uniformly separated, the Jacobian $D_{\psi}G_t$ has a bounded inverse in $E$. We can therefore estimate the sensitivity of $\psi_t$ on time by noticing that

$$ 0 = \underbrace{G_t(\mathbf{x}_t, \psi_t) - G_t(\mathbf{x}_s, \psi_s)}_{\stackrel{\text{def.}}{=} \Delta_{x,\psi} G} + \underbrace{G_t(\mathbf{x}_s, \psi_s) - G_s(\mathbf{x}_s, \psi_s)}_{\stackrel{\text{def.}}{=} \Delta_{t,s} G} $$

and estimating each term.
For the first one, from the $\mathcal{C}^1$ dependence of $\psi_t$ on $\mathbf{x}_t$ for $t$ fixed, we can apply the mean value theorem to find $(\bar{\mathbf{x}}, \bar{\psi})$ such that

$$ \Delta_{x,\psi} G = D_{\mathbf{x}} G_t(\bar{\mathbf{x}}, \bar{\psi})(\mathbf{x}_t - \mathbf{x}_s) + D_{\psi} G_t(\bar{\mathbf{x}}, \bar{\psi})(\psi_t - \psi_s). $$

For the second term, each component is controlled by the modulus of continuity in time of $\varrho_t$:

$$ \Delta_{t,s} G_i = \varrho_t(\mathrm{Lag}_i(\mathbf{x}_s, \psi_s)) - \varrho_s(\mathrm{Lag}_i(\mathbf{x}_s, \psi_s)). $$

As a result, since the Jacobian w.r.t. $\psi$ has bounded inverse and $t \mapsto \mathbf{x}_t$ is Lipschitz, we get that

$$ |\psi_t - \psi_s| \le \|D_{\psi} G_t(\bar{\mathbf{x}}, \bar{\psi})^{-1}\|_{\infty} (\mathrm{diam} \, \Omega |t-s| + \omega_{\varrho}(|t-s|)). $$

This implies that $t \mapsto \psi_t$ has the same modulus of continuity of $t \mapsto \varrho_t$ in $L^1$, and from the uniform continuity of $(\mathbf{x}, \psi) \mapsto \mathrm{Lag}_i(\mathbf{x}_t, \psi_t)$ in the strong topology of $L^1(\Omega)$, we conclude that the curve of barycenters $t \mapsto \mathbf{b}_t$ also inherits the same modulus of continuity of $\varrho_t$. $\square$

We can now synthesize these arguments into the desired convergence of the distance to barycenters (4.3).

**Theorem 4.5.** If either $\varrho_t \equiv \varrho$ is fixed, or $(\varrho_t, \mathbf{x}_t)_{t \ge 0}$ solve (4.2), then

$$ \|\mathbf{x}_t - \mathbf{b}_t\| \xrightarrow{t \to \infty} 0. $$

*Proof.* If either $\varrho_t \equiv \varrho$ is fixed, or if $(\varrho_t, \mathbf{x}_t)_{t \ge 0}$ solves (4.2), then a simple gradient flow argument gives that $\mathbf{x}_t$ has finite kinetic energy (4.4). For instance, in the latter case notice that the time derivative of the energy is given by

$$ \frac{d}{dt} \mathcal{E}(\varrho_t, \mathbf{x}_t) = -I(t) - \frac{1}{N^2} \|\mathbf{x}_t - \mathbf{b}_t\|^2, $$

where $I(t)$ corresponds to the Fisher information of $\varrho_t$ w.r.t. the Gibbs measure associated with the Kantorovich potential $\Phi_t$:

$$ I(t) = \int_{\Omega} \left| \frac{\nabla \varrho_t}{\varrho_t} + \nabla \Phi_t(x) \right|^2 d\varrho_t \ge 0. $$

Therefore integrating over $\mathbb{R}_+$ we get that

$$ \frac{1}{N^2} \int_0^\infty \|\mathbf{x}_t - \mathbf{b}_t\|^2 dt \le \mathcal{E}(\varrho_0, \mathbf{x}_0) < +\infty. $$

In addition, it also holds that $t \mapsto \varrho_t$ is uniformly continuous in $L^1(\Omega)$. In fixed density case this holds trivially and in the latter case it is a consequence of being the solution of a Fokker-Planck equation with bounded drift having therefore a bounded and Hölder continuous density, see for instance [18] and [6].

Using Lemmas 4.3 and 4.4 above, we then have from the previous discussion that

$$ \int_0^\infty \|\mathbf{x}_t - \mathbf{b}_t\|^2 dt < +\infty \quad \text{and} \quad t \mapsto \|\mathbf{x}_t - \mathbf{b}_t\| \text{ is uniformly continuous,} $$

as a consequence $\|\mathbf{x}_t - \mathbf{b}_t\| \xrightarrow{t \to \infty} 0$, and the result follows. $\square$

## 5. NUMERICAL SIMULATIONS AND CONJECTURES

In this section we describe a splitting scheme for the numerical simulation of the dynamic quantization equation (1.5). Our goal with the numerical experiments is to help elaborate conjectures on the long-time behavior of our system. The scheme consists on discretizing the time interval where one wishes to compute solutions and update alternatively the approximations for the continuous density $\varrho$ and the atoms' positions $\mathbf{x}$.

The PDE that gives the evolution of $\varrho$ is discretized on a mesh $\mathcal{T}$ that is independent of the evolution of the points $\mathbf{x}$ or the Laguerre tessellation induced by the semi-discrete optimal transportation problem. This allows for a splitted scheme; at each time step we first evolve the positions $\mathbf{x}$ with one step of an Euler scheme for their underlying ODE with the density fixed, which can be easily computed with the `pysdot` package [29]. This evolution has the advantage of being performed off the grid $\mathcal{T}$. Afterwards the evolution of $\varrho$ can be computed efficiently with general purpose solvers of Fokker-Planck type equations such as `fipy` [25].
From one hand, better integrating the evolution of Laguerre cells and the underlying mesh $\mathcal{T}$ which determines the resolution of the continuous density could lead to numerical methods that take into account the flow of mass entering and exiting each Laguerre cell, for instance with a finite volume approach [15]. On the other hand, this approach with an independent mesh and Laguerre tesselation allows for an off the grid approach for the atomic measure's evolution, since then the Laguerre tessellation is uniquely determined by the scalar Kantorovitch potentials $(\psi_i)_{i=1}^N$. This description of our method is summarized in Algorithm 1.

<caption>**Algorithm 1:** Splitting scheme for (1.5)</caption>

**Input:** initial density $\rho_0$ on a mesh $\mathcal{T}$, initial points $\mathbf{x}_0$, time step $\tau$
**Output:** sequence of approximations $(\rho_n, \mathbf{x}_n)$
**for** $n = 0, 1, 2, \dots$ **do**
    // Compute OT (semi-discrete) with pysdot
    With pysdot to compute $\psi_n = (\psi_{n,i})_{i=1}^N$ and $\mathbf{b}_n = (b_{n,i})_{i=1}^N$ corresponding to
    $W_2^2 (\rho_n, \frac{1}{N} \sum_{i=1}^N \delta_{x_{n,i}})$;
    // Compute OT map via optimality conditions
    $T_n(x) = x_{n,i^*}$ where $i^* = \underset{i=1, \dots, N}{\operatorname{argmin}} \frac{1}{2} \|x - x_{n,i}\|^2 - \psi_{n,i}$
    from optimality conditions for semi-discrete OT;
    // Point ODE update with explicit Euler
    Update atoms $\mathbf{x}_{n+1} = \mathbf{x}_n - \alpha(\mathbf{x}_n - \mathbf{b}_n)$ (e.g. $\alpha = \tau N$);
    // Fix advection over interval $[\tau n, \tau(n+1)]$
    Set advection $\nabla\Phi_n(x) = x - T_n(x)$;
    // Solve PDE with fipy over interval $[\tau n, \tau(n+1)]$
    Solve $\partial_t \rho^n = \Delta\rho^n + \operatorname{div}(\rho^n \nabla\Phi_n)$ with $\rho^n(0) = \rho_n$ and set $\rho_{n+1} = \rho^n(\tau)$;
**end**

On all experiments, the dynamics of $\mathbf{x}$ are slightly modified in order to observe the long-time behavior more quickly; more specifically we multiply their dynamics by a multiplicative factor

$$ \dot{\mathbf{x}}_t = \alpha (\mathbf{b}_t - \mathbf{x}_t). $$

Since we know that the distance to the barycenters will go to zero from Theorem 4.5, it makes sense to add the multiplicative constant in the dynamics since for sufficiently large times $|x_i(t) - b_i(t)|$ is at most of the other of the diam $\Omega_i(t)$ which, heuristically is of the order $N^{-1/d}$, the typical distance between minimizers of the uniform quantization problem. This speed parameter serves therefore to accelerate the convergence of the distance to the barycenters and observe the asymptotic behavior more quickly.

Before discussing the experiments, we briefly emphasize that the numerical results below are exploratory and primarily intended to illustrate qualitative features of the long-time behavior of (4.2).

**A crystallization phenomenon.** On Figures 4–5, we display the approximate steady states obtained for increasing numbers of points $N = 50, 100, 200, 300, 400, 500$ and a velocity scaling $\alpha_N = \sqrt{N}$. In the case of linear diffusion, shown in Figure 4, the steady configurations of the atoms tend toward a uniform triangular lattice as $N$ increases, revealing a clear dynamic crystallization effect. The color map in this figure is not uniform across subplots, in order to highlight the tendency of the diffuse density $\rho$ to concentrate according to a Gibbs-type weight with respect to the optimal Kantorovich potential inside each Laguerre cell. Formally, for large $N$ and $t$, the stationary configuration is well approximated by

$$ (5.1) \qquad \rho_N(t) \propto \sum_{i=1}^{N} e^{-(\frac{1}{2}|x-x_i(t)|^2 - \psi_i(t))} 1_{\Omega_i(t)}. $$
As $N$ grows, the Laguerre cells $\Omega_i$ shrink and $\rho_N$ becomes nearly constant within each cell, which makes (5.1) visually appear as a uniform density. This behavior is consistent with the crystallization phenomena established in [8,9] for models where the diffuse component is homogeneous. The trend is further illustrated in Figures 6 and 7, which use a common color map across time to show that the densities flatten as equilibrium is approached.

<figure><img src="image_8.png" /><figcaption>FIGURE 4. Crystallization phenomenon for the model with linear diffusion and speed $\alpha = \sqrt{N}$. Long-time configurations for $N = 50, 100, 200, 300, 400, 500$ after 200 time steps with $\tau = 0.01$. The density $\rho$ stabilizes to an almost uniform profile, leading to a static triangular lattice of Laguerre cells.</figcaption></figure>

A similar but less pronounced crystallization effect is observed for the evolution with a Porous Medium-type diffusion term, shown in Figure 5. Here we set $P(\rho) = \rho^m$ with $m = 10$, again with $\alpha = \sqrt{N}$. The nonlinear diffusion favors locally concentrated regions, producing a sharper crystalline pattern near zones of higher density.

**Discussion on the long-time behavior of the densities.** We now turn to a closer inspection of the long-time behavior of the diffuse component $\rho_t$. As suggested by the previous section, it is natural to conjecture that stationary configurations of the coupled system should approach the Gibbs-type profile (5.1). However, this correspondence does not always manifest clearly for small numbers of atoms.

Figure 6 shows the evolution of (4.2) with only three atoms under linear diffusion. While the atoms converge to the barycenters of their respective Laguerre cells, the resulting tessellation does not align perfectly with the visible structure of the density. This suggests that for small $N$, the coupling between the discrete and continuous components is too coarse for the equilibrium state to fully reflect the formal stationary form (5.1).

The situation changes markedly for larger $N$. In Figure 7, we repeat the same experiment with $N = 50$ atoms. Although the color map is kept fixed across all times (which visually exaggerates the flattening), the final configuration shows an almost constant density, consistent with the stationary structure (5.1) and with the crystalline organization observed in Figure 4. This supports the idea that as $N \to \infty$, the discrete measure of atoms becomes dense enough to recover the expected macroscopic equilibrium.
<figure><img src="image_9.png" /><figcaption>FIGURE 5. Crystallization phenomenon for the model with PME-type diffusion, i.e. $\Delta\varrho^m$ with $m = 10$ and speed $\alpha = \sqrt{N}$. Shown are steady configurations for $N = 50, 100, 200, 500, 800, 1000$ after 200 time steps with $\tau = 0.01$. The nonlinear diffusion enhances localization and produces sharper crystalline regions in high-density zones, though small fluctuations remain visible even at long times.</figcaption></figure>

Figures 8 and 9 illustrate the corresponding dynamics for the porous-medium variant of the model, in which the diffusion term is replaced by $\Delta P(\varrho)$ with $P(\varrho) = \varrho^m$ and $m = 10$. The first feature we observe in this evolution is that, even for a well-spread initial condition (here a truncated Gaussian), the early dynamics are dominated by the attraction toward the atomic measure. This leads to a rapid concentration of mass around the points $x_t$, followed by a slower diffusion-driven relaxation toward profiles typically favored by the porous-medium equation. This two-stage behavior highlights the subtle balance between aggregation and diffusion inherent to the coupled system.

Although the nonlinear diffusion case is not directly covered by the analysis of Theorem 4.5, the numerical results consistently show convergence of each atom $x_i$ toward the barycenter of its Laguerre cell. This behavior is in line with the known regularizing effects of porous-medium equations, which yield Hölder-continuous densities [13, 14, 42]. Even though rigorous results for the coupled setting are lacking, the persistence of convergence strongly suggests that this regularity continues to stabilize the underlying ODE dynamics.

A different picture emerges when the velocity scaling is reduced to $\alpha = 1$, as opposed to the accelerated $\alpha_N = \sqrt{N}$ used in the crystallization experiments. In this slower regime, especially for large $N$ (e.g. $N = 500$), the density no longer spreads uniformly but instead tends to form a single large connected component of higher concentration, as shown in Figure 9. This indicates a competing effect between nonlinear diffusion, which promotes compactly supported profiles, and the attraction to barycenters, which drives spatial homogenization.

## REFERENCES

[1] L. Ambrosio, N. Fusco, and D. Pallara. *Functions of Bounded Variation and Free Discontinuity Problems*. Oxford Mathematical Monographs. Oxford University Press, Oxford, New York, Mar. 2000.

[2] L. Ambrosio, N. Gigli, and G. Savaré. *Gradient flows in metric spaces and in the space of probability measures*. Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, 2005.
<figure><img src="image_10.png" /><figcaption>FIGURE 6. Long-time behavior for the model with linear diffusion and $N = 3$ atoms. The atoms converge to their barycenters, but the diffuse density $\rho_t$ does not coincide with the Gibbs-type stationary profile associated with $\Phi_t[\mathbf{x}, \psi]$.</figcaption></figure>

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

[4] M. Bertsch and D. Hilhorst. A density dependent diffusion equation in population dynamics: stabilization to equilibrium. *SIAM Journal on Mathematical Analysis*, 17(4):863–883, 1986.

[5] F. Blanchini and S. Miani. *Set-theoretic methods in control*, volume 78 of SCFA. Birkhäuser Cham, 2 edition, 2015.

[6] V. I. Bogachev, N. V. Krylov, M. Röckner, and S. V. Shaposhnikov. *Fokker-Planck-Kolmogorov Equations*, volume 207 of *Mathematical Surveys and Monographs*. American Mathematical Society, 2015.

[7] G. Bouchitté, C. Jimenez, and R. Mahadevan. Asymptotic analysis of a class of optimal location problems. *Journal de mathématiques pures et appliquées*, 95(4):382–419, 2011.

[8] D. P. Bourne and R. Cristoferi. Asymptotic optimality of the triangular lattice for a class of optimal location problems. *Communications in Mathematical Physics*, 387(3):1549–1602, 2021.

[9] D. P. Bourne, M. A. Peletier, and F. Theil. Optimality of the triangular lattice for a particle system with wasserstein interaction. *Communications in Mathematical Physics*, 329(1):117–140, 2014.

[10] D. P. Bourne and S. M. Roper. Centroidal power diagrams, Lloyd's algorithm, and applications to optimal location problems. *SIAM Journal on Numerical Analysis*, 53(6):2545–2569, 2015.

[11] A. Bressan and B. Piccoli. *Introduction to the mathematical theory of control*, volume 1. American institute of mathematical sciences, 2007.

[12] G. Buttazzo and F. Santambrogio. A model for the optimal planning of an urban area. *SIAM journal on mathematical analysis*, 37(2):514–530, 2005.

[13] L. A. Caffarelli and A. Friedman. Continuity of the density of a gas flow in a porous medium. *Transactions of the American Mathematical Society*, 252:99–113, 1979.

[14] L. A. Caffarelli and A. Friedman. Regularity of the free boundary of a gas flow in an n-dimensional porous medium. *Indiana University Mathematics Journal*, 29(3):361–391, 1980.

[15] C. Cances, T. O. Gallouët, and G. Todeschi. A variational finite volume scheme for wasserstein gradient flows. *Numerische Mathematik*, 146(3):437–480, 2020.

[16] G. Carlier and M. Laborde. A splitting method for nonlinear diffusions with nonlocal, nonpotential drifts. *Nonlinear Analysis: Theory, Methods & Applications*, 150:1–18, 2017.

[17] G. Carlier and C. Poon. On the total variation wasserstein gradient flow and the tv-jko scheme. *ESAIM: Control, Optimisation and Calculus of Variations*, 25:42, 2019.

[18] L. Chizat, M. Colombo, and X. Fernández-Real. Convergence of drift-diffusion pdes arising as wasserstein gradient flows of convex functions. *arXiv preprint arXiv:2507.12385*, 2025.
<figure><img src="image_11.png" /><figcaption>FIGURE 7. Long-time behavior for the model with linear diffusion and $N = 50$ atoms. The diffuse component $q_t$ becomes nearly uniform, consistent with the asymptotic configuration described in (5.1). The apparent constancy of the density is accentuated by using the same color scale for all times.</figcaption></figure>

[19] F. De Gournay, J. Kahn, and L. Lebrat. Differentiation and regularity of semi-discrete optimal transport with respect to the parameters of the discrete measure. *Numerische Mathematik*, 141:429–453, 2019.

[20] G. De Philippis, A. R. Mészáros, F. Santambrogio, and B. Velichkov. Bv estimates in optimal transportation and applications. *Archive for Rational Mechanics and Analysis*, 219:829–860, 2016.

[21] M. Di Francesco and D. Matthes. Curves of steepest descent are entropy solutions for a class of degenerate convection–diffusion equations. *Calculus of Variations and Partial Differential Equations*, 50:199–230, 2014.

[22] S. Di Marino and F. Santambrogio. JKO estimates in linear and non-linear Fokker–Planck equations, and Keller–Segel: $L^p$ and Sobolev bounds. *Ann. Inst. H. Poincaré Anal. Non Linéaire*, 39(6):1485–1517, 2022.

[23] M. Fernandes and F. Zanolin. Remarks on strongly flow-invariant sets. *Journal of Mathematical Analysis and Applications*, 128:176–188, 1987.

[24] I. Fonseca and G. Leoni. *Modern methods in the Calculus of Variations: $L^p$ spaces*. Springer Science & Business Media, 2007.

[25] J. E. Guyer, D. Wheeler, and J. A. Warren. Fipy: Partial differential equations with python. *Computing in Science & Engineering*, 11(3):6–15, 2009.

[26] J.-B. Hiriart-Urruty and C. Lemaréchal. *Convex analysis and minimization algorithms I: Fundamentals*, volume 305. Springer science & business media, 1996.

[27] 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.

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

[29] H. Leclerc. pysdot: Semi-discrete optimal transport library (python), 2025.

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

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

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

[33] Q. Mérigot and J.-M. Mirebeau. Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport. *SIAM Journal on Numerical Analysis*, 54(6):3465–3492, 2016.

[34] Q. Mérigot, F. Santambrogio, and C. Sarrazin. Non-asymptotic convergence bounds for wasserstein approximation using point clouds. *Advances in Neural Information Processing Systems*, 34:12810–12821, 2021.

[35] Q. Merigot and B. Thibert. Optimal transport: discretization and algorithms. In *Handbook of numerical analysis*, volume 22, pages 133–212. Elsevier, 2021.
<figure><img src="image_12.png" /><figcaption>FIGURE 8. Long-time behavior for the nonlinear diffusion model with $P(\rho) = \rho^{10}$, $N = 5$ atoms, and $\alpha = 1$. The convergence of atoms to barycenters persists, although the equilibrium density remains spatially nonuniform.</figcaption></figure>

[36] M. Nagumo. Über die lage der integralkurven gewöhnlicher differentialgleichungen. *Proceedings of the Physico- Mathematical Society of Japan.*, 24:551–559, 1942.

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

[38] L. Portales, E. Cazelles, and E. Pauwels. On the sequential convergence of Lloyd's algorithms. *Mathematics of Operations Research*, 2025.

[39] 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*, 2(2):395–431, 2003.

[40] F. Santambrogio. Optimal transport for applied mathematicians. *Birkhäuser, NY*, 55(58-63):94, 2015.

[41] F. Santambrogio and G. Toshpulatov. Strong $L^2H^2$ convergence of the JKO scheme for the Fokker-Planck equation. *Archive for Rational Mechanics and Analysis*, 248(6):99, 2024.

[42] J. L. Vázquez. *The porous medium equation: mathematical theory*. Oxford university press, 2007.

[43] C. Villani et al. *Optimal transport: old and new*, volume 338. Springer, 2009.
<figure><img src="image_13.png" /><figcaption>FIGURE 9. Long-time behavior for the nonlinear diffusion model with $P(\varrho) = \varrho^{10}$, $N = 500$ atoms, and $\alpha = 1$. The slower dynamics lead to the formation of a large, connected concentration region, contrasting with the uniform crystallization pattern observed for faster scaling.</figcaption></figure>