Proof. The second-variation formula along constant-speed geodesics in Wasserstein space is well known in the literature. It appears formally in [34], and under our stated assumptions we provide a rigorous proof in Lemma 6.6. The formula reads \[\label{eq-165} \frac{d^2}{dt^2}\Big \vert_{t=0}F((id+t \nabla \phi)_\# \mu^* , \nu^*) = \int_{\mathcal{X}} \bigl (|\nabla^2 \phi |^2 + \langle \nabla \phi, \nabla^2 V_{\nu^*} \nabla \phi\rangle \bigr) \, d\mu^* .\tag{13}\] Here \(V_{\nu^*}(x)=\int_{\mathcal{Y}} f(x,y)\,d\nu^*(y)\), and \(|\nabla^2\phi|\) denotes the Hilbert–Schmidt norm. At the equilibrium \(\frac{d\mu^*}{dx} \propto e^{-V_{\nu^*}}\) we are considering, it is non-negative (Lemma 12 or [35]) \[\label{eq-166}\int \bigl (|\nabla^2 \phi |^2 + \langle \nabla \phi, \nabla^2 V_{\nu^*} \nabla \phi\rangle \bigr) \, d\mu^ *= \int |\Delta \phi - \langle \nabla \phi ,\nabla V_{\nu^*}\rangle |^2 \, d\mu^* \ge 0.\tag{14}\] This shows that the Hessian is nonnegative in every gradient direction \(\nabla\phi\).

Moreover, the quadratic form in 14 does not vanish on any nontrivial direction. Indeed, if it were zero for some \(\phi\), then \(\Delta \phi - \langle \nabla \phi,\nabla V_{\nu^*}\rangle =0\). Observe this implies \(\nabla \cdot (e^{-V_{\nu^*}} \nabla \phi )=0\). Multiplying this equation by \(\phi\) and integrating by parts, we conclude that \(\int |\nabla \phi |^2 d\mu^* =0\). Thus we have \(\nabla \phi \equiv 0\), a contradiction. The remaining point is to upgrade this strict positivity to a uniform gap \(\lambda_{gap}>0\). This follows from a compactness argument which utilizes the gradient energy \(\int |\nabla^2 \phi|^2 d\mu^*\), or equivalently from a spectral theory for the elliptic operator \(L\) in 9 . Roughly, this means that the spectrum of \(L\) is contained in \([0,\infty)\) by 14 . If one thinks in terms of a discrete-spectrum picture, then the absence of \(0\) as an eigenvalue forces the smallest eigenvalue to be strictly positive.⁸, which yields 11 . Although the underlying idea is standard in the spectral theory of elliptic equations, we include a self-contained argument tailored to our setting. The same argument applies to the \(\nu\)-variable and gives 12 . See Remark 18 for characterization of optimal (the largest) \(\lambda_{gap}>0\) make the lemma hold and an alternative proof based on the spectral theory of weighted scalar Laplacian. ◻

Proof. See Appendix 7.1. ◻

Proof. In view of Lemma 4 (smoothing), \(\Vert \mu_t \Vert_{C^{1,\frac{1}{2}}}+\Vert \nu_t \Vert _{C^{1,\frac{1}{2}}} \le C=C(f)<\infty\) for all \(t\ge1\). Choose and fix some \(p>n\) and \(\lambda \in (0,\lambda_{gap})\). Let \(\varepsilon=\varepsilon(f,p,\lambda)>0\) be the constant from Proposition 3 (local contraction). By Lemma 5 (interpolation), Corollary 1 (\(C^{1,\alpha}\hookrightarrow W^{1,p}\)), and the previous observations, there exists \(\delta=\delta(f,p,\lambda)>0\) such that if \[\label{eq-profthm2-1}W_2^2(\mu_t,\mu^*)+W_2^2(\nu_t,\nu^*)\le \delta\tag{17}\] at some \(t=t'\ge1\), then ?? holds at this time \(t=t'\). Let us conclude the proof. Suppose \(W_2^2(\mu_0,\mu^*)+W_2^2(\nu_0,\nu^*)< e^{2\lambda '}\delta\), where \(\lambda'=\lambda'(\Vert f \Vert _{C^2})>-\infty\) is such that \(\nabla^2_{xx}f\ge \lambda'I\) and \(\nabla^2_{yy} f\le -\lambda'I\) at every \((x,y)\in ({ \mathbb{T}^n })^2\). By Lemma 3 (stability) and the previous observations, 17 holds for all \(t\ge1\) and ?? holds for all a.e. \(t\ge1\). In particular, this implies \[W_2^2(\mu_t,\mu^*)+W_2^2(\nu_t,\nu^*) < e^{-2\lambda'}(W_2^2(\mu_0,\mu^*)+W_2^2(\nu_0,\nu^*))e^{-2\lambda (t-1)},\] and thus 4 follows. Convergence in any \(C^{1,\alpha}\)-norm follows by Lemmas 4 and 5. ◻

Proof of Theorem 2. We first estimate the process stopped before it exits the stability neighborhood. On \(\{\tau>m\}\), the empirical measures stay in \(\mathcal{D}\) at the discrete times \(0,T,\ldots,mT\), and hence the one-step contraction can be iterated. The stopping is then removed by controlling the exit probability \(\mathbb{P}(\tau\le m)\). The immediate consequence of Proposition 4 is that, for all \(N\ge 1\) and \(t\ge 0\), taking \(m= \lfloor t/T\rfloor\)¹⁰ and using that \(\mathbf{1}_{\{\tau>k+1\}} \leqslant \mathbf{1}_{\{\tau > k\}}\) for all \(k\in\mathbb{N}\), \[\begin{align} \mathbb{E} \left ( r(\mu_{t}^N,\nu_{t}^N) \mathbf{1}_{\{ \tau > {m} \}} \right) &\le C_0 \mathbb{E} \left ( r(\mu_{mT}^N,\nu_{mT}^N) \mathbf{1}_{\{\tau > {m}\}} \right) + C' N^{-\eta} \\ &\le C_0 2^{-m} \mathbb{E} \left ( r(\mu_{0}^N,\nu_{0}^N) \mathbf{1}_{\{\tau > 0\}} \right) + 3 C' N^{-\eta} \\ &\le 2C_0 e^{-\lambda_1 t} \mathbb{E} \left ( r(\mu_{0}^N,\nu_{0}^N) \mathbf{1}_{\{\tau > 0\}} \right) + 3 C' N^{-\eta} \end{align}\] with \(\lambda_1 = \ln(2)/T\). Writing \(D_{\mathbb{T}^n}\) for the diameter of \(\mathbb{T}^n\), we immediately deduce that \[\begin{align} \mathbb{E} \left ( r(\mu_{t}^N,\nu_{t}^N) \right) &\le \mathbb{E} \left ( r(\mu_{t}^N,\nu_{t}^N) \mathbf{1}_{\{\tau > {m} \}} \right) + 2D_{\mathbb{T}^n} \mathbb{P}(\tau \le m ) \\ &\le C_1 e^{-\lambda_1 t} \mathbb{E} \left ( r(\mu_{0}^N,\nu_{0}^N) \right) + 3 C' N^{-\eta} + 2D_{\mathbb{T}^n} \mathbb{P}(\tau \le m ) \,. \end{align}\] It only remains to control the last term to conclude the proof of Theorem 2. First, \[\mathbb{P}(\tau \le m ) = \mathbb{P}(\tau \le m , \tau >0 ) + \mathbb{P}(\tau = 0 )\,.\] The last term is bounded by \(R_0 N^{-\eta}\) according to 7 (choosing \(\delta = \delta_0\) in Theorem 2). The proof is then concluded with Proposition 5 stated below. ◻

Proof of Proposition 5. Set \(\varepsilon = \delta_0/2\). For \(k \in \mathbb{N}\), consider the events \[E_k = \Big\{\sup_{t\in[kT,(k+1)T]} W_2(\mu_t^N,\mu_{k,t}^{\rm mf})+ W_2(\nu_t^N,\nu_{k,t}^{\rm mf})\le \varepsilon\Big \}, \qquad F_k = \left\{ (\mu_{kT}^N,\nu_{kT}^N) \in \mathcal{D} \right\},\] where \((\mu_{k,t}^{\rm mf},\nu_{k,t}^{\rm mf})\) is defined as previously \((\mu_{t}^{\rm mf},{\nu_{t}^{\rm mf}})\) but starting at time \(kT\) (i.e. they are the solutions of 2 initialized at \((\mu_{kT}^N,\nu_{kT}^N)\) at time \(kT\)). Under the event \(F_{k-1} \cap E_{k-1}\), applying Theorem 1 and the triangular inequality, \[\begin{align} r(\mu_{kT}^N,\nu_{kT}^N) & \leqslant r(\mu_{k-1,kT}^{\rm mf},\nu_{k-1,kT}^{\rm mf}) + W_2(\mu_{kT}^N,\mu_{k-1,kT}^{\rm mf})+ W_2(\nu_{kT}^N,\nu_{k-1,kT}^{\rm mf}) \\ &\le \frac{1}{2} r(\mu_{(k-1)T}^N,\nu_{(k-1)T}^N) + \varepsilon \;\le \;\delta_0, \end{align}\] namely \(F_k\) holds. As a consequence, using that \(\{\tau>0\} = F_0\), for any \(m\ge1\) \[\mathbb{P} \left( \tau >0, \tau \le m \right) \le \mathbb{P}\left(\tau >0 , \cup_{k=0}^{m-1} E_k^c\right) \le \sum_{k=0}^{m-1} \mathbb{P}(E_k^c) \le m e^{-bN}\] thanks to Lemma 6, which concludes the proof. ◻

Proof. Let \(\{u_m\}\) be a bounded sequence in \(W^{k+1, p}({ \mathbb{T}^n })\). By definition of the Sobolev norm, the sequence of derivatives \(\{D^\alpha u_m\}\) is bounded in \(W^{1, p}({ \mathbb{T}^n })\) for all multi-indices \(|\alpha| \le k\).

Applying Theorem 10 (specifically the case for \(L^p\) convergence) to each derivative:

1. For \(\alpha = 0\), \(\{u_m\}\) is bounded in \(W^{1,p}\), so there exists a subsequence converging strongly in \(L^p\).

2. For any \(|\alpha| \le k\), \(\{D^\alpha u_m\}\) is bounded in \(W^{1,p}\), so it also admits a subsequence converging strongly in \(L^p\).

Since there are finitely many multi-indices with \(|\alpha| \le k\), we can extract a single subsequence (via a diagonal argument) such that \(D^\alpha u_m\) converges strongly in \(L^p\) for all \(|\alpha| \le k\). By the definition of the Sobolev norm \[\|u\|_{W^{k,p}} = \left( \sum_{|\alpha| \le k} \|D^\alpha u\|_{L^p}^p \right)^{1/p},\] this immediately implies strong convergence in \(W^{k, p}({ \mathbb{T}^n })\). ◻

Proof. We divide the proof into several steps.

Step 1: (fixed point setting on path-space) Let \(\Omega:=C([0,T],M)\) and define, for \(\rho_1,\rho_2\in\mathcal{P}(\Omega)\) and \(t\in[0,T]\), \[\label{eq:W2t2} W_{2,t}^2(\rho_1,\rho_2) :=\inf_{\pi\in\Pi(\rho_1,\rho_2)} \int_{\Omega\times\Omega}\Big(\sup_{0\le s\le t} d^2(X_1(s),X_2(s))\Big)\,d\pi(X_1,X_2).\tag{23}\]

For \(\rho\in\mathcal{P}(\Omega)\) and \(t\in[0,T]\), we denote by \(\rho_t\in\mathcal{P}(M)\) the one-time marginal at time \(t\), i.e.the law of \(\omega(t)\) when \(\omega\sim\rho\). Equivalently, \[\rho_t:=(\mathrm{ev}_t)_\#\rho\in\mathcal{P}(M).\] where \(\mathrm{ev}_t:\Omega\to M\) is the evaluation map defined by \(\mathrm{ev}_t(\omega):=\omega(t)\).

Given \(\eta\in\mathcal{P}(\Omega)\) and \(\rho\in\mathcal{P}(\Omega)\), define \[b^\eta(t,x):=-\nabla_1\int_M f(x,y)\,d\eta_t(y), \qquad c^\rho(t,y):=\nabla_2\int_M f(x,y)\,d\rho_t(x).\] By 19 , for each fixed \(\eta\) the map \(x\mapsto b^\eta(t,x)\) is globally Lipschitz uniformly in \(t\) (and similarly for \(c^\rho\)). Consider the SDEs on \(M\): \[\label{eq:SDE95frozen2} dX_t=b^\eta(t,X_t)\,dt+\sqrt{2}\,dW_t,\qquad X_0\sim\mu_0,\tag{24}\] \[\label{eq:SDE95frozenY2} dY_t=c^\rho(t,Y_t)\,dt+\sqrt{2}\,d\bar W_t,\qquad Y_0\sim\nu_0.\tag{25}\] (For concreteness: one may lift to \(\mathbb{R}^n\) with periodic drifts and then project modulo \(\mathbb{Z}^n\).) Each equation has a unique strong solution, hence defines path laws \[\Phi_1(\eta):=\mathop{\mathrm{Law}}(X_\cdot)\in\mathcal{P}(\Omega),\qquad \Phi_2(\rho):=\mathop{\mathrm{Law}}(Y_\cdot)\in\mathcal{P}(\Omega),\] and the coupled map \[\Phi(\rho,\eta):=\big(\Phi_1(\eta),\Phi_2(\rho)\big)\qquad\text{on }\mathcal{P}(\Omega)\times\mathcal{P}(\Omega).\]

Step 2: (Picard fixed point) We show that the map \(\Phi\) admits a unique fixed point \((\rho,\eta)\in\mathcal{P}(\Omega)\times\mathcal{P}(\Omega)\); First, we prove the following lemma for stability estimate:

Pick any \((\rho^{(0)},\eta^{(0)})\) and define iterates \((\rho^{(k+1)},\eta^{(k+1)}):=\Phi(\rho^{(k)},\eta^{(k)})\). Set \[d_k(t):=W_{2,t}^2(\rho^{(k+1)},\rho^{(k)})+W_{2,t}^2(\eta^{(k+1)},\eta^{(k)}),\qquad t\in[0,T].\] Applying Lemma 8 with \((\rho,\eta)=(\rho^{(k)},\eta^{(k)})\) and \((\tilde{\rho},\tilde{\eta})=(\rho^{(k-1)},\eta^{(k-1)})\) yields \[d_k(t)\le C\int_0^t d_{k-1}(s)\,ds.\] Iterating this inequality gives the factorial bound \[d_k(t)\le \frac{(Ct)^k}{k!}\,d_0(t),\qquad t\in[0,T],\] hence \(\sum_{k\ge0} d_k(T)<\infty\). Therefore \((\rho^{(k)},\eta^{(k)})\) is Cauchy for \(W_{2,T}\) and converges to a limit \((\rho,\eta)\), which is easily checked to satisfy \(\Phi(\rho,\eta)=(\rho,\eta)\).

For uniqueness, let \((\rho,\eta)\) and \((\tilde{\rho},\tilde{\eta})\) be two fixed points and set \(D(t):=W_{2,t}^2(\rho,\tilde{\rho})+W_{2,t}^2(\eta,\tilde{\eta})\). Then Lemma 8 gives \(D(t)\le C\int_0^t D(s)\,ds\), so \(D(t)\equiv 0\) by Grönwall.

Step 3: (2 is uniquely solved by \((\mu_t,\nu_t)\)) Let \((\rho,\eta)\) be the unique fixed point from Step 2. Let \(X\) and \(Y\) be the corresponding solutions of 24 –25 with \(\mathop{\mathrm{Law}}(X_\cdot)=\rho\) and \(\mathop{\mathrm{Law}}(Y_\cdot)=\eta\), and define the time marginals \[\mu_t:=\mathop{\mathrm{Law}}(X_t)=(\mathrm{ev}_t)_\#\rho,\qquad \nu_t:=\mathop{\mathrm{Law}}(Y_t)=(\mathrm{ev}_t)_\#\eta.\] Then \((\mu_t,\nu_t)\) is narrowly continuous and satisfies \(\mu_{t=0}=\mu_0\), \(\nu_{t=0}=\nu_0\). For \(\varphi\in C^2(M)\), Itô’s formula for 24 gives \[d\varphi(X_t)=\nabla\varphi(X_t)\cdot b^\eta(t,X_t)\,dt+\Delta\varphi(X_t)\,dt+\sqrt{2}\,\nabla\varphi(X_t)\cdot dW_t.\] Taking expectations and using \(\nu_t=\eta_t\) yields 20 . The same argument applied to \(Y_t\) yields 21 . Thus \((\mu_t,\nu_t)\) is a distributional solution of 2 .

To show uniqueness, let \((\tilde{\mu}_t,\tilde{\nu}_t)_{t\in[0,T]}\) be another distributional solution of 2 with the same initial data. Define \[\tilde{b}_t(x):=-\nabla_1\int_M f(x,y)\,d\tilde{\nu}_t(y),\qquad \tilde{c}_t(y):=\nabla_2\int_M f(x,y)\,d\tilde{\mu}_t(x).\] Since \(f\in C^2(M\times M)\) and \(M\) is compact, we have \(\tilde{b},\tilde{c}\in L^\infty([0,T];W^{1,\infty}(M))\).

By definition of a distributional solution 20 21 , the curves \((\tilde{\mu}_t)\) and \((\tilde{\nu}_t)\) solve the linear Fokker–Planck equations \[\partial_t \tilde{\mu} = \Delta\tilde{\mu} - \nabla\cdot(\tilde{b}_t\,\tilde{\mu}),\qquad \partial_t \tilde{\nu} = \Delta\tilde{\nu} - \nabla\cdot(\tilde{c}_t\,\tilde{\nu}).\] in the sense of distributions.

By the superposition principle for Fokker-Planck equations (see e.g.[42]), there exist \(\tilde{\rho},\tilde{\eta}\in\mathcal{P}(\Omega)\) such that \[(\mathrm{ev}_t)_\#\tilde{\rho}=\tilde{\mu}_t,\qquad (\mathrm{ev}_t)_\#\tilde{\eta}=\tilde{\nu}_t \quad\text{for all }t\in[0,T],\] and \(\tilde{\rho}\) (resp.\(\tilde{\eta}\)) is the law of a weak solution to \[dX_t=\tilde{b}_t(X_t)\,dt+\sqrt2\,dW_t,\qquad X_0\sim\mu_0,\] (resp.to \(dY_t=\tilde{c}_t(Y_t)\,dt+\sqrt2\,d\bar W_t,\quad Y_0\sim\nu_0\) ).

On the other hand, since \(\tilde{b},\tilde{c}\) are bounded and globally Lipschitz in space (uniformly in \(t\)), these SDEs are well posed in law. Therefore the above path laws coincide with the unique laws defining \(\Phi_1(\tilde{\eta})\) and \(\Phi_2(\tilde{\rho})\), and hence \[\tilde{\rho}=\Phi_1(\tilde{\eta}),\qquad \tilde{\eta}=\Phi_2(\tilde{\rho}).\] Thus \((\tilde{\rho},\tilde{\eta})\) is a fixed point of \(\Phi\). By Step 2 the fixed point is unique, so \((\tilde{\rho},\tilde{\eta})=(\rho,\eta)\) and consequently \(\tilde{\mu}_t=\mu_t\) and \(\tilde{\nu}_t=\nu_t\) for all \(t\in[0,T]\).

Step 4: (instantaneous regularization for \(t>0\)) Let \((\mu_t,\nu_t)\) be the unique distributional solution from the previous steps and define \[b_t(x):=-\int_M \nabla_1 f(x,y)\,d\nu_t(y),\qquad c_t(y):=\int_M \nabla_2 f(x,y)\,d\mu_t(x).\] Since \(\mu_t,\nu_t\) are probability measures and \(f\in C^2(M\times M)\), \[\label{eq:drift95bounds95revised} \|b_t\|_{L^\infty(M)}+\|c_t\|_{L^\infty(M)}\le \|\nabla f\|_{L^\infty(M\times M)}, \qquad \|\nabla b_t\|_{L^\infty(M)}+\|\nabla c_t\|_{L^\infty(M)}\le \|\nabla^2 f\|_{L^\infty(M\times M)}.\tag{27}\]

Fix \(\tau\in(0,T)\). Consider the linear Fokker-Planck equation \(\partial_t\mu=\Delta\mu + \nabla\cdot(V_t\mu)\) with \(V\in L^\infty_tW^{1,\infty}_x\). Uniform ellipticity of \(\Delta\) and boundedness of \(b\) imply that the associated Markov semigroup admits a transition density for each \(t>0\); equivalently, for any initial probability measure, \(\mu_t\) is absolutely continuous for \(t>0\). The same holds for \(\nu_t\). Hence there exist densities \(u,v\) on \((\tau,T]\times M\) with \(\mu_t=u(t,\cdot)\,dx\) and \(\nu_t=v(t,\cdot)\,dx\).

On \((\tau,T)\times M\), the densities solve (in the weak sense) the linear parabolic equations \[\label{eq:density95form95revised} \begin{cases} \partial_t u - \Delta u = -\nabla\cdot(u\,b_t),\\ \partial_t v - \Delta v = -\nabla\cdot(v\,c_t). \end{cases}\tag{28}\] Because \(b,c\in L^\infty((\tau,T);W^{1,\infty}(M))\), standard parabolic maximal \(L^p\) regularity theory on the torus yields 22 for every \(p\in(1,\infty)\), with constants depending on \(p,\tau,T,M\) and the bounds in 27 . (One may first note that \(u(\tau,\cdot),v(\tau,\cdot)\in L^p(M)\) for every \(p<\infty\) by smoothing of the elliptic part together with the bounded drift; the resulting bounds deteriorate as \(\tau\downarrow 0\).)

From 22 with \(p\) large, Sobolev embedding gives Hölder regularity of \(u,v\) on \((\tau,T]\times M\). Under the standing assumption \(f\in C^2\), the bounds 22 are the (safe) conclusion needed here. If one additionally assumes \(f\in C^\infty(M\times M)\) (or \(C^{k+2}\) for finite \(k\)), then a standard bootstrap argument yields \(u,v\in C^\infty((\tau,T]\times M)\) (resp.\(W^{k+2,p}\)-regularity). This completes the proof of Theorem 16. ◻

Proof. Write \(\mu=\rho\,dx\) and \(\nu=\eta\,dx\), where \(\rho,\eta\in W^{k,p}({ \mathbb{T}^n })\) with \(0<c\le \rho,\eta \le C\). Since \(\nu=T_\#\mu\) with \(T=id+\nabla\psi\), the change-of-variables formula yields \[\label{eq-changeofvar} \eta = \frac{\rho\circ T^{-1}}{\det(\nabla T)\circ T^{-1}} = (\rho\circ T^{-1})\,\det(\nabla T^{-1}), \qquad \nabla T = I+\nabla^2\psi.\tag{29}\]

We decompose \[\eta-\rho = (\rho\circ T^{-1}-\rho)\,\det(\nabla T^{-1}) + \rho\,(\det(\nabla T^{-1})-1).\] Since \(p>n\), the Sobolev space \(W^{k,p}({ \mathbb{T}^n })\) is a Banach algebra, hence \[\begin{align} \|\eta-\rho\|_{W^{k,p}} &\le C\Bigl( \|\rho\circ T^{-1}-\rho\|_{W^{k,p}} + \|\rho\|_{W^{k,p}}\, \|\det(\nabla T^{-1})-1\|_{W^{k,p}} \Bigr). \end{align}\]

We first estimate the composition term. Since \(\|\nabla\psi\|_{W^{k+1,p}}\) is sufficiently small, \(T\) is a \(W^{k+1,p}\)-diffeomorphism with inverse \(T^{-1}=id+u\), where \[\|u\|_{W^{k+1,p}}\le C\|\nabla\psi\|_{W^{k+1,p}}.\] By standard composition estimates in Sobolev spaces for \(p>n\), \[\|\rho\circ(id+u)-\rho\|_{W^{k,p}} \le C\,\|u\|_{W^{k,p}}\,\|\rho\|_{W^{k+1,p}} \le C\,\|\nabla\psi\|_{W^{k+1,p}}.\]

Next, we estimate the Jacobian term. Using the chain rule, \[\nabla T^{-1} = (\nabla T)^{-1}\circ T^{-1}, \qquad \nabla T = I + A,\;\;A=\nabla^2\psi.\] The smallness of \(\|A\|_{W^{k,p}}\) implies \[(\nabla T)^{-1}=I+\widetilde{A}, \qquad \|\widetilde{A}\|_{W^{k,p}}\le C\|A\|_{W^{k,p}}.\] Since the determinant is a smooth polynomial function of the matrix entries, it is locally Lipschitz in \(W^{k,p}\), and therefore \[\|\det(\nabla T^{-1})-1\|_{W^{k,p}} \le C\|\widetilde{A}\|_{W^{k,p}} \le C\|\nabla^2\psi\|_{W^{k,p}} \le C\|\nabla\psi\|_{W^{k+1,p}}.\]

Combining the above estimates, we conclude that \[\|\mu-\nu\|_{W^{k,p}} = \|\eta-\rho\|_{W^{k,p}} \le C\,\|\nabla\psi\|_{W^{k+1,p}},\] where the constant \(C\) depends only on \(n,p,k\), \(c,C\), and \(\|\rho\|_{W^{k+1,p}}\). This proves the claim. ◻

Proof. We prove this by a perturbative inverse function theorem argument around the identity map.

For \(\phi\in W^{k+2,p}({ \mathbb{T}^n })\) with \(\int_{{ \mathbb{T}^n }}\phi\,dx=0\), consider the map \(S:=id+\nabla\phi\). We define the operator \(J(\phi)\) to be the density which is pushed forward to \(\rho\) by the map \(S\). According to the Monge-Ampère equation 29 , this corresponds to the pull-back of \(\rho\) by \(S\): \[J(\phi)(x) := \rho(S(x))\,\det(\nabla S(x)) = \rho(x+\nabla\phi(x))\,\det(I+\nabla^2\phi(x)).\] Since \(p>n\), the Sobolev space \(W^{k+1,p}({ \mathbb{T}^n })\) is a Banach algebra. The assumption \(\rho\in W^{k+1,p}({ \mathbb{T}^n })\) ensures that \(J(\phi)\in W^{k,p}({ \mathbb{T}^n })\). Standard composition and determinant estimates imply that \(J:W^{k+2,p}({ \mathbb{T}^n })\to W^{k,p}({ \mathbb{T}^n })\) is of class \(C^1\) in a neighborhood of \(\phi=0\).

By construction, the equation \(J(\phi)=\rho_1\) is equivalent to the mass transport constraint \[S_\#(\rho_1\,dx)=\rho\,dx.\] Our strategy is to find such a map \(S\) for \(\rho_1\) close to \(\rho\), and then obtain the optimal transport map \(T\) from \(\mu\) to \(\mu_1\) as the inverse map \(T=S^{-1}\).

A direct computation shows that the Fréchet derivative of \(J\) at \(\phi=0\) is given by \[dJ(0)[\varphi] = \nabla\cdot(\rho\nabla\varphi),\] for \(\varphi\in W^{k+2,p}({ \mathbb{T}^n })\) with zero mean. We claim that the linearized operator \[L: \Bigl\{\varphi\in W^{k+2,p}({ \mathbb{T}^n }): \int_{{ \mathbb{T}^n }}\varphi\,dx=0\Bigr\} \;\longrightarrow\; \Bigl\{h\in W^{k,p}({ \mathbb{T}^n }): \int_{{ \mathbb{T}^n }}h\,dx=0\Bigr\}, \quad L\varphi:=\nabla\cdot(\rho\nabla\varphi)\] is an isomorphism. Indeed, \(L\) is a self-adjoint, uniformly elliptic operator with coefficients in \(W^{k+1,p}({ \mathbb{T}^n })\). Since we restrict the domain to mean-zero functions, the Poincaré inequality implies coercivity. Thus, by the Lax-Milgram theorem and standard elliptic regularity [39], \(L\) admits a bounded inverse.

Since the target space \[M:=\Bigl\{\tilde{\rho}\in W^{k,p}({ \mathbb{T}^n }): \int_{{ \mathbb{T}^n }}\tilde{\rho}\,dx=1\Bigr\}\] is an affine subspace of \(W^{k,p}({ \mathbb{T}^n })\), the inverse function theorem applies directly. Thus there exist \(\varepsilon>0\) and \(C<\infty\) such that for every \(\rho_1\in M\) satisfying \(\|\rho_1-\rho\|_{W^{k,p}}\le\varepsilon\), there exists a unique \(\phi\in W^{k+2,p}({ \mathbb{T}^n })\) with \(\int_{{ \mathbb{T}^n }}\phi\,dx=0\) such that \(J(\phi)=\rho_1\), and \[\|\phi\|_{W^{k+2,p}} \le C\,\|\rho_1-\rho\|_{W^{k,p}}.\]

From Step 1, the solution \(\phi\) defines a map \(S=id+\nabla\phi\) satisfying \(S_\#(\rho_1\,dx)=\rho\,dx\). For \(\varepsilon\) sufficiently small, \(\|\nabla^2\phi\|_{L^\infty}\) is small, so \(S\) is the gradient of a strictly convex function. By the uniqueness of the Brenier map, \(S\) is the optimal transport map from \(\mu_1\) to \(\mu\).

Consequently, the optimal transport map \(T\) from \(\mu\) to \(\mu_1\) is given by the inverse \(T=S^{-1}\). From the properties of convex conjugates, \(T\) is also the gradient of a convex function, \(T=id+\nabla\psi\). Applying standard stability estimates for the inversion of diffeomorphisms close to the identity, we obtain \[\|\nabla\psi\|_{W^{k+1,p}} \le C'\,\|\nabla\phi\|_{W^{k+1,p}}.\] Combining this with the estimate in Step 3 yields \[\|\nabla\psi\|_{W^{k+1,p}} \le C' C\,\|\rho_1-\rho\|_{W^{k,p}}.\] This completes the proof. ◻

Proof. Since \(\Phi=\nabla\phi \in W^{2,p}({ \mathbb{T}^n })\) with \(p>n\), Morrey’s inequality (Theorem 6) guarantees that \(\Phi \in C^{1,\alpha}({ \mathbb{T}^n })\) for some \(\alpha > 0\). Consequently, \(\nabla \Phi\) is continuous and uniformly bounded on \({ \mathbb{T}^n }\), which implies for sufficiently small \(t\), thanks to the continuity of the determinant, the Jacobian determinant satisfies \(\det(I_n + t\nabla \Phi(x)) \ge 1/2>0\) for all \(x \in { \mathbb{T}^n }\). Moreover, it is straightforward to check that \(T_t(x)=x+t\Phi(x)\) is a \(C^{1,\alpha}\)-diffeomorphism for sufficiently small \(t\), since, after lifting \(T_t\) periodically to \(\mathbb{R}^n\), it holds \[|T_t(x)-T_t(y)|\ge (1-|t|\|\nabla\Phi\|_{L^\infty})|x-y| \qquad\forall x,y\in\mathbb{R}^n,\] so \(T_t\) is injective (hence bi-Lipschitz onto its image) for \(|t|\|\nabla\Phi\|_{L^\infty}<1\). Combining this with the invertibility of \(\nabla T_t\), we conclude that \(T_t\) is a \(C^{1,\alpha}\)-diffeomorphism of \(\mathbb{T}^n\).

Let \(\mu_t = (id+t\Phi)_\#\mu\). Using the change of variables formula \(y = x + t\Phi(x)\), the functional \(\widetilde{G}(t):= G(\mu_t)\) transforms to (unless otherwise stated, we omit the domain \({ \mathbb{T}^n }\) in the integral sign): \[\begin{align} \widetilde{G}(t) &= \int \rho_t(y) \log \rho_t(y) \, dy + \int V(y) \, d\mu_t(y) \\ &= \int \frac{\rho(x)}{\det(I_n + t\nabla \Phi)} \log \left( \frac{\rho(x)}{\det(I_n + t\nabla \Phi)} \right) \det(I_n + t\nabla \Phi) \, dx + \int V(x+t\Phi) \rho(x) \, dx \\ & = \int \log \rho \, d\mu - \int \log \det(I_n + t\nabla \Phi) \, d\mu + \int V(x+t\Phi) \, d\mu \end{align}\] Note that the condition \(0 < c \le \rho \le C\) ensures that the entropy term is well-defined and finite.

We show that the first and second variations exist and coincide with the stated formulas. As we saw earlier, there exists \(\bar{t}>0\) such that if \(|t|<\bar{t}\), \((I_n +t \nabla\Phi)\) is invertible (pointwise). Then for all \(|t|<\bar{t}\) we obtain \[\frac{d}{dt} \log (\det A_t) = \textrm{tr}(A_t^{-1}\dot{A}_t), \quad A_t=I_n + t\nabla\Phi.\]

Given \(t_0\) with \(|t_0|< \bar{t}\), we have \[\label{forDCT} \Big|\frac{ \log \det(I_n + t \nabla\Phi) - \log \det (I_n +t_0 \nabla\Phi ) }{t-t_0} \Big| \leq C \Vert \nabla\Phi \Vert_{L^\infty}\tag{30}\] for all \(t\) sufficiently close to \(t_0\), uniformly in \(x\in { \mathbb{T}^n }\). By the dominated convergence theorem, we have for all \(|t|<\bar{t}\) \[\frac{d}{dt} \widetilde{G}( t ) =-\int \textrm{tr}( (I_n + t \nabla\Phi)^{-1} \nabla\Phi \big) d\mu + \int \nabla V(id +t \Phi) \cdot \Phi d\mu\] In particular, we get \[\begin{align} \label{energydiff} \frac{d}{dt}\Big|_{t=0} \widetilde{G}(t) &= -\int \rho (\nabla \cdot \Phi) \, dx + \int (\nabla V \cdot \Phi) \rho \, dx =\int \nabla\rho \cdot \Phi \,dx + \int \nabla V \cdot \Phi d \mu { \nonumber }\\ &=\int \nabla \log \rho \cdot \Phi \, d\mu +\int \nabla V \cdot \Phi d \mu. { \nonumber } \end{align}\tag{31}\]

The second variation formula is derived as follows: By Morrey’s inequality, \(W^{2,p}({ \mathbb{T}^n })\hookrightarrow C^{1,\alpha}({ \mathbb{T}^n })\) for some \(\alpha\in(0,1)\) since \(p>n\). Hence \(\Phi,\Psi\in C^{1,\alpha}({ \mathbb{T}^n })\), and for \(|s|+|t|\) sufficiently small the map \(F_{s,t}=id+t\Phi+s\Psi\) is a \(C^{1,\alpha}\)-diffeomorphism of \({ \mathbb{T}^n }\).

Let \(\mu_{s,t}:=(id+t\Phi+s\Psi)_{\#}\mu\). Then the pushforward density of \((F_{s,t})_{\#}\mu\), say \(\rho_{s,t}\), satisfies \[\rho_{s,t}(F_{s,t}(x))\, \det\bigl(I_n +t\nabla\Phi(x)+s\nabla\Psi(x)\bigr) =\rho(x).\] Therefore, by change of variables, \[\begin{align} & H(s,t):=G(\mu_{s,t}) \\ &= \int \rho(x)\log\!\Bigl(\frac{\rho(x)}{\det(I_n+t\nabla\Phi(x)+s\nabla\Psi(x))}\Bigr)\,dx +\int V(x+t\Phi(x)+s\Psi(x))\,\rho(x)\,dx \\ &= \int \rho\log\rho\,dx -\int \rho(x)\log\det(I_n+t\nabla\Phi(x)+s\nabla\Psi(x))\,dx \\ &\qquad +\int V(x+t\Phi(x)+s\Psi(x))\,\rho(x)\,dx . \end{align}\] Similarly as in 30 , we see that \[\frac{\partial}{\partial t}\Big|_{t=0} \int - \rho \log \det(I_n +t \nabla\Phi + s \nabla\Psi) = - \int \textrm{tr}( (I_n+s \nabla\Psi)^{-1} \nabla\Phi) \rho dx.\] Therefore, we have \[\begin{align} &\frac{\partial^2}{\partial s\,\partial t}\Big|_{t=s=0} \int -\rho(x)\log\det(I_n +t\nabla\Phi(x)+s\nabla\Psi(x))\,dx =\int \textrm{tr}(\nabla\Phi(x) \nabla\Psi(x)) \rho(x) dx { \nonumber }\\ &= \int \rho(x)\,\nabla\Phi(x):\nabla\Psi(x)\,dx, { \nonumber } \end{align}\] where the last equality follows from the fact that \(\nabla\Phi, \nabla\Psi\) are symmetric. Next, since \(x+t\Phi(x)+s\Psi(x)\) is affine in \((s,t)\), the chain rule yields \[\frac{\partial^2}{\partial s\,\partial t}\Big|_{t=s=0} V(x+t\Phi(x)+s\Psi(x)) = \nabla^2V(x)[\Phi(x),\Psi(x)].\] Thus \[\frac{\partial^2}{\partial s\,\partial t}\Big|_{t=s=0} \int V(x+t\Phi(x)+s\Psi(x))\,\rho(x)\,dx = \int \nabla^2V(x)[\Phi(x),\Psi(x)]\,\rho(x)\,dx.\] Combining the two identities, we obtain \[\frac{\partial^2}{\partial s\,\partial t}\Big|_{t=s=0} H(s,t) = \int \Bigl( \nabla\Phi:\nabla\Psi+\nabla^2V[\Phi,\Psi] \Bigr)\,d\mu.\] Finally, by definition of \(L\) and integration by parts, \[\begin{align} \langle L\Phi,\Psi\rangle_{L^2_\mu} &= \int \Bigl( -\rho^{-1}\nabla\cdot(\rho\nabla\Phi)+\nabla^2V\cdot\Phi \Bigr)\cdot \Psi\,d\mu \\ &= -\int \nabla\cdot(\rho\nabla\Phi)\cdot \Psi\,dx +\int \rho\,\nabla^2V[\Phi,\Psi]\,dx \\ &= \int \rho\,(\nabla\Phi:\nabla\Psi)\,dx +\int \rho\,\nabla^2V[\Phi,\Psi]\,dx \\ &= \int \Bigl( \nabla\Phi:\nabla\Psi+\nabla^2V[\Phi,\Psi] \Bigr)\,d\mu. \end{align}\] Since the right-hand side is symmetric in \((\Phi,\Psi)\), we also get \[\langle L\Phi,\Psi\rangle_{L^2_\mu} = \langle \Phi,L\Psi\rangle_{L^2_\mu}.\] This completes the proof. ◻

Proof. We use the explicit variational formula for the entropy-regularized best responses.

Step 1: Entropy duality identities. For any bounded measurable function \(U\) on \({ \mathbb{T}^n }\), we have the following: \[\label{eq:entropy-dual} \max_{\nu\in{\mathcal{P}}({ \mathbb{T}^n })} \Bigl\{\int U\,d\nu - H(\nu)\Bigr\} = \log\int_{{ \mathbb{T}^n }} e^{U(y)}\,dy,\tag{35}\] with maximizer \(\nu=Z^{-1}e^{U(y)}\,dy\), where \(Z=\int e^{U(y)}dy\). Similarly, for any bounded \(V\), \[\label{eq:entropy-dual-min} \min_{\mu\in{\mathcal{P}}({ \mathbb{T}^n })} \Bigl\{\int V\,d\mu + H(\mu)\Bigr\} = -\log\int_{{ \mathbb{T}^n }} e^{-V(x)}\,dx,\tag{36}\] with minimizer \(\mu= Z^{-1} e^{-V(x)}dx\), where \(Z=\int e^{-V(x)}dx\). We prove 35 : for any \(\nu\ll dx\), write \(\nu=\rho\,dx\) and compute \[\int g\,d\nu - H(\nu) = \int \rho(x)\Bigl(g(x)-\log\rho(x)\Bigr)\,dx.\] Maximizing under the constraint \(\int\rho\,dx=1\) yields \(\rho(x)\propto e^{g(x)}\), which gives 35 . The proof of 36 is analogous.

Step 2: Closed form of the NI error. For fixed \((\mu,\nu)\) define \[U_\mu(y):=\int_{{ \mathbb{T}^n }} f(x,y)\,d\mu(x), \qquad V_\nu(x):=\int_{{ \mathbb{T}^n }} f(x,y)\,d\nu(y).\] Applying 35 –36 to the definition of \(\mathrm{NI}\) yields \[\label{eq:NI-closed-form} \mathrm{NI}(\mu,\nu) = H(\mu)+H(\nu) +\log\int_{{ \mathbb{T}^n }} e^{U_\mu(y)}\,dy +\log\int_{{ \mathbb{T}^n }} e^{-V_\nu(x)}\,dx.\tag{37}\]

Step 3: Estimate.

Now we work under the assumption of the lemma. Note we have the pointwise estimates \(|\rho_1\log \rho_1 - \rho_2\log \rho_2| \le (1+ \log C_1) |\rho_1-\rho_2|\) and \(|U_{\mu_1}-U_{\mu_2}|\le \Vert f\Vert_{L^\infty}\Vert \rho_1-\rho_2\Vert _{L^\infty}\). These estimates, together with corresponding estimates for \(\nu\), implies \[|\mathrm{NI}(\mu_1,\nu_1)-\mathrm{NI}(\mu_2,\nu_2)|\le C'(\Vert \rho_1 -\rho_2\Vert_{L^\infty}+ \Vert \eta_1-\eta_2\Vert_{L^\infty}),\] for some \(C'=C'(C_1,\Vert f\Vert_{L^\infty})\). ◻

Proof. Set \(\mathcal{E}(t):=W_2^2(\mu_t,\mu^*)+W_2^2(\nu_t,\nu^*)\).

Step 1: (showing EVI once) the \(\mu\)-part. Fix \(t\) and consider the \(\lambda'\)-displacement convex functional \(\mu\mapsto F(\mu,\nu_t)\). The above-tangent inequality [40] yields \[\label{eq:AT-mu-once} F(\mu^*,\nu_t)-F(\mu_t,\nu_t) \ge \big\langle \mathop{\mathrm{grad}}_{W_2,\mu}F(\mu_t,\nu_t),\,T_{\mu_t\to\mu^*}-id\big\rangle_{L^2_{\mu_t}} +\frac{\lambda'}{2}W_2^2(\mu_t,\mu^*).\tag{38}\] On the other hand, since \(t\mapsto\mu_t\) is absolutely continuous and \(\partial_t\mu_t+\nabla\cdot(\mu_t v_t)=0\), the first variation formula [40], namely Theorem 15, gives for a.e.\(t\), \[\label{eq:1stvar-mu-once} \frac{d}{dt}\frac{1}{2} W_2^2(\mu_t,\mu^*) =\big\langle -v_t,\,T_{\mu_t\to\mu^*}-id\big\rangle_{L^2_{\mu_t}}.\tag{39}\] Using \(v_t=-\mathop{\mathrm{grad}}_{W_2,\mu}F(\mu_t,\nu_t)\) and combining 38 –39 , we obtain the EVI \[\label{eq:EDI-mu} \frac{d}{dt}\frac{1}{2} W_2^2(\mu_t,\mu^*) \le F(\mu^*,\nu_t)-F(\mu_t,\nu_t)-\frac{\lambda'}{2}W_2^2(\mu_t,\mu^*).\tag{40}\]

Step 2: the \(\nu\)-part (same argument). Since \(\nu\mapsto -F(\mu_t,\nu)\) is \(\lambda'\)-displacement convex and \(w_t=+\mathop{\mathrm{grad}}_{W_2,\nu}F(\mu_t,\nu_t)\), applying the same argument to the descent flow of \(-F(\mu_t,\cdot)\) yields for a.e.\(t\), \[\label{eq:EDI-nu} \frac{d}{dt}\frac{1}{2} W_2^2(\nu_t,\nu^*) \le F(\mu_t,\nu_t)-F(\mu_t,\nu^*)-\frac{\lambda'}{2}W_2^2(\nu_t,\nu^*).\tag{41}\]

Summing 40 and 41 gives \[\label{eq:lyap-pre} \frac{d}{dt}\frac{1}{2}\mathcal{E}(t) \le -\frac{\lambda'}{2}\mathcal{E}(t)+\Bigl(F(\mu^*,\nu_t)-F(\mu_t,\nu^*)\Bigr).\tag{42}\]

Step 3: improved cancellation at the stationary saddle point. Apply the above-tangent inequality to the \(\lambda'\)-displacement convex map \(\mu\mapsto F(\mu,\nu^*)\) at the point \(\mu^*\) and use \(\mathop{\mathrm{grad}}_{W_2,\mu}F(\mu^*,\nu^*)=0\) to obtain \[F(\mu_t,\nu^*)-F(\mu^*,\nu^*)\ge \frac{\lambda'}{2}W_2^2(\mu_t,\mu^*), \quad\text{i.e.}\quad F(\mu^*,\nu^*)-F(\mu_t,\nu^*)\le -\frac{\lambda'}{2}W_2^2(\mu_t,\mu^*).\] Similarly, since \(\nu\mapsto F(\mu^*,\nu)\) is \(\lambda'\)-concave and \(\mathop{\mathrm{grad}}_{W_2,\nu}F(\mu^*,\nu^*)=0\), applying above-tangent to \(-F(\mu^*,\cdot)\) yields \[F(\mu^*,\nu_t)-F(\mu^*,\nu^*)\le -\frac{\lambda'}{2}W_2^2(\nu_t,\nu^*).\] Adding these inequalities gives \[F(\mu^*,\nu_t)-F(\mu_t,\nu^*) \le -\frac{\lambda'}{2}\mathcal{E}(t).\] Combining with 42 we obtain \[\frac{d}{dt}\mathcal{E}(t)\le -2\lambda'\mathcal{E}(t),\] and Grönwall’s inequality yields the claim. ◻

Proof. We prove ?? . A proof of ?? follows by considering \(\hat{F} (\nu,\mu):= -F(\mu,\nu)\) and repeat the same proof of ?? ; see the last sentence of this proof. Let us define \(V_{\nu^*}(x):= \int_{\mathcal{Y}} f(x,y) d\nu^*(y)\) and \[G( \mu ) := F(\mu , \nu^*)= H(\mu)+\int V_{\nu^*}(x) d\mu(x) - H(\nu^*).\] Let us denote \(\frac{d \mu^*}{dx}=:\rho_*\). It is known that \(\rho_* \propto e^{-V_{\nu^*} }\) [44]. More precisely, \(\rho_*= e^{-V_{\nu^*}} / \int e^{-V_{\nu^*}}dx\). Similarly, \(\frac{d \nu^*}{dy}\) satisfies \(\frac{d \nu^*}{dy} = e^{+U_{\mu^*}}/ \int e^{+U_{\mu^*}}dy\) for \(U_{\mu^*}(y):= \int _{\mathcal{X}}f(x,y) d\mu^*(x)\). We also note that the MNE \((\mu^*,\nu^*)\) has \(C^2\) densities: \[\label{eq-C942bound} C^{-1}\le \frac{d\mu^*}{dx}, \, \frac{d\nu^*}{dy} \le C, \quad \Big \Vert \frac{d\mu^*}{dx}\Big \Vert_{C^2},\Big \Vert \frac{d\nu^*}{dy}\Big \Vert_{C^2} \le C, \quad \text{ for some }C=C(\Vert f\Vert_{C^2})<\infty,\tag{43}\] which follows from the observation \(\Vert U_{\mu^*}\Vert_{C^2(\mathcal{Y})}\), \(\Vert V_{\nu^*}\Vert_{C^2(\mathcal{X})}\le\Vert f\Vert_{C^2(\mathcal{X}\times \mathcal{Y})}\).

For \(\nabla \phi \in W^{2,p}\) for \(p>n\), we have the second variation formula \[\label{eq-secondvariation}\frac{d^2}{d t ^2}\Big\vert_{t=0} G((id+t\nabla \phi )_\# \mu^* )= \int (|\nabla^2 \phi|^2 + \nabla^2 V_{\nu^*}[\nabla\phi,\nabla\phi]) \, \rho_* dx .\tag{44}\] Here, \(|\nabla ^2 \phi|^2= \sum_{i,j} \nabla^2_{ij}\phi \nabla^2_{ij}\phi\) and \(\nabla^2 V_{\nu^*}[\nabla\phi,\nabla\phi] = \sum_{i,j} \nabla^2_{ij}V_{\nu^*} \nabla _i\phi \nabla_j\phi = \langle \nabla^2 V_{\nu^*} \nabla \phi ,\nabla \phi \rangle\). This formula has been well-known in the literature: e.g., formal derivations in [34], [9]. Here the condition \(\nabla \phi \in W^{2,p}\) is assumed for a rigorous justification shown in Lemma 11. Since \(\rho_* \propto e^{-V_{\nu^*}}\), using regularities of \(\rho_*\) and \(V_{\nu^*}\), \[\label{eq-102} \nabla^2 V_{\nu^*} = \nabla (-\rho_* ^{-1} \nabla \rho_* )= \rho_*^{-2} \nabla\rho_* \otimes \nabla\rho_* -\rho_*^{-1} \nabla^2 \rho_* .\tag{45}\] After several applications of integration by parts using 45 , we have \[\label{eq-101}\int (|\nabla^2 \phi|^2 + \nabla^2 V_{\nu^*}[\nabla\phi,\nabla\phi]) \, \rho_* dx= \int |\Delta \phi -\langle \nabla\phi,\nabla V_{\nu^*}\rangle|^2 \,\rho_*dx.\tag{46}\] See Lemma 12 where we derive 46 . Later the authors learned that this observation is already well-known in the literature; e.g., [35].

Let us give a proof of ?? . We aim to prove a stronger statement that the constant \[\label{eq-lambdagap}\lambda_{gap}:=\inf_{\phi \in W^{2,2}, \nabla \phi \neq 0 } \frac{\int (|\nabla^2 \phi|^2 + \nabla^2 V_{\nu^*}[\nabla\phi,\nabla\phi]) \, \rho_* dx}{\int |\nabla \phi|^2 \,\rho_*dx }\tag{47}\] is positive and the infimum is achieved by some \(\phi_0 \in W^{2,2}\). The result ?? follows since \(W^{3,p}\subset W^{2,2}\) for \(p>n\). From its definition and 46 , \(\lambda_{gap}\) is finite non-negative number and there exists a sequence \(\{\phi_i\}_{i=1}^\infty\subset W^{2,2}\) such that the quotient on the right hand side of 47 converges to \(\lambda_{gap}\). For such a sequence we may further assume \(\int |\nabla \phi_i|^2 \rho_*dx =1\) and \(\int \phi_i dx =0\) since replacing \(\phi_i\) by \(\bar\phi_i= \beta_i (\phi_i-\alpha_i)\), for \(\alpha_i=\int \phi_i dx\), \(\beta_i^{-2}= \int |\nabla \phi_i|^2 \rho_*dx\), does not change the quotient.

From the previous observations and the choice \(\phi_i\), the numerator of quotient is bounded uniformly in \(i\). In particular, this implies \(\int |\nabla ^2 \phi_i |^2 dx \le C<\infty\). By weak-compactness Theorem (Banach-Alaoglu) and Rellich-Kondrachov compactness theorem, stated in Theorem 8 and Theorem 10 respectively, there exists a subsequence \(\{\phi_{n_j}\}\subset W^{2,2}\) and \(\phi_0 \in W^{2,2}\) such that \(\phi_{n_j}\) converges to \(\phi_0\) weakly in \(W^{2,2}\) and strongly in \(W^{1,2}\). The strong convergence implies \(\int |\nabla \phi_0|^2 \,\rho_*= \lim_{j\to \infty} \int |\nabla \phi_{n_j}|^2 \, \rho_*=1\). Due to lower-semi-continuity of \(L^2\)-norm under \(L^2\) weak convergence, \[\int (|\nabla^2 \phi_0|^2 + \nabla^2 V_{\nu^*}[\nabla\phi_0,\nabla\phi_0]) \, \rho_* \le \liminf _{j\to \infty} \int (|\nabla^2 \phi_{n_j}|^2 + \nabla^2 V_{\nu^*}[\nabla\phi_{n_j},\nabla\phi_{n_j}]) \, \rho_* =\lambda_{gap}.\] By the definition of \(\lambda_{gap}\), the previous inequality has to be an equality. This shows that \(\lambda_{gap}\ge0\) is achieved by some \(\phi_0\in W^{2,2}\). It remains to prove \(\lambda_{gap}>0\). If \(\lambda_{gap}=0\), 46 implies \[\Delta \phi_0 - \langle \nabla \phi_0 ,\nabla V_{\nu^*}\rangle =0,\] which implies \(\nabla \cdot (e^{-V_{\nu^*}} \nabla \phi_0)=0\). By multiplying with \(\phi_0\) and then integrating by parts, \[\label{eq-energyargument} 0 = \int \phi_0 \nabla \cdot (e^{-V_{\nu^*}} \nabla \phi_0)=-\int e^{-V_{\nu^*}} | \nabla \phi_0 |^2 .\tag{48}\] This contradicts \(\Vert \nabla\phi_0 \Vert_{L^2_{\mu^*}}=1\) and thereby finishes ?? .

For ?? , following similarly as suggested in the beginning, we obtain \(\bar \lambda_{gap}>0\) so that ?? holds with equality for some \(\bar \phi_0\in W^{2,2}\). We choose \(\lambda_{gap}\) as the minimum of two resulting \(\lambda_{gap}\) s obtained. ◻

Proof. Let us prove ?? as ?? similarly follows. Recall we have \(0<C^{-1}\le \frac{d\mu^*}{dx},\frac{d\nu^*}{dx} \le C\) for some \(C=C(f)<\infty\) on \({ \mathbb{T}^n }\). Observe that, for every \(\xi >0\), there exists \(\delta>0\) such that \(\mu\) and \(\nu\) in the statement satisfies \(\Vert \mu-\mu^*\Vert _{W^{1,p}}+\Vert \nu-\nu^*\Vert _{W^{1,p}}\le \xi\) (see Lemma 9). Using Morrey embedding \(W^{1,p} \hookrightarrow C^{0,\alpha}\), after possibly replacing \(\delta>0\) by smaller one, we may get \(\Vert \mu-\mu^*\Vert _{C^0}+\Vert \nu-\nu^*\Vert _{C^0}\le \xi\). Let us prove the statement. Recall 44 (Lemma 11) that, for \(V_\nu(x) = \int _{\mathcal{Y}}f(x,y) \,d\nu(y)\), there holds \[\frac{d^2}{dt^2}\Big \vert_{t=0}F((id+t \nabla \phi)_\# \mu , \nu) = \int |\nabla^2 \phi|^2 + \nabla ^2 V_{\nu}[\nabla \phi ,\nabla \phi] \, d\mu .\] Let us denote \(\frac{d\mu}{dx}=\rho\). For a fixed small \(\xi>0\), \[\begin{align} & \int ( |\nabla^2 \phi|^2 + \nabla ^2 V_{\nu}[\nabla \phi ,\nabla \phi] )\rho = (1-\xi ) \int( |\nabla^2 \phi|^2 + \nabla ^2 V_{\nu^*}[\nabla \phi ,\nabla \phi] )\rho_* \\ &+ \int (\rho -(1-\xi )\rho_*) |\nabla^2 \phi|^2 + (\rho \nabla^2 V_{\nu}-(1-\xi )\rho_*\nabla^2 V_{\nu^*} )[\nabla \phi,\nabla \phi ] \end{align}\] We rearrange the matrices in the last term \[\begin{align} \rho \nabla^2 V_{\nu}-(1-\xi )\rho_*\nabla^2 V_{\nu^*} = \rho (\nabla^2 V_{\nu}-\nabla^2 V_{\nu^*})+ (\rho-\rho_*) \nabla^2 V_{\nu^*} +\xi \rho_* \nabla^2 V_{\nu^*} .\end{align}\] Note that \(\nabla^2 V_{\nu}(x) =\int \nabla^2_xf(x,y)d\nu(y)\) and \(\nabla^2 V_{\nu}\in C^0\). Using \(f\in C^2\), we conclude that there exists \(C=C(f)<\infty\) such that for given \(\xi>0\) there exists \(\delta>0\) small such that, for all \(x\in { \mathbb{T}^n }\), \[\Big \Vert \frac{ \rho(x) \nabla^2 V_{\nu}(x)-(1-\xi )\rho_*(x)\nabla^2 V_{\nu^*}(x) }{\rho_*(x)}\Big \Vert \le C\xi ,\] and \(\rho(x) -(1-\xi )\rho_*(x) \ge 0\). Here \(\Vert \cdot \Vert\) is the operator norm of matrix. Combining above, \[\begin{align} \label{eq-1677}\int ( |\nabla^2 \phi|^2 + \nabla ^2 V_{\nu}[\nabla \phi ,\nabla \phi] )\rho \ge ((1-\xi) \lambda_{gap} - C \xi) \int |\nabla\phi|^2 \rho _* .\end{align}\tag{51}\] For \(\eta>0\), choose \(\xi>0\) and \(\delta>0\) so that \((1-\xi) \lambda_{gap} - C \xi \ge( 1-\frac{1}{2} \eta)\lambda_{gap}\) and \(\Vert 1-\frac{\rho_*}{\rho}\Vert_{L^\infty} \le \frac{1}{2} \eta\). \[\int ( |\nabla^2 \phi|^2 + \nabla ^2 V_{\nu}[\nabla \phi ,\nabla \phi] )\rho \ge (1-\frac{1}{2} \eta)(1-\frac{1}{2} \eta )\lambda_{gap} \int |\nabla \phi |^2 \rho ,\] and the lemma follows as \((1-\frac{1}{2} \eta)^2 \ge (1-\eta)\). ◻

Proof. Fix \(p>n\) and \(\lambda\in(0,\lambda_{gap})\). Choose \(\eta\in(0,1)\) such that \(\lambda \le (1-\eta)\lambda_{gap}\). Let \(\delta=\delta(f,p,\eta)\) be as in Lemma 2 (local convex–concavity).

By Lemma 10, taking \(\varepsilon>0\) small enough we make sure ?? implies: (i) \(\mu_t,\nu_t\) remain uniformly bounded away from \(0\) and \(\infty\) on \({ \mathbb{T}^n }\), and (ii) the optimal maps \(T_{\mu^*\to\mu_t}\) and \(T_{\nu^*\to\nu_t}\) satisfy \[\|T_{\mu^*\to\mu_t}-id\|_{W^{2,p}}+\|T_{\nu^*\to\nu_t}-id\|_{W^{2,p}} \le \delta.\] As a consequence, for every \(s\in[0,1]\), the displacement interpolations (a constant speed geodesic) \[\label{eq-1066} \mu_{t,s}:=\bigl(id+s(T_{\mu^*\to\mu_t}-id)\bigr)_\#\mu^*, \qquad \nu_{t,s}:=\bigl(id+s(T_{\nu^*\to\nu_t}-id)\bigr)_\#\nu^*\tag{52}\] also stay in the \(\delta\)-neighborhood required by Lemma 2.

Fix \(t\in[a,b]\). Let \(T_{\mu_t\to\mu^*}=id+\nabla\tilde{\phi}_t\) be the optimal transport map from \(\mu_t\) to \(\mu^*\). Define a constant speed geodesic from \(\mu_t\) to \(\mu^*\) by \[\tilde{\mu}_{t,s}:=(id+s\nabla\tilde{\phi}_t)_\#\mu_t,\qquad s\in[0,1],\] so that \(\tilde{\mu}_{t,0}=\mu_t\) and \(\tilde{\mu}_{t,1}=\mu^*\). Using \(T_{\mu_t \to \mu^*} = (T_{\mu^*\to \mu_t})^{-1}\) and 52 , we observe that \(\mu_{t,s}=\tilde{\mu}_{t,1-s}\) and similarly \(\nu_{t,s}=\tilde{\nu}_{t,1-s}\). Let \[\ell(s):=F(\tilde{\mu}_{t,s},\nu_t).\] In view of Lemma 11 in Appendix 6.6, it is straightforward to see \(s\mapsto \ell(s)\) is \(C^2\) and thus Taylor’s theorem gives \[\ell(1)=\ell(0)+\ell'(0)+\int_0^1 (1-s)\ell''(s)\,ds.\] A direct first-variation computation (Lemma 11) yields \[\label{eq-taylorexpandell} \ell'(0) =\Big\langle \nabla\Big(\log\frac{d\mu_t}{dx}+V_{\nu_t}\Big),\,\nabla\tilde{\phi}_t\Big\rangle_{L^2_{\mu_t}}, \qquad V_{\nu_t}(x):=\int f(x,y)\,d\nu_t(y).\tag{53}\] Next, an estimate on \(\ell''(s)\) follows by the second variation along geodesic (Lemma 2): \[\ell''(s)\ge (1-\eta)\lambda_{gap} W_2^2(\mu_{t},\mu^*)\ge \lambda W_2^2( \mu_t,\mu^*).\] Explaining this derivation more in detail, for each fixed \(s'\), observe\[id +s \nabla \tilde{\phi}_t =(id +(s-s') \nabla \tilde{\phi}_t\circ (id +s'\nabla \tilde{\phi}_t)^{-1})\circ (id +s'\nabla \tilde{\phi}_t) .\] By a direct calculation, \[\nabla \tilde{\phi}_t\circ (id +s'\nabla \tilde{\phi}_t)^{-1}=\nabla \tilde{\phi}_{t,s'},\] where \(\tilde{\phi}_{t,s'}(y):= \tilde{\phi}_t(x)+\frac{s'}{2} |\nabla \tilde{\phi}_t (x)|^2\), \(x= (id +s'\nabla \tilde{\phi}_t)^{-1}(y)\). Thus Lemma 2 implies \[\ell''(s)\ge (1-\eta)\lambda_{gap} \int |\nabla \tilde{\phi}_{t,s}|^2 d \tilde{\mu}_{t,s}=(1-\eta)\lambda_{gap}\int |\nabla \tilde{\phi}_t|^2 d \mu_t = (1-\eta)\lambda_{gap} W_2^2(\mu_t,\mu^*).\] Combining estimates above, from 53 we infer \[\label{eq:local-AT-mu} F(\mu^*,\nu_t)-F(\mu_t,\nu_t) \ge \Big\langle \nabla\Big(\log\frac{d\mu_t}{dx}+V_{\nu_t}\Big),\,\nabla\tilde{\phi}_t\Big\rangle_{L^2_{\mu_t}} +\frac{\lambda}{2}W_2^2(\mu_t,\mu^*).\tag{54}\]

Under MFL-DA, \(\mu_t\) satisfies the continuity equation with velocity \[u_t=-\nabla\Big(\log\frac{d\mu_t}{dx}+V_{\nu_t}\Big).\] By the first variation formula [40], namely Theorem 15, for a.e.\(t\), \[\frac{d}{dt}\frac{1}{2} W_2^2(\mu_t,\mu^*) = \big\langle -u_t,\,T_{\mu_t\to\mu^*}-id\big\rangle_{L^2_{\mu_t}} = \Big\langle \nabla\Big(\log\frac{d\mu_t}{dx}+V_{\nu_t}\Big),\,\nabla\tilde{\phi}_t\Big\rangle_{L^2_{\mu_t}}.\] Combining with 54 yields the Energy Variational Inequality(EVI) with respect to the equilibrium. \[\label{eq:local-EDI-mu} \frac{d}{dt}\frac{1}{2} W_2^2(\mu_t,\mu^*) \le F(\mu^*,\nu_t)-F(\mu_t,\nu_t)-\frac{\lambda}{2}W_2^2(\mu_t,\mu^*).\tag{55}\]

The same argument (using the local concavity bound in Lemma 2) gives for a.e.\(t\), \[\label{eq:local-EDI-nu} \frac{d}{dt}\frac{1}{2} W_2^2(\nu_t,\nu^*) \le F(\mu_t,\nu_t)-F(\mu_t,\nu^*)-\frac{\lambda}{2}W_2^2(\nu_t,\nu^*).\tag{56}\]

Set \(\mathcal{E}(t):=W_2^2(\mu_t,\mu^*)+W_2^2(\nu_t,\nu^*)\). Summing 55 –56 yields \[\mathcal{E}'(t) \le -\lambda \mathcal{E}(t)+2\big(F(\mu^*,\nu_t)-F(\mu_t,\nu^*)\big).\] Finally, applying the same local above-tangent argument with \((\mu,\nu)\) near \((\mu^*,\nu^*)\) and using stationarity \(\mathop{\mathrm{grad}}_{W_2,\mu}F(\mu^*,\nu^*)=\mathop{\mathrm{grad}}_{W_2,\nu}F(\mu^*,\nu^*)=0\), we get \[F(\mu_t,\nu^*)-F(\mu^*,\nu^*)\ge \frac{\lambda}{2}W_2^2(\mu_t,\mu^*), \qquad F(\mu^*,\nu_t)-F(\mu^*,\nu^*)\le -\frac{\lambda}{2}W_2^2(\nu_t,\nu^*).\] Hence \[F(\mu^*,\nu_t)-F(\mu_t,\nu^*) \le -\frac{\lambda}{2}\mathcal{E}(t).\] Plugging back gives \(\mathcal{E}'(t)\le -2\lambda \mathcal{E}(t)\) for a.e.\(t\in[a,b]\), which is exactly ?? . ◻

Proof. The proof proceeds as in [12]. We rewrite 2 as \[\label{eq-566} \partial_t \mu_t - \Delta \mu_t = \nabla\cdot(\mu_t\nabla B_t), \qquad \partial_t \nu_t - \Delta \nu_t = \nabla\cdot(\nu_t\nabla C_t),\tag{57}\] where \[B_t(x)=\int_{\mathbb{T}^n} f(x,y)\,d\nu_t(y), \qquad C_t(y)=-\int_{\mathbb{T}^n} f(x,y)\,d\mu_t(x).\] Since \(f\in C^2\), we have the uniform bound \[\|B_t\|_{C^2},\|C_t\|_{C^2}\le \|f\|_{C^2}.\]

Uniform \(L^2\) bound. Multiplying 57 by \(\mu_t\) and \(\nu_t\), respectively, and integrating by parts, we obtain \[\partial_t \|\mu_t\|_{L^2}^2 \le -\|\nabla\mu_t\|_{L^2}^2 + C\|\mu_t\|_{L^2}^2, \qquad \partial_t \|\nu_t\|_{L^2}^2 \le -\|\nabla\nu_t\|_{L^2}^2 + C\|\nu_t\|_{L^2}^2,\] for some \(C=C(f)\).

To leverage on the dissipative term, we use the Gagliardo–Nirenberg inequality [45] \[\|\rho\|_{L^2}^{2+\frac{4}{n}} \le C\Bigl( \|\nabla\rho\|_{L^2}^2 \|\rho\|_{L^1}^{\frac{4}{n}} + \|\rho\|_{L^1}^{2+\frac{4}{n}} \Bigr),\] which interpolates between the \(L^1\) mass and the \(H^1\) energy. Since \(\mu_t\) and \(\nu_t\) are probability densities, \(\|\mu_t\|_{L^1}=\|\nu_t\|_{L^1}=1\), this yields \[\|\nabla\mu_t\|_{L^2}^2 \ge \frac{1}{C}\|\mu_t\|_{L^2}^{2(1+\frac{2}{n})}-C,\] and hence \[\partial_t \|\mu_t\|_{L^2}^2 \le -\frac{1}{C}\|\mu_t\|_{L^2}^{2(1+\frac{2}{n})} + C(\|\mu_t\|_{L^2}^2+1).\] Integrating this differential inequality gives \[\label{eq-1626} \|\mu_t\|_{L^2}^2 \le C(t^{-n/2}+1),\tag{58}\] and therefore a uniform \(L^2\) bound for all \(t\ge\delta\). The same estimate holds for \(\nu_t\).

Regularization through bootstrapping. The remaining steps rely on standard regularity theory for linear parabolic equations. We invoke the \(L^p\)-maximal regularity estimate (or Calderón–Zygmund estimate) for the heat equation \(\partial_t\rho_t-\Delta\rho_t=f_t\) [46], [47]: for \(\delta\in(0,1)\), \(k\in\mathbb{Z}\), \(p,q\in(1,\infty)\), and any \(s\ge 0\), \[\label{eq:max-reg-shifted} \int_{s+\delta}^{s+1} \Bigl( \Vert \partial_t \rho_t \Vert _{W^{k,p}}^q + \|\rho_t\|_{W^{k+2,p}}^q \Bigr)\, dt \le C_{\delta,k,p,q } \int_s^{s+1} \Bigl(\|f_t\|_{W^{k,p}}^q + \|\rho_t\|_{W^{k,p}}^q \Bigr) \,dt .\tag{59}\] (Here the left-hand side is measured on the interior time interval \([s+\delta,s+1]\), while the right-hand side uses \([s,s+1]\); the constant \(C_{\delta,k,p,q}\) may blow up as \(\delta\downarrow0\).) Throughout this argument, we repeatedly apply 59 to 57 by viewing \(\nabla \cdot (\mu_t \nabla B_t)\) and \(\nabla \cdot ( \nu_t \nabla C_t)\) as the source term \(f_t\).

First, by 58 , we may view \(\nabla \cdot (\mu_t \nabla B_t) \in L^\infty_t W_x^{-1,2}([s,s+1]) \subset L^{q}_t W^{-1,2}_x([s,s+1])\) for all \(q<\infty\) for \(s>0\).¹² The estimate 59 with \(k=-1\) and \(p=2\) gives uniform \(L^q_tW^{1,2}_x([s+\delta,s+1])\) bounds of \(\mu_t\) and \(\nu_t\). Using this improved estimate, applying 59 once more with \(k=0\) and \(p=2\) yields \(L^q_tW^{2,2}_x([s+\delta,s+1])\) bounds of \(\mu_t\) and \(\nu_t\). Although \(\nabla B_t, \nabla^2 B_t \in L^\infty_t L^\infty_x \subset L^q_t L^\infty_x\) is not sufficiently regular to iterate in \(k\), the Sobolev embedding (Theorem 7) \(W^{2,2}\hookrightarrow W^{1,2^*}\) with \((2^*)^{-1}=2^{-1}-n^{-1}\) allows us to increase the spatial integrability exponent. Iterating in \(p\) while fixing \(k=0\), and applying this argument on intervals of the form \([s,s+1]\), we obtain that for every \(p\in(1,\infty)\) and \(q\in (1,\infty)\), \[\label{eq-1627} \int_{s}^{s+1}\Bigl( \|\partial_t \mu_t\| _{L^p}^q + \|\mu_t\|^q _{W^{2,p}} +\|\partial_t \nu_t\| _{L^p}^q+\|\nu_t\|^q _{W^{2,p}} \Bigr)\, dt \le C(f,\delta,p,q) \quad\text{for all } s\ge \delta .\tag{60}\]

Next, differentiating the system 2 in \(t\), yields \(\dot{\mu}_t = \partial_t \mu_t\) solves \[\label{eq-1629} \partial_t \dot{\mu}_t -\Delta \dot{\mu}_t = \nabla \cdot ( \dot{\mu}_t \nabla B_t+ \mu_t \nabla \dot{B}_t ),\tag{61}\] where \(\dot{B}_t(x)=\int f(x,y)\,d\dot{\nu}_t(y)\). Note \(\nabla \dot{B}_t, \nabla^2 \dot{B}_t \in L^q_t L^\infty_x\) for all \(q\in(1,\infty)\): indeed, for \(|\beta|\le2\), \(\| \nabla^\beta_x \dot{B}_t\|_{L^\infty_x} \le \|\nabla _x^\beta f\|_{L^\infty_{x,y}}\|\dot{\nu}_t\|_{L^1_y} \le C \|\nabla _x^\beta f\|_{L^\infty_{x,y}}\|\dot{\nu}_t\|_{L^p_y}\), and this follows by 60 . Using 60 , we repeat the previous bootstrapping argument for \(\dot{\mu}_t\) and \(\dot{\nu}_t\) on 61 . More precisely, we first apply 59 with \(k=-1\) to 61 to obtain \(\dot{\mu}_t,\dot{\nu}_t\in L^q_tW^{1,p}_x\) for all \(p,q\in(1,\infty)\); this in turn implies \(\nabla\cdot(\dot{\mu}_t\nabla B_t)\) and \(\nabla\cdot(\mu_t\nabla \dot{B}_t)\in L^q_tL^p_x\) for all \(p,q\in(1,\infty)\), so we apply the second time with \(k=0\). As a consequence, \[\label{eq-16272} \int_{s}^{s+1} \Bigl( \|\partial_t \dot{\mu}_t\| _{L^p}^q + \|\dot{\mu}_t\|^q _{W^{2,p}} +\|\partial_t \dot{\nu}_t\| _{L^p}^q+\|\dot{\nu}_t\|^q _{W^{2,p}} \Bigr) \, dt \le C(f,\delta,p,q) \quad\text{for all } s\ge \delta .\tag{62}\]

Finally, \(\mu_t \in L^2_t W^{2,p}_x([s,1+s])\) and \(\dot{\mu}_t \in L^2_t W^{2,p}_x([s,1+s])\) imply \(\mu_t \in C^{0,\frac{1}{2}}_t W^{2,p}_x([s,1+s])\). Namely \(\mu_t \in W^{1,2}_t W^{2,p}_x([s,1+s]) \hookrightarrow C^{0,\frac{1}{2}}_t W^{2,p}_x([s,1+s])\) and this in particular implies (e.g., see [36]) \[\Vert \mu _{t} \Vert_{W^{2,p}} \le C(f,\delta,p), \quad \text{ for all } t\ge \delta .\] The result follows by Morrey’s embedding (Theorem 6) \(W^{2,p}\hookrightarrow C^{1,\alpha}\) for \(p\) sufficiently large. ◻

Proof. The proof follows by compactness argument, as in [12]. Assume the statement is false. Then there exist sequences \(\rho_{0,n},\rho_{1,n}\) with \(\|\rho_{i,n}\|_{C^{1,\alpha}}\le C\) for \(i\in\{0,1\}\) and all \(n\), such that \[\label{eq-tt1} W_2(\rho_{0,n},\rho_{1,n})\le \frac{1}{n},\tag{63}\] while \(\|\rho_{0,n}-\rho_{1,n}\|_{C^{1,\beta}}\ge \varepsilon_0\) for some \(\varepsilon_0>0\). By the compact embedding \(C^{1,\alpha}\hookrightarrow C^{1,\beta}\) (Lemma 7), passing to a subsequence we may assume \(\rho_{0,n}\to\rho_0\) and \(\rho_{1,n}\to\rho_1\) in \(C^{1,\beta}\). Recall that \(W_2\)-distance is controlled by \(L^1\)-distance on compact manifold: Let \(D:=\mathrm{diam}({ \mathbb{T}^n })<\infty\). By Kantorovich–Rubinstein duality, \[W_1(\rho_0,\rho_1) =\sup_{\mathrm{Lip}(\varphi)\le 1}\int_{{ \mathbb{T}^n }}\varphi(\rho_0-\rho_1) \le D\|\rho_0-\rho_1\|_{L^1}.\] Moreover, for any coupling \(\pi\) one has \(\int d^2\,d\pi\le D\int d\,d\pi\), hence \(W_2^2(\rho_0,\rho_1)\le D\,W_1(\rho_0,\rho_1)\). Therefore \(W_2(\rho_0,\rho_1)\le D\,\|\rho_0-\rho_1\|_{L^1}^{1/2}\).

By the triangle inequality and the observation above, \[W_2(\rho_0,\rho_1) \le W_2(\rho_0,\rho_{0,n})+W_2(\rho_{0,n},\rho_{1,n})+W_2(\rho_{1,n},\rho_1) \longrightarrow 0,\] hence \(\rho_0=\rho_1\). Therefore, \[\|\rho_{0,n}-\rho_{1,n}\|_{C^{1,\beta}} \le \|\rho_{0,n}-\rho_0\|_{C^{1,\beta}}+\|\rho_{1,n}-\rho_1\|_{C^{1,\beta}} \longrightarrow 0,\] a contradiction. ◻

Proof of Proposition 4. By the Markov property and working conditionally to \((\mu_{mT}^N,\nu_{mT}^N)\), it is sufficient to prove that there exists \(\eta,C'>0\) such that for any fixed initial condition \((x_1,\dots,x_N,y_1,\dots,y_N)\) with empirical distributions \((\mu_0^N,\nu_0^N)\) satisfying \(r(\mu_0^N,\nu_0^N) \leqslant \delta_0\) (corresponding to the condition that \(\tau >0\)), all \(N \ge 1\) and \(s\in[0,T]\), \[\label{eq:stopped-expect2} \mathbb{E} \left ( r(\mu_{s}^N,\nu_{s}^N) \right) \le C_0 e^{- \lambda_0 s} r(\mu_{0}^N,\nu_{0}^N) + C' N^{-\eta}\,.\tag{72}\] In the rest of the proof, we thus consider such a fixed initial data with \(r(\mu_0^N,\nu_0^N) \leq \delta_0\). For \(s\in[0,T]\), by the triangular inequality, we bound \[\begin{align} r(\mu_s^N,\nu_s^N) & \le r( \mu_s^{\rm mf}, \nu_s^{\rm mf}) + W_2(\mu_s^N,\mu_s^{\rm mf}) + W_2(\nu_s^N, \nu_s^{\rm mf}) \\ & \le C_0 e^{-\lambda_0 s} r(\mu_0^N,\nu_0^N) + \sqrt{2\rm{e}(s)} + W_2(\tilde{\mu}_s^N, \mu_s^{\rm mf}) + W_2(\tilde{\nu}_s^N, \nu_s^{\rm mf}) , \end{align}\] where we used Theorem 1 for the first term (with the fact that \((\mu_0^{\rm mf},\nu_0^{\rm mf})= (\mu_0^N,\nu_0^N)\) by design) and 68 and 69 for the second term. Thanks to [16] (using that a moment of any order is bounded uniformly over all probability measures over the compact torus), there exists \(C,\eta>0\) (which depends only on the dimension \(n\)) such that \[\mathbb{E} \left[W_2(\tilde{\mu}_s^N, \mu_s^{\rm mf})^2 + W_2(\tilde{\nu}_s^N, \nu_s^{\rm mf})^2\right] \leq C N^{-2\eta},\] for all \(s\in[0,T]\). Taking the expectation in 70 , applying the Grönwall Lemma and finally Jensen’s inequality yields \[\begin{align} \mathbb{E} \left[ \sqrt{2\rm{e}(s)} + W_2(\tilde{\mu}_s^N, \mu_s^{\rm mf}) + W_2(\tilde{\nu}_s^N, \nu_s^{\rm mf}) \right] &\le C(T) N^{-\eta} \end{align}\] uniformly over \(s\in[0,T]\) for some constant \(C(T)\), which concludes the proof. ◻

Proof of Lemma 6. Combining 68 and 69 (and similar bound for \(W_2(\nu_t^N,\nu_t^{\rm mf})\)) yields \[\begin{align} \sup_{t\in[0,T]} \left\{ W_2(\mu_t^N,\mu_t^{\rm mf})+ W_2(\nu_t^N,\nu_t^{\rm mf}) \right\} & \le \sup_{t\in[0,T]} \left\{\sqrt{{2\rm e}(t)} + W_2(\tilde{\mu}_t^N,\mu_t^{\rm mf})+ W_2(\tilde{\nu}_t^N,\nu_t^{\rm mf})\right\} \\ & \le C(T) \sup_{t\in[0,T]} \left\{ W_2(\tilde{\mu}_t^N,\mu_t^{\rm mf})+ W_2(\tilde{\nu}_t^N,\nu_t^{\rm mf})\right\}\,, \end{align}\] for some \(C(T)>0\), where in the second line we used 70 , 71 and the Grönwall Lemma. Consequently, Lemma 6 is reduced to proving the same statement with \((\mu_t^N,\nu_t^N)\) replaced by \((\tilde{\mu}_t^N,\tilde{\nu}_t^N)\), namely for i.i.d. (time-inhomogeneous) processes. This result follows by adapting [37]; the theorem does not exactly apply stricto sensu but the proof applies mutatis mutandis. In particular, it is written for processes in \(\mathbb{R}^n\) instead of \(\mathbb{T}^n\), and thus the constants \(b\) depends on \(T\), \(\varepsilon\), the dimension, and \(\|f\|_{C^2}\), since the moments of any order are bounded uniformly over all probability measures over \(\mathbb{T}^n\). ◻