Propagation of chaos for mean field Schrödinger problems[The authors would like to thank Julio Backhoff-Veraguas for helpful discussions. Camilo Hernández acknowledges the support of a Presidential Postdoctoral Fellowship at Princeton University and the Cecilia Tanner Research Fund at Imperial College London. Ludovic Tangpi acknowledges partial support by the NSF under Grants DMS-2005832 and CAREER DMS-2143861.]

Camilo Hernández [Princeton University, ORFE department, USA<EMAIL_ADDRESS>Ludovic Tangpi [Princeton University, ORFE department, USA<EMAIL_ADDRESS>
In this work, we study the mean field Schrödinger problem from a purely probabilistic point of view by exploiting its connection to stochastic control theory for McKean-Vlasov diffusions.
Our main result shows that the mean field Schrödinger problem arises as the limit of “standard” Schrödinger problems over interacting particles.
Due to the stochastic maximum principle and a suitable penalization procedure, the result follows as a consequence of novel (quantitative) propagation of chaos results for forward-backwards particle systems.
The approach described in the paper seems flexible enough to address other questions in the theory.
For instance, our stochastic control technique further allows us to solve the mean field Schrödinger problem and characterize its solution, the mean field Schrödinger bridge, by a forward-backward planning equation.

§ INTRODUCTION

In 1932  [Schrödinger, 1931, Schrödinger, 1932] introduced the problem of determining the most likely path of a cloud of independent scattering particles, conditionally on the observation of their initial, $t=0$, and terminal, $t=T\in (0,\infty)$ distributions.
This was translated in modern probabilistic terms, see  [Föllmer, 1988], as an entropy minimization problem with marginal constraints by employing the theory of large deviations.
Let us present this heuristically.
For positive integers $d$ and $i\geq 1$, let $\Omega^i:=\Cc([0,T],\R^d)$, $X^i$ be the canonical process, $\P_i$ be the Wiener measures under which the particle $X^i$ has initial distribution $\mu_{\rm in}$, and let $\Omega$ (resp. $\P$) be the associated infinite product space (resp. probability measure).
Because $(X^i)_{i\ge1}$ are independent and identically distributed, it then follows from Sanov's theorem that for any probability measure $\Q$ on $\Cc([0,T],\R^d)$,

\[
\P\bigg( \frac{1}N \sum_{i=1}^N \delta_{X^i} =\Q\bigg)\underset{N\to \infty}{\approx} \exp\big(-N\cdot   \Hc(\Q|\P_1)\big),\; \text{where } \Hc(\Q| \P_1):=\begin{cases} \displaystyle\E\bigg[ \frac{\d \Q}{\d \P_1}\log \frac{\d \Q}{\d \P_1}\bigg] & \text{if }\Q\ll\P_1\\
\infty&  \text{otherwise}
\end{cases}
\]

denotes the relative entropy between $\P_1$ and $\Q$.
That is, the most likely behaviour of the process $(X^i)_{i\geq 1}$ is given by the measure $\Q$ which minimizes the relative entropy $\Hc(\cdot | \P_1)$ with respect to the measure $\P_1$, i.e.,

\begin{equation}
\label{eq:dyn.SP}
\Vc_{\rm S}(\muin,\mufin):= \inf \Big\{ \Hc(\Q|\P_1): \Q\in \Pc(\Omega),\, \Q\circ (X_0^1)^{-1}=\muin, \,\Q\circ (X_T^1)^{-1}=\mufin \Big\}.
\end{equation}

Problem (<ref>) is known as the (dynamic) Schrödinger problem.
Its value, the entropic transportation cost, and its optimal solution, the Schrödinger bridge, are respective probabilistic counterparts to the transportation cost and displacement interpolation.
This problem plays an important role in optimal transport theory.
Indeed, the (dynamic) Schrödinger's problem is strongly connected to the (dynamic) Monge–Kantorovich optimal transport problem,

see  [Mikami, 2004] or  [Léonard, 2012, Léonard, 2016].
We refer to  [Léonard, 2014] for an extensive (but by now somewhat outdated) survey on the topic, and to the more recent survey by  [Chen et al., 2021] and the lecture notes of  [Nutz, 2022] for expositions of the theory.
Let us nonetheless mention that the static version of Schrödinger problem, obtained by projecting[We refer to [Léonard, 2014] for a detailed statement on the connection between the static and the dynamic versions of the Schrödinger problem.] onto the endpoint marginals has received sustained attention in recent years mostly due to the striking benefit of entropy penalization in computational optimal transport,

see  [Cuturi, 2013]
and  [Altschuler et al., 2017].

In this context, but from a theoretical point of view, we also mention  [Bernton et al., 2022, Bernton et al., 2022] and  [Nutz and Wiesel, 2022, Nutz and Wiesel, 2022] who study the stability and convergence of static Schrödinger bridges and Schrödinger potentials, respectively.
This work will focus on the dynamic Schrödinger problem (<ref>).
It is less studied from the probabilistic standpoint, but has found several practical applications where the so-called interpolating flow is essential.
These include problems in economics, see  [Galichon et al., 2015], network routing, see  [Chen et al., 2017],
image interpolation, see  [Papadakis et al., 2014] and  <cit.>, and statistics, see  [Hamdouche et al., 2023], among others.[For completeness we also mention the works  [Monteiller et al., 2019], and  [Lavenant et al., 2018] and  [Solomon et al., 2015].]

We will develop a new stochastic control and purely probabilistic approach to the dynamic Schrödinger problem complementing works by   [Mikami and Thieullen, 2006, Mikami and Thieullen, 2008] and  [Chen et al., 2016, Chen et al., 2017].

Let us present the problem in which we will be interested: Let $T>0$ and a positive integer $d$ be fixed, and denote by $(\Omega, \mathcal{F},\PP)$ the infinite product space of $\Cc([0,T],\R^d)$, endowed with the topology of uniform convergence induced by the sup norm $\|\cdot\|_\infty$, carrying a sequence of independent $\mathbb{R}^d$-valued Brownian motions $(W^i)_{i\geq 1}$.
For every positive integer $N$, let $\mathcal{ F}_0$ be an initial $\sigma$-field independent of $W^1,\dots, W^N$.
We equip $\Omega$ with the filtration $\mathbb{F}^N:=(\mathcal{F}_t^N)_{t\in[0,T]}$, which is the $\PP$-completion of the filtration generated by $W^1, \dots, W^N$ and $\mathcal{F}_0$, and we further denote by $\mathbb{F}^i:=(\mathcal{F}_t^i)_{t\in[0,T]}$ the $\PP$-completion of the filtration generated by $W^i$ and $\mathcal{F}_0$.
Without further mention, we use the identifications

\begin{equation*}
W \equiv W^1;  \quad \cF \equiv \cF^1 \quad \text{and} \quad \FF \equiv \FF^1.
\end{equation*}

It is well–known from the work of <cit.> (see also  [Léonard, 2012]) that (<ref>) admits the reformulation

\begin{equation}
\label{eq.dyn.SP.control}
\Vc_{\rm S}(\muin,\mufin)=\inf\bigg\{  \E^\P\bigg[\int_0^T \frac1{2}  \|\alpha_t\|^2\d t\bigg]:  \d X_t=\alpha_t\d t+\d W_t, \, \P\circ (X_0)^{-1}= \muin,\, \P\circ (X_T)^{-1}= \mufin\bigg\},
\end{equation}

where the infimum is taken over control processes $ \alpha$ which are $\R^d$–valued, $\F$–adapted and square integrable.

This reformulation gives $\Vc_{\rm S}(\muin,\mufin)$ as the value of a (constrained)

stochastic control problem,

an idea that has been very successfully exploited by [Mikami and Thieullen, 2006, Mikami and Thieullen, 2008]

to study the asymptotic behaviour of the zero–noise limit

as well as by  [Chiarini et al., 2022] and  [Tan and Touzi, 2013] to study the case of particles following Langevin and semimartingale dynamics instead of Brownian motion, respectively.

A crucial feature of this rich theory

is the independence of the particles $(X_i)_{i\geq 0}$.
In this work, we study the most likely path of a cloud of interacting particles with given initial and final configurations.
That is, given an interaction potential $\Psi$, we consider

the problem

\begin{align}
\notag
\Vc^N_{\rm e}(\muin,\mufin):=\inf\bigg\{\E^{\P}\bigg[ \frac1{2N} \sum_{i=1}^N \int_0^T\! \|\alpha_t^i\|^2\d t\bigg]: &\;  \d X_t^i=\alpha_t^i+\sum_{j=1}^N \nabla\Psi (X_t^i- X_t^j )\d t+\d W_t^i,\\\label{eq:N.SP.intro}
&  \; \P\circ (X_0^i)^{-1} =\muin,\; \P\circ ( X_T^i)^{-1} = \mufin  ,\; 1\leq i\leq N\bigg\},
\end{align}

where the infimum is taken over vectors $(\alpha^1,\dots,\alpha^N)$ such that each $\alpha^i$ is $\R^d$–valued, $\F^N$–adapted and square integrable.

This problem was first studied by  [Backhoff-Veraguas et al., 2020].
Because of the high intractability of the problem for $N$ large,

these authors borrowed inspiration from the theory of McKean–Vlasov control to introduce the so-called mean field (entropic) Schrödinger problem

\begin{align*}
\Vc_{\rm e}(\muin,\mufin):=\inf\Big\{\mathcal{H}(\Q |\Gamma(\Q)) :&\;   \d X_t= -(\nabla \Psi \ast \Q\circ (X_t)^{-1} ) (X_t)  \d t+\d W_t, \\
&\;\Gamma(\Q):=\P\circ( X)^{-1},\, \Q \circ (X_0)^{-1}=\muin, \, \Q \circ (X_T)^{-1}=\mufin  \Big\},
\end{align*}

where $\ast$ denotes the convolution operator, and the infimum is taken over probability measures on $\Omega$. The value $\Vc_{\rm e}(\muin,\mufin)$ and its optimizers are referred to as the mean field entropic transportation cost and mean field Schrödinger bridges (MFSB for short), respectively.
In complete analogy to the independent case, [Backhoff-Veraguas et al., 2020] showed that

\begin{align}\label{MFSP:entropic}
\Vc_{\rm e}(\muin,\mufin)=\inf\bigg\{\E^\P\bigg[\int_0^T  \frac1{2}  \|\alpha_t\|^2 \diff t\bigg]:&  \; \d X_t =\big( \alpha_t - (\nabla \Psi \ast \P\circ (X_t)^{-1}) (X_t)\big)\d t +\d W_t,\notag \\
&  \; \P\circ (X_0)^{-1}=\muin,\, \P\circ  (X_T)^{-1}=\mufin  \bigg\},
\end{align}

where the infimum is taken over $\R^d$–valued, $\F$–adapted and square integrable processes $\alpha$.
In this setting, the optimal states $\hat X$ correspond, via their law, to mean field Schrödinger bridges.

The main objective of the present work is to give conditions guaranteeing that the mean field Schrödinger problem actually arises as the limiting value of the sequence of finite particle Schrödinger problems (<ref>).
That is, rigorously proving the limit	$\Vc_{\rm e}^N\to \Vc_{\rm e}$ as $N$ goes to infinity.
This is strongly connected to the so-called convergence problem in mean field game theory  [Cardaliaguet et al., 2019].
Our approach will be based on stochastic calculus arguments reminiscent of those of [Laurièrie and Tangpi, 2022] [Laurièrie and Tangpi, 2022] that allow reducing the convergence problem into a “forward-backward propagation of chaos” issue.

Due to the singularity of the terminal cost, our analysis begins with a suitable penalization argument that allows us to to remove the constraint on the law of $X_T$, thus transforming the stochastic optimal transport problem into a “standard” stochastic control problem of McKean-Vlasov type.
A map $g:\Pc_2(\R^m)\longrightarrow [0,\infty)$ is said to be a smooth penalty function if: $(i)$ $g(\mu)=0$ if and only if $\mu=\mufin$; $(ii)$ $g$ is L-differentiable.
We refer to <Ref> for details as well as a reminder of the notion of $L$-differentiability and examples of penalty functions.
Given a smooth penalty function $g$, for each $k\ge 1$, we consider the penalized problem

\begin{equation}
\label{eq:inf.entropic.problem-k.penalized}
V^{k}_{\rm e}(\muin, \mufin) := \inf \bigg\{ \EE\bigg[\int_0^T\frac12 \|\alpha_t\|^2 \d t + kg(\mathcal{L}(X^\alpha_T)) \bigg]:    \; \d X_t =\big( \alpha_t - (\nabla \Psi \ast \P\circ (X_t)^{-1}) (X_t)\big)\d t +\d W_t, \; \P\circ (X_0)^{-1}=\muin \bigg\},
\end{equation}

which, by strong duality, allows us to see $\Vc_{\rm e}$ as the limit, as $k$ goes to infinity, of the penalized problems (<ref>), i.e., $V_{\rm e}^{k}\to \Vc_{\rm e}$ as $k$ goes to infity.
By introducing the corresponding penalized version $V^{N,k}_{\rm e}$,

the stochastic maximum principle recently developed by  [Carmona and Delarue, 2015]

reduces the convergence problem to a forward-backward propagation of chaos with the caveats that for any $k$, the $N$-particle system associated to $V^{N,k}_{\rm e}$ is of McKean–Vlasov type itself!
This makes the results in [Laurièrie and Tangpi, 2022]

Moreover, using a synchronous coupling argument along with techniques developed in [Laurièrie and Tangpi, 2022] would require imposing conditions on the coefficients of the particle system that depend on $k$, which is not desirable given that $k$ needs to go to infinity eventually.
One way we propose to get around these issues is to impose a convexity property on the penalty function well-known in optimal transport theory, see  [McCann, 1997], called displacement convexity, see also  <cit.> for its equivalent formulation in terms of $L$-differentiability which we adopt here.
The relevance of this type of convexity in the convergence problem in mean field games was recently highlighted by  [Jackson and Tangpi, 2023].

In addition to investigating the large particle limit, we also derive

new existence and characterization results for mean field Schrödinger bridges.
We show that for interaction functions that are nonlinear in the law of the state value, the mean field Schrödinger problem admits an optimal Schrödinger bridge that, in addition, can be characterize as the solution of a planning FBSDE.
This extends the results first obtained in [Backhoff-Veraguas et al., 2020] using very different arguments.

While the focus of this work is on the case of interacting particles, the approach (and some results) are new even for the classical problem (<ref>).
For instance, we will consider non–quadratic cost functions as well as general drifts.
Let us illustrate our main findings on the entropic Schrödinger problem discussed so far.

*Main results for the entropic Schrödinger problem.

Our first set of results provides the existence of an MFSB and establishes a connection between the above problem and a class of McKean-Vlasov FBSDEs.

Suppose $\Vc_{\rm e}(\muin,\mufin)<\infty$ and let $\Psi: \RR^\xdim  \longrightarrow \RR$ be twice differentiable with Lipschitz–continuous derivatives.

* the mean field entropic Schrödinger problem (<ref>) admits a mean field Schrödinger bridge $\hat X$ with associated optimal control $\hat \alpha$;

* the MFSB $\hat X$ and its associated optimal control $\hat \alpha$ satisfy

$$\hat X \stackrel{\text{d}}{=} \mybar X,\; \hat\alpha_t \stackrel{\text{d}}{=} - \mybar Y_t, \text{ for almost every } t\in [0,T],$$

for $(\mybar X, \mybar Y, \mybar Z)$, defined on a probability space $(\mybar\Omega,\mybar \cF, \mybar \P)$ supporting a Brownian motion $\mybar W$, satisfying the forward backward stochastic differential equation

\begin{equation*}
\begin{cases}
\diff \mybar X_t = \mybar \alpha_t -\tilde \E [  \nabla \Psi (\mybar X_t-\tilde X_t)]  \diff t + \diff
\mybar W_t\\
\diff \mybar Y_t = \big( \tilde \E  [ \nabla^2 \Psi( \mybar X_t-\tilde X_t)\cdot \mybar Y]-\tilde \E [\nabla^2 \Psi (\tilde X_t-\mybar X_t)\cdot \tilde Y_t ]\big)  \diff t + \mybar Z_t\diff
\mybar W_t\\
\mybar X_0\sim \muin,\quad \mybar X_T \sim \mufin,
\end{cases}
\end{equation*}

where $\tilde X$ and $\tilde Y$ denote independent copies of $ \mybar X$ and $\mybar Y$, respectively, on some probability space $(\tilde \Omega, \tilde \Fc, \tilde \P)$ and $\tilde \E$ denotes the associated expectation operator. In particular, the process $M$ given by

\[
M_t:= \mybar Y_t -\int_0^t\big( \tilde \E [  \nabla^2 \Psi ( \mybar X_r-\tilde X_r) \cdot Y_r ]-\tilde \E [\nabla^2 \Psi (\tilde X_r-\mybar X_r)  \cdot \tilde Y_r] \big) \diff r
\]

is a $\mybar \P$-martingale.

* if in addition $\Psi$ is symmetric,[ This is, $\Psi(x)=\Psi(-x), x\in \R^\xdim$.] then $\hat\alpha_t \stackrel{\text{d}}{=}  \mybar Y_t $ for almost every $t\in [0,T]$, and the above system simplifies to

\begin{equation}\label{eq:fbsde.MFSP.charact}
\begin{cases}
\diff \mybar X_t = \mybar Y_t -\tilde \E [  \nabla \Psi (\mybar X_t-\tilde X_t)]  \diff t + \diff \mybar W_t\\
\diff \mybar Y_t = \tilde \E [ \nabla^2 \Psi (\mybar X_t-\tilde X_t)\cdot (\mybar Y_t-\tilde Y_t) ] \diff t + \mybar Z_t\diff \mybar W_t\\
\mybar X_0\sim \muin,\quad \mybar X_T \sim \mufin.
\end{cases}
\end{equation}

<Ref> ensures the existence of an MFSB for (<ref>).
In addition, its associated optimal control is given in terms of a solution to the FBSDE system (<ref>).
Recalling that in the case of non-interacting particles the optimal control is a martingale, see  <cit.>, we complement our results by identifying the form of the martingale process associated with the optimal control in the case the particles interact through the potential $\nabla \Psi$.
We remark that under the additional assumption of $\muin$ and $\mufin$ being absolutely continuous with respect to the Lebesgue measure, (<ref>) and the corresponding martingale associated with the optimal control had already appeared in <cit.>, see the Assumption ${\rm (H2)}$ and Equation $(4)$ therein. Here $\muin$ and $\mufin$ need not be absolutely continuous.
Observe that, similar to <cit.>, the second statement above provides a characterization of an optimizer.

However, our results go well beyond the entropic problem.
In fact, our analysis in <Ref> generalizes <Ref> to a class of Schrödinger problems beyond (<ref>), see <Ref> for details.
In particular, we show that the MFSB and the associated optimal control can be derived as limits (in law) of a sequence of solutions of McKean–Vlasov equations derived from (<ref>).
Said result holds for a fairly general class of drifts and for cost functionals that are not necessarily quadratic in the control process and are allowed to depend on the state and law of the controlled process.
Even so, <Ref> is interesting not only because of its generality but also because it ellucidates a legitimate method to compute MFSBs by means of McKean–Vlasov equations, objects for which there are known theoretical guarantees for numerical implementations, see for instance the recent results in  [Baudelet et al., 2023],  [Han et al., 2022], as well as the review of  [Lauriére, 2021].
What's more, the same reasoning reveals a procedure to numerically approximate solutions to (<ref>).
This should be contrasted with the naïve approach of trying to solve (<ref>) directly for which, to the best of our knowledge, there is no theoretical analysis or numerical implementations at this level of generality.

Our second and main result is about the convergence of the finite particle problem (<ref>) to its mean field counterpart (<ref>).

Suppose $\Vc_{\rm e}(\muin,\mufin)<\infty$ and let $\Psi: \RR^\xdim  \longrightarrow \RR$ be twice continuously differentiable with Lipschitz continuous derivatives and

$\mufin$ admits a smooth penalty function which is displacement convex.[See <Ref> for the definition of displacement convexity.]

Then, we have that

\begin{equation*}
\lim_{N\to \infty}\Vc^N_{\rm e}(\muin,\mufin) = \Vc_{\rm e} (\muin,\mufin).
\end{equation*}

As alluded to above, our arguments reduce the convergence problem to a type of forward-backward propagation of chaos result which does not follow from standard results since the terminal condition of the $N$-particle system depends on the law of the forward process $X$ due to constraint on the law of $X_T$.
<Ref> presents one scenario under which the finite particle approximation holds.
The result holds for arbitrary finite time horizons whenever the target measure $\mufin$ admits a displacement convex smooth penalty function.
As with our characterization, the analysis in <Ref> extends <Ref> to a class of Schrödinger problems beyond (<ref>), see <Ref> and <Ref> for details.

The reader familiar with mean field problems might recognize the role of the previous convexity condition but wonder whether there are other known scenarii in which the convergence is possible.
Indeed, in the study of the convergence problem in mean field games (and uniqueness of mean field games), either small time conditions or monotonicity/convexity conditions are always needed to achieve strong convergence results [Cardaliaguet et al., 2019, Jackson and Tangpi, 2023, Laurièrie and Tangpi, 2022].
In our setting, the FBSDE associated with $V^{N,k}_{\rm e}$ is also well-behaved over a small time horizon, said $T_o^k$.
Yet, the value of $T_o^k$ depends on the Lagrange multiplier $k$ (appearing in the terminal condition), which makes this approach unviable since, as explained above, $k$ needs to go to infinity eventually.
We elaborate more on this in <Ref>, and mention that this is precisely the motivation for <Ref>, our second convergence result, which we refer to as the case of regular data.

We direct the reader to see <Ref> for the presentation of this setting, and to Remarks <ref> and <ref> for its implications to the entropic problem.

*Organization of the paper.
The rest of this paper is organized as follows: we dedicate the rest of this section to fixing some frequently used notation.
<Ref> begins presenting the formulation of the mean field Schrödinger problem and the standing assumptions on the data of the problem. It also introduces the penalization argument. <Ref> provides the proofs of the existence result as well as the characterization in the general setting, see <Ref>, which we then specialize to the entropic case of <Ref>.
<Ref> introduces the finite particle Schrödinger problem and presents the strategy for the analysis. <Ref> and <Ref> present the proof of the finite particle approximation under convexity assumptions and in the case of regular data, <Ref> and <Ref>, respectively. <Ref> brings the analysis back to the entropic problem presenting the proof of <Ref>.
An Appendix section includes some technical auxiliary results.

*Notation. Let $p$ be a positive integer.

For $(\alpha,\beta) \in \R^p\times\R^p$, $\alpha\cdot \beta$ denotes their inner product, with associated norm $\|\cdot\|$, which we simplify to $|\cdot|$ when $p$ is equal to $1$.

$I_p$ denotes the identity matrix of $\R^{p\times p}$.

Let $E$ be an arbitrary Polish space with a metric $d_E$.
Given a positive integer $\ell$ and a vector $x \in E^\ell$, we will denote by

the empirical measure associated with $x$.
We denote by $\Pc(E)$ the set of probability measures on $E$ (endowed with $\Bc(E)$, the Borel $\sigma$-algebra for the topology induced by $d_E$) and by $\Pc_p(E)$ the subset of $\Pc(E)$ containing measures with finite $p$-th moment.
Notice then that for any $x\in E^\ell$, we have $L^\ell(x)\in\Pc_p(E)$.
On a product space $E_1\times E_2$, $\pi^i$ denotes the projection on the $i$-th coordinate and for $\nu\in \Pc(E_1\times E_2)$, $\pi^i_{\#}\nu$ denotes the pushforward measure, i.e. $\pi^i_{\#}\nu (A):=\nu\circ (\pi^i)^{-1}(A)\in \Pc(E_i)$, for $A\in E_i$. For $\mu, \nu\in \Pc_p(E)$, we denote by $\cW_{p}(\mu, \nu)$ the $p$-Wasserstein distance between $\mu$ and $\nu$, that is, letting $\Gamma(\mu,\nu)$ be the set of all couplings $\pi$ of $\mu$ and $\nu$, i.e., probability measures on $(E^2,\Bc(E)^{\otimes 2})$ with first and second marginals $\mu$ and $\nu$, respectively, we have that

\begin{equation*}
\cW_{p}(\mu, \nu) := \bigg(\inf_{\pi\in\Gamma(\mu,\nu)}\int_{E\times E}d_E(x,y)^p\pi(\mathrm{d}x,\mathrm{d}y) \bigg)^{1/p}.
\end{equation*}

For a random variable $X$ on $(\Omega,\Fc,\P)$, we write $X\sim \mu$ whenever $\Lc(X)=\mu$, where $\cL(X)$ denotes the law of a random variable $X$ with respect to $\PP$. Given random variables $X, Y$, we write $X\stackrel{d}{=}Y$ whenever they are equal in distribution.

§ PROBLEM STATEMENT AND PENALIZATION

As already mentioned in the introduction, thanks to the flexibility afforded by our method, we will actually prove more general results than announced.
In particular, we will allow the interaction potential and cost to be generic, possibly nonlinear functions of the measure argument.
Let us thus present the general class of problems under consideration.

Fix a nonnegative integer $m$ and a function $f_1:[0,T]\times A \to \R$, we let $\mathfrak{A}$ be the set defined as

\begin{equation*}
\mathfrak{A} := \bigg\{\alpha:[0,T]\times \Omega \to A,\,\, \FF\text{--progressive processes such that } \EE\bigg[\int_0^Tf_1(t,\alpha_t) \diff t\bigg]<\infty\bigg\},
\end{equation*}

with $A$ a given nonempty, convex, and closed subset of $\mathbb{R}^m$.[In the typical example, $A=\R^\xdim$.]

Given two probability distributions $\muin, \mufin \in \mathcal{P}_2(\mathbb{R}^m)$, we are interested in the general stochastic optimal transport problem

\begin{equation}
\label{eq:problem}
V(\muin, \mufin) := \inf\bigg\{\EE\bigg[\int_0^Tf(t, X^\alpha_t, \alpha_t,\mathcal{L}(X^\alpha_t))\diff t \bigg]: \alpha \in \mathfrak{A},\,  X_0^\alpha\sim\muin,\, X_T^\alpha \sim   \mufin \bigg\},
\end{equation}

where $X^\alpha$ is the unique strong solution of the controlled stochastic differential equation

\begin{equation}
\label{eq:control.SDE.MF}
X^\alpha_t = X_0^\alpha + \int_0^t\alpha_u+b(u, X^\alpha_u, \mathcal{L}(X^\alpha_u))\diff u + \sigma W_t,\; t\in [0,T],
\end{equation}

where $\sigma\in  \R^{\xdim\times d}$. The assumptions on the functions $b$, $f_1$ and $ f$ as well as the matrix $\sigma$ will be introduced below.

Given $\muin, \mufin \in \cP_2(\RR^\xdim)$ we call $V(\muin,\mufin)$, with some abuse of language, the value of the mean field Schrödinger problem and when the cost $f$ is given by $f\equiv \|a\|^2/2$, then $V(\muin,\mufin)$ is called the mean field entropic transportation cost.
Moreover, the optimal states of $V(\muin,\mufin)$ are mean field Schrödinger bridges.

Problem (<ref>) means that we are looking for the “best trajectory”, in the sense of minimal cost $f$, to transport the probability distribution $\muin$ to $\mufin$. We will assume $V(\muin, \mufin)<\infty$, i.e. (<ref>) is feasible, and refer to [Backhoff-Veraguas et al., 2020] for conditions on $\muin$ and $\mufin$ guaranteeing it.

We make a few comments on the nature of Problem (<ref>)$:$

* In light of the presence of the terminal condition, we may interpret $V(\muin,\mufin)$ as a planning McKean–Vlasov stochastic control problem, i.e., the cooperative counterpart of the mean field game planning problem introduced in  [Lions, 2007].

* Choosing $\mu_{\rm in}=\delta_{x_0}$, $\mu_{\rm fin}=\delta_{0}$ we may reinterpret $V(\muin,\mufin)$ as a $($mean field$)$ optimal liquidation problem where the goal is to optimally steer the total inventory $X$ of a trader starting with $x_0$ to the value $X_T=0$, see for instance

[Fu et al., 2021],  [Graewe et al., 2015]

and the references therein.

* When $f(t, x, a,\mu)= \|a\|^2/2$, and, $b(t, x, \mu) = \int_{\R^m}\nabla\Psi(x-\tilde x)\mu(\d \tilde x)$, then Problem (<ref>) becomes the mean field entropic Schrödinger problem (<ref>).
If in addition $b=0$, it becomes Problem (<ref>).
These two settings will serve as our standard examples.

We now present the standing assumptions used in this article.
To that end, let us recall the notion of $L$-derivative on the space of probability measures, see for instance  [Lions, 2007] or  <cit.> for further details.
The definition is based on the notion of the lifting of a map $\phi:\mathcal{P}_2(\mathbb{R}^m)\longrightarrow \mathbb{R}$ to functions $\Phi$ defined on the Hilbert space $\L^2(\tilde \Omega ,\tilde \Fc ,\tilde \P;\mathbb{R}^m)$, over some atomless probability space $(\tilde \Omega,\tilde \Fc,\tilde \P)$, by setting $\Phi:\L^2(\tilde \Omega,\tilde \Fc,\tilde \P;\mathbb{R}^m)\longrightarrow \mathbb{R}$, $\Phi(\xi)=\phi(\Lc(\xi))$.

A function $\phi:\mathcal{P}_2(\mathbb{R}^m)\longrightarrow \mathbb{R}$ is said to be continuously $L$-differentiable if the lifted function $\Phi$ is continuously Fréchet differentiable with derivative $\partial_\xi \Phi:\L^2( \Omega,  \Fc, \P;\mathbb{R}^m)\longrightarrow \L^2( \Omega,  \Fc, \P;\mathbb{R}^m)$.
In this case, there is a measurable function $\partial_\mu \phi(\mu): \mathbb{R}^m\longrightarrow \mathbb{R}^m$, which we call the L-derivative of $\phi$, which for $\xi, \zeta\in \L^2(  \Omega,  \Fc,  \P;\mathbb{R}^m)$, $\Lc(\xi)=\mu$, satisfies, $\partial_\mu \phi(\mu)(\xi)=\partial_\xi \Phi(\xi)$, and

\[
\phi(\Lc(\xi+\zeta))-\phi(\mu)= \E[ \partial_\mu \phi(\mu)(\xi)\cdot \zeta]+o(\|\zeta\|_{\L^2}).
\]

In particular, the derivative is uniquely defined up to an additive constant.
To ease the notation, the above definition of $L$-differentiability assumes the lifting is defined on the same probability space in which the mean field Schrödinger problem was defined, but any probability space supporting a random variable with uniform distribution on $[0, 1]$ can be used.

We assume the matrix $\sigma$ is nonzero as well as:

* The function $b:[0,T]\times \RR^\xdim \times  \cP_1(\RR^\xdim)\longrightarrow \RR^m$ is Borel–measurable, continuously differentiable and Lipschitz–continuous in its last two arguments uniformly in $t$.
That is,

\begin{align*}
\|b(t, x,\mu) - b(t, x', \mu')\| \le \ell_b\big( \|x-x'\| + \cW_1(\mu,\mu')\big),\; \|b(t, x, \mu)\| &\le \ell_b\Big(1+ \|x\| + \Big(\int_{\RR^\xdim}\|x\|\mu(\d x)\Big)\Big),
\end{align*}

for all $t \in [0,T]$, $(x,\mu), (x', \mu') \in \RR^\xdim \times \cP_1(\RR^\xdim)$ and some constant $\ell_b>0$.

*  The function $f:[0,T]\times \RR^\xdim \times \R^d \times \cP_1(\RR^\xdim)\longrightarrow \RR$ is Borel–measurable and of quadratic growth. That is,

\begin{equation*}
|f(t, x, a,\mu)|	\le \ell_f\Big(1 +  \|x\|^2 + \|a\|^2 + \int_{\RR^\xdim}\|x\|^2\mu(\d x)\Big),
\end{equation*}

for all $(t,x,a,\mu) \in [0,T] \times \RR^\xdim \times A\times \cP_2(\RR^\xdim)$ and some constant $\ell_f>0$.
The function $f$ admits the decomposition

\begin{equation*}
f(t, x, a, \mu) = f_1(t, a) + f_2(t, x, \mu),
\end{equation*}

where, $f_2$ is bounded from below and continuously differentiable in the last two arguments, and, $f_1$ is continuously differentiable in the second argument, convex, positive, and satisfies $f_1(t,0)=0$. Moreover, there is $p>1$ and $C_1\in\R,C_2>0$, such that

$$f_1(t,a)\geq C_1+C_2\|a\|^p.$$

There is a Borel–measurable function $\Lambda:[0,T] \times \mathbb{R}^\xdim  \longrightarrow A$, that is Lipschitz–continuous and with linear growth in the last arguments with Lipschitz–constant $\ell_\Lambda$, such that

\begin{equation}
\label{eq:def.H}
\Lambda(t, y) \in \partial_y H_1(t, y), \text{ with}\;  H_1 (t, y): = \inf_{a\in A} \big\{ f_1 (t, a ) + a\cdot y\big\},
\end{equation}

where $\partial_y H_1$ denotes the set of subgradients of $H_1$, and such that the function $F:[0,T]\times \RR^\xdim\times \RR^\xdim\times \mathcal{P}_2(\RR^\xdim\times \RR^\xdim)\longrightarrow \RR^m$ given by

\begin{align}
\label{eq:def.F}
F(t,x,y,\nu) :=  \partial_x H_2 (t, x, y,\pi^1_{\#} \nu)
+\int_{\RR^\xdim\!\times \R^m}\partial_\mu H_2 (t, \tilde x  ,\tilde y,\pi^1_{\#} \nu )(x) \nu (\d \tilde x\times \d \tilde y) ,
% \\
%&+\int_{\RR^\xdim\times\RR^\xdim}\partial_\mu b\big(t, u, v, \Lambda(t, u, v, \xi),\varphi(t,\xi)\big)(x)y\xi(du, dv)
\end{align}

$$H_2(t,x,y,\mu):=f_2(t,x,\mu)+ b (t, x,\mu )\cdot y$$
is Lipschitz-continuous and with linear growth in its last three arguments with Lipschitz-constant $\ell_F$.

Let us comment on the previous set of assumptions. <Ref><ref> imposes Lipschitz and growth assumptions of the interaction term $b$ so that (<ref>) admits a unique strong solution $X^\alpha$. The quadratic growth condition of $f$ in <Ref><ref> is motivated by our baseline examples, see <Ref>. Our choice of the $\Wc_1$ norm for the measure term in $b$, as well as the convexity and coercivity conditions of $f_1$, would allow us to obtain the existence of optimizers and MFSBs without having to impose compactness on the set $A$ where the controls take values. Finally, as we will make use of the maximum principle for McKean-Vlasov control problems, <Ref><ref> introduces the Hamiltonian functional, the associated first-order sensitivity functional $F$, and imposes standard linear growth and Lipschitz conditions. Let us mention that (<ref>) is equivalent to $\partial_a f_1(t,\Lambda(t,y))=-y$,[Indeed, $\partial_y H_1(t,y)=\partial_a( f(t,a)+a \cdot y)^\star(0)$, where the superscript $\star$ denotes the Legendre transform, and $a\longmapsto f_1(t,a)$ is differentiable.] and if $f_1(t,\cdot)$ is strongly convex, the subgradient is a singleton (the gradient) and the Lipschitz–continuity assumption on $\Lambda(t,\cdot)$ follows from the inverse mapping theorem. These assumptions are easily verified in our benchmark example introduced in <Ref>, see <Ref>.

Let us now explain our strategy to analyse Problem (<ref>). In the language of stochastic control, the terminal cost of the stochastic optimal transport problem (<ref>) is the function $\chi_{\mufin}(\cL(X^\alpha_T))$, i.e., the convex indicator of the set $\{\mufin\}$ which is defined as $0$ when $\cL(X^\alpha_T) = \mufin$ and infinity otherwise, see <cit.>.
This terminal cost is clearly not differentiable, and consequently the problem is not amenable to standard techniques of stochastic control theory.
Thus, the starting point of our analysis is a penalization procedure allowing to write the optimal transport problem (<ref>) as a “standard” McKean–Vlasov control problem.

We say that $g:\Pc_2(\R^m)\longrightarrow [0,\infty)$ is a smooth penalty function if

* $g(\mu)=0$ if and only if $\mu=\mufin;$

* $g$ is L-differentiable.

Given a smooth penalty function $g$, for each $k\ge 1$, we consider the penalized problem

\begin{equation}
\label{eq:inf.problem-k.penalized}
V^{\infty,k}(\muin, \mufin) := \inf \bigg\{ \EE\bigg[\int_0^Tf(t, X^\alpha_t, \alpha_t,\mathcal{L}(X^\alpha_t))\diff t + kg(\mathcal{L}(X^\alpha_T)) \bigg]: \alpha \in \mathfrak{ A},\, X_0^{\alpha}\sim \muin\bigg\},
\end{equation}

where $X^\alpha$ is the unique strong solution of the controlled stochastic differential equation (<ref>).

Let us now discuss Problem (<ref>) and <Ref>; and provide some examples of penalty functions.
To begin with, note that (<ref>) is a standard stochastic control problem of McKean–Vlasov type (also called mean–field control problem).

Regarding the notion of smooth penalty function, our choice of <Ref><ref> is motivated by the above discussion on the terminal cost of the stochastic optimal transport problem (<ref>), but also by the following duality argument. Letting

\[
J(\alpha,k):=\EE\bigg[\displaystyle \int_0^Tf(t, X^\alpha_t, \alpha_t,\mathcal{L}(X^\alpha_t))\diff t + kg (\mathcal{L}(X^\alpha_T))\bigg],
\]

it is clear from <Ref><ref> that

\[
V(\muin, \mufin)=\displaystyle \inf_{\alpha\in \mathfrak{A}} \sup_{k \geq 1} J(\alpha,k), \text{ because }\,  \sup_{k \geq 1} J(\alpha,k ) =
\begin{cases}
\displaystyle \EE\bigg[\int_0^Tf(t, X^\alpha_t, \alpha_t,\mathcal{L}(X^\alpha_t))\diff t \bigg] & \text{if } \mathcal{L}(X^\alpha_T)=\mufin \\
\infty & \text{otherwise.}
\end{cases}
\]

Moreover, by weak duality

\begin{align}\label{eq:weakduality}
V(\muin, \mufin)\geq \displaystyle \sup_{k\geq 1} V^{\infty, k}(\muin,\mufin),
\end{align}

and it is natural that as $k$ goes to infinity, the penalized problems converge to the stochastic optimal transportation problem, i.e., equality in (<ref>) holds $($see <Ref> for details$)$. 		This will be crucial in our construction of solutions.
On the other hand, the smooth property of a penalty function, <Ref><ref>, will allow us to employ the stochastic maximum principle to analyze Problem (<ref>).

Regarding the choice of penalty function, at first glance, and by looking at the distance-like condition in <Ref><ref>, it seems appropriate to consider the mapping $g_1:\cP_1(\RR^\xdim)\ni \mu \longmapsto\mathcal{W}_1(\mu , \mufin)$.

this is not is a smooth function on the Wasserstein space, thus violating <Ref><ref>.
Alternatively, we could have considered the (perhaps more natural) penalty function $g_2:\cP_2(\RR^\xdim)\ni \mu \longmapsto\mathcal{W}_2^2(\mu , \mufin)$.
However, as noted by  [Alfonsi and Jourdain, 2020], this function is differentiable only when $\mufin$ is a Dirac mass, which would be quiet a restrictive condition but still amenable to applications in optimal liquidation, see <Ref><ref>.
One possible avenue to circumvent this would be to make use of the mollification theory for functions on the Wasserstein space $\cP_1(\RR^\xdim)$ recently developed by  <cit.>, see also  [Cosso and Martini, 2023].
Indeed, since the map $g_1$ is continuous and 1-Lipschitz on $\cP_1(\RR^\xdim)$, there exists a sequence $g^n_1:\cP_1(\mathbb{R}^\xdim) \to \RR$ of smooth and 1-Lipschitz functions that converges to $g_1$ locally uniformly.
This would have nevertheless led to an extra approximation level in our approach.
A candidate that bypasses this approximation is given by the Fourier-Wasserstein metric on $(\Pc_1(\R^\xdim),\Wc_1)$,
recently introduced in  [Soner and Yan, 2022] in the study of viscosity solutions to the dynamic programming equation.
At the intuitive level, the Fourier-Wasserstein metric $\hat \rho$ corresponds to a smooth equivalent of $\Wc_1$ and leads to the penalty function $g_3:\cP_1(\mathbb{R}^\xdim) \ni\mu \to \hat \rho(\mu,\mufin)$.
Another typical example of a penalty function
are Gauge-type functions on $(\Pc_2(\R^\xdim),\Wc_2)$, introduced notably by  [Cosso et al., 2022], (see also  [Bayraktar et al., 2022]) in the study of viscosity solutions of the master equation.
Similar to the previous example, one may think of a Gauge function $\rho$ as a smooth equivalent of $\Wc_2$ suggesting the choice $g_4:\cP_2(\mathbb{R}^\xdim) \ni\mu \to \rho(\mu,\mufin)$.

§ MEAN FIELD SCHRÖDINGER PROBLEM: EXISTENCE AND CHARACTERIZATION

Given the diverse nature of the examples presented in the previous section and the growing literature on PDEs over the space of probability measures, in this section we work with a generic penalty function $g$ in the sense of <Ref>.

We now state the main result of this section which covers <Ref> as a particular case.

Let <Ref> be satisfied and let $g$ be a penalty function in the sense of <Ref>.
Then, the mean field Schrödinger problem (<ref>) admits a MFSB $\hat X$ with associated optimal control $\hat \alpha$. Moreover,

* $\hat X$ is the limit in law of $(\hat X^k)_{k\ge1}$,

where for each $k\ge1$, the process $(\hat X^k, \hat Y^k, \hat Z^k)$ solves the McKean–Vlasov equation

\begin{equation}
\label{eq:fbsde.reg.MFSP.Theorem}
\begin{cases}
\diff \hat X^{k}_t = B\big(t, \hat X^{k}_t, \hat Y^k_t, \cL(\hat X^{k}_t)\big)\diff t + \sigma \diff W_t\\
\diff \hat Y^{k}_t = -F\big(t,  \hat X^{k}_t, \hat Y^{k}_t, \cL(\hat X^{k}_t,\hat Y^{k}_t)\big) \diff t + \hat Z^{k}_t\diff W_t\\
\hat X^{k}_0\sim \muin,\; \hat Y^{k}_T = k\partial_\mu	g(\cL(X_T^{k}))(X^{k}_T),
\end{cases}
\end{equation}

whereby $F$ is the function given in <Ref> and $B$ is the function defined as

\begin{equation}
\label{eq:def.B}
B(t, x, y, \mu) := \Lambda(t,  y) + b(t,x, \mu ).
\end{equation}

* If in addition $a\longmapsto f_1(t,a)$ and $(x,\mu)\longmapsto f_2(t,x,\mu)$ are Lipschitz smooth uniformly in $t\in [0,T]$,[Given a function $f:[0,T]\times\R^\xdim \longmapsto \R$ which is continuously differentiable in the second argument, the mapping $x\longmapsto f(t,x)$ is Lipschitz smooth uniformly in $t\in [0,T]$, if there exists $\ell_{f}>0$ such that $\|\partial_x f(t,x)-\partial_x f(t,y)\ |\leq \ell_{f} \|x-y\|$, for all $x,y\in \R^\xdim$.] and $y\longmapsto \Lambda(t,y)$ is linear, then the MFSB $\hat X$ and its associated optimal control $\hat \alpha$ satisfy

\begin{align}\label{eq.charactirizationlimit}
\hat X \stackrel{\text{d}}{=} \mybar X, \text{ and, } {\color{black} \hat\alpha_t \stackrel{\text{d}}{=} \Lambda(t,\mybar Y_t)},\text{ for almost all } t\in [0,T],
\end{align}

for $( \mybar X, \mybar  Y,  \mybar Z)$, defined on a probability space $(\mybar \Omega, \mybar \Fc, \mybar \P)$ supporting a Brownian motion $\mybar W$, satisfying

\begin{equation}
\label{eq:fbsde.MFSP.charact.corollary}
\begin{cases}
\diff \mybar X_t =  \big(\Lambda(t,\mybar Y_t) -b(t,  \mybar X_t, \cL(\mybar X_t)) \big) \diff t + \sigma \diff \mybar W_t\\
\diff   \mybar Y_t = -F(t,\mybar X_t,\mybar Y_t,\Lc(\mybar X_t,\mybar Y_t))  \diff t +  \hat Z_t\diff \mybar W_t\\
\mybar X_0\sim \muin,\; \mybar X_T \sim \mufin.
\end{cases}
\end{equation}

In particular,

\[
\mybar M_t:= \mybar Y_t +\int_0^t F(r,\mybar X_r,\mybar Y_r,\Lc(\mybar X_r,\mybar Y_r))  \diff r,
\]

is a $\mybar \P$-martingale.

The above result complements the existence statement of a MFSB by: $(i)$ an asymptotic characterization in the case of general cost functionals $f$, $(ii)$ a characterization of the MFSB in terms of a so-called planning FBSDEs.
On the one hand, let us highlight that all the additional assumptions in the second part of the statement of <Ref> are satisfied in the entropic case, see <Ref>.
On the other hand, we remark that the asymptotic characterization is in itself interesting, for instance, for numerical simulation.
In fact, while estimating an MFSB seems to be a very difficult task, the above result suggests that one could obtain an MFSB by numerically solving the FBSDE (<ref>) for increasingly large values of $k$.
Further notice that the present characterization contrasts with <Ref> in that <Ref> does not seem amenable to numerical approximation because it does not have boundary conditions, whereas <Ref> is a standard McKean–Vlasov FBSDE for which numerical schemes exist.
This observation will be explored thoroughly in a separate work to present a numerical scheme for the computation of Schrödinger bridges.

I think we should say something about PDEs, at least in the entropic case.

[inline]Can we show that $\hat\alpha^k_t = v^k(t, \hat X^k_t, \Lc(\hat X^k_t))$? and thus, perhaps that $\hat\alpha_t = v(t, \hat X_t, \Lc(\hat X_t))$ for $t\in [0,T)$? We can under well-posedness of the FBSDE (or T small)

Let $d_2$ be a metric on $\Pc_2(\R^m)$ for which $(\Pc_2(\R^m),d_2)$ is complete. A map $\rho:\Pc_2(\R^m)\times \Pc_2(\R^m)\longrightarrow [0,+\infty)$ is said to be a gauge-type function if: $\rho(\mu,\mu)=0$ for all $\mu\in \Pc_2(\R^m)$; $\rho$ is continuous on $\Pc_2(\R^m)\times \Pc_2(\R^m)$; and for any $\eps>0$ there exists $\delta>0$ such that, for all $\mu,\nu\in \Pc_2(\R^m)$, $\rho(\mu,\nu)\leq \delta$ implies $d_2(\mu,\nu)\leq \eps$.

By <cit.>, for any $\varrho>0$

\[
\rho_{2,\varrho}(\mu,\nu):=c_m \sum_{\hat n\geq 0} 2^{2\hat n} \sum_{\hat \ell\geq 0 } 2^{-2\hat \ell} \sum_{B\in \Pc_{\hat\ell}} \bigg( \sqrt{ \big|(\mu\ast \Nc_\varrho)\big( (2^{\hat n} B)\cap B_{\hat n}\big) -(\mu\ast \Nc_\varrho)\big( (2^{\hat n} B)\cap B_{\hat n}\big) \big|^2 +\delta_{\hat n, \hat \ell}^2} -\delta_{\hat n, \hat \ell}\Big),
\]

is a gauge-type function on the complete metric space $\big(\Pc_2(\R^m),\Wc_2^{(\varrho)}(\mu,\nu)\big)$, see <cit.>, where

* $\Wc_2^{(\varrho)}(\mu,\nu):=\Wc_2(\mu\ast \Nc_\varrho,\nu\ast\Nc_\varrho)$ and $\Nc_\varrho:=\Nc(0,\varrho I_m)$ denotes an independent Normal random variable

Moreover, we note that the following holds

\[
\rho_{2,\varrho}(\mu,\nu)=c_m \sum_{\hat n\geq 0} 2^{2\hat n} \sum_{\hat \ell\geq 0 } 2^{-2\hat \ell} \sum_{B\in \Pc_{\hat\ell}} \bigg( \sqrt{ \big|\langle \psi^B_{\hat n} , \mu-\nu  \rangle\big|^2 +\delta_{\hat n, \hat \ell}^2} -\delta_{\hat n, \hat \ell}\Big)
\]

where .

We thus let $\varrho>0$ be fixed and consider the function

$$g(\mu):=\rho_{2,\varrho}(\mu, \mufin),\quad \mu \in \cP_2(\RR^\xdim).$$

§.§ Preparatory lemmas

As mentioned above, our analysis takes full advantage of the fact that Problem (<ref>) is a standard stochastic control problem of McKean–Vlasov type.
The literature on McKean–Vlasov control problems has quickly developed in recent years, notably due to its link to mean–field games.
We refer the reader for instance to  <cit.>;  [Carmona and Delarue, 2015];  [Djete et al., 2022],  [Cosso et al., 2022];  [Soner and Yan, 2022] and  [Pham and Wei, 2017].
The proof of <Ref> will need the following two intermediate lemmas:

Let <Ref> hold.
For every $k \ge 1$ the control problem (<ref>) with value $V^{\infty, k}(\muin,\mufin)$ admits an optimizer $\hat\alpha^{k}$ satisfying $\hat\alpha^{k}_ t= \Lambda (t, \hat Y^{k}_t )$, $\d t\otimes \d \P\text{\rm --a.e.}$, where $\Lambda$ is a function satisfying <Ref>, and $(\hat X^{k}, \hat Y^{k}, \hat Z^{k})$ solving the McKean–Vlasov equation

\begin{equation}
\label{eq:fbsde.reg.MFSP}
\begin{cases}
\diff\hat X^{k}_t = B(t, \hat X^{k}_t,\hat Y^{k}_t, \cL(\hat X^{k}_t))\diff t + \sigma \diff W_t\\
\diff \hat Y^{k}_t = -F(t, \hat X^{k}_t, \hat Y^{k}_t, \cL(X^{k}_t,Y^{k}_t)) \diff t + \hat Z^{k}_t\diff W_t\\
\hat X^{k}_0\sim \muin,\quad \hat Y^{k}_T = k\partial_\mu	g(\cL(\hat X_T^{k}))(\hat X^{k}_T).
\end{cases}
\end{equation}

In a first step, we extract a candidate optimizer. For every positive integer $n$ and $k\geq 1$, there is $\alpha^{k,n} \in \mathfrak{A}$ such that

\begin{equation}
\label{eq:lower.bound.Vinfkn}
V^{\infty,k}(\muin, \mufin) \ge \EE\bigg[\int_0^Tf(t, X^{\alpha^{k,n}}_t, \alpha^{k,n}_t,\cL(X^{\alpha^{k,n}}_t))\diff t + kg(\cL(X^{\alpha^{k,n}}_T)) \bigg] - \frac1n.
\end{equation}

Since $V^{\infty,k}(\muin,\mufin)\le V(\muin, \mufin)$ and the functions $f_2$ and $g$ are bounded from below, there is $C>0$ independent of $n$ such that

\begin{equation}
\label{eq:bound.f1}
\EE\bigg[\int_0^Tf_1(t, \alpha^{k,n}_t)\diff t\bigg] \le C.
\end{equation}

Thus, since $f_1$ is convex and coercive, and (<ref>) holds, by <Ref>,
the sequence $(\int_0^\cdot \alpha^{k,n}_t\diff t)_{n\geq 1}$ is tight and converges to $\int_0^\cdot \hat \alpha_t^k\diff t$ in law in $C([0,T],\RR^m)$ and by convexity of $f_1$ it holds

\begin{equation*}
\liminf_{n\to \infty}\EE\bigg[\int_0^Tf_1(t, \alpha^{k,n}_t)\diff s\bigg]\ge \EE\bigg[\int_0^Tf_1(t, \hat \alpha_t^k)\diff t\bigg].
\end{equation*}

Note that this last inequality shows that $\hat\alpha^k\in \mathfrak{A}$.With this, we now proceed to show the weak convergence of the associated controlled processes. We claim that $(X^n)_{n\geq 1}$, where $X^n:=X^{\alpha^{k,n}}$, is tight. For this note that there is $p>1$ and $C>0$ such that

\begin{align*}
\|X_t^n\|^p& \leq  C\bigg(\|X_0\|^p+ \int_0^T \| \alpha_r^{k,n}\|^p\diff r + \int_0^t \| b(r,X_r^n,\Lc(X_r^n))\|^p \diff r+ \|\sigma W_t\|^p\bigg)\\
& \leq   C\bigg(1+ \|X_0\|^p+ \int_0^T f_1(r, \alpha_r^{k,n}) \diff r + \int_0^t \big(  \|X_r^n\|^p +\E[\|X_r^n\|] ^p\big) \diff r+ \|\sigma W_t\|^p\bigg),
\end{align*}

where we use $f_1(t,a)\geq C_1+C_2\|a\|^p$. It then follows from (<ref>) and Gronwall's inequality that there is $C>0$ such that

\begin{align}\label{eq:boundXntight}
\sup_{n\geq 1} \E\bigg[\sup_{t\in [0,T]} \|X_t^n\|^p\bigg]\leq C.
\end{align}

Similarly, for the same $p>1$, we may find that there is $C>0$ such that

\begin{align*}
\|X_t^n-X_s^n\|^p \leq  C\bigg( |t-s|^{p-1}  \int_0^T \Big( f_1(r, \alpha_r^{k,n})  +\|X_r^n\|^p+\E[\|X_r^n\|^p]\Big) \diff r+ |t-s|^{\frac p2}\|\sigma\|^p\bigg),
\end{align*}

which, together with (<ref>), (<ref>) and Aldous criterion, see <cit.>, allows us to conclude that $(X^n)_{n\geq 1}$ is tight.We are now in a position to establish the weak convergence of $(X^{n})_{n\geq 1}$.

By Skorokhod's extension theorem, we can find a common probability space $(\mybar \Omega, \mybar \cF, \mybar \P)$ supporting $(\mybar X^n,\int_0^\cdot \mybar \alpha_t^{k,n} \diff t, \mybar W^n) \stackrel{d}{=} (X^n,\int_0^\cdot \alpha_t^{k,n} \diff t, W)$, for all $n\geq 1$, and $(\mybar X,\int_0^\cdot \mybar \alpha_t^{k} \diff t,\mybar W)\stackrel{d}{=} (X,\int_0^\cdot \hat \alpha_t^{k} \diff t, W)$, and on which $( \mybar X^{n}, \int_0^\cdot \mybar\alpha^{k,n}_t\diff t , \mybar W^n) \longrightarrow  (\mybar X, \int_0^\cdot \mybar{\alpha}_t^k\diff t,\mybar W)$, $\mybar \P\text{--a.s.}$
Let us note that by dominated convergence, $(\int_0^\cdot \mybar \alpha^{k,n}_t\diff t)_{n\geq 1}$ converges to $\int_0^\cdot {\mybar{\alpha}}_t^k\diff t$ in $\mathbb{L}^1(\mybar  \Omega,\mybar  \cF, \mybar \P)$.We now note that since $b$ is Lipschitz–continuous and $(\mybar X^n,\int_0^\cdot \mybar \alpha_t^{k,n} \diff t, \mybar W^n) \stackrel{d}{=} (X^n,\int_0^\cdot \alpha_t^{k,n} \diff t, W)$, $\mybar X^n$ solves, under $\mybar \P$, the controlled SDE

\begin{equation*}
\diff \mybar X^{n}_t = \mybar \alpha_t^{k,n} + b(t, \mybar X^{n}_t, \cL(\mybar X^{n}_t))\diff t + \sigma\diff \mybar W_t^n,\quad \mybar X^{n}_0 \sim \muin.
\end{equation*}

Letting $\mybar \E$ denote the expectation under $\mybar \P$, the Lipschitz continuity of $b$ and Gronwall's inequality leads to

\begin{equation*}\label{eq:estimateX.alpha}
\mybar \EE\bigg[\sup_{t\in [0,T]} \|\mybar X^{n}_t - \mybar X^{\mybar{\alpha}^k}_t \|\bigg] \le \mathrm{e}^{2 \ell_bT}\mybar \EE \bigg[\|\mybar X_0^n-\mybar X_0 \|+\sup_{t\in [0,T]} \bigg\|\int_0^t\mybar \alpha^{k,n}_s\diff s - \int_0^t\mybar {\alpha}_s^k\diff s\bigg\|+  \|\sigma\| \sup_{t\in [0,T]}\|\mybar W^n-\mybar W \|  \bigg].
\end{equation*}

Consequently, $\Wc_1\big(\Lc\big(\mybar X^{n}\big),\Lc\big(\mybar X^{\bar \alpha^k}\big)\big)=\Wc_1\big(\Lc\big( X^{n}\big),\Lc\big( X^{ \alpha^{k}}\big)\big)$ converges to zero as $n$ goes to infinity.

Thus, $(X^{\alpha^{k,n}})_{n\ge0}$ converges in law to $X^{\hat \alpha^k}$, where $X^{\hat \alpha^k}$ is uniquely given by the controlled SDE

\begin{equation*}
\diff X^{\hat \alpha^k}_t = \hat \alpha_t^k + b(t, X^{\hat \alpha^k}_t, \cL(X^{\hat \alpha^k}_t))\diff t + \sigma\diff W_t,\quad X^{\hat \alpha^k}_0 \sim \muin.
\end{equation*}

Therefore, taking the limit in $n$ in (<ref>), it follows by continuity of $f$ and $g$ that

\begin{equation*}
V^{\infty, k}(\muin, \mufin) \ge \EE\bigg[\int_0^Tf(t, X^{\hat \alpha^k}_t, \hat \alpha_t^k,\cL(X^{\hat \alpha^k}_t))\diff t + kg(\cL(X^{\hat \alpha^k}_T)) \bigg],
\end{equation*}

showing that $\hat \alpha^k$ is optimal since it is admissible.

At last, we note that the characterization follows from the maximum principle.
In fact, it follows by <cit.> that an optimal control $\hat\alpha^k$ satisfies

\begin{equation*}
\big(\partial_af_1(t, \hat\alpha_t^{k}) +  \hat Y^{k}_t\big )\cdot (\hat\alpha^{k}_t - a)\leq 0,\;  \text{for all } a \in A,\;  \d t\otimes \d \PP \text{--a.e.,}
\end{equation*}

where $\hat X^k:=X^{\hat \alpha^k}$ and $(\hat Y^{k}, \hat Z^k)$

satisfy (<ref>). Since $a\longmapsto f_1(t,a)+a\cdot y$ is convex, we see that $0$ is a subgradient, and consequently $\hat\alpha^{k}_t = \Lambda (t,\hat Y^{k}_t)$, $\d t\otimes \d \PP \text{--a.e.,}$ for a function $\Lambda$ as in (<ref>).

Under <Ref> we have

\begin{equation*}
\lim_{k \to \infty}V^{\infty,k}(\muin, \mufin)= V(\muin, \mufin).
\end{equation*}

Moreover, strong duality holds, i.e., there is equality in (<ref>).

Note that the sequence $(V^{\infty, k}(\muin, \mufin))_{k\ge0}$ is monotone increasing. Thus, the second statement follows from the first.
Let us show the first result. By weak duality, see (<ref>), since $V^{\infty,k}(\muin,\mufin)\le V(\muin, \mufin)<\infty$, we have that

\begin{equation*}
\lim_{k \to \infty}V^{\infty, k}(\muin, \mufin) \le V(\muin, \mufin).
\end{equation*}

The sequence $(V^{\infty, k}(\muin, \mufin))_{k\ge0}$ is monotone increasing and for all $k$ we have $V^{\infty, k}(\muin, \mufin) \le V(\muin, \mufin)$.
In fact,

\begin{align*}
V^{\infty, k}(\muin, \mufin) &\le \sup_{k \in \NN}\inf_{\alpha \in \mathfrak{A}}\EE\bigg[\int_0^Tf(t, X^{\alpha}_t, \alpha_t,\cL(X^{\alpha}_t))\diff t + kg(\cL(X^{\alpha}_T)) \bigg]\\
&\le \inf_{\alpha \in \mathfrak{A}}\sup_{k \in \NN}\EE\bigg[\int_0^Tf(t, X^{\alpha}_t, \alpha_t,\cL(X^{\alpha}_t))\diff t + kg(\cL(X^{\alpha}_T)) \bigg]\\
& = V(\muin, \mufin).
\end{align*}

where in the last step we used that, by definition of $g$, $g(\cL(X^{\alpha}_T))=0$  implies $\Wc_2^{(\varrho)}(\Lc(X^{\alpha}_T),\mufin)=0$, and the fact that $\varrho\longmapsto \Wc_2^{(\varrho)}(\mu,\nu)$ is non-increasing, see <cit.>. Thus, we have

Let us show the reverse inequality.
For every positive $k$, there is $\alpha^k\in \mathfrak{A}$ such that

\begin{equation}
\label{eq:upper.bound.lim.k}
V^{\infty,k}(\muin, \mufin) \ge \EE\bigg[\int_0^Tf(t, X^{\alpha^k}_t, \alpha^k_t,\cL(X^{\alpha^k}_t))\diff t + kg(\cL(X^{\alpha^k}_T)) \bigg] - \frac1k.
\end{equation}

Since $V^{\infty,k}(\muin,\mufin)\le V(\muin, \mufin)<\infty$ and the functions $g$ and $f_2$ are bounded from below, as in the proof of <Ref>, the sequence $(\int_0^\cdot \alpha^{k}_t\diff t)_{k\geq 1}$ is tight and converges to $\int_0^\cdot \hat \alpha_t\diff t$ in law, with $\hat \alpha\in \mathfrak A$, and consequently, $X^{\alpha^k}$ is tight and converges to $X^{\hat \alpha}$ in law. Moreover, there is $C>0$ independent of $k$ such that

\begin{align*}
g(\cL(X^{\alpha^k}_T)) &\le \frac1k\bigg(V^{\infty, k}(\muin,\mufin)- \EE\bigg[\int_0^Tf(t, X^{\alpha^k}_t, \alpha^k_t,\cL(X^{\alpha^k}_t))\diff t  \bigg] \bigg) + \frac{1}{k^2}\notag \\
&\le \frac1k\Big(V(\muin,\mufin)- C \Big) + \frac{1}{k^2}\xrightarrow [k\to \infty]{}0.
\end{align*}

Thus, by definition of $g$, $ X_T^{\hat \alpha}\sim\mufin$. Now, since $g$ is non-negative we have, thanks to the continuity of $f$, back in (<ref>) that

\begin{align*}
\lim_{k\to \infty}  V^{\infty,k}(\muin, \mufin) & \ge \lim_{k\to \infty} \EE\bigg[\int_0^Tf(t, X^{\alpha^k}_t, \alpha^k_t,\cL(X^{\alpha^k}_t))\diff t  \bigg] - \frac1k\\
&\geq  \EE\bigg[\int_0^Tf(t, X_t^{\hat \alpha}, \hat \alpha_t,\cL(X^{\hat \alpha}_t))\diff t  \bigg] .
\end{align*}

Since $\hat\alpha\in \mathfrak A$ and $X_T^{\hat \alpha}\sim \mufin$, we find that

\[
\lim_{k\to \infty} V^{\infty,k}(\muin,\mufin)\geq
\]

which yields the result.

This shows that there is a constant $C>0$, independent of $k$, such that $\alpha^k\in {
\mathfrak A^k}:=\{\alpha\in \mathfrak A: g(\Lc(X_T^\alpha))\leq C/k\}$.
Since $g$ is nonnegative, it then follows from (<ref>) that

\begin{align*}
\inf_{\alpha\in \mathfrak A^k}\EE\bigg[\int_0^Tf(t, X^{\alpha}_t, \alpha_t,\cL(X^{\alpha}_t))\diff t \bigg] - \frac1k.
\end{align*}

Now note that for any $k\geq k^\prime>0$, $\mathfrak A^k\subseteq \mathfrak A^{k^\prime}$, and $\mathfrak A^\infty:=\{\alpha \in \mathfrak A: \Lc(X^{\alpha}_T)=\mufin\} \subseteq \mathfrak A^k $ thanks to <Ref><ref>. Thus

\[
\lim_{k\to \infty} V^{\infty,k}(\muin,\mufin)\geq
\]

which yields the result.

§.§ The general case: Proof of Theorem <ref>

Let us argue the first statement and $(i)$ together.
By <Ref> for every $k\geq 0$ there exists $\hat\alpha^k\in \mathfrak{A}$ such that $\hat \alpha^k_t=\Lambda(t ,\hat Y^k_t )$, $\d t\otimes \d \PP \text{--a.e.,}$ with associated optimal state process $\hat X^{k}$ and adjoint process $\hat Y^{k}$ satisfying (<ref>) which is an optimizer of $V^{\infty,k}$, i.e.,

\begin{equation}
\label{eq:optimizer.Vk}
V^{\infty,k}(\muin, \mufin) = \EE\bigg[\int_0^Tf(t,\hat X^{k}_t, \hat\alpha^{k}_t, \cL(\hat X^{k}_t))\diff t + kg(\cL(\hat X^k_T)) \bigg].
\end{equation}

Once again, since $f_2, g$ are bounded from below and $V^{\infty,k}(\muin, \mufin)\le V(\muin,\mufin) <\infty$ we find that

\begin{equation}\label{eq:bound.f1.alphak}
\EE\bigg[\int_0^Tf_1(t, \hat\alpha^k_t)\diff t\bigg] \le C,
\end{equation}

for some constant $C>0$.
Thus, from <Ref> we obtain the convergence of $\hat A^k(\cdot):=\int_0^\cdot\hat \alpha^{k}_t\diff t $ to $\hat A(\cdot):=\int_0^\cdot \hat \alpha_t\diff t$ in law in $C([0,T],\R^\xdim)$, which in turn, arguing as in the proof of <Ref>, allows us to derive that $(\hat X^{k})_{k\ge1}$ converges in law to $\hat X^{\hat \alpha}$. Moreover, as in the proof of <Ref>, since $V^{\infty,k}(\muin,\mufin)\le V(\muin, \mufin)<\infty$ and the functions $g$ and $f$ are bounded from below, there is $C>0$ independent of $k$ such that

\begin{align*}
g(\cL(\hat X^{k}_T)) \le \frac1k\Big(V(\muin,\mufin)- C \Big) ,
\end{align*}

from which we deduce that $\Lc(\hat X_T)=\mufin$. All in all, $(\hat X^{k})_{k\ge1}$ converges in law to $\hat X^{\hat \alpha}$ satisfying

\begin{equation}\label{eq.existencethm.limitX}
\diff \hat X_t = \hat\alpha_t+ b(t,  \hat X_t, \cL(\hat X_t))\diff t + \sigma \diff W_t, \quad \hat X_0 \sim \muin, \quad \hat X_T \sim \mufin.
\end{equation}

This establishes $(i)$ provided that we verify the optimality of $\hat \alpha$. Back in (<ref>), we see that

\begin{align}
\label{eq:ine.lim.Vk}
\lim_{k\to \infty}V^{\infty,k}(\muin,\mufin)\ge \EE\bigg[\int_0^Tf(t, \hat X_t, {\hat \alpha}_t,\cL(\hat X_t))\diff t \bigg].
\end{align}

Since the reverse inequality in (<ref>) is clear by definition of $V^{\infty, k}(\muin, \mufin)$ and the fact that $\hat \alpha\in \mathfrak{A}$, it follows that

\begin{equation*}
\lim_{k\to \infty}V^{\infty,k}(\muin,\mufin) =
\EE\bigg[\int_0^Tf(t, \hat X_t, {\hat \alpha}_t,\cL(\hat X_t))\diff t \bigg]
\end{equation*}

and thus, by  <Ref> and uniqueness of the limit, we have

\begin{equation*}
V(\muin,\mufin) = \EE\bigg[\int_0^Tf(t, \hat X_t, {\hat \alpha}_t,\cL(\hat X_t))\diff t \bigg].
\end{equation*}

That is, $\hat\alpha$ is optimal. We now argue $(ii)$. We argue in three steps.Step 1. To begin with, recall that by $(i)$, as in <Ref>, the sequence $(\hat X^k,\int_0^\cdot \hat \alpha_t^{k} \diff t,  W)_{k\geq 1}$ is tight.
Thus, by Skorokhod extension theorem, there is a common probability space $(\bar \Omega, \bar \cF, \bar \P)$ supporting $(\mybar X^k,\bar A^k , \mybar W^k) \stackrel{d}{=} (X^k,  \hat A^k, W)$, for all $k\geq 1$, and $(\mybar X,\bar A  ,\mybar W)\stackrel{d}{=} (X,  \hat A, W)$, with $\bar A^k(\cdot):=\int_0^\cdot \mybar\alpha^{k}_t\diff t$, $\bar A(\cdot):=\int_0^\cdot \mybar{\alpha}_t\diff t$, and on which $(\mybar X^{k}, \bar A^k , \mybar W^k) \longrightarrow (\mybar X,  \bar A, \mybar W)$, $\mybar \P\text{--a.s.}$ as $k$ goes to infinity.
Moreover, $\mybar X$ satisfies (<ref>), $\mybar \P\text{--a.s.}$ We will write $\mybar \E$ for the expectation under $\mybar \P$ and let $\mybar \F^k:=\F^{\overline W^k}$ and $\mybar \F:=\F^{\overline W}$, e.g., $\mybar \F$ denotes the filtration generated by $\mybar W$.  To obtain the result, i.e., to arrive at (<ref>), we construct $(\mybar X^k, \mybar Y^k, \mybar Z^k)_{k\geq 1}$ on $(\mybar \Omega, \mybar \cF, \mybar \P)$, analog to (<ref>), and study its limit. Note that the previous paragraph defines $(\mybar X^k)_{k\geq 1}$ and identifies its limit. We are left to introduce the pair $(\mybar Y^k, \mybar Z^k)$.On the one hand, since <Ref> is equivalent to $y=-\partial_a f_1(t,\Lambda(t,y))$, we let

$$\mybar Y^k:=-\partial_a f_1(\cdot,\mybar \alpha_\cdot^k).$$

Notice that for any $k\geq 1$ and any bounded continuous $G:\R^\xdim \longrightarrow \R$ and almost every $t\in [0,T]$, $\E[G(h^{-1}(\hat A^k(t+h)-\hat A^k(t)))]=\E[G(h^{-1}(\bar A^k(t+h)-\bar A^k(t)))]\longrightarrow\E[G(\mybar \alpha_t^k)]$, as $h\longrightarrow 0$.
Thus, we have that $\mybar \alpha^k_t\stackrel{d}{=} \hat \alpha^k_t$ for almost every $t\in [0,T]$.
Thus, $\mybar Y_t^k\stackrel{d}{=} \hat Y^k_t$ for almost every $t\in [0,T]$, $k\geq 1$.
We now claim that, for $p>1$ given by <Ref><ref>, there exists $\mybar Y\in\L^p([0,T]\times \mybar \Omega)$ such that $\mybar Y^k$ converges to $\mybar Y$ weakly in $\L^p([0,T]\times \mybar \Omega)$.
Indeed, since $f_1$ satisfies  $f_1(t,a)\geq C_1+C_2\|a\|^p$, and $a\longmapsto f_1(t,a)$ is Lipschitz smooth, we have that

\[
\mybar \E \bigg[\int_0^T \|\mybar Y_t^k\|^p \diff t\bigg]=\mybar \E \bigg[\int_0^T \|\partial_a f_1(t,\mybar \alpha_t^k)\|^p \diff t\bigg]\leq \ell_{f_1}\mybar  \E \bigg[\int_0^T(1+\| \mybar \alpha_t^k \|^p )\diff t\bigg]\leq C,
\]

where the last inequality follows from (<ref>).
Consequently, $(\mybar Y^k)_{k\geq 1}$ (and $(\mybar \alpha^k)_{k\geq 1}$), being a bounded sequence in a reflexive normed space, contain weakly convergent subsequences, see  <cit.>.
That is, $(\mybar Y^k)_{k\geq 1}$ converges weakly in $\L^p([0,T]\times \mybar \Omega)$ to $\mybar Y\in\L^p([0,T]\times \mybar \Omega)$. In particular, $(\mybar Y^k_t)_{k\geq 1}$ converges to $\mybar Y_t$ for almost every $t\in [0,T]$. We now claim that $\mybar \alpha_t =\Lambda(t,\mybar Y_t)$, for almost every $t\in [0,T]$, thus establishing (<ref>). Indeed, note that the weak convergence of $\mybar Y^k_t$ implies, by Mazur's lemma that, up to a subsequence, there exists $(\gamma^{k})_{n\leq k \leq N(n)}$, $n\geq 1$, $\sum_{k=n}^{N(n)} \gamma_k=1$ such that $\sum_{k=n}^{N(n)} \gamma_k\mybar Y^k_t$
converges $\P$–a.s. to $\mybar Y_t$. Therefore, the linearity of $y\longmapsto \Lambda(t,y)$ gives the almost sure convergence of

\[
\sum_{k=n}^{N(n)} \gamma_k\mybar \alpha_t^k=  \Lambda\Big(t, \sum_{k=n}^{N(n)} \gamma_k\mybar Y_t^k\Big),
\]

which in turn implies, by the weak convergence of $\mybar \alpha^k_t$, that $\mybar \alpha_t=\Lambda(t, \mybar Y_t)$.
On the other hand, since $(\hat X^k,\hat Y^k,\hat Z^k)$ satisfies (<ref>), and $(\Omega, \|\cdot\|_\infty)$ is a Polish space, we have thanks to <cit.> that there exists measurable functions $\Phi^k:[0,T]\times C([0,T],\R^\xdim)\longrightarrow \R^{\xdim\times n}$, such that $\hat Z_t^k=\Phi^k(t,W_{\cdot \wedge t})$ on $(\Omega,\Fc)$ where $W_{\cdot \wedge t}$ denotes the path of $W$ up to time $t$.
Letting $\mybar Z_t^k:=\Phi^k(t,\mybar W^k_{\cdot \wedge t})$ we claim that $(\mybar X^k,\mybar Y^k,\mybar Z^k)$ satisfies (<ref>) on $(\mybar \Omega, \mybar \cF, \mybar \P)$.
Indeed, by construction $(\mybar X^k,\mybar Y^k , \mybar Z^k) \stackrel{d}{=} (\hat X^k,  \hat Y^k, \hat Z^k)$ meaning that their finite dimensional distributions coincide. Thus, by discretization we have that

\begin{align*}
& \mybar \E\bigg[\bigg\| \mybar Y_t^k-\mybar Y_0^k+\int_0^t F(r,\mybar X_r^k,\mybar Y_r^k,\Lc(\mybar X_r^k,\mybar Y_r^k))\diff r-\int_0^t \mybar Z_r^k \diff \mybar W_r^k\bigg\|^2\bigg]\\
&=\E \bigg[\bigg\| Y_t^k-Y_0^k+\int_0^t F(r,X_r^k,Y_r^k,\Lc(X_r^k,Y_r^k))\diff r-\int_0^t Z_r^k \diff W_r\bigg\|^2\bigg]=0.
\end{align*}

Consequently, we have that the sequence $(\mybar M^k)_{k\geq 1}$ given by

\[
\mybar M^k_t:= \mybar{Y_t}^{k}+\int_0^t\Big(  \partial_x H_2(r,\mybar X_r^k, \mybar Y_r^k ,\Lc(\mybar X_r^k))   +  \tilde \E \big[ \partial_\mu H_2 (r,\widetilde{ X}_r^k, \widetilde{Y}_r^k ,\Lc(\mybar X_r^k))(\mybar  X_r^k)  \big] \Big)   \diff r,
\]

is a $(\mybar \P,\mybar \F^k)$–martingale for every $k$, where $\widetilde{ X}^k$ and $\widetilde Y^k$ denote i.i.d. copies of $\mybar X^k$ and $\mybar Y^k$.

Step 2. We claim that for every $t\in [0,T]$, $(\mybar M^k_t)_{k\geq 1}$ converges weakly in $\L^1(\mybar\Omega)$ to

\[
\mybar M_t:= \mybar Y_t +\int_0^t \Big( \partial_x H_2(r,\mybar X_r,\mybar Y_r,\Lc(\mybar X_r))   +  \tilde \E \big[\partial_\mu H_2 (r,\widetilde{ X}_r, \widetilde{ Y}_r,\Lc(\mybar X_r))(\mybar X_r) \big] \Big)  \diff r.
\]

To show the claim, we prove that each of the summands converges weakly to its respective limit.

Let $\mybar K$ be a bounded random variable on $(\mybar \Omega, \mybar \Fc,\mybar \P)$ and note that there exists $C>0$ and $q>0$, $1/p+1/q=1$, such that

\begin{align*}
&\mybar \E\bigg[ K \int_0^t   \partial_x H_2(r,\mybar X_r^k, \mybar Y_r^k ,\Lc(\mybar X_r^k))  \diff r- K \int_0^t  \partial_x H_2(r,\mybar X_r,\mybar Y_r,\Lc(\mybar X_r))  \diff r\bigg]\\
& =\mybar \E \bigg[ K \int_0^t  \Big(\partial_x b(r,\mybar X_r^k,\Lc(\mybar X_r^k))- \partial_x b(r,\mybar X_r,\Lc(\mybar X_r))\Big)  \mybar Y_r^k\diff r \bigg]+ \mybar\E \bigg[  \int_0^t K \partial_x b(r,\mybar X_r,\Lc(\mybar X_r))   ( \mybar Y_r^k- \mybar Y_r ) \diff r\bigg] \\
&\;\; + \mybar\E \bigg[ K \int_0^t  \big(\partial_x f_2(r,\mybar X_r^k,\Lc(\mybar X_r^k))- \partial_x f_2 (r,\mybar X_r,\Lc(\mybar X_r))\big)   \diff r \bigg] \\
&\leq  C  \mybar\E  \bigg[\int_0^t  \| \partial_x b(r,\mybar X_r^k,\Lc(\mybar X_r^k))- \partial_x b(r,\mybar X_r,\Lc(\mybar X_r))\|^q\diff r\bigg]^{\frac1q} + \mybar\E \bigg[  \int_0^t K \partial_x b(r,\mybar X_r,\Lc(\mybar X_r))  ( \mybar Y_r^k- \mybar Y_r ) \diff r\bigg] \\
&\;\; + \mybar\E \bigg[ K \int_0^t  \big(\partial_x f_2(r,\mybar X_r^k,\Lc(\mybar X_r^k))- \partial_x f_2 (r,\mybar X_r,\Lc(\mybar X_r))\big)   \diff r \bigg],
\end{align*}

where the inequality follows from Cauchy-Schwartz, the boundedness of $K$, and the fact that $(\mybar Y^k)_{k\geq 1}$ is bounded in $\L^p([0,T]\times \mybar \Omega)$.
Now, recall that under <Ref><ref> the mapping $(x,\mu)\longmapsto b (t,x,\mu)$ is continuously differentiable and Lipschitz, and under the extra assumptions in $(ii)$ $(x,\mu)\longmapsto f_2 (t,x,\mu)$ is Lipschitz smooth uniformly in $t\in[0,T]$.
Thus, by dominated convergence, the almost sure convergence of $\mybar X^k$ to $\mybar X$ and of $\Lc(\mybar X_t^k)$ to $\Lc(\mybar X_t)$ for $t\in [0,T]$, implies that the first and third terms tend to zero.
The second term tends to zero in light of the boundedness of $\partial_x b$ and $H$ and the convergence of $\mybar Y^k$ to $\mybar Y$ weakly in $\L^1([0,T]\times \mybar \Omega)$.

We are left to argue for the last term in $\mybar M^k$. Again, for $\mybar K$ bounded, there exists $C>0$ and $q>0$, $1/p+1/q=1$, such that

\begin{align*}
&\mybar \E  \bigg[ K \int_0^t \widetilde{ \E} \big[ \partial_\mu    H_2(r,\widetilde{ X}_r^k, \widetilde{ Y}_r^k ,\Lc(\mybar X_r^k))(\mybar X_r^k) \big] \d r- K\int_0^t    \widetilde \E \big[\partial_\mu   H_2(r,\widetilde{ X}_r,\widetilde{ Y}_r,\Lc(\mybar  X_r))(\mybar X_r^k)\big] \d r\bigg]\\
&=   \mybar \E  \bigg[  \int_0^t K \widetilde \E \Big[ \big(  \partial_\mu b(r,\widetilde{ X}_r^k,\Lc(\mybar X_r^k))(\mybar X_r^k)- \partial_\mu b(r,\widetilde{ X}_r,\Lc(\mybar X_r))(\mybar X_r)\big)   \widetilde{ Y}_r^k  \Big]  \d r\bigg]  +\mybar\E \bigg[ K \int_0^t \widetilde \E\big[  \partial_\mu b(r,\widetilde{ X}_r,\Lc(\mybar X_r))(\mybar X_r)   ( \widetilde{ Y}_r^k- \widetilde{ Y}_r )\big] \d r\bigg]\\
&\;\;  +\mybar\E  \bigg[ K \int_0^t   \widetilde \E\big[ \partial_\mu f_2(r,\widetilde{ X}_r^k,\Lc(\mybar X_r^k))(\mybar X_r^k)- \partial_\mu f_2 (r,\widetilde{ X}_r,\Lc(\mybar X_r))(\mybar X_r) \big]  \d r \bigg] \\
&\leq  C \mybar\E \bigg[  \int_0^t  \widetilde \E \big[ \|   \partial_\mu b(r,\widetilde{ X}_r^k,\Lc(\mybar X_r^k))(\mybar X_r^k)- \partial_\mu b(r,\widetilde{ X}_r,\Lc(\mybar X_r))(\mybar X_r)\|^q \big] \d r\bigg]^{\frac{1}q}   +\widetilde \E\bigg[  \int_0^t  \mybar\E \big[K  \partial_\mu b(r,\widetilde{ X}_r,\Lc(\mybar X_r))(\mybar X_r)\big]   ( \widetilde{ Y}_r^k- \widetilde{ Y}_r ) \diff r\bigg]\\
&\;\;  + \mybar\E \bigg[ K \int_0^t   \widetilde \E\big[ \partial_\mu f_2(r,\widetilde{ X}_r^k,\Lc(\mybar X_r^k))(\mybar X_r^k)- \partial_\mu f_2 (r,\widetilde{ X}_r,\Lc(\mybar X_r))(\mybar X_r) \big]  \d r \bigg] ,
\end{align*}

where we used the boundedness of $K$ and of $(\widetilde Y^k)_{k\geq 1}$ in $\L^p([0,T]\times \widetilde \Omega)$ in the first term, and Fubini's theorem to exchange the order of the expectations in the second term.
Arguing as above, the convergence of the second term follows from that of $(\widetilde{Y}^k)_{k\geq 1}$, whereas the first and third terms exploit the convergence of $(\mybar X^k)_{k\geq 1}$, $(\widetilde X^k)_{k\geq 1}$ and $(\Lc(\mybar X_t^k))_{k\geq 1}$ for $t\in [0,T]$ with the assumptions on $b$ and $f_2$. With this, the claim is proved.

Step 3. We show that $\mybar M$ is a $(\mybar \P,\mybar \F)$-martingale for almost every $t\in [0,T]$. Let $0\leq s\leq t\leq T$ and $K$ be a bounded $\mybar \Fc_s$–measurable random variable and note that since $M^k$ is a $(\P,\F^k)$–martingale we have that

\begin{align*}
\mybar \E\big[ K\big( \mybar \E[\mybar M_t|\mybar\Fc_s]-\mybar M_s\big)\big]&=\mybar \E\big[ K\big( \mybar \E[\mybar M_t|\mybar\Fc_s]-\mybar \E[\mybar M_t|\mybar\Fc_s^k]\big)\big]+\mybar \E\big[ K\big( \mybar \E[\mybar M_t|\mybar\Fc_s^k]-\mybar \E[\mybar M_t^k|\mybar\Fc_s^k]\big)\big]+\mybar \E\big[ K\big( \mybar \E[\mybar M_t^k|\mybar\Fc_s^k]-\mybar M_s\big)\big]\\
&= \mybar \E\big[ K\big( \mybar \E[\mybar M_t|\mybar\Fc_s]-\mybar \E[\mybar M_t|\mybar\Fc_s^k]\big)\big]+\mybar \E\big[ K\big( \mybar \E[\mybar M_t- \mybar M_t^k|\mybar\Fc_s^k]\big)\big]+\mybar \E\big[ K\big( \mybar M_s^k -\mybar M_s\big)\big].
\end{align*}

Now, note that the convergence of $(\mybar M^k_t)_{k\geq 1}$ weakly in $\L^1(\mybar \Omega)$ for every $t\in [0,T]$, implies the last term goes to zero. Let us now look at the first term. Since $\mybar W^k$ converges to $\mybar W$, $\mybar \P\text{--a.s.}$, $\mybar W^k$ converges to $\mybar W$ in probability in the $J_1(\R^\bmdim)$ topology. Thus, by definition of $\mybar \F^k$ we have thanks to  <cit.> that $\mybar \F^k$ converges weakly to $\mybar \F$ in the sense of <cit.>, and consequently, $\mybar \E[\mybar M_t|\mybar\Fc_\cdot^k]$ converges to $\mybar \E[\mybar M_t|\mybar\Fc_\cdot]$ in probability. Since $\mybar\E[\|\mybar M_t\|^p]<C$ for almost every $t\in [0,T]$, by de la Vallée-Poussin's theorem, we have that the family $\mybar \E[\mybar M_t|\mybar\Fc_\cdot^k]$ is uniformly integrable, which in light of the convergence in probability implies that $\mybar \E[\mybar M_t|\mybar\Fc_\cdot^k]$ converges to $\mybar \E[\mybar M_t|\mybar\Fc_\cdot]$ in $\L^1(\mybar \Omega)$ and thus the first term goes to zero for almost every $t\in [0,T]$. We are only left to study the convergence of the second term. For this, we notice that since $\mybar \F^k$ converges weakly to $\mybar \F$, we have that $\mybar \Fc^k_s$ converges weakly to $\mybar \Fc^k_s$ in the sense of <cit.>, i.e., $K^k:=\E[K|\mybar \Fc^k_s]$ converges to $K$ in probability. Since $K$ is bounded, by de la Vallée-Poussin's theorem, $K^k$ converges to $K$ in $\L^1(\mybar \Omega)$, but, by the boundedness of $K$, also in $\L^q(\mybar \Omega)$, $1/q+1/p=1$. Thus,

\begin{align*}
\mybar \E\big[ K\big( \mybar \E[\mybar M_t- \mybar M_t^k|\mybar\Fc_s^k]\big)\big]&=\mybar \E\big[ (K-K^k)\big( \mybar \E[\mybar M_t- \mybar M_t^k|\mybar\Fc_s^k]\big)\big]+\mybar \E\big[ K^k \big( \mybar M_t- \mybar M_t^k\big)\big]\\
&\leq \|K-K^k\|_{\L^q(\overline \Omega)}\| \mybar M_t- \mybar M_t^k\|_{\L^p(\overline \Omega)}+\mybar \E\big[ K^k \big( \mybar M_t- \mybar M_t^k\big)\big],
\end{align*}

where we used the fact that $K^k$ is $\mybar \Fc^k_s$-measurable. Thus, the first term goes to zero since $\mybar\E[\|\mybar M_t\|^p],\mybar\E[\|\mybar M_t^k\|^p]<C$, for all $k\geq 1$ and the convergence of the sequence $(K^k)_{k\geq 1}$. The second term goes to zero in light of the convergence of $(\mybar M^k_t)_{k\geq 1}$ weakly in $\L^1(\mybar \Omega)$. In light of the continuity of $\mybar M$, we conclude that $\mybar M$ is a $(\mybar \P,\mybar \F)$-martingale establishing the last statement in $(ii)$. Since $\mybar\F=\F^{\overline W}$, by martingale representation, we find that there exists a predictable process $\mybar Z$ such that

\[
\mybar Y_t =\mybar Y_s +\int_t^s F(r,\mybar X_r,\mybar Y_r,\Lc(\mybar X_r,\mybar Y_r))   \diff r -\int _t^s \mybar Z_r\diff \mybar W_r,
\]

where we recall $F$ is given by (<ref>). This implies that $(\mybar X, \mybar Y, \mybar Z)$ solve the FBSDE (<ref>).

§.§ The entropic case: Proof of Theorem <ref>

At last, we bring ourselves back to the case of the entropic transportation cost problem presented in the introduction.

In the setting of the mean field entropic transportation cost problem we have that $f(t,x,a,\mu)=f_1(t,a)=\|a\|^2/2$, $H_1(t,y)=-|y|^2/2$, which readily verifies <Ref><ref> and the smooth Lipschitz assumption in <Ref><ref>, and $b(t,x,a,\mu)=a-\nabla \Psi \ast \mu (x)$ so that
the Lipschitz assumption on $\nabla\Psi$ guarantees that <Ref><ref> holds. In addition,  and $\Lambda(t, y)=-y$

\[
F(t,x,y,\nu) =-\int_{\RR^\xdim\times \R^m} \big(  \nabla^2 \Psi ( x -\tilde x)\cdot y -\nabla^2 \Psi (\tilde x-x) \cdot  \tilde  y \big) \nu (\d \tilde x\times \d \tilde y) ,
\]

the Lipschitz assumption on $\nabla^2\Psi$ guarantees that  <Ref><ref> holds. Thus, $(i)$ and $(ii)$ follow directly from <Ref><ref>.
Lastly, note that $(iii)$ is a direct consequence of the symmetry of $\Psi$.

§ FINITE PARTICLE APPROXIMATION OF THE MEAN FIELD SCHRÖDINGER PROBLEM

We now turn to the approximation of the mean field Schrödinger problem discussed in the previous section by the Schrödinger problem over finite, interacting particles.
Let us begin by introducing the generic finite particle problem we will consider.

Fix $N \in \mathbb{N}$.
Let $\cA$ be the set of admissible controls defined as

\begin{equation*}
\cA := \bigg\{\alpha:[0,T]\times \Omega \to A,\,\, \FF^N\text{--progressive processes such that } \EE\bigg[\int_0^Tf_1(t ,\alpha_t)  \diff t\bigg]<\infty\bigg\}.
\end{equation*}

We use the notation

\begin{equation*}
\balpha^N:= (\alpha^{1,N},\dots,\alpha^{N,N}),\quad \text{and}\quad \X^{N,\balpha}:= (X^{1,\alpha^N},\dots, X^{N,\alpha^N}),
\end{equation*}

for any control processes $\alpha^{i,N} \in \cA$, $i\in \{1,\dots,N\}$, and controlled diffusion

\begin{equation}
\label{eq:controlled.SDE.NP}
X^{i,\balpha^N}_t = X_0^{i} + \int_0^t\alpha_u^{i,N} + b(u, X^{i,\balpha^N}_u,L^N(\X^{N,\balpha}_u))\diff u + \sigma W^i_t,
\end{equation}

recall the notation $L^N(\x)$ for the empirical measure of the vector $\x$ as introduced at the end of <Ref>.
We consider the stochastic optimal transport problem

\begin{align}\label{eq:N.problem}
V^N(\muin, \mufin) := \inf\bigg\{\E\bigg[\frac1N\sum_{i=1}^N\int_0^Tf(t, X^{i,\balpha^N}_t, \alpha_t^{i,N},L^N(\X^{N,\balpha}_t))\diff t\bigg]: \balpha^N  \in \cA^N, \,   X_0^{i,\alpha^N}\sim \muin , \, X_T^{i,\alpha^N} \sim \mufin  , \,1\leq i\leq N  \bigg\}.
\end{align}

This problem can be seen as a natural generalization of the entropic stochastic optimal transport problem.
Here, we are optimally steering interacting particles with initial configuration $\muin$ to the terminal configuration $\mufin$.

Let us remark that, once again, appropriate choices of $f$ and $b$ make $V^N(\muin, \mufin)$ coincide with the value $\Vc^N_{\rm e}(\muin,\mufin)$ of (<ref>) and $\Vc_{\rm S}(\muin,\mufin)$ of (<ref>).

The main results of this section consist of convergence results that give conditions guaranteeing that the value of the finite particle problem converges to the value of the mean field Schrödinger problem, i.e., the convergence of the sequence $(V^{N}(\muin,\mufin))_{N\geq 1}$ to $V(\muin,\mufin)$.

These exploit additional convexity or extra regularity to establish the convergence for arbitrary finite time horizons.

Here again, the argument starts with the introduction of an appropriate penalization of the stochastic optimal transport problem (<ref>).
Given a smooth penalty function $g$ in the sense of <Ref>, we thus consider the penalized control problems

\begin{equation}
\label{eq:N.problem-k.penalized}
V^{N,k}(\muin, \mufin) := \inf \bigg\{ \frac1N\sum_{i=1}^N\EE\bigg[\int_0^Tf(t, X^{i,\alpha^N}_t\!\!, \alpha_t^{i,N},L^N(\X^{N, \balpha}_t))\diff t + kg(\mathcal{L}(X^{i,\alpha^N}_T)) \bigg]: \alpha^N \in \mathcal{A}^N,\,   X_0^{i,\alpha^N}\sim \muin ,\, 1\leq i\leq N\bigg\},
\end{equation}

where $X^{i,\alpha^N}$ is the unique strong solution of the controlled stochastic differential equation (<ref>).

The main ingredient for proving the convergence of the sequence $(V^{N, k}(\muin,\mufin))_{N\geq 1}$ to $V^{\infty, k}(\muin,\mufin)$ consists of a propagation of chaos result for coupled FBSDEs.
Indeed, in complete analogy to the analysis in <Ref>, the next lemma, whose proof is deferred to <Ref>, introduces the FBSDE system associated with Problem (<ref>).
To simplify the presentation, in the remainder of this section, we drop the superscript $k$ from the FBSDEs in the analysis unless otherwise stated.

Let $N,k\ge1$ be fixed.
Under the conditions of <Ref>, the control problem (<ref>) with value $V^{N,k}(\muin, \mufin)$ admits an optimal control $\hat\balpha^N:= (\alpha^{1,N},\dots, \alpha^{N,N})$ such that for each $i=1,\dots,N$ it holds

\begin{equation*}
\hat\alpha^{i,N}_t = \Lambda\big(t, Y^{i,N}_t \big), \quad \P\otimes \d t\text{--a.s.},
\end{equation*}

where $\Lambda$ is a function satisfying <Ref> and $(X^{i,N}, Y^{i,N},Z^{i,j,N} )$
solves the FBSDE

\begin{equation}
\label{eq:N-FBSDE-game.proof.main}
\begin{cases}
\diff X_t^{i,N} = B\big(t, X_t^{i,N}, Y^{i,N}_t, L^N(\X^{N}_t)\big)\diff t + \sigma \diff W^i_t\\
\diff Y^{i,N}_t = -F\big(t, X^{i,N}_t, Y_t^{i,N}, L^N(\X^{N}_t,\Y^{N}_t) \big)\diff t + \sum_{j=1}^NZ^{i,j,N}_t\diff W_t^j\\
X_0^{i,N} \sim \muin, \quad Y_T^{i,N} = k\partial_\mu g ( \cL(X^{i,N}_T))(X^{i,N}_T),
\end{cases}
\end{equation}

with the functions $B$ and $F$ defined respectively in <Ref> and <Ref>.

Recall from <Ref> that for each $k \ge1$, an optimizer $\hat\alpha$ of the control problem (<ref>) with value $V^{\infty, k}(\muin,\mufin)$ is given by $\hat\alpha = \Lambda (t ,  Y_t)$ where $( X,  Y,  Z)$ solves the McKean–Vlasov equation

\begin{equation}
\label{eq:fbsde.reg.MFSP.proof}
\begin{cases}
\diff X_t = B(t, X_t, Y_t, \cL( X_t))\diff t + \sigma \diff W_t\\
\diff  Y_t = -F(t,  X_t, Y_t, \cL(X_t,Y_t)) \diff t + Z_t\diff W_t\\
X_0\sim \muin,\quad Y_T = k\partial_\mu	g(\cL(X_T))( X_T).
\end{cases}
\end{equation}

The goal now is to show that the particle system (<ref>) evolving forward and backward in time converges to the McKean–Vlasov system (<ref>) in an appropriate sense.
This propagation of chaos would follow from results developed in [Laurièrie and Tangpi, 2022, Peng and Tangpi, 2022] if not for the fact that, for $k$ fixed, the $N$-particle system is actually a McKean–Vlasov FBSDE itself, since the terminal value of the value processes $Y^{i,N}$ depends on the law of the forward process.
We thus need new arguments to obtain convergence in the present case.
Indeed, establishing this result (for $k$ fixed) dictates the assumptions behind our finite particle approximation results <Ref> and <Ref> below.
Moreover, when we bring ourselves back to the convergence of $V^N(\muin,\mufin)_{N\geq 1}$ the assumptions under which the previous propagation of chaos result holds should not depend on the value of $k$.
Naturally, the analysis comes with extra requirements on the terminal condition, which, as a consequence of our penalization approach, will lead us to tighten to choice of penalty functions.
We elaborate more on this below.

§.§ Finite particle approximation under displacement convexity

Our first result exploits the notion of displacement convexity, which we present next, to show the finite particle approximation.

An $L$-differentiable function $\phi:\Pc(\R^\xdim)\longmapsto \R$ is said to be displacement convex if there is $\lambda_\phi\geq 0$ such that for any two $\mu_0,\mu_1\in \Pc_2(\R^\xdim)$ and $\mu\in \Gamma(\mu_0,\mu_1)$ we have that

\[
\int_{\R^\xdim \times \R^\xdim} (x-y) \cdot \big( \partial_\mu \phi(\mu_1)(x)-\partial_\mu \phi(\mu_0)(y) \big)\mu(\d x, \d y )\geq \lambda_\phi \int_{\R^\xdim \times \R^\xdim} \| x-y\|^2\mu(\d x, \d y).
\]

Let <Ref> hold with the strengthening that $f_1(t,\cdot)$ is strongly convex uniformly in $t\in [0,T]$.
Suppose that $g$ is a displacement convex penalty function.
Then we have that

\begin{equation*}
\lim_{N\to \infty}V^N(\muin,\mufin) = V (\muin,\mufin).
\end{equation*}

Let us comment on the additional assumptions in <Ref>, namely, $(i)$ $g$ is a displacement convex penalty function; and $(ii)$ the strong convexity of $f_1(t,\cdot)$.
With regards to $(i)$, displacement convexity, which is strongly linked to the notion of displacement monotonicity, is known to play a relevant role in the analysis of mean field problems. In particular, when $\mufin=\delta_{x_0}$ and we take the squared Wasserstein–2 distance as a penalty function, i.e.,

\[
g_2:\cP_2(\RR^\xdim)\longrightarrow [0,\infty),\; \mu \longmapsto \mathcal{W}_2^2(\mu , \delta_{x_0}),
\]

it follows from <cit.>, see also Remark 3.3 therein, that $g_2$ is $L$-differentiable and for any $(x,\mu) \in \R^\xdim\times \Pc_2(\R^\xdim)$

\[
\partial_\mu g_2(\mu)(x)=2(x-x_0).
\]

Consequently, $g_2$ is a displacement convex penalty function with $\lambda_{g_2}=2$.
Recall that as mentioned in <Ref><ref> this example is still meaningful from a practical point of view.
As for $(ii)$, recall if $f_1(t,\cdot)$ is strongly convex, first, the map $\Lambda(t,y)$ is unique, and second, since (<ref>) is equivalent to $\partial_a f_1(t,\Lambda(t,y))=-y$, the strong convexity of $f_1(t,\cdot)$[This is, there exists $\lambda_{f_1}>0$ such that for any $(t,a,a^\prime)\in [0,T]\times A^2$, we have that
f_1(t,a)-f_1(t,a^\prime)-(a-a^\prime)\cdot \partial_a f_1(t,a) \geq \lambda_{f_1} \|a-a^\prime\|^2,
$] is equivalent to

\begin{align}\label{eq:displaassumpLambda}
(y - y^\prime )\cdot(\Lambda(t,y) - \Lambda(t,y^\prime)) \le - \lambda_{f_1} \|y-y^\prime\|^2.
\end{align}

for any $(t,y,y^\prime)\in [0,T]\times (\R^\xdim)^2$. A typical example is given by the entropic case, see <Ref>.

As mentioned already, the crux of the finite particle approximation lies in the following propagation of chaos result. For this, we introduce the space $(\Hc^2,\|\cdot\|_{\Hc^2})$ of adapted $\R^\xdim\times\R^\xdim$–valued processes $(X,Y)$ equipped with the norm

\[
\| (X,Y)\|_{\Hc^2}^2:=\E\bigg[\sup_{t\in [0,T]} \| X _t \|^2 \bigg]+\E\bigg[\int_0^T\|Y_t\|^2\diff t \bigg].
\]

Let <Ref> hold with the strengthening that $f_1(t,\cdot)$ is strongly convex uniformly in $t\in [0,T]$. Suppose that $g$ is a displacement convex penalty function.
Let $(  X^{i,N},  Y^{i,N})$ be part of the solution to (<ref>), and let $( \widetilde X^i, \widetilde Y^i)$, $i=1,\dots,N$, denote i.i.d. copies of $( X, Y, Z)$ solution to (<ref>). Then, there is $C>0$ such that for $i=1,\dots,N$, we have that

\begin{align}\label{eq.propagationresults}
\| (X^{i,N}-\widetilde X^i,Y^{i,N}-\widetilde Y^i)\|_{\Hc^2}^2  \leq C \E \big[\epsilon^N\big], \text{ where, } \epsilon^N:= \int_0^T  \cW_1^2\big(L^N(\widetilde \X^N_t, \widetilde \Y^N_t), \cL(X_t, Y_t)\big) \diff t .
\end{align}

In particular, $( X^{i,N},Y^{i,N})$ converges to $(  \widetilde X^i,  \widetilde Y^i)$ in $\Hc^2$.

Let $( \widetilde X^i, \widetilde Y^i, \widetilde Z^{i})$, $i=1,\dots,N$ be i.i.d. copies of $( X, Y, Z)$ solving the McKean-Vlasov FBSDE (<ref>).

Let $\delta X^i_t:= X^{i,N}_t-\widetilde X^i_t$ and $\delta Y^i_t := Y^{i,N}_t - \widetilde Y^{i}_t$.
Applying Itô's formula, the dynamics of $\delta X^i_t\cdot \delta Y^i_t$ are, for a martingale $M^i$, given by

\begin{align*}
\delta X^i_T \cdot \delta Y^i_T- \delta X^i_0 \cdot \delta Y^i_0 &= \int_0^T \delta X^i_t\cdot\Big(-F\big(t, X^{i,N}_t, Y^{i,N}_t, L^N(\X_t^N, \Y^N_t)\big) + F\big(t, \widetilde X^i_t, \widetilde Y^i_t, \cL(X_t, Y_t)\big)\Big)\d t\\
&\quad  +\int_0^T  \delta Y^i_t\cdot\Big(B\big(t, X^{i,N}_t, Y^{i,N}_t, L^N(\XX_t^N)\big) - B\big(t, \widetilde X^i_t, \widetilde Y^i_t, \cL(X_t)\big)\Big) \diff t + M_T^i
\end{align*}

Observe that by displacement convexity of $g$, we have

\begin{align*}
\E[\delta X^i_T \cdot \delta Y^i_T] &= k\lambda_g \E\Big[\delta X^i_T\cdot \big( \partial_\mu g(\cL(X^{N,i}_T))(X^{N,i}_T) - \partial_\mu g(\cL(X_T))(X_T)\big) \Big]\\
&\ge k\lambda_g  \E\big[\|\delta X_T^i\|^2\big].
\end{align*}

Hence, using Lipschitz–continuity of $F$ and $b$ (and recalling that $B(t,x,y, \mu) =\Lambda (t,y) + b(t,x,\mu)$) and then Young's inequality, for every $\varepsilon>0$ it holds that

\begin{align*}
k\lambda_g  \E\big[\|\delta X^i_T\|^2\big] &\le \E\bigg[\int_0^T \ell_F \|\delta X_t^i\|\Big( \|\delta X_t^i\|+\|\delta Y_t^i\|+\frac{1}N\sum_{j=1}^N \|\delta X_t^j\|+\frac{1}N\sum_{j=1}^N \|\delta Y_t^j\| + \cW_1\big(L^N(\widetilde \X^N_t, \widetilde \Y^N_t), \cL(X_t, Y_t)\big) \Big)\d t\bigg]\\
&\quad + \E\bigg[\int_0^T  \ell_b \|\delta Y_t^i\|\Big( \|\delta X_t^i\|+\frac{1}N\sum_{j=1}^N \|\delta X_t^j\| + \cW_1\big(L^N(\widetilde \X^N_t), \cL(X_t)\big)\Big)\d t\bigg]+  \E\bigg[\int_0^T\delta Y^i_t\cdot\big( \Lambda(t, Y^{i,N}_t) - \Lambda (t, \widetilde Y^i_t) \big)\diff t \bigg]\\
&\le  \E\bigg[\int_0^TC_1^\eps \|\delta X^i_t\|^2 + C_2^\eps \frac1N\sum_{j=1}^N\|\delta X^j_t\|^2\diff t \bigg] + \varepsilon\E\bigg[ \int_0^T\|\delta Y^i_t\|^2 + \frac1N\sum_{j=1}^N\|\delta Y^j_t\|^2\diff t \bigg] - \lambda_{f_1} \E\bigg[\int_0^T\|\delta Y^i_t\|^2\diff t \bigg]\\
&\quad + \E\bigg[\int_0^T\eps \cW_1^2\big(L^N(\widetilde \X^N_t, \widetilde \Y^N_t), \cL(X_t, Y_t)\big) + \frac{\ell_b^2}{\eps} \cW_1^2\big(L^N(\widetilde \X^N_t), \cL(X_t)\big)\diff t \bigg]
\end{align*}

where $C_1^\eps:=\ell_F + (2 \ell_F^2+\ell_b^2)\varepsilon^{-1}$, $C_2^\eps:=\eps \ell_F + \ell_b^2\varepsilon^{-1}$, in the first inequality we use triangular inequality for the Wasserstein distance, and in the second inequality we use (<ref>).

Thus, rearranging terms, picking $\varepsilon<\lambda_{f_1} /2$ and letting $\gamma:=\lambda_{f_1} -2\varepsilon>0$ we obtain

\begin{align}
\notag
\gamma \E\bigg[\int_0^T\|\delta Y^i_t\|^2\diff t \bigg] &\le \E\bigg[\int_0^TC_1^\eps \|\delta X^i_t\|^2 + C_2^\eps \frac1N\sum_{j=1}^N\|\delta X^j_t\|^2\diff t \bigg] + \varepsilon\E\bigg[ \int_0^T \frac1N\sum_{j=1}^N\|\delta Y^j_t\|^2\diff t \bigg]\\
&\quad + \E\bigg[\int_0^T\eps \cW_1^2\big(L^N(\widetilde \X^N_t, \widetilde \Y^N_t), \cL(X_t, Y_t)\big) +\frac{\ell_b^2}{\eps}   \cW_1^2\big(L^N(\widetilde \X^N_t), \cL(X_t)\big)\diff t \bigg].\label{eq:estim1.disp.conv1}
\end{align}

Now, applying Itô's formula to $\|\delta X^i_t\|^2$ and Young's inequality, it follows by the Lipschitz–continuity of $B$ that

\begin{align}\label{eq:estim3.displ.conv}
\|\delta X^i_t\|^2 &= \int_0^t2\delta X^i_s\cdot\Big( B(s, X^{i,N}_s, Y^{i,N}_s, L^N(\XX_s^N) - B(t, \widetilde X^i_s, \widetilde Y^i_s, \cL(X_s))\Big) \diff s\notag \\
&\le \ell_B\int_0^t5 \|\delta X^i_s\|^2 + \|\delta Y^i_s\|^2 + \frac1N\sum_{j=1}^N\|\delta X^j_s\|^2 + \cW^2_1\big(L^N(\XX^N_s), \cL(X_s)\big)\diff s.
\end{align}

Applying Gronwall's inequality (resp. averaging and then applying Gronwall's inequality) we find that

\begin{align}
\E\big[ \|\delta X^i_t\|^2\big] & \le \e^{5\ell_B t}\E\bigg[\int_0^t \|\delta Y^i_s\|^2+ \frac1N\sum_{j=1}^N\|\delta X^j_s\|^2  + \cW^2_1\big(L^N(\widetilde\XX^N_s), \cL(X_s)\big)\diff s\bigg]  \label{eq:estim5.displ.conv1}\\
\E\Big[\frac1N\sum_{j=1}^N\|\delta X^j_t\|^2\Big] &\le \e^{6 \ell_B t}\E\bigg[\int_0^t\frac1N\sum_{j=1}^N\|\delta Y^j_s\|^2 + \cW^2_1\big(L^N(\widetilde\XX^N_s), \cL(X_s)\big)\diff s\bigg]  .\label{eq:estim4.displ.conv1}
\end{align}

Let us now average out on both sides of (<ref>) and use (<ref>) to obtain

\begin{align*}
( \gamma-\eps) \E\bigg[\int_0^T \frac1N\sum_{i=1}^N \|\delta Y^i_t\|^2\diff t \bigg] &\le   \e^{5\ell_B T} ( C_1^\eps + C_2^\eps)  \E\bigg[ \int_0^T  \int_0^t  \frac1N\sum_{j=1}^N\|\delta Y^j_s\|^2  \diff s \diff t \bigg]    \\
&\quad + \E\bigg[\int_0^T\eps \cW_1^2\big(L^N(\widetilde \X^N_t, \widetilde \Y^N_t), \cL(X_t, Y_t)\big) +\big( T \e^{5\ell_B T} ( C_1^\eps + C_2^\eps)+ \frac{\ell_b^2}{\eps}\big)   \cW_1^2\big(L^N(\widetilde \X^N_t), \cL(X_t)\big)\diff t \bigg].
\end{align*}

It then follows from Gronwall's inequality, updating $\eps$ so that $\gamma-\eps>0$ and with $\epsilon^N$ as in (<ref>), that there is $C>0$ such that

\begin{align}\label{eq:estim2.displ.conv}
\E\bigg[\int_0^T \frac1N\sum_{i=1}^N \|\delta Y^i_t\|^2\diff t \bigg] &\le   C \E\big[\epsilon^N\big] ,
\end{align}

where we used the fact that  $ \cW_1\big(L^N(\widetilde \X^N_t), \cL(X_t)\big)\leq \cW_1\big(L^N(\widetilde \X^N_t, \widetilde \Y^N_t), \cL(X_t, Y_t)\big)$. Consequently, back in (<ref>) we obtain that there is $C>0$ such that

\begin{align}\label{eq:estim5.displ.conv}
\E\bigg[\frac1N\sum_{j=1}^N\|\delta X^j_t\|^2\bigg]&\le   C \E\big[\epsilon^N\big].
\end{align}

Now we use (<ref>) back in (<ref>) to derive that there is $C>0$ such that

\begin{align*}
\gamma \E\bigg[\int_0^T\|\delta Y^i_t\|^2\diff t \bigg] &\le \E\bigg[\int_0^TC_1^\eps  \e^{5\ell_B t} \int_0^t \|\delta Y^i_s\|^2\diff s   + C \frac1N\sum_{j=1}^N\|\delta X^j_t\|^2\diff t \bigg] + \varepsilon\E\bigg[ \int_0^T \frac1N\sum_{j=1}^N\|\delta Y^j_t\|^2\diff t \bigg] + C \E\big[\epsilon^N\big],
\end{align*}

which after using Grownwall's inequality shows, in light of (<ref>) and (<ref>), that there is $C>0$ such that

\begin{align*}
\E\bigg[\int_0^T\|\delta Y^i_t\|^2\diff t \bigg] &\le  C \E\big[\epsilon^N\big].
\end{align*}

Going back to (<ref>) now allows to conclude that $\E\big[\sup_{t\in [0,T]} \|\delta X^i_t\|^2 \big]$ satisfies the same estimate.
This shows (<ref>).
The second assertion follows since by the law of large numbers, see e.g. <cit.>, $ \E[\epsilon^N] \longrightarrow 0$, as $N\longrightarrow \infty$.

Having established the convergence of the associated FBSDEs, we have all the elements to present the proof of the finite particle approximation result.

We argue in three steps.

Step 1: We show that for $k\geq 0$ fixed we have that

\[
V^{N,k}(\muin, \mufin) \longrightarrow V^{\infty,k}(\muin, \mufin), \text{ as } N\longrightarrow \infty.
\]

Note that in light of Lemmata <ref> and <ref>, $\hat\alpha_t = \Lambda(t, Y_t)$ is optimal for the control problem (<ref>) with value $V^{\infty,k}:=V^{\infty, k}(\muin,\mufin)$ and,  by uniquenes of $\Lambda$, so is $\hat\balpha^{N}=(\hat\alpha^{1,N},\dots,\hat\alpha^{N,N})$ with $\hat\alpha^{i,N}_t = \Lambda(t,   Y^{i,N}_t)$ for the control problem (<ref>) with value $V^{N,k}:=V^{N, k}(\muin,\mufin)$, i.e.,

\begin{align*}
V^{N,k}  = \frac1N\sum_{i=1}^N\EE\bigg[\int_0^Tf(t,  X^{i,N}_t, \hat \alpha_t^{i,N},L^N( \X^{N}_t))\diff t + kg(\cL(X^{i,N}_T)) \bigg],\;
V^{\infty,k}  = \EE\bigg[\int_0^Tf(t,  X_t, \hat \alpha_t,\cL( X_t))\diff t + kg(\cL(X_T)) \bigg].
\end{align*}

Thus, since $f$ is locally Lipschitz-continuous, there exists a linearly growing function $\beta$ such that

\begin{align*}
|V^{N,k} -V^{\infty,k} |& \leq \frac1N\sum_{i=1}^N\E\bigg[ \int_0^T \beta( t,X^{i,N}_t, X_t , \hat \alpha^{i,N} ,  \hat \alpha_t) \big( \|  X^{i,N}_t- X_t \|+  \|\hat \alpha^{i,N}_t - \hat \alpha_t\| +\Wc_2 (L^N( \X^{N}_t) , \mathcal{L}(X_t))\big) \d t\bigg]\\
&\; +k \bigg( \sup_{i\leq N}  g(\cL(X^{i,N}_T)) -g(\mathcal{L}(X_T))\bigg).
\end{align*}

We now note that under <Ref> there exists $C>0$, independent of $N$, such that

\[
M:=\sup_{i\leq N} \E \Big[ \sup_{t\in [0,T]}\Big \{ 1+\|X^{i,N}_t\|^2+ \|Y^{i,N}_t\|^2+ \| X_t\|^2+\| Y_t\|^2\Big \} \Big]\leq C.
\]

Thus, using the fact that $\Lambda$ is Lipschitz–continuous and applying Cauchy–Schwarz inequality, we find that for a constant $C>0$ independent of $N$ it holds that

\begin{align*}
|V^{N,k} - V^{\infty,k} | \leq &\, C   \frac1N\sum_{i=1}^N \E \bigg[\int_0^T(1 + \ell_\Lambda)\|X^{i,N}_t - X_t\|^2 + \ell_\Lambda \|Y^{i,N}_t - Y_t\|^2 + (1 + \ell_\Lambda)\Wc_2^2(L^N(\X^N_t), \cL(X_t)) \Big) \diff t\bigg]^{\frac{1}2}\\
& +k  \bigg( \sup_{i\leq N}  g(\cL(X^{i,N}_T)) -g(\mathcal{L}(X_T))\bigg).
\end{align*}

The claim now follows since under the assumptions of <Ref> we have that <Ref> holds, and consequently, $(X^{i,N},Y^{i,N})$ converges to $(X,Y)$ in $\Hc^2$ as $N\longrightarrow \infty$.

Step 2: We show that

\[
\limsup_{N\to \infty} V^N (\muin,\mufin) \leq V (\muin,\mufin).
\]

Let $\alpha\in \mathfrak{A}$, and let $\tilde \alpha^{i}$, $i=1,\dots,N$, be i.i.d. copies such that $\tilde \alpha^{i}$ is $\F^i$-adapted and $X_0^{i,N,\tilde \alpha^i}\sim \muin$, $X_T^{i,N,\tilde \alpha^i}\sim \mufin$, for $X^{i,N,\tilde \alpha^i}$ solution to (<ref>), i.e.,

\[
X^{i,N,\tilde \alpha^i}_t = X_0^{i,N,\tilde \alpha^i} + \int_0^t\tilde \alpha_u^{i} + b(u, X^{i,N,\tilde \alpha^i}_u,L^N(\X^{N,\tilde \alpha^N}_u))\diff u + \sigma W^i_t,
\]

where, as before, $\X^{N,  \tilde \alpha^N }:=(X^{1,N,\tilde \alpha^1},\dots,X^{N,N,\tilde \alpha^N})$ and $ \tilde \alpha^N :=(\tilde \alpha^{1},\dots,\tilde \alpha^{N})$. Note that $ \tilde \alpha^N \in \mathcal{A}^N$. It is then clear that

\begin{align}\label{eq:prop.con.aux1}
V^N (\muin,\mufin) \leq \frac1N\sum_{i=1}^N\EE\bigg[\int_0^Tf(t, X^{i,N, \tilde \alpha^i}_t, \tilde \alpha_t^{i},L^N(\X^{N, \tilde \alpha^N}_t))\diff t  \bigg].
\end{align}

Now, let $X^{i,\tilde \alpha^{i}}$ denote the solution to

\[
X^{i,\tilde \alpha^i}_t = X_0^{i,\tilde \alpha^i}+ \int_0^t\tilde \alpha_u^{i} + b(u, X^{i,\tilde \alpha^i}_u,\Lc(X^{i,\tilde \alpha^i}_u))\diff u + \sigma W^i_t,\; X_0^{i,\tilde \alpha^i}\sim \muin,
\]

and note that by a standard propagation of chaos argument $L^N(\X^{N,\tilde \alpha^N })$ converges to $\Lc(X^{i,\tilde \alpha^{i}})$. In particular, $X_T^{i,\tilde \alpha^i}\sim \mufin$. Consequently, the estimate

\[
\E \bigg[\sup_{t\in[0,T]} \| X^{i,N,\tilde \alpha^{i}}_t-X^{i,\tilde \alpha^{i}}_t\|\bigg] \leq C \int_0^T \E \big[   \Wc_1(L^N(\X^{N,{\tilde \alpha}^N }_r) , \mathcal{L}(X^{i,\tilde  \alpha^{i}}_r)) \big] \diff r,
\]

which holds for some $C>0$, implies that $\Lc(X^{i,N,\tilde \alpha^{i}})$ converges to $\Lc(X^{i,\tilde \alpha^{i}})$. Thus, by the continuity of $f$, the term on the right in (<ref>) converges and thus

\[\limsup_{N\to \infty} V^N (\muin,\mufin)\leq
\EE\bigg[\int_0^Tf(t, X^{i, \tilde \alpha^{i}}_t, \tilde \alpha_t^{i},\mathcal{L}(X^{i,\tilde \alpha^{i}}_t))\,d t   \bigg].
\]

Since $\tilde \alpha^{i}\sim\alpha \in \mathfrak A$, $X_0^{i,\tilde \alpha^{i}}\sim \muin$ and $X_T^{i,\tilde \alpha^{i}}\sim \mufin$ the claim follows.

Step 3: We show that

\[
\liminf_{N\to \infty} V^N (\muin,\mufin) \geq V (\muin,\mufin).
\]

Recall that by weak duality we have that $V^N (\muin,\mufin) -V^{N,k} (\muin,\mufin)\geq 0$ for any $k\geq 0$ and thus

\begin{align*}
V^N (\muin,\mufin) - V (\muin,\mufin)& =  V^N (\muin,\mufin) -V^{N,k} (\muin,\mufin) +V^{N,k} (\muin,\mufin) -V^{\infty, k} (\muin,\mufin)+V^{\infty, k} (\muin,\mufin)- V (\muin,\mufin)\\
&\geq V^{N,k} (\muin,\mufin) -V^{\infty, k} (\muin,\mufin)+V^{\infty,k} (\muin,\mufin)- V (\muin,\mufin)
\end{align*}

The claim follows since the first term on the right goes to $0$ by Step 1 and the second term vanishes by <Ref>. Step 2 and Step 3 establish the result.

We are left to show the convergence of the Shrödinger Bridges. We now make explicit the dependence on $k$ and $N$. This is $\hat X^{k,N}$, $\hat X^{k}$, $\hat X$, denote the optimal controlled processes associated to the problems (<ref>), (<ref>), (<ref>). Let us also introduce the process $\hat X^{N}$, obtained letting $k\longrightarrow \infty$ in (<ref>). What's more, arguing as in <Ref> and <Ref>, we can show that $\hat X^{k,N}$ converges in law to $\hat X^N$, as $k\longrightarrow \infty$. The convergence now follows from

\[
\hat X^{N}- \hat X= \hat X^{N}-\hat X^{k,N}+\hat X^{k,N} - \hat X^{k}+\hat X^{k}-\hat X,
\]

since by <Ref>, for every $k\geq 1$, we have that $\hat X^{k,N}$ converges in law to $\hat X^k$, as $N\longrightarrow \infty$, for $\hat X^{k,N}$ and $\hat X^k$ given by <Ref> and <Ref> respectively. Moreover, by <Ref> so does $\hat X^k$ to $\hat X$, as $k\longrightarrow \infty$. Lastly, arguing as in <Ref> and <Ref> we obtain that $\hat X^{k,N}$ converges in law to $\hat X^N$, as $k\longrightarrow \infty$.

First observe that as a direct consequence of the propagation of chaos <Ref> we have $V^{N, k,n}(\muin, \mufin)\rightarrow V^{\infty, k,n}(\muin, \mufin)$ uniformly in $n$.
In fact,
By optimality of $\hat\balpha^{N}$ with $\hat\alpha^{i,N}_t = \Lambda(t, X^{i,N}_t, Y^{i,N}_t, L^N(\XX^N_t))$ for the control problem (<ref>) with value $V^{N, k,n}(\muin,\mufin)$ and optimality of $\hat\alpha^i_t = \Lambda(t, X_t, Y_t, \cL(X_t))$ for the control problem  (<ref>) with value $V^{\infty, k,n}(\muin,\mufin)$, we have

\begin{equation*}
\begin{cases}
V^{N,k,n}(\muin, \mufin) = \frac1N\sum_{j=1}^N\EE\bigg[\int_0^Tf(t,  X^{j,N}_t, \hat \alpha_t^{j,N},L^N( \X^{N}_t))\diff t + kg^n(\cL(X^{j,N}_T)) \bigg]\\
V^{\infty,k,n}(\muin, \mufin) = \EE\bigg[\int_0^Tf(t,  X_t, \hat \alpha_t,\cL( X_t))\diff t + kg^n(\cL(X_T)) \bigg].
\end{cases}
\end{equation*}

Thus, since $g^n$ is $1$–Lipschitz for all $n$ and $f$ is locally Lipschitz–continuous, we have

\begin{align*}
|V^{N,k,n}(\muin, \mufin) - V^{\infty,k,n}(\muin, \mufin) | &\le \frac1N\sum_{j=1}^N\EE\bigg[\int_0^T\beta(X^{j,N}_t, X_t, \hat\alpha^{i,N}_t, \hat\alpha_t)\Big(\|X^{j,N}_t - X_t\| + \|\hat\alpha^{j,N}_t - \hat\alpha_t\| + \W_2(L^N(\X^N_t), \cL(X_t)) \Big) \diff t\bigg]
\end{align*}

for a linearly growing function $\beta$.
Thus, using the fact that $\Lambda$ is Lipschitz–continuous and applying Cauchy–Schwarz inequality yields

\begin{align}
\notag
&\sup_{n\ge1}|V^{N,k,n}(\muin, \mufin) - V^{\infty,k,n}(\muin, \mufin) |\\
\notag
&\le C\EE\Big[\sup_{t\in [0,T]}\|X^{j,N}_t\|^2+ \| X_t\|^2+1\Big]^{1/2}\EE\bigg[\int_0^T(1 + \ell_\Lambda)\|X^{j,N}_t - X_t\|^2 + \ell_\Lambda|Y^{i,N}_t - Y_t\| + (1 + \ell_\Lambda)\W_2^2(L^N(\X^N_t), \cL(X_t)) \Big) \diff t\bigg]^{1/2} \\
\label{eq.conv.V.N}
&\quad \xrightarrow [N\to \infty]{}0
\end{align}

which converges to zero by <Ref>.

On the other hand, by Lemma <ref>, the control problem (<ref>) with value $V^{\infty, k,n}(\muin,\mufin)$ admits an optimizer $\hat\alpha^{k,n}$ satisfying $\hat\alpha^{k,n} = \Lambda_t\big(\hat X^{k,n}_t, \hat Y^{k,n}_t,\cL(\hat X^{k,n}_t)\big)$ with $(\hat X^{k,n}, \hat Y^{k,n}, \hat Z^{k,n})$ solving the McKean–Vlasov equation (<ref>) and therefore

\begin{equation*}
V^{\infty,k,n}(\muin, \mufin) = \EE\bigg[\int_0^Tf(t,  \hat X^{k,n}_t, \hat \alpha_t^{k,n},\cL( \hat X^{k,n}_t))\diff t + kg^n(\cL(\hat X^{k,n}_T)) \bigg].
\end{equation*}

To prove this theorem it suffices to put the above lemmas together.
In fact, since $V^N(\muin, \mufin) \ge V^{N,k}(\muin, \mufin)$ for all $k,N$, we have

\begin{align}
\notag
|V^N(\muin, \mufin) - V(\muin, \mufin)| &=  |V^N(\muin, \mufin) - V^{N,k}(\muin, \mufin)| + |V^{N,k}(\muin, \mufin) -  V(\muin, \mufin)|\\
\label{eq:estim.V.conv.proof}
&\le \big(V^N(\muin, \mufin) - V^{N,k}(\muin, \mufin)\big)^++ |V^{N,k}(\muin, \mufin) -  V(\muin, \mufin)|.
\end{align}

We now show that the right hand side above converges to zero.

There also exists $\balpha^k = (\alpha^{1,k},\dots, \alpha^{N,k}) \in \cA^N$ such that

\begin{equation*}
V^{N,k}(\muin, \mufin) \ge \frac1N\sum_{i=1}^N\EE\bigg[\int_0^Tf(t, X^{i,\alpha^{i,k}}_t, \alpha_t^{i,k},L^N(\X^{N,\balpha^k}_t))\diff t + kg^n(\cL(X^{i,\alpha^{i,k}}_T)) \bigg] - \varepsilon.
\end{equation*}

Arguing exactly as in the proof of Lemma <ref> allows to construct a weak limit $(\int_0^\cdot \alpha^1_s\diff s, \dots, \int_0^\cdot \alpha^N_s\diff s)$ of $(\int_0^\cdot \alpha^{1,k}_s\diff s, \dots, \int_0^\cdot \alpha^{N,k}_s\diff s)$ such that
$X^{i,\alpha^{i,k}}$ converges to $X^{i, \alpha^i}$ (on a possibly new probability space) and

\begin{align*}
\liminf_{k\to \infty}V^{N,k}(\muin,\mufin)&\ge
\liminf_{k\to \infty}\frac1N\sum_{i=1}^N\EE\bigg[\int_0^Tf(t, X^{i,\alpha^{i,k}}_t, \alpha_t^{i,k},L^N(\X^{N,\balpha^k}_t))\diff t + kg^n(\cL(X^{i,\alpha^{i,k}}_T)) \bigg] -\varepsilon\\
&\ge V^N(\muin, \mufin) - \varepsilon.%\frac1N\sum_{i=1}^N\EE\bigg[\int_0^Tf(t, X^{i,\alpha^i}_t, \alpha^i_t, L^N(\X_t^{N,\balpha})\diff t  \bigg].
\end{align*}

Thus, we have

\begin{align*}
\notag
\limsup_{k\to \infty}\sup_{N\in \NN}\big(V^N(\muin, \mufin) - V(\muin,\mufin)\big)^+ &\le \sup_{N\ge1}\Big( V^N(\muin,\mufin) - \liminf_{k\to \infty}V^{N,k}(\muin, \mufin) \Big)\\
\label{eq:conv.V.k}
&\le \varepsilon.
\end{align*}

Therefore, coming back to (<ref>), we have

\begin{align*}
\big|V^{N}(\muin, \mufin) -  V(\muin, \mufin)\big| &\le \sup_{N\in \NN}\big(V^{N}(\muin, \mufin) -  V^{N,k}(\muin, \mufin)\big) + \big|V^{N,k}(\muin, \mufin) - V^{N,k,n}(\muin, \mufin)\big|\\
&\quad + \sup_{n\ge1}\big|V^{N,k,n}(\muin, \mufin) - V^{\infty, k,n}(\muin, \mufin)\big| \\
&\quad + \big|V^{\infty,k,n}(\muin, \mufin) - V^{\infty, k}(\muin, \mufin)\big| + \big|V^{\infty,k}(\muin, \mufin) - V(\muin, \mufin)\big|.
\end{align*}

Subsequently taking the limit in $n, N$ and $k$, it follows from Lemma <ref> and <Ref> that

\begin{equation*}
\lim_{N\to \infty}\big|V^{N}(\muin, \mufin) -  V(\muin, \mufin)\big| \le 2\varepsilon.
\end{equation*}

Dropping $\varepsilon$ yields the result.

§.§ Finite particle approximation under regularity of the data

Our second result holds for regular data. It allows us to work with an arbitrary terminal distribution $\mufin\in \Pc_2(\R^\xdim)$ in exchange for some boundedness assumptions on the solution to the limiting system. These can ultimately be verified under regularity assumptions on the data, see <Ref> for details.

In addition to <Ref>, $\sigma\in \R^{\xdim\times d}$ is nonsingular and the mappings $F$ and $B$ defined by Equations (<ref>) and (<ref>) satisfy that for $(t,x,y,\mu,\nu) \in [0,T]\times \RR^\xdim \times \RR^\xdim \times \cP_2(\RR^\xdim) \times \cP_2(\RR^\xdim\times \RR^\xdim)$

\begin{align}\label{eq:growth.ass.b.f.convergence}
\| F(t,x,y,\nu)\| \leq \ell_F \bigg(1+\|y\|+\int_{\R^\xdim} \|y\| (\pi^2_\# \nu)(\d y)  \bigg), \; B(t,x,y,\mu) \leq \ell_B \big(1+\|y\|  \big),
\end{align}

as well as the following regularity conditions

* $B$ $($resp. $F)$ is continuously differentiable in $(t,w, \mu)$, for $w=(x,y)$ $(resp.\, (t,w, \nu))$, with bounded and Lipschitz-continuous derivatives.

Moreover, for any $X,Y\in \L^2(\Omega)$, the families $(\partial_\mu B(t,x,y,\nu) (X))_{(t,x,y)}$ and $(\partial_\mu F(t,x,y,\nu) (X,Y))_{(t,x,y)}$ have uniformly integrable second moments$;$

* for any $\mu$ $($resp. $\nu)$ the mapping $w\longmapsto \partial_w B(w,\mu)$ $($resp. $w\longmapsto \partial_w F(w,\nu))$ is continuously differentiable with bounded and Lipschitz-continuous derivatives. In addition, $(w,x)\longmapsto \partial_\mu B(w,\mu)(x)$ and $(w,w^\prime)\longmapsto \partial_\mu F(w,\mu)(w^\prime)$ are continuously differentiable with bounded and Lipschitz-continuous derivatives.

Let <Ref> hold with the strengthening that $\sigma\in \R^{\xdim\times d}$ is nonsingular, and the mappings $F$ and $B$ defined by (<ref>) and (<ref>) satisfy that for $(t,x,y,\mu,\nu) \in [0,T]\times \RR^\xdim \times \RR^\xdim \times \cP_2(\RR^\xdim) \times \cP_2(\RR^\xdim\times \RR^\xdim)$

\begin{align}\label{eq:growth.ass.b.f.convergence}
\| F(t,x,y,\nu)\| \leq \ell_F \bigg(1+\|y\|+\int_{\R^\xdim} \|y\| (\pi^2_\# \nu)(\d y)  \bigg), \; B(t,x,y,\mu) \leq \ell_B \big(1+\|y\|  \big).
\end{align}

Suppose $g$ is a smooth penalty function for which $(x,\mu)\longmapsto \partial_\mu g(\mu)(x)$ is bounded and Lipschitz continuous with Lipschitz constant $\ell_g>0$. If there exists $C_{\tilde z}>0$ such that $Z^k$ solution to (<ref>) satisifes

\begin{align}\label{eq.assumptionzbounded}
\|Z_t^k  \|_\infty \leq k C_{\tilde z}, \; \text{for almost every } \; t\in [0,T].
\end{align}

Then, we have that

\begin{equation*}
\lim_{N\to \infty}V^N(\muin,\mufin) = V (\muin,\mufin).
\end{equation*}

The above assumptions can be explained by our strategy to obtain the convergence of the particle systems in this case.
We will translate the convergence of the FBSDE (<ref>) into that of a quadratic BSDE.
As such, we need to have control of the norm of $($part of$)$ the solution $(Y,Z)$, and to perform appropriate changes of measure.
The former we obtain from the growth assumptions on $F$ in (<ref>) and condition (<ref>), which makes explicit the dependence on $k\geq 1$ and on which elaborate more in <Ref>.
The latter we obtain from the assumptions on $B$ in (<ref>) and on $\sigma$.
Ultimately, back to the entropic problem, the above assumption reduces to some extra regularity on the interaction potential together with Condition <ref>, see <Ref> and <Ref>.Lastly, we note that the above result imposes extra assumptions on $\partial_\mu g(\mu)(x)$.
Since $\partial_\mu g(\mu)(x)$ shows up in the analysis due to our penalization approach and despite our list of examples of smooth penalty functions in <Ref>, it is not immediate that there is a smooth penalty function satisfying those extra assumptions.
Fortunately, the Gauge function introduced in [Cosso et al., 2022] can be modified to provide a penalty function verifying them, although it will no longer be a Gauge function per se.

There exists a smooth penalty function $g$ whose derivative $\partial_\mu g(\mu)(x)$ satisfies:

* the map $x \longmapsto\partial_\mu g(\mu)(x)$ is continuous differentiable, with bounded and Lipschitz continuous derivative $\partial_x \partial_\mu g(\mu)(x)$,

* the map $\mu\longmapsto \partial_\mu g(\mu)(x)$ is $L$-differentiable, with bounded and Lipschitz-continuous derivative $\partial_{\mu\mu}^2 g(\mu)(x)(y)$,

* the map $(x,y) \longmapsto\partial_{\mu\mu}^2 g(\mu)(x)(y)$ is continuously differentiable with bounded and Lipschitz-continuous derivatives,

* there is $\ell_g>0$ depending only on $m$, $\delta$ and $\varrho$ such that for all $(\mu,x,y)\in \Pc_2(\R^m)\times \R^m\times \R^m$

\begin{align}\label{eq:boundsgauge}
|\partial_\mu g(\mu)(x)|+|\partial_x \partial_\mu g(\mu)(x)|+|\partial_{\mu\mu}^2 g(\mu)(x)(y)|+|\partial_x \partial_{\mu\mu}^2 g(\mu)(x)(y)|+|\partial_y \partial_{\mu\mu}^2 g(\mu)(x)(y)|\leq \ell_g .
\end{align}

Let the assumptions in <Ref> hold.

Let $(  X^{i,N},  Y^{i,N})$ be part of the solution to (<ref>), and let $( \widetilde X^i, \widetilde Y^i)$, $i=1,\dots,N$, denote i.i.d. copies of $( X, Y, Z)$ solution to (<ref>). Then $( X^{i,N},Y^{i,N})$ converges to $(  \widetilde X^i,\widetilde Y^i)$ in $\Hc^2$.

Let us comment on the small time horizon assumption. The reader familiar with the literature of FBSDEs might have considered the direct analysis of systems (<ref>) and (<ref>) over a small time horizon without any regularity assumptions. Yet, the crux of the problem with this approach is that the resulting upper bound on $T$ turns out to be of order $1/k$. Thus, even though the result holds for $k$, i.e., $T^k$, fixed, we cannot gain much out of it when we bring ourselves back to the proof of <Ref> where we will have to let $k$ go to infinity. The regularity assumptions in <Ref> allow us to bypass this issue.

On the one hand, we remark that the argument we present below can not be iterated, without extra assumptions, to arbitrary length intervals as it is usual in the analysis of FBSDEs. This is because the solution to the master equation given by <cit.> fails to have derivatives which are globally Lipschitz with respect to the measure argument, see <cit.>. Thus, to extend the analysis to arbitrary time horizons one has to assume, see <cit.>, the well-posedness of the associated FBSDE, i.e., of (<ref>), on the entire time interval, and show an apriori estimate on $Y$ which is: $(i)$ independent of the terminal condition, i.e., on the penalty function and hence on the target distribution$;$ $(ii)$ uniform in the length of the interval. These rather strong assumptions only hold on a case-by-case basis, see <cit.>.

We argue in four steps.

Step 1: In this step we introduce two auxiliary BSDEs. Let $( X, Y, Z)$ be a solution to the McKean-Vlasov FBSDE (<ref>). Notice that, by assumption, we have that for $t\in [0,T]$

\|Y_t\| \leq \|Y_T\| +\int_t^T \ell_F (1 +   \| Y_r\|+\E[\| Y_r\|] ) \d r -\int_t^T Z_r \d W_r.

It then follows from taking expectation and applying Gronwall's lemma that
$\E[ \|Y_t\| ] \leq ( k\ell_g +T\ell_F)\e^{ 2 \ell_F T}$.
Let us now note that under assumption <Ref> there exists a process $\alpha$ bounded by $\ell_F$ such that

\[
Y_t=Y_T+\int_t^T (\alpha_r Y_r + F(r,X_r,0,\Lc(X_r,Y_r))\big)  \d r-\int_t^T Z_r \d W_r.
\]

where the product is coordinate-wise. It then follows from Itô's formula,

(<ref>), and taking conditional expectation that

\[
\|Y_t\|\leq \e^{\ell_F T} \|Y_T\| +\e^{\ell_F T} \int_t^T \E[\| Y_r\|] \d r \leq \e^{\ell_F T}( k\ell_g+ T ( k\ell_g +T\ell_F)\e^{ 2 \ell_F T}).\]

This is, together with (<ref>) we obtain that there exists $C_{\tilde y},C_{\tilde z}>0$ such that

\[
\|Y_t^k\|_\infty \leq k C_{\tilde y},\; t\in [0,T], \text{ and, } \|Z_t^k  \|_\infty \leq k C_{\tilde z}, \; \text{for almost every } \; t\in [0,T].
\]

Let $( \widetilde X^i, \widetilde Y^i, \widetilde Z^{i})$, $i=1,\dots,N$, denote i.i.d. copies of $( X, Y, Z)$ solving the McKean-Vlasov FBSDE (<ref>) and consider the probability

\[
\frac{\d \P^B}{\d \P}:=\Ec \bigg(- \sum_{i=1}^N \int_0^T \sigma^{-1} B(t, \widetilde  X_t^i, \widetilde  Y_t^i,\Lc( X_t))\d W_t^i\bigg),
\]

where $\Ec(\cdot)$ denotes the Doleans-Dade stochastic exponential. Notice that $\P^B$ is well-defined in light of the growth assumption on $B$ and the boundedness of the $\widetilde Y^i$'s. By Girsanov's theorem we have that for $i=1,\dots,N$, $W^{i,B}:=W^i+\int_0^\cdot \sigma^{-1} B(t,  \widetilde  X_t^i,  \widetilde  Y_t^i, \cL(  X_t))\d t$ is a $\P^{B}$–Brownian motion and $ \widetilde  X^i_t = \widetilde X_0^i+ \sigma W_t^{i,B}$. Moreover, $(\widetilde Y^i,\widetilde Z^i)$ solves the McKean–Vlasov BSDE under $ \P^{B}$ given by

\begin{align}\label{eq:convergence:BSDEgirsanov}
\begin{cases}
\d  \widetilde Y_t^i =- \big( F(t,  \widetilde X_t^i, \widetilde  Y_t^i, \cL( X_t, Y_t)) + \sigma^{-1} B(t, \widetilde X_t^i, \widetilde Y_t^i, \cL(  X_t)) \widetilde Z_t^i  \big)  \diff t + \widetilde  Z_r^i\diff W_r^{i,B},\\
\widetilde Y_T^i = k\partial_\mu	g(\cL( X_T))( \widetilde X_T^i) .
\end{cases}
\end{align}

The second auxiliary BSDEs is given by $(\hat Y^N,\hat Z^N)$ solution under $\P^{B}$ to the multidimensional coupled BSDE

\begin{align}\label{eq:convergence:Npart:BSDEgirsanov}
\begin{cases}
\d  \hat  Y_t^{i,N} =-\big( F(t, \widetilde  X_t^i, \hat Y_t^{i,N},L^N(\widetilde \X^{N}_t,\hat \Y^{N}_t) ) + \sigma^{-1} B(t, \widetilde X_t^i, \hat Y_t^{i,N}, L^N(\widetilde \X^{N}_t) )\hat Z_t^{i,i,N}  \big)  \diff t + \sum_{j=1}^N  \hat Z_t^{i,j,N}\diff W_t^{j,B},\\
\hat Y_T^{i,N}  = k\partial_\mu	g(\cL( X_T))( \widetilde X_T^i) .
\end{cases}
\end{align}

The well-posedness of a global solution to (<ref>) follows from <cit.> provided we verify (H1)-(H4) therein.
Noticing the generator of (<ref>) is diagonally quadratic, we see that (H1) and (H2) hold.
To verify (H3) recall that the terminal condition is bounded by assumption, whereas the boundedness of the process $\alpha$ appearing in [Fan et al., 2020] follows from the growth assumptions of $F$ and $B$.
Lastly, (H4) follows from the diagonal quadratic structure of the generator.
We also notice that, analogous to the estimate for $Y$ at the beginning of this step, there is $C_{\hat y}>0$ independent of $N$ such that for any $i=1,\dots,N$, $\| \hat  Y^{i,N}_t\|\leq  k C_{\hat y}, \; t\in [0,T]$.

Step 2: We derive an estimate for the BSDEs (<ref>) and (<ref>).
For this let us set $\delta Y^i: = \widetilde Y^i -  \hat Y^{i,N}$ and $\delta Z^{i,j}: = \widetilde Z^{i}1_{\{i=j\}} - \hat Z^{i,j,N}$, so that $\delta Y^i_T=0$ and, by assumption, we have that $\d \P^B\otimes \d t\text{--a.e.}$

\begin{align*}
\|F(r,  \widetilde X_r^i, \widetilde Y_r^i, \cL( X_r, Y_r))- F(r, \widetilde  X_r^i, \hat Y_r^{i,N},L^N( \widetilde \X^{N}_t,\hat \Y^{N}_t) )\| & \leq \ell_F \big( \|\delta Y^i_t\| + \Wc_1(L^N(\widetilde \X_t^N,\hat \Y_t^N), \Lc(X_t, Y_t)) \big),\\
\| B(r, \widetilde X_r^i,\widetilde  Y_r^i, \cL(  X_r))   Z_r^i   - B(r,\widetilde  X_r^i, \hat Y_r^{i,N}, L^N(\widetilde  \X^{N}_t) )\hat Z_r^{i,i,N} \|  & \leq k C_{\tilde z}   \big(   \|\delta Y^i_t\| + \Wc_1(\Lc(X_t),L^N(\widetilde \X_t^N)) \big) +  \ell_B(1+ k C_{\hat y})  \|\delta Z^{i,i}_t\|.
\end{align*}

Moreover, we have that

\begin{align}
\| \delta Y^i_t\|^2 &\leq\int_t^T 2\| \delta Y^i_r\| \Big( \big(\ell_F+\frac{k C_{\tilde z}}\sigma\big)\big( \| \delta Y^i_r\|   + \Wc_1(L^N(\widetilde \X_t^N,\hat  \Y_t^N), \Lc(X_t, Y_t)) \big) + \frac{\ell_B}\sigma (1+k C_{\hat y}) \| \delta Z^{i,i}_r \|  \Big) \d r\notag -\int_t^T \|\delta Z_r^{i,i}\|^2\d r+ M_t^i\notag \\
%& \leq \int_t^T\big( C^k \| \delta Y^i_r\|^2 +\Wc_1^2(\Lc(  X_r),L^N(\widetilde  \X_r^N)) +  \Wc_1^2(\Lc(  Y_r),L^N(  \hat \Y_r^N)) \big) \d r+M_t^i\notag \\
& \leq \int_t^T\big( C^k \| \delta Y^i_r\|^2+  \frac{1}N \sum_{j=1}^N \|\delta Y_r^j\|^2   +   \Wc_1^2(L^N(\widetilde \X_t^N,\widetilde  \Y_t^N), \Lc(X_t, Y_t)) \big)  \d t+ M_t^i  \label{eq:aux:convergence1}
\end{align}

where, $M_t^i:= \sum_{j=1}^N \int_t^T 2\| \delta Y^i_r\| \|\delta Z^{i,j}_r\| \d W_r^{j,B}$, in the second inequality, we used Young's inequality and introduced the constant $C^k:=2(\ell_F+kC_{\tilde z}\sigma^{-1})+ (\ell_F+ kC_{\tilde z}\sigma^{-1})^2+ \ell_F^2 +  \ell_B^2 \sigma^{-2} (1+kC_{\hat y})^2$.
Averaging over $i$ we find that

\begin{align*}
\frac{1}N \sum_{i=1}^N \|\delta Y_t^i\|^2\leq \int_t^T\big(\Wc_1^2(L^N(\widetilde \X_t^N,\widetilde  \Y_t^N), \Lc(X_t, Y_t)) + (1+C^k) \frac{1}N \sum_{i=1}^N \|\delta Y_r^i\|^2 \big) \d t+ \frac{1}N  \sum_{i=1}^N M_t^i ,
\end{align*}

which in turn, taking expectation and using Gronwall's inequality and with $\epsilon^N$ as in (<ref>), leads to

\begin{align}\label{eq:aux:convergence2}
\E^{\P^B}\bigg[ \frac{1}N \sum_{i=1}^N \|\delta Y_t^i\|^2\bigg] \leq  \e^{ (1+C^k)T} \E^{\P^B}\big[\epsilon^N\big] .
\end{align}

We now come back to (<ref>). Taking expectation, plugging (<ref>), and using Gronwall's inequality we find that

\begin{align}\label{eq.estimate.quadraticBSDEs}
\E^{\P^B} \bigg[\int_0^T \| \delta Y^i_t\|^2\d t \bigg] \leq T \e^{ C^k T} (\e^{ (1+C^k)T}+1) \E^{\P^B}\big[\epsilon^N\big] .
\end{align}

Step 3: We now show that there exists a common probability space $(\mybar \Omega, \mybar \Fc,\mybar \P)$, supporting $(\mybar X^{i,N},\mybar Y^{i,N},\mybar W^{i,N})\stackrel{d}{=} ( X^{i,N}, Y^{i,N},W^i)$, and $(\mybar X^i,\mybar Y^i, \mybar W^i)\stackrel{d}{=} ( \widetilde X^i,  \widetilde Y^i, W^i)$, on which $(\mybar X^{i,N},\mybar Y^{i,N})$ converges to $(\mybar X^i,\mybar Y^i)$, $\mybar \P\text{--a.s.}$

For this we consider (<ref>) as well as the measure

\[
\frac{\d \P^{B^N}}{\d \P}:=\Ec \bigg(- \sum_{i=1}^N \int_0^T \sigma^{-1} B(t, X_t^{i,N}, Y_t^{i,N},L^N(\X_t^N))\d W_t^{i,B}\bigg).
\]

Since, as we argued for $\hat Y^{i,N}$ in Step 1, $Y^{i,N}$ is bounded, the growth condition on $B$ guarantees, via Girsanov's theorem, that for $i=1,\dots,N$, $W^{i,B^N}:=W^i+\int_0^\cdot \sigma^{-1} B(t, X_t^{i,N}, Y_t^{i,N},L^N(\X^{N}_t))\d t$ is a $\P^{B^N}$–Brownian motion, $X^{i,N}_t =  X_0^i+ \sigma W_t^{i,B^N}$ and $ \P^{B^N}\text{--a.s.}$

\begin{equation}\label{eq:convergence:OriginalNpart:BSDEgirsanov}
\begin{cases}
\d  Y_t^{i,N} =-\big( F(r,  X_r^{i,N},  Y_r^{i,N},L^N( \X^{N}_t, \Y^{N}_t) ) + \sigma^{-1} B(r, X_r^{i,N},  Y_r^{i,N}, L^N(\X^{N}_t) ) Z_r^{i,i,N}  \big)  \diff t + \sum_{j=1}^N   Z_r^{i,j,N}\diff W_r^{j,B^N},\\
Y_T^{i,N}  = k\partial_\mu	g(\cL( X_T))( X_T^{i,N}).
\end{cases}
\end{equation}

We claim that $(X^{i,N},Y^{i,N})$ converges in law to $(\widetilde X^{i},\widetilde Y^{i})$. Indeed, let us first note that $\P^{B^N}\circ (X^{i,N})^{-1} = \P^{B}\circ (\widetilde X^{i})^{-1}$. Second, we verify that $\P^{B}\circ (\widetilde Y^{i})^{-1}$ does not depend on $N$. To see this, let

\[
\frac{\d \P^{i}}{\d \P}:=\Ec (Z^i):=\Ec \bigg( -  \int_0^T \sigma^{-1} B(t, \widetilde X_t^i,\widetilde Y_t^{i},\Lc( X_t))\d  W_t^{i}\bigg),
\]

and note that for any continuous and bounded function $f:\Cc([0,T];\R^\xdim)\longrightarrow \R$, the independence of $(\widetilde X^i,\widetilde Y^i,W^i)$ implies

\[
\E^{\P^B}[f(\widetilde Y^i)]=\E^{\P}\bigg[f(\widetilde Y^i)\prod_{k=1}^N \Ec(Z^k)\bigg]=\E^{\P}\bigg[f(\widetilde Y^i)\Ec(Z^i) \prod_{k\neq i}^N \Ec(Z^k)\bigg]=\E^{\P}\big[f(\widetilde Y^i)\Ec(Z^i)\big]=\E^{\P^i}[f(\widetilde Y^i)].
\]

Lastly, noticing that (<ref>) and (<ref>) are both weak solutions of the same BSDE, and by <cit.> this BSDE has a unique weak solution, it follows from (<ref>) and the law of large numbers that

\[
\P^{B^N}\circ (X^{i,N},Y^{i,N})^{-1} \xrightarrow{N\rightarrow \infty} \P^{B}\circ (\widetilde X^{i},\widetilde Y^{i})^{-1}.
\]

The conclusion of Step 3 follows from Skorodhod representation theorem.

Step 4: We now consider (<ref>) and (<ref>) under $\P$.

With a slight abuse of notation we let $\delta X^i: =  \widetilde X^i -  X^{i,N}$, $\delta Y^i: =  \widetilde Y^i -  Y^{i,N}$, and $\delta Z^{i,j}: =  \widetilde Z^{i}1_{\{i=j\}} - Z^{i,j,N}$, so that for $N_t^i:= \sum_{j=1}^N \int_t^T 2\| \delta Y^i_r\|\|\delta Z^{i,j}_r\| \d W_r^{j}$, we have that

\begin{align*}
\| \delta Y^i_t\|^2 \leq&\, 2 k\ell_g ( \|\delta X_T^i\|^2  +\Wc_1^2(\Lc( \widetilde X_T^i),\Lc(X_T^{i,N}))  )
-\int_t^T \|\delta Z_r^{i,i}\|^2\d r+N_t^i\\
&  + \ell_F \int_t^T 2 \| \delta Y^i_r\|  \big( \| \delta X^i_r\| + \| \delta Y^i_r\|   +\Wc_1(L^N( \X_r^N,\Y_r^N),\Lc( X_r,Y_r))  \big) \d r\\
\leq &\, 2 k\ell_g \big( \|\delta X_T^i\|^2  +\Wc_1^2(\Lc( \widetilde X_T^i),\Lc(X_T^{i,N}))  \big)
+ N_t^i \\
&  + \ell_F \int_t^T \big( 6 \| \delta Y^i_r\|^2+ \| \delta X^i_r\|^2 +  \frac{1}N \sum_{i=1}^N \|\delta Y_r^i\|^2  +  \frac{1}N \sum_{i=1}^N \|\delta X_r^i\|^2  +\Wc_1^2(L^N( \widetilde \X_r^N,\widetilde \Y_r^N),\Lc( X_r, Y_r))  \big) \d r.
\end{align*}

where we used the Lipschitz continuity of $(x,\mu)\longmapsto \partial_\mu g(\mu)(x)$. As in Step 3, taking expectation, averaging and using Gronwal's inequality, we find that there is a constant $C>0$, depending only $T$, $\ell_g$ and $\ell_F$, such that

\begin{align*}
\E\bigg[ \int_0^T \| \delta Y^i_r\|^2 \d r\bigg] &\leq C  \bigg(  k \E\big[ \|\delta X_T^i\|^2\big] +k\Wc_1^2(\Lc( \widetilde X_T^i),\Lc(X_T^{i,N}))  +  \int_0^T\Big( \E\big[  \| \delta X^i_r\|^2 \big] +  \frac{1}N \sum_{i=1}^N \E[ \|\delta X_r^i\|^2]  \Big) \d r+\E\big[\epsilon^N\big]\bigg)
\end{align*}

We now claim the expression on the right side goes to zero as $N$ tends to infinity.
The last term goes to zero by the law of large numbers.
We now recall that the laws of $X^i$ and of $X^{i,N}$ are uniquely determined by the controlled SDE satisfied by $X^i$ and $X^{i,N}$, respectively.
Thus, since $\mybar X^{i}$ and $\mybar X^{i,N}$ satisfied the corresponding controlled SDEs under $\mybar \P$, we have that $\Wc_1^2(\Lc( \tilde X_T^i),\Lc(X_T^{i,N}))=\Wc_1^2(\Lc_{\bar \P}(\bar X_T^i),\Lc_{\bar \P}(\mybar X_T^{i,N}))$.
Likewise, the expectation under $\P$ and $\mybar \P$ of $ \|\delta X_t^i\|^2$ and $\frac{1}N \sum_{i=1}^N  \|\delta X_t^i\|^2$ agree. The claim then follows from Step 3.At last, we note that the convergence of $\E\displaystyle \bigg[ \int_0^T \| \delta Y^i_r\|^2 \d r\bigg]$ implies that of $\E\big[\sup_{t\in [0,T]} \| \delta X^i_r\|^2 \big]$. This concludes the proof.

Coming back to <Ref>, we remark that the argument follows identically as in the proof of <Ref> noticing that the conclusion at the end of Step 1 now follows thanks to <Ref>. We close this subsection with some comments on the unviability of other strategies and the assumption (<ref>).

Let us now further discuss the unviability of two approaches based on other known scenarii in which the study of the convergence problem in mean field games (and uniqueness of mean field games) is often achieved. Indeed, in the literature, either small time conditions or monotonicity/convexity conditions are always needed to achieve strong convergence results [Cardaliaguet et al., 2019, Jackson and Tangpi, 2023, Laurièrie and Tangpi, 2022].
In our setting, for $k$ fixed, the associated FBSDE is also well-posed over a small time horizon, said $T_o^k$. This suggests two frameworks for the study of the convergence result using the same proof idea as in <Ref>.

* A first approach would be based on the results of  [Chassagneux et al., 2014] which provide conditions for the well-posedness of bounded solutions to the so-called master equation in two scenarii. The first one is a small time horizon $T_o^k$, whereas the second scenario iterates the previous one to be able to obtain a result over arbitrary time horizons. Though this second result seems viable at first sight, the conditions under which the iteration can be carried up impose assumptions on the type of interaction between the particles, which are more restrictive than those in <Ref>, it also requires extra regularity assumptions on the data and, perhaps more critically, it imposes the same convexity assumptions on the penalty function prescribed by <Ref>.

* Alternatively, the theory of Malliavin calculus allows us to identify another scenario under which one can show the existence of $T_o^k <\infty$, such that for $T\leq T_o^k$, there is a bounded solution $(Y^k,Z^{k})$ to (<ref>).

In either setting, the particle system (<ref>) converges to (<ref>) over $T_0^k$, for every $k\geq 1$. However, to establish $V^N(\muin,\mufin)\longrightarrow V (\muin,\mufin)$ we see from Step 3 in the proof of Theorem <ref> that there needs to exist $T_o$ under which the convergence of the associated FBSDEs holds uniformly in $k$. Unfortunately, in either $(i)$ or $(ii)$ above $T^k_o$ decreased to $0$ as we let $k$ go to infity.

Building upon the previous remark, we now discuss condition in (<ref>) in more detail. First, it should be clear that (<ref>) allows us to bypass the issues raised in <Ref> since it is an estimate on $[0,T]$. Second, let us describe a simple setting, familiar to those acquainted with the Malliavin analysis of FBSDEs, under which (<ref>) can be verified. It develops the scenario in <Ref><ref> which consists on assuming (<ref>), that $(x,\mu)\longmapsto \partial_\mu g(\mu)(x)$ is bounded and Lipschitz continuous, that $(x,y)\longmapsto (F,B)(t,x,y,\mu)$ is continuously differentiable, and that $X_0\sim\muin$ is Malliavin differentiable with bounded Malliavin derivative. If, in addition, the condition: for any $k\geq 1$ there is

\begin{align}\label{eq.decoupling} \tag{D}
U^k:[0,T]\times\R^\xdim\times \Pc_2(\R^\xdim)\longmapsto \R^\xdim,\text{ bounded, s.t. } \P[Y_t^k=U^k(t,X_t^k,\Lc(X_t^k))]=1, \; t\in [0,T],
\end{align}

for which $x\longmapsto U^k(t,x,\mu)$ is differentiable with bounded derivative, holds, then we can verify (<ref>). We highlight that $U^k$ is to the so-called decoupling field for the FBSDE (<ref>).

§.§ The entropic case: Proof of Theorem <ref>

The result follows from specializing <Ref> to the entropic case in which $f$ and $b$ are as in <Ref>.
In this case, $\Lambda(t,y)=-y$ and the Hamiltonian (<ref>) and the functionals $F$ and $B$ takes the form

\begin{align*}
H^{N}(t,\bx, {\bf a}, \by,,L^N({\bf x})) &= \sum_{i=1}^N \bigg[  \frac 12 \|a^i\|^2
+ \bigg( -\frac1N\sum_{j=1}^N\nabla \Psi( x^i - x^j) +  a^i  \bigg) \cdot y^{i} \bigg],\\
F(t,x,y,L^N({\bf x},{\bf y}))&=- \frac1N\sum_{j=1}^N\nabla^2 \Psi(x - x^{j}) \cdot y + \frac1N\sum_{j=1}^N\nabla^2\Psi(x^{j} - x)\cdot y^{j},\; B(t,x,y,L^N({\bf x}))=y-\frac1N\sum_{j=1}^N\nabla \Psi( x - x^j).
\end{align*}

Let us argue the result. We will use <Ref> for which we must verify <Ref> and the strong convexity of $f_1(t,\cdot)$. The displacement convex penalty function is given by assumption. As in the proof of <Ref>, <Ref> holds in light of the assumptions on $\Psi$. Noticing that $f_1(t,a)=\|a\|^2/2$, we see that the strong convexity assumption hold. In fact, since $\Lambda(t,y)=-y$ we see that (<ref>) holds with equality and $\lambda_{f_1}=1$. The result follows from <Ref>.

Let us comment on the implications of <Ref> in the entropic case.
Indeed, suppose that $\Psi$ is three times differentiable with bounded and Lipschitz–continuous derivatives, that $X_0\sim\muin$ is Malliavin differentiable with bounded derivatives and that Condition <ref> holds.
We claim that the finite particle approximation holds.
Recall that the condition on the smooth penalty function $g$ is satisfied by the modification of the gauge function alluded to in <Ref>.
From the form of $F$ and $B$ above, it's clear that the growth assumptions on $F$ and $B$ are satisfied thanks to the boundedness of the derivatives of $\Psi$. Likewise, the differentiability of $(x,y)\longmapsto (F,B)(t,x,y,\mu)$ follows from that of $\Psi$. Since by assumption we have that Condition <ref> holds, as argued in <Ref>, we have that (<ref>) is satisfied. Thus, <Ref> holds and establishes the claim.

[inline]Camilo: the initial statement included... “if for each $N$ the $N$-particle Schrödinger problems admit a Schrödinger bridge $(\hat X^{i,\np})_{i\in \{1,\dots,\np\}}$, then for each $i$ the sequence $(X^{i,\np})_{\np\ge 1}$ converges in law to a mean field Schrödinger bridge $\hat X$”... I don't think we can show that (again the issue is with $k$). Writing

\[
\hat X^{N}- \hat X= \hat X^{N}-\hat X^{k,N}+\hat X^{k,N} - \hat X^{k}+\hat X^{k}-\hat X,
\]

we see that once again, we will need an estimate on $\hat X^{k,N} - \hat X^{k}$ which is uniform in $k$.

§ SMALL NOISE AND LONG TIME ASYMPTOTIC PROPERTIES FOR THE MFGP

§.§ Small noise asymptotics

Let us consider the problem

\begin{equation}
\label{eq:problem.eps}
V_\eps(\muin, \mufin) = \inf\bigg\{\EE\bigg[\int_0^Tf(t, X^{\alpha,\eps}_t, \alpha_t,\mathcal{L}(X^{\alpha,\eps}_t))\diff t \bigg]: \alpha \in \mathfrak{A}, \; X^{\alpha,\eps}_0\sim  \muin \text{, and } X^{\alpha,\eps}_T\sim  \mufin\bigg\},
\end{equation}

where $X^{\alpha,\eps}$ is the unique strong solution of the controlled stochastic differential equation

\begin{equation}
\label{eq:control.SDE.MF.eps}
X^{\alpha,\eps}_t = X_0^{\alpha,\eps} + \int_0^t\alpha_u+b(u, X^{\alpha,\eps}_u, \mathcal{L}(X^{\alpha,\eps}_u))\diff u + \eps\sigma W_t.
\end{equation}

The goal of this part is to derive the limiting problem when $\eps\downarrow 0$.

For this, we introduce the following set of assumptions

Let <Ref> be satisfied. We further assume that $f_2(t,\cdot,\cdot)$ is bounded and that there is $C>0$ such that $\langle y,\sigma^\top \sigma y\rangle> C\|y\|$, for any $y\in \R^\bmdim$. Then, $(V_\eps(\muin,\mufin))_{\eps>0}$ converge to $V_0(\muin,\mufin)$ satisfying

\begin{equation*}
V_0(\muin, \mufin) = \inf_{(\phi,\mu)}\, \int_0^T\int_{\RR^\xdim}f\big(t,x, \phi(t,x,\mu_t),\mu_t\big)\mu_t(\d x)\diff t,
\end{equation*}

where the infimum is over Borel–measurable maps $\phi:[0,T]\times \R^m\times\cP_2(\R^\xdim)\longrightarrow \R^\xdim$ and absolutely continuous curves $\mu_t:[0,T]\longrightarrow \Pc_2(\R^m)$ such that $\phi(t,\cdot,\mu_t)\in \mathbb{L}^p(\RR^d,\mu_t)$ and with $\mu$ satisfying the continuity equation

\begin{equation}
\label{eq:contin.eqt}
\partial_t\mu  + \nabla\cdot \Big( \mu  b(t, \cdot,\mu ) + \mu \sigma\phi(t, \cdot,\mu )\Big) = 0, \quad \mu_0 = \muin, \quad \mu_T = \mufin,
\end{equation}

in the sense of distributions.

If the function $b$ does not depend on $\mu$, i.e. $b(t, x, \mu) = b(t,x)$, then it holds

\begin{equation*}
V_0(\muin,\mufin) := \inf_{(\phi,\mu)}\,  \int_0^T\int_{\R^\xdim} f\big(t,x, \phi(t, x), \mu_t\big)\mu_t(\d  x)\diff t ,
\end{equation*}

where the supremum is over Borel–measurable maps $\phi:[0,1]\times \R^\xdim\longrightarrow \R^\xdim$ and absolutely continuous curves $\mu_t:[0,T]\longrightarrow \Pc_2(\R^m)$ such that $\phi(t,\cdot) \in \L^p(\R^\bmdim,\mu_t)$ and

\begin{equation*}
\partial_t\mu + \nabla\cdot \big( \mu  b(t, \cdot) + \mu  \sigma\phi(t, \cdot) \big) = 0, \quad \mu_0 = \muin,\quad \mu_T = \mufin,
\end{equation*}

in the sense of distributions.

The zero–noise control problem with value $V_0(\muin,\mufin)$ obtained above should be reminiscent to the Benamou–Brenier formulation of (static) optimal transport in fact, taking $b=0$, $\sigma=1=T$ and $f(t, x, a, \mu) = \|a\|^2$ above $V_0(\muin, \mufin)$ becomes

\begin{align*}
V_0(\muin,\mufin) &= \inf\bigg\{ \int_0^1\int_{\R^\xdim} \| \phi(t, x)\|^2\mu^\phi_t(\diff x)\diff t: 	\quad	\partial_t\mu^\phi + \nabla\cdot \big( \mu^\phi \phi(t, \cdot) \big) = 0, \quad \mu^\phi_0 = \muin,\quad \mu^\phi_T = \mufin \bigg\} \\
&= \inf\bigg\{\iint_{\R^\xdim\times\R^\xdim}\|x-y\|^2\pi(\diff x, \diff y), \quad \pi \in \Gamma(\muin, \mufin) \bigg\}
\end{align*}

where the second inequality follows by  <cit.>.

The proof of this result starts again with the penalization argument used above.
In fact, let us consider the McKean–Vlasov control problem

\begin{equation}
\label{eq:problem.eps.penal}
V_\eps^{k}(\muin, \mufin) := \inf_{\alpha \in \mathfrak{A}}\EE\bigg[\int_0^Tf(t, X^{\alpha,\eps}_t, \alpha_t,\mathcal{L}(X^{\alpha,\eps}_t))\diff t + kg(\cL(X^{\alpha,\eps}_T)) \bigg]
\end{equation}

and its zero–noise analogue

\begin{equation*}
V_0^k(\muin,\mufin) := \inf_{(\phi,\mu)} \bigg(\int_0^T\int_{\R^\xdim} f\big(t,x, \phi(t,x,\mu_t),\mu_t\big)\mu_t(\d  x)\diff t +kg(\mu_T)\bigg)
\end{equation*}

with the infimum taken over  $\phi$ and $\mu$ satisfying

\begin{equation}
\label{eq:contin.eqt.k}
\partial_t\mu_t + \nabla\cdot \Big( \mu_t b(t, \cdot,\mu_t) + \mu_t\sigma\phi(t, \cdot,\mu_t)\Big) = 0, \quad \mu_0 = \muin %, \quad \mu^\phi_T = \mufin
\end{equation}

in the sense of distributions.

We begin noticing that by <Ref> we know that

\begin{equation*}
\lim_{k\to \infty} V^k_\eps(\muin,\mufin) =V_\eps(\muin,\mufin) .