# A LAGRANGE-GALERKIN SCHEME FOR FIRST ORDER MEAN FIELD GAMES SYSTEMS

E Carlini, Francisco José Silva, Ahmad Zorkot

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

## To cite this version:

E Carlini, Francisco José Silva, Ahmad Zorkot. A LAGRANGE-GALERKIN SCHEME FOR FIRST ORDER MEAN FIELD GAMES SYSTEMS. 2023. hal-04046845

HAL Id: hal-04046845

https://hal.science/hal-04046845

Preprint submitted on 26 Mar 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A LAGRANGE-GALERKIN SCHEME FOR FIRST ORDER MEAN FIELD GAMES
SYSTEMS

ELISABETTA CARLINI<fnref n="1" />, FRANCISCO J. SILVA<fnref n="2" />, AND AHMAD ZORKOT<fnref n="3" />

ABSTRACT. In this work, we consider a first order mean field games system with non-local couplings. A Lagrange-Galerkin scheme for the continuity equation, coupled with a semi-Lagrangian scheme for the Hamilton-Jacobi-Bellman equation, is proposed to discretize the mean field games system. The convergence of solutions to the scheme towards a solution to the mean field game system is established in arbitrary space dimensions. The scheme is implemented to approximate two mean field games systems in dimension one and two.

AMS-Subject Classification: 91A16, 49N80, 35Q89, 65M12.

**Keywords:** First order mean field games, Lagrange-Galerkin schemes, semi-Lagrangian schemes, convergence results, numerical experiences.

# 1. INTRODUCTION

In view of its applications in Economics, Physics, and Social Sciences, the study of optimal control problems and differential games with a large number of agents has attracted the attention of several researchers during the last two decades. An important step in this direction has been achieved with the introduction of the theory of Mean Field Games (MFGs) by J.-M. Lasry-Lions [38, 39, 40] and, independently, by M. Huang, R.P. Malhamé, and P.E. Caines [37]. The main purpose of this theory is to characterize Nash equilibria for a class of symmetric differential games with a continuum of agents. One of the main applications of MFGs theory is that such equilibria can be used to provide approximate equilibria for the corresponding games with a large, but finite, number of players. In its standard form, MFGs are described by a system of two Partial Differential Equations (PDEs); a Hamilton-Jacobi-Bellman (HJB) equation, describing the optimal cost of a typical player in the game, and a Fokker-Planck (FP) equation, describing the evolution of the initial distribution when all the players act optimally. We refer the reader to the monographs [33, 23, 24], the survey [34], and the lectures [4] for a throughout overview on MFGs.

The numerical approximation of MFGs with nonlocal couplings has been an active research topic in recent years. In the case where the MFGs system includes nondegenerate second order terms, finite-difference schemes have been studied in [3, 35, 1, 6, 7], semi-Lagrangian scheme where investigated in [22], and machine learning methods such as deep learning and reinforcement learning have been analyzed in [25, 26, 10]. In the case where the dynamics of the underlying differential games are deterministic, the resulting MFGs system is of first order and several numerical methods have been proposed to approximate its solutions; see e.g. [20, 17] for semi-Lagrangian discretizations, [36, 31] for the approximation by discrete-time finite state space MFGs (see [32]), and [44, 42] for Fourier analysis techniques. We refer the reader to [5, 41], and the references therein, for an overview on numerical methods to approximate MFGs equilibria including also the case of local couplings and variational methods.

<footnote n="1">¹Dipartimento di Matematica Guido Castelnuovo, Sapienza Università di Roma, 00185 Rome, Italy (carlini@mat.uniroma1.it).</footnote>

<footnote n="2">²Institut de recherche XLIM-DMI, UMR-CNRS 7252, Faculté des Sciences et Techniques, Université de Limoges, 87060 Limoges, France (francisco.silva@unilim.fr).</footnote>

<footnote n="3">³Institut de recherche XLIM-DMI, UMR-CNRS 7252, Faculté des Sciences et Techniques, Université de Limoges, 87060 Limoges, France (ahmad.zorkot@unilim.fr).</footnote>
In this paper we focus our attention on the approximation of first order MFGs systems. Namely, we consider the PDE system

$$ \begin{array}{l} -\partial_t v + H(x, D_x v) = F(x, m(t)) \quad \text{in } ]0, T[ \times \mathbb{R}^d, \\ \text{(MFG)} \qquad \qquad \qquad \qquad \begin{aligned} v(T, x) &= G(x, m(T)) \quad \text{in } \mathbb{R}^d, \\ \partial_t m - \operatorname{div}(D_p H(x, D_x v)m) &= 0 \quad \text{in } \mathbb{R}^d \times ]0, T[, \\ m(0) &= m_0^* \end{aligned} \end{array} $$

where $H: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ is convex with respect to its second argument and, denoting by $\mathcal{P}_1(\mathbb{R}^d)$ the set of probability measures over $\mathbb{R}^d$ with finite first order moment, $F: \mathbb{R}^d \times \mathcal{P}_1(\mathbb{R}^d) \to \mathbb{R}^d$, $G: \mathbb{R}^d \times \mathcal{P}_1(\mathbb{R}^d) \to \mathbb{R}^d$, and $m_0^* \in \mathcal{P}_1(\mathbb{R}^d)$. In the article [17], the authors propose a convergent semi-discrete scheme to approximate solutions to (MFG). A fully-discrete version has been proposed in [20]. In the proposed scheme, the HJB equation is discretized by using a semi-Lagrangian approximation of the HJB equation (see e.g. [29]), while the FP equation, or continuity equation, is approximated by a scheme which is dual to a linearized version of the scheme for the HJB equation. The existence of solutions to this approximation is shown and a convergence result to a solution to (MFG) is established when the dimension $d$ of the space variable is equal to one. An extension of this scheme to the case where (MFG) involves non-local and fractional diffusions terms has been studied in [27]. If the resulting system has non-smooth solutions, the convergence of solutions of the scheme is also shown when the space dimension is equal to one.

In order to obtain a convergent scheme for general state dimensions, the key point is to provide a scheme which preserves, under standard conditions on the data (see Section 2.2 below), the main properties of solutions to both equations in (MFG). Namely, the boundedness, Lipschitzianity, and semiconcavity of the solution to the HJB equation and a uniform compact support, equicontinuity, and uniform bounds in $L^p$ spaces for solutions to the continuity equation. As shown in [2, 20, 27], a standard semi-Lagrangian scheme for the HJB equation enjoys the former properties under suitable assumptions on the discretization parameters. In order to treat the continuity equation, we consider the Lagrange-Galerkin (LG) scheme introduced in [43] and recalled in Section 4 below. As we show, it turns out that, for a specific choice of the basis functions, the resulting scheme for the continuity equation coincides with the one introduced in [45] and further studied in [48] for Lipschitz velocity fields. The desired properties for the solutions to this scheme are established in Section 4.2. In particular, we provide a uniform $L^p$ estimate, not available in the schemes considered in [20, 27] in arbitrary space dimensions, which will play a key role in our main convergence result. Combining the semi-Lagrangian scheme for the HJB equation and the LG scheme for the continuity equation, we obtain a discretization of (MFG) for which the existence of solutions is established and, using stability and compactness arguments, the convergence to a solution to (MFG) is established.

The rest of this article is organized as follows. In Section 2 we fix some standard notation, and we state our main assumptions on the data of (MFG). Some important results about solutions to HJB and continuity equations are recalled, as well as existence and uniqueness results for solutions to (MFG). The next two sections deal with the discretization of the HJB and continuity equations in (MFG) separately. Section 3 recalls a standard semi-Lagrangian scheme to approximate the solution to the HJB equation in (MFG). Several important properties of this scheme are reviewed and a new semiconcavity estimate for the solution to the scheme is provided in Proposition 3.5. This estimate will play a crucial role in Section 4, which is devoted to the study of a LG scheme to approximate the continuity equation in (MFG). Notice that, in general, this continuity equation is driven by a non-smooth velocity field. We show that the solutions to the LG scheme inherit the equicontinuity and $L^p$-stability of the solution to the original equation and we establish in Proposition 4.6 a convergence result as the discretization steps tend to zero. In Section 5 we couple the schemes studied in the previous sections to obtain a discretization of (MFG). The existence of a solution to the discretized MFG system is provided in Theorem 5.1 and the convergence result, valid in arbitrary dimensions, is shown in Theorem 5.2. Finally, Section 6 is devoted to the numerical implementation of the scheme for the MFGs system. Since the LG scheme for the continuity equation involves some integrals depending on the discrete characteristics of the equation, we approximate them by numerical quadrature and by the so-called area-weighting method introduced in [43]. The performances of these two approximations are compared in a one-dimensional example with an explicit solution, and the area-weighting method is implemented to approximate the solution to a MFGs in a two-dimensional space.
**Acknowledgements.** F. J. Silva and A. Zorkot where partially supported by l'Agence Nationale de la Recherche (ANR), project ANR-22-CE40-0010. For the purpose of open access, the authors have applied a CC-BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.

The three authors were partially supported by KAUST through the subaward agreement ORA-2021-CRG10-4674.6.

## 2. PRELIMINARIES

**2.1. Notation.** Let $d \in \mathbb{N}$. In what follows, $\langle \cdot, \cdot \rangle$ and $| \cdot |$ denote the standard scalar product in $\mathbb{R}^d$ and its induced norm, respectively. We set $| \cdot |_{\infty}$ for the maximum norm in $\mathbb{R}^d$ and $B_{\infty}(0, C)$ and $\bar{B}_{\infty}(0, C)$ for the associated open and closed balls, centered at 0 and of radius $C > 0$, respectively. Let $\mathcal{P}(\mathbb{R}^d)$ be the set of probability measures on $\mathbb{R}^d$. For every $\nu \in \mathcal{P}(\mathbb{R}^d)$ we denote by supp($\nu$) its support. Let $\mathcal{P}_1(\mathbb{R}^d) = \{\nu \in \mathcal{P}(\mathbb{R}^d) | \int_{\mathbb{R}^d} |x| d\nu(x) < \infty\}$, and, for every $\nu_1, \nu_2 \in \mathcal{P}_1(\mathbb{R}^d)$, set

$$ (2.1) \qquad d_1(\nu_1, \nu_2) = \inf_{\gamma \in \Pi(\nu_1, \nu_2)} \int_{\mathbb{R}^d \times \mathbb{R}^d} |x-y| d\gamma(x,y), $$

where $\Pi(\nu_1, \nu_2)$ denotes the set of probabilities measures on $\mathbb{R}^d \times \mathbb{R}^d$ with first and second marginals given by $\nu_1$ and $\nu_2$, respectively. By the Kantorovich-Rubinstein theorem (see e.g. [9, Section 7.1]) we have

$$ (2.2) \qquad d_1(\nu_1, \nu_2) = \sup \left\{ \int_{\mathbb{R}^d} \varphi(x) d(\nu_1 - \nu_2)(x) \mid \varphi \in \text{Lip}_1(\mathbb{R}^d) \right\}, $$

where $\text{Lip}_1(\mathbb{R}^d)$ denotes the set of all nonexpansive functions on $\mathbb{R}^d$. Given $\nu \in \mathcal{P}(\mathbb{R}^d)$ and a Borel function $\Psi: \mathbb{R}^d \to \mathbb{R}^q$ ($q \in \mathbb{N}$), the push-forward measure $\Psi \sharp \nu$, defined on the $\sigma$-algebra of Borel sets $\mathcal{B}(\mathbb{R}^q)$, is defined by

$$ (2.3) \qquad \Psi \sharp \nu(A) = \nu(\Psi^{-1}(A)) \quad \text{for all } A \in \mathcal{B}(\mathbb{R}^q), $$

or, equivalently (see e.g. [14, Theorem 3.6.1]), for every $\varphi: \mathbb{R}^d \to \mathbb{R}$ such that $\varphi \circ \Psi$ is integrable with respect to $\nu$, one has

$$ (2.4) \qquad \int_{\mathbb{R}^q} \varphi(x) d(\Psi \sharp \nu)(x) = \int_{\mathbb{R}^d} \varphi(\Psi(x)) d\nu(x). $$

**2.2. Assumptions.** Our hypothesis on the data of (MFG) are the following:

(H1) It holds that

$$ (2.5) \qquad H(x,p) = \sup_{a \in \mathbb{R}^d} (\langle a, p \rangle - L(x,a)) \quad \text{for all } x, p \in \mathbb{R}^d, $$

where $L: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ is of class $C^2$, bounded from below, and, for every $x, a \in \mathbb{R}^d$, we have

$$ (2.6) \qquad L(x,a) \leq C_{L,1}|a|^2 + C_{L,2}, $$

$$ (2.7) \qquad |D_x L(x,a)| \leq C_{L,3}(1 + |a|^2), $$

$$ (2.8) \qquad C_{L,4}|b|^2 \leq \langle D_{aa}^2 L(x,a)b, b \rangle \quad \text{for all } b \in \mathbb{R}^d, $$

$$ (2.9) \qquad \langle D_{xx}^2 L(x,a)y,y \rangle \leq C_{L,5}(1 + |a|^2)|y|^2 \quad \text{for all } y \in \mathbb{R}^d, $$

for some constants $C_{L,i} > 0$ ($i = 1, \dots, 5$).

(H2) The functions $F$ and $G$ are continuous and, for every $x, y \in \mathbb{R}^d$ and $\nu \in \mathcal{P}_1(\mathbb{R}^d)$, we have

$$ (2.10) \qquad |F(x,\nu)| \leq C_{F,1}, $$

$$ (2.11) \qquad |G(x,\nu)| \leq C_{G,1}, $$

$$ (2.12) \qquad |F(x,\nu) - F(y,\nu)| \leq C_{F,2}|x - y|, $$

$$ (2.13) \qquad |G(x,\nu) - G(y,\nu)| \leq C_{G,2}|x - y|, $$

$$ (2.14) \qquad F(x+y,\nu) - 2F(x,\nu) + F(x-y,\nu) \leq C_{F,3}|y|^2, $$

$$ (2.15) \qquad G(x+y,\nu) - 2G(x,\nu) + G(x-y,\nu) \leq C_{G,3}|y|^2, $$

for some constants $C_{F,i} > 0, C_{G,i} > 0$ ($i = 1, 2, 3$).
(H3) The initial condition $m_0^*$ is absolutely continuous with respect to the Lebesgue measure and satisfies:

(i) There exists $C^* > 0$ such that $\text{supp}(m_0^*) \subset \bar{B}_\infty(0, C^*)$.

(ii) There exists $\mathbf{p} \in ]1, \infty[$ such that the density of $m_0^*$, still denoted by $m_0^*$, belongs to $L^p(\mathbb{R}^d)$.

**Remark 2.1.** Since $L$ is bounded from below, the strong convexity assumption (2.8) on $L(x, \cdot)$, which is uniform with respect to $x \in \mathbb{R}^d$, and (2.6), imply the existence of $C_{L,6} > 0$ and $C_{L,7} > 0$ such that

$$ (2.16) \qquad L(x, a) \geq C_{L,6} |a|^2 - C_{L,7} \quad \text{for all } x, a \in \mathbb{R}^d. $$

It follows from (2.5), (2.6), and (2.16), that there exist $C_{H,i} > 0$ ($i = 1, 2, 3, 4$) such that

$$ (2.17) \qquad C_{H,1}|p|^2 - C_{H,2} \leq H(x,p) \leq C_{H,3}|p|^2 + C_{H,4} \quad \text{for all } x, p \in \mathbb{R}^d. $$

Moreover, by (2.5), (2.16), and Danskin's theorem (see e.g. [15, Theorem 4.13]), we deduce that $H$ is of class $C^1$ and, for every $x, p \in \mathbb{R}^d$, the following equalities hold

$$ (2.18) \qquad D_a L(x, D_p H(x, p)) = p, $$

$$ (2.19) \qquad D_x H(x, p) = -D_x L(x, D_p H(x, p)). $$

Since $D_p H(x, p)$ is the unique maximizer of $\sup_{a \in \mathbb{R}^d} (\langle a, p \rangle - L(x, a))$, (2.6), and (2.16), yield the existence of $C_{H,5} > 0$ such that

$$ (2.20) \qquad |D_p H(x, p)| \leq C_{H,5}(1 + |p|) \quad \text{for all } x, p \in \mathbb{R}^d. $$

Finally, since $L$ is of class $C^2$, by (2.8) and the implicit function theorem applied to (2.18), it follows that $D_p H$ is of class $C^1$ and hence, by (2.19), we obtain that $H$ is of class $C^2$.

A typical example of a function $H$ satisfying (H1) is given by $H(x, p) = a(x)|p|^2 + \langle b(x), p \rangle$, where $a: \mathbb{R}^d \to \mathbb{R}$ is of class $C^2$, with bounded first and second order derivatives, there exist $\underline{a}, \bar{a} \in ]0, \infty[$ such that $\underline{a} \le a(x) \le \bar{a}$ for all $x \in \mathbb{R}^d$, and $b: \mathbb{R}^d \to \mathbb{R}^d$ is bounded, of class $C^2$, with bounded first and second order derivatives.

**2.3. The first order mean field games system.** Given $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, consider the HJB equation

$$ (2.21) \qquad \begin{aligned} & -\partial_t v(t, x) + H(x, D_x v(t, x)) = F(x, \mu(t)) && \text{for } (t, x) \in ]0, T[ \times \mathbb{R}^d, \\ & v(x, T) = G(x, \mu(T)) && \text{for } x \in \mathbb{R}^d. \end{aligned} $$

It follows from [12, 28] that (2.21) admits a unique viscosity solution $v[\mu]$ and, for every $t \in [0, T[$, $x \in \mathbb{R}^d$, and $\alpha \in L^2([t, T]; \mathbb{R}^d)$, setting $X^{t,x,\alpha}(\cdot) = x - \int_t^{\cdot} \alpha(s)ds$ and

$$ (2.22) \qquad J^{t,x}[\mu](\alpha) = \int_t^T \left( L(X^{t,x,\alpha}(s), \alpha(s)) + F(X^{t,x,\alpha}(s), \mu(s)) \right) ds + G(X^{t,x,\alpha}(T), \mu(T)), $$

we have

$$ (2.23) \qquad v[\mu](t, x) = \inf \left\{ J^{t,x}[\mu](\alpha) \mid \alpha \in L^2([t, T]; \mathbb{R}^d) \right\}. $$

The proof of the following result follows from standard arguments (see e.g. [18]). However, for the sake of completeness, we provide its proof in the appendix of this work.

**Proposition 2.1.** Assume (H1)-(H2) and let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. Then the following hold:

(i) [Existence of an optimal control] For every $(t, x) \in [0, T[ \times \mathbb{R}^d$, there exists $\alpha^{t,x} \in L^\infty([t, T]; \mathbb{R}^d)$ such that $v[\mu](t, x) = J^{t,x}[\mu](\alpha^{t,x})$. Moreover, there exists $C_b > 0$, independent of $(\mu, t, x)$, such that $\|\alpha^{t,x}\|_{L^\infty([0,T];\mathbb{R}^d)} \le C_b$.

(ii) [Uniform bound] We have

$$ (2.24) \qquad |v[\mu](t, x)| \leq C_v \quad \text{for all } (t, x) \in [0, T] \times \mathbb{R}^d, $$
where $C_v > 0$ is independent of $\mu$.

(iii) [Lipschitz property] We have

$$ (2.25) \qquad |v[\mu](t, x) - v[\mu](t, y)| \leq C_{\text{Lip}} |x - y| \quad \text{for all } t \in [0, T], x, y \in \mathbb{R}^d, $$
where $C_{\text{Lip}} > 0$ is independent of $\mu$.
(iv) [Semi-concavity] We have

$$
(2.26) \quad v[\mu](t, x+y) - 2v[\mu](t, x) + v[\mu](t, x-y) \le C_{sc}|y|^2 \quad \text{for all } t \in [0, T], x, y \in \mathbb{R}^d,
$$
where $C_{sc} > 0$ is independent of $\mu$.

**Remark 2.2.** Assertion (i) in Proposition 2.1 implies that, for every $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, we have

$$
(2.27) \quad v[\mu](t,x) = \inf \left\{ J^{t,x}[\mu](\alpha) \mid \alpha \in L^{\infty}([0,T]; \mathbb{R}^d), \|\alpha\|_{L^{\infty}([0,T]; \mathbb{R}^d)} \le C_b \right\}
$$

for all $(t, x) \in [0, T[ \times \mathbb{R}^d$. In particular, $v[\mu]$ is also characterized by the HJB equation

$$
-\partial_t v(t, x) + H_b(x, D_x v(t, x)) = F(x, \mu(t)) \quad \text{for } (t, x) \in ]0, T[ \times \mathbb{R}^d,
$$

$$
(2.28) \qquad v(x, T) = G(x, \mu(T)) \quad \text{for } x \in \mathbb{R}^d,
$$

where

$$
(2.29) \quad H_b(x, p) = \sup_{a \in \bar{B}(0, C_b)} \{\langle a, p \rangle - L(x, a)\} \quad \text{for all } x, p \in \mathbb{R}^d.
$$

Consider the set-valued map $D_x^+ v[\mu]: [0, T] \times \mathbb{R}^d \to 2^{\mathbb{R}^d}$ defined by

$$
D_x^+ v[\mu](t, x) = \left\{ p \in \mathbb{R}^d \mid \limsup_{y \to x} \frac{v[\mu](t, y) - v[\mu](t, x) - \langle p, y-x \rangle}{|y-x|} \le 0 \right\} \quad \text{for all } (t,x) \in [0,T] \times \mathbb{R}^d.
$$

It follows from Proposition 2.1(iv) and [18, Proposition 3.1.5 and Proposition 3.3.4] that $D_x^+ v[\mu]$ takes nonempty and closed values and its graph is closed. In particular, since Proposition 2.1(iv) and [18, Theorem 3.3.6] imply that $D_x^+ v[\mu](t, x) \subset \bar{B}(0, C_{\text{Lip}})$ for all $(t, x) \in [0, T] \times \mathbb{R}^d$, by [11, Chapter 1, Corollary 1] we have that $D_x^+ v[\mu]$ is upper-semicontinuous, i.e. for every $M \subset \mathbb{R}^d$ closed, $D_x^+ v[\mu]^{-1}(M)$ is closed. Therefore, $D_x^+ v[\mu]$ is a Borel measurable set-valued map and hence admits a Borel measurable selection (see e.g. [47, Corollary 14.6]). Notice that Proposition 2.1(iii), Rademacher's theorem, and [18, Proposition 3.1.5] imply that all the measurable selections of $D_x^+ v[\mu]$ coincide almost everywhere in $[0, T] \times \mathbb{R}^d$ and hence, hereafter, we will denote likewise by $D_x v[\mu]$ any choice among them.

Let $\mathbf{p} \in ]1, \infty[$ be as in (H3). We say that $m \in L^p([0, T] \times \mathbb{R}^d)$ solves the continuity equation

$$
\partial_t m - \operatorname{div} (D_p H(x, D_x v[\mu]) m) = 0 \quad \text{in } ]0, T[ \times \mathbb{R}^d,
$$

$$
(2.30) \qquad m(0) = m_0^* \quad \text{in } \mathbb{R}^d,
$$

if, for every $\varphi \in C_0^\infty(\mathbb{R}^d)$ and $t \in [0, T]$, we have

$$
(2.31) \quad \int_{\mathbb{R}^d} \varphi(x)m(t,x)dx = \int_{\mathbb{R}^d} \varphi(x)m_0^*(x)dx - \int_0^t \int_{\mathbb{R}^d} \langle D_p H(x, D_x v[\mu](s,x)), D\varphi(x) \rangle m(s,x)dxds.
$$

**Proposition 2.2.** Assume (H1)-(H3) and let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. Then (2.31) admits a solution $m \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d)) \cap L^p([0, T] \times \mathbb{R}^d)$ and there exists $\tilde{C} > 0$ such that

$$
(2.32) \qquad \|m(t, \cdot)\|_{L^p(\mathbb{R}^d)} \le \tilde{C} \|m_0^*\|_{L^p(\mathbb{R}^d)} \quad \text{for all } t \in [0, T].
$$

If, in addition, for every $t \in [0, T]$, the functions $F(\cdot, \mu(t))$ and $G(\cdot, \mu(T))$ are differentiable, then the
solution $m$ to (2.30) is unique.

Proof. The first assertion in the statement follows from Proposition 4.6 below, while the second one
follows by arguing as in the proof of [19, Lemma 1.10]. The crucial steps in the latter are the use of the
superposition principle in [8] for solutions to (2.30) and the fact that, under the differentiability assump-
tions over $F(\cdot, \mu(t))$ and $G(\cdot, \mu(T))$, the optimal control problem inf $\{J^{0,x}[\mu](\alpha) | \alpha \in L^2([0,T]; \mathbb{R}^d)\}$
admits a unique solution for almost every $x \in \mathbb{R}^d$. $\square$

Finally, we say that $(v^*, m^*)$, with $m^* \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d)) \cap L^p([0, T] \times \mathbb{R}^d)$, solves (MFG) if $v^* = v[m^*]$ and $m^*$ solves (2.30) with $\mu = m^*$.

**Proposition 2.3.** Assume (H1)-(H3). Then system (MFG) admits a solution $(v^*, m^*)$.

*Proof.* A proof of this result, under slightly different assumptions, can be found, for instance, in [19, Section 1.3.4]. In the present context, the result follows from Theorem 5.2 below. $\square$
A uniqueness result for solutions to (MFG) can be shown under additional assumptions on the coupling terms F and G. A sufficient condition is the so-called Lasry-Lions monotonicity condition which states that, for $h = F, G$, it holds

$$ (2.33) \qquad \int_{\mathbb{R}^d} (h(x, m_1) - h(x, m_2)) d(m_1 - m_2)(x) \geq 0 \quad \text{for all } m_1, m_2 \in \mathcal{P}_1(\mathbb{R}^d). $$

**Proposition 2.4.** Assume (H1)-(H3), the monotonicity condition (2.33) and that, for all $\nu \in \mathcal{P}_1(\mathbb{R}^d)$, the functions $F(\cdot, \nu)$ and $G(\cdot, \nu)$ are differentiable. Then system (MFG) admits a unique solution.

*Proof.* The existence of a solution to (MFG) follows from Proposition 2.3 while, under (2.33) and the differentiability assumptions on $F(\cdot, \nu)$ and $G(\cdot, \nu)$, the proof of the uniqueness of the solution follows by arguing as in the proof of [19, Theorem 1.8]. □

### 3. A SEMI-LAGRANGIAN SCHEME FOR THE HJB EQUATION

Let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. In this section we recall a standard semi-Lagrangian scheme to approximate the viscosity solution $v[\mu]$ to (2.21). Most of the results for the semi-Lagrangian scheme that will be needed in the remainder of the article follow similarly to those in the monograph [29] and the contributions [20, 21, 27]. The principal differences come from our assumptions on L in (H1), which allow us to consider cost functionals not covered in these references (see e.g. the example in the last paragraph of Remark 2.1). Therefore, we confine ourselves to explain the main changes in the proofs of the aforementioned properties. On the other hand, the estimate in Proposition 3.5 below seems to be new and will play a key role later in this article.

In order to define the scheme, let $N \in \mathbb{N}^*$ be the number of time steps, let $\Delta t = T/N$ be the time step, let $\mathcal{I}_{\Delta t} = \{0, \dots, N\}$, let $\mathcal{I}_{\Delta t}^* = \mathcal{I}_{\Delta t} \setminus \{N\}$, let $t_k = k\Delta t$ for all $k \in \mathcal{I}_{\Delta t}$, and set $\mathcal{G}_{\Delta t} = \{t_k | k \in \mathcal{I}_{\Delta t}\}$. Given a space step $\Delta x > 0$ and $i = (i_1, \dots, i_d) \in \mathbb{Z}^d$, define $\beta_i^1: \mathbb{R}^d \to \mathbb{R}$ as

$$ (3.1) \qquad \beta_i^1(z) = \prod_{l=1}^{d} \hat{\beta}\left(\frac{z_l}{\Delta x} - i_l\right) \quad \text{for all } z = (z_1, \dots, z_d) \in \mathbb{R}^d, $$

where

$$ (3.2) \qquad \hat{\beta}(\xi) = \max\{0, 1 - |\xi|\} \quad \text{for all } \xi \in \mathbb{R}. $$

Notice that $\beta_i^1 \ge 0$, $\sum_{i \in \mathbb{Z}^d} \beta_i^1(x) = 1$ for all $x \in \mathbb{R}^d$ and, setting $x_i = i\Delta x$, we have $\beta_i^1(x_j) = 1$, if $i = j$, and $\beta_i^1(x_j) = 0$, otherwise. Let $\mathcal{G}_{\Delta x} = \{i\Delta x | i \in \mathbb{Z}^d\}$ be the uniform grid and, given $\phi: \mathcal{G}_{\Delta x} \to \mathbb{R}$, define its interpolate as

$$ I^1[\phi](x) = \sum_{i \in \mathbb{Z}^d} \beta_i^1(x) \phi_i \quad \text{for all } x \in \mathbb{R}^d, $$

where, for notational simplicity, we have set $\phi_i = \phi(x_i)$. For every $\varphi: \mathbb{R}^d \to \mathbb{R}$ denote by $\varphi|_{\mathcal{G}_{\Delta x}}$ its restriction to $\mathcal{G}_{\Delta x}$. If $\varphi$ is of class $C^2$ and has bounded second order derivatives, it follows from [46, Remark 3.4.2] that

$$ (3.3) \qquad \|\varphi(x) - I[\varphi|_{\mathcal{G}_{\Delta x}}](x)\|_{\infty} \le C_{\varphi}(\Delta x)^2, $$

where $C_{\varphi} > 0$ depends only on $\varphi$.

We consider the following fully-discrete semi-Lagrangian scheme: find $\{v_k: \mathcal{G}_{\Delta x} \to \mathbb{R} | k \in \mathcal{I}_{\Delta t}\}$ such that

$$ (3.4) \qquad \begin{aligned} v_{k,i} &= S_{k,i}^{\text{fd}}[\mu](v_{k+1}) && \text{for all } k \in \mathcal{I}_{\Delta t}^*, i \in \mathbb{Z}^d, \\ v_{N,i} &= G(x_i, \mu(T)) && \text{for all } i \in \mathbb{Z}^d, \end{aligned} $$

where, for every $\phi: \mathcal{G}_{\Delta x} \to \mathbb{R}$, bounded, $k \in \mathcal{I}_{\Delta t}$, and $i \in \mathbb{Z}^d$,

$$ (3.5) \qquad S_{k,i}^{\text{fd}}[\mu](\phi) = \inf_{a \in \bar{B}(0, C_b)} [\Delta t L(x_i, a) + I^1[\phi](x_i - \Delta t a)] + \Delta t F(x_i, \mu(t_k)). $$

Notice that, being explicit, scheme (3.4) admits a unique solution. By definition, $S^{\text{fd}}[\mu]$ is monotone, i.e. for every $\phi^1, \phi^2: \mathcal{G}_{\Delta x} \to \mathbb{R}$, bounded, with $\phi_i^1 \le \phi_i^2$ for all $i \in \mathbb{Z}^d$, we have that

$$ (3.6) \qquad S_{k,i}^{\text{fd}}[\mu](\phi^1) \le S_{k,i}^{\text{fd}}[\mu](\phi^2) \quad \text{for all } k \in \mathcal{I}_{\Delta t}^*, i \in \mathbb{Z}^d. $$
Moreover, using (H1) and (H2), standard arguments (see e.g. [29, Section 5.2.3]) yield the following consistency property for $S^{\text{fd}}[\mu]$ : let $(\mu_n)_{n \in \mathbb{N}} \subset C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, $((\Delta t_n, \Delta x_n))_{n \in \mathbb{N}} \subset [0, \infty[$, $((t_{k_n}, x_{i_n}))_{n \in \mathbb{N}} \subset \mathcal{G}_{\Delta t_n} \times \mathcal{G}_{\Delta x_n}$, and $(t, x) \in ]0, T[ \times \mathbb{R}^d$ such that, as $n \to \infty$, $\mu_n \to \mu$, $(\Delta t_n, \Delta x_n) \to 0$, $(\Delta x_n)^2/\Delta t_n \to 0$, and $(t_{k_n}, x_{i_n}) \to (t, x)$. Then, recalling the definition of $H_b$ in (2.29), for every $\varphi: [0, T] \times \mathbb{R}^d \to \mathbb{R}$ of class $C^1$, with bounded derivatives, we have

$$ (3.7) \quad \lim_{n \to \infty} \frac{1}{\Delta t_n} \left( \varphi(x_{i_n}, t_{k_n}) - S_{k_n, i_n}^{\text{fd}}[\mu_n] (\varphi(t_{k_n+1}, \cdot)|_{\mathcal{G}_{\Delta x_n}}) \right) = -\partial_t \varphi(x, t) + H_b(x, D_x \varphi(t, x)) - F(x, \mu(t)). $$
Given $(\Delta t, \Delta x) \in ]0, \infty[$, let us set

$$ (3.8) \quad v^{\Delta t, \Delta x}[\mu](t_k, x) = I^1[v_k](x) \quad \text{for all } k \in \mathcal{I}_{\Delta t}, x \in \mathbb{R}^d, $$

where, for every $k \in \mathcal{I}_{\Delta t}$, $v_k: \mathcal{G}_{\Delta x} \to \mathbb{R}$ is computed with (3.5). We extend this definition to $[0, T] \times \mathbb{R}^d$, by setting

$$ (3.9) \quad v^{\Delta t, \Delta x}[\mu](t, x) = v^{\Delta t, \Delta x}[\mu](t_k, x) \quad \text{if } t \in [t_k, t_{k+1}, [k, k \in \mathcal{I}_{\Delta t}^*]. $$

The following result provides properties for $v^{\Delta t, \Delta x}[\mu]$ that are analogous to those in Proposition 2.1(ii)-(iv) for $v[\mu]$.

**Proposition 3.1.** Assume (H1)-(H2), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $(\Delta t, \Delta x) \in ]0, \infty[$. Then the following hold:

(i) [Stability] We have

$$ (3.10) \quad |v^{\Delta x, \Delta t}[\mu](t, x)| \le \tilde{C}_v \quad \text{for all } (t, x) \in [0, T] \times \mathbb{R}^d, $$
where $\tilde{C}_v > 0$ is independent of $(\mu, \Delta t, \Delta x)$.

(ii) [Lipschitz property] We have

$$ (3.11) \quad |v^{\Delta t, \Delta x}[\mu](t, x) - v^{\Delta t, \Delta x}[\mu](t, y)| \le \tilde{C}_{\text{Lip}} |x - y| \quad \text{for all } t \in [0, T], x, y \in \mathbb{R}^d, $$
where $\tilde{C}_{\text{Lip}} > 0$ is independent of $(\mu, \Delta t, \Delta x)$.

(iii) [Discrete semi-concavity] We have

$$ (3.12) \quad v^{\Delta t, \Delta x}[\mu](t, x + x_i) - 2v^{\Delta t, \Delta x}[\mu](t, x) + v^{\Delta t, \Delta x}[\mu](t, x - x_i) \\ \le \tilde{C}_{\text{sc}} |x_i|^2 \quad \text{for all } t \in [0, T], x \in \mathbb{R}^d, i \in \mathbb{Z}^d, $$
where $\tilde{C}_{\text{sc}} > 0$ is independent of $(\mu, \Delta t, \Delta x)$.

*Proof.* (i): This follows directly from (3.4), (2.6), (2.16), (2.10), (2.11), and iteration.

(ii): It follows from (2.7) that

$$ (3.13) \quad |L(x, a) - L(y, a)| \le C_{L,3} (1 + C_b^2) |x - y| \quad \text{for all } x, y \in \mathbb{R}^d, a \in \bar{B}(0, C_b). $$

Using this inequality, (2.12), and (2.13), the result follows from the same arguments than those in [21, Lemma 3.1(i)] (see also the proof of [27, Lemma 5.3(a)]).

(iii): It follows from (2.9) that

$$ (3.14) \quad L(x + y, a) - 2L(x, a) + L(x - y, a) \le C_{L,5} (1 + C_b^2) |y|^2 \quad \text{for all } x, y \in \mathbb{R}^d, a \in \bar{B}(0, C_b). $$

In turn, using (2.14) and (2.15), the result follows by arguing as in the proof of [2, Lemma 4.1] (see also the proof of [27, Lemma 3.2(ii)]). $\square$

Using the monotonicity of $S^{\text{fd}}[\mu]$, the consistency property in (3.7), and the stability result in Proposition 3.1(i), the Barles-Souganidis relaxed limit method (see [13]) yields the following convergence result (see [20, Theorem 3.3] for a detailed proof).

**Proposition 3.2.** Assume (H1)-(H2), let $(\mu_n)_{n \in \mathbb{N}} \subset C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $((\Delta t_n, \Delta x_n))_{n \in \mathbb{N}} \subset [0, \infty[$. Suppose that, as $n \to \infty$, $\mu_n \to \mu$, for some $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, $(\Delta t_n, \Delta x_n) \to 0$, and $(\Delta x_n)^2/\Delta t_n \to 0$. Then $(v^{\Delta t_n, \Delta x_n}[\mu_n])_{n \in \mathbb{N}}$ converges to $v[\mu]$ uniformly over compact subsets of $[0, T] \times \mathbb{R}^d$.
Given $\varepsilon > 0$, consider the mollifier $\mathbb{R}^d \ni x \mapsto \rho_\varepsilon(x) = \rho(x/\varepsilon)/\varepsilon^d \in \mathbb{R}^d$, where $\rho \in C^\infty(\mathbb{R}^d)$ has bounded derivatives of any order and satisfies $\rho(\mathbb{R}^d) \subset [0, \infty[$ and $\int_{\mathbb{R}^d} \rho(x)dx = 1$. Given $\varphi \in W^{1,\infty}(\mathbb{R}^d)$, a standard computation shows that

$$
(3.15) \qquad \sup_{x \in \mathbb{R}^d} |(\rho_\varepsilon * \varphi)(x) - \varphi(x)| \le \varepsilon \|D\varphi\|_{L^\infty(\mathbb{R}^d)}
$$

$$
(3.16) \qquad \sup_{x \in \mathbb{R}^d} \| D^\ell (\rho_\varepsilon * \varphi)(x) \| \le c_\ell \varepsilon^{1-\ell} \quad \text{for all } \ell \in \mathbb{N},
$$

where $\|D^\ell(\rho_\varepsilon * \varphi)(x)\|$ denotes the operator norm of $D^\ell(\rho_\varepsilon * \varphi)(x)$ and $c_\ell > 0$ depends only on $\ell$. Let us set $\Delta = (\Delta t, \Delta x, \varepsilon)$ and define

$$
(3.17) \qquad v^{\Delta}[\mu](t, \cdot) = \rho_{\varepsilon} * v^{\Delta t, \Delta x}[\mu](t, \cdot) \quad \text{for all } t \in [0, T].
$$

The function $v^{\Delta}[\mu]$ satisfies similar properties than $v^{\Delta t, \Delta x}[\mu]$, as the following proposition shows.

**Proposition 3.3.** Assume (H1)-(H2), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[^3$.
Then the following holds:

(i) [Stability] We have

$$
(3.18) \qquad |v^\Delta[\mu](t, x)| \le \tilde{C}_v \quad \text{for all } (t, x) \in [0, T] \times \mathbb{R}^d,
$$
with $\tilde{C}_v > 0$ being as in Proposition 3.1(i).

(ii) [Lipschitz property] We have

$$
(3.19) \qquad |v^\Delta[\mu](t, x) - v^\Delta[\mu](t, y)| \le \tilde{C}_{\text{Lip}} |x - y| \quad \text{for all } t \in [0, T], x, y \in \mathbb{R}^d,
$$
with $\tilde{C}_{\text{Lip}} > 0$ being as in Proposition 3.1(ii).

(iii) [Approximate semi-concavity] We have

$$
(3.20) \quad v^{\Delta}[\mu](t, x+y) - 2v^{\Delta}[\mu](t, x) + v^{\Delta}[\mu](t, x-y) \\
\le \tilde{C}_{\text{asc}} \left( |y|^2 + (\Delta x)^2 + \frac{(\Delta x)^2}{\varepsilon} \right) \quad \text{for all } t \in [0, T], x, y \in \mathbb{R}^d,
$$
where $\tilde{C}_{\text{asc}} > 0$ is independent of $(\mu, \Delta)$.

*Proof.* Assertions (i) and (ii) follow directly from (3.17) and the corresponding assertions in Proposition 3.1. The proof of (iii) follows from Proposition 3.1(iii) and arguing as in the proof of [2, Lemma 4.2] (see also the proof of [27, Lemma 5.5(b)]).
□

In the following, given $A \subset \mathbb{R}^d$, we denote by $I_A$ the indicator function of $A$. The convergence result in the following proposition will play an important role in the next section.

**Proposition 3.4.** Assume (H1)-(H2), let $(\mu_n)_{n \in \mathbb{N}} \subset C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $((\Delta t_n, \Delta x_n, \varepsilon_n))_{n \in \mathbb{N}} \subset ]0, \infty[^3$. Set $\Delta_n = (\Delta t_n, \Delta x_n, \varepsilon_n)$ and let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. Suppose that, as $n \to \infty$, $\mu_n \to \mu$, $\Delta_n \to 0$, and $(\Delta x_n)^2/\Delta t_n \to 0$. Then the following hold:

(i) $(v^{\Delta_n}[\mu_n])_{n \in \mathbb{N}}$ converges to $v[\mu]$ uniformly over compact subsets of $[0, T] \times \mathbb{R}^d$.

(ii) If, in addition, $\Delta x_n/\varepsilon_n \to 0$, then, for every $K \subset [0, T] \times \mathbb{R}^d$ compact and $q \in [1, \infty[$,

$$
(3.21) \qquad \mathbb{I}_K D_p H(\cdot, D_x v^{\Delta_n}[\mu_n]) \to \mathbb{I}_K D_p H(\cdot, D_x v[\mu]) \quad \text{in } L^q([0, T] \times \mathbb{R}^d).
$$

*Proof.* (i): This follows from Proposition 3.3(ii), estimate (3.15), and Proposition 3.2.

(ii): It follows from Proposition 3.3(iii) and [2, Remark 6] that

$$
(3.22) \quad \left\langle D_x v^{\Delta}[\mu](t, y) - D_x v^{\Delta}[\mu](t, x), y - x \right\rangle \le c \left( |y - x|^2 + \frac{(\Delta x)^2}{\varepsilon^2} \right) \quad \text{for all } t \in [0, T], x, y \in \mathbb{R}^d.
$$

Using this inequality and arguing as in the proof of [20, Theorem 3.5], one deduces that, as $n \to \infty$, $D_x v^{\Delta_n}[\mu_n] \to D_x v[\mu]$ almost everywhere in $[0, T] \times \mathbb{R}^d$. In turn, since $H$ is of class $C^2$, we get that $D_p H(\cdot, D_x v^{\Delta_n}[\mu_n]) \to D_p H(\cdot, D_x v[\mu])$ almost everywhere in $[0, T] \times \mathbb{R}^d$. Thus, (3.21) follows from Proposition 3.3(ii) and Lebesgue's dominated convergence theorem.
$\square$

We conclude this section with a useful estimate for $D^2 v^\Delta[\mu]$.
**Proposition 3.5.** Assume (H1)-(H2), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[^{2\times} ]0, 1[$. Then it holds that

$$ (3.23) \qquad \langle D_x^2 v^\Delta[\mu](t, x)y, y \rangle \le \tilde{C}_{hb} \left( 1 + \frac{(\Delta x)^2}{\varepsilon^4} \right) |y|^2 \quad \text{for all } t \in [0, T], x, y \in \mathbb{R}^d, $$

where $\tilde{C}_{hb} > 0$ is independent of $(\mu, \Delta)$.

*Proof.* Let us fix $t \in [0, T]$ and $x, y \in \mathbb{R}^d$. If $y=0$ the result is true and, hence, let us assume that $y \ne 0$ and set $\tau = \Delta x/(\sqrt{\varepsilon}|y|)$. In what follows, $C > 0$ denotes a constant, independent of $(\mu, \Delta)$ which may change from line to line. It follows from Proposition 3.3(iii) that

$$ (3.24) \qquad v^\Delta[\mu](t, x + \tau y) - 2v^\Delta[\mu](t, x) + v^\Delta[\mu](t, x - \tau y) \le C\tau^2|y|^2. $$

On the other hand, a Taylor expansion of order 4 and (3.16) imply that

$$ v^\Delta[\mu](t, x + \tau y) \ge v^\Delta[\mu](t, x) + \langle D_x v^\Delta[\mu](t, x), \tau y \rangle + \frac{1}{2} \langle D_x^2 v^\Delta[\mu](t, x) \tau y, \tau y \rangle \\ + \frac{1}{6} D_x^3 v^\Delta[\mu](t, x) (\tau y)^3 - \frac{1}{\varepsilon^3} C |\tau y|^4, $$

$$ v^\Delta[\mu](t, x - \tau y) \ge v^\Delta[\mu](t, x) - \langle D_x v^\Delta[\mu](t, x), \tau y \rangle + \frac{1}{2} \langle D_x^2 v^\Delta[\mu](t, x) \tau y, \tau y \rangle \\ - \frac{1}{6} D_x^3 v^\Delta[\mu](t, x) (\tau y)^3 - \frac{1}{\varepsilon^3} C |\tau y|^4, $$

Adding both inequalities, using (3.24) and the relation $\tau|y| = \Delta x/\sqrt{\varepsilon}$, we get

$$ (3.25) \qquad \langle D_x^2 v^\Delta[\mu](t, x)\tau y, \tau y \rangle \le C \left( 1 + \frac{(\Delta x)^2}{\varepsilon^4} \right) \tau^2 |y|^2. $$

Dividing by $\tau^2$ yields (3.23). $\square$

### 4. A LAGRANGE-GALERKIN TYPE SCHEME FOR THE CONTINUITY EQUATION

Given $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ and $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[^{3}$, let $v^\Delta[\mu]$ be defined as in (3.17). Consider the continuity equation

$$ (4.1) \qquad \begin{aligned} \partial_t m - \operatorname{div} (D_p H(x, D_x v^\Delta[\mu]) m) &= 0 && \text{in } ]0, T[ \times \mathbb{R}^d, \\ m(0) &= m_0^* && \text{in } \mathbb{R}^d. \end{aligned} $$

Since $H$ is of class $C^2$ and $D_x v^\Delta[\mu]$ is bounded and Lipschitz, by [9, Proposition 8.1.8] equation (4.1) admits a unique solution $m^\Delta[\mu] \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, which can be represented as

$$ (4.2) \qquad m^\Delta[\mu](t) = \Phi^\Delta[\mu](0, t, \cdot)\sharp m_0^* \quad \text{for all } t \in [0, T], $$

where, for all $s \in [0, T)$ and $x \in \mathbb{R}^d$, $\Phi^\Delta[\mu](s, \cdot, x)$ denotes the unique solution to

$$ (4.3) \qquad \begin{aligned} \dot{X}(t) &= -D_p H(X(t), D_x v^\Delta[\mu](t, X(t))) && \text{for a.e. } t \in (s, T), \\ X(s) &= x. \end{aligned} $$

Relations (4.2) and (4.3) imply that

$$ (4.4) \qquad m^\Delta[\mu](t) = \Phi^\Delta[\mu](s, t, \cdot)\sharp m^\Delta[\mu](s) \quad \text{for all } s, t \in [0, T], s \le t. $$

On the other hand, notice that (2.20) and the uniform bound $\|D_x v^\Delta[\mu](t, \cdot)\|_{L^\infty(\mathbb{R}^d)} \le \tilde{C}_{Lip}$ for all $t \in [0, T]$, which follows from Proposition 3.3(ii), yield the existence of $C_{bf} > 0$, independent of $(\mu, \Delta)$, such that

$$ (4.5) \qquad \sup_{(t,x) \in [0,T] \times \mathbb{R}^d} |D_p H(x, D_x v^\Delta[\mu](t, x))| \le C_{bf}. $$

Relation (4.4) and estimate (4.5) have two straightforward consequences. The first one is that, by (H3)(i), we have

$$ (4.6) \qquad \operatorname{supp} (m^\Delta[\mu](t)) \subset \bar{B}(0, C^* + TC_{bf}) \quad \text{for all } t \in [0, T]. $$
The second one, is the uniform equicontinuity of the family $\{m^\Delta[\mu] \mid \Delta \in ]0, \infty[\}^n$ in $C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. More precisely, using (2.2), an easy computation shows that

$$ (4.7) \qquad d_1(m^\Delta[\mu](t), m^\Delta[\mu](s)) \le C_{bf}|t-s| \quad \text{for all } s, t \in [0, T]. $$

The purpose of the following two propositions is to provide some stability estimates for $m^\Delta[\mu]$, which, together with (4.6) and (4.7), motivate the forthcoming analysis for a LG discretization of (4.1).

**Proposition 4.1.** Assume (H1)-(H3), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[$. Then for every $c > 0$ there exists $\tilde{c} > 0$, independent of $\mu$, such that, if $\Delta x \le c\varepsilon^2$, it holds that

$$ (4.8) \qquad \|m^\Delta[\mu](t, \cdot)\|_{L^p(\mathbb{R}^d)} \le \tilde{c} \|m_0^*\|_{L^p(\mathbb{R}^d)} \quad \text{for all } t \in [0, T]. $$

*Proof.* Let $c > 0$ and suppose that $\Delta x \le c\varepsilon^2$. Proposition 3.5 implies that

$$ (4.9) \qquad \langle D_x^2 v^\Delta[\mu](t, x)y, y \rangle \le \tilde{C}_{hb} (1 + c^2) |y|^2 \quad \text{for all } t \in [0, T], x, y \in \mathbb{R}^d. $$

Using this inequality, the proof of (4.8) follows from exactly the same arguments than those in the proof of [30, Proposition 4.1]. □

**Proposition 4.2.** Assume (H1)-(H3), let $(\mu_n)_{n \in \mathbb{N}} \subset C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $((\Delta t_n, \Delta x_n, \varepsilon_n))_{n \in \mathbb{N}} \subset (0, \infty)^3$. Set $\Delta_n = (\Delta t_n, \Delta x_n, \varepsilon_n)$ and let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. Suppose that, as $n \to \infty$, $\mu_n \to \mu$, $\Delta_n \to 0$, $(\Delta x_n)^2/\Delta t_n \to 0$, and $\Delta x_n = O(\varepsilon_n^2)$. Then, up to some subsequence, the following hold:

(i) $(v^{\Delta_n}[\mu^n])_{n \in \mathbb{N}}$ converges to $v[\mu]$, uniformly over compact subsets of $[0, T] \times \mathbb{R}^d$, and, for every $K \subset [0, T] \times \mathbb{R}^d$ compact and $q \in [1, \infty[$, $(\mathbb{I}_K D_x v^{\Delta_n}[\mu_n])_{n \in \mathbb{N}}$ converges to $\mathbb{I}_K D_x v[\mu]$ in $L^q([0, T] \times \mathbb{R}^d)$.

(ii) $(m^{\Delta_n}[\mu_n])_{n \in \mathbb{N}}$ converges in $C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ towards a solution to (2.30). Moreover, the convergence also hold weakly in $L^p([0, T] \times \mathbb{R}^d)$, if $p < \infty$, and weakly* in $L^\infty([0, T] \times \mathbb{R}^d)$, if $p = \infty$.

*Proof.* For every $n \in \mathbb{N}$, let us set $v^n = v^{\Delta_n}[\mu_n]$ and $m^n = m^{\Delta_n}[\mu_n]$.

(i): This follows from Proposition 3.4.

(ii): It follows from (4.6), [9, Proposition 7.1.5], (4.7), and the Ascoli-Arzelà theorem, that there exists $m^* \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ such that, as $n \to \infty$ and up to some subsequence, $m^n \to m^*$ in $C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. By Proposition 4.1, the convergence also hold weakly in $L^p([0, T] \times \mathbb{R}^d)$, if $p < \infty$, and weakly* in $L^\infty([0, T] \times \mathbb{R}^d)$, if $p = \infty$. Since $m^n$ solves (4.1), for every $t \in [0, T]$, and $\varphi \in C_0^\infty(\mathbb{R}^d)$, we have

$$ (4.10) \qquad \int_{\mathbb{R}^d} \varphi(x)m^n(t, x)dx = \int_{\mathbb{R}^d} \varphi(x)m_0^*(t, x)dx \\ \qquad \qquad + \int_0^t \int_{\mathbb{R}^d} \langle D\varphi(x), D_pH(x, D_xv^n(s, x)) \rangle m^n(s, x)dxd\bar{s}. $$

Thus, by (i), we can pass to the limit in the previous expression to deduce that $m^*$ solves (2.30). □

**4.1. The Lagrange-Galerkin approximation.** The main purpose of this section is to provide some results in the vein of Propositions 4.1 and 4.2 for solutions $m^\Delta[\mu]$ to a LG approximation of (4.1) that we proceed to construct.

Let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ and let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[$. For every $k \in \mathcal{I}_{\Delta t}^*$ and $x \in \mathbb{R}^d$, let $\Phi_k^\Delta[\mu](x)$ be the explicit Euler approximation of $\Phi^\Delta[\mu](t_k, t_{k+1}, x)$, i.e.

$$ (4.11) \qquad \Phi_k^\Delta[\mu](x) = x - \Delta t D_p H(x, D_x v^\Delta[\mu](t_k, x)) \quad \text{for all } x \in \mathbb{R}^d. $$

As in [17], we consider the following semi-discrete approximation of (4.4):

$$ (4.12) \qquad \begin{aligned} m_{k+1} &= \Phi_k^\Delta[\mu] \sharp m_k \quad \text{for all } k \in \mathcal{I}_{\Delta t}^*, \\ m_0 &= m_0^*, \end{aligned} $$

or, equivalently, for every $k \in \mathcal{I}_{\Delta t}^*$ and every Borel function $\varphi: \mathbb{R}^d \to \mathbb{R}$, such that $\varphi(\Phi_k^\Delta[\mu](\cdot))$ is integrable with respect to $m_k$,

$$ (4.13) \qquad \int_{\mathbb{R}^d} \varphi(x) dm_{k+1}(x) = \int_{\mathbb{R}^d} \varphi(\Phi_k^\Delta[\mu](x)) dm_k(x). $$
Following [43], which mainly deals with a LG approximation of the dual (or transport) equation associated to (4.1), let us formally deduce from (4.13) a time-space approximation of (4.1). For every $i \in \mathbb{Z}^d$, set

$$ (4.14) \qquad E_i = \{x \in \mathbb{R}^d \mid |x - x_i|_{\infty} \le \Delta x/2\} $$

and define

$$ (4.15) \qquad m_{0,i} = \frac{1}{(\Delta x)^d} \int_{E_i} m_0^*(x)dx. $$

Given the regular mesh defined by $\{E_i | i \in \mathbb{Z}^d\}$, let $\{\beta_i | i \in \mathbb{Z}^d\}$ be a finite element basis. In the following, we look for an approximation $m^\Delta[\mu]$ of the solution $m^\Delta[\mu]$ to (4.1) such that

$$ (4.16) \qquad m^\Delta[\mu](t_k, x) = \sum_{j \in \mathbb{Z}^d} m_{k,j} \beta_j(x) \quad \text{for all } k \in \mathcal{I}_{\Delta t}, x \in \mathbb{R}^d, $$

for some constants $\{m_{k,j} | k \in \mathcal{I}_{\Delta t}, j \in \mathbb{Z}^d\}$. In order to determine the latter, we replace $m_k$ and $m_{k+1}$ in (4.13) by $m^\Delta(t_k, \cdot)$ and $m^\Delta(t_{k+1}, \cdot)$, respectively, and, given $i \in \mathbb{Z}^d$, we take $\varphi = \beta_i$ to obtain the following equations

$$ (4.17) \qquad \sum_{j \in \mathbb{Z}^d} m_{k+1,j} \int_{\mathbb{R}^d} \beta_j(x) \beta_i(x) dx = \sum_{j \in \mathbb{Z}^d} m_{k,j} \int_{\mathbb{R}^d} \beta_i(\Phi_k^\Delta[\mu](x)) \beta_j(x) dx \quad \text{for all } k \in \mathcal{I}_{\Delta t}^*. $$

In the context of second-order Fokker-Planck equations, scheme (4.17) has already been considered in [16] to provide a high-order accurate LG scheme to solve second-order mean field games problems with smooth solutions. In this reference, the authors consider symmetric Lagrangian basis of odd order which preserve the mass but not the positivity of the initial condition $\{m_{0,i} | i \in \mathbb{Z}^d\}$. Since we aim to approximate solutions to (4.1), which in general are not smooth, from now on we take $\beta_i = \beta_i^0 := \mathbb{I}_{E_i}$ for all $i \in \mathbb{Z}^d$. Under this choice, (4.17) and (4.15) yield the following LG scheme for (4.1):

$$ (4.18) \qquad m_{k+1,i} = \frac{1}{(\Delta x)^d} \sum_{j \in \mathbb{Z}^d} m_{k,j} \int_{E_j} \beta_i^0(\Phi_k^\Delta[\mu](x)) dx \quad \text{for all } k \in \mathcal{I}_{\Delta t}^*, i \in \mathbb{Z}^d, $$

$$ (4.19) \qquad m_{0,i} = \frac{1}{(\Delta x)^d} \int_{E_i} m_0^*(x)dx \quad \text{for all } i \in \mathbb{Z}^d. $$

The scheme above is explicit and hence admits a unique solution. Interestingly, the scheme (4.18)-(4.19) coincides with the one proposed in [45] (see also [48]) to approximate solutions to continuity equations. Indeed, we have

$$ (4.20) \qquad \int_{E_j} \beta_i^0(\Phi_k^\Delta[\mu](x))dx = \int_{\mathbb{R}^d} \mathbb{I}_{E_j \cap \Phi_k^\Delta[\mu]^{-1}(E_i)}(x)dx = \mathcal{L}^d(E_j \cap \Phi_k^\Delta[\mu]^{-1}(E_i)), $$

where $\mathcal{L}^d$ denotes the Lebesgue measure in $\mathbb{R}^d$. Plugging this expression in (4.18) yields the scheme in [45, Section 2.2]. Notice that our main results for solutions to (4.18)-(4.19), contained in Propositions 4.5 and 4.6 below, do not follow from the results in [45, 48]. Therefore, the analysis in this section provides a complementary study to the one in [45, 48] for the approximation (4.18)-(4.19) of continuity equations.

## 4.2. Properties of LG scheme.
We begin with a preliminary result stating that the solution to (4.18)-(4.19) is supported on a compact set, which is independent of the discretization parameters provided that $\Delta x$ is of the order of $\Delta t$.

**Lemma 4.1.** Assume (H1)-(H3), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[$, and let $\{m_{k,i} | k \in \mathcal{I}_{\Delta t}, i \in \mathbb{Z}^d\}$ be the solution to (4.18)-(4.19). Then for every $c > 0$ there exists $\tilde{C}^* > 0$, independent of $\mu$, such that, if $\Delta x \le c\Delta t$, for every $k \in \mathcal{I}_{\Delta t}$ we have $m_{k,i} = 0$ if $x_i \notin \bar{B}_\infty(0, \tilde{C}^*)$.

*Proof.* Let $c > 0$ and suppose that $\Delta x \le c\Delta t$. For every $k \in \mathcal{I}_{\Delta t}^*$, set $r_k = \sup\{|x_i|_\infty | m_{k,i} \ne 0, i \in \mathbb{Z}^d\} \in [0, \infty]$. By (4.5), (4.18), and (H3)(i), we have

$$ (4.21) \qquad r_{k+1} \le r_k + \Delta t C_{\text{bf}} + \frac{\Delta x}{2} \le r_k + \Delta t \left( C_{\text{bf}} + \frac{c}{2} \right) \le C^* + \Delta x N_{\Delta t} \Delta t \left( C_{\text{bf}} + \frac{c}{2} \right) = C^* + T \left( C_{\text{bf}} + \frac{c}{2} \right), $$

for all $k \in \mathcal{I}_{\Delta t}^*$. The result follows by letting $\tilde{C}^* = C^* + T(C_{\text{bf}} + c/2)$. $\square$
Let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ and let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[^3$. As a consequence of the previous result, in (4.18) it suffices to compute $m_{i,k+1}$ for $i \in \mathcal{I}_{\Delta x}$, where

$$ (4.22) \qquad \mathcal{I}_{\Delta x} := \{i \in \mathbb{Z}^d | x_i \in \bar{B}_{\infty}(0, \tilde{C}^*)\}. $$

Given the constants $\{m_{k,i} | k \in \mathcal{I}_{\Delta t}, i \in \mathbb{Z}^d\}$, computed with (4.18)-(4.19), we extend $\mathbf{m}^{\Delta}[\mu]$, given by (4.16), to $[0, T] \times \mathbb{R}^d$ as follows:

$$ (4.23) \qquad \mathbf{m}^{\Delta}[\mu](t, x) = \left(\frac{t_{k+1} - t}{\Delta t}\right) \mathbf{m}^{\Delta}[\mu](t_k, x) + \left(\frac{t - t_k}{\Delta t}\right) \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) \\ \qquad \qquad \qquad \qquad \text{for all } k \in \mathcal{I}_{\Delta t}^*, t \in [t_k, t_{k+1}), x \in \mathbb{R}^d. $$

In the following proposition, we state, for later use, some simple properties of the solution to (4.18)-(4.19).

**Proposition 4.3.** Assume (H1)-(H3), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[^3$, and let $\{m_{k,i} | k \in \mathcal{I}_{\Delta t}, i \in \mathbb{Z}^d\}$ be the solution to (4.18)-(4.19). Then the following hold:

(i) $\mathbf{m}^{\Delta}[\mu](t, x) \ge 0$ for all $t \in [0, T]$, $x \in \mathbb{R}^d$.

(ii) Let $c > 0$. If $\Delta x \le c\Delta t$ and $\tilde{C}^* > 0$ is given by Lemma 4.1, we have

$$ (4.24) \qquad \operatorname{supp} (\mathbf{m}^{\Delta}[\mu](t, \cdot)) \subseteq \bar{B}_{\infty}(0, \tilde{C}^*) \quad \text{for all } t \in [0, T]. $$

(iii) Let $a = (a_i)_{i \in \mathbb{Z}^d} \subset \mathbb{R}$ and set $\varphi_a(x) = \sum_{i \in \mathbb{Z}^d} a_i \beta_i^0(x)$ for all $x \in \mathbb{R}^d$. Then we have

$$ (4.25) \qquad \int_{\mathbb{R}^d} \varphi_a(x) \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) dx = \int_{\mathbb{R}^d} \varphi_a(\Phi_k^{\Delta}[\mu](x)) \mathbf{m}^{\Delta}[\mu](t_k, x) dx \quad \text{for all } k \in \mathcal{I}_{\Delta t}^*. $$

$$ (4.26) \qquad \int_{\mathbb{R}^d} \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) dx = \sum_{i \in \mathbb{Z}^d} \sum_{j \in \mathbb{Z}^d} a_j m_{k+1,i} \int_{\mathbb{R}^d} \beta_i^0(x) \beta_j^0(x) dx = \sum_{i \in \mathbb{Z}^d} a_i m_{k+1,i} (\Delta x)^d \\ = \sum_{i \in \mathbb{Z}^d} a_i \sum_{j \in \mathbb{Z}^d} m_{k,j} \int_{E_j} \beta_i^0(\Phi_k^{\Delta}[\mu](x)) dx = \sum_{j \in \mathbb{Z}^d} m_{k,j} \int_{\mathbb{R}^d} \varphi_a(\Phi_k^{\Delta}[\mu](x)) \beta_j^0(x) dx \\ = \int_{\mathbb{R}^d} \varphi_a(\Phi_k^{\Delta}[\mu](x)) \mathbf{m}^{\Delta}[\mu](t_k, x) dx. $$

*Proof.* (i): Using that $m_0^* \ge 0$ and $\beta_i^0 \ge 0$ for all $i \in \mathbb{Z}^d$, this assertion follows directly from (4.18)-(4.19).

(ii): This follows from Lemma 4.1 and (4.23).

(iii): For every $k \in \mathcal{I}_{\Delta t}^*$, by (4.18), we have

$$ (4.27) \qquad \int_{\mathbb{R}^d} \varphi_a \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) dx = \sum_{i \in \mathbb{Z}^d} \sum_{j \in \mathbb{Z}^d} a_j m_{k+1,i} \int_{\mathbb{R}^d} \beta_i^0(x) \beta_j^0(x) dx = \sum_{i \in \mathbb{Z}^d} a_i m_{k+1,i} (\Delta x)^d \\ = \sum_{i \in \mathbb{Z}^d} a_i \sum_{j \in \mathbb{Z}^d} m_{k,j} \int_{E_j} \beta_i^0(\Phi_k^{\Delta}[\mu](x)) dx = \sum_{j \in \mathbb{Z}^d} m_{k,j} \int_{\mathbb{R}^d} \varphi_a(\Phi_k^{\Delta}[\mu](x)) \beta_j^0(x) dx \\ = \int_{\mathbb{R}^d} \varphi_a(\Phi_k^{\Delta}[\mu](x)) \mathbf{m}^{\Delta}[\mu](t_k, x) dx. $$

Notice that the changes of the order of summation above are justified by (ii).

(iv): By (iii), with $\varphi_a(x) := \sum_{i \in \mathbb{Z}^d} \beta_i^0(x) = 1$, we obtain the result for $t = t_k$, with $k \in \mathcal{I}_{\Delta t}$. The result for every $t \in [0, T]$ follows from (4.23). $\square$

In what follows, given $\varphi: \mathbb{R}^d \to \mathbb{R}$ we set

$$ (4.28) \qquad I^0[\varphi](x) = \sum_{i \in \mathbb{Z}^d} \varphi(x_i) \beta_i^0(x) \quad \text{for all } x \in \mathbb{R}^d. $$

We will need the following estimate in some of the proofs below.

**Lemma 4.2.** Assume (H1)-(H3), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[^3$, and, given $L > 0$, let $\varphi: \mathbb{R}^d \to \mathbb{R}$ be L-Lipschitz. Then, for every $k \in \mathcal{I}_{\Delta t}^*$, we have

$$ (4.28) \qquad \left| \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) dx - \int_{\mathbb{R}^d} \varphi(\Phi_k^{\Delta}[\mu](x)) \mathbf{m}^{\Delta}[\mu](t_k, x) dx \right| \le L \sqrt{d} \Delta x. $$
Proof. Since $\sum_{i \in \mathbb{Z}^d} \beta_i^0(x) = 1$ for all $x \in \mathbb{R}^d$, we have

$$ (4.29) \quad \begin{aligned} & |\varphi(x) - I^0[\varphi](x)| = \left| \sum_{i \in \mathbb{Z}^d} (\varphi(x) - \varphi(x_i)) \beta_i^0(x) \right| \\ & \le \sum_{i \in \mathbb{Z}^d} |\varphi(x) - \varphi(x_i)| \beta_i^0(x) \le \frac{L\sqrt{d}}{2} \Delta x \quad \text{for all } x \in \mathbb{R}^d. \end{aligned} $$

It follows that $\|\varphi - I^0[\varphi]\|_{\infty} \le (L\sqrt{d}/2)\Delta x$ and, hence,

$$ (4.30) \quad \left| \frac{1}{(\Delta x)^d} \int_{E_i} \varphi(x) dx - \varphi(x_i) \right| = \left| \frac{1}{(\Delta x)^d} \int_{E_i} (\varphi(x) - \varphi(x_i)) dx \right| \le \frac{L\sqrt{d}}{2} \Delta x, $$

from which we deduce that, for every $k \in \mathcal{I}_{\Delta t}^*$ and $j \in \mathbb{Z}^d$,

$$ (4.31) \quad \begin{aligned} & \left| \sum_{i \in \mathbb{Z}^d} \frac{1}{(\Delta x)^d} \int_{E_i} \varphi(x) dx \int_{E_j} \beta_i^0(\Phi_k^\Delta[\mu](y)) dy - \int_{E_j} I^0[\varphi](\Phi_k^\Delta[\mu](x)) dx \right| \\ & = \left| \sum_{i \in \mathbb{Z}^d} \left( \frac{1}{(\Delta x)^d} \int_{E_i} \varphi(x) dx - \varphi(x_i) \right) \int_{E_j} \beta_i^0(\Phi_k^\Delta[\mu](y)) dy \right| \\ & \le \frac{L\sqrt{d}}{2} \Delta x \sum_{i \in \mathbb{Z}^d} \int_{E_j} \beta_i^0(\Phi_k^\Delta[\mu](y)) dy \\ & = \frac{L\sqrt{d}}{2} (\Delta x)^{d+1}. \end{aligned} $$

Therefore, by (4.31) and (4.29), we obtain

$$ (4.32) \quad \left| \sum_{i \in \mathbb{Z}^d} \frac{1}{(\Delta x)^d} \int_{E_i} \varphi(x) dx \int_{E_j} \beta_i^0(\Phi_k^\Delta[\mu](y)) dy - \int_{E_j} \varphi(\Phi_k^\Delta[\mu](x)) dx \right| \le L\sqrt{d}(\Delta x)^{d+1}. $$

Finally, from (4.18), (4.32), and Proposition 4.3(iv), we get

$$ (4.33) \quad \begin{aligned} & \left| \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^\Delta[\mu](t_{k+1}, x) dx - \int_{\mathbb{R}^d} \varphi(\Phi_k^\Delta[\mu](x)) \mathbf{m}^\Delta[\mu](t_k, x) dx \right| \\ & = \left| \sum_{i \in \mathbb{Z}^d} \int_{E_i} \varphi(x) dx \frac{1}{(\Delta x)^d} \sum_{j \in \mathbb{Z}^d} m_{k,j} \int_{E_j} \beta_i^0(\Phi_k^\Delta[\mu](y)) dy - \sum_{i \in \mathbb{Z}^d} m_{k,i} \int_{E_i} \varphi(\Phi_k^\Delta[\mu](x)) dx \right| \\ & = \left| \sum_{j \in \mathbb{Z}^d} m_{k,j} \left( \sum_{i \in \mathbb{Z}^d} \frac{1}{(\Delta x)^d} \int_{E_i} \varphi(x) \int_{E_j} \beta_i^0(\Phi_k^\Delta[\mu](y)) dy - \int_{E_j} \varphi(\Phi_k^\Delta[\mu](x)) dx \right) \right| \le L\sqrt{d}(\Delta x), \end{aligned} $$

which shows (4.28). $\square$

In the next result, we study the equicontinuity of the family $\{\mathbf{m}^\Delta[\mu] \mid \Delta \in ]0, \infty[\}^3$, under the condition that $\Delta x$ is, at most, of the order of $\Delta t$.

**Proposition 4.4.** Assume (H1)-(H3), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[\cdot$.
Then, for every $c > 0$, if $\Delta x \le c\Delta t$, we have

$$ (4.34) \quad d_1(\mathbf{m}^\Delta[\mu](t), \mathbf{m}^\Delta[\mu](s)) \le (C_{\text{bf}} + c\sqrt{d})|t-s| \quad \text{for all } t, s \in [0, T], $$
where, for every $t \in [0, T]$, $\mathbf{m}^\Delta[\mu](t) \in \mathcal{P}_1(\mathbb{R}^d)$ denotes the measure $d\mathbf{m}^\Delta[\mu](t)(x) = \mathbf{m}^\Delta[\mu](t, x)dx$.

*Proof.* Let $\varphi \in \text{Lip}_1(\mathbb{R}^d)$ and define $\psi_\varphi: [0, T] \to \mathbb{R}$ by

$$ (4.35) \quad \psi_\varphi(t) = \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^\Delta[\mu](t, x) dx \quad \text{for all } t \in [0, T]. $$

It follows from (4.23) that $\psi_\varphi$ is continuous and affine on every interval $[t_k, t_{k+1}]$ ($k \in \mathcal{I}_{\Delta t}^*$). Thus, $\psi_\varphi \in W^{1,\infty}([0, T])$ and

$$ (4.36) \quad \left\| \frac{d}{dt} \psi_\varphi \right\|_{\infty} = \frac{1}{\Delta t} \max_{k \in \mathcal{I}_{\Delta t}^*} \left| \int_{\mathbb{R}^d} \varphi(x) \left( \mathbf{m}^\Delta[\mu](t_{k+1}, x) - \mathbf{m}^\Delta[\mu](t_k, x) \right) dx \right| $$
In order to estimate the right-hand side of (4.36), fix $k \in \mathcal{I}_{\Delta t}^*$ and notice that, by Lemma 4.2, (4.11), (4.5), Proposition 4.3(iv), and $\Delta x \le c\Delta t$, we have

$$ (4.37) \qquad \begin{aligned} \int_{\mathbb{R}^d} \varphi(x) \left( \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) - \mathbf{m}^{\Delta}[\mu](t_k, x) \right) dx &\le \int_{\mathbb{R}^d} \left( \varphi(\mathbf{\Phi}_k^{\Delta}[\mu](x)) - \varphi(x) \right) \mathbf{m}^{\Delta}[\mu](t_k, x) dx + \sqrt{d}\Delta x \\ &\le \int_{\mathbb{R}^d} \left| \mathbf{\Phi}_k^{\Delta}[\mu](x) - x \right| \mathbf{m}^{\Delta}[\mu](t_k, x) dx + \sqrt{d}\Delta x \\ &\le (C_{\text{bf}} + c\sqrt{d})\Delta t. \end{aligned} $$

Changing $\varphi$ by $-\varphi$ in the previous computation, (4.37) implies that

$$ (4.38) \qquad \left| \int_{\mathbb{R}^d} \varphi(x) \left( \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) - \mathbf{m}^{\Delta}[\mu](t_k, x) \right) dx \right| \le (C_{\text{bf}} + c\sqrt{d})\Delta t. $$

and hence, by (4.36), $\|\frac{d}{dt}\psi_{\varphi}\|_{\infty} \le (C_{\text{bf}} + c\sqrt{d})$. Thus, we deduce that

$$ (4.39) \qquad \int_{\mathbb{R}^d} \varphi(x) \left( \mathbf{m}^{\Delta}[\mu](t, x) - \mathbf{m}^{\Delta}[\mu](s, x) \right) dx \le (C_{\text{bf}} + c\sqrt{d})|t - s| \quad \text{for all } t, s \in [0, T] $$

and (4.34) follows from (2.2). $\square$

The following result states a stability property for $\mathbf{m}^{\Delta}[\mu]$ which is analogous to the one in Proposition 4.1 for $\mathbf{m}^{\Delta}[\mu]$.

**Proposition 4.5.** Assume (H1)-(H3), let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[$. Then for every $c_1, c_2 > 0$, there exists $\tilde{C} > 0$, independent of $\mu$, such that, if $\Delta$ is small enough, $\Delta x \le c_1 \Delta t$, and $\Delta t \le c_2 \varepsilon^2$, we have $\mathbf{m}^{\Delta}[\mu](t, \cdot) \in L^p(\mathbb{R}^d)$ for all $t \in [0, T]$ and

$$ (4.40) \qquad \|\mathbf{m}^{\Delta}[\mu](t, \cdot)\|_{L^p(\mathbb{R}^d)} \le \tilde{C} \|m_0^*\|_{L^p(\mathbb{R}^d)} \quad \text{for all } t \in [0, T]. $$

*Proof.* Let $c_1, c_2 > 0$, suppose that $\Delta x \le c_1 \Delta t$, $\Delta t \le c_2 \varepsilon^2$, fix $k \in \mathcal{I}_{\Delta t}^*$, and let $\tilde{C}^* > 0$ be as in Lemma 4.1. Then, by (4.5), there exists $R > 0$, independent of $(\mu, \Delta, k)$, such that $(\mathbf{\Phi}_k^{\Delta}[\mu])^{-1}(\tilde{B}_{\infty}(0, \tilde{C}^*)) \subset B_{\infty}(0, R)$. The regularity of $H$ and estimate (3.16), with $\ell = 2$, yield the existence of $C_R > 0$, independent of $(\mu, \Delta, k)$, such that, for every $x \in B_{\infty}(0, R)$, the norm of the matrix

$$ (4.41) \qquad D_x(D_pH(x, D_x v^{\Delta}[\mu](t_k, x))) = D_p^2 H(x, D_x v^{\Delta}[\mu](t_k, x)) D_x^2 v^{\Delta}[\mu](t_k, x) + D_{x_p}^2 H(x, D_x v^{\Delta}[\mu](t_k, x)), $$

induced by the 2-norm in $\mathbb{R}^d$, is bounded by $C_R/\varepsilon$. In particular, $D_pH(\cdot, D_x v^{\Delta}[\mu](t_k, \cdot))$ is $(C_R/\varepsilon)$-Lipschitz on $B_{\infty}(0, R)$. Thus, expression (4.11) and the inequality $\Delta t/\varepsilon \le c_2\varepsilon$ imply that, if $\varepsilon$ is small enough, there exists $C_1 > 0$, independent of $(\mu, \Delta, k)$, such that

$$ (4.42) \qquad |\mathbf{\Phi}_k^{\Delta}[\mu](x) - \mathbf{\Phi}_k^{\Delta}[\mu](y)| \ge C_1|x - y| \quad \text{for all } x, y \in B_{\infty}(0, R), $$

which implies that $\mathbf{\Phi}_k^{\Delta}[\mu]$ is injective on $B_{\infty}(0, R)$, and, denoting by $I_d$ the $d \times d$ identity matrix,

$$ (4.43) \qquad D_x \mathbf{\Phi}_k^{\Delta}[\mu](x) = I_d - \Delta t D_x [D_p H(x, D_x v^{\Delta}[\mu](t_k, x))] $$

is invertible for $x \in B_{\infty}(0, R)$. In particular, $\mathbf{\Phi}_k^{\Delta}[\mu]$ is a diffeomorphism of $B_{\infty}(0, R)$ onto $\mathbf{\Phi}_k^{\Delta}[\mu](B_{\infty}(0, R))$. Let us suppose first that $p \in ]1, \infty[$. By the change of variable formula, we have

$$ (4.44) \qquad \int_{\mathbb{R}^d} \left( \mathbf{m}^{\Delta}[\mu](t_{k+1}, \mathbf{\Phi}_k^{\Delta}[\mu](x)) \right)^p dx = \int_{B_{\infty}(0, R)} \left( \mathbf{m}^{\Delta}[\mu](t_{k+1}, \mathbf{\Phi}_k^{\Delta}[\mu](x)) \right)^p dx \\ = \int_{\mathbf{\Phi}_k^{\Delta}[\mu](B_{\infty}(0, R))} \left( \mathbf{m}^{\Delta}[\mu](t_{k+1}, y) \right)^p \left| \det \left( D_x \mathbf{\Phi}_k^{\Delta}[\mu] \left( \mathbf{\Phi}_k^{\Delta}[\mu]^{-1}(y) \right) \right) \right|^{-1} dy. $$

Using again that the norm of $D_x(D_pH(\cdot, D_x v^{\Delta}[\mu](t_k, \cdot)))$ is bounded by $C_R/\varepsilon$ on $B_{\infty}(0, R)$, relation (4.43) and a Taylor expansion for the determinant yield the existence of $C_2 > 0$, independent of $(\mu, \Delta, k)$, such that

$$ \left| \det \left( D_x \mathbf{\Phi}_k^{\Delta}[\mu](x) \right) - \left( 1 - \Delta t \operatorname{Tr} \left( D_x \left( D_p H(x, D_x v^{\Delta}[\mu](t_k, x)) \right) \right) \right) \right| \le C_2 (\Delta t / \varepsilon)^2 \quad \text{for all } x \in B_{\infty}(0, R), $$
where, given $B \in \mathbb{R}^{d \times d}$, $\mathrm{Tr}(B)$ denotes its trace. In turn, we get the existence of $C_3 > 0$, independent of $(\mu, \Delta, k)$, such that

$$ (4.45) \qquad \left| \det (D_x \Phi_k^\Delta[\mu](x)) \right|^{-1} - \left(1 + \Delta t \mathrm{Tr} (D_x (D_p H(x, D_x v^\Delta[\mu](t_k, x))))\right) \le C_3 (\Delta t / \epsilon)^2, $$

for all $x \in B_\infty(0, R)$. Since $\Delta x \le c_1 c_2 \epsilon^2$, Proposition 3.5 implies that $D_x^2 v^\Delta[\mu](t_k, x) - \tilde{C}_{hb}(1 + (c_1 c_2)^2) I_d$ is negative semidefinite. Using that $H(\cdot, \cdot)$ is of class $C^2$ and convex with respect to its second argument, it follows from (4.41) and [18, Lemma 1.6.4] that there exists $C_4 > 0$, independent of $(\mu, \Delta, k)$, such that

$$ (4.46) \qquad \mathrm{Tr}(D_x [D_p H(x, D_x v^\Delta[\mu](t_k, x))]) \le C_4 \quad \text{for all } x \in B_\infty(0, R), $$

which, together with (4.45), yields

$$ (4.47) \qquad \left| \det (D_x \Phi_k^\Delta[\mu](x)) \right|^{-1} \le 1 + C_5 \Delta t \quad \text{for all } x \in B_\infty(0, R), $$

where $C_5 = C_4 + C_3 c_2$. Therefore, by (4.44), we get

$$ (4.48) \qquad \int_{\mathbb{R}^d} \left( \mathbf{m}^{\Delta}[\mu](t_{k+1}, \Phi_k^{\Delta}[\mu](x)) \right)^{\mathbf{p}} dx \le (1 + C_5 \Delta t) \int_{\mathbb{R}^d} \left( \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) \right)^{\mathbf{p}} dx. $$

Setting $\mathbf{p}^* = \mathbf{p}/(\mathbf{p}-1)$, it follows from (4.25) and Hölder's inequality that

$$ \begin{align} \|\mathbf{m}^{\Delta}[\mu](t_{k+1}, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)} &= \int_{\mathbb{R}^d} (\mathbf{m}^{\Delta}[\mu](t_{k+1}, x))^{\mathbf{p}-1} \mathbf{m}^{\Delta}[\mu](t_{k+1}, x) dx \nonumber \\ &= \int_{\mathbb{R}^d} (\mathbf{m}^{\Delta}[\mu](t_{k+1}, \Phi_k^{\Delta}[\mu](x)))^{\mathbf{p}-1} \mathbf{m}^{\Delta}[\mu](t_k, x) dx \nonumber \\ &\le \left( \int_{\mathbb{R}^d} (\mathbf{m}^{\Delta}[\mu](t_{k+1}, \Phi_k^{\Delta}[\mu](x)))^{\mathbf{p}} dx \right)^{\frac{1}{\mathbf{p}*}} \|\mathbf{m}^{\Delta}[\mu](t_k, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)} \nonumber \\ &\le (1 + C_5 \Delta t)^{\frac{1}{\mathbf{p}*}} \|\mathbf{m}^{\Delta}[\mu](t_{k+1}, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)}^{\frac{\mathbf{p}}{\mathbf{p}*}} \|\mathbf{m}^{\Delta}[\mu](t_k, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)} \nonumber \\ &\le (1 + C_5 \Delta t) \|\mathbf{m}^{\Delta}[\mu](t_{k+1}, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)}^{\mathbf{p}-1} \|\mathbf{m}^{\Delta}[\mu](t_k, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)}. \tag{4.49} \end{align} $$

In turn, we deduce that

$$ (4.50) \qquad \|\mathbf{m}^{\Delta}[\mu](t_{k+1}, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)} \le (1 + C_5 \Delta t) \|\mathbf{m}^{\Delta}[\mu](t_k, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)}. $$

By (4.19) and Jensen's inequality, we have $\|\mathbf{m}^{\Delta}[\mu](0, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)} \le \|m_0^*\|_{L^{\mathbf{p}}(\mathbb{R}^d)}$, and hence

$$ \begin{align} \|\mathbf{m}^{\Delta}[\mu](t_{k+1}, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)} &\le (1 + C_5 \Delta t)^N \|\mathbf{m}^{\Delta}[\mu](0, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)} \nonumber \\ &\le e^{C_5 T} \|m_0^*\|_{L^{\mathbf{p}}(\mathbb{R}^d)}, \tag{4.51} \end{align} $$

which, by (4.23), shows (4.40), with $\tilde{C} = e^{C_5 T}$. If $\mathbf{p} = \infty$, then (4.40) holds for every $\mathbf{p}' \in ]1, \in[$. Noticing that $\tilde{C}$ is independent of $\mathbf{p}'$ and that, by Proposition 4.3(ii), for every $t \in [0, T]$, the support of $\mathbf{m}^{\Delta}[\mu](t, \cdot)$ is contained in $\bar{B}_{\infty}(0, \tilde{C}^*)$, (4.40) for $\mathbf{p} = \infty$ follows by letting $\mathbf{p}' \to \infty$. $\square$

The next result provides the analogous for $\mathbf{m}^{\Delta}[\mu]$ of Proposition 4.2 for $m^{\Delta}[\mu]$.

**Proposition 4.6.** Assume (H1)-(H3), let $(\mu_n)_{n \in \mathbb{N}} \subset C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$, and let $((\Delta t_n, \Delta x_n, \varepsilon_n))_{n \in \mathbb{N}} \subset (0, \infty)^3$. Set $\Delta_n = (\Delta t_n, \Delta x_n, \varepsilon_n)$ and let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. Suppose that, as $n \to \infty$, $\mu_n \to \mu$, $\Delta_n \to 0$, $\Delta x_n = o(\Delta t_n)$, and $\Delta t_n = O(\varepsilon_n^2)$. Then, up to some subsequence, the following hold:

(i) $(v^{\Delta_n}[\mu^n])_{n \in \mathbb{N}}$ converges to $v[\mu]$, uniformly over compact subsets of $[0, T] \times \mathbb{R}^d$, and, for every $K \subset [0, T] \times \mathbb{R}^d$ compact and $q \in [1, \infty[, (\mathbb{I}_K D_x v^{\Delta_n}[\mu_n])_{n \in \mathbb{N}}$ converges to $\mathbb{I}_K D_x v[\mu]$ in $L^q([0, T] \times \mathbb{R}^d)$.

(ii) $(\mathbf{m}^{\Delta_n}[\mu_n])_{n \in \mathbb{N}}$ converges in $C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ towards a solution $m \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d)) \cap L^p([0, T] \times \mathbb{R}^d)$ to (2.30). Moreover, the convergence also hold weakly in $L^p([0, T] \times \mathbb{R}^d)$, if $\mathbf{p} < \infty$, and weakly* in $L^\infty([0, T] \times \mathbb{R}^d)$, if $\mathbf{p} = \infty$. In addition, there exists $\tilde{C} > 0$ such that

$$ (4.52) \qquad \|m(t, \cdot)\|_{L^{\mathbf{p}}(\mathbb{R}^d)} \le \tilde{C} \|m_0^*\|_{L^{\mathbf{p}}(\mathbb{R}^d)} \quad \text{for all } t \in [0, T]. $$
Proof. Let us set $v^n := v^{\Delta_n, \varepsilon_n}[\mu^n]$, $\mathbf{m}^n = \mathbf{m}^{\Delta_n, \varepsilon_n}[\mu_n]$, and $\Phi_k^n = \Phi^{\Delta_n, \varepsilon_n}[\mu_n]$ for all $k \in \mathcal{I}_{\Delta t_n}^*$.

(i): This corresponds to Proposition 4.2(i).

(ii): It follows from Proposition 4.3(ii), [9, Proposition 7.1.5], Proposition 4.4, and the Arzelá-Ascoli theorem, that there exists $m \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ such that, as $n \to \infty$ and up to some subsequence, $\mathbf{m}^n \to m$ in $C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. Moreover, by Proposition 4.5, the convergence holds weakly, if $p < \infty$, and weakly* in $L^\infty([0, T] \times \mathbb{R}^d)$, if $p = \infty$, and $m$ satisfies (4.52). It remains to show that $m$ solves (2.30). Let $t \in [0, T]$ and let $k(n) \in \mathcal{I}_{\Delta t_n}^*$ be such that $t \in [t_{k(n)}, t_{k(n)+1}]$. For every $\varphi \in C_0^\infty(\mathbb{R}^d)$, we have

$$ (4.53) \quad \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^n(t_{k(n)}, x) \,dx = \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^n(0, x) \,dx + \sum_{k=0}^{k(n)-1} \int_{\mathbb{R}^d} \varphi(x) (\mathbf{m}^n(t_{k+1}, x) - \mathbf{m}^n(t_k, x)) \,dx. $$

Let $k \in \mathcal{I}_{\Delta t_n}^*$. Since (4.11) and (4.5) yield $|\Phi_k^n(x) - x| = O(\Delta t_n)$ for all $x \in \text{supp}(\varphi)$, by Lemma 4.2, a Taylor expansion, and Proposition 4.3(iv), we have

$$ \begin{align} \int_{\mathbb{R}^d} \varphi(x) (\mathbf{m}^n(t_{k+1}, x) - \mathbf{m}^n(t_k, x)) dx &= \int_{\mathbb{R}^d} (\varphi(\Phi_k^n(x)) - \varphi(x)) \mathbf{m}^n(t_k, x) dx + O(\Delta x_n) \nonumber \\ &= -\Delta t_n \int_{\mathbb{R}^d} \langle D\varphi(x), D_p H(x, D_x v^n(t_k, x)) \rangle \mathbf{m}^n(t_k, x) dx + O(\Delta x_n) \nonumber \\ &\quad + O((\Delta t_n)^2), \tag{4.54} \end{align} $$

which, combined with (4.53), yields

$$ \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^n(t_{k(n)}, x) dx = \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^n(0, x) dx - \Delta t_n \sum_{k=0}^{k(n)-1} \int_{\mathbb{R}^d} \langle D\varphi(x), D_p H(x, D_x v^n(t_k, x)) \rangle \mathbf{m}^n(t_k, x) dx \\ + O \left( \frac{\Delta x_n}{\Delta t_n} + \Delta t_n \right). \quad (4.55) $$

Since $\varphi$ has a compact support, it follows from (3.16), with $\ell = 2$, that there exists $C_\varphi > 0$ such that $\langle D\varphi(\cdot), D_p H(\cdot, D_x v^n(t_k, \cdot)) \rangle$ is $(C_\varphi/\varepsilon_n)$-Lipschitz. Thus, by Proposition 4.4, for every $k \in \mathcal{I}_{\Delta t_n}^*$, we have

$$ (4.56) \quad \left| \int_{\mathbb{R}^d} \langle D\varphi(x), D_p H(x, D_x v^n(t_k, x)) \rangle (\mathbf{m}^n(s, x) - \mathbf{m}^n(t_k, x)) dx \right| = O \left( \frac{\Delta t_n}{\varepsilon_n} \right) \quad \text{for all } s \in [t_k, t_{k+1}]. $$

Recalling that $D_x v^n(s, x) = D_x v^n(t_k, x)$ for all $s \in [t_k, t_{k+1}]$ and $x \in \mathbb{R}^d$, we obtain

$$ \begin{align} \int_{\mathbb{R}^d} \Delta t_n \langle D\varphi(x), D_p H(x, D_x v^n(t_k, x)) \rangle \mathbf{m}^n(t_k, x) dx & \nonumber \\ &= \int_{t_k}^{t_{k+1}} \int_{\mathbb{R}^d} \langle D\varphi(x), D_p H(x, D_x v^n(s, x)) \rangle \mathbf{m}^n(s, x) dx ds + O \left( \frac{(\Delta t_n)^2}{\varepsilon_n} \right) \end{align} $$

and hence, in view of (4.55), we deduce that, for $n$ large enough,

$$ (4.58) \quad \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^n(t_{k(n)}, x) dx \\ = \int_{\mathbb{R}^d} \varphi(x) \mathbf{m}^n(0, x) dx - \int_0^T \int_{\mathbb{R}^d} \mathbb{I}_{[0, t_{k(n)}]} \langle D\varphi(x), D_p H(x, D_x v^n(s, x)) \rangle \mathbf{m}^n(s, x) dx ds + O \left( \frac{\Delta x_n}{\Delta t_n} + \frac{\Delta t_n}{\varepsilon_n} \right). $$

Finally, by (i),

$$ (4.59) \quad \mathbb{I}_{[0, t_{k(n)}]}(\cdot) \langle D\varphi(\cdot), D_p H(x, D_x v^n(\cdot, \cdot)) \rangle \underset{n \to \infty}{\longrightarrow} \mathbb{I}_{[0, t]}(\cdot) \langle D\varphi(\cdot), D_p H(x, D_x v[m](\cdot, \cdot)) \rangle, $$

in $L^q([0, T] \times \mathbb{R}^d)$, for every $q \in [1, \infty[$, and, hence, we can pass to the limit in (4.58) to obtain that $m$ satisfies (2.30). $\square$
5. A LAGRANGE-GALERKIN SCHEME FOR THE MEAN FIELD GAMES SYSTEM

In this section, we combine the schemes discussed in Sections 3 and 4 to obtain a scheme for system (MFG) and we provide a convergence result.

Let $\Delta = (\Delta t, \Delta x, \varepsilon) \in ]0, \infty[$, let $\tilde{C}^* > 0$ be as in Lemma 4.1, and define

$$ (5.1) \qquad \mathcal{D}^{\Delta t, \Delta x} = \left\{ \mu = (\mu_{k,i}) \mid \mu_{k,i} \ge 0, \sum_{j \in \mathbb{Z}^d} \mu_{k,j} (\Delta x)^d = 1 \text{ for all } k \in \mathcal{I}_{\Delta t}, i \in \mathcal{I}_{\Delta x} \right\}, $$

where $\mathcal{I}_{\Delta x}$ is defined in (4.22). Notice that $\mathcal{D}^{\Delta t, \Delta x}$ is a convex and compact subset of $\mathbb{R}^{(N_{\Delta t}+1) \times (2N_{\Delta x}+1)^d}$.
Given $\mu \in \mathcal{D}^{\Delta t, \Delta x}$ define $\tilde{\mu} \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ as

$$ (5.2) \quad \mathrm{d}\tilde{\mu}(t)(x) = \left(\frac{t-t_k}{\Delta t}\right) \sum_{i \in \mathcal{I}_{\Delta x}} \mu_{k+1,i} \beta_i^0(x) \mathrm{d}x + \left(\frac{t_{k+1}-t}{\Delta t}\right) \sum_{i \in \mathcal{I}_{\Delta x}} \mu_{k,i} \beta_i^0(x) \mathrm{d}x $$
$$ \text{for all } k \in \mathcal{I}_{\Delta t}^*, t \in [t_k, t_{k+1}.] $$

The discretization of (MFG) that we propose is the following: find $\mu \in \mathcal{D}^{\Delta t, \Delta x}$ such that

$$ (\mathrm{MFG}^{\Delta}) \qquad \mu_{k,i} = \mathbf{m}^{\Delta}[\tilde{\mu}](t_k, x_i) \quad \text{for all } k \in \mathcal{I}_{\Delta t}^{*}, i \in \mathcal{I}_{\Delta x}, $$

where we recall that $\mathbf{m}^{\Delta}[\tilde{\mu}]$ is defined in (4.23).

**Theorem 5.1.** Assume that (H1)-(H3) hold. Then, if $\Delta t/\varepsilon$ is small enough, system (MFG$^{\Delta}$) admits at least one solution.

*Proof.* Consider the application $T: \mathcal{D}^{\Delta t, \Delta x} \to \mathbb{R}^{(N_{\Delta t}+1) \times (2N_{\Delta x}+1)^d}$ defined by

$$ (T(\mu))_{k,i} = \mathbf{m}^{\Delta}[\tilde{\mu}](t_k, x_i) \quad \text{for all } k \in \mathcal{I}_{\Delta t}^{*}, i \in \mathcal{I}_{\Delta x}. $$

It follows from Proposition 4.3(i),(ii),(iv) that $T(\mathcal{D}^{\Delta t, \Delta x}) \subseteq \mathcal{D}^{\Delta t, \Delta x}$. Moreover, if $(\mu_n)_{n \in \mathbb{N}} \subset \mathcal{D}^{\Delta t, \Delta x}$ converges to $\mu$ then the continuity of $L$, $F$, and $G$, imply that, as $n \to \infty$, $v_{k,i}^{\Delta t, \Delta x}[\mu_n] \to v_{k,i}^{\Delta t, \Delta x}[\mu]$ for all $k \in \mathcal{I}_{\Delta t}$ and $i \in \mathbb{Z}^d$. Thus, $(v_{k,i}^{\Delta t, \Delta x}[\mu_n])_{n \in \mathbb{N}}$, defined in (3.9), converges to $v_{k,i}^{\Delta t, \Delta x}[\mu]$ pointwise and hence, by Lebesgue's dominated convergence, the sequence $(v_{k,i}^{\Delta}[\mu_n])_{n \in \mathbb{N}}$, defined in (3.17), satisfies that $v_{k,i}^{\Delta}[\mu_n] \to v_{k,i}^{\Delta}[\mu]$ and $D_x v_{k,i}^{\Delta}[\mu_n] \to D_x v_{k,i}^{\Delta}[\mu]$ pointwise. Consequently, given $k \in \mathcal{I}_{\Delta t}^*$, it follows from (4.11) that $\Phi_k^{\Delta}[\mu_n] \to \Phi_k^{\Delta}[\mu]$ pointwise. In particular, $\beta_i^0(\Phi_k^{\Delta}[\mu_n](x)) \to \beta_i^0(\Phi_k^{\Delta}[\mu](x))$ for all $x \in \mathbb{R}^d \setminus (\Phi_k^{\Delta}[\mu]^{-1}(\partial E_i))$. If $R > 0$ is as in the proof of Proposition 4.5 and $\Delta t/\varepsilon$ is small enough, we have that $\Phi_k^{\Delta}[\mu]$ is a diffeomorphism of $B_\infty(0, R)$ onto $\Phi_k^{\Delta}[\mu](B_\infty(0, R))$. Therefore, since $\mathcal{L}^d(\partial E_i) = 0$, we have $\mathcal{L}^d(\Phi_k^{\Delta}[\mu]^{-1}(\partial E_i)) = 0$ and hence $\beta_i^0(\Phi_k^{\Delta}[\mu_n](x)) \to \beta_i^0(\Phi_k^{\Delta}[\mu](x))$ for almost every $x \in \mathbb{R}^d$. Therefore, by Lebesgue's dominated convergence,

$$ \int_{E_j} \beta_i^0(\Phi_k^{\Delta}[\mu_n](x)) \mathrm{d}x \xrightarrow{n \to \infty} \int_{E_j} \beta_i^0(\Phi_k^{\Delta}[\mu](x)) \mathrm{d}x \quad \text{for all } k \in \mathcal{I}_{\Delta t}^*, i, j \in \mathbb{Z}^d. $$

Altogether, it follows from (4.18) that, as $n \to \infty$, $T(\mu_n) \to T(\mu)$, i.e. $T$ is continuous. Finally, the existence of a solution to (MFG$^{\Delta}$), i.e. of a fixed point of $T$, follows from Brouwer's fixed-point theorem. $\square$

In the next result we provide our main result, which shows the convergence, up to some subsequence, of solutions to (MFG$^{\Delta}$) towards a solution to (MFG).

**Theorem 5.2.** Assume (H1)-(H3), let $((\Delta t_n, \Delta x_n, \varepsilon_n))_{n \in \mathbb{N}} \subset ]0, \infty[$, and set $\Delta_n = (\Delta t_n, \Delta x_n, \varepsilon_n)$. Suppose that, as $n \to \infty$, $\Delta_n \to 0$, $\Delta x_n = o(\Delta t_n)$, and $\Delta t_n = O(\varepsilon_n^2)$. For every $n$, large enough, let $m^n \in \mathcal{S}^{\Delta_n}$ be a solution to (MFG$^{\Delta_n}$), define $\tilde{m}^n$ by (5.2), and set $v^n = v^{\Delta_n}[\tilde{m}^n]$. Then there exists a solution $(v^*, m^*)$ to (MFG) such that, up to some subsequence, the following hold:

(i) $(v^n)_{n \in \mathbb{N}}$ converges to $v^*$, uniformly over compact subsets of $[0, T] \times \mathbb{R}^d$.

(ii) $(\tilde{m}^n)_{n \in \mathbb{N}}$ converges in $C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ towards $m^*$. Moreover, the convergence also hold weakly in $L^p([0, T] \times \mathbb{R}^d)$, if $p < \infty$, and weakly* in $L^\infty([0, T] \times \mathbb{R}^d)$, if $p = \infty$. In addition, there exists $\tilde{C} > 0$ such that

$$ (5.3) \qquad \|m^*(t, \cdot)\|_{L^p(\mathbb{R}^d)} \le \tilde{C} \|m_0^*\|_{L^p(\mathbb{R}^d)} \quad \text{for all } t \in [0, T]. $$
Proof. For all $n \in \mathbb{N}$, large enough, we have $\tilde{m}^n = m^{\Delta_n}[\tilde{m}^n]$. Thus, Proposition 4.3(ii), [9, Proposition 7.1.5], Proposition 4.4, and the Arzelá-Ascoli theorem, imply the existence of $m^* \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ and subsequence, still labelled by $n$, such that $(\tilde{m}^n)_{n \in \mathbb{N}}$ converges to $m^*$ in $C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$. It follows from Proposition 4.6 that $m^*$ solves (2.30), with $\mu = m^*$, i.e. $(v[m^*], m^*)$ solves (MFG). Therefore, assertions (i)-(ii) follow from the corresponding assertions in Proposition 4.6. $\square$

**Remark 5.1.** (i) Theorem 5.2 shows, in particular, that system (MFG) admits at least one solution $(v^*, m^*)$. If the solution to (MFG) is unique, then the entire sequence $(v^n, m^n)$ converges to $(v^*, m^*)$ (see Theorem 2.4 for a sufficient condition ensuring uniqueness).

(ii) The condition on $\Delta x_n = o(\Delta t_n)$ in Theorem 5.2 is stronger than the condition $(\Delta x_n)^2 = o(\Delta t_n)$ needed for convergence, when the space dimension is equal to one, in the scheme studied in [20] (see also [27]). This can be explained by the estimate (4.28) in Lemma 4.2, which seems difficult to improve, even if $\varphi$ is smooth, and it is in compliance with Assumption 3.1 in [48], which plays an important role in the LG approximation of continuity equations with Lipschitz vector fields.

## 6. NUMERICAL RESULTS

In this section, given $\Delta = (\Delta t, \Delta x, \varepsilon)$, we use $(\text{MFG}^\Delta)$ to approximate the solutions to two first order MFGs systems. In order to obtain an implementable version, we need to approximate the integrals in the LG scheme (4.18)-(4.19). We consider two methods. In the first one, the integrals are approximated by numerical quadrature, while, in the second one, we use the so-called *area weighting* technique, introduced in [43] and recalled in Section 6.1 below.

In the first example, the state dimension is equal to one and the data of the MFGs system does not satisfy some of the assumptions in Section 2.2. On the other hand, the PDE system admits an explicit solution, which allows to compare the quadrature and area weighting methods to solve $(\text{MFG}^\Delta)$. For comparable accuracies, the area weighting method is less expensive than the quadrature method and, hence, we use the former in order to treat the second example, where the state dimension is equal to two and no explicit solution is known. Let us point out that the data of the second example fulfills all the assumptions in Section 2.2.

We solve $(\text{MFG}^\Delta)$ heuristically by fixed point iterations that are stopped as soon as the uniform norm of the difference between two consecutive iterates is smaller than a given threshold $\tau$, which in the simulations is set to $10^{-3}$. In particular, we use the classical Picard iterations in the first test, as in [20], and Picard iterations with damping parameter 0.5 in the second test, as in [27].

**6.1. Area-weighted LG approximation.** Let $\mu \in C([0, T]; \mathcal{P}_1(\mathbb{R}^d))$ and consider the continuity equation (4.1). The main idea of the area-weighting technique is to replace, for each $k \in \mathcal{I}_{\Delta x}^*$, the local nonlinear discrete flow $E_i \mapsto x \mapsto \Phi_k^\Delta[\mu](x) \in \mathbb{R}^d$, defined by (4.11), by the local affine approximation

$$ (6.1) \qquad E_i \mapsto x \mapsto \bar{\Phi}_k^\Delta[\mu](x) = x - \Delta t D_p H(x_i, D_x v^\Delta[\mu](t_k, x_i)) \in \mathbb{R}^d. $$

Notice that $\bar{\Phi}_k^\Delta[\mu](x) = x - x_i + \Phi_k^\Delta[\mu](x_i)$ for all $x \in E_i$. Under this approximation, we can compute the integrals in (4.18)-(4.19) explicitly. Indeed, for all $i = (i_1, \dots, i_d) \in \mathbb{Z}^d$ and $l = 1, \dots, d$, let us set $T_{i,l} = [(x_i)_l - \Delta x/2, (x_i)_l + \Delta x/2]$, and observe that

$$ (6.2) \qquad \beta_i^0(y) = \prod_{l=1}^{d} \mathbb{I}_{T_{i,l}}(y_l) \quad \text{for all } y = (y_1, \dots, y_d) \in \mathbb{R}^d. $$

It follows from (6.1) and (6.2) that, for every $i, j \in \mathbb{Z}^d$, we have

$$ \begin{align} \int_{E_j} \beta_i^0 (\bar{\Phi}_k^\Delta [\mu](y)) dy &= \int_{E_j} \beta_i^0 (y - x_j + \Phi_k^\Delta [\mu](x_j)) dy \tag{6.3} \\ &= \prod_{l=1}^d \int_{(x_j)_l - \Delta x/2}^{(x_j)_l + \Delta x/2} \mathbb{I}_{T_{i,l}} (y_l - (x_j)_l + (\Phi_k^\Delta [\mu](x_j))_l) dy_l = \prod_{l=1}^d \int_{(\Phi_k^\Delta [\mu](x_j))_l - \Delta x/2}^{(\Phi_k^\Delta [\mu](x_j))_l + \Delta x/2} \mathbb{I}_{T_{i,l}} (y_l) dy_l \nonumber \\ &= \prod_{l=1}^d \mathcal{L}^1 \left( [(x_i)_l - \Delta x/2, (x_i)_l + \Delta x/2] \cap [(\Phi_k^\Delta [\mu](x_j))_l - \Delta x/2, (\Phi_k^\Delta [\mu](x_j))_l + \Delta x/2] \right). \nonumber \end{align} $$
On the other hand, for every $l = 1, \dots, d$, it follows from (3.2) that

$$ \begin{aligned} & \mathcal{L}^1 \left( [(x_i)_l - \Delta x/2, x_i)_l + \Delta x/2] \cap [(\Phi_k^\Delta[\mu](x_j))_l - \Delta x/2, (\Phi_k^\Delta[\mu](x_j))_l + \Delta x/2] \right) \\ &= \begin{cases} \Delta x + (\Phi_k^\Delta[\mu](x_j))_l - (x_i)_l & \text{if } (\Phi_k^\Delta[\mu](x_j))_l \in [(x_i)_l - \Delta x, (x_i)_l], \\ \Delta x + (x_i)_l - (\Phi_k^\Delta[\mu](x_j))_l & \text{if } (\Phi_k^\Delta[\mu](x_j))_l \in (x_i)_l, (x_i)_l + \Delta x, \\ 0 & \text{otherwise,} \end{cases} \\ &= \Delta x \hat{\beta} ((\Phi_k^\Delta[\mu](x_j))_l / \Delta x - i_l), \end{aligned} $$

which, combined with (3.1) and (6.3), yields

$$ (6.4) \qquad \frac{1}{(\Delta x)^d} \int_{E_j} \beta_i^0 (\bar{\Phi}_k^\Delta[\mu](y)) dy = \beta_i^1 (\Phi_k^\Delta[\mu](x_j)). $$

Thus, replacing $\Phi_k^\Delta[\mu]$ by $\bar{\Phi}_k^\Delta[\mu]$ in (4.18) and, for every $i \in \mathbb{Z}^d$, denoting by $m_{0,i}^*$ any approximation of $\int_{E_i} m_0^*(x)dx/(\Delta x)^d$, we obtain the following area-weighted LG version of (4.18)-(4.19):

$$ \begin{aligned} \bar{m}_{k+1,i} &= \sum_{j \in \mathbb{Z}^d} \bar{m}_{k,j} \beta_i^1 (\Phi_k^\Delta[\mu](x_j)) & \text{for all } k \in \mathcal{I}_{\Delta t}^*, i \in \mathbb{Z}^d, \\ \bar{m}_{0,i} &= m_{0,i}^* & \text{for all } i \in \mathbb{Z}^d. \end{aligned} \qquad (6.5) $$

**Remark 6.1.** Notice that (6.5) corresponds to the scheme proposed in [20] for the continuity equation (4.1). Therefore, the latter can be seen as an area-weighted version of the LG scheme of Section 4.

**6.2. Non-local MFG with analytical solution.** We consider system (MFG) with a quadratic Hamiltonian $H(x, p) = \frac{\nu^2}{2}$, coupling terms

$$ (6.6) \qquad F(x, \nu) = \frac{1}{2} \left( x - \int_{\mathbb{R}^d} y \, d\nu(y) \right)^2, \quad G(x, \nu) = 0, $$

and initial data $m_0^*$ given by the distribution of a $d$-dimensional Gaussian random variable with mean $\mu^* \in \mathbb{R}^d$ and covariance matrix $\Sigma_0 \in \mathbb{R}^{d \times d}$ assumed, for simplicity, to be diagonal. Notice that the coupling term $F$ in (6.6) and the initial distribution $m_0^*$ do not satisfy assumptions (H2) and (H3), respectively. On the other hand, the MFG system admits in this case an explicit solution, which allows to compare the performance of quadrature and area-weighting methods to approximate the continuity equation. Indeed, setting

$$ \Pi(t) = \left( \frac{e^{2T-t} - e^t}{e^{2T-t} + e^t} \right) I_d, \quad s(t) = -\Pi(t)\mu^*, \quad c(t) = \frac{1}{2} \langle \Pi(t)\mu^*, \mu^* \rangle \quad \text{for all } t \in [0, T], $$

and arguing as in [16, Section 5.2], one finds that (MFG) admits a unique solution $(v^*, m^*)$ given by

$$ v^*(t, x) = \frac{1}{2} \langle \Pi(t)x, x \rangle + \langle s(t), x \rangle + c(t) \quad \text{for all } (t, x) \in [0, T] \times \mathbb{R}^d $$

and, for every $t \in [0, T]$, $m^*(t)$ is the joint distribution of $d$ independent Gaussian random variables $\{X_l(t) | l = 1, \dots, d\}$ with means $\mu_l(t)$ and variances $\sigma_l^2(t)$ ($l = 1, \dots, d$), given by

$$ \mu_l(t) = \mu_l^* \quad \text{and} \quad \sigma_l^2(t) = \left( \frac{e^{2T-t} + e^t}{e^{2T} + 1} \right)^2 (\Sigma_0)_{l,l}. $$

In the numerical test, we take $T = 0.25$, $d = 1$, $\mu^* = 0.1$, and $\Sigma_0 = 0.105$. Since the exact solution $m^*(t)$ does not have a compact support, we approximate the system on the bounded domain $\mathcal{O} = ]-2, 2[$ and we impose Dirichlet boundary conditions at $x = -2$ and $x = 2$, which are equal to the values of the exact solution at these points. In order to implement the latter, we proceed as in [29, Section 5.1.5].

We test our scheme for different values of $\Delta x$, the time step is chosen as $\Delta t = (\Delta x)^{2/3}/2$, and the mollifier in (3.17) is defined with $\mathbb{R} \ni x \mapsto \rho(x) = e^{-x^2/2}/\sqrt{2\pi} \in \mathbb{R}$ and $\varepsilon = \sqrt{\Delta t}$. We denote by $(v^\Delta, m^\Delta)$ and $(\bar{v}^\Delta, \bar{m}^\Delta)$ the approximations of solutions to (MFG$^\Delta$) obtained by estimating the integrals in (4.18) by numerical quadrature and by the area-weighting method, respectively. For the numerical quadrature of the integrals in the computation of $(v^\Delta, m^\Delta)$, we divide each interval $E_i$ into $\lfloor 4/\Delta x \rfloor$ subintervals
and we use the midpoint rule on each one of them. The initial condition $m_0^*$ being smooth, we use the midpoint rule to approximate the integrals in (4.19).

Setting $\mathcal{G}_{\Delta x}(\bar{O}) := \mathcal{G}_{\Delta x} \cap \bar{O}$, Tables 1 and 2 below show the uniform and $L^2$ relative discrete errors

$$ (6.7) \qquad E_{\infty}(h^{\Delta}) = \frac{\max_{x_i \in \mathcal{G}_{\Delta x}(\Omega)} |h^{\Delta}(x_i) - h(x_i)|}{\max_{x_i \in \mathcal{G}_{\Delta x}(\Omega)} |h(x_i)|}, \quad E_2(h_{\Delta}) = \left( \frac{\sum_{x_i \in \mathcal{G}_{\Delta x}(\Omega)} |h^{\Delta}(x_i) - h(x_i)|^2}{\sum_{x_i \in \mathcal{G}_{\Delta x}(\Omega)} |h(x_i)|^2} \right)^{\frac{1}{2}}, $$

for $(h, h^{\Delta}) = (m^*(T, \cdot), m^{\Delta}(T, \cdot))$, $(m^*(T, \cdot), \bar{m}^{\Delta}(T, \cdot))$, $(v^*(0, \cdot), v^{\Delta}(0, \cdot))$, and $(v^*(0, \cdot), \bar{v}^{\Delta}(0, \cdot))$. Our results show smaller errors for the approximations computed with numerical quadrature, specially in the uniform norm. However, such precision is achieved at the expense of a high computational cost compared with the area-weighted approximation. Table 2 shows that the higher precision obtained by computing an approximation of $m^*$ by numerical quadrature does not significantly affect the approximation of the value function $v^*$.

<table><caption>TABLE 1. Errors for the approximation of $m^*(T, \cdot)$.</caption><thead><tr><th>Δx</th><th>E<sub>∞</sub>(m<sup>Δ</sup>(T, ·))</th><th>E<sub>∞</sub>(m<sup>Δ</sup>(T, ·))</th><th>E<sub>2</sub>(m<sup>Δ</sup>(T, ·))</th><th>E<sub>2</sub>(m<sup>Δ</sup>(T, ·))</th></tr></thead><tbody><tr><td>4.80 · 10<sup>-2</sup></td><td>8.41 · 10<sup>-3</sup></td><td>3.69 · 10<sup>-2</sup></td><td>6.30 · 10<sup>-3</sup></td><td>1.09 · 10<sup>-2</sup></td></tr><tr><td>2.40 · 10<sup>-2</sup></td><td>6.91 · 10<sup>-3</sup></td><td>3.25 · 10<sup>-2</sup></td><td>4.39 · 10<sup>-3</sup></td><td>1.05 · 10<sup>-2</sup></td></tr><tr><td>1.20 · 10<sup>-2</sup></td><td>3.94 · 10<sup>-3</sup></td><td>2.62 · 10<sup>-2</sup></td><td>2.77 · 10<sup>-3</sup></td><td>6.77 · 10<sup>-3</sup></td></tr><tr><td>6.00 · 10<sup>-3</sup></td><td>1.83 · 10<sup>-3</sup></td><td>2.44 · 10<sup>-2</sup></td><td>6.89 · 10<sup>-4</sup></td><td>2.67 · 10<sup>-3</sup></td></tr></tbody></table>

<table><caption>TABLE 2. Errors for the approximation of $v^*(0, \cdot)$.</caption><thead><tr><th>Δx</th><th>E<sub>∞</sub>(v<sup>Δ</sup>(0, ·))</th><th>E<sub>∞</sub>(v<sup>Δ</sup>(0, ·))</th><th>E<sub>2</sub>(v<sup>Δ</sup>(0, ·))</th><th>E<sub>2</sub>(v<sup>Δ</sup>(0, ·))</th></tr></thead><tbody><tr><td>4.80 · 10<sup>-2</sup></td><td>7.02 · 10<sup>-3</sup></td><td>7.11 · 10<sup>-3</sup></td><td>6.20 · 10<sup>-3</sup></td><td>6.31 · 10<sup>-3</sup></td></tr><tr><td>2.40 · 10<sup>-2</sup></td><td>5.74 · 10<sup>-3</sup></td><td>5.82 · 10<sup>-3</sup></td><td>4.90 · 10<sup>-3</sup></td><td>5.12 · 10<sup>-3</sup></td></tr><tr><td>1.20 · 10<sup>-2</sup></td><td>4.34 · 10<sup>-3</sup></td><td>4.37 · 10<sup>-3</sup></td><td>3.70 · 10<sup>-3</sup></td><td>3.75 · 10<sup>-3</sup></td></tr><tr><td>6.00 · 10<sup>-3</sup></td><td>3.30 · 10<sup>-3</sup></td><td>3.36 · 10<sup>-3</sup></td><td>2.95 · 10<sup>-3</sup></td><td>3.01 · 10<sup>-3</sup></td></tr></tbody></table>

### 6.3. A two-dimensional example.
In this test, we consider system (MFG) with $d = 2$, a quadratic Hamiltonian $H(x, p) = |p|^2/2$, and coupling terms having the form

$$ (6.8) \qquad F(x, \nu) = \gamma \min\{|x - \bar{x}|^2, R\} + (r_{\sigma} * \nu)(x) \quad \text{and} \quad G(x, \nu) = 0, $$

where $\gamma > 0$, $\bar{x} \in \mathbb{R}^2$, $R > 0$, and, for $\sigma > 0$, $r_{\sigma}(x) = e^{-|x|^2/2\sigma^2}/(2\pi\sigma^2)$ for all $x \in \mathbb{R}^2$. Given $\ell > 0$, $x_0^* \in [0, \ell]^2$, and $\sigma_0 > 0$, we consider the initial density

$$ (6.9) \qquad m_0^*(x) = \frac{\chi(x)}{\int_{[0,\ell]^2} \chi(y) dy} \quad \text{with} \quad \chi(x) = e^{-|x-x_0^*|^2/2\sigma_0^2} \mathbb{I}_{[0,\ell]^2}(x) \quad \text{for all } x \in \mathbb{R}^d. $$

Notice that the data above satisfy (H1), (H2), and (H3), with $p = \infty$. In our tests below, we choose $T = 1$, $\ell = 2$, $x_0^* = (0.75, 0.75)$, $\sigma_0 = 0.07$ in the initial condition, $\bar{x} = (1.75, 1.75)$, $R = 5$, $\sigma = 0.25$, and two values $\gamma = 0.5$ and $\gamma = 3$ in the running cost $F$. Since in this two-dimensional example the computational cost to solve (MFG$^{Δ}$) is important, in view of the discussion in Section 6.2 we implement the area-weighting method of Section 6.1 to approximate the integrals in (4.18). The integrals in (4.19), to approximate the initial condition $m_0^*$, are computed by using the midpoint rule. We set $\Delta x = 0.025$, $\Delta t = \Delta x^{2/3}$, and the mollifier in (3.17) is defined with $\mathbb{R}^2 \ni x \mapsto \rho(x) = e^{-|x|^2/2}/2\pi \in \mathbb{R}$ and $\varepsilon = \sqrt{\Delta t}/2$. Figure 1 shows the approximation $\bar{m}^\Delta$ of the exact distribution $m^*$ in the $x_1-x_2$ plane obtained after solving (MFG$^{Δ}$) for $\gamma = 0.5$ and $\gamma = 3$. On the left, we display the evolution of the initial distribution, concentrated around $x_0^*$, by overlaying the distributions $\bar{m}^\Delta(t_k, \cdot)$ for $k \in \mathcal{I}_{\Delta t}$. On the right, we display only the final distribution $\bar{m}^\Delta(T, \cdot)$. The simulation shows the effect of the positive constant $\gamma$, which weights the importance of reaching the target point $\bar{x}$. If $\gamma = 0.5$, the aversion to crowed regions, modeled
by the second term in the definition of $F$, has a more relevant impact on the distribution of the players than the term penalizing the distance to $\bar{x}$, while, if $\gamma = 3$, the latter term has a preponderant role in the evolution of the distribution of the agents.

<figure><img src="image_3.png" /><figcaption>(A) Time evolution $(\bar{m}^{\Delta}(t_k, \cdot))_{k \in \mathcal{I}_{\Delta t}}$ for $\gamma = 0.5$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(B) Final distribution $\bar{m}^{\Delta}(T, \cdot)$ for $\gamma = 0.5$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(C) Time evolution $(\bar{m}^{\Delta}(t_k, \cdot))_{k \in \mathcal{I}_{\Delta t}}$ for $\gamma = 3$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(D) Final distribution $\bar{m}^{\Delta}(T, \cdot)$ for $\gamma = 3$.</figcaption></figure>

FIGURE 1. Approximation of $m^*$ in both cases $\gamma = 0.5$ and $\gamma = 3$.

## APPENDIX A

*Proof of Proposition 2.1.* Let us fix $(t, x) \in [0, T \times \mathbb{R}^d$. The existence of $\alpha^{t,x} \in L^2([t, T]; \mathbb{R}^d)$, such that $v[\mu](t, x) = J^{t,x}[\mu](\alpha^{t,x})$, follows from (2.16), the continuity assumption on $F$ and $G$ in (H2), and the direct method in the calculus of variations. Setting $\alpha_0(s) = 0$ for all $s \in [t, T]$, the inequalities $J^{t,x}[\mu](\alpha^{t,x}) \le J^{t,x}[\mu](\alpha_0(s))$, (2.6), (2.10), (2.11), and (2.16), imply that

$$ (A.1) \qquad \int_t^T |\alpha^{t,x}(s)|^2 ds \le \tilde{C} := \frac{T(C_{L,2} + 2C_{F,1} + C_{L,7}) + 2C_{G,1}}{C_{L,6}}. $$

In particular, setting $\mathcal{A}^t = \{\alpha \in L^2([t, T]; \mathbb{R}^d), \int_t^T |\alpha(s)|^2 ds \le \tilde{C}\}$, we have

$$ (A.2) \qquad v[\mu](t, x) = \inf \{J^{t,x}[\mu](\alpha) \mid \alpha \in \mathcal{A}^t\}. $$
Thus, assertion (ii) follows from (2.6), (2.16), (A.1), (2.10), and (2.11). Moreover, it follows from conditions (2.7), (2.12), (2.13), and expression (A.2) that, for every $y \in \mathbb{R}^d$, we have

$$ egin{aligned} |v[\mu](t, x) - v[\mu](t, y)| &\le \sup_{\alpha \in \mathcal{A}^t} |J^{t,x}[\mu](\alpha) - J^{t,y}[\mu](\alpha)| \\ &\le \sup_{\alpha \in \mathcal{A}^t} \left\{ \int_t^T \left( C_{L,3}(1 + |\alpha(s)|^2) + C_{F,2} \right) |X^{t,x,\alpha}(s) - X^{t,y,\alpha}(s)|ds \right. \\ &\qquad \left. + C_{G,2} |X^{t,x,\alpha}(T) - X^{t,y,\alpha}(T)| \right\} \\ &\le (T(C_{L,3} + C_{F,2}) + C_{L,3}\tilde{C} + C_{G,2}) |x - y|, \end{aligned} $$

which shows (iii). Let us set $\bar{X} = X^{t,x,\alpha^{t,x}}$ and let $s \in [t, T]$. Since $v[\mu]$ satisfies the dynamic programming inequality

$$ v[\mu](s, \bar{X}(s)) \le \int_s^{s+h} \left( L(X^{s,\bar{X}(s),\alpha}(r), \alpha(r)) + F(X^{s,\bar{X}(s),\alpha}(r), \mu(r)) \right) dr + v[\mu](s+h, X^{s,\bar{X}(s),\alpha}(s+h)) $$

for all $h \in [0, T - s[$ and $\alpha \in L^2([t, T]; \mathbb{R}^d)$, by taking $\alpha = \alpha_0$, the equality

$$ (A.3) \qquad v[\mu](s, \bar{X}(s)) = J^{s,\bar{X}(s)}[\mu](\alpha^{t,x}|_{[s,T]}), $$

the estimates (2.16), (2.6), (2.10), the equality $X^{s,\bar{X}(s),\alpha_0}(s+h) = \bar{X}(s)$, and (2.25), imply that

$$ egin{aligned} C_{L,6} \int_s^{s+h} |\alpha^{t,x}(r)|^2 dr &\le h(C_{F,1} + C_{L,2}) + h(C_{L,7} + C_{F,1}) \\ &\quad + v[\mu](s+h, \bar{X}(s)) - v[\mu](s+h, \bar{X}(s+h)) \\ &\le h(2C_{F,1} + C_{L,2} + C_{L,7}) + C_{Lip} \int_s^{s+h} |\alpha^{t,x}(r)|dr. \end{aligned} $$

By Young's inequality, we get the existence of $C > 0$, independent of $(\mu, t, x)$, such that

$$ \int_s^{s+h} |\alpha^{t,x}(r)|^2 dr \le Ch $$

and, hence, by the Lebesgue differentiation theorem (see e.g. [14]), we have $\alpha^{t,x} \in L^\infty([0, T]; \mathbb{R}^d)$ and $\|\alpha^{t,x}\|_{L^\infty([0,T];\mathbb{R}^d)} \le \sqrt{C}$, which shows (i).

Finally, in order to show (iv), notice that, for every $y \in \mathbb{R}^d$, (2.9) implies that

$$ (A.4) \qquad L(x+y,a) - 2L(x,a) + L(x-y,a) \le C_{L,5}(1+|a|^2)|y|^2 \quad \text{for all } a \in \mathbb{R}^d. $$

Estimates (A.4), (2.14), (2.15), and (A.1), imply

$$ egin{aligned} v(t, x+y) + v(t, x-y) &\le \int_t^T \left( L(X^{t,x+y,\alpha^{t,x}}(s), \alpha^{t,x}(s)) + L(X^{t,x-y,\alpha^{t,x}}(s), \alpha^{t,x}(s)) \right. \\ &\qquad \left. + F(X^{t,x+y,\alpha^{t,x}}(s), \mu(s)) + F(X^{t,x-y,\alpha^{t,x}}(s), \mu(s)) \right) ds \\ &\quad + G(X^{t,x+y,\alpha^{t,x}}(T), \mu(T)) + G(X^{t,x-y,\alpha^{t,x}}(T), \mu(T)) \\ (A.5) \qquad &\le 2 \int_t^T \left( L(X^{t,x,\alpha^{t,x}}(s), \alpha^{t,x}(s)) + F(X^{t,x,\alpha^{t,x}}(s), \mu(s)) \right) ds \\ &\quad + 2G(X^{t,x,\alpha^{t,x}}(T), \mu(T)) \\ &\quad + \int_t^T \left( C_{L,5}(1 + |\alpha^{t,x}(s)|^2) + C_{F,3} \right) |y|^2 ds + C_{G,3} |y|^2 \\ &\le 2v[\mu](t, x) + (T(C_{L,5} + C_{F,3}) + C_{L,5}\tilde{C} + C_{G,3}) |y|^2, \end{aligned} $$

from which (2.26) follows. $\square$

<figure><img src="image_7.png" /></figure>
REFERENCES

[1] Y. Achdou, F. Camilli, and I. Capuzzo-Dolcetta. Mean field games: convergence of a finite difference method. *SIAM J. Numer. Anal.*, 51(5):2585-2612, 2013.

[2] Y. Achdou, F. Camilli, and L. Corrias. On numerical approximation of the Hamilton-Jacobi-transport system arising in high frequency approximations. *Discrete Contin. Dyn. Syst. Ser. B*, 19(3):629-650, 2014.

[3] Y. Achdou and I. Capuzzo-Dolcetta. Mean field games: numerical methods. *SIAM J. Numer. Anal.*, 48(3):1136-1162, 2010.

[4] Y. Achdou, P. Cardaliaguet, F. Delarue, A. Porretta, and F. Santambrogio. *Mean field games*, volume 2281 of *Lecture Notes in Mathematics*. Springer, Cham; Centro Internazionale Matematico Estivo (C.I.M.E.), Florence, 2020.

[5] Y. Achdou and M. Laurière. Mean field games and applications: numerical aspects. In *Mean field games*, volume 2281 of *Lecture Notes in Math.*, pages 249-307. Springer, Cham, 2020.

[6] Y. Achdou and A. Porretta. Convergence of a finite difference scheme to weak solutions of the system of partial differential equations arising in mean field games. *SIAM J. Numer. Anal.*, 54(1):161-186, 2016.

[7] N. Almulla, R. Ferreira, and D. Gomes. Two numerical approaches to stationary mean-field games. *Dyn. Games Appl.*, 7(4):657-682, 2017.

[8] L. Ambrosio. Transport equation and Cauchy problem for BV vector fields. *Invent. Math.*, 158(2):227-260, 2004.

[9] L. Ambrosio, N. Gigli, and G. Savaré. *Gradient flows in metric spaces and in the space of probability measures*. Second edition. Lecture notes in Mathematics ETH Zürich. Birkhäuser Verlag, Bassel, 2008.

[10] A. Angiuli, J.-P. Fouque, and M. Laurière. Unified reinforcement Q-learning for mean field game and control problems. *Math. Control Signals Systems*, 34(2):217-271, 2022.

[11] J.-P. Aubin and A. Cellina. *Differential inclusions*, volume 264 of *Grundlehren der mathematischen Wissenschaften*. Springer-Verlag, Berlin, 1984.

[12] M. Bardi and I. Capuzzo Dolcetta. *Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations*. Birkauser, 1996.

[13] G. Barles and P. E. Souganidis. Convergence of approximation schemes for fully nonlinear second order equations. *Asymptotic Anal.*, 4(3):271-283, 1991.

[14] V. I. Bogachev. *Measure theory. Vol. I, II*. Springer-Verlag, Berlin, 2007.

[15] J.F. Bonnans and A. Shapiro. *Perturbation analysis of optimization problems*. Springer-Verlag, New York, 2000.

[16] E. Calzola, E. Carlini, and F. J. Silva. A high-order scheme for mean field games. *arXiv:2207.08463*, 2022.

[17] F. Camilli and F. J. Silva. A semi-discrete in time approximation for a first order-finite mean field game problem. *Network and Heterogeneous Media*, 7-2:263-277, 2012.

[18] P. Cannarsa and C. Sinestrari. *Semiconcave Functions, Hamilton-Jacobi Equations, and Optimal Control*. Progress in Nonlinear Differential Equations and Their Applications. Birkauser, 2004.

[19] P. Cardaliaguet and A. Porretta. An introduction to mean field game theory. In *Mean field games*, volume 2281 of *Lecture Notes in Math.*, pages 1-158. Springer, Cham, 2020.

[20] E. Carlini and F. J. Silva. A fully discrete semi-Lagrangian scheme for a first order mean field game problem. *SIAM J. Numer. Anal.*, 52(1):45-67, 2014.

[21] E. Carlini and F. J. Silva. A semi-Lagrangian scheme for a degenerate second order mean field game system. *Discrete and Continuous Dynamical Systems*, 35(9):4269-4292, 2015.

[22] E. Carlini and F. J. Silva. On the discretization of some nonlinear Fokker-Planck-Kolmogorov equations and applications. *SIAM J. Numer. Anal.*, 56(4):2148-2177, 2018.

[23] R. Carmona and F. Delarue. *Probabilistic theory of mean field games with applications. I*, volume 83 of *Probability Theory and Stochastic Modelling*. Springer, Cham, 2018.

[24] R. Carmona and F. Delarue. *Probabilistic theory of mean field games with applications. II*, volume 84 of *Probability Theory and Stochastic Modelling*. Springer, Cham, 2018.

[25] R. Carmona and M. Laurière. Convergence analysis of machine learning algorithms for the numerical solution of mean field control and games I: The ergodic case. *SIAM J. Numer. Anal.*, 59(3):1455-1485, 2021.

[26] R. Carmona and M. Laurière. Convergence analysis of machine learning algorithms for the numerical solution of mean field control and games: II—The finite horizon case. *Ann. Appl. Probab.*, 32(6):4065-4105, 2022.

[27] I. Chowdhury, O. Ersland, and E.R. Jakobsen. On numerical approximations of fractional and nonlocal mean field games. *Found Comput Math*, 2022.

[28] F. Da Lio and O. Ley. Convex Hamilton-Jacobi equations under superlinear growth conditions on data. *Appl. Math. Optim.*, 63(3):309-339, 2011.

[29] M. Falcone and R. Ferretti. *Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations*. MOS-SIAM Series on Optimization, 2013.

[30] M. Fischer and F. J. Silva. On the asymptotic nature of first order mean field games. *Appl. Math. Optim.*, 84(2):2327-2357, 2021.

[31] J. Gianatti and F. J. Silva. Approximation of deterministic mean field games with control-affine dynamics. *Preprint*, 2022.

[32] D. A. Gomes, J. Mohr, and R. Souza. Discrete time, finite state space mean field games. *Journal de Mathématiques Pures et Appliquées*, 93:308-328, 2010.

[33] D. A. Gomes, E. A. Pimentel, and V. Voskanyan. *Regularity theory for mean-field game systems*. SpringerBriefs in Mathematics. Springer, Cham, 2016.

[34] D. A. Gomes and J. Saúde. Mean field games models-A brief survey. *Dyn. Games Appl.*, 4(2):110-154, 2014.
[35] O. Guéant. Mean field games equations with quadratic Hamiltonian: a specific approach. *Math. Models Methods Appl. Sci.*, 22(9):1250022, 37, 2012.

[36] S. Hadikhanloo and F. J. Silva. Finite mean field games: fictitious play and convergence to a first order continuous mean field game. *J. Math. Pures Appl. (9)*, 132:369–397, 2019.

[37] M. Huang, R. P. Malhamé, and P. E. Caines. Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle. *Commun. Inf. Syst.*, 6(3):221–251, 2006.

[38] J.-M. Lasry and P.-L. Lions. Jeux à champ moyen I. Le cas stationnaire. *C. R. Math. Acad. Sci. Paris*, 343:619–625, 2006.

[39] J.-M. Lasry and P.-L. Lions. Jeux à champ moyen II. Horizon fini et contrôle optimal. *C. R. Math. Acad. Sci. Paris*, 343:679–684, 2006.

[40] J.-M. Lasry and P.-L. Lions. Mean field games. *Jpn. J. Math.*, 2:229–260, 2007.

[41] M. Laurière. Numerical methods for mean field games and mean field type control. In *Mean field games*, volume 78 of Proc. Sympos. Appl. Math., pages 221–282. Amer. Math. Soc., Providence, RI, 2021.

[42] S. Liu, M. Jacobs, W. Li, L. Nurbekyan, and S. J. Osher. Computational methods for first-order nonlocal mean field games with applications. *SIAM J. Numer. Anal.*, 59(5):2639–2668, 2021.

[43] K. W. Morton, A. Priestley, and E. Süli. Stability of the Lagrange-Galerkin method with non-exact integration. *RAIRO, Modélisation Math. Anal. Numér.*, 22(4):625–653, 1988.

[44] L. Nurbekyan and J. Saúde. Fourier approximation methods for first-order nonlocal mean-field games. *Port. Math.*, 75(3-4):367–396, 2018.

[45] B. Piccoli and A. Tosin. Time-evolving measures and macroscopic modeling of pedestrian flow. *Arch. Ration. Mech. Anal.*, 199(3):707–738, 2011.

[46] A. Quarteroni and A. Valli. *Numerical approximation of partial differential equations*. Springer Verlag, 1994.

[47] R. T. Rockafellar and R. J.-B. Wets. *Variational analysis*, volume 317 of *Grundlehren der mathematischen Wissenschaften*. Springer-Verlag, Berlin, 1998.

[48] A. Tosin and P. Frasca. Existence and approximation of probability measure solutions to models of collective behaviors. *Netw. Heterog. Media*, 6(3):561–596, 2011.