\section{Introduction}

The approximation of the equilibria is at the heart of the game theory. The classic literature introduces a natural relation connecting the equilibria of the games and the optima of the sequential decision problems.  This leads to a fruitful research on the topics such as approachability, regret and calibration, see the survey by V. Perchet \cite{perchet14} and the books by N. Cesa-Bianchi and G. Lugosi \cite{CL06} and by D. Fudenberg and D. K. Levine \cite{FL98}. In particular, the gradient-based strategy often plays a crucial rule in approximating the equilibria. In the present paper we study the analog to the gradient-based strategy in the continuous-time setting, namely, we aim at approximating the equilibria of games using the diffusion processes encoded with the gradients of potential functions.

Consider a game with $n$ players. The mixed Nash equilibrium is defined to be a collection of probability measures $(\nu^{*,i})_{i=1,\cdots, n}$ such that
\begin{equation}\label{nashclassic}
\nu^{*, i} \in \argmin_{\nu^i} \int f^i(x^1 ,\cdots, x^n) \nu^i(dx^i) \prod_{j\neq i}\nu^{*,j}(dx^j),
\end{equation}
which means that each player can no longer improve his performance  by making a unilateral change of strategy. Note that in this classical setting, the potential function of each player
\begin{equation*}
\nu^i\mapsto  F^i\big(\nu^i, (\nu^j)_{j\neq i}\big):=\int f^i(x^1 ,\cdots, x^n) \nu^i(dx^i) \prod_{j\neq i}\nu^{j}(dx^j)
\end{equation*}
is linear. In this paper we shall allow the potential function to be {\it nonlinear} in view of the applications, in particular, to the neural networks (see \cref{sec:app}).
As another generalization to  the classic theory, we consider games on a {\it random environment}. Introduce a space of environment $\mathbb{Y}$ and fix a probability measure ${\mathsf m}$ on it. We urge each player to choose a strategy among the probability measures $\nu^i$ on the product space $\mathbb{R}^{n^i}\times \mathbb{Y}$ such that the marginal law of $\nu^i$ on $\mathbb{Y}$, $\nu^i_Y$, matches the fixed distribution ${\mathsf m}$. Typically, in our framework we  consider the game of which the Nash equilibrium is  a collection of probability measures $(\nu^{*,i})_{i=1,\cdots, n}$ on the product spaces such that
\begin{equation}\label{nashcor}
\nu^{*, i} \in \argmin_{\nu^i : \nu^i_Y={\mathsf m}} \int f^i(x^1 ,\cdots, x^n, {\mathsf y}) \nu^i(dx^i|{\mathsf y}) \prod_{j\neq i}\nu^{*,j}(dx^j|{\mathsf y}) {\mathsf m}(d{\mathsf y}) +\frac{\sigma^2 }{2} H(\nu^i | {\rm Leb}\times {\mathsf m} ),
\end{equation}
where $\nu(\cdot|{\mathsf y})$ denotes the conditional probability given ${\mathsf y}$, and we add the relative entropy $H$ as a regularizer.
In contrast to the conventional definition of the Nash equilibrium \cref{nashclassic}, where the players' strategies are uncorrelated, in our setting the strategies of the players are allowed to be coupled through the environment. Moreover,  the general framework of the present paper goes beyond the particular game \cref{nashcor}, by allowing the cost function to be nonlinear in $(\nu^i)_{i=1,\cdots, n}$. As an application,  we observe (\cref{eg:dynamicgame}) that relaxed dynamic games can be viewed as games on random environment, where the environment $\mathbb{Y}$ is the time horizon.

One of our main contributions is the first order condition of the optimization on the probability space given a marginal constraint (\cref{thm:FOC}), which naturally provides a necessary condition for being a Nash equilibrium of a game on  random environment (\cref{cor:NC_Nash}). This result is a generalization to the first order condition in Proposition 2.4 in \cite{HRSS19} for the optimization on the probability space without marginal constraint. The key ingredient for this analysis is the linear functional derivative $\frac{\delta F}{\delta \nu}$, first introduced for the variational calculus and recently popularized by the study on the mean-field games, see e.g. Cardaliaguet et al. \cite{cardaliaguet2015master}, Delarue et al. \cite{DLR19, DLR19II}, Chassagneux et al. \cite{CLTse19}. Roughly speaking, we prove that
\begin{multline}\label{intro:FOC}
\mbox{if $\nu^*\in \argmin_{\nu: \nu_Y ={\mathsf m}} F(\nu) + \frac{\sigma^2}{2}H(\nu| {\rm Leb}\times {\mathsf m})$,}\\
\mbox{then $\nabla_x  \frac{\delta F}{\delta \nu}(\nu^*, x, {\mathsf y}) + \frac{\sigma^2 }{2} \nabla_x \ln \nu^*(x|{\mathsf y}) = 0$, for all $x$, ${\mathsf m}$-a.s. ${\mathsf y}$. }
\end{multline}
Besides the first order condition for the Nash equilibrium, we also provide sufficient conditions on the linear functional derivative so that the game on a random environment admits a (unique) equilibrium.

The first order equation in \cref{intro:FOC}  clearly links the minimizer $\nu^*$ (or the Nash equilibrium in the context of games) to the invariant measure  of a system of  diffusion processes, see \cref{eq:MFL} below. Since the dynamics of the diffusion processes depends on their marginal distributions (in other word, McKean-Vlasov diffusion, see \cite{Meleard96, Sznit89}) and involves the gradients of the potential functions, we name the system mean-field Langevin  (MFL) system. Further, we study the different settings where the marginal laws of the MFL system converge to the unique invariant measure, which, due to the first order condition, must coincide with the Nash equilibrium of the game on random environment.
In the case with small dependence on the marginal laws, we prove the exponential ergodicity of the MFL system, using the reflection coupling (see \cite{eberle11, eberle2019quantitative}).
In the  case of one-player game (in other word, optimization), once the potential function is convex, we use an argument, similar to that in \cite{HRSS19, HKR19}, based on the Lasalle's invariant principle to prove the (non-exponential) ergodicity of the MFL system under  mild conditions on the coefficients.

In view of applications, our result can be used to justify the applicability of the gradient descent algorithm for training (deep) neural networks. As mentioned in \cite{JSS19, mei2018mean, mei2019mean, HRSS19, HKR19}, the supervised learning with (deep) neural networks can be viewed as a minimization problem (or optimal control problem in the context of deep learning) on the space of probability measures, and the gradient descent algorithm is approximately a discretization of the corresponding mean-field Langevin dynamics. The present paper provides a more general framework for such studies. In particular, it is remarkable that in \cref{sec:app} we provide a theoretically convergent numerical scheme for the generative adversarial networks (GAN) as well as a way to characterize the training error.

The rest of the paper is organized as follows. In \cref{sec:notation} we introduce the definitions of a game on a random environment and the corresponding MFL system.  In \cref{sec:main}  the main theorems of the paper are stated without proofs. In \cref{sec:app} we present the applications to the dynamic games and the GAN. Finally in \cref{sec:proof} we present the proofs of the main theorems.

\section{Notation and definitions}\label{sec:notation}

\subsection{Preliminary}

Denote by $\mathbb{Y}$ the space of environment, and assume $\mathbb{Y}$ to be Polish. Throughout the paper, we fix a probability measure ${\mathsf m}$ on $\mathbb{Y}$. Define the product space $\bar \mathbb{R}^d:= \mathbb{R}^d\times \mathbb{Y}$ for $d\in \mathbb{N}$.  In this paper we consider a game in which the players choose strategies among $\Pi:=\{\bar\pi\in {\cal P}_p(\bar\mathbb{R}^d): \bar\pi(\mathbb{R}^d,\cdot) = {\mathsf m} \}$, where ${\cal P}_p(\bar\mathbb{R}^d)$ is the space of probability measures on $\bar\mathbb{R}^d$ with finite $p$-moments for some $p\ge 1$.  We say that a function $F:\Pi\rightarrow \mathbb{R}$ is in ${\cal C}^1$ if there exists a  function $\frac{\delta F}{\delta \nu}: \Pi\times \bar\mathbb{R}^d \rightarrow \mathbb{R}$ such that for all $\nu,\nu'\in \Pi$
\begin{equation}\label{C1derv}
F(\nu') - F(\nu) = \int_0^1 \int_{\bar\mathbb{R}^d} \frac{\delta F}{\delta \nu}\big((1-\lambda)\nu+\lambda \nu',  \bar x\big)(\nu'-\nu)(d \bar x)d\lambda.
\end{equation}
We will refer to $\frac{\delta F}{\delta \nu}$ as the linear functional derivative.
There is at most one $\frac{\delta F}{\delta \nu}$, modulo a constant, satisfying \cref{C1derv}.

 Here is the basic assumption we apply throughout the paper.

\begin{assum}[basic assumption]\label{assum:base}
Assume that for some $p\ge 1$, the function $F: \Pi\rightarrow \mathbb{R}$ belongs to ${\cal C}^1$ and
\begin{itemize}
\item $F$ is ${\cal W}_p$-continuous, where ${\cal W}_p$ stands for the $p$-Wasserstein distance;

\item $\frac{\delta F}{\delta \bar\pi}: (\bar\pi, x, {\mathsf y})\in\Pi\times\mathbb{R}^d\times\mathbb{Y}\mapsto \frac{\delta F}{\delta \bar\pi}(\bar\pi, x, {\mathsf y}) \in \mathbb{R}$ is ${\cal W}_p$-continuous in $\bar\pi$ and  continuously differentiable in $x$;

\item $\frac{\delta F}{\delta \bar\pi}$ is of $p$-polynomial growth in $\bar x= (x, {\mathsf y})$, that is, $\sup_{\bar\pi\in \Pi}\big| \frac{\delta F}{\delta \bar\pi}(\bar\pi, \bar x) \big| \le C(1+|\bar x|^p)$.
\end{itemize}
\end{assum}

\begin{remark}
 Since in our setting the law on the environment $\mathbb{Y}$ is fixed, by disintegration we may identify a distribution
$\bar \pi\in \Pi$
with the probability measures $\big(\pi(\cdot | {\mathsf y})\big)_{{\mathsf y}\in\mathbb{Y}}\subset {\cal P}_p(\mathbb{R}^d)$ such that
$\bar \pi(d\bar x) = \pi(dx | {\mathsf y}) {\mathsf m}(d{\mathsf y})$.
\end{remark}

\subsection{Game on random environment}

In this paper, we consider a particular game in which the strategies of the $n$ players are correlated through the random environment (or signal) $\mathbb{Y}$.
 Let $n^i\in \mathbb{N}$ for $i=1,\cdots, n$ and $N:=\sum_{i=1}^n n^i$.
 As mentioned before, the $i$-th player chooses his strategy (a probability measure) among $\Pi^i:=\{ \nu\in {\cal P}_p(\bar\mathbb{R}^{n^i}) : \nu(\mathbb{R}^{n^i} ,\cdot) = {\mathsf m} \}$, while the joint distribution of the other players' strategies belongs to the space $\Pi^{-i}:=\{ \nu\in {\cal P}_p(\bar\mathbb{R}^{N-n^i}) : \nu(\mathbb{R}^{N-n^i} ,\cdot) = {\mathsf m} \}$.
 The $i$-th player aims at optimizing his objective function $F^i: \Pi^i\times\Pi^{-i}\rightarrow \mathbb{R}$. More precisely, he faces the optimization:
\begin{equation*}
\mbox{given $\mu\in \Pi^{-i}$},\quad\q\quad
\mbox{solve}\quad\inf_{\nu\in \Pi^i} F^i(\nu, \mu).
\end{equation*}
In this paper, we are more interested in solving a regularized version of the game above. We use the relative entropy with respect to ${\rm Leb}^{n^i}\times{\mathsf m}$, denoted by $H^i$, as the regularizer. Namely, given $\mu\in \Pi^{-i}$ the $i$-th player solves:
\begin{equation*}
\quad\inf_{\nu\in \Pi^i} V^i(\nu, \mu), \quad V^i(\nu, \mu):= F^i(\nu, \mu) +\frac{\sigma^2}{2} H^i(\nu)\quad\mbox{for some $\sigma>0$.}
\end{equation*}
For $\bar\pi\in \Pi$, we denote by $\bar\pi^i\in \Pi^i$ its marginal distribution on $\bar\mathbb{R}^{n^i}$, and by $\bar\pi^{-i}\in \Pi^{-i}$ the marginal distribution on  $\bar\mathbb{R}^{N-n^i}$.

\begin{definition}\label{defn:nash}
A probability measure $\bar\pi \in \Pi$ is a Nash equilibrium of this game, if
\begin{equation*}
\bar\pi^i \in \arg\min_{\nu\in \Pi^i} V^i(\nu, \bar\pi^{-i}), \quad\q\mbox{for all}\quad i=1,\cdots, n.
\end{equation*}
\end{definition}

\begin{eg}\label{eg:dynamicgame}
To have a concrete example of games on random environment, we refer to the dynamic games, both discrete-time and continuous-time models. In the discrete-time case, let $\mathbb{Y}:=\{1, \cdots, T\}$ for some $T\in \mathbb{N}$ and  ${\mathsf m}$ be the uniform distribution on $\mathbb{Y}$. Define the controlled dynamics:
\begin{equation*}
\Theta_{\mathsf y}^i =   \varphi^i_{\mathsf y} (\Theta_{{\mathsf y}-1}^i, \pi^i(\cdot| {\mathsf y}), \bar \pi^{-i}), \quad\mbox{where}\quad \bar\pi^i(\cdot, {\mathsf y}) =  \pi^i(\cdot|{\mathsf y}) {\mathsf m}({\mathsf y}), \quad\mbox{for}\quad {\mathsf y}\in \mathbb{Y}.
\end{equation*}
If the $n$ players minimize the objective functions of the form $f^i\big((\Theta^i_{\mathsf y})_{{\mathsf y}\in \mathbb{Y}}\big)$ by choosing the strategy $\bar\pi^i$, then the game fits the framework of this paper.

Similarly for the continuous-time model, consider the space $\mathbb{Y}:=[0,T]$ for $T\in \mathbb{R}$ and let ${\mathsf m}$ be the uniform distribution on the interval. Define the continuous-time dynamics:
\begin{equation*}
d\Theta^i_{\mathsf y} = \varphi^i \big( \pi^i(\cdot|{\mathsf y}), \bar \pi^{-i}, \Theta_{\mathsf y}^i, {\mathsf y}\big)d{\mathsf y}, \quad\mbox{where}\quad \bar\pi^i (\cdot, d {\mathsf y}) = \pi^i(\cdot|{\mathsf y}) {\mathsf m}(d{\mathsf y}), \quad\mbox{for}\quad {\mathsf y}\in \mathbb{Y}.
\end{equation*}
If the $n$ players  minimize the objective functions of the form $f^i\big((\Theta^i_{\mathsf y})_{{\mathsf y}\in \mathbb{Y}}\big)$ by choosing the strategy $\bar\pi^i$, then this game also fits in the  framework discussed above.
\end{eg}

\subsection{Mean-field Langevin system}

For any fixed $\mu\in \Pi^{-i}$, we assume that $F^i(\cdot, \mu): \nu\in \Pi^i \mapsto F^i(\nu,\mu)\in  \mathbb{R}$ satisfies \cref{assum:base}. The linear derivative is denoted by $ \frac{\delta F^i}{\delta \nu}(\cdot, \mu, \cdot):  (\nu, \bar x^i) \mapsto  \frac{\delta F^i}{\delta \nu} (\nu , \mu, \bar x^i)$, with $\bar x^i= (x^i, {\mathsf y})\in \bar\mathbb{R}^{n^i}$.  In order to compute Nash equilibria of the game on random environment,  we are interested in the following mean-field Langevin (MFL) dynamics:
\begin{equation}\label{eq:MFL}
d X^i_t = - \nabla_{x^i} \frac{\delta F^i}{\delta \nu} (\bar\pi^i_t , \bar\pi^{-i}_t, X^i_t, Y)dt + \sigma dW^i_t, \quad\mbox{for}\quad i=1,\cdots,n,
\end{equation}
where $W=(W^i)_i$ is an $N$-dimension Brownian motion, $Y$ is a random variable taking values in $\mathbb{Y}$ and satisfying the law ${\mathsf m}$, and $\bar\pi_t := {\rm Law}\bar X_t$ with $\bar X_t :=(X^1_t, \cdots, X^n_t, Y)$.
 In this paper we will discuss the relation between the MFL dynamics and the Nash equilibrium of the game on the random environment.

\begin{remark}
Here are some important observations:
\begin{itemize}

\item The random variable $Y$ plays the role of parameter in the MFL system. This leads us to study the system:
\begin{equation}\label{eq:MFL_sys}
d X^{\mathsf y}_t = - \left(\nabla_{x^i} \frac{\delta F^i}{\delta \nu} (\bar\pi^i_t , \bar\pi^{-i}_t, X^{{\mathsf y},i}_t, {\mathsf y})\right)_{i=1,\cdots,n}dt
		 + \sigma dW_t, \quad\mbox{for}\quad  {\mathsf m}\mbox{-a.s.}~{\mathsf y}\in \mathbb{Y}.
\end{equation}
Formally, the marginal laws of the MFL system above with a fixed $y\in \mathbb{Y}$ satisfy the following system of Fokker-Planck equations:
\begin{multline}
\label{eq:FP}
\partial_t \pi^i(\cdot|{\mathsf y}) = \nabla_{x^i}\cdot \left(\nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\bar\pi^i, \bar\pi^{-i}, \cdot, {\mathsf y})\pi^i(\cdot|{\mathsf y})   +\frac{\sigma^2}{2} \nabla_{x^i} \pi^i(\cdot|{\mathsf y})\right), \quad \\
\mbox{for all $i=1, \cdots,n$, ${\mathsf m}$-a.s. ${\mathsf y}\in \mathbb{Y}$.}
\end{multline}

\item For fixed ${\mathsf y}\in \mathbb{Y}$, the dynamic systems for $\big(X^i({\mathsf y})\big)_i$ are only weakly coupled through the marginal distributions.

\item Although we name the system after Langevin, the drift term of the dynamics of the aggregated vector $ (X^i)_{i=1,\cdots, n}$ is in general not in the form of the gradient of a potential function.
\end{itemize}
\end{remark}

\section{Main results}\label{sec:main}

\subsection{Optimization with marginal constraint}

One of our observations  is the following first order condition of the optimization over the probability measures with marginal constraint.

\begin{theorem}[first order condition]\label{thm:FOC}
 Let $F: \Pi\rightarrow \mathbb{R}$ satisfy \cref{assum:base}. Define $V(\bar\pi):=F(\bar\pi) + \eta H(\bar\pi)$ for some $\eta>0$.  If $\bar\pi^*\in \argmin_{\bar\pi\in \Pi} V(\bar\pi)$, then
\begin{equation}\label{eq:necessary}
\nabla_x \frac{\delta F}{\delta \bar\pi}(\bar\pi^*, \cdot,{\mathsf y}) + \eta \nabla_x\ln\big(\pi^*(\cdot|{\mathsf y})\big) =0\quad\mbox{for ${\mathsf m}$-a.s. ${\mathsf y}$.}
\end{equation}
Conversely, if we additionally assume that $F$ is convex, then $\bar\pi^*\in \Pi$ satisfying \cref{eq:necessary} implies
$\bar\pi^*\in \argmin_{\bar\pi\in \Pi} V(\bar\pi)$.
\end{theorem}

\begin{remark}
We remark that
\begin{itemize}
\item  the regularizer $H(\bar\pi)$ plays an important role for the proof of the necessary condition. Without it, for  $\bar\pi^*\in \argmin_{\bar\pi\in \Pi} V(\bar\pi)$  we can only conclude that there is a measurable function $f:\mathbb{Y}\rightarrow\mathbb{R}$ such that
\begin{equation*}
\frac{\delta F}{\delta \bar\pi}(\bar\pi^*, x,{\mathsf y}) =f({\mathsf y}) ,\quad \bar\pi^*\mbox{-a.s.};
\end{equation*}

\item for the readers more interested in the minimization of the unregularized potential function $F$, by standard argument (see e.g. \cite[Proposition 2.3]{HRSS19}) one may prove that under mild conditions the minimum of $F+\eta H$ converges to the minimum of $F$ as $\eta\rightarrow 0$.

\end{itemize}
\end{remark}

\subsection{Equilibria of games on random environment}

 Applying the first order condition above to the context of the game on  random environment, we immediately obtain the following necessary condition for the Nash equilibria.

\begin{corollary}[Necessary condition for Nash equilibria]
\label{cor:NC_Nash}
For $i=1,\cdots, n$ and $\mu\in \Pi^{-i}$, let $F^i(\cdot, \mu): \nu \in\Pi^i\mapsto F^i(\nu, \mu)\in \mathbb{R}$  satisfy \cref{assum:base}. If  $\bar\pi\in \Pi$ is a Nash equilibrium, we have for  $i=1,\cdots,n,$
\begin{equation}\label{eq:FOCeq}
\nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\bar\pi^i, \bar\pi^{-i}, x^i,{\mathsf y}) + \frac{\sigma^2}{2} \nabla_{x^i}\ln\big(\pi^i(x^i|{\mathsf y})\big) =0\quad\mbox{for all $x^i\in \mathbb{R}^{n^i}$ and ${\mathsf m}$-a.s. ${\mathsf y}\in \mathbb{Y}$}.
\end{equation}
\end{corollary}

We shall use  the first order equation \cref{eq:FOCeq} to show the following sufficient condition for the uniqueness of Nash equilibrium.

\begin{corollary}[Uniqueness of Nash equilibrium: Monotonicity]
\label{cor:uniq_monot}
The functions $(F^i)_{i=1,\cdots, n}$ satisfy the monotonicity condition, if for $\bar\pi, \bar\pi'\in \Pi$ we have
\begin{equation*}
\label{eq:monotonicity}
\sum_{i=1}^n \int \left( \frac{\delta F^i}{\delta \nu}(\bar\pi^i, \bar\pi^{-i}, \bar x^i) - \frac{\delta F^i}{\delta \nu}(\bar\pi'^i, \bar\pi'^{-i}, \bar x^i)  \right)(\bar\pi - \bar\pi')(d\bar x) \ge 0.
\end{equation*}
We have the following results:
\begin{itemize}
\item[{\rm (i)}]  for $n=1$,  if a function $F$ satisfies the monotonicity condition then it is convex on $\Pi$. Conversely, if $F$ is convex and  satisfies \cref{assum:base}, then $F$ satisfies the monotonicity condition.

\item[{\rm (ii)}] in general ($n\ge 1$), for $i=1,\cdots, n$  and any $\mu\in\Pi^{-i}$, let $F^i(\cdot, \mu): \nu \in\Pi^i \mapsto F^i(\nu, \mu)\in \mathbb{R}$  satisfy \cref{assum:base} and $ (F^i)_{i=1, \cdots, n} $ satisfy the monotonicity condition. Then for any two Nash equilibria $ \bar \pi^* , \bar\pi'^* \in \Pi $ we have $(\bar\pi^*)^i =(\bar\pi'^*)^i$ for all $i=1, \cdots,n$.
\end{itemize}
\end{corollary}

\begin{remark}
Similar monotonicity conditions are common assumptions to ensure the uniqueness of equilibrium in the game theory, in particular in the literature of mean-field games, see e.g. \cite{lasry2007}.
\end{remark}

As for the existence of Nash equilibria, we obtain the following result following the classical argument based on the fixed point theorem.

\begin{theorem}[Existence of equilibria]\label{thm:exist_nash}
Assume  that  for $i=1,\cdots, n$,  and $\mu\in \Pi^{-i}$
\begin{itemize}
\item[{\rm (i)}]  the set $\argmin_{\nu\in \Pi^i} V^i(\nu, \mu) $ is  non-empty and convex;

\item[{\rm (ii)}]  the function $\tilde F^i (\bar\pi):= F^i(\bar\pi^i, \bar\pi^{-i})$ is  ${\cal W}_p$-continuous on $\Pi$;

\item[{\rm (iii)}]  the function $F^i(\cdot, \mu): \nu \in\Pi^i \mapsto F^i(\nu, \mu)\in \mathbb{R}$  satisfies \cref{assum:base},  and there exist some $q\ge q'>0,~ C,C'>0 \in \mathbb{R}$ such that for all $\bar\pi\in \Pi$ we have
\begin{eqnarray}\label{basic_dissipative}
C'| \bar x^i |^{q'} - C\le \frac{\delta F^i}{\delta \nu}(\bar\pi^i, \bar\pi^{-i},  \bar x^i) \le C| \bar x^i |^q + C.
\end{eqnarray}
\end{itemize}
Then there exists at least one Nash equilibrium $\bar\pi^*\in \Pi$ for the game on random environment.
\end{theorem}

\begin{remark}
There are various sufficient conditions so that  the set $\argmin_{\nu\in \Pi^i} V^i(\nu, \mu)$ is convex, for example, the function $\nu\mapsto V^i(\nu, \mu)$ is quasi-convex, or  $V^i(\nu, \mu)$ has a unique minimizer. That is why we leave the assumption (i) in the abstract form.
\end{remark}

\subsection{Invariant measure of the MFL system}

In view of the Fokker-Planck equation \cref{eq:FP}, the first order equation \cref{eq:FOCeq} appears to be a sufficient condition for $ \bar \pi$ being an invariant measure of the MFL system \cref{eq:MFL}. That is why we consider the MFL dynamics as a reasonable tool to compute the Nash equilibria of the game on random environment.

The following \cref{thm:well} suggests that proving the {\it existence} of Nash equilibria  and the {\it uniqueness} of the invariant measure, we can establish the equivalence between the invariant measure of  \cref{eq:MFL} and one Nash equilibrium.
While the existence of Nash equilibria has been discussed in\cref{thm:exist_nash},
the uniqueness of invariant measure of mean-field dynamics is more complicated and is indeed a long-standing problem in probability and analysis. We are going to use the coupling argument in order to obtain the contraction result in \cref{thm:contraction}.

Define the average Wasserstein distance:
\begin{equation*}
\overline{\cal W}_p(\bar\pi, \bar\pi') : =\Big( \int_\mathbb{Y} {\cal W}_p^p\big(\pi(\cdot|{\mathsf y}), \pi'(\cdot|{\mathsf y})\big) {\mathsf m}(d{\mathsf y}) \Big)^\frac1p,
\end{equation*}
and  the spaces of flow of probability measures:
\begin{eqnarray*}
& C_p([0,T], \Pi) : =\left\{(\bar\pi_t)_{t\in [0,T]}: ~ \mbox{for each $t$,}~ \bar\pi_t\in \Pi,~\mbox{and}~ t \mapsto \bar\pi_t ~\mbox{is continuous in ${\cal W}_p$}\right\},\\
& \overline C_p([0,T], \Pi) : =\left\{(\bar\pi_t)_{t\in [0,T]}: ~ \mbox{for each $t$,}~ \bar\pi_t\in \Pi,~\mbox{and}~ t \mapsto \bar\pi_t ~\mbox{is continuous in $\overline{\cal W}_p$}\right\}.
\end{eqnarray*}

\begin{theorem}\label{thm:well}
For $i=1,\cdots, n$ and $\mu\in \Pi^{-i}$, let $F^i(\cdot, \mu): \nu \in\Pi^i\mapsto F^i(\nu, \mu)\in \mathbb{R}$  satisfy \cref{assum:base}. Further assume that
\begin{itemize}
\item the initial distribution $\bar\pi_0={\rm Law}(\bar X_0)\in \Pi$;
\item for each $i=1,\cdots,n$, the function $\nabla_{x^i} \frac{\delta F^i}{\delta \nu}$ is Lipschitz continuous in the following sense
\begin{multline*}
\left| \nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\nu, \mu,  x^i, {\mathsf y}) -\nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\nu', \mu', x'^i, {\mathsf y})\right|\\
\le C\left({\cal W}_p(\nu, \nu') + {\cal W}_p (\mu, \mu')  + | x^i -  x'^i |\right)
+ C_0 {\cal W}_p \big(\mu(\cdot|{\mathsf y}), \mu'(\cdot|{\mathsf y})\big),
\end{multline*}
and satisfies
\begin{equation}\label{eq:uniform0}
\sup_{\nu\in \Pi^i, \mu\in \Pi^{-i}, {\mathsf y}\in \mathbb{Y}} \left| \nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\nu, \mu,  0, {\mathsf y}) \right| <\infty.
\end{equation}
\end{itemize}
Then the MFL system \cref{eq:MFL} admits a unique strong solution in $\overline C_p([0,T], \Pi)$ for all $T>0$. In particular, if $C_0=0$ then the unique solution lies in $C_p([0,T], \Pi)$ for all $T>0$. Moreover,  each Nash equilibrium $\bar\pi^*$ defined in \cref{defn:nash} is an invariant measure of  \cref{eq:MFL}.
\end{theorem}

\begin{remark}
{\rm (i)}\quad The dependence on $\mu(\cdot|{\mathsf y})$ of the function $\nabla_{x^i}  \frac{\delta F^i}{\delta \nu}(\nu, \mu,  x^i, {\mathsf y}) $ is inevitable for some interesting examples such as \cref{nashcor} in the introduction, where under some mild conditions we may compute
\begin{equation*}
 \nabla_{x^i}  \frac{\delta F^i}{\delta \nu} \big(\nu^i, (\nu^j)_{j\neq i} , x^i, {\mathsf y}\big):=\int \nabla_{x^i} f^i(x^1 ,\cdots,x^i, \cdots, x^n, {\mathsf y})  \prod_{j\neq i}\nu^{j}(dx^j|{\mathsf y}).
\end{equation*}
Note that when there is only one player, there is no such dependence.

\vspace{3mm}
\noindent{\rm (ii)}\quad The Lipschitz condition with respect to ${\cal W}_p(\nu, \nu')$ and ${\cal W}_p(\mu, \mu')$ can be replaced by the one with respect to  $\overline{\cal W}_p(\nu, \nu')$ and $\overline{\cal W}_p(\mu, \mu')$. The latter is weaker.  Under such assumption we cannot prove the particular case that the unique solution  lies in $C_p([0,T], \Pi)$ for all $T>0$ when $C_0=0$. In the following analysis of the one-player problem it is crucial for us that the solution is in  $C_p([0,T], \Pi)$, so we prefer to state the Lipschitz condition in its current form.
\end{remark}

\begin{theorem}[Uniqueness of invariant measure: Contraction]\label{thm:contraction}
For $i=1,\cdots, n$ and $\mu\in \Pi^{-i}$, let $F^i(\cdot, \mu): \nu \in\Pi^i\mapsto F^i(\nu, \mu)\in \mathbb{R}$  satisfy \cref{assum:base}.
   Assume that
\begin{itemize}
\item  for each $i=1,\cdots,n$, the function $\nabla_{x^i} \frac{\delta F^i}{\delta \nu}$ is Lipschitz continuous in the following sense:
\begin{multline*}
\left| \left(\nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\bar\pi^i, \bar\pi^{-i},  x^i, {\mathsf y}) \right)_{i=1,\cdots,n}
- \left(\nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\bar\pi'^i, \bar\pi'^{-i}, x'^i, {\mathsf y}) \right)_{i=1,\cdots,n} \right|\\
\le \gamma \Big(  \overline{\cal W}_1(\bar\pi,\bar\pi') +{\cal W}_1\big(\pi(\cdot|{\mathsf y}),\pi'(\cdot|{\mathsf y})\big)   \Big) + C|x^i -x'^i|,
\end{multline*}
and \cref{eq:uniform0} holds true;
\item  there is a continuous function $ \kappa:(0,+\infty)\rightarrow\mathbb{R} $ s.t. $ \limsup\limits_{r\rightarrow+\infty}\kappa(r)<0$, $\int_{0}^{1}r\kappa(r) d r<+\infty  $ and for any $(\bar\pi,{\mathsf y})\in \Pi\times\mathbb{Y}$ we have for all $x, x'\in \mathbb{R}^N$ ($x\neq x'$)
\begin{equation*}
\sum_{i=1}^n (x^i -  x'^i)\cdot\Big(-\nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\bar\pi^i, \bar\pi^{-i}, x^i, {\mathsf y}) +\nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\bar\pi^i, \bar\pi^{-i},  x'^i, {\mathsf y})\Big)
\leq \kappa\left(|x -x'|\right)\left|x - x'\right|^2.
\end{equation*}
\end{itemize}
Let $\bar\pi_0, \bar\pi'_0 \in {\cal P}_q(\bar\mathbb{R}^d)\cap \Pi$ for some $q>1$ be two initial distributions of the MFL system \cref{eq:MFL}. Then we have
\begin{equation}\label{eq:contractionresult}
 \overline{\cal W}_1 (\bar\pi_t, \bar\pi'_t)
\le e^{(2\gamma-c\sigma^2) t}  \frac{2}{\varphi(R_1)}  \overline{\cal W}_1 (\bar\pi_0, \bar\pi'_0),
\end{equation}
where the coefficients read
\begin{eqnarray*}
	& \varphi(r) = \exp\left( -\frac12\int_{0}^{r}\frac{u\kappa^+(u)}{\sigma^2} d u \right),\quad
	c^{-1} = \int_{0}^{R_2}\Phi(r)\varphi(r)^{-1} d r,
	\quad \Phi(r) = \int_{0}^{r}\varphi(s) d s,\\
	& R_1 := \inf\{ R\geq 0: \kappa(r)\leq0 \text{ for all }r\geq R \},\\
&R_2 := \inf\{ R\geq R_1: \kappa(r)R(R-R_1)\leq-4\sigma^2\text{ for all }r\geq R \}.
\end{eqnarray*}
In particular, if $2\gamma < c\sigma^2$ (i.e. the MFL system has small dependence on the marginal laws), there is a unique invariant measure in $\cup_{q>1}{\cal P}_q(\bar\mathbb{R}^N)\cap\Pi$.
\end{theorem}

\subsection{Special case: one player}

When the problem degenerates to the case of a single player, the MFL dynamics becomes a gradient flow and the function $V=F+\frac{\sigma^2}{2}H$ is a natural Lyapunov  function for the dynamics.

\begin{theorem}[Gradient flow]\label{thm:GD}
Consider a function $F$ satisfying \cref{assum:base} with $p=2$. Let the assumption of \cref{thm:well} hold true, and further assume that
\begin{itemize}
\item  there is $\varepsilon>0$ such that  for all $\bar\pi\in \Pi$ and ${\mathsf y}\in \mathbb{Y}$
\begin{equation}\label{assum:dissipative}
 x \cdot \nabla_x \frac{\delta F}{\delta \nu}(\bar\pi, x, {\mathsf y}) \ge  \varepsilon |x|^2, \quad\mbox{for $|x|$  big enough} ;
 \end{equation}

\item for all $\bar\pi\in \Pi$ and ${\mathsf y}\in \mathbb{Y}$, the mapping $x\mapsto  \nabla_x \frac{\delta F}{\delta \nu}(\bar\pi, x, {\mathsf y})$ belongs to $C^3$;

\item for all ${\mathsf y}\in \mathbb{Y}$, the function $(\bar\pi, x)\mapsto \nabla_x \frac{\delta F}{\delta \nu}(\bar\pi, x, {\mathsf y}), \nabla_x^2 \frac{\delta F}{\delta \nu}(\bar\pi, x, {\mathsf y})$ are jointly continuous.
\end{itemize}
Then we have for $s'>s>0$
\begin{equation}\label{eq:GD}
 V(\bar\pi_{s'}) -V(\bar\pi_s) =-\int_s^{s'}  \int_{\bar\mathbb{R}^N} \left| \nabla_{x} \frac{\delta F}{\delta \nu}(\bar\pi_t,  x,{\mathsf y}) + \frac{\sigma^2}{2} \nabla_{x}\ln\big(\pi_t(x |{\mathsf y})\big) \right|^2  \bar\pi_t (d\bar x) dt
\end{equation}
\end{theorem}

Using an argument, similar to that in \cite{HRSS19, HKR19},  based on the Lasalle's invariant principle, we can show  the following theorem.
\begin{theorem}\label{thm:invariant}
Consider the following statements:
\begin{itemize}
\item[{\rm (i)}] $\bar\pi_0\in \cup_{q>2}{\cal P}_q(\bar\mathbb{R}^N)$;
\item[{\rm (ii.a)}] $\mathbb{Y}$ is countable;
\item[{\rm (ii.b)}] $\mathbb{Y} = \mathbb{R}^m$, ${\mathsf m}$ is absolutely continuous with respect to the Lebesgue measure  and the function $e^{-\frac{2}{\sigma^2}\frac{\delta F }{\delta \nu}(\bar\pi, \cdot)}{\mathsf m}$ is semiconvex in $\bar x$ for any given $\bar\pi\in \Pi$.
\end{itemize}
Let  the assumptions of \cref{thm:GD} hold true. Further assume {\rm (i)}, {\rm (ii.a)} or {\rm (i)}, {\rm (ii.b)}.
Then all the ${\cal W}_2$-cluster points of the marginal laws $(\bar\pi_t)_{t\ge 0}$ of the MFL system \cref{eq:MFL} belong to  the set
\begin{equation}\label{eq:limitset}
{\cal I}:=\left\{ \bar\pi \in\Pi:~ \nabla_{x} \frac{\delta F}{\delta \nu}(\bar\pi,  \cdot,{\mathsf y}) + \frac{\sigma^2}{2} \nabla_{x}\ln\big(\pi(\cdot|{\mathsf y})\big) =0\quad\mbox{${\mathsf m}$-a.s.} \right\}.
\end{equation}
\end{theorem}

\begin{remark}[The limit set and the mean-field equilibria on the environment]
Consider the case where the probability measure on the environment ${\mathsf m}$ is atomless. In particular, for a fixed ${\mathsf y}\in \mathbb{Y}$ the probability $\bar\pi\in \Pi$ does not depend on $\pi(\cdot |{\mathsf y} )$. Therefore the equation in \cref{eq:limitset} is a sufficient and necessary condition for
\begin{equation}\label{nash_x0}
\pi(\cdot|{\mathsf y}) = \argmin_{\nu}  \left( \int  \frac{\delta F}{\delta \nu}(\bar\pi, x ,{\mathsf y}) \nu(dx) + \frac{\sigma^2}{2} H(\nu|{\rm Leb} )\right),
\quad\mbox{for ${\mathsf m}$-a.s. ${\mathsf y}$},
\end{equation}
where $H(\cdot|{\rm Leb})$ is the relative entropy with respect to the Lebesgue measure.
If we view the variable ${\mathsf y}$ as the index of the `players', \cref{nash_x0} indicates that all $\bar\pi\in {\cal I}$ are (mean-field) Nash equilibria of the game where the ${\mathsf y}$-player aims at:
\begin{equation*}
\inf_{\nu}  \left( \int  \frac{\delta F}{\delta \nu}(\bar\pi, x ,{\mathsf y}) \nu(dx) + \frac{\sigma^2}{2} H(\nu|{\rm Leb})\right).
\end{equation*}
\end{remark}

\begin{corollary}\label{thm:convex_1P}
If the function $V$ is convex, the limit set ${\cal I}$ is a singleton and thus the marginal laws $(\bar\pi_t)_{t\ge 0}$ converge in ${\cal W}_2$ to the minimizer of $V$.
\end{corollary}

\section{Applications}\label{sec:app}

\subsection{Dynamic games and deep neural networks}

As mentioned in \cref{eg:dynamicgame}, both discrete-time and continuous-time dynamic games  can be viewed as games on the random environment.

Take the continuous-time dynamic game as an example, in particular $\mathbb{Y}=[0,T]$. Consider the controlled process of the $i$-th player
\begin{equation}\label{eq:controlledproc}
d\Theta^i_{\mathsf y} = \int \varphi^i \big( x^i, \bar \pi^{-i}, \Theta_{\mathsf y}^i, {\mathsf y}\big) \pi^i(dx^i | {\mathsf y}) d{\mathsf y}.
\end{equation}
and his objective function
\begin{equation*}
F^i(\bar\pi^i, \bar\pi^{-i}) := \int_0^T c^i \big( x^i, \bar \pi^{-i}, \Theta_{\mathsf y}^i, {\mathsf y}\big) \pi^i(dx^i|{\mathsf y}) d{\mathsf y} + g^i(\Theta^i_T).
\end{equation*}
Define the Hamiltonian function $H^i(x^i, \mu, \theta^i, {\mathsf y}, p) := c^i(x^i, \mu, \theta^i, {\mathsf y}) + p \cdot \varphi^i (x^i, \mu, \theta^i, {\mathsf y}) $.  Assume that
\begin{itemize}
\item the coefficients $\varphi^i, c^i $ are uniformly Lipschitz in $(x^i, \theta^i)$;

\item $\varphi^i, c^i, g^i$ are continuously differentiable in $\theta^i$;

\item  $\nabla_{\theta^i} \varphi^i, \nabla_{\theta^i} c^i $ are uniformly Lipschitz in $(x^i, \theta^i)$, and $\nabla_{\theta^i} g^i$ is uniformly Lipschitz in $x^i$.
\end{itemize}
\noindent It follows from a standard variational calculus that
\begin{equation*}
\frac{\delta F^i}{\delta \nu}(\nu, \mu,x^i, {\mathsf y}) = H^i(x^i, \mu, \Theta^i_{\mathsf y}, {\mathsf y}, P^i_{\mathsf y}),
\end{equation*}
where $\Theta^i$ follows the dynamics \cref{eq:controlledproc} and $P^i$ is the solution to the linear ODE:
\begin{equation*}
P^i_{\mathsf y} = \nabla_\theta g(\Theta^i_T) + \int_{\mathsf y}^T \nabla_\theta H^i(x^i, \mu, \Theta^i_{\mathsf y}, {\mathsf y}, P^i_{\mathsf y})  \pi^i(dx^i | {\mathsf y}) d{\mathsf y}.
\end{equation*}
Therefore, according to \cref{thm:contraction}, the Nash equilibrium of this dynamic game can be approximated by the marginal law of the MFL system.

In case the number of players $n=1$, the marginal laws of the MFL system approximates the minimizer of the optimization.  There is a rising interest in modeling the forward propagation of the deep neural networks using a controlled dynamics and in connecting the deep learning to the optimal control problems, see e.g. \cite{CLT19, CRBD19, EHL19,JSS19,  LM19, LMLLY20}.  For the controlled processes in the particular form \cref{eq:controlledproc},  we refer to Section 4 in \cite{HKR19} for the connection between the optimal control problem and the deep neural networks. In particular, we remark that the backward propagation algorithm is simply a discretization of the corresponding MFL dynamics.

\subsection{Linear-convex zero-sum game and GAN}

Consider the zero-sum game between two players, i.e. $F^1 = -F^2 (=:F)$. For $F$ satisfying the assumption in \cref{thm:contraction} so that the contraction result \eqref{eq:contractionresult} holds true, we may use the following MFL system to approximate the unique Nash equilibrium:
\begin{eqnarray}
\begin{cases}\label{eq:zerosumMFL2}
 d X^1_t = - \nabla_{x^1} \frac{\delta F}{\delta \nu} (\bar\pi^1_t , \bar\pi^2_t, X^1_t, Y)dt + \sigma dW^1_t,\\
 d X^2_t =  \nabla_{x^2} \frac{\delta F}{\delta \nu} (\bar\pi^2_t , \bar\pi^1_t, X^2_t, Y)dt + \sigma dW^2_t.
\end{cases}
\end{eqnarray}

Now we consider a particular subclass of the zero-sum games. Assume that $F:(\nu, \mu)\in \Pi^1\times \Pi^2\mapsto F(\nu, \mu)\in \mathbb{R}$ is linear in $\mu$ and convex in $\nu$,  and define $\tilde V(\nu,\mu):= F(\nu,\mu) + \frac{\sigma^2}{2} \big(H(\nu) -H(\mu)\big)$. In particular,
\begin{itemize}
\item $\frac{\delta F}{\delta \mu}$ does not depend on $\mu$;
\item the function $\Phi:\nu\mapsto\max_{\mu\in \Pi^{2}} \tilde V(\nu, \mu) - \frac{\sigma^2}{2} H(\nu)$ is convex.
\end{itemize}
If the Nash equilibrium exists, denoted by $\bar\pi^*$, by the standard argument we have

\begin{equation*}
\min_{\nu\in \Pi^{1}}\max_{\mu\in \Pi^{2}} \tilde V(\nu, \mu) = \tilde  V(\bar\pi^{*,1}, \bar\pi^{*,2})=\max_{\mu\in \Pi^{2}} \min_{\nu\in \Pi^{1}} \tilde V (\nu, \mu).
\end{equation*}
It follows from \cref{thm:FOC} that $\mu^*[\nu] := \argmax_{\mu\in \Pi^{2}} \tilde V(\nu, \mu)$ has the explicit density
\begin{equation*}
\frac{ \mu^*[\nu]  (x^2 , d{\mathsf y})}{{\mathsf m}(d{\mathsf y})}= C(\nu, {\mathsf y}) e^{-\frac{2}{\sigma^2} \frac{\delta F}{\delta \mu}(\nu, x^2, {\mathsf y})},
\end{equation*}
where $C(\nu, {\mathsf y})$ is the normalization constant. Further, assume that the function $\Phi(\nu) = \tilde V(\nu, \mu^*[\nu])$ satisfies the assumption of \cref{thm:invariant} and recall that $\Phi$ is convex, it follows from \cref{thm:convex_1P} that we may approximate the minimizer $\bar\pi^{*,1}$ using the dynamics:
\begin{equation}\label{zerosumMFL}
dX_t = -\nabla_x\frac{\delta \Phi}{\delta \nu} (\bar\pi_t, X_t, Y) dt + \sigma dW_t.
\end{equation}
Compared to the dynamics \eqref{eq:zerosumMFL2}, the dynamics \cref{zerosumMFL} enjoys the natural Lyapunov function $\tilde V$, i.e. $ t \mapsto \tilde V( \bar \pi_t ) $ decreases monotonically.
\vspace{3mm}

As an application, the generative adversarial networks (GAN) can be viewed as a linear-convex game. Given a bounded, continuous, non-constant activation function $\varphi$, consider the parametrized functions
\begin{equation}\label{eq:parafunGAN}
\{z\mapsto \mathbb{E}[\varphi(X,z)]: ~{\rm Law}(X) = \nu\in {\cal P}_2(\mathbb{R}^{n^2}) \}
\end{equation}
as the options of the discriminators. The regularized GAN aims at computing the Nash equilibrium of the game:
\begin{eqnarray*}
&\mbox{Define}\quad \tilde V(\nu, \mu):= - \int   \mathbb{E}[\varphi(X,z)]  (\mu-\hat\mu)(dz) - \frac12 \lambda \Big(\int |z|^2 \mu(dz) - \mathbb{E}[|X|^2] \Big) - \frac{\sigma^2}{2}\big(H(\mu)-H(\nu)\big)&\\
&\begin{cases}
{\rm Generator:} & \sup_{\mu\in {\cal P}_2(\mathbb{R}^{n^1})} \tilde V(\nu, \mu) \\
{\rm Discriminator:} & \inf_{\nu\in {\cal P}_2(\mathbb{R}^{n^2})} \tilde V(\nu, \mu)
\end{cases},&
\end{eqnarray*}
where $\hat\mu\in {\cal P}_2(\mathbb{R}^{n^1})$ is the distribution of interest. Indeed, in order to compute the Nash equilibrium of the game, it is appealing to sample the MFL dynamics \eqref{eq:zerosumMFL2} and approximate its invariant measure. However, in order that the contraction result in  \cref{thm:contraction} holds true, it is crucial that the MFL system has small dependence on the marginal laws, which is not necessarily true in the context of GAN.  Here we present another approach, which exploits the particular structure of the linear-convex game. As discussed before, the optimizer of the generator given $\nu\in {\cal P}_2(\mathbb{R}^{n^2})$ is explicit and has the density:
\begin{equation}
\label{generator_explicit}
\mu^*[\nu](z) = C(\nu) e^{-\frac{2}{\sigma^2} \big( \mathbb{E}[\varphi(X,z)] +\frac{\lambda}{2} |z|^2 \big)}.
\end{equation}
Further for the potential function $\Phi(\nu):= \tilde V(\nu, \mu^*[\nu]) - \frac{\sigma^2}{2} H(\nu)$ we have
\begin{equation}
\label{discrim_formula}
\frac{\delta \Phi}{\delta \nu}(\nu, x)= -\int \varphi(x, z) (\mu^*[\nu]-\hat\mu)(dz) + \frac{\lambda}{2}|x|^2.
\end{equation}
Then the strategy of the discriminator in the Nash equilibrium can be approximated by the MFL dynamics \cref{zerosumMFL}.  In the perspective of numerical realization, note that the law $\mu^*[\nu]$ can be simulated by the MCMC algorithms such as Metropolis-Hastings.

\vspace{3mm}

In order to illustrate the advantage of the algorithm using the MFL dynamics \cref{zerosumMFL}, here we present the numerical result for a toy example. We are going to use the GAN to generate the samples of the exponential distribution with intensity $1$. In this test, the optimal response of the generator, $\mu^*[\nu]$, is computed via Metropolis algorithm with Gaussian proposal distribution with zero mean and variance optimised according to \cite{metropolis_efficient}. The discriminator chooses parametrized functions among \cref{eq:parafunGAN}, where
\begin{equation*}
\varphi( X, z ) = C ( A z + b )^+,\quad  \text{ with } X = ( C, A, b ).
\end{equation*}
When we numerically run the MFL dynamics  \cref{zerosumMFL} to train the discriminator, we use a $3000$-particle system, that is, the network is composed of $3000$ neurones, and set its initial distribution to be standard Gaussian.  The other parameters are chosen as follows: $ \sigma = 0.4 $, $ dt = 0.01 $, $ \lambda = 0.2 $. \Cref{fig:exponential} shows the training result after $ 60 $ iterations.
\begin{figure}[h]
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.9\linewidth]{histogram_exp_2.png}
  \caption{Histogram of GAN's sampling and the target density}
  \label{fig:sub1}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.9\linewidth]{training_error_2.png}
  \caption{Training error}
  \label{fig:sub2}
\end{subfigure}
\caption{Learning via MCMC-GAN: histogram of learned distribution and training error.}
\label{fig:exponential}
\end{figure}
 In particular, we see that the training error decreases monotonically as suggested by our theoretical results.

\section{Proofs}\label{sec:proof}

\subsection{Optimization with marginal constraint}

\begin{proof}[Proof of \cref{thm:FOC}] \quad  \underline{Necessary condition} ~ {\it Step 1.}\quad Let $\bar\pi^*\in \Pi$ be a minimizer of $V$. Since $H(\bar\pi^*)<\infty$, the probability measure $\bar\pi^*$ is absolutely continuous wrt ${\rm Leb}\times {\mathsf m}$. Take any probability measure $\bar\pi\in \Pi $ such that $H(\bar\pi)<\infty$, in particular $\bar\pi$ is also absolutely continuous wrt ${\rm Leb}\times {\mathsf m}$. Denote the convex combination by $\bar\pi^\varepsilon: = \varepsilon \bar\pi + (1-\varepsilon) \bar\pi^* \in \Pi$. Define the function $h(z): = z\ln z$ for $z\in\mathbb{R}^+$ and $h(0)=0$. Then
\begin{eqnarray*}
0\le \frac{V(\bar\pi^\epsilon)-V(\bar\pi^*)}{\epsilon} &=& \frac{F(\bar\pi^\epsilon)-F(\bar\pi^*)}{\epsilon} + \eta\int\frac{h(\pi^\epsilon( x|{\mathsf y}))  - h(\pi^*( x|{\mathsf y}))}{\varepsilon} dx {\mathsf m}(d{\mathsf y})  \\
&= &\frac{1}{\epsilon}\int_{0}^{\epsilon}\int\frac{\delta F}{\delta \bar\pi}(\bar\pi^\lambda,\bar x)(\bar\pi-\bar\pi^*)(d\bar x) d\lambda
	+ \eta\int\frac{h(\pi^\epsilon( x|{\mathsf y}))  - h(\pi^*( x|{\mathsf y}))}{\varepsilon} dx {\mathsf m}(d{\mathsf y}) .
\end{eqnarray*}
Since $\sup_{\lambda\in [0,\varepsilon]}|\frac{\delta F}{\delta \bar\pi}(\bar\pi^\lambda,\bar x)| \le C(1+|\bar x|^p)$ and $\bar\pi,\bar\pi^*\in \Pi$, by the dominated convergence theorem
\begin{equation*}
\lim_{\varepsilon\rightarrow 0} \frac{1}{\epsilon}\int_{0}^{\epsilon}\int\frac{\delta F}{\delta \bar\pi}(\bar\pi^\lambda,\bar x)(\bar\pi-\bar\pi^*)(d \bar x) d\lambda
= \int\frac{\delta F}{\delta \bar\pi}(\bar\pi^*,\bar x)(\bar\pi-\bar\pi^*)(d \bar x).
\end{equation*}
Since the function $h$ is convex, we have $\frac{h(\pi^\epsilon( x|{\mathsf y}))  - h(\pi^*( x|{\mathsf y}))}{\varepsilon} \le h(\pi( x|{\mathsf y}))  - h(\pi^*( x|{\mathsf y}))$. Note that $\int \big(h(\pi( x|{\mathsf y}))  - h(\pi^*( x|{\mathsf y}))\big) d x{\mathsf m}({\mathsf y}) = H(\bar\pi)-H(\bar\pi^*)<\infty$. By Fatou lemma, we obtain
\begin{eqnarray*}
\limsup_{\varepsilon\rightarrow 0}\int\frac{h(\pi^\epsilon( x|{\mathsf y}))  - h(\pi^*( x|{\mathsf y}))}{\varepsilon} d  x {\mathsf m}(d{\mathsf y})
&\le& \int \lim_{\varepsilon\rightarrow 0} \frac{h(\pi^\epsilon( x|{\mathsf y}))  - h(\pi^*( x|{\mathsf y}))}{\varepsilon}  d  x {\mathsf m}(d{\mathsf y}) \\
&=&\int \ln\pi^*(x|{\mathsf y}) (\bar\pi-\bar\pi^*)(d\bar x).
\end{eqnarray*}
Therefore we have
\begin{equation}\label{FOCineq}
0\le \limsup_{\varepsilon\rightarrow 0}\frac{V(\bar\pi^\epsilon)-V(\bar\pi^*)}{\epsilon}
\le \int \left( \frac{\delta F}{\delta \bar\pi}(\bar\pi^*,\bar x) + \eta\ln\pi^*(x|{\mathsf y})  \right)(\bar\pi-\bar\pi^*)(d\bar x).
\end{equation}

\vspace{3mm}

\noindent {\it Step 2.}\quad  We are going to show that for ${\mathsf m}$-a.s. ${\mathsf y}$
\begin{eqnarray}\label{eq:nece_constant}
\Xi_{{\mathsf y}}(x):= \frac{\delta F}{\delta \bar\pi}(\bar\pi^*,\bar x) + \eta \ln\pi^*( x|{\mathsf y}) \quad\mbox{is equal to a constant $\pi^*(\cdot |{\mathsf y})$-a.s.}
\end{eqnarray}
Define the mean value $\overline c({\mathsf y}) : = \int_{\mathbb{R}^d} \Xi_{{\mathsf y}}(x) \pi^*(dx|{\mathsf y})$ and let $\varepsilon, \varepsilon'>0$. Consider the probability measure $\bar\pi\in \Pi$ absolutely continuous  wrt $\bar\pi^*$ such that
\begin{eqnarray*}
\frac{d\pi (\cdot|{\mathsf y})}{d\pi^*(\cdot|{\mathsf y})} =
\begin{cases}
1, & \mbox{for ${\mathsf y}$ such that} ~\pi^*\big(\Xi_{{\mathsf y}}\le \overline c({\mathsf y}) -\varepsilon\big| {\mathsf y}\big) <\varepsilon' \\
\frac{1_{\Xi_{{\mathsf y}} \le \overline c({\mathsf y}) -\varepsilon}}{\pi^*\big(\Xi_{{\mathsf y}} \le \overline c({\mathsf y}) -\varepsilon\big| {\mathsf y}\big) }, & \mbox{otherwise}
\end{cases}.
\end{eqnarray*}
Since $\frac{d\bar\pi}{d\bar\pi^*}$ is bounded, we have that $\bar\pi\in \Pi$ and $H(\bar\pi)<\infty$. In particular \cref{FOCineq} holds true for this $\bar\pi$.
Also note that $\Xi_{{\mathsf y}} \le \overline c({\mathsf y})-\varepsilon$, $\pi(\cdot| {\mathsf y})$-a.s. for ${\mathsf y}$ such that $\pi^*\big(\Xi_{{\mathsf y}} \le \overline c({\mathsf y}) -\varepsilon\big| {\mathsf y}\big) \ge \varepsilon'$. So we have
\begin{eqnarray*}
0&\le &\int_\mathbb{Y} \int_{\mathbb{R}^d} \Xi_{{\mathsf y}}(x) \big(\pi(dx|{\mathsf y}) -\pi^*(dx |{\mathsf y})\big){\mathsf m}(d{\mathsf y}) \\
&=&\int_{\pi^*\big(\Xi_{{\mathsf y}}\le \overline c({\mathsf y}) -\varepsilon\big| {\mathsf y}\big) \ge \varepsilon'} \int_{\mathbb{R}^d} \Xi_{{\mathsf y}}(x) \big(\pi(dx|{\mathsf y}) -\pi^*(dx |{\mathsf y})\big){\mathsf m}(d{\mathsf y}) \\
& = & \int_{\pi^*\big(\Xi_{{\mathsf y}}\le \overline c({\mathsf y}) -\varepsilon\big| {\mathsf y}\big) \ge \varepsilon'}\left( \int_{\mathbb{R}^d} \Xi_{{\mathsf y}}(x) \pi(dx|{\mathsf y}) -\overline c({\mathsf y})\right){\mathsf m}(d{\mathsf y}) \\
&\le& -\varepsilon ~{\mathsf m}\left\{{\mathsf y}:\pi^*\big(\Xi_{{\mathsf y}}\le \overline c({\mathsf y}) -\varepsilon\big| {\mathsf y}\big) \ge \varepsilon'\right\}.
\end{eqnarray*}
Therefore we conclude that
$\pi^*\big(\Xi_{{\mathsf y}}\le \overline c({\mathsf y}) -\varepsilon\big| {\mathsf y}\big) <\varepsilon'$ for ${\mathsf m}$-a.s. ${\mathsf y}$.
Since this is true for arbitrary $\varepsilon', \varepsilon>0$, we obtain \eqref{eq:nece_constant}.

\vspace{3mm}
\noindent{\it Step 3}.\quad We are going to show that $\bar\pi^*$ is equivalent to  ${\rm Leb}\times {\mathsf m}$, so that $\Xi_{\mathsf y}$ does not depend on $x$, ${\rm Leb}\times {\mathsf m}$-a.s. and the first order equation \cref{eq:necessary}  holds true.
Suppose the opposite, i.e. there is a set ${\cal K}\in \bar\mathbb{R}^d$ such that $\bar\pi^*({\cal K})=0$ and ${\rm Leb}\times{\mathsf m}({\cal K})>0$. In particular,  $\ln\pi^*(x|{\mathsf y}) = -\infty$ on ${\cal K}$.
Denote ${\cal K}_{\mathsf y}:=\{x\in \mathbb{R}^d: (x,{\mathsf y})\in {\cal K}\}$. We may assume that there exist $K>\varepsilon>0$ such that $ {\rm Leb}({\cal K}_{\mathsf y}) \in [\varepsilon, K]$ for all ${\mathsf y}\in \mathbb{Y}$.
Define a probability measure $\bar\pi\in \Pi$ such that for all Borel-measurable $A\subset \mathbb{R}^d$
\begin{equation*}
\pi (A|{\mathsf y}):= \frac12 \pi^*(A|{\mathsf y}) + \frac{1}{2 {\rm Leb}({\cal K}_{\mathsf y})}\int_{A\cap {\cal K}_{\mathsf y}} dx.
\end{equation*}
It is easy to verify that  $H(\bar\pi)<\infty$, so \cref{FOCineq} holds true and it implies
\begin{eqnarray*}
0&\le& \frac12 \int_{\cal K} \left( \frac{\delta F}{\delta \bar\pi}(\bar\pi^*,\bar x) + \eta\ln\pi^*(x|{\mathsf y})  \right)\bar\pi(d\bar x) - \frac12\int \left( \frac{\delta F}{\delta \bar\pi}(\bar\pi^*,\bar x) + \eta\ln\pi^*(x|{\mathsf y})  \right) \bar\pi^*(d\bar x)\\
&\le & -\infty + \int C(1+|\bar x|^p)  \bar\pi^*(d\bar x) -\frac{\eta}{2} H(\bar\pi^*) =-\infty.
\end{eqnarray*}
It is a contradiction, so $\bar\pi^*$ is equivalent to  ${\rm Leb}\times {\mathsf m}$.

\vspace{5mm}

\noindent\underline{Sufficient condition}~ Assume that $F$ is convex. Let $ \bar\pi^*\in \Pi $  satisfy the first order equation \cref{eq:necessary}, in particular, $\bar\pi^*$ is equivalent to ${\rm Leb}\times {\mathsf m}$.
Take any $ \bar\pi\in \Pi  $ absolutely continuous wrt ${\rm Leb}\times {\mathsf m}$ (otherwise $ V(\bar\pi)=+\infty $), and thus absolutely continuous wrt the measure $\bar\pi^*$.
Define $\bar\pi^\epsilon := (1-\epsilon)\bar\pi^* + \epsilon \bar\pi$ for $\epsilon>0$.
By the convexity of $F$ we obtain
\begin{eqnarray*}
F(\bar\pi) - F(\bar\pi^*)
&\ge& \lim_{\varepsilon\rightarrow 0} \frac{1}{\varepsilon} \big(F(\bar\pi^\varepsilon) - F(\bar\pi^*)\big) \\
&=&\lim_{\varepsilon\rightarrow 0} \frac{1}{\varepsilon}\int_0^\varepsilon \int_{\bar\mathbb{R}^d} \frac{\delta F}{\delta \bar\pi} (\bar\pi^\lambda, \bar x)  (\bar\pi -\bar\pi^*) (d\bar x)d\lambda
=\int_{\bar\mathbb{R}^d} \frac{\delta F}{\delta \bar\pi}(\bar\pi^*,\bar x ) (\bar\pi -\bar\pi^*) (d\bar x).
\end{eqnarray*}
The last equality is due to the dominated convergence theorem. On the other hand, by convexity of the function $ h $,
\begin{equation*}
H(\bar\pi) - H(\bar\pi^*) \ge \int_{\bar\mathbb{R}^d} \ln\pi^*(x|{\mathsf y}) (\bar\pi -\bar\pi^*)(d\bar x).
\end{equation*}
Hence
\begin{equation*}
V(\bar\pi)-V(\bar\pi^*)
\geq \int_{\bar\mathbb{R}^d} \bigg(\frac{\delta F}{\delta \bar\pi}(\bar\pi^*,\cdot) +\eta \ln \pi^*(x|{\mathsf y})  \bigg)(\bar\pi -\bar\pi^*) (d\bar x) = 0,
\end{equation*}
so $\bar\pi^*$ is a minimizer.
\end{proof}

\subsection{Equilibria of game}

\begin{proof}[Proof of \cref{cor:uniq_monot}] \quad
{\rm (i)}\quad Let $n=1$. Take three probability measures $\bar\pi,\bar\pi', \bar\pi''\in \Pi$ such that $\bar\pi =\frac12(\bar\pi'+\bar\pi'')$. Denote $\bar\pi'^{\lambda}:= \lambda \bar\pi+ (1-\lambda)\bar\pi'$ and $ \bar\pi''^{\lambda}:=\lambda \bar\pi+ (1-\lambda) \bar\pi''$. By the definition of the linear derivative of $F$ we obtain
\begin{equation*}
F(\bar\pi') - 2 F(\bar\pi)+F(\bar\pi'') = \int_0^1 \int_{\bar\mathbb{R}^N} \frac{\delta F}{\delta \nu} (\bar\pi'^\lambda, \bar x ) (\bar\pi'-\bar\pi)(d\bar x) d\lambda
					+ \int_0^1 \int_{\bar\mathbb{R}^N} \frac{\delta F}{\delta \nu} (\bar\pi''^\lambda, \bar x ) (\bar\pi''-\bar\pi)(d\bar x) d\lambda.
\end{equation*}
Note that $\bar\pi' - \bar\pi = \bar\pi-\bar\pi'' =\frac{1}{2-2\lambda} \big( \bar\pi'^\lambda-\bar\pi''^\lambda\big)$. Therefore we have
\begin{equation*}
F(\bar\pi') - 2 F(\bar\pi)+F(\bar\pi'') =  \int_0^1 \frac{1}{2-2\lambda} \int_{\bar\mathbb{R}^N} \left( \frac{\delta F}{\delta \nu} (\bar\pi'^\lambda, \bar x )-\frac{\delta F}{\delta \nu} (\bar\pi''^\lambda, \bar x ) \right) (\bar\pi'^\lambda-\bar\pi''^\lambda) (d\bar x) d\lambda \geq 0.
\end{equation*}
Finally note that $\lambda\mapsto (F(\bar\pi'^\lambda), F(\bar\pi''^\lambda))$ is continuous. So $F$ satisfying the monotonicity condition must be convex on $\Pi$.

On the other hand, suppose $F$ is convex on $\Pi$.  Following a similar computation, we  obtain
\begin{equation*}
0\le \frac{1}{\varepsilon} \Big( F(\bar\pi')- F(\bar\pi'^\varepsilon) - F(\bar\pi''^\varepsilon) + F(\bar\pi'') \Big)\\
= \frac{1}{\varepsilon}\int_0^\varepsilon \frac{1}{2} \int_{\bar\mathbb{R}^N}  \left( \frac{\delta F}{\delta \nu} (\bar\pi'^\lambda, \bar x )-\frac{\delta F}{\delta \nu} (\bar\pi''^\lambda, \bar x ) \right) (\bar\pi'-\bar\pi'') (d\bar x) d\lambda.
\end{equation*}
Let $\bar\pi', \bar\pi''\in \Pi$. It follows from the dominated convergence theorem that
\begin{eqnarray*}
0 &\le& \lim_{\varepsilon\rightarrow 0}  \frac{1}{\varepsilon}\int_0^\varepsilon \frac{1}{2} \int_{\bar\mathbb{R}^N}  \left( \frac{\delta F}{\delta \nu} (\bar\pi'^\lambda, \bar x )-\frac{\delta F}{\delta \nu} (\bar\pi''^\lambda, \bar x ) \right) (\bar\pi'-\bar\pi'') (d\bar x) d\lambda\\
&=& \frac12 \int_{\bar\mathbb{R}^N}  \left( \frac{\delta F}{\delta \nu} (\bar\pi', \bar x )-\frac{\delta F}{\delta \nu} (\bar\pi'', \bar x ) \right) (\bar\pi'-\bar\pi'') (d\bar x).
\end{eqnarray*}
So $F$ satisfies the monotonicity condition on $\Pi$.

\vspace{3mm}

\noindent {\rm (ii)}\quad Let $ \bar\pi^* \in \Pi$ be a Nash equilibrium of the game.
Then, by \cref{cor:NC_Nash} we have that for every $ i $ there exists a function $ f^i : \mathbb{Y} \rightarrow \mathbb{R}$
\begin{equation*}
\frac{ \delta F^i }{ \delta \nu } ( ( \bar \pi^* )^{i}, ( \bar \pi^* )^{-i}, x^i, {\mathsf y} ) + \frac{ \sigma^2 }{ 2 } \ln( ( \pi^* )^i( x^i | {\mathsf y} )  ) = f^i( {\mathsf y} ),
\quad\mbox{for ${\mathsf m}$-a.s. ${\mathsf y}$.}
\end{equation*}
Let $ \bar \pi^*, \bar \pi'^* \in \Pi$ be Nash equilibriums. Then monotonicity condition \cref{eq:monotonicity} implies
\begin{equation*}
\sum_{ i = 1 }^n \int\limits_{ \bar \mathbb{R}^N } \left( - \frac{ \sigma^2 }{ 2 } \ln( (\pi^*)^i( x^i | {\mathsf y} ) ) + \frac{ \sigma^2 }{ 2 } \ln( (\pi'^*)^{i}( x^i | {\mathsf y} ) )  \right) ( \bar \pi^* - \bar \pi'^* )( d \bar x ) \geq 0,
\end{equation*}
because the marginal distributions of $ \bar \pi^* $ and $ \bar \pi'^* $ on $ \mathbb{Y} $ coincide.
The latter inequality can be rewritten
\begin{eqnarray*}
0 &\le & - \sum_{ i = 1 }^n \int\limits_{ \bar \mathbb{R}^{ n^i } } \ln \left( \frac{ ( \pi^* )^i( x^i | {\mathsf y} )  }{ ( \pi'^* )^i( x^i | {\mathsf y} ) } \right) ( ( \bar \pi^*)^i - ( \bar \pi'^* )^i )( d \bar x^i )\\
& = & -\sum_{ i = 1 }^n \left( H( ( \bar \pi^* )^i | ( \bar\pi'^* )^i ) + H( ( \bar \pi'^* )^i | ( \bar \pi^* )^i ) \right).
\end{eqnarray*}
This is only possible if $ ( \bar \pi^* )^i = ( \bar \pi'^* )^i $ for all $ i = 1, \ldots, n $.
\end{proof}

\vspace{3mm}

\begin{proof}[Proof of \cref{thm:exist_nash}] \quad  For $\bar\pi\in \Pi$ denote ${\cal R}^i(\bar\pi):=\argmin_{\nu\in \Pi^i} V^i(\nu, \bar\pi^{-i})$ for $i=1,\cdots,n$, and define
\begin{equation*}
{\cal R}(\bar\pi):=\{\bar\pi'\in \Pi : ~\bar\pi'^i \in {\cal R}^i(\bar\pi)~\mbox{for $i=1, \cdots, n$}\}.
\end{equation*}

\noindent{\it Step 1.}\quad First we prove that ${\cal R}(\bar\pi)$ is ${\cal W}_p$-compact. For any $\nu^i \in {\cal R}^i(\bar\pi)$, it follows from the first order equation \cref{eq:necessary} that
\begin{equation*}
\nu^i( \bar x^i) = C( \nu^i, \bar\pi^{-i}) e^{-\frac{2}{\sigma^2} \frac{\delta F^i}{\delta \nu}(\nu^i, \bar\pi^{-i}, \bar x^i)},
\end{equation*}
where $C( \nu^i, \bar\pi^{-i})>0$ is the normalization constant so that $\nu^i$ is a probability measure. Take a $p'>p$. The condition \eqref{basic_dissipative} implies that $C( \nu^i, \bar\pi^{-i})$ is uniformly bounded as well as
\begin{equation*}
\int_{\bar\mathbb{R}^{n^i}} |\bar x^i |^{p'} \nu^i(d\bar x^i)
\le C( \nu^i, \bar\pi^{-i}) \int_{\bar\mathbb{R}^{n^i}} |\bar x^i |^{p'} e^{- C'|\bar x^i|^q +C} dx^i {\mathsf m}(d{\mathsf y}).
\end{equation*}
So we have
\begin{equation*}
\overline C^i:= \sup_{\nu^i \in {\cal R}^i(\bar\pi) }\int_{\bar\mathbb{R}^{n^i}} |\bar x^i |^{p'} \nu^i(d\bar x^i) <\infty.
\end{equation*}
Therefore, ${\cal R}(\bar\pi) \subset  {\cal E} :=\{\nu\in \Pi: \nu^i\in {\cal E}^i \mbox{ for $i=1,\cdots,n$}\}$, with ${\cal E}^i:= \{\nu^i\in \Pi^i: \int_{\bar\mathbb{R}^{n^i}} |\bar x^i |^{p'} \nu^i(d\bar x^i) \le  \overline C^i\}$, and note that ${\cal E}$ is  ${\cal W}_p$-compact.

\vspace{3mm}
\noindent{\it Step 2.}\quad We are going to show that the graph of $\bar\pi \in {\cal E} \mapsto {\cal R}(\bar\pi) $ is ${\cal W}_p$-closed,  i.e. given $\bar\pi_m, \bar\pi_\infty\in \Pi\cap {\cal E}$,  $\bar\pi'_m\in  {\cal R}(\bar\pi_m)$ such that $\bar\pi_m\rightarrow\bar\pi_\infty$ in ${\cal W}_p$ and $\bar\pi'_m\rightarrow\bar\pi'_\infty\in  \Pi$, we want to show that $\bar\pi'_\infty\in {\cal R}(\bar\pi_\infty)$.

Denote the concatenation of two probability measures $\nu^i\in \Pi^i, \mu^{-i}\in \Pi^{-i}$ by
\begin{equation*}
\nu^i\otimes \mu^{-i} (dx, d{\mathsf y}) = \nu^i(x^i|{\mathsf y}) \mu^{-i}(x^{-i}|{\mathsf y})dx {\mathsf m} (d{\mathsf y}).
\end{equation*}
Note that for $\bar\pi, \bar\pi'\in {\cal E}$ we have $\bar\pi'^{i}\otimes \bar\pi^{-i}\in {\cal E}$. Since ${\cal E}$ is ${\cal W}_p$-compact, there is a subsequence, still denoted by $(\bar\pi_m ,\bar\pi'_m)$, and $ \bar\pi^*\in \Pi\cap{\cal E}$ such that $\bar\pi'^{i}_m\otimes \bar\pi^{-i}_m \rightarrow \bar\pi^*$  in ${\cal W}_p$ and $\bar\pi^{*,i} = \bar\pi'^{i}_\infty$, $\bar\pi^{*,-i} = \bar\pi^{-i}_\infty$. By the lower-semicontinuity of the mapping: $\bar\pi\mapsto V^i(\bar\pi^i, \bar\pi^{-i})$, we have
\begin{equation}\label{eq:estimatelow}
V^i(\bar\pi'^i_\infty, \bar\pi^{-i}_\infty) = V^i(\bar\pi^{*, i}, \bar\pi^{*,-i}) \leq \liminf_{m\rightarrow\infty} V^i (\bar\pi'^i_m, \bar\pi^{-i}_m).
\end{equation}
Further, fix $\nu^i\in \Pi^i\cap{\cal E}^i$.
 Again by the compactness of ${\cal E}$, there is a subsequence, still denoted by $(\nu^{i} ,\bar\pi'_m)$, and $ \bar\pi^\nu\in \Pi\cap{\cal E}$ such that $\nu^{i}\otimes \bar\pi^{-i}_m \rightarrow \bar\pi^\nu$  in ${\cal W}_p$ and $\bar\pi^{\nu,i} = \nu^i$, $\bar\pi^{\nu,-i} = \bar\pi^{-i}_\infty$. Therefore
\begin{equation*}
\liminf_{m\rightarrow\infty} V^i (\bar\pi'^i_m, \bar\pi^{-i}_m)
\le \liminf_{m\rightarrow\infty} V^i (\nu^{i}, \bar\pi^{-i}_m)
= V^i (\bar\pi^{\nu, i}, \bar\pi^{\nu,-i})
= V^i (\nu^i, \bar\pi^{-i}_\infty).
\end{equation*}
Together with \cref{eq:estimatelow}, we conclude that $\bar\pi^{'i}_\infty \in {\cal R}^i(\bar\pi_\infty)$ for all $i$, and thus $\bar\pi'_\infty\in {\cal R}(\bar\pi_\infty)$.

\vspace{3mm}
\noindent{\it Step 3.}\quad From the condition of the theorem and the result of {\it Step 1\&2}, we conclude that for any $\bar\pi\in \Pi$ the set ${\cal R}(\bar\pi)$ is non-empty, convex, that the set $\cup_{\bar\pi\in \Pi} {\cal R}(\bar\pi) $ is a subset of a ${\cal W}_p$-compact set, and that the graph of the mapping $\bar\pi\in {\cal E} \mapsto {\cal R}(\bar\pi) \subset {\cal E}$ is ${\cal W}_p$-closed. Therefore, it follows from the Kakutani fixed-point theorem that the mapping $\bar\pi \mapsto {\cal R}(\bar\pi)$ has a fixed point, which is a Nash equilibrium.
\end{proof}

\subsection{Invariant measure of MFL system}

\begin{proof}[Proof of \cref{thm:well}] \quad
The proof is based on the Banach fixed point theorem.  Given fixed $(\bar\pi_t), (\bar\pi'_t) \in \overline C_p([0,T], \Pi) $, the SDEs \cref{eq:MFL} apparently have unique strong solutions, denoted by $\bar X, \bar X'$ respectively. Denote by
 $\hat\bar\pi_t: = {\rm Law}(\bar X_t)$ and $\hat\bar\pi'_t: = {\rm Law}(\bar X'_t)$. We are going to show that the mapping $\Psi: \bar\pi\mapsto \hat\bar\pi$ is a contraction as $T$ is small enough.

First, by the standard SDE estimate we know $\mathbb{E}\big[\sup_{t\in [0,T]} |X_t|^p \big]<\infty$, and this implies $(\hat\bar\pi_t), (\hat\bar\pi'_t) \in \overline C_p([0,T], \Pi) $. Note the SDEs for $\tilde X, \tilde X'$ share the same Brownian motion.  Denoting $\delta X = X -X'$, we have
\begin{eqnarray*}
 |\delta X^i_s|^p &=& \left|\int_0^s \left( - \nabla_{x^i} \frac{\delta F^i}{\delta \nu} (\bar\pi^i_t , \bar\pi^{-i}_t, X^i_t, Y)    + \nabla_{x^i} \frac{\delta F^i}{\delta \nu} (\bar\pi'^i_t , \bar\pi'^{-i}_t, X'^i_t, Y)dt \right)dt \right|^p\\
 &\le & T^{1-\frac1p} \Big(\int_0^s C\Big({\cal W}_p^p(\bar\pi^i_t, \bar\pi'^i_t)+  {\cal W}_p^p(\bar\pi^{-i}_t, \bar\pi'^{-i}_t) + |\delta X^i_{t'}|^p\Big)dt +\int_0^s C_0  {\cal W}_p^p(\bar\pi^{-i}_t(\cdot|Y), \bar\pi'^{-i}_t(\cdot|Y)) dt \Big).
\end{eqnarray*}
Taking expectation on both sides, by the Gronwall inequality we obtain that
\begin{equation*}
\mathbb{E}\big[|\delta X^{\mathsf y}_s|^p \big] \le C e^{CT}\int_0^T \Big( {\cal W}^p_p(\bar\pi_t, \bar\pi'_t) + {\cal W}^p_p(\bar\pi_t(\cdot|{\mathsf y}), \bar\pi'_t(\cdot|{\mathsf y}))\Big)dt.
\end{equation*}
 Note that
$\mathbb{E}\big[ |\delta X^{\mathsf y}_s|^p \big]
\ge{\cal W}_p^p (\hat\bar\pi_s(\cdot|{\mathsf y}), \hat\bar\pi'_s(\cdot|{\mathsf y}))$.
Therefore
\begin{equation*}
{\cal W}_p^p (\hat\bar\pi_s(\cdot|{\mathsf y}), \hat\bar\pi'_s(\cdot|{\mathsf y})) \le   C e^{CT}\int_0^T \Big( {\cal W}^p_p(\bar\pi_t, \bar\pi'_t) + {\cal W}^p_p(\bar\pi_t(\cdot|{\mathsf y}), \bar\pi'_t(\cdot|{\mathsf y}))\Big)dt.
\end{equation*}
Integrating both sides with respect to ${\mathsf m}$, we obtain
\begin{equation*}
\overline{\cal W}_p^p (\hat\bar\pi_t, \hat\bar\pi'_t) \le   C e^{CT}\int_0^T 2\overline {\cal W}^p_p(\bar\pi_t, \bar\pi'_t) dt
\le 2CT e^{CT}\sup_{t\in [0,T]} \overline {\cal W}^p_p(\bar\pi_t, \bar\pi'_t) ,
\end{equation*}
 and $\Psi$ is a contraction whenever $T$ is small enough. In case $C_0=0$ the result can be deduced similarly.
\end{proof}

\vspace{3mm}

In order to prove \cref{thm:contraction}, the main ingredient is the reflection coupling  in Eberle \cite{eberle11}. For this mean-field system, we shall adopt the reflection-synchronous coupling as in \cite{eberle2019quantitative}.

We first recall the reflection-synchronous coupling. Fix a parameter $\varepsilon>0$. Introduce the Lipschitz functions ${\rm rc}: \mathbb{R}^N \times \mathbb{R}^N\rightarrow [0,1]$ and ${\rm sc}: \mathbb{R}^N \times \mathbb{R}^N\rightarrow [0,1]$ satisfying
\begin{eqnarray*}
{\rm sc}^2 (x, x') + {\rm rc}^2(x,x') =1, \quad
{\rm rc}(x,x') =1 ~~\mbox{for $|x-x'|\ge \varepsilon$}, \quad
{\rm rc}(x,x') =0 ~~\mbox{for $|x-x'|\le \varepsilon/2$}.
\end{eqnarray*}
 Let $\bar\pi_0,\bar\pi'_0\in \Pi\cap {\cal P}_q(\bar\mathbb{R}^N)$ with some $q>1$ be two initial distributions of the MFL system \cref{eq:MFL}, and $W, \tilde W$ be two independent $N$-dimensional Brownian motions. It follows from \cref{thm:well} the two MFL systems have strong solutions, and denote the marginal laws by $\bar\pi_t, \bar\pi'_t$. Denote the drift of the dynamics \cref{eq:MFL_sys} by
 \begin{equation}\label{eq:drift}
 b^{\mathsf y}(t,x) := \left(- \nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\bar\pi^i_t, \bar\pi^{-i}_t, x^i, {\mathsf y})\right)_{i=1,\cdots,n},
 \quad \tilde b^{\mathsf y}(t,x) := \left(- \nabla_{x^i} \frac{\delta F^i}{\delta \nu}(\bar\pi'^i_t, \bar\pi'^{-i}_t, x^i, {\mathsf y})\right)_{i=1,\cdots,n}
 \end{equation}
 Further, for a fixed ${\mathsf y}\in \mathbb{Y}$,  define the coupling $\Sigma^{\mathsf y}=(X^{\mathsf y}, X'^{\mathsf y})$ as the strong solution to the  SDE
 \begin{eqnarray*}
 && dX^{\mathsf y}_t = b^{\mathsf y}(t, X^{\mathsf y}_t ) dt + {\rm rc} (\Sigma^{\mathsf y}_t) \sigma dW_t + {\rm sc}(\Sigma^{\mathsf y}_t)\sigma d\tilde W_t,\\
 && dX'^{\mathsf y}_t  = \tilde b^{\mathsf y}(t, X'^{\mathsf y}_t) dt + {\rm rc} (\Sigma^{\mathsf y}_t) \big({\rm Id} - 2 e_t \langle e_t, \cdot\rangle \big)\sigma dW_t + {\rm sc}(\Sigma^{\mathsf y}_t)\sigma d\tilde W_t,
 \end{eqnarray*}
where $e_t:= \frac{X^{\mathsf y}_t -X'^{\mathsf y}_t}{|X^{\mathsf y}_t -X'^{\mathsf y}_t| }$ for $X^{\mathsf y}_t \neq X'^{\mathsf y}_t $, otherwise $e_t:=\hat e$ some arbitrary fixed unit vector in $\mathbb{R}^N$.

In the following proof, we will use the concave increasing function $f$ constructed in  \cite[Theorem 2.3]{eberle2019quantitative}:
\begin{equation*}
	f(r) := \int_{0}^{r}\varphi(s)g(s\wedge R_2) d s,\quad
	\mbox{where}\quad g(r):= 1-\frac{c}{2}\int_0^r \Phi(s)\varphi(s)^{-1}ds,
\end{equation*}
and the function $ \varphi  $ and  the constants $R_2, c$ are defined  in the statement of \cref{thm:contraction}.
In particular, on  $(0,R_2)\cup(R_2,+\infty)$ the function $f$ satisfies
\begin{eqnarray}
& 2\sigma^2 f''(r)  \le -r \big(\kappa(r)+2\gamma \eta(r)\big) f'(r) - c\sigma^2 f(r) \label{eq:fODE}
\end{eqnarray}
and for $ r\in\mathbb{R}^+ $
\begin{eqnarray}
& r\varphi(R_1) \le \Phi (r) \le 2 f(r) \le 2\Phi (r) \le 2r. \label{eq:bound_f}
\end{eqnarray}

\vspace{3mm}
\begin{proof}[Proof of \cref{thm:contraction}] \quad Note that $\pi_0(\cdot|{\mathsf y}),\pi'_0(\cdot|{\mathsf y}) \in {\cal P}_q(\mathbb{R}^N)$ with $q>1$ for ${\mathsf m}$-a.s. ${\mathsf y}\in \mathbb{Y}$. For such ${\mathsf y}$, we may choose the coupling $ ( X^{\mathsf y}_0, X'^{\mathsf y}_0 ) $ so that
\begin{eqnarray}\label{eq:fWlower}
{\cal W}_1 \big(\pi_0(\cdot|{\mathsf y}),\pi'_0(\cdot|{\mathsf y}) \big) = \mathbb{E} \big[|X^{\mathsf y}_0 - X'^{\mathsf y}_0|\big]
\ge \mathbb{E} \Big[f\big(|X^{\mathsf y}_0-X'^{\mathsf y}_0|\big)\Big]  .
\end{eqnarray}
The last inequality is due to \eqref{eq:bound_f}. On the other hand, for all $t\ge 0$ we have
\begin{eqnarray}\label{eq:fWupper}
{\cal W}_1 \big(\pi_t(\cdot|{\mathsf y}),\pi'_t(\cdot|{\mathsf y}) \big) \le  \mathbb{E} \big[|X^{\mathsf y}_t - X'^{\mathsf y}_t|\big]
\le  \frac{2}{\varphi(R_1)}\mathbb{E} \Big[f\big(|X^{\mathsf y}_t - X'^{\mathsf y}_t|\big)\Big].
\end{eqnarray}
Denote $\delta X^{\mathsf y}_t : = X^{\mathsf y}_t - X'^{\mathsf y}_t $. By the definition of the coupling above, we have
\begin{eqnarray*}
d \delta X^{\mathsf y}_t = \Big(b^{\mathsf y}(t, X^{\mathsf y}_t) - \tilde b^{\mathsf y}(t, X'^{\mathsf y}_t)\Big) dt + 2 {\rm rc}(\Sigma^{\mathsf y}_t) \sigma e_t d\bar W_t,
\end{eqnarray*}
where $\bar W_t:=\int_0^t  e_s \cdot dW_s $ is a one-dimensional Brownian motion. Denote $r_t := |\delta X^{\mathsf y}_t|$ and note that by the definition of ${\rm rc}$ we have ${\rm rc} (\Sigma^{\mathsf y}_t) = 0$ whenever $r_t\le \varepsilon/2$. Therefore,
\begin{eqnarray*}
d r_t = e_t \cdot \Big( b^{\mathsf y}(t, X^{\mathsf y}_t)- \tilde b^{\mathsf y}(t, X'^{\mathsf y}_t)\Big) dt + 2 {\rm rc}(\Sigma^{\mathsf y}_t) \sigma d \bar W_t.
\end{eqnarray*}
Then it follows from the It\^o-Tanaka formula and the concavity of $f$ that
\begin{eqnarray*}
 f(r_t) - f(r_0)
&\le & \int_0^t \left( f'(r_s) e_s \cdot \big( b^{\mathsf y}(s,X^{\mathsf y}_s)- \tilde b^{\mathsf y}(s, X'^{\mathsf y}_s)\big) + 2{\rm rc}(\Sigma^{\mathsf y}_s)^2 \sigma^2 f''(r_s) \right) ds
+ M_t,
\end{eqnarray*}
where $ M_t:= 2 \int_0^t {\rm rc}(\Sigma^{\mathsf y}_s) f' (r_s) \sigma d\bar W_s $ is a martingale. Now  note that
\begin{eqnarray*}
  e_s \cdot \big( b^{\mathsf y}(s,X^{\mathsf y}_s)- \tilde b^{\mathsf y}(s, X'^{\mathsf y}_s)\big)
 \le 1_{\{ r_s\ge \varepsilon\}} r_s \kappa(r_s) + 1_{\{r_s<\varepsilon\}} C \varepsilon + \gamma \big(\overline{\cal W}_1(\bar\pi_s, \bar\pi'_s) + {\cal W}_1(\pi_s(\cdot|{\mathsf y}), \pi'_s(\cdot|{\mathsf y})) \big).
\end{eqnarray*}
Further, since $f'' \le 0$ and ${\rm rc}(\Sigma^{\mathsf y}_s) = 1$ whenever $r_s \ge\varepsilon$, we have
\begin{eqnarray*}
f(r_t) - f(r_0)
&\le& \int_0^t \Big( 1_{\{ r_s\ge \varepsilon\}}\big( f'(r_s) r_s \kappa(r_s)  + 2 \sigma^2 f''(r_s)\big)  +  1_{\{ r_s < \varepsilon\}} C \varepsilon \\
        && \quad\q\quad\q\quad\q\quad + \gamma f'(r_s) \big(\overline{\cal W}_1(\bar\pi_s, \bar\pi'_s) + {\cal W}_1(\pi_s(\cdot|{\mathsf y}), \pi'_s(\cdot|{\mathsf y})) \big)\Big) ds  + M_t.
\end{eqnarray*}
By taking expectation on both sides, we obtain
\begin{eqnarray*}
&&\int \mathbb{E} \Big[f(r_t) - f(r_0)\Big] {\mathsf m}(d{\mathsf y}) \\
&\le& \int_0^t \Big( \int \mathbb{E}\Big[ 1_{\{ r_s\ge \varepsilon\}}\big( f'(r_s) r_s \kappa(r_s)  + 2 \sigma^2 f''(r_s)\big) +  C \varepsilon  \Big] {\mathsf m}(d{\mathsf y})  + 2 \gamma  f'(r_s)\overline {\cal W}_1(\bar\pi_s, \bar\pi'_s) \Big) ds\\
&\le & \int_0^t \Big( \int  \mathbb{E}\Big[  -c \sigma^2 f(r_s)1_{\{ r_s\ge \varepsilon\}}  +  C \varepsilon   \Big]  {\mathsf m}(d{\mathsf y}) + 2 \gamma  \overline {\cal W}_1(\bar\pi_s, \bar\pi'_s)  \Big)ds\\
&\le & \int_0^t \Big( \int  \mathbb{E}\Big[  -c \sigma^2 f(r_s) +  (C + c\sigma^2) \varepsilon   \Big]  {\mathsf m}(d{\mathsf y}) + 2 \gamma  \overline {\cal W}_1(\bar\pi_s, \bar\pi'_s)  \Big)ds.
\end{eqnarray*}
The second  last inequality is due to \eqref{eq:fODE}.
Together with  \eqref{eq:fWlower} and \eqref{eq:fWupper}, we obtain
\begin{eqnarray*}
\frac{\varphi(R_1)}{2}e^{c\sigma^2 t} \overline{\cal W}_1 \big(\bar\pi_t,\bar\pi'_t \big) - \overline{\cal W}_1\big(\bar\pi_0,\bar\pi'_0\big)
&\le & e^{c\sigma^2t}\int \mathbb{E} \big[f(r_t) \big] {\mathsf m}(d{\mathsf y})-\int \mathbb{E} \big[ f(r_0)\big] {\mathsf m}(d{\mathsf y})\\
&\le& \int_0^t e^{c\sigma^2 s} \big( ( C + c \sigma^2 ) \varepsilon  +2 \gamma  \overline {\cal W}_1(\bar\pi_s, \bar\pi'_s) \big)ds .
\end{eqnarray*}
This holds true for all $\varepsilon>0$, so finally we obtain \cref{eq:contractionresult} by Gronwall's inequality.
\end{proof}

\subsection{One player case}

Throughout this subsection, we suppose that the assumptions in \cref{thm:GD} hold true. Recall the drift function $b^{\mathsf y}$ defined in \cref{eq:drift} with $n=1$, i.e.
\begin{equation*}
b^{\mathsf y}(t,x):=- \nabla_{x} \frac{\delta F}{\delta \nu}(\bar\pi_t, x, {\mathsf y}).
\end{equation*}
Under the assumption of \cref{thm:GD} the function $b^{\mathsf y}$ is continuous in $(t,x)$ and $C^3$ in $x$ for all $t\in [0,T]$. Due to a classical regularity result in the theory of linear PDEs (see e.g. \cite[p.14-15]{JKO98}), we obtain the following result.
\begin{lemma}\label{lem:FP}
Under the assumption of \cref{thm:GD},  the marginal laws $(\pi_t(\cdot | {\mathsf y}))_{t\ge 0}$ of the solution to \cref{eq:MFL_sys}  are weakly continuous solutions  to the Fokker-Planck equations:
\begin{eqnarray}\label{eq:FoP}
\partial_t \nu = \nabla_x \cdot (-b^{\mathsf y} \nu  + \frac{\sigma^2}{2} \nabla_x \nu) \quad\mbox{for}\quad {\mathsf y}\in \mathbb{Y}.
\end{eqnarray}
In particular, we have that $(t,x)\mapsto \pi_t(x|{\mathsf y})$ belongs to $C^{1,2}\big((0,\infty)\times \mathbb{R}^N)\big)$.
\end{lemma}

The following results can be proved with the same argument as in Lemma 5.5-5.7 in \cite{HKR19}, so the proof is omitted.

\begin{lemma}\label{lem:finit_ent}
 Fix a ${\mathsf y}\in \mathbb{Y}$ and assume $ \mathbb{E}\big[ |X^{\mathsf y}_0|^2 \big] < \infty $, where we recall the $X^{\mathsf y}$ defined in \cref{eq:MFL_sys}. Denote by $ \mathbb{Q}^{\sigma}_{\mathsf y}$ the scaled Wiener measure\footnote{Let $B$ be the canonical process of the Wiener space and $\mathbb{Q}$ be the Wiener measure, then the scaled Wiener measure $ \mathbb{Q}^{\sigma}: = \mathbb{Q}\circ (\sigma B)^{-1} $.} with initial distribution $ \pi_0(\cdot|{\mathsf y})$ and by $({\cal F}_t)_{t\ge 0}$ the canonical filtration of the Wiener space. Then
 \begin{enumerate}
 	\item[{\rm (i)}] for any finite horizon $T>0$, the law of the solution to \cref{eq:MFL}, $ \mathbf{\Pi}(\cdot|{\mathsf y}):={\rm Law}(X^{\mathsf y})$, is equivalent to $ \mathbb{Q}^{\sigma}_{\mathsf y}$  on $ {\cal F}_T $ and the relative entropy
 	\begin{equation*}
 		\int \ln \Big( \frac{d\mathbf{\Pi}(\cdot|{\mathsf y})}{d\mathbb{Q}^{\sigma}_{\mathsf y}}\Big|_{{\cal F}_T} \Big) d \nu_t
		= \mathbb{E}\Big[ \int_0^T \big| b^{\mathsf y}(t, X^{\mathsf y}_t) \big|^2 dt \Big]
		<+\infty.
 	\end{equation*}

 	\item[{\rm (ii)}] the marginal law $ \pi_t(\cdot|{\mathsf y})$ admits a density s.t. $\pi_t(\cdot|{\mathsf y})>0 $ and $ H\big(\pi_t(\cdot|{\mathsf y})\big)<+\infty $ for $t>0$.

	\item[{\rm (iii)}]  the function $\ln \pi_t(x|{\mathsf y})$ is continuous differentiable in $x$ for $t>0$, and for any $t_0\in (0,t]$ it satisfies
	\begin{equation*}
 			\nabla_x \ln \pi_t(x|{\mathsf y}) = -\frac1{t_0}\mathbb{E}\left[ \int_{0}^{t_0}\big(1-s \nabla_x b^{\mathsf y}(s,X^{\mathsf y}_{t-t_0+s})\big)dW^{t-t_0}_s \Big| X^{\mathsf y}_t=x \right],
 		\end{equation*}
 		where $ W^{t-t_0}_s := W_{t-t_0+s} - W_{t-t_0} $ and $W$ is the Brownian motion in \cref{eq:MFL_sys}.
In particular,  for any $ t_*>0 $ we have
\begin{eqnarray*}
C:=\sup_{s\ge t_*}\int_{\mathbb{R}^N } \big| \nabla_x \ln \pi_s(x|{\mathsf y})  \big|^2 \pi_s(x|{\mathsf y})  dx <+\infty,
\end{eqnarray*}
and $C$ only depends on $t_*$  and the Lipschitz constant of $\nabla_x b^{\mathsf y}$ with respect to $ x $.

	\item[{\rm (iv)}] we have the estimates
\begin{eqnarray*}
		&\int_{\mathbb{R}^N} |\nabla_x \ln \pi_t(x|{\mathsf y}) | dx <+\infty, \quad \int_{\mathbb{R}^N}|x\cdot\nabla_x \ln \pi_t(x|{\mathsf y})| dx<+\infty\quad\mbox{for all}\quad t>0,\\
		&\mbox{and}\quad \int_t^{t'} \int_{\mathbb{R}^N}|\Delta_{xx} \pi_s(x|{\mathsf y})| dx ds <+\infty\quad\mbox{for all}\quad t'>t>0,
\end{eqnarray*}
	and together with the integration by parts we obtain for all $t'>t>0$
\begin{eqnarray}\label{lem:IPP}
		&\int_{\mathbb{R}^N} \Delta_{xx}  \frac{\delta F}{\delta \nu}(\bar\pi_t, x, {\mathsf y}) \pi_t(x|{\mathsf y}) dx
		= -\int_{\mathbb{R}^N} \nabla_{x} \frac{\delta F}{\delta \nu}(\bar\pi_t, x, {\mathsf y}) \cdot \nabla_x \pi_t(x|{\mathsf y})dx,\\
		&\int_t^{t'}\int_{\mathbb{R}^N}\Delta_{xx} \big(\ln \pi_s(x|{\mathsf y}) \big)\pi_s(x|{\mathsf y}) dxds = -\int_t^{t'} \int_{\mathbb{R}^N}\left|  \nabla_x \ln \pi_s(x|{\mathsf y}) \right|^2 \pi_s(x|{\mathsf y}) dxds.\notag
\end{eqnarray}
 \end{enumerate}
 \end{lemma}

 The proof of \cref{thm:GD} is based on the previous lemma and  It\^o calculus.

 \vspace{3mm}

\begin{proof}[Proof of \cref{thm:GD}] \quad
It follows from  the It\^o-type formula \cite[Theorem 4.14]{Carmona+Rene} and  \eqref{lem:IPP} that
\begin{eqnarray}
\mathrm{d}F(\bar\pi_t)
&=& \int_{\mathbb{R}^N} \left(- \big|\nabla_x \frac{\delta F}{\delta \nu}(\bar\pi_t,\bar x) \big|^2 + \frac{\sigma^2}{2}{\rm Tr}( \nabla^2_x \frac{\delta F}{\delta \nu}(\bar\pi_t,\bar x))\right) \bar\pi_t (d\bar x) d t \nonumber \\
&=& \int_{\mathbb{R}^N} \left(-|\nabla_x \frac{\delta F}{\delta \nu}(\bar\pi_t, \bar x)|^2 - \frac{\sigma^2}{2}\nabla_x \frac{\delta F}{\delta \nu}(\bar\pi_t,\bar x)\cdot \frac{\nabla_x \pi_t(x|{\mathsf y})}{\pi_t(x|{\mathsf y})}\right) \bar\pi_t(dx)d t. \label{FreeEnergyEq01}
\end{eqnarray}
On the other hand, by It\^o's formula and  the Fokker-Planck equation \eqref{eq:FoP}, we have
\begin{equation*}
 \mathrm{d} \log\pi_t(X_t|{\mathsf y}) = \left(\sigma^2\nabla_x \cdot \frac{\nabla_x \pi_t(X_t|{\mathsf y})}{\pi_t(X_t|{\mathsf y})} + \frac{\sigma^2}{2}\left| \frac{\nabla_x \pi_t(X_t|{\mathsf y})}{\pi_t(X_t|{\mathsf y})}\right|^2 - \nabla_x\cdot \nabla_x\frac{\delta F}{\delta\nu}(\bar\pi_t, X_t, {\mathsf y}) \right) dt+ \mathrm{d} M_t .
\end{equation*}
where $M$ is a martingale on $[t,T]$ for any $0<t<T$. By taking expectation on both sides and using  \eqref{lem:IPP}, we obtain for $t>0$:
\begin{eqnarray}	\label{FreeEnergyEq02}
\mathrm{d}H(m_t)
&=&\mathbb{E} \left[ - \frac{\sigma^2}{2}\left|\frac{\nabla_x \pi_t(X_t|{\mathsf y})}{\pi_t(X_t|{\mathsf y})} \right|^2 +  \nabla_x\frac{\delta F}{\delta\nu}(\bar\pi_t, X_t, {\mathsf y})\cdot \frac{\nabla_x \pi_t(X_t|{\mathsf y})}{\pi_t(X_t|{\mathsf y})} \right]\mathrm{d}t \notag \\
&=& \int_{\mathbb{R}^N} \left( - \frac{\sigma^2}{2}\left|  \frac{\nabla_x \pi_t(x|{\mathsf y})}{\pi_t(x|{\mathsf y})} \right|^2 - \nabla_x\frac{\delta F}{\delta\nu}(\bar\pi_t, x, {\mathsf y}) \cdot \left( \frac{\nabla_x \pi_t(x|{\mathsf y})}{\pi_t(x|{\mathsf y})}\right) \right)\bar\pi(\bar x) dt.
\end{eqnarray}
Summing up equations \eqref{FreeEnergyEq01} and \eqref{FreeEnergyEq02}, we obtain \cref{eq:GD}.
\end{proof}

 \vspace{3mm}

   Let $ (\bar\pi_t)_{t\ge 0} $ be the flow of marginal laws of the solution of \cref{eq:MFL} given an initial law $ \bar\pi_0$. Define a dynamic system $ {\cal S}(t)\left[\bar\pi_0\right] := \bar\pi_t $.
Due to the result of \cref{thm:GD}, we can view the function $V$ as a Lyapunov function of the dynamic system $({\cal S}_t)_{t\ge 0}$, and then it is natural to prove \cref{thm:invariant} using LaSalle's invariance principle (see the following \cref{prop:Lasalle}). However, $V^\sigma$ is not continuous (only lower-semicontinuous), in particular, the mapping $t\mapsto V^\sigma(\bar\pi_t)$ is not a priori continuous at $+\infty$, which makes the proof non-trivial. Here we follow the strategy first developed in \cite{HRSS19} to overcome the difficulty.
 Define the $ \omega $-limit set:
\begin{equation*}
	\omega(\bar\pi_0) := \Big\{ \bar\pi\in\Pi:\quad \text{there exists }t_n\rightarrow+\infty\text{ such that } {\cal W}_2\left({\cal S}(t_n)\left[\bar\pi_0\right], \bar\pi\right)\rightarrow 0  \Big\}.
\end{equation*}

\begin{lemma}
Assume that   $ \bar\pi^0\in \Pi\cap{\cal P}_q(\bar\mathbb{R}^N)$ with $q\ge 2$. Then for the solution  to the MFL system \cref{eq:MFL} we have
\begin{equation}\label{eq:uniformmoment}
\sup_{t\ge 0} \int_{\bar\mathbb{R}^N} |\bar x|^q  \bar\pi_t(d \bar x) <\infty.
\end{equation}
\end{lemma}
\begin{proof}
By It\^o formula, we obtain
\begin{equation*}
d |X^{\mathsf y}_t|^q = q | X^{\mathsf y}_t |^{q-2}\left(- 2 X^{\mathsf y}_t \cdot \frac{\delta F}{\delta \nu}(\bar\pi_t, X^{\mathsf y}_t , {\mathsf y} ) +\sigma^2q(q-2+N) \right) dt + \sigma q|X^{\mathsf y}_t|^{q-2}X^{\mathsf y}_t \cdot dW_t.
\end{equation*}
By the linear growth assumption on $\frac{\delta F}{\delta \nu}$ and the dissipative condition \cref{assum:dissipative}, there is a constant $M$ such that
\begin{eqnarray*}
d |X^{\mathsf y}_t|^q  &\le & q | X^{\mathsf y}_t |^{q-2}\left( C - \varepsilon |X^{\mathsf y}_t|^2 1_{\{|X^{\mathsf y}_t|\ge M\}} \right) dt + \sigma q|X^{\mathsf y}_t|^{q-2}X^{\mathsf y}_t \cdot dW_t\\
			&\le & q | X^{\mathsf y}_t |^{q-2}\left( (C+\varepsilon M^2) - \varepsilon |X^{\mathsf y}_t|^2  \right) dt + \sigma q|X^{\mathsf y}_t|^{q-2}X^{\mathsf y}_t \cdot dW_t.
\end{eqnarray*}
In particular, the constant $C$ above does not depend on ${\mathsf y}$. In case $q =2$, by taking expectation on both sides and using the Gronwall inequality, we obtain
\begin{equation}\label{eq:tight}
\sup_{t\ge 0} \mathbb{E}[|X^{\mathsf y}_t|^q] \le C\big(1+  \mathbb{E}[|X^{\mathsf y}_0|^q]\big),
\end{equation}
with $q=2$.
For $q>2$, a similar inequality follows from the induction. The desired result \cref{eq:uniformmoment} follows from integrating with respect to ${\mathsf m}$ on both sides of the inequality.
\end{proof}

\begin{proposition}[Invariance Principle]\label{prop:Lasalle}
	Assume that   $ \bar\pi^0\in \Pi\cap{\cal P}_q(\bar\mathbb{R}^N)$ with $q>2$. Then the set $ \omega(\bar\pi_0) $ is non-empty, compact and invariant, that is
	\begin{enumerate}
		\item[{\rm (a)}] for any $ \bar\pi\in\omega(\bar\pi_0) $, we have $ S(t)\left[\bar\pi\right]\in\omega(\bar\pi_0) $ for all $ t\ge 0 $;

		\item[{\rm (b)}] for any $ \bar\pi\in\omega(\bar\pi_0) $ and  $t\ge 0$, there exists $ \bar\pi'\in\omega(\bar\pi_0) $ such that $ S(t)\left[\bar\pi'\right]=\bar\pi $.
	\end{enumerate}
\end{proposition}

\begin{proof}
Note that
\begin{itemize}
\item the mapping $\bar\pi_0 \mapsto  {\cal S}(t)\left[\bar\pi_0\right]$ is ${\cal W}_2$-continuous due to the stability on the initial law;
\item the mapping $t\mapsto  {\cal S}(t)\left[\bar\pi_0\right]$ belongs to $C_2\big(\mathbb{R}^+, \Pi\big)$, due to \cref{thm:well};
\item the set $\big\{ {\cal S}(t)[\bar\pi_0],~t\ge 0 \big\}$ belongs to a ${\cal W}_2$-compact set, due to \cref{eq:uniformmoment}.
\end{itemize}
The rest follows the standard argument for Lasalle's invariance principle (see e.g. \cite[Theorem 4.3.3]{Henry81} or \cite[Proposition 6.5]{HRSS19}).
\end{proof}

\vspace{3mm}
\begin{proof}[Proof of \cref{thm:invariant}] \quad
{\it Step 1.}\quad Recall the set ${\cal I}$ defined in \cref{eq:limitset}. We first prove the existence of a converging subsequence towards a  member in ${\cal I}$. Since $\omega (\bar\pi_0)$ is ${\cal W}_2 $-compact and $V$ is ${\cal W}_2$-lower semicontinuous, there is $\bar\pi^*\in \argmin_{\bar\pi\in \omega(\bar\pi_0)} V(\bar\pi)$. By the backward invariance (b), given $t>0$ there is $\nu\in \omega(\bar\pi_0)$ such that ${\cal S} (t)[\nu] = \bar\pi^*$. By \cref{thm:GD}, we have
\begin{equation*}
V\big({\cal S}(t+s)[\nu]\big) \le V\big({\cal S}(t)[\nu]\big) = V(\bar\pi^*), \quad\mbox{for all $s>0$}.
\end{equation*}
Further by the forward invariance (a), we know ${\cal S}(t+s)[\nu] \in \omega(\bar\pi_0)$, and by the optimality of $\bar\pi^*$ we obtain $V\big({\cal S}(t+s)[\nu]\big) = V(\bar\pi^*)$. Again by \cref{thm:GD}, we get
\begin{equation*}
0= \frac{d V\big({\cal S}(t)[\nu]\big) }{dt}
=-  \int_{\bar\mathbb{R}^N} \left| \nabla_{x} \frac{\delta F}{\delta \nu}(\bar\pi^*,  x,{\mathsf y}) + \frac{\sigma^2}{2} \nabla_{x}\ln\big(\pi^*(x |{\mathsf y})\big) \right|^2  \bar\pi^* (d\bar x).
\end{equation*}
Since $\bar\pi^*={\cal S} (t)[\nu] $ is equivalent to ${\rm Leb}\times {\mathsf m}$ according to \cref{lem:finit_ent} (ii), we have $\bar\pi^*\in {\cal I}$. By the definition of $\omega(\bar\pi_0)$, there is a subsequence of $(\bar\pi_t)$ converging towards $\bar\pi^*$.

\vspace{3mm}

\noindent{\it Step 2 (a).}\quad We first prove the result under the assumption (ii.a). Let $(\bar\pi_{t_n})_n$ be a sequence converging to $\bar\pi^*$ in ${\cal W}_2$.
Due to the estimate \cref{eq:tight} and the fact that $\mathbb{Y}$ is countable, there is subsequence, still denoted by $(t_n)_n$ such that for each ${\mathsf y}\in \mathbb{Y}$, $\pi_{t_n}(\cdot|{\mathsf y})$ converges to a probability measure $\pi^{\mathsf y}$ in ${\cal W}_2$.
Then clearly $\pi^{\mathsf y}=\pi^*(\cdot|{\mathsf y})$ for ${\mathsf m}$-a.s. ${\mathsf y}$, and thus $\bar\pi_{t_n}(\cdot|{\mathsf y})\rightarrow \bar\pi^*(\cdot|{\mathsf y})$ in ${\cal W}_2$ for ${\mathsf m}$-a.s. ${\mathsf y}$. Note that
\begin{equation*}
\bar\pi^*(x|{\mathsf y})= C \exp\left( -\frac{2}{\sigma^2}\frac{\delta F}{\delta \nu}(\bar\pi^*, x,{\mathsf y})\right),
\end{equation*}
in particular, $\bar\pi^*(\cdot|{\mathsf y})$ is log-semiconcave. By the HWI inequality (see \cite[Theorem 3]{OV00}) we have
\begin{equation}\label{HWI}
\int  \big(\ln \bar\pi_{t_n}(x|{\mathsf y}) -\ln  \bar\pi^*(x|{\mathsf y}) \big) \bar\pi_{t_n}(dx|{\mathsf y})
\le
{\cal W}_2 \big( \bar\pi_{t_n}(x|{\mathsf y})  , \bar\pi^*(x|{\mathsf y}) \big) \left(\sqrt{{\cal I}_n^{\mathsf y}} + C {\cal W}_2 \big(\bar\pi_{t_n}(x|{\mathsf y})  , \bar\pi^*(x|{\mathsf y}) \big) \right),
\end{equation}
where ${\cal I}^{\mathsf y}_n$ is the relative Fisher information defined as
\begin{eqnarray*}
{\cal I}^{\mathsf y}_n &: =& \int \Big| \nabla_x \ln\bar\pi_{t_n}(x|{\mathsf y})  - \nabla_x \ln\bar\pi^*(x|{\mathsf y})   \Big|^2 \bar\pi_{t_n}(dx|{\mathsf y}) \\
&=&   \int \Big| \nabla_x \ln\bar\pi_{t_n}(x|{\mathsf y})  +\frac{2}{\sigma^2}\nabla_x \frac{\delta F}{\delta \nu}(\bar\pi^*, x,{\mathsf y})  \Big|^2 \bar\pi_{t_n}(dx|{\mathsf y})  \\
&\le & 2 \int \big| \nabla_x \ln \bar\pi_{t_n}(x|{\mathsf y})  \big|^2 \bar\pi_{t_n}(d x|{\mathsf y}) + C \Big(1+\int |x|^2 \bar\pi_{t_n}(d x|{\mathsf y}) \Big),
\end{eqnarray*}
where the last inequality is due to the linear growth of $\nabla_x \frac{\delta F}{\delta \nu}$ in $x$.
It follows from \cref{lem:finit_ent} (iii) that $\sup_{n, {\mathsf y}} {\cal I}^{\mathsf y}_n<\infty$. Integrate both sides of \cref{HWI} with respect to ${\mathsf m}$, and obtain
\begin{equation}\label{HWIest}
H\big(\bar\pi_{t_n} \big|\bar\pi^*\big) \le C \overline{\cal W}_2(\bar\pi_{t_n}, \bar\pi^*)\Big(1+ \overline{\cal W}_2(\bar\pi_{t_n}, \bar\pi^*) \Big).
\end{equation}
The right hand side converges to $0$ as $n\rightarrow\infty$ by the dominated convergence theorem. Therefore,
\begin{eqnarray*}
&& \limsup_{n\rightarrow\infty} H(\bar\pi_{t_n}| {\rm Leb}\times {\mathsf m}) - H(\bar\pi^*| {\rm Leb}\times {\mathsf m})\\
& = & \limsup_{n\rightarrow\infty} \int \ln\left(\frac{\pi_{t_n}(x|{\mathsf y})}{ \pi^*(x|{\mathsf y})}\right) \bar\pi_{t_n}(d\bar x)+ \int \ln\big(\pi^*(x|{\mathsf y})\big) (\bar\pi_{t_n} -\bar\pi^*)(d\bar x) \\
 & =& \limsup_{n\rightarrow\infty}H\big(\bar\pi_{t_n} \big|\bar\pi^*\big)  ~\le~ 0,
\end{eqnarray*}
where the last equality is due to the dominated convergence theorem and the last inequality is due to \cref{HWIest}.
Since $H$ is ${\cal W}_2$-lower-semicontinuous, we have $ \lim\limits_{n\rightarrow\infty} H(\bar\pi_{t_n}| {\rm Leb}\times {\mathsf m}) = H(\bar\pi^*| {\rm Leb}\times {\mathsf m}) $. Together with the fact that $F$ is ${\cal W}_2$-continuous, we have $\lim\limits_{t\rightarrow\infty} V(\bar\pi_t) = V (\bar\pi^*)$. Further by the ${\cal W}_2$-lower-semicontinuity of $V$, we obtain
\begin{equation*}
V(\bar\pi) \le \lim_{t'_n\rightarrow\infty} V(\bar\pi_{t'_n})=V(\bar\pi^*), \quad\mbox{for all $\bar\pi\in \omega (\bar\pi_0)$}.
\end{equation*}
Together with the optimality of $\bar\pi^*$, we have $V(\bar\pi) =V(\bar\pi^*)$ for all $\bar\pi\in \omega (\bar\pi_0)$. Finally by the invariant principle and \cref{eq:GD}, we conclude that $\omega(\bar\pi_0)\subset{\cal I}$.

\vspace{3mm}

\noindent{\it Step 2 (b).}\quad Similarly we can prove the result under the assumption (ii.b). Let $(\bar\pi_{t_n})_n$ be a sequence converging to $\bar\pi^*$ in ${\cal W}_2$. Note that
\begin{equation*}
\bar\pi^*(x, y)= C \exp\left( -\frac{2}{\sigma^2}\frac{\delta F}{\delta \nu}(\bar\pi^*, x,{\mathsf y})\right){\mathsf m}({\mathsf y})
\end{equation*}
is log-semiconcave due to the assumption. Due to the HWI inequality, we have
\begin{equation*}
\int  \big(\ln \bar\pi_{t_n}(x|{\mathsf y}) -\ln  \bar\pi^*(x|{\mathsf y}) \big) \bar\pi_{t_n}(d\bar x)
\le
{\cal W}_2 \big( \bar\pi_{t_n} , \bar\pi^* \big) \left(\sqrt{{\cal I}_n} + C {\cal W}_2 \big(\bar\pi_{t_n} , \bar\pi^* \big) \right),
\end{equation*}
where ${\cal I}_n$ is the relative Fisher information defined as
\begin{eqnarray*}
{\cal I}_n &: =& \int \Big| \nabla_x \ln\bar\pi_{t_n}(x|{\mathsf y})  - \nabla_x \ln\bar\pi^*(x|{\mathsf y})   \Big|^2 \bar\pi_{t_n}(d\bar x) \\
&=&   \int \Big| \nabla_x \ln\bar\pi_{t_n}(x|{\mathsf y})  +\frac{2}{\sigma^2}\nabla_x \frac{\delta F}{\delta \nu}(\bar\pi^*, x,{\mathsf y})  \Big|^2 \bar\pi_{t_n}(d\bar x)  \\
&\le & 2 \int \big| \nabla_x \ln \bar\pi_{t_n}(x|{\mathsf y})  \big|^2 \bar\pi_{t_n}(d \bar x) + C \Big(1+\int |x|^2 \bar\pi_{t_n}(d \bar x) \Big).
\end{eqnarray*}
Again by \cref{lem:finit_ent} (iii) we have $\sup_{n} {\cal I}_n<\infty$. For the rest, we may follow the same lines of arguments in Step 2 (a) to conclude the proof.
\end{proof}

\bibliographystyle{siamplain}