Proof. Let \(\zeta_0\) be a Brenier potential from \(f_0\) to \(f_1\), then by the discussion above \(\zeta_{t}:=(1-t)\frac{|x|^{2}}{2}+t\zeta_0\) is a Brenier potential from \(f_0\) to \(f_{t}\) for any \(t\in [0, 1]\), which belongs to \(\mathcal{C}^{0,\alpha}(\overline{\Omega})\) by Lemma 1 and solves 2 . Re-writing 2 and using that \(\nabla \zeta_t\) is invertible, we have \[f_{t}(z)=\frac{f_0\left((\nabla\zeta_{t})^{-1}(z)\right)}{\det D^{2}\zeta_{t}\left((\nabla\zeta_{t})^{-1}(z)\right)}.\] and again we make use of Lemma 1, which yields the uniform ellipticity bounds \[\Lambda^{-1}\leq D^{2}\zeta_{t}\leq\Lambda,\] where \(\Lambda\) depends on \(\Omega\), \(a\), and \(M\). This proves that \(f_{t}\) has a lower bound depending only on \(a\), \(M\), and \(\Omega\). Finally, we note that \[D\left[(\nabla\zeta_{t})^{-1}\right]=\left(D^{2}\zeta_{t}\right)^{-1}\circ(\nabla\zeta_{t})^{-1}\] is bounded by \(\Lambda\) , so \((\nabla\zeta_{t})^{-1}(z)\) is \(\Lambda\)-Lipschitz and \(f_0\left((\nabla\zeta_{t})^{-1}(z)\right)\) has a \(\mathcal{C}^{0,\alpha}\) norm controlled in terms of \(a\), \(M\), and \(\Omega\). Then the proof is finished after applying Lemmas 11, and 13. ◻

Proof. Write \(T_t(x):=(1-t)x+t\nabla \zeta_0(x)\), recalling that \(\zeta_0\) is a Brenier potential from \(f_0\) to \(f_1\), then under our assumptions, as mentioned above \(T_t\) is an optimal map from \(f_0\) to \(f_t\), smooth on \(\overline{\Omega}\), and is invertible on \(\overline{\Omega}\) for each \(t\in [0, 1]\). For a fixed \(x\in \overline{\Omega}\), define the \(\mathcal{C}^{1, \alpha}\) map \(F_x: [0, 1]\times \overline{\Omega}\to \overline{\Omega}\) by \(F_x(t, y):=T_t(y)-x\). We calculate the Jacobian as \[\begin{align} D_{t, y}F_x(t, y)= \begin{pmatrix} \nabla\zeta_0(y)-y& (1-t)\mathop{\mathrm{Id}}+tD^2\zeta_0(y) \end{pmatrix} \end{align}\] which has full rank for all \(t\in [0, 1]\), thus by the implicit function theorem, \(T_t^{-1}(x)\) is \(\mathcal{C}^1\) in \(t\in (0, 1)\) for any \(x\in \overline{\Omega}\). Differentiating the relation \(F_x(t, T_t^{-1}(x))=0\) in \(t\), we have \[\begin{align} \partial_tT_t^{-1}(x) &=-[(1-t)\mathop{\mathrm{Id}}+tD^2\zeta_0(T_t^{-1}(x))]^{-1}(\nabla \zeta_0(T_t^{-1}(x))-T_t^{-1}(x)) \end{align}\] which we see is continuous and bounded uniformly in \((t, x)\in [0, 1]\times \overline{\Omega}\). Note that by differentiating 2 , first order derivatives of \(\zeta_0\) satisfy a linear, uniformly elliptic (by Lemma 1) equation, with boundary condition that is uniformly oblique by [7]. Thus by the assumed regularity of \(\partial\Omega\) and \(\partial\Omega_*\), we can apply [37] to see \(\zeta_0\in \mathcal{C}^3(\overline{\Omega})\). Since \(f_t(x)=f_0(T_t^{-1}(x))/\det((1-t)\mathop{\mathrm{Id}}+tD^2\zeta_0(T_t^{-1}(x)))\), by differentiating we find for a fixed \(x\in \overline{\Omega}\) \[\begin{align} \partial_tf_t &=f_t\left(\frac{\left\langle \nabla f_0\circ T_t^{-1}, \partial_tT_t^{-1}\right\rangle}{f_0\circ T_t^{-1}}\right.\\ -&\mathop{\mathrm{tr}}\Bigg([(1-t)\mathop{\mathrm{Id}}+t(D^2\zeta_0\circ T_t^{-1})]^{-1}(D^2\zeta_0\circ T_t^{-1})-\mathop{\mathrm{Id}}+\sum_{i=1}^dD^2(\partial_i\zeta_0)\circ T_t^{-1}(\partial_tT_t^{-1})_i\Bigg), \end{align}\] which is bounded uniformly and continuous in \((t, x)\in [0, 1]\times \overline{\Omega}\). Similar calculations hold for \(g_t\).

Thus we may apply Corollary 1, and we have that \(\phi_{t}(x)\) is differentiable with respect to \(t\) for each \(x\) and \(\xi_t\) is a weak solution of \[\begin{cases} {\displaystyle -\relax\left[f_t\left(D^{2}\phi_t\right)^{-1}\nabla\xi_t\right]=\frac{\partial_{t}g_t(\nabla\phi_t)}{g_t(\nabla\phi_t)}f_t-\partial_{t}f_t,} & \textrm{in }\Omega,\\ \left\langle f_t\left(D^{2}\phi_t\right)^{-1}\nabla\xi_{t}, \vec{n}\right\rangle =0, & \textrm{on }\partial\Omega. \end{cases}\] The term \(\partial_{t}f_t\) can be written in divergence form as \(f_{t}\) solves the continuity equation 6 . To show that \(\frac{\partial_{t}g_t(\nabla\phi_t)}{g_t(\nabla\phi_t)}f_t\) can be written as a divergence, fix \(t\in [0, 1]\), consider \(\theta\in\mathcal{C}^{\infty}(\overline{\Omega}_{*})\), and use \(\tilde{\theta}(x):=\theta\left(\nabla\phi_t(x)\right)\) as a test function in 6 to obtain \[\begin{align} \int_{\Omega}\frac{\partial_{t}g_t\left(\nabla\phi_t(x)\right)}{g_t\left(\nabla\phi_t(x)\right)}\tilde{\theta}(x)f_t(x)\,dx & =\int_{\Omega_{*}}\partial_{t}g_t(y)\theta(y)\,dy\\ & =\int_{\Omega_{*}}\left\langle \nabla v_t(y), \nabla\theta(y)\right\rangle g_t(y)\,dy\\ & =\int_{\Omega}\left\langle \nabla v_t\left(\nabla\phi_t(x)\right), \nabla\theta\left(\nabla\phi_t(x)\right)\right\rangle f_t(x)\,dx\\ & =\int_{\Omega}\left\langle \left(D^{2}\phi_t(x)\right)^{-1}\nabla v_t\left(\nabla\phi_t(x)\right), \nabla\tilde{\theta}(x)\right\rangle f_t(x)\,dx. \end{align}\] Since \(x\mapsto \nabla \phi_t(x)\) is a diffeomorphism of \(\overline{\Omega}\) with \(\overline{\Omega}_*\), this shows that in the weak sense, \[\frac{\partial_{t}g_t\left(\nabla\phi_t(x)\right)}{g_t\left(\nabla\phi_t(x)\right)}f_t(x)=-\relax\left[f_{t}\left(D^{2}\phi_{t}\right)^{-1}\nabla v_{t}(\nabla\phi_{t})\right],\] thus \(\xi_t\) is a weak solution of ?? .

Now since \(\Omega\) and \(\Omega_*\) are uniformly convex, by Lemmas 2 and 1, we conclude that the coefficient matrix from ?? is uniformly elliptic, that is \[\Lambda^{-1}\mathop{\mathrm{Id}}\leq\left(D^{2}\phi_{t}\right)^{-1}\leq\Lambda\mathop{\mathrm{Id}},\] where \(\Lambda>0\) only depends on \(\left\Vert f_i\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega})},\left\Vert g_i\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega}_*)}\), \(a\), \(\Omega\), and \(\Omega_{*}\). Since \(f_t\) is also bounded away from zero and infinity on \(\Omega\), by the smoothness of the expression on the right hand side of ?? , we can apply Lemma 6 to obtain that \(\xi_t\in \mathcal{C}^1(\overline{\Omega})\). ◻

Proof. Take \(\theta\in\mathcal{C}^{1}(\overline{\Omega\times\Omega_{*}})\) and let \(\phi_t\) be the Brenier potential from \(f_0\) to \(f_t\) with zero average on \(\Omega\) so that \(\pi_t=(\mathop{\mathrm{Id}}\times \nabla \phi_t)_\sharp f_t\). By Lemma 3 and the discussion within its proof, we have that \(\partial_t(\theta\left(x,\nabla\phi_{t}(x)\right)f_{t}(x))\) is bounded uniformly in \((t, x)\in [0, 1]\times \Omega\), and Clairaut’s theorem applies to yield \(\partial_t\nabla \phi_t=\nabla \xi_t\). Then we may differentiate under the integral below to see \[\begin{align} \frac{d}{dt}\int_{\Omega\times\Omega_{*}}\theta(x,y)\,d\pi_{t}(x,y) & =\frac{d}{dt}\int_{\Omega}\theta\left(x,\nabla\phi_{t}(x)\right)f_{t}(x)\,dx\\ & =\int_{\Omega}\left\langle \nabla\xi_{t}(x), \nabla_{y}\theta\left(x,\nabla\phi_{t}(x)\right)\right\rangle f_{t}(x)\,dx\\ &+\int_{\Omega}\theta\left(x,\nabla\phi_{t}(x)\right)\partial_{t}f_{t}(x)\,dx. \end{align}\] By the regularity of \(\phi_t\) on \(\overline{\Omega}\), we can use the continuity equation 6 for \(f_t\) to obtain \[\begin{align} \int_{\Omega}\theta\left(x,\nabla\phi_{t}(x)\right)\partial_{t}f_{t}(x)\,dx & =\int_{\Omega}\left\langle \nabla_{x}\theta\left(x,\nabla\phi_{t}(x)\right), \nabla u_t(x)\right\rangle f_{t}(x)\,dx \\ & +\int_{\Omega}\left\langle D^{2}\phi_{t}(x)\nabla_{y}\theta\left(x,\nabla\phi_{t}(x)\right), \nabla u_t(x)\right\rangle f_{t}(x)\,dx, \end{align}\] thus combining, we have \[\frac{d}{dt}\int_{\Omega\times\Omega_{*}}\theta\,d\pi_{t}= \int_{\Omega\times\Omega}\left\langle \big(\nabla_{x}\theta,\nabla_{y}\theta\big), \left(\nabla u_t,\nabla\xi_{t}+D^{2}\phi_{t}\nabla u_t\right)\right\rangle d\pi_{t},\] an equivalent weak formulation for the continuity equation by [38]. ◻

Proof of 1. First we prove the estimate when \(f_0\), \(f_1\), \(g_0\), \(g_1\) are \(\mathcal{C}^\infty\), then the result will be extended to Hölder continuous densities by approximation. Note that, due to Lemmas 2 and 1, we have a uniform ellipticity bound on \(D^2\phi_t\). By Lemma 3, for each fixed \(t\) we may use \(\xi_t\) as a test function in ?? to obtain the second line below, followed by Young’s inequality with an \(\varepsilon\) to see, \[\begin{align} \int_{\Omega}\left|\nabla\xi_t\right|^{2}f_t & \leq\Lambda\int_{\Omega}\left\langle \left(D^{2}\phi_t\right)^{-1}\nabla\xi_t, \nabla\xi_t\right\rangle f_t \\ & =\Lambda\int_{\Omega}\left\langle \left( \nabla u_t-\left(D^{2}\phi_t\right)^{-1}\nabla v_t(\nabla\phi_t)\right), \nabla\xi_t\right\rangle f_t\\ & \leq\frac{1}{2}\int_{\Omega}\left|\nabla\xi_t\right|^{2}f_t+\Lambda^{2}\int_{\Omega}\left|\nabla u_t\right|^{2}f_t+\Lambda^{4}\int_{\Omega}\left|\nabla v_t(\nabla\phi_t)\right|^{2}f_t\\ & =\frac{1}{2}\int_{\Omega}\left|\nabla\xi_t\right|^{2}f_t+\Lambda^{2}\int_{\Omega}\left|\nabla u_t\right|^{2}f_t+\Lambda^{4}\int_{\Omega_{*}}\left|\nabla v_t\right|^{2}g_t. \end{align}\] Rearranging and using the Benamou–Brenier formula ([34]) implies \[\label{eq:d295estimate95xi} \int_{0}^{1}\int_{\Omega}\left|\nabla\xi_t\right|^{2}f_tdt\leq C\left[\mathrm{d}_2(f_0,f_1)^2+\mathrm{d}_2(g_0,g_1)^2\right],\tag{8}\] where \(C\) depends on \(\left\Vert f_i\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega})},\left\Vert g_i\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega}_*)}\), \(a\), \(\Omega\), and \(\Omega_{*}\). Since \(\nabla \phi_t(\Omega)\subset \Omega_*\), we have \(\phi_t(x)\) is bounded uniformly for \((t, x)\in [0, 1]\times \overline{\Omega}\), then by Lemma 5 applied to ?? , we see that \(t\mapsto \phi_t(x)\) is Sobolev on \((0, 1)\) for each \(x\), hence absolutely continuous on \([0, 1]\). From here the bound on the gradient is immediate since by Minkowski’s integral inequality, \[\begin{align} \left\Vert \nabla\phi_{1}-\nabla\phi_{0}\right\Vert _{L^{2}(\Omega)} & =\left\Vert \int_{0}^{1}\nabla\xi_t\,dt\right\Vert _{L^{2}(\Omega)}\\ & \leq\int_{0}^{1}\left\Vert \nabla\xi_t\right\Vert _{L^{2}(\Omega)}dt\leq C\left[\mathrm{d}_2(f_0,f_1)+\mathrm{d}_2(g_0,g_1)\right]. \end{align}\]

We now turn to the proof of the estimate ?? . Since by Lemma 4 we have an explicit vector field that drives a curve connecting \(\pi_{0}\) and \(\pi_{1}\), again the Benamou–Brenier formula ([34]) gives us an upper bound on the distance between \(\pi_{0}\), \(\pi_{1}\), yielding \[\begin{align} \mathrm{d}_2(\pi_{0},\pi_{1})^2 & \leq\int_{0}^{1}\int_{\Omega\times\Omega_{*}}\left(|\nabla u_t(x)|^{2}+\left|\nabla\xi_{t}(x)+D^{2}\phi_{t}(x)\nabla u_t(x)\right|^{2}\right)\,d\pi_{t}(x,y)\,dt \\ & \leq\left(2\left\Vert D^{2}\phi_{t}\right\Vert _{L^{\infty}(\Omega\times[0,1])}^{2}+1\right)\int_{0}^{1}\int_{\Omega}|\nabla u_t(x)|^{2}f_{t}(x)\,dx\,dt+2\int_{0}^{1}\int_{\Omega}|\nabla\xi_{t}(x)|^{2}f_{t}(x)\,dx\,dt\\ & =\left(2\left\Vert D^{2}\phi_{t}\right\Vert _{L^{\infty}(\Omega\times[0,1])}^{2}+1\right)\mathrm{d}_2(f_0,f_1)^2+2\int_{0}^{1}\int_{\Omega}|\nabla\xi_{t}(x)|^{2}f_{t}(x)\,dx\,dt. \end{align}\] Now recall that Lemma 2 together with Lemma 1 shows that \[\left\Vert D^{2}\phi_{t}\right\Vert _{L^{\infty}(\Omega\times[0,1])}\leq C\] where \(C\) depends on \(\Omega\), \(\Omega_{*}\), \(a\), \(\left\Vert f_i\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega})}\), and \(\left\Vert g_i\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega}_{*})}\). Then, the proof in this case is finished combining the above inequality with 8 .

Finally, we note that the constants obtained in the estimates ?? and ?? do not depend on any measure of smoothness of the densities other than \(\left\Vert f_i\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega})}\) and \(\left\Vert g_i\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega}_{*})}\). Thus we may approximate \(f_i\) and \(g_i\) by smooth densities in the \(\mathcal{C}^{0,\alpha}\) norm, then taking limits in the resulting estimates finishes the proof: here we can note that by a proof similar to [28], uniform convergence of densities on a bounded domain implies convergence in \(\mathrm{d}_2\), then we use [15] for ?? , and [15] combined with weak lower semicontinuity of \(\mathrm{d}_2\) ([15]) to take the limit in ?? . ◻

Proof of 2. Again the proof will follow from the intuition outlined in [sec:sec:intuition]. In this case, we consider the linear interpolation of the densities \[f_t=(1-t)f_0+tf_1,\qquad g_t=(1-t)g_0+tg_1.\] First assume the densities are \(\mathcal{C}^{1, \alpha}\) regular, then \(t\mapsto f_t\) and \(t\mapsto g_t\) are clearly smooth as curves valued in \(\mathcal{C}^{1, \alpha}(\overline{\Omega})\) and \(\mathcal{C}^{1, \alpha}(\overline{\Omega}_*)\) respectively, so that we can apply Theorem 5. In particular, if \(\tilde{\phi}_t\) is the Brenier potential from \(f_{t}\) to \(g_{t}\) with \(\int_\Omega \tilde{\phi}_t=0\) for all \(t\in [0, 1]\) we have that \(\tilde{\phi}_t\in C^{2, \alpha}(\overline{\Omega})\), is differentiable with respect to \(t\), and \(\tilde{\xi}_t:=\partial_t\tilde{\phi}_t\) solves, in the weak sense (using that \(f_t\) and \(g_t\) are linear interpolations), \[\label{eq:PDE95xi95weaker95metrics} \begin{cases} {\displaystyle-\relax\left[f_t\big(D^{2}\tilde{\phi}_t\big)^{-1}\nabla \tilde{\xi}_t\right]=\frac{g_1(\nabla\tilde{\phi}_t)-g_0(\nabla\tilde{\phi}_t)}{g_{t}(\nabla\tilde{\phi}_t)}f_{t}-(f_1-f_0),} & \textrm{in }\Omega,\\ \left\langle f_t\big(D^{2}\tilde{\phi}_t\big)^{-1}\nabla\tilde{\xi}_t, \vec{n}\right\rangle =0, & \textrm{on }\partial\Omega. \end{cases}\tag{10}\] Now define \(\phi_t:=\tilde{\phi}_t+\log\left(\int_\Omega e^{-\tilde{\phi}_t}\right)\), then we see \(\int_\Omega e^{-\phi_t(x)}\,dx=1\) and \(\nabla\phi_t=\nabla\tilde{\phi}_t\). Since \(\partial_t\tilde{\phi}_t\) is bounded on \(\overline{\Omega}\) and \(\int_\Omega e^{-\tilde{\phi}_t}\neq 0\), we can justify differentiating under the integral to see that \(\phi_t\) is also differentiable in \(t\), and if \(\xi_t:=\partial_t\phi_t\) we have \(\nabla \xi_t=\nabla \tilde{\xi}_t\). Thus 10 holds with \(\xi_t\) and \(\phi_t\) replacing \(\tilde{\xi}_t\) and \(\tilde{\phi}_t\).

Next, since \(\nabla\phi_t(\Omega)=\Omega_{*}\), we see \(\phi_t\) is uniformly Lipschitz independent of \(t\), thus \(\left\Vert \phi_t\right\Vert _{L^{\infty}(\Omega)}\) is bounded in terms of \(\Omega\) and \(\Omega_{*}\); in particular there are \(c_{1}\) and \(c_{2}\) depending on \(\Omega\) and \(\Omega_{*}\) such that \(c_1\leq \rho_t(x):=e^{-\phi_t(x)}\leq c_2\). For the rest of the computation \(C>0\) will denote a constant depending on \(d\), \(\textrm{diam}(\Omega)\), \(\textrm{diam}(\Omega_*)\), \(a\), and \(A\) which may change from line to line. For each \(t\in [0, 1]\) we can apply 9 with \(F=\phi_t\) and \(u=\xi_t\) to obtain \[\begin{align} \int_{\Omega}\left\langle \left(D^{2}\phi_t\right)^{-1}\nabla\xi_t, \nabla\xi_t\right\rangle f_t & \geq C^{-1}\int_{\Omega}\left\langle \left(D^{2}\phi_t\right)^{-1}\nabla\xi_t, \nabla\xi_t\right\rangle \,d\rho_t\\ & \geq C^{-1}\int_{\Omega}\left(\xi_t-\int_\Omega \xi_t d\rho_t\right)^{2}\,d\rho_t\geq C^{-1}\int_{\Omega}\left(\xi_t-\int_\Omega \xi_t d\rho_t\right)^{2}\,dx. \end{align}\] Taking \(\xi_t-\int_\Omega \xi_t d\rho_t\) as a test function in 10 , and combining with the above, we obtain \[\begin{align} &\int_{\Omega}\left(\xi_t-\int_\Omega \xi_t d\rho_t\right)^{2}\,dx \leq C\int_{\Omega}\left\langle \left(D^{2}\phi_t\right)^{-1}\nabla\xi_t, \nabla\xi_t\right\rangle f_t \\ & =-C\int_{\Omega}(f_1-f_0)\left(\xi_t-\int_\Omega \xi_t d\rho_t\right)+C\int_{\Omega}\frac{g_1(\nabla\phi_t)-g_0(\nabla\phi_t)}{g_{t}(\nabla\phi_t)}f_{t}\left(\xi_t-\int_\Omega \xi_t d\rho_t\right) \\ & \leq C\int_{\Omega}(f_1-f_0)^{2}+C\int_{\Omega_{*}}\left(g_1-g_0\right)^{2}+\frac{1}{2}\int_{\Omega}\left(\xi_t-\int_\Omega \xi_t d\rho_t\right)^{2} \end{align}\] which shows \[\label{eqn:32variance32bound} \left\Vert \xi_t-\int_\Omega \xi_t d\rho_t\right\Vert _{L^{2}(\Omega)}\leq C\left(\left\Vert f_1-f_0\right\Vert _{L^{2}(\Omega)}+\left\Vert g_1-g_0\right\Vert _{L^{2}(\Omega_*)}\right).\tag{11}\] For a fixed \(x\in \overline{\Omega}\) we calculate, \[\begin{align} \int_0^1\left(\xi_t(x)-\int_\Omega \xi_t d\rho_t\right)dt &=\phi_1(x)-\phi_0(x)- \int_0^1\int_\Omega\partial_t\phi_t(y)e^{-\phi_t(y)}dydt\\ &=\phi_1(x)-\phi_0(x)+ \int_\Omega\int_0^1\partial_t\left(e^{-\phi_t(y)}\right)dtdy\\ &=\phi_1(x)-\phi_0(x)+ \int_\Omega\left(e^{-\phi_1(y)}-e^{-\phi_0(y)}\right)dy=\phi_1(x)-\phi_0(x). \end{align}\] Now the stability bound is a consequence of Minkowski’s integral inequality combined with 11 : \[\left\Vert \phi_{1}-\phi_{0}\right\Vert _{L^{2}(\Omega)}=\left\Vert \int_{0}^{1}\left(\xi_t-\int_\Omega \xi_t d\rho_t\right)\,dt\right\Vert _{L^{2}(\Omega)}\leq C\left(\left\Vert f_1-f_0\right\Vert _{L^{2}(\Omega)}+\left\Vert g_1-g_0\right\Vert _{L^{2}(\Omega_*)}\right).\] This finishes the proof of 2 for smooth densities.

Now note that the regularity of \(f_i\), \(g_i\) was only used to justify the differentiability of \(\phi_t\) with respect to \(t\) and the bound obtained at the end is independent of such regularity. Therefore, an approximation argument allows us to extend the inequality to discontinuous densities: for general \(f_i, g_i\), \(i=0,1\) as in the hypotheses, let \(\{f_{i, j}\}_{j\in \mathbb{N}}\) and \(\{g_{i, j}\}_{j\in \mathbb{N}}\) be smooth densities with the same upper and lower bounds converging to \(f_i\) and \(g_i\) in \(L^2\), with corresponding Brenier potentials \(\{\phi_{i, j}\}_{j\in \mathbb{N}}\) which have average zero. Since \(\nabla \phi_{i, j}(\Omega)\subset \Omega_*\), we can use Arzelà–Ascoli to obtain subsequences that converge uniformly on \(\Omega\) to \(\phi_i\), necessarily convex; note the Legendre transforms of \(\phi_{i, j}\) also uniformly converge to the Legendre transform of \(\phi_i\). In particular, we can take a limit in the Kantorovich dual problem (see [15]) to see \(\phi_i\) is also a Brenier potential from \(f_i\) to \(g_i\). Thus we can take a limit in the corresponding estimates for \(\phi_{i, j}\) to finish the proof of ?? .

Finally, ?? follows immediately from ?? combined with [22]. ◻

Proof of Theorem 3. As before, we consider \[f_{t}=(1-t)f_0+tf_1,\qquad g_{t}=(1-t)g_0+tg_1\] and interpolate between \(\phi_0\), \(\phi_1\) by letting \(\phi_t\) be the Brenier potential from \(f_{t}\) to \(g_{t}\) with \(\int_\Omega \phi_t=0\) for all \(t\in [0, 1]\). Then Theorem 5 shows that \(\xi_{t}:=\partial_{t}\phi_{t}\) is a weak solution of \[\begin{cases} {\displaystyle -\relax\left(f_{t}\left(D^{2}\phi_{t}\right)^{-1}\nabla\xi_{t}\right)=\frac{g_1(\nabla\phi_{t})-g_0(\nabla\phi_{t})}{g_{t}(\nabla\phi_{t})}f_{t}-(f_1-f_0),} & \textrm{in }\Omega,\\ \left\langle f_t \left(D^{2}\phi_{t}\right)^{-1}\nabla\xi_t,\vec{n}\right\rangle =0, & \textrm{on }\partial\Omega. \end{cases}\] As done previously, we use Lemma 1 to obtain a uniform ellipticity estimate \[C^{-1} \mathop{\mathrm{Id}}\leq f_t(D^2\phi_t)^{-1} \leq C \mathop{\mathrm{Id}},\] then combining this with Lemma 13 shows \[\Vert f_{t}(D^{2}\phi_{t})^{-1}\Vert_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq C.\] Therefore, we can apply Lemma 6 followed by Lemma 5 to obtain \[\begin{align} \left\Vert \xi_t\right\Vert _{\mathcal{C}^{1,\alpha}(\overline{\Omega})} & \leq C\left(\left\Vert \xi_t\right\Vert _{L^{\infty}(\Omega)}+\left\Vert \frac{g_1(\nabla\phi_{t})-g_0(\nabla\phi_{t})}{g_{t}(\nabla\phi_{t})}f_{t}-(f_1-f_0)\right\Vert _{L^{\frac{d}{1-\alpha}}(\Omega)}\right) \\ &\leq C\left\Vert \frac{g_1(\nabla\phi_{t})-g_0(\nabla\phi_{t})}{g_{t}(\nabla\phi_{t})}f_{t}-(f_1-f_0)\right\Vert _{L^{\frac{d}{1-\alpha}}(\Omega)} \\ & \le C\left(\left\Vert f_1-f_0\right\Vert _{L^{\frac{d}{1-\alpha}}(\Omega)}+\left\Vert g_1-g_0\right\Vert _{L^{\frac{d}{1-\alpha}}(\Omega_{*})}\right). \end{align}\] By the above bound, we can justify exchanging integration and differentiation below to obtain \[\begin{align} \left\Vert \phi_{1}-\phi_{0}\right\Vert _{\mathcal{C}^{1,\alpha}(\overline{\Omega})} &=\left\Vert \int_{0}^{1}\xi_{t}\,dt\right\Vert _{\mathcal{C}^{1,\alpha}(\overline{\Omega})} \leq \int_{0}^{1}\left\lVert \xi_{t}\right\rVert _{\mathcal{C}^{1,\alpha}(\overline{\Omega})}\,dt\\ &\leq C\left(\left\Vert f_1-f_0\right\Vert _{L^{\frac{d}{1-\alpha}}(\Omega)}+\left\Vert g_1-g_0\right\Vert _{L^{\frac{d}{1-\alpha}}(\Omega_{*})}\right). \qedhere \end{align}\] ◻

Proof of Remark 1. Let us first show the one dimensional bound \[\left\Vert \phi_{1}'-\phi_{0}'\right\Vert _{L^{\infty}(I)}\leq\frac{1}{\varepsilon}\left(\left\Vert f_1-f_0\right\Vert _{L^{1}(I)}+\left\Vert g_1-g_0\right\Vert _{L^{1}(J)}\right),\] which we stated in (?? ). We use that in one dimension the optimal transport map can be simply expressed as \[\phi_{i}'=G_i^{-1}\circ F_i,\] where \(F_i\), \(G_i\) are the cumulative distribution functions of \(f_i\) and \(g_i\). From here the bound is an elementary computation. Indeed, for \(x_1\leq x_2\) in \(J\), \[\begin{align} \lvert G_i(x_1)-G_i(x_2)\rvert &=\int_{x_1}^{x_2}g_i\geq a\lvert x_2-x_1\rvert, \end{align}\] hence the \(G_i^{-1}\) are \(\frac{1}{a}\)-Lipschitz and we have \[\begin{align} \left|G_1^{-1}\left(F_1(x)\right)-G_0^{-1}\left(F_0(x)\right)\right| & \leq\left|G_1^{-1}\left(F_1(x)\right)-G_1^{-1}\left(F_0(x)\right)\right|+\left|G_1^{-1}\left(F_0(x)\right)-G_0^{-1}\left(F_0(x)\right)\right|\\ & \le\frac{1}{a}\left\Vert F_1-F_0\right\Vert _{L^{\infty}(I)}+\left\Vert G_1^{-1}-G_0^{-1}\right\Vert _{L^{\infty}([0,1])}\\ & \leq\frac{1}{a}\left\Vert f_1-f_0\right\Vert _{L^{1}(I)}+\left\Vert G_1^{-1}-G_0^{-1}\right\Vert _{L^{\infty}([0,1])}. \end{align}\] For the second term, we note that given \(y\in[0,1]\) \[\left|G_1^{-1}(y)-G_0^{-1}(y)\right|=\left|G_0^{-1}\left(G_0\left(G_1^{-1}(y)\right)\right)-G_0^{-1}(y)\right|\leq\frac{1}{a}\left|G_0\left(G_1^{-1}(y)\right)-y\right|.\] Therefore \[\begin{align} \left\Vert G_1^{-1}-G_0^{-1}\right\Vert _{L^{\infty}([0,1])} & \leq\frac{1}{a}\sup_{y\in[0,1]}\left|G_0\left(G_1^{-1}(y)\right)-y\right|\\ & =\frac{1}{a}\sup_{x\in J}\left|G_0(x)-G_1(x)\right|\leq\frac{1}{a}\left\Vert g_1-g_0\right\Vert _{L^{1}(J)}, \end{align}\] which proves (?? ).

Next, we show that, if we remove the lower bound hypothesis \(g_i\geq a\) from Theorem 3, there exists no pair of constants \(0<\eta <1\) and \(C>0\) such that \[\left\Vert \nabla\phi_{1}-\nabla\phi_{0}\right\Vert _{L^{\infty}(\Omega)}\leq C\left(\left\Vert f_0-f_1\right\Vert _{L^{\infty}(\Omega)}^\eta+\left\Vert g_0-g_1\right\Vert _{L^{\infty}(\Omega_{*})}^\eta\right).\] We prove this with a one dimensional example on \(\Omega = \Omega_*=(0,1)\). We set \(f\equiv 1\) and \[g_{a}(x)= \begin{cases} (p+1)a^{-p}(a-x)^{p}, & 0\leq x\leq a\\ (p+1)(1-a)^{-p}(x-a)^{p}, & a\leq x\leq1 \end{cases}, \qquad 0<a<1.\] In this case, the optimal transport map between \(f\) and \(g_a\) is given by \(G_a^{-1}\), where \(G_a\) is the cumulative function of \(g_a\). We will consider \(a\) in a small neighborhood of \(\frac{1}{2}\) so that, for each fixed \(p>0\), the functions \(g_a\) have a uniform \(\mathcal{C}^{0,\min(p,1)}(\Omega)\) bound. A simple computation shows that \[G_{a}^{-1}(y)= \begin{cases} a\left(1-\left(1-\frac{y}{a}\right)^{\frac{1}{p+1}}\right), & 0\leq y\leq a,\\ a+(1-a)\left(\frac{y-a}{1-a}\right)^{\frac{1}{p+1}}, & a\leq y\leq1. \end{cases}\] Since the \(g_a\) are defined so that \(G_a(a)=a\), setting \(a=\frac{1}{2(1-\varepsilon)}=\frac{1}{2}+\frac{\varepsilon}{2}+O(\varepsilon^{2})\) we have \[\left\Vert G_{a}^{-1}-G_{1/2}^{-1}\right\Vert _{L^{\infty}([0,1])}\geq\left|G_{a}^{-1}(1/2)-G_{1/2}^{-1}(1/2)\right|=\left|\frac{1-\varepsilon^{\frac{1}{p+1}}}{2(1-\varepsilon)}-\frac{1}{2}\right|\gtrsim\varepsilon^{\frac{1}{p+1}}.\] On the other hand, a simple computation shows that for each fixed \(p>1\) it holds \[\left\Vert g_{a}-g_{1/2}\right\Vert _{L^{\infty}([0,1])}\lesssim \varepsilon.\] Therefore, for any \(\frac{1}{p+1}<\eta<1\) we have \[\limsup_{\varepsilon\to 0}\frac{\left\Vert G_{a}^{-1}-G_{1/2}^{-1}\right\Vert _{L^{\infty}([0,1])}}{\left\Vert g_{a}-g_{1/2}\right\Vert _{L^{\infty}([0,1])}^{\eta}}=+\infty.\qedhere\] ◻

Proof. Recall by Kantorovich duality (see [15]) it holds for each \(t\), \[\frac{1}{2}\mathrm{d}_2(f_t, g_t)^2=\sup\left\{ \int_\Omega\left(\frac{\lvert x\rvert^2}{2}-\phi(x)\right)f_t(x)\,dx+\int_{\Omega_*}\left(\frac{\lvert y\rvert^2}{2}-\psi(y)\right)g_t(y)\,dy\right\},\] where the supremum is taken over \(\phi(x)\), \(\psi(y)\) satisfying \(\phi(x)+\psi(y)\leq \left\langle x, y\right\rangle\). Then the pair \((\phi_t, \phi_t^*)\) is a maximizer in the problem above, where \(\phi_t^*\) is the Legendre transform of \(\phi_t\) defined by \[\begin{align} \phi_t^*(y):=\sup_{x\in \Omega}(\left\langle x, y\right\rangle-\phi_t(x)). \end{align}\] Moreover, if \(T_t=\nabla\phi_t(x)\) is the optimal map from \(f_t\) to \(g_t\), then under our conditions \(T_t\) is invertible and it holds for all \(y\in \Omega_*\) that \[\begin{align} \label{eqn:32dual32potential32equality} \phi_t^*(y)=\left\langle T_t^{-1}(y), y\right\rangle-\phi_t(T_t^{-1}(y)). \end{align}\tag{12}\] Regarding time differentiability, Corollary 1 shows that \(\nabla\phi_t\) is \(\mathcal{C}^1\) with respect to \(t\), and using the same argument as in the beginning of the proof of Lemma 3, we see that \(T_t^{-1}(y)\) is differentiable in \(t\) for any \(y\in \Omega_*\) and \(\lvert t\rvert\) small enough, with \(\partial_tT_t^{-1}(y)\) uniformly bounded for such \((t, y)\). Thus 12 yields that \(\phi_t^*(y)\) is differentiable in \(t\) for \(y\in \Omega_*\), with uniformly bounded derivative, and we can justify differentiating under the integral to obtain \[\begin{align} 0 &=\partial_h\vert_{h=0}\left[\int_\Omega\left(\frac{\lvert x\rvert^2}{2}-\phi_{t+h}(x)\right)f_t(x)\,dx+\int_{\Omega_*}\left(\frac{\lvert y\rvert^2}{2}-\phi_{t+h}^*(y)\right)g_t(y)\,dy\right]\\ &=-\int_\Omega\partial_t\phi_t(x)f_t(x)\,dx-\int_{\Omega_*}\partial_t\phi_t^*(y)g_t(y)\,dy, \end{align}\] where we have used that quantity being differentiated above achieves its maximum at \(h=0\). Then we calculate, \[\begin{align} \frac{1}{2}\frac{d}{dt}\mathrm{d}_2(f_{t},g_{t})^2 & =\frac{d}{dt}\int_{\Omega}\left(\frac{|x|^2}{2}-\phi_{t}(x)\right)f_{t}(x)\,dx +\frac{d}{dt}\int_{\Omega_{*}}\left(\frac{|y|^2}{2}-\phi_{t}^*(y)\right)g_{t}(y)\,dy\\ &=\int_{\Omega}\left(\frac{|x|^2}{2}-\phi_{t}(x)\right)\partial_{t}f_{t}(x)\,dx +\int_{\Omega_{*}}\left(\frac{|y|^2}{2}-\phi_{t}^*(y)\right)\partial_{t}g_{t}(y)\,dy\\ & -\int_{\Omega}\partial_{t}\phi_{t}(x)f_{t}(x)\,dx-\int_{\Omega_{*}}\partial_{t}\phi_{t}^*(y)g_{t}(y)\,dy\\ & =\int_{\Omega}\left(\frac{|x|^2}{2}-\phi_{t}(x)\right)h(x)f(x)\,dx +\int_{\Omega_{*}}\left(\frac{|y|^2}{2}-\phi_{t}^*(y)\right)k(y)g(y)\,dy. \end{align}\] We can see the last expression above is also differentiable in \(t\), proving the claim. ◻

Proof of 4. First we assume \(f\), \(g\), \(h\), \(k\) are all \(\mathcal{C}^\infty\). Note that for \(|t|<\delta\) small enough we have that \(f_t = (1+th)f\), \(g_t = (1+tk)g\) are bounded away from zero, so we can apply Theorem 5 and conclude that for any \(0<\beta <\alpha\) the mapping \(t\mapsto \phi_t\in \mathcal{C}^{2,\beta}(\overline{\Omega})\) is \(\mathcal{C}^2\). This also implies that \(\phi\in \mathcal{C}^{2}\big((-\delta,\delta)\times \overline{\Omega}\big)\) as a function of \(t\), \(x\) and that \[t\mapsto \nabla \phi_t \in \mathcal{C}^{1,\beta}(\overline{\Omega}), \qquad t\mapsto D^2\phi_t\in \mathcal{C}^{0,\beta}(\overline{\Omega};\mathbb{R}^{d\times d})\] are twice differentiable as curves on Banach spaces. This justifies the application of the dominated convergence theorem in the subsequent computations.

To obtain the exact expression in ?? , we compute the discrete difference \[\begin{align} \frac{\mathrm{d}_2(f_{t},g_{t})^2+\mathrm{d}_2(f_{-t},g_{-t})^2-2\mathrm{d}_2(f,g)^2}{t^{2}} & =\frac{1}{t^{2}}\int_{\Omega}\left|\nabla\phi_{t}-x\right|^{2}f_{t}\,dx \\ & +\frac{1}{t^{2}}\int_{\Omega}\left|\nabla\phi_{-t}-x\right|^{2}f_{-t}\,dx-\frac{2}{t^{2}}\int_{\Omega}\left|\nabla\phi-x\right|^{2}f\,dx\\ & =\int_{\Omega}\underbrace{\frac{\left|\nabla\phi_{t}-x\right|^{2}+\left|\nabla\phi_{-t}-x\right|^{2}-2\left|\nabla\phi-x\right|^{2}}{t^{2}}}_{I}\,f\,dx\\ & +\int_{\Omega}\frac{\left|\nabla\phi_{t}-x\right|^{2}-\left|\nabla\phi_{-t}-x\right|^{2}}{t}hf\,dx. \end{align}\] After some algebraic manipulations we have \[I=\left\langle \frac{(\nabla\phi_{t}-\nabla \phi)+(\nabla \phi-\nabla\phi_{-t})}{t}, \frac{\nabla\phi_{t}-\nabla\phi}{t}\right\rangle+\left\langle \nabla\phi_{-t}+\nabla\phi-2x, \frac{\nabla\phi_{t}+\nabla\phi_{-t}-2\nabla\phi}{t^{2}}\right\rangle,\] and we similarly obtain \[\frac{\left|\nabla\phi_{t}-x\right|^{2}-\left|\nabla\phi_{-t}-x\right|^{2}}{t} =\left\langle \nabla\phi_{t}+\nabla\phi_{-t}-2x, \frac{(\nabla\phi_{t}-\nabla \phi)+(\nabla \phi-\nabla\phi_{-t})}{t}\right\rangle.\] Therefore, letting \(t\to 0\) and denoting \(\xi=\partial_{t}|_{t=0}\phi_t\), we deduce \[\begin{align} \nonumber \left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\mathrm{d}_2(f_{t},g_{t})^2 & = 4\int_{\Omega}\left\langle \nabla\phi(x)-x, \nabla\xi(x)\right\rangle h(x)f(x)\,dx + 2\int_{\Omega}\left|\nabla\xi\left(x\right)\right|^{2}f(x)\,dx \\ & +2\int_{\Omega}\left\langle \nabla\phi(x)-x, \left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\nabla\phi_{t}(x)\right\rangle f(x)\,dx \label{eq:missing95term952var95OT} \end{align}\tag{13}\] Next, we compute \[\int_{\Omega}\left\langle \nabla\eta\left(\nabla\phi(x)\right), \left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\nabla\phi_{t}(x)\right\rangle f(x)\,dx,\] where here \(\eta\in \mathcal{C}^2(\overline{\Omega}_*)\) is a general test function. To do this we note that \[\begin{align} 0 & =\int_{\Omega_*}\frac{g_{t}(y)+g_{-t}(y)-2g(y)}{t^{2}}\eta(y)\,dy \\ & =\frac{1}{t^{2}}\int_{\Omega}\eta\left(\nabla\phi_{t}(x)\right)f_{t}(x)\,dx +\frac{1}{t^{2}}\int_{\Omega}\eta\left(\nabla\phi_{-t}(x)\right)f_{-t}(x)\,dx -\frac{2}{t^{2}}\int_{\Omega}\eta\left(\nabla\phi(x)\right)f(x)\,dx\\ & =\int_{\Omega}\frac{\eta\left(\nabla\phi_{t}(x)\right)+\eta\left(\nabla\phi_{-t}(x)\right)-2\eta\left(\nabla\phi(x)\right)}{t^{2}}f(x)\,dx \\ &\qquad\qquad+\int_{\Omega}\frac{\big(\eta\left(\nabla\phi_{t}(x)\right)-\eta\left(\nabla\phi(x)\right)\big) +\big(\eta\left(\nabla\phi(x)\right)-\eta\left(\nabla\phi_{-t}(x)\right)\big)}{t}h(x)f(x)\,dx. \end{align}\] Letting \(t\to0\) yields \[\begin{align} 0 & =\int_{\Omega}\left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\left[\eta\left(\nabla\phi_{t}(x)\right)\right]f(x)\,dx +2\int_{\Omega}\left.\frac{d}{dt}\right|_{t=0}\left[\eta\left(\nabla\phi_{t}(x)\right)\right]h(x)f(x)\,dx\\ & =\int_{\Omega}\left\langle D^{2}\eta\left(\nabla\phi(x)\right)\nabla\xi\left(x\right),\nabla\xi\left(x\right)\right\rangle f(x)\,dx +\int_{\Omega}\left\langle \nabla\eta\left(\nabla\phi(x)\right), \left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\nabla\phi_{t}(x)\right\rangle f(x)\,dx\\ & \,+2\int_{\Omega}\left\langle \nabla\eta\left(\nabla\phi(x)\right), \nabla\xi\left(x\right)\right\rangle h(x)f(x)\,dx, \end{align}\] which shows that \[\begin{align} \int_{\Omega}\left\langle \nabla\eta\left(\nabla\phi(x)\right), \left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\nabla\phi_{t}(x)\right\rangle f(x)\,dx & =-\int_{\Omega}\left\langle D^{2}\eta\left(\nabla\phi(x)\right)\nabla\xi\left(x\right),\nabla\xi\left(x\right)\right\rangle f(x)\,dx \\ & \,-2\int_{\Omega}\left\langle \nabla\eta\left(\nabla\phi(x)\right), \nabla\xi\left(x\right)\right\rangle h(x)f(x)\,dx. \end{align}\] Now we can set \({\eta(y)= \frac{1}{2}|y|^2-\psi(y)}\) where \(\psi\) is a conjugate Brenier potential from \(g\) to \(f\) satisfying \({\nabla \psi = (\nabla \phi)^{-1}}\) to obtain \[\begin{align} \int_{\Omega}\left\langle \nabla\phi(x)-x, \left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\nabla\phi_{t}(x)\right\rangle f(x)\,dx & =-\int_{\Omega}\left\langle \left(\mathop{\mathrm{Id}}-D^{2}\psi\left(\nabla\phi(x)\right)\right)\nabla\xi(x),\nabla\xi(x)\right\rangle f(x)\,dx \\ & \,-2\int_{\Omega}\left\langle \nabla\phi(x)-\nabla\psi\left(\nabla\phi(x)\right), \nabla\xi(x)\right\rangle h(x)f(x)\,dx\\ & =-\int_{\Omega}\left\langle \left(\mathop{\mathrm{Id}}-\big(D^{2}\phi(x)\big)^{-1}\right)\nabla\xi(x),\nabla\xi(x)\right\rangle f(x)\,dx\\ & \,-2\int_{\Omega}\left\langle \nabla\phi(x)-x, \nabla\xi(x)\right\rangle h(x)f(x)\,dx. \end{align}\] Plugging this into (13 ) yields \[\left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\mathrm{d}_2(f_{t},g_{t})^2 =2\int_{\Omega}\left\langle \left(D^{2}\phi(x)\right)^{-1}\nabla\xi(x), \nabla\xi(x)\right\rangle f(x)\,dx,\] which proves Theorem 4 for smooth densities.

Next, we prove Theorem 4 assuming only \(f\), \(h\in \mathcal{C}^{0,\alpha}(\overline{\Omega})\), \(g\), \(k\in \mathcal{C}^{0,\alpha}(\overline{\Omega}_*)\) through an approximation argument. There exist linear and bounded extension operators \[\mathcal{E}:\mathcal{C}^{0,\alpha}(\overline{\Omega})\longrightarrow C^{0,\alpha}(\mathbb{R}^d), \qquad \mathcal{E}_*:\mathcal{C}^{0,\alpha}(\overline{\Omega}_*)\longrightarrow C^{0,\alpha}(\mathbb{R}^d),\] constructed as in \(\mathcal{E}_0\) in [45]. Therefore, without loss of generality, we may assume that \(f,g,h,k\in\mathcal{C}^{0,\alpha}(\mathbb{R}^{d})\) with bounds of the form \({\left\Vert f\right\Vert _{\mathcal{C}^{0,\alpha}(\mathbb{R}^{d})}\leq c_d\left\Vert f\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega})}}\). Then given an approximation kernel \(m^\varepsilon\in \mathcal{C}^\infty_c(\mathbb{R}^d)\), we set \[f^{\varepsilon}=c_{\varepsilon,f}(f*m^{\varepsilon}),\quad g^{\varepsilon}=c_{\varepsilon,g}(g*m^{\varepsilon}),\quad h^{\varepsilon}=h*m^{\varepsilon}-c_{\varepsilon,h},\quad k^{\varepsilon}=k*m^{\varepsilon}-c_{\varepsilon,k},\] where the constants are chosen so that \(f^{\varepsilon},g^{\varepsilon}\) are probability densities on \(\Omega,\Omega_{*}\) and \[\int_{\Omega}h^{\varepsilon}f^{\varepsilon}=\int_{\Omega_{*}}k^{\varepsilon}g^{\varepsilon}=0.\] Choosing \(\varepsilon>0\) small enough it is easy to see that \(f^{\varepsilon},g^{\varepsilon}\) will be bounded from below by \(a/2\), in which case we can define the curve of densities \[f_{t}^{\varepsilon}:=(1+th^{\varepsilon})f^{\varepsilon},\qquad g_{t}^{\varepsilon}:=(1+tk^{\varepsilon})g^{\varepsilon}\] for \(\left|t\right|<\delta\) sufficiently small. Let \(\phi^\varepsilon_t\) be the Brenier potential from \(f^\varepsilon_t\) to \(g^\varepsilon_t\) with average zero and \(\xi^\varepsilon_t:=\partial_t\phi^\varepsilon_t\). Since the densities are smooth, we can apply the previous computations to \(\phi_t\) replaced by \(\phi^\varepsilon_{t+u}\) for \(\lvert u\rvert\) sufficiently small to obtain \[\left.\frac{d^{2}}{dt^{2}}\right|_{t=u}\mathrm{d}_2(f_{t}^{\varepsilon},g_{t}^{\varepsilon})^2=2\int_{\Omega}\left\langle \left(D^{2}\phi_{u}^{\varepsilon}\right)^{-1}\nabla\xi_{u}^{\varepsilon},\nabla\xi_{u}^{\varepsilon}\right\rangle f_{u}^{\varepsilon}.\] Now we let \(\varepsilon\to 0\) after re-writing the previous expression in integral form. In general, if \(\varphi:[-t,t]\longrightarrow\mathbb{R}\) is a \(\mathcal{C}^{2}\) function, the second order quotient can be re-written as \[\frac{\varphi(t)+\varphi(-t)-2\varphi(0)}{t^{2}} =2\int_{0}^{1}\int_{0}^{1}s\varphi''\left((2\tau-1)st\right)\,d\tau\,ds.\] Denoting \(u_{t}(\tau,s):=(2\tau-1)st\in[-t,t]\), we see \(\lvert u_t(\tau, s)\rvert\) can be made small uniformly in \(\tau\), \(s\in [0, 1]\) by taking \(\lvert t\rvert\) small, thus applying the above integral identity to \(t\mapsto \mathrm{d}_2(f_{t}^{\varepsilon},g_{t}^{\varepsilon})^2\) we obtain \[\begin{align} &\frac{\mathrm{d}_2(f_{t}^{\varepsilon},g_{t}^{\varepsilon})^2+\mathrm{d}_2(f_{-t}^{\varepsilon},g_{-t}^{\varepsilon})^2-2\mathrm{d}_2(f_0^{\varepsilon},g_0^{\varepsilon})^2}{t^{2}}\\ &=4\int_{0}^{1}\int_{0}^{1}\int_{\Omega}s\left\langle \left(D^{2}\phi_{u_{t}(\tau,s)}^{\varepsilon}\right)^{-1}\nabla\xi_{u_{t}(\tau,s)}^{\varepsilon},\nabla\xi_{u_{t}(\tau,s)}^{\varepsilon}\right\rangle f_{u_{t}(\tau,s)}^{\varepsilon}\,d\tau\,ds. \end{align}\] In order to pass to the limit we obtain estimates uniform in \(\varepsilon\) and \(t\). First, Lemma 1 implies that for \(\beta<\alpha\), \[\left\Vert \phi_{t}^{\varepsilon}\right\Vert _{\mathcal{C}^{2,\beta}(\overline{\Omega})}\leq C,\] where \(C\) depends on \(\left\Vert f_{t}^{\varepsilon}\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\), \(\left\Vert g_{t}^{\varepsilon}\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega}_{*})}\) , and the lower bounds of \(f^{\varepsilon}\) and \(g^{\varepsilon}\). Restricting the time interval if necessary and choosing \(\varepsilon\) small enough, all these quantities are controlled by the bounds of \(f_{t}\) and \(g_{t}\). Consequently, we obtain \[\left\Vert \phi_{t}^{\varepsilon}\right\Vert _{\mathcal{C}^{2,\beta}(\overline{\Omega})}\leq C, \qquad \Lambda^{-1}\mathop{\mathrm{Id}}\leq D^{2}\phi_{t}^{\varepsilon}\leq\Lambda\mathop{\mathrm{Id}}\] with constants independent of \(t\), \(\varepsilon\). Then classical elliptic regularity estimates (see Lemmas 5 and 6) give uniform bounds on \(\Vert \xi_t^\varepsilon\Vert_{\mathcal{C}^{1,\beta}(\overline{\Omega})}\) which allows us to take \(\varepsilon\to 0\) in the weak formulation of the PDE satisfied by \(\xi_t^\varepsilon\) using the Arzelà–Ascoli Theorem. In summary, as \(\varepsilon\to0\) we have the following convergence, uniformly for \(\lvert u\rvert\) small: \[\left(D^{2}\phi_{u}^{\varepsilon}\right)^{-1}\to\left(D^{2}\phi_{u}\right)^{-1},\qquad f_{u}^{\varepsilon}\to f_{u},\qquad\nabla\xi_{u}^{\varepsilon}\to\nabla\xi_{u}.\] Therefore, letting \(\varepsilon\to0\) and using the dominated convergence theorem gives \[\begin{align} &\frac{\mathrm{d}_2(f_{t},g_{t})^2+\mathrm{d}_2(f_{-t},g_{-t})^2-2\mathrm{d}_2(f_0,g_0)^2}{t^{2}} \\ &= 4\int_{0}^{1}\int_{0}^{1}\int_{\Omega}s\left\langle \left(D^{2}\phi_{u_{t}(\tau,s)}\right)^{-1}\nabla\xi_{u_{t}(\tau,s)},\nabla\xi_{u_{t}(\tau,s)}\right\rangle f_{u_{t}(\tau,s)}\,d\tau\,ds. \end{align}\] Given the uniform bounds on \(\left\Vert \phi_{t}\right\Vert _{\mathcal{C}^{2,\beta}(\overline{\Omega})}\), \(\left\Vert \xi_{t}\right\Vert _{\mathcal{C}^{1,\beta}(\overline{\Omega})}\), a simple compactness argument shows that \(D^{2}\phi_{t}\), \(\nabla\xi_{t}\) are continuous in \(t\), so we can let \(t\to0\) and the dominated convergence theorem gives \[\begin{align} \left.\frac{d^{2}}{dt^{2}}\right|_{t=0}\left[\mathrm{d}_2(f_{t},g_{t})^2\right] =4\int_0^1sds \int_{\Omega}\left\langle \left(D^{2}\phi\right)^{-1}\nabla\xi,\nabla\xi\right\rangle f\,dx =2\int_{\Omega}\left\langle \left(D^{2}\phi\right)^{-1}\nabla\xi,\nabla\xi\right\rangle f\,dx, \end{align}\] which finishes the proof of Theorem 4. ◻

Proof. The formula for \(D_\phi\Gamma\) is well known and can be easily obtained by setting \(\phi_s = \phi+s\xi\) and differentiating with respect to \(s\). Thus it is sufficient to show \(\Gamma\) is \(m\) times continuously differentiable.

We first show the domain and co-domain of \(\Gamma\) are as claimed. Given a fixed \(\phi\in \mathcal{C}^{2,\beta}_{u}(\overline{\Omega})\), there is a constant \(0<\Lambda\) such that \(\Lambda^{-1} < D^2\phi(x)< \Lambda\) for all \(x\in \overline{\Omega}\). Also note that the function \(\log \det\) belongs to \(\mathcal{C}^\infty(\mathcal{M}_{\Lambda})\) with all derivatives bounded when restricted to the convex set \[\mathcal{M}_{\Lambda}:= \left\{ A\in\mathbb{R}^{d\times d}\mid \Lambda^{-1}\,\mathop{\mathrm{Id}}<A<\Lambda\,\mathop{\mathrm{Id}}\right\},\] therefore \(\log \det D^2\phi\in \mathcal{C}^{0,\beta}(\overline{\Omega})\) as a consequence of Lemma 11. Similarly, using that \(f_t\), \(g_t\) are bounded from above and below, Lemma 11 shows that \(\log g_t(\nabla \phi)\in \mathcal{C}^{1,\beta}(\overline{\Omega})\), \(\log f_t\in \mathcal{C}^{m,\beta}(\overline{\Omega})\) and \(\omega_*(\nabla\phi)\in\mathcal{C}^{1,\beta}(\partial\Omega)\).

Let us now show regularity of \(\Gamma\). In order to verify the hypothesis of Lemma 12, we consider the linear and bounded extension operator \[\mathcal{E}_0:\mathcal{C}^{m,\alpha}(\overline{\Omega}_*)\longrightarrow C^{m,\alpha}(\mathbb{R}^d).\] as constructed in [45]. We can then view \(\Gamma=\mathcal{G}\circ \mathcal{F}\) as a composition of two mappings defined by \[\begin{align} \mathcal{F}: I\times\mathcal{C}_{u}^{2,\beta}(\overline{\Omega}) & \to\mathcal{C}^{1,\beta}(\overline{\Omega};\mathbb{R}^{d}) \times\mathcal{C}^{0,\beta}(\overline{\Omega};\mathbb{R}^{d\times d}_{\textrm{sym}})\times \mathcal{C}^{1,\alpha}(\overline{\Omega})\times \mathcal{C}^{m,\alpha}(\overline{\Omega_*^\delta})\notag\\ \mathcal{F}(t,\phi) :&=(\nabla\phi,D^{2}\phi,f_{t},\mathcal{E}_0 g_{t}|_{\Omega_*^\delta}), \end{align}\] where \(\mathbb{R}^{d\times d}_{\textrm{sym}}\) denotes the set of symmetric matrices, and \(\Omega_*^\delta\) is a convex neighborhood of \(\overline{\Omega}_*\) small enough to guarantee that \(\mathcal{E}_0g_t > \frac{a}{2}\) on \(\Omega_*^\delta\), and \[\begin{align} &\mathcal{G}:\mathcal{U} \times\mathcal{C}_{>0}^{0,\beta}(\overline{\Omega};\mathbb{R}^{d\times d}_{\textrm{sym}}) \times\mathcal{C}_{>0}^{1,\alpha}(\overline{\Omega}) \times\mathcal{C}_{>0}^{m,\alpha}(\overline{\Omega_*^\delta}) \to\mathcal{C}^{0,\beta}(\overline{\Omega})\times\mathcal{C}^{1,\beta}(\partial\Omega)\\ &\mathcal{G}(p,H,f,g) :=\left(\log\det H+\log(g\circ p)-\log f,\omega_{*}(p)\vert_{\partial\Omega}\right), \end{align}\] where \(\mathcal{C}_{>0}^{m,\alpha}(\overline{\Omega_*^\delta})\) is the subset of \(\mathcal{C}^{m,\alpha}(\overline{\Omega_*^\delta})\) of functions bounded away from \(0\), \(\mathcal{C}_{>0}^{0,\beta}(\overline{\Omega};\mathbb{R}^{d\times d}_{\textrm{sym}})\) is the subset of \(\mathcal{C}^{0,\beta} (\overline{\Omega};\mathbb{R}^{d\times d}_{\textrm{sym}})\) which are positive definite uniformly over \(\overline{\Omega}\), and \(\mathcal{U}\) is the open subset of \(p\in \mathcal{C}^{1,\beta}(\overline{\Omega};\mathbb{R}^{d})\) with \(p(\overline{\Omega})\subset \Omega_*^{\delta/2}\).

It is easy to see that \(\mathcal{F}\) is \(m\) times differentiable, since \(t\mapsto f_t,g_t\) are differentiable by hypothesis and the maps \(\phi\mapsto (\nabla \phi, D^2\phi)\), \(g\mapsto \mathcal{E}_0g|_{\Omega_*^\delta}\) are linear (hence infinitely Frechet differentiable). Note that \(\mathcal{F}\) takes values in the set \[\mathcal{U} \times\mathcal{C}_{>0}^{0,\beta}(\overline{\Omega};\mathbb{R}^{d\times d}_{\textrm{sym}})\times \mathcal{C}_{>0}^{1,\alpha}(\overline{\Omega})\times \mathcal{C}_{>0}^{m,\alpha}(\overline{\Omega_*^\delta}),\] therefore, by Lemma 11, it is sufficient to show that \(\mathcal{G}\) is \(m\) times differentiable, which we can do term by term.

To this end, we now show that the mapping \[\label{eq:log95map} (p,g)\in\left[\mathcal{U}\subset\mathcal{C}^{1,\beta}(\overline{\Omega};\mathbb{R}^{d})\right]\times\mathcal{C}_{>0}^{m,\alpha}(\overline{\Omega_{*}^{\delta}}) \mapsto\log(g\circ p)\in\mathcal{C}^{0,\beta}(\overline{\Omega})\tag{14}\] is \(m\) times Frechet differentiable. To do so, we view it as the composition \(\mathcal{L}\circ \mathcal{H}\), where \[\begin{align} \mathcal{H}&:\left[\mathcal{U}\subset\mathcal{C}^{1,\beta}(\overline{\Omega};\mathbb{R}^{d})\right]\times\mathcal{C}_{>0}^{m,\alpha}(\overline{\Omega_{*}^{\delta}}) \longrightarrow\mathcal{C}_{>0}^{0,\beta}(\overline{\Omega}), \qquad \mathcal{H}(p,g)=g\circ p,\\ \mathcal{L}&: \mathcal{C}_{>0}^{0,\beta}(\overline{\Omega})\longrightarrow\mathcal{C}^{0,\beta}(\overline{\Omega}), \qquad \mathcal{L}(h)=\log h. \end{align}\] Note that \(\mathcal{H}\) is Frechet differentiable due to Lemma 12 statement [item:comp95reg] applied with \(E=F=\mathbb{R}^d\), \(G=\mathbb{R}\), \(U=\Omega\), \(V=\Omega_*^\delta\), \(r=1+\beta\), \(t=\beta\), \(s=m+\alpha\). As for \(\mathcal{L}\), we note that we can decompose \[\mathcal{C}_{>0}^{0,\beta}(\overline{\Omega}) = \cup_{\varepsilon>0}\mathcal{V}_\varepsilon,\qquad \mathcal{V}_\varepsilon:=\{f\in \mathcal{C}^{0,\beta}(\overline{\Omega})\mid \varepsilon<f<\varepsilon^{-1}\},\] and it is enough to show that \(\mathcal{L}\) is \(m\) times Frechet differentiable on each open set \(\mathcal{V}_\varepsilon\). For \(\varepsilon>0\) fixed, we may restrict the logarithm to a sub-interval of \((0,\infty)\) where all its derivatives are bounded. Then Lemma 12, statement [item:g4295comp] with \(E=\mathbb{R}^d\), \(F=G=\mathbb{R}\), \(U=\Omega\), \(V=(\frac{\varepsilon}{2},\frac{2}{\varepsilon})\), \(\mathcal{U}\) replaced by \(\mathcal{V}_\varepsilon\), \(r=t=\beta\), and \(s=\infty\) shows that \(\mathcal{L}\) is infinitely Frechet differentiable. All in all, this shows that 14 is \(m\) times Frechet differentiable

The terms \(p\mapsto \omega_*(p)\), \(H\mapsto \log \det H\) from \(\mathcal{G}\) are treated analogously. ◻

Proof. According to Lemma 8, \(D_\phi \Gamma(t, \phi_t)\xi=(p,q)\) corresponds to the problem \[\begin{cases} \displaystyle{\textrm{tr}\left[(D^{2}\phi_t)^{-1}D^{2}\xi\right]+\frac{\left\langle \nabla g_t(\nabla\phi_t), \nabla\xi\right\rangle}{g_t(\nabla\phi_t)} = p}, & \textrm{in }\Omega, \\ \left\langle \vec{n}_*(\nabla\phi_t), \nabla\xi\right\rangle = q, & \textrm{on }\partial\Omega. \end{cases}\] Let us first re-write the boundary condition. Consider the function \(H(x):=\omega_{*}\left(\nabla\phi_{t}(x)\right)\), which verifies \(H(x)<0\) for \(x\in\Omega\) and \(H\equiv0\) on \(\partial\Omega\). Then, if we fix \(x_{0}\in\partial\Omega\), for any direction \(\tau\) tangential to \(\partial\Omega\) at \(x_{0}\) we have \[0=\partial_{\tau}H(x_{0})=\left\langle D^{2}\phi_{t}(x_{0})\vec{n}_{*}\left(\nabla\phi_{t}(x_{0})\right), \tau\right\rangle.\] As a consequence, we see that \[D^{2}\phi_{t}(x_{0})\vec{n}_{*}\left(\nabla\phi_{t}(x_{0})\right)=c\vec{n}(x_{0})\] for some \(c\in\mathbb{R}\), \(c\neq 0\) depending on \(x_{0}\). This implies \[\vec{n}_{*}\left(\nabla\phi_{t}(x)\right)=\frac{\left(D^{2}\phi_{t}(x)\right)^{-1}\vec{n}(x)}{\left|\left(D^{2}\phi_{t}(x)\right)^{-1}\vec{n}(x)\right|},\] which shows the equivalence between the boundary conditions. Next, observe that \[\label{eq:magic95eq} f_t\,\left(\textrm{tr}\left[(D^{2}\phi_t)^{-1}D^{2}\xi\right]+\frac{\left\langle \nabla g_t(\nabla\phi_t), \nabla\xi\right\rangle}{g_t(\nabla\phi_t)}\right)=\relax\left[f_t(D^{2}\phi_t)^{-1}\nabla\xi\right],\tag{15}\] which is shown as a simple calculus exercise after differentiating the Monge–Ampère equation in the \(x\) variable. Since \(\phi_t\in \mathcal{C}^3(\overline{\Omega})\), there are no regularity issues. ◻

Proof. The previous Lemma 9 tells us that \(D_\phi \Gamma(t,\phi_t)\xi=(p,q)\) corresponds to the problem \[\begin{cases} \displaystyle{\relax\left[f_t(D^{2}\phi_t)^{-1}\nabla\xi\right]=pf_t}, & \textrm{in }\Omega \\ \displaystyle{\left\langle (D^{2}\phi_t)^{-1}\nabla\xi,\vec{n}\right\rangle =\left|(D^{2}\phi_t)^{-1}\vec{n}\right|q}, & \textrm{on }\partial\Omega. \end{cases}\] It is easy to see that (?? ) is a necessary condition, as \[\begin{align} \int_{\Omega}pf_{t}\,dx & =\int_{\Omega}\relax\left[f_{t}\left(D^{2}\phi_{t}\right)^{-1}\nabla\xi\right]\,dx \\ & =\int_{\partial\Omega}\left\langle \left(D^{2}\phi_{t}\right)^{-1}\nabla\xi, \vec{n}\right\rangle f_{t}\,dx =\int_{\partial\Omega}\left|\left(D^{2}\phi_{t}\right)^{-1}\vec{n}\right|qf_{t}\,dx. \end{align}\] Existence of solutions is then a simple consequence of the Lax-Milgram Theorem and elliptic regularity theory. ◻

Proof of Theorem 5. The previous analysis shows that \(D_\phi\Gamma (t, \phi_t)\) defines an isomorphism \(\mathcal{X}\longrightarrow\mathcal{Y}_t\) where \[\begin{align} \mathcal{X} & :=\left\{ \xi\in\mathcal{C}^{2,\beta}(\overline{\Omega})\mid\int_{\Omega}\xi =0\right\} \\ \mathcal{Y}_{t} & :=\left\{ (p,q)\in\mathcal{C}^{0,\beta}(\overline{\Omega})\times\mathcal{C}^{1,\beta}(\partial\Omega)\mid\int_{\Omega}pf_{t}=\int_{\partial\Omega}\left|(D^{2}\phi_{t})^{-1}\vec{n}\right|qf_{t}\right\} . \end{align}\] This is not enough to apply the implicit function theorem directly, but we are able to apply it to a modified operator. Without loss of generality, let us show that \(t\mapsto \phi_t\) is a smooth map in a neighborhood of \(t=0\). To do this we define \[\tilde{\Gamma}:[0,1]\times\left(\mathcal{C}_{>0}^{2,\beta}(\overline{\Omega})\cap\mathcal{X}\right)\longrightarrow\mathcal{Y}_{0}, \quad \tilde{\Gamma}^{(1)}(t, \phi)=\Gamma_{t}^{(1)}(\phi)+C_{\phi}, \quad \tilde{\Gamma}^{(2)}(t, \phi)=\Gamma_{t}^{(2)}(\phi),\] where \(\tilde{\Gamma}=(\tilde{\Gamma}^{(1)}, \tilde{\Gamma}^{(2)})\), and the constant \(C_\phi\in \mathbb{R}\) is chosen so \(\tilde{\Gamma}(t, \phi)\in\mathcal{Y}_{0}\). In other words, \[\tilde{\Gamma}^{(1)}(t,\phi)= \Gamma_{t}^{(1)}(\phi)+\int_{\partial\Omega}\left|(D^{2}\phi_{0})^{-1}\vec{n}\right| \Gamma_{t}^{(2)}(\phi)f_0-\int_{\Omega}\Gamma_{t}^{(1)}(\phi)f_0.\] This functional is now different from the original \(\Gamma\), so we must check if it verifies all the hypothesis of the implicit function theorem. Its \(\mathcal{C}^m\)-smoothness is a simple consequence of the \(\mathcal{C}^m\)-smoothness of \(\Gamma\). Also, given \(\xi\in\mathcal{X}\), we have \[D_{\phi}\tilde{\Gamma}^{(1)}(t, \phi)\xi=D_{\phi}\Gamma_{t}^{(1)}(\phi)\xi+ \int_{\partial\Omega}\left|(D^{2}\phi_{0})^{-1}\vec{n}\right|\big(D_{\phi}\Gamma_{t}^{(2)}(\phi)\xi\big) f_0 -\int_{\Omega}\big(D_{\phi}\Gamma_{t}^{(1)}(\phi)\xi\big) f_0.\] Note that \(D_\phi \Gamma (t, \phi_t)\xi\in \mathcal{Y}_t\), so the last two expressions cancel out if we set \(\phi=\phi_0\), \(t=0\). Therefore, we have \[D_\phi\tilde{\Gamma}(0, \phi_0) = D_\phi\Gamma(0, \phi_0)|_{\mathcal{X}},\] which shows that \(D_\phi\tilde{\Gamma}(0, \phi_0)\) is an isomorphism and that the implicit function theorem can be applied. This implies the existence of some \(\delta>0\) and a smooth mapping \[t\in (-\delta,\delta) \mapsto \tilde{\phi}_t\in \mathcal{C}^{2,\beta}_{>0}(\overline{\Omega}),\] such that \(\tilde{\Gamma}(t,\tilde{\phi}_t)\equiv (0, 0)\). As a last step, we show that we actually have \(\tilde{\phi}_t = \phi_t\). We know that for each \(t\), \[\det D^{2}\tilde{\phi}_t(x)=C_t\frac{f_{t}(x)}{g_{t}\big(\nabla\tilde{\phi}_{t}(x)\big)},\qquad\nabla\tilde{\phi}_t(\Omega)=\Omega_*.\] for some constant \(C_t>0\). However, using the change of variables \(y=\nabla\tilde{\phi}_t(x)\) \[C_t=\int_{\Omega}C_tf_t(x)\,dx=\int_{\Omega}g_t\left(\nabla\tilde{\phi}_t(x)\right)\det D^{2}\tilde{\phi}_t(x)\,dx=\int_{\Omega_{*}}g_t(y)\,dy=1,\] which implies \(\tilde{\phi}_t = \phi_t\) by uniqueness of solutions, finishing the proof of Theorem 5. ◻

Proof of Corollary 1. First let the linear and bounded extension operators \[\mathcal{E}:\mathcal{C}^{0,\gamma}(\overline{\Omega})\longrightarrow\mathcal{C}^{0,\gamma}(\mathbb{R}^{d}),\qquad\mathcal{E}_{*}:\mathcal{C}^{0,\gamma}(\overline{\Omega}_{*})\longrightarrow\mathcal{C}^{0,\gamma}(\mathbb{R}^{d})\] correspond to \(\mathcal{E}_0\) as constructed in [45]; these operators are defined independent of \(\gamma\in (0, 1]\). Additionally, the extended function can be expressed as a finite sum outside of the original domain, hence differentiation in \(t\) commutes with either extension operator. Thus we may assume \(f_t\), \(g_t\) are defined on \(\mathbb{R}^d\) with \(f_t\), \(g_t\geq \frac{a}{2}\) on a neighborhood of \(\overline{\Omega}\), \(\overline{\Omega}_*\), such that \[\begin{align} \label{eqn:32extended32bounded} \sup_{t\in (0, 1)}\max\left\{ \Vert f_t\Vert_{\mathcal{C}^{0,\alpha}(\mathbb{R}^d)},\Vert \partial_tf_t\Vert_{L^\infty(\mathbb{R}^d)}, \Vert g_t\Vert_{\mathcal{C}^{0,\alpha}(\mathbb{R}^d)}, \Vert \partial_tg_t\Vert_{L^\infty(\mathbb{R}^d)}\right\}<\infty, \end{align}\tag{16}\] and with \(\partial_tf_t\), \(\partial_tg_t\) are continuous on \((0, 1)\times \mathbb{R}^d\). We further extend \(f_t\), \(g_t\) to \(t\in \mathbb{R}\) by setting \(f_t=f_0\) for \(t<0\) and \(f_t=f_1\) for \(t>1\), and similarly for \(g_t\). Next, let \({\eta\in\mathcal{C}_{c}^{\infty}\left(\mathbb{R}\times\mathbb{R}^{d}\right)}\), be nonnegative, define its rescaling \(\eta^{\varepsilon}(t, x)=\varepsilon^{-(d+1)}\eta(t/\varepsilon, x/\varepsilon)\), and let \[\begin{align} f_{t}^{\varepsilon}(x):&=\left(\int_\Omega (f*\eta^{\varepsilon})(t, y)dy\right)^{-1}(f*\eta^{\varepsilon})(t, x),\\ g_{t}^{\varepsilon}(x):&=\left(\int_{\Omega_*} (g*\eta^{\varepsilon})(t, y)dy\right)^{-1}(g*\eta^{\varepsilon})(t, x), \end{align}\] so that \(f_{t}^{\varepsilon}\), \(g_{t}^{\varepsilon}\) are probability densities on \(\Omega\), \(\Omega_{*}\) for each \(t\in [0, 1]\); above the convolutions are taken in both time and space. Then for each \(0<t<1\) we have that \(f_{t}^{\varepsilon}\), \(g_{t}^{\varepsilon}\) are bounded from below by \(\frac{a}{2}\) on \(\Omega\) and \(\Omega_*\) respectively, and belong to \(\mathcal{C}^{\infty}\) in \((t, x)\in \mathbb{R}\times \mathbb{R}^d\).

Then, one can see that \(t\mapsto f_t^\varepsilon\), \(t\mapsto g_t^\varepsilon\) are infinitely Frechet differentiable as curves in any space \(\mathcal{C}^{k,\gamma}\) (see Proposition 6 for the case \(k=0\), and note the argument can be extended to any \(k\in \mathbb{N}\).) In particular, we satisfy the hypotheses of Theorem 5 with \(m=1\). Thus if \(\phi_{t}^{\varepsilon}\) is the Brenier potential from \(f^\varepsilon_t\) to \(g^\varepsilon_t\) with integral zero, it is differentiable with respect to \(t\). Fixing \(t\in (0, 1)\), by integration by parts with ?? , \(\xi_{t}^{\varepsilon}:=\partial_{t}\phi_{t}^{\varepsilon}\) satisfies for any test function \(\theta\in \mathcal{C}^\infty(\overline{\Omega})\), \[\begin{align} \begin{aligned} \int_\Omega \left\langle (D^{2}\phi^\varepsilon_t)^{-1}\nabla\xi^\varepsilon_t, \nabla \theta\right\rangle f^\varepsilon_tdx&= \int_\Omega\theta\left(\frac{\partial_{t}g^\varepsilon_t\left(\nabla\phi^\varepsilon_t\right)}{g^\varepsilon_t\left(\nabla\phi^\varepsilon_t\right)}-\frac{\partial_{t}f^\varepsilon_t}{f^\varepsilon_t}\right)f^\varepsilon_tdx\\ &=\int_\Omega\left(\partial_{t}g^\varepsilon_t(\theta\circ\left(\nabla\phi^\varepsilon_t\right)^{-1})-\theta\partial_tf^\varepsilon_t\right)dx \end{aligned} \label{eqn:32mollified32pde} \end{align}\tag{17}\] Let \(\psi^\varepsilon_t\) be the Brenier potential from \(g^\varepsilon_t\) to \(f^\varepsilon_t\) which is zero at some designated point in \(\Omega_*\) for all \(\varepsilon\), then by Lemma 1 we obtain \[\begin{align} \left\lVert \phi_{t}^{\varepsilon}\right\rVert _{\mathcal{C}^{2,\alpha}(\overline{\Omega})},\;\left\lVert \psi_{t}^{\varepsilon}\right\rVert _{\mathcal{C}^{2,\alpha}(\overline{\Omega}_*)}\leq C, \end{align}\] where \(C\) depends on \(\left\Vert f_{t}^{\varepsilon}\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\), \(\left\Vert g_{t}^{\varepsilon}\right\Vert _{\mathcal{C}^{0,\alpha}(\overline{\Omega}_{*})}\), and the lower bound of \(f_{\varepsilon},g_{\varepsilon}\), which by construction are independent of \(\varepsilon>0\). Thus we can apply Arzelà–Ascoli to find subsequences which, for any \(0<\beta<\alpha\), converge as \(\varepsilon\to 0\) in \(\lVert\cdot\rVert_{\mathcal{C}^{2,\beta}(\overline{\Omega})}\) to \(\phi_t\) and \(\lVert\cdot\rVert_{\mathcal{C}^{2,\beta}(\overline{\Omega}_{*})}\) to \(\psi_t\) respectively, which are clearly Brenier potentials from \(f_t\) to \(g_t\) and \(g_t\) to \(f_t\) respectively (for example, use that \(\phi_t\) and \(\psi_t\) are convex as a limit of convex functions, then [15] combined with the uniqueness in [16]). Since \(\nabla \psi^\varepsilon_t=(\nabla \phi^\varepsilon_t)^{-1}\) and \(D^2 \psi^\varepsilon_t=(D^2 \phi^\varepsilon_t)^{-1}\) by, this also yields \[\begin{align} \Lambda^{-1}\mathop{\mathrm{Id}}\leq (D^{2}\phi_{t}^{\varepsilon})^{-1}\leq\Lambda\mathop{\mathrm{Id}} \end{align}\] for \(\Lambda>0\) independent of \(\varepsilon\). Then by Lemmas 5 and 6 with 17 , we can find a subsequence such that \(\xi^\varepsilon_t\) converges in \(\lVert\cdot\rVert_{\mathcal{C}^{1,\beta}(\overline{\Omega})}\) for any \(\beta <\alpha\), and the limit can be seen as \(\xi_t:=\partial_t\phi_t\) by taking \(\varepsilon\to 0\) in the expression \[\int_0^1\int_{\Omega}\xi^{\varepsilon}(t, x)\theta(t,x)\,dx\,dt =-\int_0^1\int_{\Omega}\phi^{\varepsilon}(t, x)\partial_{t}\theta(t, x)\,dx\,dt,\] where \(\theta\in\mathcal{C}_{c}^{\infty}\left((0, 1)\times\Omega\right)\) is an arbitrary test function. Thus the left hand side of 17 converges to the left hand side of ?? as \(\varepsilon\to 0\). By the regularity of \(\partial_tf_t\), the sequence \(\partial_tf^\varepsilon_t\) converges pointwise to \(\partial_tf_t\) on \(\overline{\Omega}\) (viewing \(\{t\}\times \overline{\Omega}\) as a compact subset of \((0, 1)\times \mathbb{R}^d\)), and similarly for \(\partial_tg^\varepsilon_t\). Thus by the bound 16 , we may apply dominated convergence and let \(\varepsilon\to 0\) in 17 , which yields the weak formulation of ?? as \(\varepsilon\to 0\), finishing the proof. ◻

Proof. The clear candidate for differential \(d\mathcal{F}(t):\mathbb{R}\longrightarrow\mathcal{C}^{0,\beta}(\overline{\Omega})\) is \(\left\langle d\mathcal{F}(t),h\right\rangle =\partial_{t}f(t,\cdot)h\), which we can identify with \(\mathcal{F}'(t):=\partial_{t}f(t,\cdot)\in\mathcal{C}^{0,\beta}(\overline{\Omega})\). We must show that for \(t\in (0,1)\), \[\label{eq:defn95Frechet95diff} \lim_{h\to0}\left\Vert \frac{\mathcal{F}(t+h)-\mathcal{F}(t)}{h}-\mathcal{F}'(t)\right\Vert _{\mathcal{C}^{0, \beta}(\overline{\Omega})}=0.\tag{18}\] Using the inequality \(\Vert g\Vert_{\mathcal{C}^{0, \beta}(\overline{\Omega})}\leq C\Vert g\Vert_{\mathcal{C}^{0}(\overline{\Omega})}^{1-\frac{\beta}{\alpha}}\Vert g\Vert_{\mathcal{C}^{0, \alpha}(\overline{\Omega})}^{\frac{\beta}{\alpha}}\), we obtain when \(0<t+h<1\), \[\begin{align} \left\Vert \frac{\mathcal{F}(t+h)-\mathcal{F}(t)}{h}-\mathcal{F}'(t)\right\Vert _{\mathcal{C}^{0, \beta}(\overline{\Omega})} & =\left\Vert \int_{0}^{1}\Big(\partial_{t}f(t+sh,\cdot)-\partial_{t}f(t,\cdot)\Big)\,ds\right\Vert _{\mathcal{C}^{0, \beta}(\overline{\Omega})}\\ & \leq CM^{\frac{\beta}{\alpha}}\int_{0}^{1}\left\Vert \partial_{t}f(t+sh,\cdot)-\partial_{t}f(t,\cdot)\right\Vert _{\mathcal{C}^{0}(\overline{\Omega})}^{1-\frac{\beta}{\alpha}}\,ds. \end{align}\] Note that \({\displaystyle \lim_{h\to0}}\left\Vert \partial_{t}f(t+sh,\cdot)-\partial_{t}f(t,\cdot)\right\Vert _{\mathcal{C}^{0}(\overline{\Omega})}=0\) for each \(s\) due to the uniform continuity of \(\partial_{t}f\), so the dominated convergence theorem shows 18 as desired. The continuity of the differential follows from a similar calculation, as \[\left\Vert \mathcal{F}'(t+h)-\mathcal{F}'(t)\right\Vert _{\mathcal{C}^{0, \beta}(\overline{\Omega})}\leq CM^{\frac{\beta}{\alpha}}\left\Vert \partial_{t}f(t+h,\cdot)-\partial_{t}f(t,\cdot)\right\Vert _{\mathcal{C}^{0}(\overline{\Omega})}^{1-\frac{\beta}{\alpha}} \overset{h\to0}{\longrightarrow}0.\qedhere\] ◻

Proof of Lemma 1. Let \(\phi\) be the solution of ?? and fix \(\alpha<\beta\). By Lemma 15, there exists \(\delta>0\) such that \(\phi \in \mathcal{C}^{1,\delta}\) with an estimate on \(\Vert \phi\Vert_{\mathcal{C}^{1,\delta}(\overline{\Omega})}\). Therefore, letting \(\tilde{\alpha}:=\min(\alpha, \delta\beta)<\beta\), Lemmas 11 and 13 show that \[\left\Vert f/g(\nabla\phi)\right\Vert _{\mathcal{C}^{0,\tilde{\alpha}}(\overline{\Omega})} \leq C\left(d,a,\Omega,\Omega_{*},\left\Vert f\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega})},\left\Vert g\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega}_*)}\right)=:M.\] Now note that \(\phi=\phi_{f/g(\nabla\phi)}\in C^{2, \alpha}(\overline{\Omega})\) in the notation of Lemma 14. Since the set \[\left\{ h\in\mathcal{C}^{0,\tilde{\alpha}}(\overline{\Omega})\mid h\geq a/\lVert g\rVert_{\mathcal{C}^{0}(\overline{\Omega}_*)},\;\left\Vert h\right\Vert _{\mathcal{C}^{0,\tilde{\alpha}}(\overline{\Omega})}\leq M\right\}\] is compact as a subset of \(\mathcal{C}^{0, \tilde{\alpha}-\varepsilon}(\overline{\Omega})\), for \(0<\varepsilon<\tilde{\alpha}\) arbitrary, applying Lemma 14 shows \[\label{eq:workaround95estimate95only95f} \left\Vert \phi\right\Vert _{\mathcal{C}^{2,\tilde{\alpha}-\varepsilon}(\Omega)} =\left\Vert \phi_{f/g(\nabla\phi)}\right\Vert _{\mathcal{C}^{2,\tilde{\alpha}-\varepsilon}(\Omega)} \leq C\left(d,a,\Omega,\Omega_{*},\left\Vert f\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega})}, \left\Vert g\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega}_*)},\alpha,\beta\right).\tag{20}\] In particular this implies a \(\mathcal{C}^2\) bound on \(\phi\), thus we can bootstrap and use Lemmas 11, 13 again to conclude \[\left\Vert f/g(\nabla\phi)\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega})}\leq C\left(d,a,\Omega,\Omega_{*},\left\Vert f\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega})},\left\Vert g\right\Vert _{\mathcal{C}^{0,\beta}(\overline{\Omega}_*)},\alpha,\beta\right),\] which can be combined with Lemma 14 again to produce the desired estimate. ◻

Proof of Lemma 14. The idea is to apply the implicit function theorem to the functional \[\Gamma:\mathcal{C}_{>0}^{0,\alpha}(\overline{\Omega})\times\mathcal{C}_{>0}^{2,\alpha}(\overline{\Omega})\longrightarrow\mathcal{C}^{0,\alpha}(\Omega),\qquad\Gamma(f,\phi):=\left(\det D^{2}\phi-|\Omega_{*}|\left\langle f\right\rangle _{\Omega}^{-1}f,\omega_{*}(\nabla\phi)\right),\] where recall that \(\omega_{*}\) is a convex defining function of \(\Omega_{*}\). The domain of \(\Gamma\) is defined by the open sets \[\begin{align} \mathcal{C}_{>0}^{0,\alpha}(\overline{\Omega}) & :=\left\{ f\in\mathcal{C}^{0,\alpha}(\overline{\Omega})\mid\inf_{\Omega}f>0\right\}, \\ \mathcal{C}_{u}^{2,\alpha}(\overline{\Omega}) & :=\left\{ \phi\in\mathcal{C}^{2,\alpha}(\overline{\Omega})\mid\inf_{x\in\Omega,|e|=1}\left\langle D^{2}\phi(x)e, e\right\rangle >0\right\} . \end{align}\] From Lemma 12 we immediately see that \(\Gamma\) is of class \(\mathcal{C}^1\) with differential \[D_{\phi}\Gamma(f,\phi)\xi=\left(\textrm{tr}\left[A[\phi]D^{2}\xi\right],\left\langle \vec{n}_{*}\left(\nabla\phi\right), \nabla\xi\right\rangle\right), \qquad A[\phi] = (\det D^2\phi) \left(D^2\phi\right)^{-1},\] where \(A[\phi]\) is the cofactor matrix of \(D^2\phi\). Therefore, showing that \(D_\phi \Gamma (f,\phi)\) is a diffeomorphism is equivalent to showing that for each \(p\in\mathcal{C}^{0,\alpha}(\overline{\Omega})\), \(q\in\mathcal{C}^{1,\alpha}(\partial\Omega)\), there is a unique solution \(\xi \in \mathcal{C}^{2,\alpha}(\overline{\Omega})\) of \[\label{eq:PDE195workaround} \begin{cases} \textrm{tr}\left[A[\phi]D^{2}\xi\right]=p, & \textrm{in }\Omega\\ \left\langle \vec{n}_{*}(\nabla\phi), \nabla\xi\right\rangle=q, & \textrm{on }\partial\Omega. \end{cases}\tag{21}\] To study this problem, we will re-write its boundary condition and note that it is equivalent to the same problem written in divergence form. This would be immediately obvious if \(\phi\) were \(\mathcal{C}^{3}\), since it is known that cofactor matrices are divergence free.

Next, we consider a fixed \(f_0\in \mathcal{C}_{>0}^{0,\alpha}(\overline{\Omega})\), its corresponding potential \(\phi_{f_0}\in \mathcal{C}^{2,\alpha}(\overline{\Omega})\) and will show that the mapping \(f\mapsto\phi_{f}\) is \(\mathcal{C}^{1}\) in a neighborhood of \(f_0\). To do this, first note the previous claim shows that \(D_{\phi}\Gamma(f_0,\phi_{f_0})\) is an isomorphism from \(\mathcal{X}\) to \(\mathcal{Y}_{\phi_0}\) defined by \[\begin{align} \mathcal{X} & :=\left\{ \xi\in\mathcal{C}^{2,\alpha}(\overline{\Omega})\mid\int_{\Omega}\xi=0\right\},\\ \quad\mathcal{Y}_{\phi_{0}}& := \left\{ (p,q)\in\mathcal{C}^{0,\alpha}(\overline{\Omega})\times\mathcal{C}^{1,\alpha}(\partial\Omega)\mid\int_{\Omega}p+\int_{\partial\Omega}\left|A[\phi_{f_0}]\vec{n}\right|q=0\right\} . \end{align}\] Then define, writing \(\Gamma(f, \phi)=(\Gamma^{(1)}(f,\phi), \Gamma^{(2)}(f,\phi))\), \[\tilde{\Gamma}:\mathcal{C}_{>0}^{0,\alpha}(\overline{\Omega})\times\left(\mathcal{X}\cap\mathcal{C}_{u}^{2,\alpha}(\overline{\Omega})\right)\longrightarrow\mathcal{Y}_{\phi_{0}},\qquad\tilde{\Gamma}(f,\phi)=\left(\Gamma^{(1)}(f,\phi)-C_{f,\phi}\left\langle f\right\rangle _{\Omega}^{-1}f,\Gamma^{(2)}(f,\phi)\right),\] where \(C_{f,\phi}\) is a constant chosen to ensure \(\tilde{\Gamma}(f,\phi)\in\mathcal{Y}_{\phi_{0}}\), i.e. \[C_{f,\phi}=\int_{\Omega}\Gamma^{(1)}(f,\phi)+\int_{\partial\Omega}\left|A[\phi_{0}]\vec{n}\right|\Gamma^{(2)}(f,\phi).\] We now check that \(\tilde{\Gamma}\) satisfies the hypothesis of the implicit function theorem. The smoothness of \(\tilde{\Gamma}\) is a consequence of the smoothness of \(\Gamma\), and we calculate \[D_{\phi}\tilde{\Gamma}(f,\phi)\xi=D_{\phi}\Gamma(f,\phi)\xi+\left(\int_{\Omega}D_{\phi}\Gamma^{(1)}(f,\phi)\xi+\int_{\partial\Omega}\left|A[\phi_{0}]\vec{n}\right|D_{\phi}\Gamma^{(2)}(f,\phi)\xi,0\right)\left\langle f\right\rangle _{\Omega}^{-1}f.\] Note that the second term is \(0\) at \((f_0,\phi_{0})\) because \(D_{\phi}\Gamma(f_0,\phi_{0})\xi\in\mathcal{Y}_{\phi_{0}}\), hence \[D_{\phi}\tilde{\Gamma}(f_0,\phi_{0})=D_{\phi}\Gamma(f_0,\phi_0)|_{\mathcal{X}},\] showing that \(D_{\phi}\tilde{\Gamma}(f_0,\phi_{0}):\mathcal{X}\longrightarrow\mathcal{Y}_{\phi_{0}}\) is a diffeomorphism and the implicit function theorem applies. Thus there is a neighborhood \(\mathcal{U}\subset\mathcal{C}_{>0}^{0,\alpha}(\Omega)\) of \(f_0\) and a \(\mathcal{C}^1\) mapping \[f\in\mathcal{U}\mapsto\tilde{\phi}_{f}\in\mathcal{C}_{u}^{2,\alpha}(\Omega)\cap\mathcal{X},\quad\textrm{such that }\tilde{\Gamma}(f,\tilde{\phi}_{f})=0.\] Finally, we verify that \(\tilde{\phi}_f = \phi_f\). From \(\tilde{\Gamma}(f,\tilde{\phi}_{f})=0\) we conclude that for some \(C\in \mathbb{R}\), \[\begin{cases} \det D^{2}\tilde{\phi}_{f}=\left\langle f\right\rangle _{\Omega}^{-1}\left(|\Omega_{*}|-C\right)f, & \textrm{in }\Omega\\ \nabla\tilde{\phi}_{f}(\Omega)=\Omega_{*}. \end{cases}\] Since \(\tilde{\phi}_{f}\) is strictly convex, \(\nabla\tilde{\phi}_f\) is a one-to-one mapping and we can perform the change of variables \(y=\nabla\tilde{\phi}_{f}(x)\), which shows \[|\Omega_{*}|=\int_{\nabla\tilde{\phi}_{f}(\Omega)}1\,dy=\int_{\Omega}\det D^{2}\tilde{\phi}_{f}(x)\,dx=\left\langle f\right\rangle _{\Omega}^{-1}\left(|\Omega_{*}|-C\right)\int_{\Omega}f=|\Omega_{*}|-C.\] Therefore, we have \(C=0\) and \(\tilde{\phi}_f = \phi_f\), finishing the proof. ◻