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

HAL
open science

# A variational formulation of a Multi-Population Mean Field Games with non-local interactions

Luigi de Pascale, Luca Nenna

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

## To cite this version:

Luigi de Pascale, Luca Nenna. A variational formulation of a Multi-Population Mean Field Games with non-local interactions. 2025. hal-04668323v3

HAL Id: hal-04668323

https://inria.hal.science/hal-04668323v3

Preprint submitted on 31 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A variational formulation of a Multi-Population Mean Field Games with non-local interactions

Luigi De Pascale,*, Luca Nenna†

October 31, 2025

We propose a MFG model with quadratic Hamiltonian involving $N$ populations. This results in a system of $N$ Hamilton-Jacobi-Bellman and $N$ Fokker-Planck equations with non-local interactions. As in the classical case we introduce an Eulerian variational formulation which, despite the non convexity of the interaction, still gives a weak solution to the MFG model. The problem can be reformulated in Lagrangian terms and solved numerically by a Sinkhorn-like scheme. We present numerical results based on this approach, these simulations exhibit different behaviors depending on the nature (repulsive or attractive) of the non-local interaction.

*Keywords*. Mean field games, optimal transport, multi-marginal optimal transport, entropic regularization, convex analysis.

2020 Mathematics Subject Classification. Primary: 49L25; Secondary: 49K40, 49N15, 65K10, 91A14.

## 1. Introduction

The aim of this paper is to establish a variational formulation (both Lagrangian and Eulerian), as well as a suitable numerical methods, for quadratic second-order Mean Field Games which involve $N$ different populations interacting through a given non-local functional. Let $d$ be the dimension of the space; then for $i = 1, \dots, N$ we consider

<footnote>*Dipartimento di Matematica ed Informatica, Università di Firenze, Viale Morgagni, 67/a - 50134 Firenze, Italy.</footnote>

<footnote>†Université Paris-Saclay, CNRS, Laboratoire de mathématiques d'Orsay, ParMA, Inria Saclay, 91405, Orsay, France. email: luca.nenna@universite-paris-saclay.fr</footnote>
the following coupled PDE system.

$$
\left\{
\begin{array}{l}
-\partial_t u^i - \frac{1}{2} \Delta u^i + \frac{1}{2} |\nabla u^i|^2 = \sum_{j \neq i} \int_{\mathbb{R}^d} V^{i,j}(x-y) \rho^j dy, \\
\\
\partial_t \rho^i - \frac{1}{2} \Delta \rho^i + \text{div}(\nabla u^i \rho^i) = 0, \\
\\
\rho^i(0,x) = \rho_0^i(x), \quad u^i(T,x) = g^i(x);
\end{array}
\right.
\qquad (1.1)
$$

where $\rho^i(0, x) \in \mathcal{P}_{ac}(\mathbb{R}^d)$ and $g^i \in \mathbb{C}^0(\mathbb{R}^d)$ and $V^{i,j}$ is a given positive and lower semi-
continuous potential, the full set of assumptions is given at the end of this introduction.
In the paper, we will mostly use the more compact notation $\rho_t$ for $\rho(t, x)$.

Mean-field games involving several populations have been studied recently by [19, 20, 15, 3, 9, 22, 31, 26]; for some more details about the theory of MFG we refer the reader to the seminal work by Lasry and Lions [24], the book [17] or the lecture notes by Cardaliaguet [13]. Notice that the interaction of the populations can be expressed via a more general functional, but we will not discuss here the extra difficulties, or via some Optimal Transport coupling as done in [3]. We will show that the above system can be seen as optimality conditions for the following variational problem *Eulerian*.

$$
\inf \left\{ \mathcal{J}(\rho^1, v^1; \dots; \rho^N, v^N) \mid \partial_t \rho_t^i - \frac{1}{2} \Delta \rho_t^i + \text{div}(\rho_t^i v_t^i) = 0, \rho^i(0, x) = \rho_0^i, \forall i \right\}
$$

where

$$
\begin{align}
\mathcal{J}(\rho^1, v^1; \dots; \rho^N, v^N) &:= \sum_i \int_{[0,T] \times \mathbb{R}^d} \frac{|v^i|^2}{2} d\rho_t^i dt + \nonumber \\
&\quad \sum_{1 \le j \le N \atop j \ne i} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V^{i,j}(x-y) d\rho_t^i \otimes d\rho_t^j dt + \int_{\mathbb{R}^d} g^i(x) d\rho^i(T,x). \tag{1.2}
\end{align}
$$

Moreover, we will also relate the above minimization problem with a *Lagrangian* relative entropy minimization problem, that is

$$
\min\{J(Q^1, \dots, Q^N) : e_{\#}^0 Q^i = \rho_0^i\} \qquad (1.3)
$$

where

$$
\begin{align*}
J(Q^1, \dots, Q^N) &:= \sum_i \mathcal{H}(Q^i | R^i) + \sum_{\substack{1 \le i \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) de_{\#}^t Q^i \otimes e_{\#}^t Q^j dt \\
&\quad + \sum_i \int_{\Omega} g^i(\omega(T)) dQ^i,
\end{align*}
$$
where $Q^i \in \mathcal{P}(\Omega)$ with $\Omega = \mathcal{C}([0, T]; \mathbb{R}^d)$ and and $\mathcal{H}(P|R)$ is the Boltzmann-Shannon entropy, that is

$$ \mathcal{H}(\gamma | \pi) = \int \rho \log \rho d\pi, \text{ if } \gamma = \rho\pi. $$

When the reference measure $\pi$ in the entropy is not indicated, is intended the Lebesgue measure.

We denote the $N$-uple by $\underline{\rho} := (\rho^1, \dots, \rho^N)$.

The existence and uniqueness of weak solutions for this system will be discussed mainly in the full space $\mathbb{R}^d$, but all of our results can be easily adapted to the periodic setting of $\mathbb{T}^d$. For $\mathbb{R}^d$ we will need some condition at $\infty$ then we will use the spaces:

$$ (i) \quad \mathcal{C}_0^2([0, T] \times \mathbb{R}^d) = \{\phi \in \mathcal{C}^2([0, T] \times \mathbb{R}^d) : \lim_{|x| \to +\infty} \phi(t, x) = 0, \forall t\}; $$

$$ (ii) \quad \mathrm{Lip}_0(\mathbb{R}^d) = \{u \in \mathrm{Lip}(\mathbb{R}^d) : \lim_{|x| \to +\infty} u(x) = 0\}. $$

We will make the following assumptions.

$$ (A1) \quad \rho_0^i \in \mathrm{Lip}_0(\mathbb{R}^d) \text{ (or } \rho_0^i \in \mathrm{Lip}(\mathbb{T}^d)); $$

$$ (A2) \quad \mathcal{H}(\rho_0^i) < +\infty; $$

$$ (A3) \quad g^i \in \mathcal{C}_0^2(\mathbb{R}^d), $$

$$ (A4) \quad V^{i,j} * \rho \in \mathcal{C}_b(\mathbb{R}^d), \text{ for all } \rho \text{ such that } \sqrt{\rho} \in H^1(\mathbb{R}^d). $$

**Remark 1.1.** Assumption (A4) above is always satisfied if $V^{i,j} \in \mathcal{C}_b(\mathbb{R}^d)$ but also, for example, if $V^{i,j}(x) = \frac{1}{|x|^\alpha}$ with $\alpha$ such that $\frac{1}{|x|^\alpha} \in L_{\mathrm{loc}}^q$ for a $\frac{d}{2d-2} \le q \le \infty$. So, if $d=3$, the Coulomb cost is allowed. The condition $\sqrt{\rho} \in H^1(\mathbb{R}^d)$ may look unnatural, however, since $\rho$ is a probability measure, it corresponds to the fact that the Fisher information $I(\rho)$ is finite. Lemma 4.1 of [7] shows that this is always the case in this setting. The relevance of the condition $\sqrt{\rho} \in H^1(\mathbb{R}^d)$ is also given by the fact that such $\rho$'s are the electron densities associated to wave functions in quantum mechanics [25] and, in this respect, the Coulomb type cost discussed above is also relevant (notice that in [24] the authors linked the MFG system with Hartree type equation in which case the Coulomb interaction arises quite naturally).

Also, assumption (A4) is related to Th. 2.6. Choosing to use duality with different functions spaces would, of course, requires different assumptions.

The aim of this work is twofold. On the one hand, we introduce a non-local interaction term, which is not only interesting for the applications (see for instance [19, 20, 15, 3, 9, 22, 31] and the references therein) but also introduces a slight non-convexity. On the other hand, we try to strip naked the structure of the problem so to be as accessible as possible to non-specialists and students.

Compared to [5] we allow a different kind of non-locality and non-convexity in the interaction term and, also, aside with the periodic setting, we consider problems on the
entire space $\mathbb{R}^d$. Working on the entire space requires some care in the use of Fenchel-Rockafellar duality since the duality between spaces of measures and spaces of continuous functions is different. Separate convexity, which we use to obtain the necessary conditions for optimality, also appears in [11]. Non-local interaction of the type considered here also appears, for example, in [10]. Here we show that this kind of non-local interaction can be studied also in the present variational setting.

The paper is organized as follows: Section 2 is devoted to introduce the Eulerian variational formulation (in the same flavor as [24, 14]) as well as the analysis of minimizers and the duality. Notice here that since the interaction term is slightly non-convex we have to introduce and study the linearized functional. In Section 3 we introduce the Lagrangian variational formulation which, as in [6], turns out to correspond to a minimization of a relative entropy. Section 4 is devoted to the time discretization of the two formulations and the $\Gamma$-convergence of the discrete problem to the continuous one. The time discretization, as well as the linearization functional introduced in section 2 are then useful in order to introduce a suitable numerical method in section 5 based on generalized Sinkhorn algorithm.

## 2. Variational formulation

Let $\mathcal{P}_2 = \{\mu \in \mathcal{P} \mid \mu \ll \text{Leb}^d, \int |x|^2 d\mu < +\infty\}$ the space of probability measure which are absolutely continuous with respect to the $d$-dimensional Lebesgue measure and have finite second moment (we recall that the ambient space is either $\mathbb{R}^d$ or $\mathbb{T}^d$). We will consider the metric space $(\mathcal{P}_2, \mathcal{W}_2)$ constituted by $\mathcal{P}_2$ equipped with the Wasserstein distance $\mathcal{W}_2$. This is a length metric space (see, for example, [2]) and absolutely continuous curve in this metric space will play a role. Since the functional (1.2) is not convex in the couple $(\rho_t^i, v^i)$ we have to introduce the momentum variable $m_t^i = \rho_t^i v_t^i$ and re-write the functional in the following convex formulation

$$ \mathcal{J}(\rho^1, m^1; \dots; \rho^N, m^N) := \sum_i \int_{[0,T] \times \mathbb{R}^d} \frac{|m_t^i|^2}{2\rho_t^i} dxdt + \sum_{\substack{1 \le j \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V^{i,j}(x-y) d\rho_t^i \otimes d\rho_t^j dt + \int_{\mathbb{R}^d} g^i(x) d\rho^i(T,x); \quad (2.1) $$

if $\rho^j \in \mathcal{P}_2$, $m^i \ll \rho^i$, with a little abuse of notations, $m_t^i = \frac{\partial m_t^i}{\partial \rho_t^i}$ denotes the Radon-Nicodym derivative of measures, and $+\infty$, otherwise. In the previous problem, then, $v^i = \frac{\partial m^i}{\partial \rho^i}$.
We will be interested in

$$ \inf \left\{ \mathcal{J}(\rho^1, m^1; \dots; \rho^N, m^N) \mid \partial_t \rho_t^i - \frac{1}{2} \Delta \rho_t^i + \operatorname{div} m_t^i = 0, \rho^i(0, x) = \rho_0^i, \forall i = 1, \dots, N \right\}. \quad (2.2) $$

Absolutely continuous curves in $\mathcal{P}_2$ are, by now, well characterized. We report here part of the characterization which may be read, for example, in Th. 5.14 of [30]. For a curve $\gamma$ in a metric space $X$ we denote by $|\gamma'|$ the metric derivative of $\gamma$.

**Theorem 2.1.** Let $\rho_t$ be an absolutely continuous curve in $\mathcal{P}_2(\mathbb{R}^d)$ then for a.e. $t \in [0, T]$, there exists a vector field $v_t \in L_{\rho_t}^2$ such that

$$ \partial_t \rho_t + \operatorname{div}(v_t \rho_t) = 0, \quad (2.3) $$

and for a.e. $t$ we have $\|v_t\|_{L_{\rho_t}^2} \le |\rho'|$, where $|\rho'|$ denotes the metric derivative of $\rho$. Conversely, if $\rho_t$ is a curve in $\mathcal{P}_2(\mathbb{R}^d)$ and for a.e. $t \in [0, T]$ there exists $v_t \in L_{\rho_t}^2(\mathbb{R}^d, \mathbb{R}^d)$ such that (2.3) holds, then $\rho_t$ is an absolutely continuous curve in $\mathcal{P}_2$ and for a.e. $t$ we have $|\rho'| \le \|v_t\|_{L_{\rho_t}^2}$.

If $\{\rho_t\}_{t \in [0, T]}$ is a curve of probability measures which solves

$$ \partial_t \rho_t - \frac{1}{2} \Delta \rho_t + \operatorname{div}(v_t \rho_t) = 0, $$

with

$$ \int_0^T \int_{\mathbb{R}^d} |v_t|^2 d\rho_t dt < +\infty $$

and $\rho(0) = \rho_0$ then, by Lemma 4.1 of [7], $\rho_t$ is absolutely continuous since for $w_t := v_t - \frac{1}{2} \nabla \log \rho_t$ it holds

$$ \partial_t \rho - \operatorname{div}(w_t \rho_t) = 0, $$

and

$$ \int_0^T \int |w_t|^2 \rho_t dx dt = \int_0^T \int |v_t|^2 \rho_t dx dt + \frac{1}{4} \int_0^T I(\rho_t) dt + \mathcal{H}(\rho(T)) - \mathcal{H}(\rho_0) < +\infty, $$

where $I$ denotes the Fischer information, that is

$$ I(\rho) = 4 ||\nabla \sqrt{\rho}||_{L^2}^2 $$

Actually Lemma 4.1 of [7] is richer and more detailed than we need and, for more, we refer the reader to the original paper. We stress that the finiteness of the average kinetic energy associated to the curve $\rho_t$ and the vector field $v_t$ implies, by (iv) of Lemma 4.1 in [7], the finiteness of the Fisher information and the entropy of $\rho(T, x)$. In particular, if $(\rho_1, \dots, \rho_N)$ is admissible for problem (2.2) then each component $\rho_i$ satisfies the assumptions above and so it is an absolutely continuous curve in $\mathcal{P}_2$.
Existence of at least one minimiser of problem (2.2) will be proved in the next section. It is well known that the first term of $\mathcal{J}$ is not convex in the couple $(\rho^i, v^i)$. However, as we have mentioned above, it can be rewritten as

$$ \int_{[0,T]\times\mathbb{R}^d} f_2(\rho_t^i, m_t^i)dxdt, $$

where

$$ f_2(t, z) := \begin{cases} \frac{z^2}{t} & \text{if } t > 0, \\ 0 & \text{if } t = 0, \text{ and, } z = 0. \\ +\infty & \text{otherwise.} \end{cases} $$

and the functional $(\rho^i, m^i) \mapsto \int f_2(\rho^i, m^i)$ is convex and lower semi-continuous (see, for example, Th. 5.18 of [30]).

The second term is only separately convex in the $N$-tuple $(\rho^1, \dots, \rho^N)$ so, obtaining the optimality conditions goes through a directional linearisation process which will also be used in the following section.

Let $(\bar{\rho}^1, \bar{m}^1; \dots; \bar{\rho}^N, \bar{m}^N)$ be a minimiser of (2.2), define, for $i = 1, \dots, N$

$$ H^i(t, x) = \sum_{j \neq i} \int (V^{i,j}(x-y) + V^{j,i}(x-y)) d\bar{\rho}_t^j(y), \quad (2.4) $$

and consider the functionals

$$ \mathcal{J}^i(\rho, m) := \int_{[0,T]\times\mathbb{R}^d} \frac{|m_t|^2}{2\rho_t} dxdt + \int_{[0,T]\times\mathbb{R}^d} H^i d\rho_t dt + k^i + \int_{\mathbb{R}^d} g^i(x) d\rho^i(T, x), $$

where the constant is given by

$$ k^i := \sum_{j \neq i} \int_{[0,T]\times\mathbb{R}^d} \frac{|\bar{m}_t^j|^2}{2\bar{\rho}_t^j} dxdt \\ + \int_{\mathbb{R}^d} g^j(x) d\bar{\rho}^j(T, x) + \sum_{\substack{1 \le k \le N \\ j \neq k \neq i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V^{j,k}(x-y) d\rho_t^k \otimes d\rho_t^j dt. $$

The following proposition is slightly more than a remark.

**Proposition 2.2.** For $i = 1, \dots, N$, the couple $(\bar{\rho}^i, \bar{m}^i)$ is a minimizer, among curves starting at $\rho_0^i$, of the functional

$$ \mathcal{J}^i(\rho, m), $$

which is, moreover, a convex functional.
*Proof.* We observe that

$$ \mathcal{J}^i(\rho, m) = \mathcal{J}(\bar{\rho}^1, \bar{m}^1; \dots; \rho, m; \dots; \bar{\rho}^N, \bar{m}^N) \geq \mathcal{J}(\bar{\rho}^1, \bar{m}^1; \dots; \bar{\rho}^N, \bar{m}^N) = \mathcal{J}^i(\bar{\rho}_t^i, \bar{m}_t^i). $$

The convexity of $\mathcal{J}^i$ is due to the separate convexity of the only term of $\mathcal{J}$ which is not convex. $\square$

We will then make a careful analysis of a functional of the type of $\mathcal{J}^i$.

## 2.1. Analysis of the minimizers

In this subsection, for an $H$ which has the properties as the one defined as in (2.4), we will write necessary conditions for the minimization of

$$ \mathcal{F}(\rho, m) := \frac{1}{2} \int_{[0,T] \times \mathbb{R}^d} \frac{|m_t|^2}{\rho_t} dxdt + \int_{[0,T] \times \mathbb{R}^d} H d\rho_t dt + \int_{\mathbb{R}^d} g(x) d\rho(T, x), $$

among solutions of

$$ \left\{ \begin{array}{l} \partial_t \rho_t - \frac{1}{2} \Delta \rho_t + \operatorname{div} m_t = 0; \\ \rho(0, x) = \rho_0(x). \end{array} \right. $$

As frequently in this settings, the convex analysis, via the Fenchel-Rockafellar duality, will be the main tool. In the next two lemmas we will write $\mathcal{F}$ as the sum of two convex conjugates. Since $(\rho, m)$ is couple of finite measures, the natural space for duality is the space of bounded continuous functions $C_b([0, T] \times \mathbb{R}^d) \times C_b([0, T] \times \mathbb{R}^d; \mathbb{R}^d)$. In the following will denote by $C_b^2([0, T] \times \mathbb{R}^d)$ the set of functions which are $C_b^1$ in time and $C_b^2$ in space.

**Lemma 2.3.** Let $(\alpha, \beta) \in C_b([0, T] \times \mathbb{R}^d) \times C_b([0, T] \times \mathbb{R}^d; \mathbb{R}^d)$, $H$ defined as in (2.4) and

$$ \mathcal{K}_1(\alpha, \beta) := \begin{cases} 0 & \text{if } \alpha + \frac{|\beta|^2}{2} \le H, \\ +\infty, & \text{otherwise.} \end{cases} $$

Then

$$ \mathcal{K}_1^*(\rho, m) = \begin{cases} \int \frac{|m|^2}{\rho} + H d\rho & \text{if } m \ll \rho, \\ +\infty & \text{otherwise.} \end{cases} $$

*Proof.* This is part of Prop. 5.18 in [30]. $\square$

**Lemma 2.4.** Let $(\alpha, \beta) \in C_b([0, T] \times \mathbb{R}^d) \times C_b([0, T] \times \mathbb{R}^d; \mathbb{R}^d)$ and

$$ \mathcal{K}_2(\alpha, \beta) = \begin{cases} \int \phi(0, x) \rho_0 dx & \text{if } \alpha = -\partial_t \phi - \frac{1}{2} \Delta \phi \text{ and } \beta = -\nabla \phi, \\ +\infty & \text{otherwise;} \end{cases} $$
for some $\phi \in \mathbb{C}_b^2([0, T] \times \mathbb{R}^d)$, such that $\phi(T, x) \geq -g(x)$. Then

$$ \mathcal{K}_2^*(\rho, m) = \begin{cases} \displaystyle \int_{+\infty} g(x)\rho(T, x)dx & \text{if } \partial_t\rho - \frac{1}{2}\Delta\rho + \text{div } m = 0, \rho(0) = \rho_0 \\ \\ & \text{otherwise;} \end{cases} $$

*Proof.* By the finiteness condition on $\mathcal{K}_2$ we can write

$$ \mathcal{K}_2^*(\rho, m) = \sup\{\langle -\partial_t\phi - \frac{1}{2}\Delta\phi, \rho \rangle - \langle \nabla\phi, m \rangle - \int \phi(0, x)\rho_0(x)dx \mid \\ \phi \in \mathbb{C}_b^2([0, T] \times \mathbb{R}^d), \phi(T, x) \geq -g(x)\} $$

Since we will need to integrate by parts, we first show that the sup can be considered on $\mathbb{C}_0^\infty$ functions. Let $\phi$ be admissible for the sup above and let $\beta$ be a positive constant and $R$ such that $\text{spt}(g) \subset \subset B(0, R)$. The function $\phi^\beta := \phi + \beta$ still satisfies the regularity constraint but, in addition, it satisfies $\phi^\beta(T, x) \geq -g(x) + \beta$. Letting $\beta$ go to 0 the value of

$$ \langle -\partial_t \phi^\beta - \frac{1}{2} \Delta \phi^\beta, \rho \rangle - \langle \nabla \phi^\beta, m \rangle - \int \phi^\beta(0, x) \rho_0(x) dx, $$

tends to the corresponding value for $\phi$ from which it differs only for $\beta$. To adjust the support we consider a smooth cut-off function $\chi_R$ such that $0 \leq \chi_R \leq 1$, $\chi_R \equiv 1$ in $B(0, R)$ and $\chi_R \equiv 0$ in $\mathbb{R}^d \setminus B(0, R)$, the gradient and the Hessian of $\chi_R$ are bounded independently of $R$. We then consider $\phi_R^\beta := \chi_R \phi^\beta$. By the choice of $R$ we have that $\phi_R^\beta \geq -g$ and $\phi_R^\beta \geq -g + \beta$ in a neighborhood of the support of $g$. Letting $R$ go to $+\infty$ the value of

$$ \langle -\partial_t \phi_R^\beta - \frac{1}{2} \Delta \phi_R^\beta, \rho \rangle - \langle \nabla \phi_R^\beta, m \rangle - \int \phi_R^\beta(0, x) \rho_0(x) dx, $$

tends to the corresponding value for $\phi^\beta$ by the dominated convergence theorem. Finally, we need a smoothing of the function. Since it has compact support, $g$ is uniformly continuous. So that we may choose $\delta$ smaller then $\text{dist}(\text{spt}(g), \partial B(0, R))$ and such that if $|x-y| < \delta$ then $|g(x)-g(y)| < \beta$. Choose a smoothing parameter $\varepsilon < \delta$ and a standard spatial mollifier $\varphi_\varepsilon$. Consider $\phi_R^{\beta, \varepsilon} := \varphi_\varepsilon * \phi_R^\beta$. By the uniform convergence, letting $\varepsilon$ go to 0, the value of

$$ \langle -\partial_t \phi_R^{\beta, \varepsilon} - \frac{1}{2} \Delta \phi_R^{\beta, \varepsilon}, \rho \rangle - \langle \nabla \phi_R^{\beta, \varepsilon}, m \rangle - \int \phi_R^{\beta, \varepsilon}(0, x) \rho_0(x) dx, $$

tends to the corresponding value for $\phi_R^\beta$. So we only need to check that $\phi_R^{\beta, \varepsilon}(T, \cdot)$ satisfies the constraint. Let $x$ be in a neighborhood of $\text{spt}g$ such that $\phi_R^\beta > -g + \beta$ then we have,
by the choice of $\delta$

$$ \phi_{R}^{\beta, \varepsilon}(T, x) = \int_{B(x, \varepsilon)} \varphi_{\varepsilon}(y-x) \phi_{R}^{\beta}(y) dy \geq - \int_{B(x, \varepsilon)} \varphi_{\varepsilon}(y-x) g(y) dy + \beta \\ \geq - \int_{B(x, \varepsilon)} \varphi_{\varepsilon}(y-x) (g(x) + \beta) dy + \beta = g(x). $$

In points $x$ far away from $\text{spt}(g)$ the inequality is true since the functions remain non-negative. We may now assume that $\phi$ is regular enough. Integrating by parts gives,

$$ \sup\{\langle\phi, \partial_t \rho - \frac{1}{2}\Delta\rho + \text{div}(m)\rangle - \int \phi(T, x)\rho(T, x)dx + \int \phi(0, x)\rho(0, x)dx - \int \phi_0(x)d\rho_0(x)\}. $$

We consider, now, the sup over suitable subsets of $\{\phi \in C_0^2([0, T] \times \mathbb{R}^d), \phi(T, x) \geq -g(x)\}$. First, let $\hat{\phi}$ be admissible and such that $\hat{\phi}(0, x) \equiv 0$ and $\hat{\phi}(T, x) > -g(x)$. For every $\varphi \in C_c^\infty((0, T) \times \mathbb{R}^d)$ the function $\hat{\phi} + \varphi$ is admissible. Considering the sup over the (affine) space $\hat{\phi} + C_c^\infty$ gives that unless

$$ \partial_t \rho_t - \frac{1}{2} \Delta \rho_t + \text{div} m_t = 0 $$

in the sense of distributions the sup is $+\infty$.

and

$$ \rho(0, x) = \rho_0 $$

the sup is $+\infty$. If these last two conditions are satisfied then

$$ \mathcal{K}_2^*(\rho, \mu) = \int g(x)\rho(T, x)dx. $$

□

By Fenchel-Rockafellar duality

**Proposition 2.5.** *Weak duality holds*

$$ \min\{\mathcal{K}_1^*(\rho, m) + \mathcal{K}_2^*(\rho, m)\} \geq \sup\{-\mathcal{K}_1(\alpha, \beta) - \mathcal{K}_2(-\alpha, -\beta)\} $$

According to the expression of the functionals obtained above

$$ \min\{\mathcal{K}_1^*(\rho, m) + \mathcal{K}_2^*(\rho, m)\} = \min \left\{ \int \frac{|m|^2}{\rho} + \int H d\rho + \int g(x)\rho(T, x)dx \right. \\ \left. \text{if } m \ll \rho \text{ and } \partial_t \rho - \frac{1}{2}\Delta\rho + \text{div} m = 0, \right\} $$
which is our original problem. On the other end the right-hand side is

$$ \sup \left\{ \int \phi(0,x) d\rho_0(x) \mid \phi \in \mathcal{C}_b^2([0,T] \times \mathbb{R}^d), \right. \\ \left. -\partial_t \phi - \frac{1}{2} \Delta \phi + \frac{|\nabla \phi|^2}{2} \le H, \phi(T,x) \le g(x). \right\} \quad (2.5) $$

Concerning the existence of a solution for this last problem, we have

**Theorem 2.6.** There exists a $\mathcal{C}_b^2$ solution $\psi$ of

$$ \left\{ \begin{array}{l} -\partial_t \psi - \frac{1}{2} \Delta \psi + \frac{|\nabla \psi|^2}{2} = H(t,x), \\ \psi(T,x) = g(x) \end{array} \right. \\ \int \psi(0,x) \rho_0(x) dx = (2.5). $$

and such $\psi$ is a maximizer for problem (2.5) above.

*Proof.* Setting $u(t,x) := \psi(T-t,x)$ the problem above is transformed in

$$ \left\{ \begin{array}{l} \partial_t u - \frac{1}{2} \Delta u + \frac{|\nabla u|^2}{2} = H(T - t, x), \\ u(0, x) = g(x) \end{array} \right. $$

Then the Hopf-Cole transform $v(x,t) = e^{-\frac{u}{2}}$ gives the further simplification

$$ \left\{ \begin{array}{l} \partial_t v - \frac{1}{2} \Delta v = -v \frac{1}{2} H(T - t, x), \\ v(0, x) = e^{-\frac{g(x)}{2}} \end{array} \right. $$

So, setting $h(x,t) := -\frac{H(T-t,x)}{2}$ we can look at the solutions of the Heat equation. By assumptions (A1), (A4), Remark 1.1 and Lemma 4.1 of [7], $h$ is continuous and bounded. So that, according also to our assumptions on the functions $g^i$ we can apply Theorem 5.1 of [23] (page 320) which gives a $\mathcal{C}^2$ solution of the problem. Since $H$, defined by (2.4) is bounded by assumption (A4), one get, by the comparison principle $|\psi| \le (1+T)C$ (see also [13]. The maximum principle implies, then, the maximality of $\psi$. $\square$

We conclude the section by connecting everything back to the system (1.1),

**Theorem 2.7.** Let $(\rho^1, m^1, \dots, \rho^N, m^N)$ be a minimizer for (2.2), define $v^i = \frac{\partial m^i}{\partial \rho^i}$, and $+\infty$, otherwise and, for all $i = 1, \dots, N$, let $u^i$ be a maximizer for (2.5) with $H = H^i$ (defined in equation (2.4)), $\rho_0 = \rho_0^i$ and $g = g^i$. Then $v^i = -\nabla u^i$ and $(\rho^1, \dots, \rho^N, u^1, \dots, u^N)$ is a solution of (1.1).
Proof. Let us consider a maximizer $u^i$ to (2.5), it is easy to show that associated to the vector field $-\nabla u^i$, there exists an absolutely continuous curve $\tilde{\rho}^i$ satisfying the Fokker-Planck equation with initial condition $\rho_0^i$. Thus, the couple $(\tilde{\rho}^i, \tilde{m}^i)$ with $\tilde{m}^i = -\nabla u^i \tilde{\rho}^i$ is now a competitor for the primal problem. Indeed, we have, for every minimizer $(\rho_t^i, m^i)$ and by 2.5

$$ 
\begin{aligned}
& \int_{[0,T]\times\mathbb{R}^d} \frac{|\tilde{m}^i|^2}{2\tilde{\rho}_t^i} dxdt + \int_{[0,T]\times\mathbb{R}^d} H^i d\tilde{\rho}_t^i dt + \int_{\mathbb{R}^d} g^i(x) d\tilde{\rho}^i(T,x) \ge \\
& \int_{[0,T]\times\mathbb{R}^d} \frac{|m^i|^2}{2\rho_t^i} dxdt + \int_{[0,T]\times\mathbb{R}^d} H^i d\rho_t^i dt + \int_{\mathbb{R}^d} g^i(x) d\rho^i(T,x) \ge \int_{\mathbb{R}^d} u^i(0,x) d\rho_0^i(x).
\end{aligned} 
\quad (2.6) $$

Since $u^i$ is a maximizer, which saturates the constraint in (2.5) by Theorem 2.6, then

$$ 
\begin{aligned}
& \int_{[0,T]\times\mathbb{R}^d} H^i d\tilde{\rho}_t^i dt + \int_{\mathbb{R}^d} g^i(x) d\tilde{\rho}^i(T,x) = \\
& \int_{[0,T]\times\mathbb{R}^d} \left( -\partial_t u^i - \frac{1}{2}\Delta u^i + \frac{|\nabla u^i|^2}{2} \right) d\tilde{\rho}_t^i + \int_{\mathbb{R}^d} u^i(T,x) d\tilde{\rho}^i(T,x) = \\
& \int_{[0,T]\times\mathbb{R}^d} u^i(\partial_t \tilde{\rho}_t^i - \frac{1}{2}\Delta \tilde{\rho}_t^i) dxdt + \int_{[0,T]\times\mathbb{R}^d} \frac{|\nabla u^i|^2}{2} d\tilde{\rho}_t^i dt + \int_{\mathbb{R}^d} u^i(0,x) d\rho_0^i.
\end{aligned} 
\quad (2.7) $$

where we have integrated by part and exploit the fact that $\tilde{\rho}^i$ solves the Fokker-Plank equation with initial datum $\rho_0^i$ and vector field $-\nabla u^i$. Moreover, the last line turns out to be equal to

$$ \int_{[0,T]\times\mathbb{R}^d} u^i(-\mathrm{div}(\tilde{m}^i))dxdt + \int_{[0,T]\times\mathbb{R}^d} \frac{|\nabla u^i|^2}{2} d\tilde{\rho}_t^i dt + \int_{\mathbb{R}^d} u^i(0,x) d\rho_0^i. $$

Combining with (2.6), integrating again by part the first term in the above equation and that $\tilde{m}_t^i = -\nabla u^i \tilde{\rho}_t^i$ we get

$$ 
\begin{aligned}
& \int_{[0,T]\times\mathbb{R}^d} |\nabla u^i|^2 \tilde{\rho}_t^i dxdt - \int_{[0,T]\times\mathbb{R}^d} |\nabla u^i|^2 \tilde{\rho}_t^i dxdt + \int_{\mathbb{R}^d} u^i(0,x) d\rho_0^i(x) \\
& \ge \int_{[0,T]\times\mathbb{R}^d} \frac{|m^i|^2}{2\rho_t^i} dxdt + \int_{[0,T]\times\mathbb{R}^d} H^i d\rho_t^i dt + \int_{\mathbb{R}^d} g^i(x) d\rho^i(T,x) \\
& \ge \int_{\mathbb{R}^d} u^i(0,x) d\rho_0^i(x).
\end{aligned} 
\quad (2.8) $$

It easily follows now that $v^i = -\nabla u^i$. is optimal for the primal problem. $\square$

<figure><img src="image_4.png" /></figure>
### 3. The Entropic problem

In this section we focus on a Lagrangian formulation to (2.2) based on a minimization of a relative entropy. In particular by proving the existence of a solution to the Lagrangian formulation we deduce the existence also for the Eulerian one.

Let $\Omega := (C([0, T], \mathbb{R}^d), \| \cdot \|_{\infty})$ and let $R^i \in \mathcal{P}(\Omega)$ be Wiener measures defined by

$$R^i = \int_{\mathbb{R}^d} \rho_0^i \delta_{x+B^i(t)(x)} dx$$

where the $B_i$ are the classical Brownian motions. For $Q^i \in \mathcal{P}(\Omega)$ consider the following variational problem to which we will refer as Lagrangian

$$\min\{J(Q^1, \dots, Q^N) : e_{\#}^0 Q^i = \rho_0^i\} \qquad (3.1)$$

where

$$J(Q^1, \dots, Q^N) := \sum_i \mathcal{H}(Q^i|R^i) + \sum_{\substack{1 \le i \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) de_{\#}^t Q^i \otimes de_{\#}^t Q^j dt \\ + \sum_i \int_{\mathbb{R}^d} g^i(x) de_{\#}^T Q^i,$$

where $e_{\#}^t: \mathcal{P}(\Omega) \to \mathcal{P}(\mathbb{R}^d)$ is the evaluation map at time $t$; furthermore, notice that for every $t \in [0, T]$, $e_{\#}^t$ is continuous for the narrow convergences.

**Remark 3.1.** Taking $Q^i = R^i$ gives at least one point at which $J$ is finite. In fact, in that case, the entropy terms are equal to 0. For the other terms, it holds:

$$\int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) de_{\#}^t R^i \otimes de_{\#}^t R^j dt = \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) (p_t * \rho_0^i)(x) (H_t * \rho_0^j)(y) dx dy dt$$

where $H_t$ is the Gaussian kernel, that is

$$H_t(z) := \frac{1}{(2\pi t)^d} \exp \left( -\frac{|z|^2}{2t} \right), \quad t > 0, z \in \mathbb{R}^d.$$

This integral is finite by assumption (A4). Analogously,

$$\int_{\mathbb{R}^d} g^i(x) de_{\#}^T R^i = \int_{\mathbb{R}^d} g^i(x) (H_T * \rho_0^i)(x) dx,$$

which is also finite by the assumptions on $g^i$ and $\rho_0^i$.

**Theorem 3.2.** *The minimum in (3.1) is finite and there exists a minimizer $(\overline{Q}^1, \dots, \overline{Q}^N)$.*

*Proof.* By Remark 3.1 above the infimum of $J$ is finite. Let $\{(Q_n^1, \dots, Q_n^N)\}_{n \in \mathbb{N}}$ be a minimizing sequence, so that $J(Q_n^1, \dots, Q_n^N) \le C$. Since $V$ is bounded from below we
have

$$
\begin{aligned}
& \sum_i \mathcal{H}(Q_n^i | R^i) + \boldsymbol{k} + \sum_i \int_{\Omega} g^i(\omega(T)) \mathrm{d}Q_n^i + \max_{\varphi} \left\{ \sum_i \int_{\mathbb{R}^d} \varphi(0,x) \mathrm{d}\rho_0^i - \int_{\Omega} \varphi(0,\omega(0)) \mathrm{d}Q_n^i \right\} \\
& \le \sum_i \mathcal{H}(Q^i | R^i) + \sum_{\substack{1 \le i,j \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) \mathrm{d}e_{\#}^t Q^i \otimes e_{\#}^t Q^j \mathrm{d}t \\
& \quad + \sum_i \int_{\Omega} g^i(\omega(T)) \mathrm{d}Q^i + \max_{\varphi} \left\{ \sum_i \int_{\mathbb{R}^d} \varphi(0,x) \mathrm{d}\rho_0^i - \int_{\Omega} \varphi(0,\omega(0)) \mathrm{d}Q^i \right\} \\
& \le C,
\end{aligned}
$$

where we remind that $\Omega$ is the set of continuous paths $\mathcal{C}([0, T], \mathbb{R}^d)$. The above inequality implies, since the other addenda of the first term have linear growth

$$ \mathcal{H}(Q_n^i | R^i) \leq C_1, \quad \forall i. $$

Thus (see appendix), up to subsequences, there exists $Q^i$ with

$$ Q_n^i \stackrel{*}{\to} Q^i $$

and since all the terms constituting $J$ are lower semicontinuous, $(Q^1, \dots, Q^N)$ is a minimizer.

The first term of $J$ is the entropy which is lower-semicontinuous (see Lemma 9.4.3 of [2] and the appendix), the second and the third terms are given by continuous functionals composed with the continuous operator $Q \mapsto e_{\#}^t Q$ finally the last term is the sup of continuous functionals and l.s.c. as such. $\square$

From now on, we denote $Q_t^i = e_{\#}^t Q^i \quad \forall t \in [0, T]$.

The next theorem from [7] relates the Lagrangian problem (3.1) to the Eulerian formulation of system (1.1) that we discussed in the previous section. This will allow us to obtain the existence of a minimizer $(\bar{\rho}^1, \bar{m}^1, \dots, \bar{\rho}^N, \bar{m}^N)$ for problem (2.2). Since only the first term of the energy in (2.2) depends on $m$, given a curve $\rho \in \mathcal{C}([0, T]; (\mathcal{P}_2, \mathcal{W}_2))$ we define

$$ \mathcal{E}(\rho) := \inf_v \left\{ \frac{1}{2} \int_0^T \int_{\mathbb{R}^d} |v_t|^2 \mathrm{d}\rho_t \mathrm{d}t : \partial_t \rho_t - \frac{1}{2} \Delta \rho_t + \mathrm{div}(v_t \rho_t) = 0 \right\}, $$

and then, we decompose the minimization (2.2) in a two steps minimization writing it
as

$$
\inf \left\{ \sum_i \mathcal{E}(\rho_t^i) + \sum_{\substack{1 \le j \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V^{i,j}(x-y) d\rho_t^i \otimes d\rho_t^j dt \right. \\
\left. + \int_{\mathbb{R}^d} g^i(x) d\rho^i(T,x) : \rho^i \in \mathcal{C}([0,T]; (\mathcal{P}_2, \mathcal{W}_2)), \rho^i(0,x) = \rho_0^i, \forall i = 1, \dots, N \right\}. \quad (3.2)
$$

This last problem only depends on the curve $\rho$ since the role of $m$ (or $v$) has already been
encoded in $\mathcal{E}$. In an analogous way we can decompose also the minimization problem
(3.1) as

$$
\inf \left\{ \sum_i \mathcal{S}(\rho^i) + \sum_{\substack{1 \le j \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V^{i,j}(x-y) d\rho_t^i \otimes d\rho_t^j dt \right. \\
\left. + \int_{\mathbb{R}^d} g^i(x) d\rho^i(T,x) : \rho^i \in \mathcal{C}([0,T]; (\mathcal{P}_2, \mathcal{W}_2)), \rho^i(0,x) = \rho_0^i, \forall i = 1, \dots, N \right\}, \quad (3.3)
$$

where

$$
S(\mu^i) = \inf\{\mathcal{H}(Q^i | R^i) \mid Q_t^i = \mu_t^i, \forall t \in [0, T], \} \quad (3.4)
$$

so obtaining a problem which depends only on the curve $\mu^i$. Then using the following
result of [7] we have the existence of a minimizer of (2.2) from that of a minimizer of
(3.1).

**Proposition 3.3** (see corollary 4.7 of [7]). If $\rho_0 \in \mathcal{P}_2$ and $\mathcal{H}(\rho_0) < +\infty$ then

$$
S(\rho) = \mathcal{E}(\rho) + \mathcal{H}(\rho_0).
$$

**Theorem 3.4.** There exists *minimizers of problems (3.2) and (2.2).*

Proof. If $\overline{Q} = (\overline{Q}^1, \dots, \overline{Q}^N)$ is a minimizer of problem (3.1), the curve $\overline{Q}_t := (\overline{\rho}^1, \dots, \overline{\rho}^N)$ is a minimizer of problem (3.3). By Proposition 3.3 above, problems (3.2) and (3.3) have the same minimizers (while the minimal values differs by a constant). So $(\overline{\rho}^1, \dots, \overline{\rho}^N)$ is also a minimizer for problem (3.2). If we choose, then, for each $i \in \{1, \dots, N\}$, $v^i$ a vector field as in the definition of $\mathcal{E}(\overline{\rho}^i)$ and consider $\overline{m}^i := v^i \overline{\rho}^i$ we have that $(\overline{\rho}^1, \overline{m}^1, \dots, \overline{\rho}^N, \overline{m}^N)$ is a minimizer for problem (2.2). $\square$

4. Time discretization and Γ-convergence

Before introducing a suitable discretizations we shortly recall the two equivalent prob-
lems we studied above. The main player in those problems is a vector curve of probability
measures

$$
\rho \in (\mathbb{C}([0, T]; (\mathcal{P}_2, \mathcal{W}_2)))^N.
$$
Given a positive integer $K$ the discrete version of the $\rho$ above is a $N$-tuple of $(K+1)$-vectors of probability measures

$$ \boldsymbol{\rho}_K \in \mathcal{P}_2^{K+1} \times \dots \times \mathcal{P}_2^{K+1}. $$

So the $i$-th components $\rho_0^i, \dots, \rho_K^i$ of $\rho_K^i$ is a discrete version of the $i$-th curve. To this $(K+1)$-tuple of probability measures we can associate the piece-wise constant curve

$$ \underline{\rho}^i(t) = \rho_j^i, \quad \text{for } t \in \left[ \frac{jT}{K}, \frac{(j+1)T}{K} \right). $$

The ambient space for the Eulerian versions of the problems was

$$ \mathcal{A}_0 = \{ \rho : [0, T] \to (\mathcal{P}_2)^N \mid \text{abs. cont. and s.t. } \rho(0) = \rho_0 \}. $$

We may first introduce the discretized space

$$ \mathcal{A}_0^K := \{ \boldsymbol{\rho}_K \in (\mathcal{P}_2^{(K+1)})^N : \rho_0^i = \rho_0^i, i = 1, \dots, N \} $$

We say that $\rho_K \to \rho$ as $K \to +\infty$ if for all $i = 1, \dots, N$

$$ \sup_{t \in [0,T]} \mathcal{W}_2(\underline{\rho}^i(t), \rho^i(t)) \to 0. $$

Given $N$ continuous curves of measures $\rho^i \in \mathcal{C}([0, T], (\mathcal{P}_2(\mathbb{R}^d), \mathcal{W}_2))$, that is $\rho^i : t \in [0, T] \mapsto \rho_t^i \in \mathcal{P}_2(\mathbb{R}^d)$, we defined the minimal energy $\mathcal{E}(\rho)$, the minimal entropic cost $\mathcal{S}(\rho)$, as well as the cost

$$ \mathcal{F}(\rho^1, \dots, \rho^N) = \sum_{\substack{1 \le i,j \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) d\rho_t^i \otimes d\rho_t^j dt + \sum_i \int_{\mathbb{R}^d} g^i(x) d\rho^i(T, x). $$

Thus, the minimization problem (2.2) and (3.1) can now be re-written in the following way

$$ \inf \left\{ \sum_{i=1}^{N} \mathcal{E}(\rho^i) + \mathcal{F}(\rho^1, \dots, \rho^N) \mid \rho \in \mathcal{A}_0 \right\}, \quad (4.1) $$

and

$$ \inf \left\{ \sum_{i=1}^{N} \mathcal{S}(\rho^i) + \mathcal{F}(\rho^1, \dots, \rho^N) \mid \rho \in \mathcal{A}_0 \right\}. \quad (4.2) $$

We define the time discretization of (4.1) as

$$ \inf \left\{ \sum_{i=1}^{N} \mathcal{E}^{K}(\boldsymbol{\rho}^i) + F^{K}(\boldsymbol{\rho}^1, \dots, \boldsymbol{\rho}^N) \mid \boldsymbol{\rho}_K \in \mathcal{T}^K \right\}, \quad (4.3) $$
where

$$ \mathcal{E}^K(\boldsymbol{\rho}^i) := \sum_{k=0}^{K-1} \mathcal{E}_{\frac{T}{K}}(\boldsymbol{\rho}_k^i, \boldsymbol{\rho}_{k+1}^i), $$

with

$$ \mathcal{E}_{\frac{T}{K}}(\mu, \nu) := \inf\{\frac{1}{2} \int_{0}^{\frac{T}{K}} \int_{\mathbb{R}^d} |v_t|^2 d\rho_t dt \mid \partial_t \rho - \frac{1}{2} \Delta \rho_t + \text{div}(\rho_t v_t) = 0, \rho_0 = \mu, \rho_{\frac{T}{K}} = \nu\}, $$

and

$$ \mathcal{T}^K := \{(\boldsymbol{\rho}^1, \dots, \boldsymbol{\rho}^N) \in ((\mathcal{P}_2(\mathbb{R}^d), \mathcal{W}_2)^{(K+1)})^N \mid \boldsymbol{\rho}_0^i = \rho_0^i, \forall i = 1, \dots, N\}. $$

In a similar way one can discretize in time the Lagrangian counterpart

$$ \mathcal{S}^K(\boldsymbol{\rho}^i) := \inf\{\mathcal{H}(Q^i | R^i) \mid Q \in \mathcal{P}(\Omega), Q_{j \frac{T}{K}}^i = \boldsymbol{\rho}_j^i, j = 0, \dots, K\} \quad (4.4) $$

as well as the interaction term and the final cost

$$ \mathcal{F}^K(\boldsymbol{\rho}^1, \dots, \boldsymbol{\rho}^N) = \frac{T}{K} \sum_{k=1}^{K-1} \sum_{\substack{1 \le i, j \le N \\ j \ne i}} \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) d\boldsymbol{\rho}_k^i \otimes d\boldsymbol{\rho}_k^j + \sum_{i=1}^{N} \int_{\mathbb{R}^d} g^i(x) d\boldsymbol{\rho}_K^i, \quad (4.5) $$

where, we recall, $\boldsymbol{\rho}^i$ stands now for a vector of measures, that is $\boldsymbol{\rho}^i \in \mathcal{P}_2(\mathbb{R}^d)^{(K+1)}$ which discretize a curve of measures. Notice that (4.4) can be equivalently reformulate as a classical multi-marginal problem; that is for $i = 1, \dots, N$ we have

$$ \mathcal{S}^K(\boldsymbol{\rho}^i) := \inf\{\mathcal{H}(\pi_K^i | R_K^i) \mid \pi_K^i \in \Pi(\rho_0^i, \dots, \rho_K^i)\}, \quad (4.6) $$

where $\Pi(\rho_0^i, \dots, \rho_K^i)$ is the set of probability measures on $(\mathbb{R}^d)^{K+1}$ having $\rho_0^i, \dots, \rho_K^i$ as marginals and

$$ R_K^i := (e^0, e^{\frac{T}{K}}, \dots, e^T)_{\#} R^i. $$

Then the discretized (4.2) takes the form

$$ \inf\{\sum_{i=1}^{N} \mathcal{S}^K(\boldsymbol{\rho}^i) + \mathcal{F}^K(\boldsymbol{\rho}^1, \dots, \boldsymbol{\rho}^N) \mid (\boldsymbol{\rho}^1, \dots, \boldsymbol{\rho}^N) \in \mathcal{T}^K\}. \quad (4.7) $$

**Theorem 4.1.** As $K \to +\infty$,

$$ \sum_{i=1}^{N} \mathcal{E}^{K}(\boldsymbol{\rho}^{i}) + \mathcal{F}^{K}(\boldsymbol{\rho}^{1}, \dots, \boldsymbol{\rho}^{n}) \xrightarrow{\Gamma} \sum_{i}^{N} \mathcal{E}(\rho^{i}) + \mathcal{F}(\rho^{1}, \dots, \rho^{n}). $$

*Proof.* $\Gamma$ – $\limsup$ **inequality** Let $\rho \in \mathcal{A}_0$ be such that

$$ \sum_{i} \mathcal{E}(\rho^i) + \mathcal{F}(\rho) < +\infty. $$
We consider the discretization $\rho_K$ of $\rho$ given by $\rho_j = \rho(\xi_j^K)$ where the times $\xi_j^K \in [j \frac{T}{K}, (j+1) \frac{T}{K})$ for $j = 1, \dots, K-1$ are chosen according to Remark 4.2 below, $\xi_K^0 = 0$ and $\xi_K^K = T$. Since $t \mapsto \rho_t$ is uniformly continuous, the convergence,

$$ \rho_K \to \rho, $$

is verified. We check the convergence of each term of the functional, starting with $\mathcal{F}^K$. Since $t \mapsto \rho_t$ is $W_2$ continuous, it is $w^*$ continuous and since $(x, y) \mapsto V(x-y)$ is lower semi-continuous, the same holds for

$$ t \mapsto \sum_{1 \le i, j \le N \atop i \ne j} \int V(x-y) d\rho_t^i(x) \otimes d\rho_t^j(y). $$

And this last map can be considered as $g$ in Remark 4.2. By the choice of discretizing times

$$ \int g^i(x) d\rho_k^i = \int g^i(x) d\rho^i(T, x). $$

Concerning the energy term, it is enough to study the convergence for each $i$. In the addendum $\mathcal{E}_{\frac{T}{K}}(\rho_k^i, \rho_{k+1}^i)$ we may take $\rho(t-k\frac{T}{K})$ and the corresponding optimal vector field $v$ as test function so obtaining

$$ \mathcal{E}^K(\rho_K^i) = \sum_{k=0}^{K-1} \mathcal{E}_{\frac{T}{K}}(\rho_k^i, \rho_{k+1}^i) \le \mathcal{E}(\rho^i), $$

which, passing to the $\limsup$ as $K \to +\infty$ concludes the proof of the $\Gamma - \limsup$ inequality.

Let $\rho \in \mathcal{A}_0$ and let $\rho_K \to \rho$. We have

$$ \frac{T}{K} \sum_{k=0}^{K-1} \sum_{1 \le i, j \le N \atop i \ne j} \int V(x-y) d\rho_k^i \otimes d\rho_k^j = \int_0^T \sum_{1 \le i, j \le N \atop i \ne j} \int V(x-y) d\underline{\rho}_t^i \otimes d\underline{\rho}_t^j dt. $$

The very definition of convergence $\rho_K \to \rho$ implies that $\underline{\rho}_t^i \otimes \rho_t^j \stackrel{*}{\rightharpoonup} d\rho_t^i \otimes d\rho_t^j$ for all $t, i, j$. Since $V$ is lower semi-continuous

$$ \mu \otimes \nu \mapsto \int V(x-y) d\mu \otimes d\nu, $$

is $w^*$ lower semi-continuous, by the lower semi-continuity of the integral

$$ \liminf_{K \to +\infty} \int_0^T \sum_{1 \le i, j \le N \atop i \ne j} \int V(x-y) d\underline{\rho}_t^i \otimes d\underline{\rho}_t^j dt. $$
Concerning the last addendum of the functional

$$ \int g^i(x) d\boldsymbol{\rho}_K^i \to \int g^i(x) d\rho^i, $$

since $g^i$ is continuous and bounded by assumptions. So let us look at the energy term $\mathcal{E}^K$. This term may be studied as in [7] by usual methods in Calculus of Variations. Also in this case it is enough to study $\mathcal{E}^K(\boldsymbol{\rho}^i)$. Since we may assume that the total energy of $\boldsymbol{\rho}_K$ is bounded, the same holds for $\mathcal{E}^K(\boldsymbol{\rho}^i)$. Let $\tilde{\rho}_K, \tilde{v}_K$ be an almost minimizer for $\mathcal{E}^K(\boldsymbol{\rho}^i)$, i.e. $\tilde{\rho}_K : [0, T] \to \mathcal{P}_2$ is absolutely continuous, $\tilde{\rho}_K(j \frac{T}{K}) = \boldsymbol{\rho}_j^i$ for $j = 0, \dots, K$ $\tilde{v}_{K,t} \in L_{\tilde{\rho}_t}^2,$

$$ \partial_t \tilde{\rho}_K - \frac{1}{2} \Delta \tilde{\rho}_K + \operatorname{div} \tilde{v}_K \tilde{\rho}_K = 0 $$

and

$$ C \geq \mathcal{E}^K(\boldsymbol{\rho}_K^i) \geq \int_0^T \int \frac{|\tilde{v}_k|^2}{2} d\tilde{\rho}_{K,t} dt - \frac{1}{K}. $$

We estimate the total variation of the measure $\tilde{m}_K := \tilde{v}_K \tilde{\rho}_K$ in $[0, T] \times \mathbb{R}^d$. Using

$$ C \geq \int_0^T \int \frac{|\tilde{m}_K|^2}{2\tilde{\rho}_K} dxdt, $$

and the Young's inequality we obtain

$$ \|\tilde{m}_K\|_{[0,T] \times \mathbb{R}^d} \leq C + \|\tilde{\rho}\|_{[0,T] \times \mathbb{R}^d}. $$

Since $\{\tilde{\rho}_K\}$ is also bounded in total variation we have that, up to subsequences

$$ (\tilde{\rho}_K, \tilde{m}_K) \stackrel{*}{\rightharpoonup} (\rho, m), $$

with $(\rho, m)$ solution of

$$ \partial_t \rho_t - \frac{1}{2} \Delta \rho_t + \operatorname{div} m_t = 0. $$

Since $m = v\rho$ for a suitable $v$ and applying the lower semi-continuity Theorem 2.34 of [1], originally from [12], we have,

$$ \liminf_{K \to +\infty} \mathcal{E}^K(\boldsymbol{\rho}_K^i) \geq \liminf_{K \to +\infty} \int_0^T \int \frac{|\tilde{v}_k|^2}{2} d\tilde{\rho}_{K,t} dt = \\ \liminf_{K \to +\infty} \int_0^T \int \frac{|\tilde{m}_k|^2}{2\tilde{\rho}_{K,t}} dxdt \geq \int_0^T \int \frac{|m|^2}{2\rho_t} dxdt \geq \mathcal{E}(\rho^i). $$

□

Remark 4.2. Let $g : [0, T] \to [0, +\infty)$ be a lower-semicontinuous functions and let $K \in \mathbb{N}$.
There exists points $\xi_K^j \in \left[j \frac{T}{K}, (j+1) \frac{T}{K}\right)$ such that

$$ \frac{T}{K} \sum_{j=1}^{K-1} g(\xi_K^j) \rightarrow \int_0^T g(t)dt. $$

## 5. Numerical Approximation

We now present a numerical scheme in order to solve the discretized in time problem (4.7). In particular the scheme is based on a variant of the Sinkhorn algorithm, successfully used to solve many variational problems involving optimal transport [21, 6, 4, 8, 28, 29, 18, 27] and it is an adaptation of the scheme introduced in [7, 3]. More numerical schemes for multi-population MFG have been developed, for example, in [16] based on some previous works of the same authors.

We recall that (4.7) reads as:

$$ \inf \left\{ \sum_{i=1}^{N} S^K(\boldsymbol{\rho}^i) + \mathcal{F}^K(\boldsymbol{\rho}^i, \dots, \boldsymbol{\rho}^i) \mid (\boldsymbol{\rho}^i, \dots, \boldsymbol{\rho}^i) \in \mathcal{F}^K \right\} $$

where $S^K$ is itself defined by (4.6) which is an entropy minimization with multi-marginal constraints.

Denoting $P^k : (\mathbb{R}^d)^{K+1} \rightarrow (\mathbb{R}^d)$ the $k$-th canonical projection we can obviously rewrite (4.7) as an optimization problem over plans $\pi_N^i$ only:

$$ \inf \left\{ \sum_{i=1}^{N} \mathcal{H}(\pi_K^i | R_K^i) + i_{\rho_0}^i (P_{\#}^0 \pi_K^i) + G(P_{\#}^K \pi_K^i) \right. \\ \left. + \mathcal{F}^K(P_{\#} \pi_K^1, \dots, P_{\#} \pi_K^N) : (\pi_K^i)_{i=1}^N \in (\mathcal{P}((\mathbb{R}^d)^{K+1}))^N \right\}, \quad (5.1) $$

where $P_{\#} \pi_K^i = (P_{\#}^0 \pi_K^i, \dots, P_{\#}^K \pi_K^i)$, $G$ is the final cost and is of the form $G(\rho) = \int g d\rho$, and

$$ i_{\rho_0}(\rho) = \begin{cases} 0 & \text{if } \rho = \rho_0 \\ +\infty & \text{otherwise} \end{cases} $$

is the indicator function in the convex analysis sense and is used to enforce the initial condition. We recall that for all $i = 1, \dots, N$, the static reference measure $R_K^i$ in the relative entropy term is defined as follows

$$ R_K^i := (e^0, e^{\frac{T}{K}}, \dots, e^T)_{\#} R^i. $$

Moreover, since we are considering the reversible Wiener measure, it turns out that $R_K^i$
can be decomposed by using the heat kernel as

$$R_K^i := \left( \prod_{k=1}^{K} H_{\frac{T}{K}}(x_k - x_{k-1}) \right) d\rho_0(x_0), \cdots dx_K,$$

where

$$H_t(z) := \frac{1}{(2\pi t)^d} \exp \left( -\frac{|z|^2}{2t} \right), \quad t > 0, z \in \mathbb{R}^d,$$

and $|\cdot|$ denotes the standard euclidean norm. We also need a discretization in space, for instance we use a $M$ grid points to discretize $\mathbb{R}^d$, then $\pi_K^i$ and $R_K^i$ become tensors in $\mathbb{R}^{MN}$. For sake of simplicity we will keep the continuous in space notation, but form now on integral must be understood as finite sums and $x_0, \cdots, x_K$ as $M$ vectors. Notice that thanks to the euclidean norm, the heat kernel $H_t(z)$ can be decomposed as a product along the dimension of the one dimensional kernel, that is

$$H_t(z) = \prod_{j=1}^{d} h_t(z_j),$$

where $h_t(z_j)$ is the heat kernel in dimension one. This implies that, instead of storing a matrix $H_t \in \mathbb{R}^{M \times M}$, one can just store $d$ small matrices belonging to $\mathbb{R}^{\frac{d}{M} \times \frac{d}{M}}$. One can now try to generalize the algorithm introduced in [18] and its multi-marginal variant [6] to the multi-population case in the same flavour as [3]. However, since the interaction term between populations is non-convex, it happens that we are out of the domain of application of Sinkhorn algorithm. A way to overcome this difficulty is through a semi-implicit approach in order to treat the interaction term, that is at step $n+1$ the $i$-th plan $\pi_K^{i,(n+1)}$ is computed as the optimal solution of a linearized problem obtained by injecting the $j$-th, with $j \neq i$, plans $\pi_K^{j,(n)}$ computed at the previous step: for all $i=1, \cdots, N$

$$\pi_K^{i,(n+1)} := \operatorname{argmin}_{\pi_K^i \in \mathcal{P}((\mathbb{R}^d)^{K+1})} \left\{ \begin{aligned} & \mathcal{H}(\pi_K^i | R_K^i) + i_{\rho_0^i} (P_{\#}^0 \pi_K^i) + G(P_{\#}^K \pi_K^i) \\ & + \mathcal{F}_i^K (P_{\#} \pi_K^i) \end{aligned} \right\}, \quad (5.2)$$

where

$$\mathcal{F}_i^K (P_{\#} \pi_K^i) := \mathcal{F}^K (P_{\#} \pi_K^{1,(n+1)}, \cdots, P_{\#} \pi_K^i, \cdots, P_{\#} \pi_K^{N,(n)})$$

We now have to solve $N$ finite-dimensional strictly convex minimization problems. Then, for each problem strong duality holds and (5.2) can be re-written as follows

$$\sup_{(u_0^i, \cdots, u_K^i)} -F_0^*(-u_0^i) - G^*(-u_K^i) - \mathcal{F}_i^*(-u_1^i, \cdots, -u_K^i - 1) \\ - \int \left( \exp \left( \oplus_{j=0}^K u_j^i \right) - 1 \right) R_K^i, \quad (5.3)$$
where with a slight abuse of notation $F_i^*$ denotes the sum of the Fenchel-Legendre transform of each term in $\mathcal{F}_i^K$. Denoting by $\pi_K^{i,(n+1)}$ and $u_j^{i,(n+1)}$ the optimal solution to (5.2) and (5.3), respectively, it follows that the unique solution to (5.2) has the form

$$ \pi_K^{i,(n+1)}(\boldsymbol{x}) := \left( \otimes_{k=0}^{K} e^{u_k^{i,(n+1)}(x_k)} \right) R_K^i(\boldsymbol{x}), $$

where $\boldsymbol{x} = (x_0, \dots, x_K)$.

Remark 5.1 (Structure of the optimal solution). By definition of the linearized term $\mathcal{F}_i^K$ it follows that the $\mathcal{F}_i^*$ is just a sum of indicator function in the convex analysis sense; this implies that for all $k = 1, \dots, K-1$

$$ u_k^{i,(n+1)}(x_k) = \sum_{j \neq i} \int V^{i,j}(x_k - y_k) \rho_k^{j,(n)}(y_k) dy_k, $$

where $\rho_k^{j,(n)}(y_k) := P_{k,\#} \pi_K^{j,(n)}$ is the marginal of the solution computed at the previous step. In the same way if the final cost of the form $G(\rho) = \int gd\rho$ then $u_K^{i,*}(x_K) = g^i(x_K)$. For sake of clarity we consider always these kinds of functional, even if the algorithm can be defined with more complex functional (we refer the reader to [6, 3] for some examples with 1 or 2 populations; the extension to $N$ populations is then straightforward).

Notice that thanks to the remark above the generalised Sinkhorn algorithm takes now the following form.

In the following numerical results we take a space $M \times M$ discretization of $[0, 1]^2$ with $M = 100$ and a time discretization of $[0, 1]$ with $K = 32$ time step. Let us, firstly, consider the 2 densities case: we have always considered the same initial data

$$ \rho_0^1 = \exp(-50(x^1 - .2)^2 - 50(y^1 - .5)^2), \quad \rho_0^2 = \exp(-50(x^2 - .8)^2 - 50(y^2 - .5)^2) $$

and the same final costs

$$ g^1 = 50((x^1 - 0.8)^2 + (y^1 - 0.45)^2), \quad g^2 = 50((x^2 - 0.2)^2 + (y^2 - 0.5)^2), $$

such that we expect the two Gaussians to switch position. As for the interaction potential we have considered in Figure 5.1 a strong repulsion given by $V(x, y) = 120\chi_{||x-y||<0.2}(x, y)$ and in Figure 5.2 a truncated Coulomb repulsion $V(x, y) = \min \left( 1000, \frac{1}{||x-y||} \right)$. Notice in both cases the effect of entropic term which oblige the densities to spread but at the same time the effect of the repulsive interaction forbid them to touch each other (the distance between them depends on the kind of the repulsion term).

It is now straightforward to extend the theory and the numerical method to a slightly
<caption>**Algorithm 1** Multi-population Sinkhorn</caption>

Require: $u_k^{i,(0)} = 0$
1: **while** $\sum_{i=1}^N ||\rho_0^{i,(n)} - \rho_0^i|| < \text{tol}$ **do**
2: &nbsp;&nbsp;&nbsp;**for** $i = 1:N$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $k = 0:K$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $k = 0$ **then**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u_0^{i,(n+1)} = \log(\rho_0^i) - \log\left(\int\left(\exp\left(\bigoplus_{j=1}^K u_j^{i,(n)}\right)\right) R_K^i\right)$
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $k \neq 0, K$ **then**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u_K^{i,(n+1)} = \sum_{j \neq i} \int V^{i,j}(x_k - y_k) \rho_k^{j,(n)}(y_k) dy_k$
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $k = K$ **then**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u_K^{i,(n+1)} = g^i$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**
13: **end while**

<figure><img src="image_5.png" /><figcaption>Figure 5.1: Support of $\rho^1$ and $\rho^2$ for $V(x, y) = 120\chi_{||x-y||<0.2}(x, y)$.</figcaption></figure>

more general model with a viscosity parameter $\epsilon$. The Mean Field Game system (1.1) now takes the form :
<figure><img src="image_6.png" /><figcaption>Figure 5.2: Support of $\rho^1$ and $\rho^2$ for $V(x, y) = \min(1000, \frac{1}{||x - y||})$.</figcaption></figure>

$$ egin{cases} -\partial_t u^i - \frac{\varepsilon}{2}\Delta u^i + \frac{1}{2}|\nabla u^i|^2 = \sum_{j \neq i} \int_{\mathbb{R}^d} V^{i,j}(x-y)\rho_t^j dy, \\ \partial_t \rho_t^i - \frac{\varepsilon}{2}\Delta \rho_t^i + \text{div}(\nabla u^i \rho_t^i) = 0, \\ \rho^i(0,x) = \rho_0^i(x), \quad u^i(T,x) = g^i(x); \end{cases} \tag{5.4} $$

The Lagrangian formulation (3.1) we have proposed becomes

$$ \min\{J_{\varepsilon}(Q^1, \dots, Q^N) : e_{\#}^0 Q^i = \rho_0^i\} \tag{5.5} $$

where

$$ J_{\varepsilon}(Q^1, \dots, Q^N) := \sum_i \mathcal{H}(Q^i | R_{\varepsilon}^i) + \sum_{\substack{1 \le i \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) de_{\#}^t Q^i \otimes de_{\#}^t Q^j dt \\ + \sum_i \int_{\mathbb{R}^d} g^i(x) de_{\#}^T Q^i. $$

where $R_{\varepsilon}^i$ are the reversible Wiener measure induced by a Brownian motion with variance $\varepsilon$. Notice that one could choose different $\varepsilon$ for each population. In particular, if we discretize the problem in time, we have that the reference measure $R_{\varepsilon,K}^i$ can be still decomposed by using the heat kernel $H_{\varepsilon t}(z)$.
Notice that we can still use the algorithm we have introduced above, but the performance, in terms of iterations to converge, will be affected by small values of $\epsilon$. At least formally, when the viscosity is small, (5.5) is an approximation of the following Lagrangian formulation of first-order variational mean-field games (see [5] for the one population case)

$$ \min\{\mathcal{K}(Q^1, \dots, Q^N) : e_{\#}^0 Q^i = \rho_0^i\} \qquad (5.6) $$

where

$$ \mathcal{K}(Q^1, \dots, Q^N) := \sum_i K(Q_i) + \sum_{\substack{1 \le i \le N \\ j \ne i}} \int_0^T \int_{\mathbb{R}^d \times \mathbb{R}^d} V(x-y) d e_{\#}^t Q^i \otimes d e_{\#}^t Q^j dt \\ + \sum_i \int_{\mathbb{R}^d} g^i(x) d e_{\#}^T Q^i. $$

where

$$ K(Q) := \frac{1}{2} \int_{\Omega} \int_{0}^{T} |\dot{\omega}(t)|^{2} dt dQ(\omega). \qquad (5.7) $$

This also implies that we can use the Sinkhorn algorithm, with small $\epsilon$, in order to approximate the solution to first-order MFGs. In Figures 5.3 and 5.4 we have considered the same data as above but with $\epsilon = 0.005$; notice now the effect of a weaker diffusion term which prevents the densities from spreading. Finally, we consider a 3 populations

<figure><img src="image_7.png" /><figcaption>Figure 5.3: Support of $\rho^1$ and $\rho^2$ for $\epsilon = .005$ and $V(x, y) = 120\chi_{||x-y||<0.2}(x, y)$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5.4: Support of $\rho^1$ and $\rho^2$ for $\varepsilon = .005$ and $V(x, y) = \min \left( 1000, \frac{1}{||x - y||} \right)$.</figcaption></figure>

example with initial data

$$ \begin{align*} \rho_0^1 &= \exp(-50(x^1 - .2)^2 - 50(y^1 - .5)^2), \quad \rho_0^2 = \exp(-50(x^2 - .8)^2 - 50(y^2 - .5)^2), \\ \rho_0^3 &= \exp(-80(x^3 - .5)^2 - 80(y^3 - .1)^2), \end{align*} $$

and final costs

$$ g^1 = 50((x^1 - 0.8)^2 + (y^1 - 0.5)^2), \quad g^2 = 50((x^2 - 0.5)^2 + (y^2 - 0.1)^2), \quad g^3 = 50((x^2 - 0.2)^2 + (y^2 - 0.5)^2) $$

which induces a rotation of the populations. In Figure 5.5 we plot the evolution of support of the 3 densities considering as interaction term the strong repulsion we have taken above.

## A. Entropy and the De la Vallée Poussin Theorem

The existence of minimizers for $J$ relies on the compactness of minimizing sequences that follows from the super-linearity of the entropy functional. This is a classical fact which we shortly report here for the reader's convenience.

Let $\mu$ and $\sigma$ be two probability measures on a metric space $X$. We say that $\mu$ is absolutely continuous with respect to $\sigma$ if there exists a function $f \in L_\sigma^1$ such that

$$ \mu = f(x)\sigma, $$
<figure><img src="image_9.png" /><figcaption>Figure 5.5: Support of $\rho^1$, $\rho^2$ and $\rho^3$ for $V(x, y) = 120\chi_{||x-y||<0.2}(x, y)$.</figcaption></figure>

and in this case we write $\mu \ll \sigma$ and we use the classical notation $f = \frac{d\mu}{d\sigma}$.

The relative entropy of $\mu$ with respect to $\sigma$ is defined as

$$ \mathcal{H}(\mu|\sigma) = \begin{cases} \int_X \frac{d\mu}{d\sigma} \log \left( \frac{d\mu}{d\sigma} \right) d\sigma & \text{if } \mu \ll \sigma, \\ +\infty & \text{otherwise.} \end{cases} $$

Following [2] (Example 9.3.6) one can introduce the function

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

which is nonnegative, lower semi-continuous, strictly convex and super-linear at $+\infty$. Then it holds

$$ \mathcal{H}(\mu|\sigma) = \int H\left(\frac{d\mu}{d\sigma}\right) d\sigma; $$

and

$$ \mathcal{H}(\mu|\sigma) = 0 \Leftrightarrow \mu = \sigma. $$

This way to rewrite the relative entropy is handy to make a connection with the De la Vallée Poussin theorem. In fact, let $\{\mu_i\}_{i \in I} \subset \mathcal{P}(X)$ be such that $\mathcal{H}(\mu_i|\sigma) < C$ then the
family $\left\{\frac{d\mu_i}{d\sigma}\right\}_{i \in I} \subset L_{\sigma}^{1}(X)$ is weakly compact thanks to the theorem

**Theorem A.1** (De la Vallée Poussin). Let $\{f_i\}_{i \in I} \subset L_{\sigma}^1$ then the following are equivalent

* the functions $\{f_i\}_{i \in I}$ are uniformly integrable (and then weakly compact in $L_{\sigma}^1$ by the Dunford-Pettis theorem,

* there exists a function $\varphi : \mathbb{R}_+ \to \mathbb{R}_+$ non-decreasing and such that

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

$$ \text{such that} \quad \int_X \varphi(|f_i|)d\sigma < C. $$

A second advantage of writing the entropy using the function $H$ is the lower semi-continuity with respect to the weak $L_{\sigma}^1$-convergence of the densities, in the sense that if

$$ f_n \xrightarrow{w-L_{\sigma}^1} f, $$

then

$$ \liminf_{n \to \infty} \int_X H(f_n)d\sigma \geq \int_X H(f)d\sigma. $$

**Acknowledgment:**

The work of the first author is partially financed by the “Fondi di ricerca di ateneo” of the University of Firenze and partially financed by the EU-Next Generation EU, (Missione 4, Componente 2, Investimento 1.1 Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN), CUPB53D23009310006 - (2022J4FYNJ). The first author is member of the research group GNAMPA of INdAM. The second author is partially on academic leave at Inria (team Matherials) for the year 2023-2024 and acknowledges the hospitality of this institution during this period. His work benefited from the support of the FMJH Program PGMO, from H-Code, Université Paris-Saclay and from the ANR project GOTA (ANR-23-CE46-0001).

Both author warmly thanks the hospitality of *equipe MoKaPlan of the INRIA Paris*, where this paper has been written.

**References**

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

[2] L. Ambrosio, N. Gigli, and G. Savaré. *Gradient Flows in Metric Spaces and in the Space of Probability Measures*. Lectures in Mathematics. Birkhäuser Boston, 2nd edition, 2004.
[3] César Barilla, Guillaume Carlier, and Jean-Michel Lasry. A mean field game model for the evolution of cities. *Journal of Dynamics & Games*, 8(3):299, 2021.

[4] J.-D. Benamou, G. Carlier, and L. Nenna. A numerical method to solve multi-marginal optimal transport problems with coulomb cost. In *Splitting Methods in Communication, Imaging, Science, and Engineering*, pages 577-601. Springer International Publishing, 2016.

[5] J.-D. Benamou, G. Carlier, and F. Santambrogio. Variational mean field games. In *Active Particles, Volume 1*, pages 141-171. Springer, 2017.

[6] Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré. Iterative Bregman projections for regularized transportation problems. *SIAM J. Sci. Comput.*, 37(2):A1111–A1138, 2015.

[7] Jean-David Benamou, Guillaume Carlier, Simone Di Marino, and Luca Nenna. An entropy minimization approach to second-order variational mean-field games. *Mathematical Models and Methods in Applied Sciences*, 29(08):1553–1583, 2019.

[8] Jean-David Benamou, Guillaume Carlier, and Luca Nenna. Generalized incompressible flows, multi-marginal transport and sinkhorn algorithm. to appear on *Numerische Mathematik*, 2018.

[9] Alain Bensoussan, Tao Huang, and Mathieu Lauriere. Mean field control and mean field game models with several populations. *arXiv preprint arXiv:1810.00783*, 2018.

[10] Chiara Bernardini and Annalisa Cesaroni. Ergodic mean-field games with aggregation of choquard-type. *Journal of Differential Equations*, 364:296-335, 2023.

[11] Ariela Briani and Pierre Cardaliaguet. Stable solutions in potential mean field game systems. *Nonlinear Differential Equations and Applications NoDEA*, 25(1):1, 2018.

[12] Giuseppe Buttazzo and Lorenzo Freddi. Functionals defined on measures and applications to non equi-uniformly elliptic problems. *Annali di matematica pura ed applicata*, 159(1):133-149, 1991.

[13] Pierre Cardaliaguet. A short course on mean field games. *Unpublished lecture notes*, 2018.

[14] Pierre Cardaliaguet, P Jameson Graber, Alessio Porretta, and Daniela Tonon. Second order mean field games with degenerate diffusion and local coupling. *Nonlinear Differential Equations and Applications NoDEA*, 22(5):1287-1317, 2015.

[15] Pierre Cardaliaguet and Alessio Porretta. An introduction to mean field game theory. In *Mean Field Games: Cetraro, Italy 2019*, pages 1-158. Springer, 2021.

[16] Elisabetta Carlini and Francisco J Silva. A fully-discrete scheme for systems of nonlinear fokker-planck-kolmogorov equations. In *PDE models for multi-agent phenomena*, pages 195-218. Springer, 2018.
[17] René Carmona, François Delarue, et al. *Probabilistic theory of mean field games with applications I-II*. Springer, 2018.

[18] L. Chizat, G. Peyré, B. Schmitzer, and F.-X. Vialard. Scaling algorithms for unbalanced transport problems. Technical report, http://arxiv.org/abs/1607.05816, 2016.

[19] Marco Cirant. Multi-population mean field games systems with neumann boundary conditions. *Journal de Mathématiques Pures et Appliquées*, 103(5):1294-1315, 2015.

[20] Marco Cirant and Gianmaria Verzini. Bifurcation and segregation in quadratic two-populations mean field games systems. *ESAIM: Control, Optimisation and Calculus of Variations*, 23(3):1145-1177, 2017.

[21] M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In *Advances in Neural Information Processing Systems (NIPS) 26*, pages 2292-2300, 2013.

[22] Gokce Dayanikli and Mathieu Lauriere. Multi-population mean field games with multiple major players: Application to carbon emission regulations. *arXiv preprint arXiv:2309.16477*, 2023.

[23] Olga Aleksandrovna Ladyzhenskaia, Vsevolod Alekseevich Solonnikov, and Nina N. Ural'tseva. *Linear and quasi-linear equations of parabolic type*, volume 23. American Mathematical Soc., 1988.

[24] Jean-Michel Lasry and Pierre-Louis Lions. Mean field games. *Japanese journal of mathematics*, 2(1):229-260, 2007.

[25] Elliott H Lieb. Density functionals for coulomb systems. *International Journal of Quantum Chemistry*, 24(3):243-277, 1983.

[26] Alpár Richárd Mészáros and Francisco J Silva. On the variational formulation of some stationary second-order mean field games systems. *SIAM Journal on Mathematical Analysis*, 50(1):1255-1277, 2018.

[27] Luca Nenna. *Numerical methods for multi-marginal optimal transportation*. PhD thesis, PSL Research University, 2016.

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

[29] Gabriel Peyré and Marco Cuturi. Computational optimal transport. Technical report, 2017.

[30] Filippo Santambrogio. *Optimal transport for applied mathematicians*. Progress in Nonlinear Differential Equations and their Applications, 87. Birkhäuser/Springer, Cham, 2015. Calculus of variations, PDEs, and modeling.
[31] Li-Hsien Sun. Mean field games with heterogeneous groups: Application to banking systems. *Journal of Optimization Theory and Applications*, 192(1):130–167, 2022.