Fix an integer $d \geq 2$. Let $K \subset \mathbb{R}^d$ be a compact convex set with nonempty interior, and let $\rho$ denote the uniform probability measure
$$
\rho = \frac{1_K(x) \, dx}{|K|}.
$$
Fix a nonempty compact set $Y \subset \mathbb{R}^d$; write $\mathrm{Prob}(Y)$ for the Borel probability measures on $Y$.

For $\mu \in \mathrm{Prob}(Y)$, the quadratic-cost optimal transport map from $\rho$ to $\mu$ is the minimizer $T_\mu$ of
$$
\int_K |x - T(x)|^2 \, d\rho(x)
$$
over all Borel $T : K \to Y$ with $T_{\#} \rho = \mu$. The minimizer exists, is unique up to $\rho$-null sets, and equals the gradient of a convex function on $K$. Both the existence and the gradient structure go back to Brenier.
pdfrefs_brenier/mccann-1997-advances.pdf
pdfrefs_brenier/carlier-galichon-santambrogio-2010-siam-j-math-anal.pdf

For compactly supported $\mu$ and $\nu$, the quadratic Wasserstein distance is
$$
W_2(\mu,\nu) = \left( \inf_{\pi \in \Pi(\mu,\nu)} \int_{\mathbb{R}^d \times \mathbb{R}^d} |y-z|^2 \, d\pi(y,z) \right)^{1/2},
$$
where $\Pi(\mu,\nu)$ is the set of couplings of $\mu$ and $\nu$.

The stability question: how strongly does $T_\mu$ depend on $\mu$? The answer below is one-third Holder in $W_2$, with a constant that depends only on $K$ and $Y$.

There is a finite constant $C(K,Y)$ such that for all $\mu,\nu \in \mathrm{Prob}(Y)$,
$$
\| T_\mu - T_\nu \|_{L^2(\rho)} \leq C(K,Y) \, W_2(\mu,\nu)^{1/3}.
$$
The targets may be atomic, singular, or absolutely continuous. The exponent $1/3$ cannot be increased: if $K = Y = [-1,1]^d$ and $\alpha > 1/3$, the ratio $\| T_\mu - T_\nu \|_{L^2(\rho)} / W_2(\mu,\nu)^\alpha$ is unbounded over pairs of measures supported on three atoms.
pdfrefs_brenier/arxiv_2103.05934v2-delalande-merigot-2023.pdf
pdfrefs_brenier/arxiv_2404.02855v1-divol-nilesweed-pooladian-2025.pdf

The bound is uniform in atom count, in positive atom masses, in site separation, and in density bounds. Neither strict convexity nor smoothness of $\partial K$ enters, and $Y$ need not be convex.

An explicit constant appears in the interpolation step later. With $K \subset B(0,R)$, $Y \subset B(0,L)$ and $P_K$ the sum of the $(d-1)$-volumes of the coordinate projections of $K$ onto the hyperplanes orthogonal to the coordinate axes, one admissible value is
$$
C_* = \left( 12 \, d \, R^2 \left( 1 + \sqrt{C_0} \right)^2 + \frac{28 L^2 P_K}{|K|} \right)^{1/2},
$$
with $C_0 = 162$.

Earlier results, in roughly historical order. Near smooth data and along target curves, Ambrosio proved a local half-Holder estimate for the map; the statement appears as Corollary 3.4 in the survey of Gigli. Berman proved Holder stability for varying sources against a fixed target whose density is bounded below, on a convex domain. Merigot, Delalande, and Chazal combined Berman's estimates with dualization to obtain a dimension-independent $W_1^{2/15}$ map bound for a convex source. Delalande and Merigot lowered the exponent further, to $W_1^{1/6}$, for source densities bounded above and below, and introduced the gradient-interpolation step that underlies the final part of the argument below. Letrouit and Merigot extended potential stability to John domains, obtaining a corresponding map estimate when the boundary is rectifiable with finite $(d-1)$-dimensional measure. Divol, Niles-Weed, and Pooladian obtained a $W_2^{1/3}$ map bound for regular source densities and nondegenerate finite targets, with a constant depending on atom count, lower mass bounds, and angular geometry. Cazelles, Pauwels, and Portales proved a $W_1^{1/5}$ bound for bounded source densities. Han and Zhu quote a $W_1^{1/4}$ bound of Merigot. Chewi, Eichinger, and Pooladian established near-Lipschitz stability for the Kim-Milman flow map, giving a comparable Lipschitz-order bound in a related framework. Caja-Lopez, Delgadino, and Kitagawa proved Lipschitz $L^2$ stability for uniformly positive Holder-continuous densities on suitably regular convex domains.
pdfrefs_brenier/merigot-delalande-chazal-2020-pmlr-v108.pdf
pdfrefs_brenier/arxiv_2103.05934v2-delalande-merigot-2023.pdf
pdfrefs_brenier/arxiv_2411.04908v3-letrouit-merigot-2024.pdf
pdfrefs_brenier/arxiv_2404.02855v1-divol-nilesweed-pooladian-2025.pdf
pdfrefs_brenier/arxiv_2604.22366v2-cazelles-pauwels-portales-2026.pdf
pdfrefs_brenier/merigot-2026-hal.pdf
pdfrefs_brenier/arxiv_2606.23037v1-han-zhu-2026.pdf
pdfrefs_brenier/arxiv_2606.23383v1-chewi-eichinger-pooladian-2026.pdf
pdfrefs_brenier/arxiv_2605.24232v1-caja-lopez-delgadino-kitagawa-2026.pdf

A three-atom construction on a fixed cube gives the sharpness assertion, disproving every exponent above $1/3$. Then the upper bound is built from four pieces: finite-optimal-map stability, cell calculus for prescribed masses, a uniform moment estimate that removes any dependence on the inverse graph-Laplacian, and a coordinate interpolation inequality that converts potential stability into map stability. Taking limits over finite nets extends everything to all Borel targets.

The exponent $1/3$ is sharp. The construction uses three affine functions on a fixed cube; a small tilt of one slope moves a thin central cell while the target displacement stays second order.

Begin with an optimality criterion for gradients of maxima of affine functions. Let $u(x) = \max_{1 \leq i \leq N} \{ x \cdot y_i + c_i \}$ on a convex $K$, with distinct slopes $y_i \in \mathbb{R}^d$ and constants $c_i \in \mathbb{R}$. Let $T(x) = y_i$ at the smallest maximizing $i$, and let $\lambda = T_\# \rho$.

On the complement of the finite union of hyperplanes where two distinct slopes tie for the maximum, $u$ equals a single affine function in a neighborhood, so $\nabla u(x) = y_i$; the tie set is null for $\rho$. Extend $c(y) = c_i$ to the slopes $y_i$, and $c(y) = -\infty$ otherwise. Then $x \cdot y + c(y) \leq u(x)$ for all $(x,y)$, with equality if and only if $y = T(x)$ for $\rho$-almost every $x$ (equal slopes are merged to the largest intercept). For any coupling $\pi$ of $\rho$ and $\lambda$,
$$
\int x \cdot y \, d\pi(x,y) \leq \int_K u \, d\rho - \int_Y c(y) \, d\lambda(y),
$$
with equality for the graph plan of $T$. The cross term is the only part of the quadratic cost not fixed by the marginals, since
$$
\int |x - y|^2 \, d\pi(x,y) = \int_K |x|^2 \, d\rho(x) + \int_Y |y|^2 \, d\lambda(y) - 2 \int x \cdot y \, d\pi(x,y).
$$
Hence the graph of $T$ is the unique optimal plan, and $T = \nabla u$ almost everywhere.
pdfrefs_brenier/mccann-1997-advances.pdf

Specialize to $K = Y = [-1,1]^d$, so that $\rho$ is uniform on $K$. Fix $0 < a < 1/2$ and put $b = a/2$. Define
$$
u_a(x) = \max\{ -x_1, x_1, a \}, \qquad \tilde u_a(x) = \max\{ -x_1, x_1, a + b x_2 \}.
$$
The slopes lie in $Y$: $-e_1, e_1, 0$ for $u_a$; $-e_1, e_1, b e_2$ for $\tilde u_a$. Set $\mu_a = (\nabla u_a)_\# \rho$ and $\nu_a = (\nabla \tilde u_a)_\# \rho$.

The central cells are $A_a = \{ x \in K : |x_1| < a \}$ and $B_a = \{ x \in K : |x_1| < a + b x_2 \}$. Since $|b x_2| \leq b = a/2 < a$ throughout $K$, the quantity $a + b x_2$ stays strictly between $0$ and $3 a / 2 < 1$. Both cells therefore live strictly inside $K$ in the first coordinate, away from the null boundaries where two affine pieces tie, and both have mass $a$:
$$
\rho(A_a) = \rho(B_a) = a.
$$
Away from $A_a$ the first map equals $\mathrm{sgn}(x_1) e_1$, and away from $B_a$ the second map equals $\mathrm{sgn}(x_1) e_1$. The two targets are three-atom measures:
$$
\mu_a = \frac{1-a}{2} (\delta_{-e_1} + \delta_{e_1}) + a \, \delta_0, \qquad
\nu_a = \frac{1-a}{2} (\delta_{-e_1} + \delta_{e_1}) + a \, \delta_{b e_2}.
$$
pdfrefs_brenier/arxiv_2404.02855v1-divol-nilesweed-pooladian-2025.pdf

The two targets share the outer atoms. Matching the outer atoms and sending $0$ to $b e_2$ costs $a b^2$. Every point of $\mathrm{supp}(\mu_a)$ has second coordinate zero, while the atom $b e_2 \in \mathrm{supp}(\nu_a)$ has second coordinate $b$ and mass $a$, so every coupling satisfies
$$
\int |y - z|^2 \, d\pi(y,z) \geq \int |y_2 - z_2|^2 \, d\pi(y,z) \geq a b^2.
$$
Consequently
$$
W_2(\mu_a, \nu_a)^2 = a b^2.
$$

For the map distance, break the pointwise identity
$$
|T_{\mu_a}(x) - T_{\nu_a}(x)|^2
$$
into four regions. On $A_a \cap B_a$ the maps are $0$ and $b e_2$, and the squared difference is $b^2$. On $A_a \setminus B_a$ the maps are $0$ and $\mathrm{sgn}(x_1) e_1$, and the squared difference is one. On $B_a \setminus A_a$ the maps are $\mathrm{sgn}(x_1) e_1$ and $b e_2$, and the squared difference is $1 + b^2$. Outside $A_a \cup B_a$ the two maps coincide.

Fix $x_2 = t$. On that slice, $A_a$ occupies an interval of length $2a$, while $B_a$ occupies an interval of length $2(a + b t)$. Elementary integration gives
$$
\rho(A_a \setminus B_a) = \rho(B_a \setminus A_a) = \frac{b}{4}, \qquad
\rho(A_a \cap B_a) = a - \frac{b}{4}.
$$
Combining,
$$
\int_K |T_{\mu_a} - T_{\nu_a}|^2 \, d\rho
= b^2 \left( a - \frac{b}{4} \right) + \frac{b}{4} + \frac{b}{4} (1 + b^2)
= \frac{b}{2} + a b^2.
$$

Substituting $b = a/2$,
$$
W_2(\mu_a, \nu_a) = \frac{a^{3/2}}{2}, \qquad
\| T_{\mu_a} - T_{\nu_a} \|_{L^2(\rho)} = \frac{\sqrt{a}}{2} \sqrt{1 + a^2}.
$$
For $\alpha > 1/3$,
$$
\frac{\| T_{\mu_a} - T_{\nu_a} \|_{L^2(\rho)}}{W_2(\mu_a, \nu_a)^\alpha}
= 2^{\alpha - 1} \sqrt{1 + a^2} \, a^{(1 - 3\alpha)/2}.
$$
The exponent is negative, so the ratio diverges as $a \downarrow 0$. The source $K$, the target container $Y$, and the number of atoms are fixed; only the geometry of three affine pieces varies. This shows the exponent $1/3$ cannot be increased.

Letrouit proved an obstruction to exponents above $1/3$ for a nonconvex source, and formulated uniform square-root stability in $W_2$ for a uniform convex source as a conjecture. The three-atom family above already rules out every exponent above $1/3$ in the uniform convex setting, so the sharp positive exponent is exactly $1/3$. The conjecture is recorded in Letrouit's study of unstable optimal transport maps, cited as reference 42 of Letrouit-Merigot.
pdfrefs_brenier/arxiv_2411.04908v3-letrouit-merigot-2024.pdf

Turn to the upper bound. Begin with finite targets. Let $\mu = \sum_{i=1}^N m_i \delta_{y_i}$ with distinct $y_i \in Y$, positive masses $m_i > 0$, and $\sum_i m_i = 1$. The optimal map $T_\mu = \nabla u$ for a finite maximum
$$
u(x) = \max_{1 \leq i \leq N} \{ x \cdot y_i + c_i \},
$$
where the intercepts $c_i$ are the unique choice that puts exactly the mass $m_i$ on cell $i$. Cell calculus and a unique normalized vector of intercepts are established first, then a moment estimate that controls variations of the first moments $B_i$ by the mass variations $M_i$.

pdfrefs_brenier/arxiv_1603.05579v2-kitagawa-merigot-thibert-2019.pdf

Cell calculus. Fix distinct sites $y_1, \ldots, y_N \in \mathbb{R}^d$ and intercepts $c = (c_1, \ldots, c_N)$. The closed Laguerre cells
$$
P_i(c) = K \cap \bigcap_{j \neq i} \{ x \in \mathbb{R}^d : x \cdot (y_i - y_j) + c_i - c_j \geq 0 \}
$$
partition $K$ up to the null tie set. Write $M_i(c) = \rho(P_i(c))$ and $B_i(c) = \int_{P_i(c)} x \, d\rho(x)$.

For each $i \neq j$ the equality plane is $H_{ij} = \{ x \in \mathbb{R}^d : x \cdot (y_i - y_j) + c_i - c_j = 0 \}$, and the common face is $F_{ij} = P_i \cap P_j$. Introduce the surface coefficients
$$
g_{ij} = \frac{1}{|K| |y_i - y_j|} \int_{F_{ij}} 1 \, dS, \qquad
b_{ij} = \frac{1}{|K| |y_i - y_j|} \int_{F_{ij}} x \, dS,
$$
where $dS$ denotes the $(d-1)$-dimensional surface measure on $H_{ij}$. Faces of lower dimension, or faces that meet the boundary of $K$ only on a null set, contribute zero.

On the open set of configurations with distinct sites and $M_i > 0$ for every $i$, the maps $M_i$ and $B_i$ are jointly $C^1$ in intercepts and sites. For an intercept direction $r = (r_i)$ and a site direction $v = (v_i)$,
$$
D M_i[r,v] = \sum_{j \neq i} \bigl( g_{ij} (r_i - r_j) + b_{ij} \cdot (v_i - v_j) \bigr),
$$
$$
D B_i[r,v] = \sum_{j \neq i} b_{ij} (r_i - r_j) + \frac{1}{|K|} \sum_{j \neq i} \frac{1}{|y_i - y_j|} \int_{F_{ij}} x \bigl( x \cdot (v_i - v_j) \bigr) \, dS(x).
$$
The intercept-only operator
$$
(G q)_i = D M_i[q,0] = \sum_{j \neq i} g_{ij} (q_i - q_j)
$$
is a symmetric weighted graph-Laplacian: its kernel is the constant vectors, and its range is $\{ z \in \mathbb{R}^N : \sum_i z_i = 0 \}$.

The derivative formula is obtained by telescoping the difference of indicator products one factor at a time. For a fixed pair $i, j$ write $\ell_{ij}(x) = x \cdot (y_i - y_j) + c_i - c_j$ and perturb the configuration by $h(r, v)$ with $h \downarrow 0$. On $K$, the region where the indicator $\mathbf{1}_{\ell_{ij} \geq 0}$ changes has $|\ell_{ij}(x)| \leq C h$ for a constant $C$ that depends on the configuration and direction. Coordinatizing by $(z, t) \mapsto z + h t \, n_{ij}$ with $n_{ij} = (y_i - y_j) / |y_i - y_j|$ and $z \in H_{ij}$, the volume factor $h$ cancels the difference quotient. Outside the null set $(\partial K \cap H_{ij}) \cup \bigcup_{k \neq i, j} (H_{ik} \cap H_{ij})$, all remaining indicators converge to the values at $z$, and the changed factor converges to $\mathbf{1}_{\{ |y_i - y_j| t + r_i - r_j + z \cdot (v_i - v_j) \geq 0 \}} - \mathbf{1}_{\{ t \geq 0 \}}$. Its integral in $t$ is $(r_i - r_j + z \cdot (v_i - v_j)) / |y_i - y_j|$. Dominated convergence yields the formula. Applying the same argument in the reverse direction gives two-sided derivatives, and continuity of the coefficients in $(c, y)$ follows from a similar uniform convergence of the face integrals and the Jacobian factors.

Connectivity of the graph on $\{ 1, \ldots, N \}$ with edges $\{ g_{ij} > 0 \}$ follows from a segment argument. Choose strictly active points $p, q$ in two cells with $p$ outside all equality planes and $q$ outside all equality planes and the affine spans of $p$ with every triple-equality set. The segment $[p, q]$ stays in the interior of $K$ and does not meet any triple equality set. Along the segment the maximizing index changes only by crossing equality planes one at a time. At each crossing exactly two affine functions are maximal; a small relatively open disk in the common equality plane lies inside the common face, so the corresponding $g_{ij} > 0$. The successive maximizing indices give a graph path. Symmetry of the face coefficients gives $q \cdot G q = \sum_{i < j} g_{ij} (q_i - q_j)^2$, which identifies the kernel of $G$ with the constant vectors and its range with the orthogonal complement of the constant vectors.

Prescribed cell masses. Fix positive weights $m_1, \ldots, m_N$ with $\sum_i m_i = 1$. For every distinct site configuration there is a unique intercept vector $c$ with
$$
M_i(c) = m_i \quad (1 \leq i \leq N), \qquad \sum_{i=1}^N m_i c_i = 0,
$$
and this normalized vector depends $C^1$ on the sites.

The proof minimizes the convex function
$$
\mathcal{F}_y(c) = \int_K \max_i \{ x \cdot y_i + c_i \} \, d\rho(x) - \sum_{i=1}^N m_i c_i
$$
on the hyperplane $H_m = \{ c \in \mathbb{R}^N : \sum_i m_i c_i = 0 \}$. On $H_m$, an upper bound on $\max_i c_i$ together with the weighted zero-mean condition bounds every coordinate from below, so $\mathcal{F}_y(c) \to \infty$ as $|c| \to \infty$; coercivity provides a minimizer. Distinct sites make the maximizing index unique almost everywhere, and differentiating the maximum under the integral (justified by bounded difference quotients) gives $\nabla \mathcal{F}_y(c) = (M_i(c) - m_i)_{i=1}^N$. At a minimum on $H_m$ this gradient lies in the span of $m$. Its coordinate sum is zero, since the masses and weights each sum to one. As $\sum_i m_i = 1$, the gradient must vanish, so $M_i(c) = m_i$. Conversely, any such $c$ minimizes $\mathcal{F}_y$ by convexity.

Local uniqueness as the sites vary follows from the implicit function theorem applied to the mass residual on $H_m$. At any solution all cells have positive mass, so the calculus above applies and the intercept derivative of the mass residual is the restriction
$$
G : H_m \longrightarrow \left\{ z \in \mathbb{R}^N : \sum_i z_i = 0 \right\}.
$$
Both spaces have dimension $N-1$ and the only constant vector in $H_m$ is zero, so the restriction is an isomorphism. The implicit function theorem gives a locally unique $C^1$ branch. Global uniqueness follows: two normalized solutions would produce a segment of minimizers of $\mathcal{F}_y$ by convexity, contradicting local isolation. The $C^1$ branches glue by uniqueness. Only invertibility of the restricted $G$ is used; no quantitative bound on its inverse is needed.

pdfrefs_brenier/arxiv_1603.05579v2-kitagawa-merigot-thibert-2019.pdf

A geometric coupling lemma. Let $U, V \subset \mathbb{R}^d$ be compact convex sets of positive volume. There is a coupling $(X, Z)$ of the uniform probability measures on $U$ and $V$ with
$$
\mathrm{Law}\!\left( \frac{X + Z}{2} \right) \leq \frac{dx}{\sqrt{|U| |V|}}.
$$
The inequality is an inequality of measures on $\mathbb{R}^d$.

The construction is Knothe's triangular rearrangement, adapted to averaging the map with its counterpart. Alternative derivations come from displacement concavity and from the arithmetic-geometric mean inequality on the diagonal entries of the triangular map.
pdfrefs_brenier/carlier-galichon-santambrogio-2010-siam-j-math-anal.pdf
pdfrefs_brenier/mccann-1997-advances.pdf

For $0 \leq k < d$ set
$$
p_k(\xi) = \int_{\mathbb{R}^{d-k}} 1_{\mathrm{int}(U)}(\xi, \eta) \, d\eta, \qquad p_d = 1_{\mathrm{int}(U)}.
$$
Thus $p_0 = |U|$. Let $D_k$ be the projection of $\mathrm{int}(U)$ onto its first $k$ coordinates. Each $p_k$ is positive and continuous on $D_k$. Continuity follows by dominated convergence: at a prefix in $D_k$, the sublevel sets of the fibers converge almost everywhere as the prefix varies, because a hyperplane meeting the interior of a convex body meets it in a face whose outer boundary is null; the boundedness of $U$ provides a common integrable bound. Fiber integrals are Borel functions on the full domains.

Given $s = (s_1, \ldots, s_d) \in (0,1)^d$, construct $X(s)$ recursively. Given a prefix $\xi = (X_1, \ldots, X_{k-1}) \in D_{k-1}$, let $X_k$ be the quantile at level $s_k$ of the density
$$
t \mapsto \frac{p_k(\xi, t)}{p_{k-1}(\xi)}.
$$
Fubini makes this a probability density on a bounded open interval, and its cumulative distribution function is strictly increasing there. The quantile is therefore unique, keeps the new prefix in $D_k$, and satisfies
$$
\partial_{s_k} X_k = \frac{p_{k-1}(X_1, \ldots, X_{k-1})}{p_k(X_1, \ldots, X_k)} > 0.
$$
The construction is measurable: the conditional cumulative distribution functions are Borel in the arguments, and quantiles of Borel functions are Borel. With $s$ uniform on $(0, 1)^d$, successive conditioning gives the uniform law on $U$. The diagonal partial derivatives telescope:
$$
\prod_{k=1}^d \partial_{s_k} X_k = \frac{p_0}{p_d} = |U|,
$$
which is constant on the interior of $U$.

The same construction for $V$ with the same parameter $s$ yields $Z(s)$ and $\prod_{k=1}^d \partial_{s_k} Z_k = |V|$. For $H = (X+Z)/2$, the arithmetic-geometric mean inequality applied to the diagonal derivatives yields
$$
\prod_{k=1}^d \partial_{s_k} H_k
\geq \sqrt{\left( \prod_{k=1}^d \partial_{s_k} X_k \right) \left( \prod_{k=1}^d \partial_{s_k} Z_k \right)}
= \sqrt{|U| |V|}.
$$

Now apply a change of variables to $H$. For a nonnegative Borel $g$ on $\mathbb{R}^d$ and fixed $a = (s_1, \ldots, s_{d-1})$, one-dimensional substitution in the strictly increasing $C^1$ map $s_d \mapsto H_d(a, s_d)$, followed by enlargement of its image to all of $\mathbb{R}$, gives
$$
\int_0^1 g\bigl(H_1(a), \ldots, H_{d-1}(a), H_d(a,t)\bigr) \partial_{s_d} H_d(a,t) \, dt \leq \int_{\mathbb{R}} g\bigl(H_1(a), \ldots, H_{d-1}(a), z\bigr) \, dz.
$$
The last integral is itself a nonnegative Borel function of $a$. Multiplying by the earlier derivative factors (which depend only on $a$) and applying Tonelli's theorem, then iterating this argument backwards over $k = d-1, \ldots, 1$, gives
$$
\int_{(0,1)^d} g(H(s)) \prod_{k=1}^d \partial_{s_k} H_k(s) \, ds \leq \int_{\mathbb{R}^d} g(x) \, dx.
$$
Combining with the pointwise lower bound on the product gives
$$
\mathbb{E} g\!\left( \frac{X + Z}{2} \right) \leq \frac{1}{\sqrt{|U| |V|}} \int_{\mathbb{R}^d} g(x) \, dx,
$$
which is the asserted measure domination.

The moment estimate. Fix distinct sites and intercepts with $M_i > 0$ for every $i$, and let $K \subset B(0, R)$ for some $R > 0$. Then every intercept direction $q \in \mathbb{R}^N$ satisfies
$$
\sum_{i=1}^N \frac{|D B_i[q, 0]|^2}{M_i} \leq C_0 \, R^2 \sum_{i=1}^N \frac{|D M_i[q, 0]|^2}{M_i},
\qquad C_0 = 162.
$$
The estimate is uniform over the number of cells, the cell masses, and the site separations.

Fix the direction $q$ and hold the sites fixed. For small $s > 0$ let $P_i^-$ and $P_i^+$ be the cells at $c - s q$ and $c + s q$. Set
$$
\alpha_i = \rho(P_i^-), \qquad \beta_i = \rho(P_i^+), \qquad \gamma_i = \sqrt{\alpha_i \beta_i}.
$$
Continuity in the cell-calculus lemma makes every $\alpha_i, \beta_i$ positive for small $s$. Couple the uniform measures on $P_i^{\pm}$ by the geometric coupling lemma, writing $X_i \in P_i^-$, $Z_i \in P_i^+$ and $D_i = \mathbb{E} |X_i - Z_i|^2$.

For any competing index $j$, the definition of $P_i^{\mp}$ gives
$$
X_i \cdot (y_i - y_j) + c_i - c_j - s (q_i - q_j) \geq 0, \qquad
Z_i \cdot (y_i - y_j) + c_i - c_j + s (q_i - q_j) \geq 0.
$$
Averaging,
$$
\frac{X_i + Z_i}{2} \cdot (y_i - y_j) + c_i - c_j \geq 0.
$$
Convexity of $K$ keeps the midpoint in $K$, so the midpoint lies in $P_i$. Since $|P_i^{\pm}| = |K| \alpha_i^{\pm}$, the coupling lemma applied to $U = P_i^-$ and $V = P_i^+$ gives the measure domination
$$
\gamma_i \, \mathrm{Law}\!\left( \frac{X_i + Z_i}{2} \right) \leq \rho\big|_{P_i}.
$$

Apply this domination with the nonnegative test function $\Phi(x) = R^2 - |x|^2$. Set $A_i^{\pm} = \mathbb{E} \Phi(X_i^{\pm})$, where $X_i^-$ and $X_i^+$ denote the two uniform samples. The algebraic identity
$$
\Phi\!\left( \frac{X_i + Z_i}{2} \right)
= \frac{\Phi(X_i) + \Phi(Z_i)}{2} + \frac{|X_i - Z_i|^2}{4}
$$
together with the measure domination gives
$$
\sum_i \gamma_i \left( \frac{A_i^- + A_i^+}{2} + \frac{D_i}{4} \right)
\leq \int_K \Phi \, d\rho
= \frac{1}{2} \sum_i \left( \alpha_i A_i^- + \beta_i A_i^+ \right).
$$
The rightmost equality uses that the two partitions at $c \pm s q$ have null overlaps.

Set $\delta_i = \alpha_i - \beta_i$ and
$$
\mathcal{A} = \sum_i \gamma_i D_i, \qquad
\mathcal{B} = \sum_i \frac{\delta_i^2}{\gamma_i}.
$$
Subtracting the averaged endpoint terms in the display above yields
$$
\frac{\mathcal{A}}{4}
\leq \sum_i \left[ \left( \frac{\alpha_i + \beta_i}{2} - \gamma_i \right) \frac{A_i^- + A_i^+}{2} + \frac{\delta_i}{4} (A_i^- - A_i^+) \right].
$$
Note that $0 \leq A_i^{\pm} \leq R^2$, and
$$
\frac{\alpha_i + \beta_i}{2} - \gamma_i
= \frac{(\alpha_i - \beta_i)^2}{2 (\sqrt{\alpha_i} + \sqrt{\beta_i})^2}
\leq \frac{\delta_i^2}{8 \gamma_i}.
$$
Since $|X_i^{\pm}| \leq R$,
$$
|A_i^- - A_i^+| \leq 2 R \, \mathbb{E} |X_i^- - X_i^+| \leq 2 R \sqrt{D_i}.
$$
Cauchy-Schwarz therefore gives
$$
\frac{\mathcal{A}}{4} \leq \frac{R^2 \mathcal{B}}{8} + \frac{R}{2} \sqrt{\mathcal{A} \mathcal{B}}.
$$
Using $\frac{R}{2} \sqrt{\mathcal{A} \mathcal{B}} \leq \frac{\mathcal{A}}{8} + \frac{R^2 \mathcal{B}}{2}$ and rearranging,
$$
\mathcal{A} \leq 5 R^2 \mathcal{B}.
$$

The first-moment difference over the symmetric change is
$$
\Delta B_i = B_i(c - s q) - B_i(c + s q) = \alpha_i \mathbb{E} X_i - \beta_i \mathbb{E} Z_i.
$$
Hence
$$
|\Delta B_i| \leq R |\delta_i| + \alpha_i \sqrt{D_i}.
$$
For the fixed configuration and the fixed direction $q$, continuity allows $s$ small enough that $\alpha_i, \beta_i, \gamma_i$ all lie in $[M_i/2, 2 M_i]$. Set $\mathcal{E} = \sum_i \delta_i^2 / M_i$. These comparisons give $\alpha_i^2 / M_i \leq 8 \gamma_i$ and $\mathcal{B} \leq 2 \mathcal{E}$. Consequently
$$
\sum_i \frac{|\Delta B_i|^2}{M_i}
\leq 2 R^2 \mathcal{E} + 2 \sum_i \frac{\alpha_i^2}{M_i} D_i
\leq 2 R^2 \mathcal{E} + 16 \mathcal{A}
\leq 162 R^2 \mathcal{E}.
$$
Divide by $(2 s)^2$ and let $s \downarrow 0$. The cell-calculus lemma gives
$$
\frac{\Delta B_i}{2 s} \to -D B_i[q, 0], \qquad
\frac{\delta_i}{2 s} = \frac{M_i(c - s q) - M_i(c + s q)}{2 s} \to -D M_i[q, 0],
$$
so the claimed estimate follows.

The threshold for $s$ may depend on the direction $q$ and on the configuration, including its smallest mass. Only the constant $162$ survives the derivative limit.
pdfrefs_brenier/arxiv_1603.05579v2-kitagawa-merigot-thibert-2019.pdf

From the moment estimate to potential stability along an optimal coupling. Let $\mu, \nu \in \mathrm{Prob}(Y)$ be finitely supported. Choose an optimal coupling
$$
\pi = \sum_{i=1}^N m_i \delta_{(y_i^0, y_i^1)}, \qquad m_i > 0, \quad \sum_i m_i = 1,
$$
after merging identical pairs and discarding zero-mass ones. Then
$$
W_2(\mu, \nu)^2 = \sum_{i=1}^N m_i |y_i^0 - y_i^1|^2.
$$
Write $v_i = y_i^1 - y_i^0$; the quantity $\sum_i m_i |v_i|^2 = W_2(\mu, \nu)^2$.

Optimality implies the monotonicity inequality
$$
(y_i^0 - y_j^0) \cdot (y_i^1 - y_j^1) \geq 0 \qquad (1 \leq i, j \leq N).
$$
If the scalar product is negative, swapping equal positive portions of the mass between pairs $i$ and $j$ reduces the total quadratic cost by $2 \varepsilon (y_i^0 - y_j^0) \cdot (y_i^1 - y_j^1)$ for $0 < \varepsilon \leq \min(m_i, m_j)$.

Set $y_i(t) = (1-t) y_i^0 + t y_i^1$. For $0 < t < 1$, the sites $y_1(t), \ldots, y_N(t)$ are distinct. Indeed, with $a = y_i^0 - y_j^0$ and $b = y_i^1 - y_j^1$,
$$
|y_i(t) - y_j(t)|^2 = (1-t)^2 |a|^2 + t^2 |b|^2 + 2 t (1-t) a \cdot b,
$$
and the right-hand side is strictly positive unless $a = b = 0$, in which case the pairs $(y_i^0, y_i^1)$ and $(y_j^0, y_j^1)$ coincide and were already merged.

By the prescribed-mass lemma, for $0 < t < 1$ there are unique intercepts $c_i(t)$ with $M_i(t) = m_i$ and $\sum_i m_i c_i(t) = 0$, and these are $C^1$ in $t$. Let $r_i = c_i'(t)$; then $\sum_i m_i r_i = 0$.

Differentiating the mass constraints gives $D M_i[c'(t), y'(t)] = 0$ for every $i$. Combining the two derivative formulas from cell calculus, for every intercept direction $q$,
$$
\sum_i q_i D M_i[r, v] = \sum_i r_i D M_i[q, 0] + \sum_i v_i \cdot D B_i[q, 0].
$$
Indeed, expanding
$$
\sum_i q_i D M_i[r, v] = \sum_{i < j} g_{ij} (q_i - q_j)(r_i - r_j) + \sum_{i < j} b_{ij} \cdot (q_i - q_j)(v_i - v_j),
$$
$$
\sum_i v_i \cdot D B_i[q, 0] = \sum_{i < j} b_{ij} \cdot (q_i - q_j)(v_i - v_j),
$$
and the difference of the first expression and the second is $\sum_{i < j} g_{ij} (q_i - q_j)(r_i - r_j) = \sum_i r_i D M_i[q, 0]$ by symmetry of $g_{ij}$. Setting $D M_i[r, v] = 0$ gives the duality relation
$$
\sum_i r_i D M_i[q, 0] = - \sum_i v_i \cdot D B_i[q, 0].
$$

The range of $G$ is $\{ z \in \mathbb{R}^N : \sum_i z_i = 0 \}$. Together with $\sum_i m_i r_i = 0$, this permits a choice of intercept direction $q$ with
$$
D M_i[q, 0] = m_i r_i \qquad (1 \leq i \leq N).
$$
No quantitative bound on $q$ is needed.

Set $E_r = \sum_i m_i |r_i|^2$. Using the duality relation, weighted Cauchy-Schwarz, the equality $\sum_i m_i |v_i|^2 = W_2(\mu, \nu)^2$, and the moment estimate,
$$
E_r = \sum_i r_i (m_i r_i) = \sum_i r_i D M_i[q, 0] = - \sum_i v_i \cdot D B_i[q, 0]
$$
$$
\leq \left( \sum_i m_i |v_i|^2 \right)^{1/2} \left( \sum_i \frac{|D B_i[q, 0]|^2}{m_i} \right)^{1/2}
\leq W_2(\mu, \nu) \left( C_0 R^2 \sum_i \frac{|D M_i[q, 0]|^2}{m_i} \right)^{1/2}
$$
$$
= W_2(\mu, \nu) \left( C_0 R^2 \sum_i m_i |r_i|^2 \right)^{1/2}
= \sqrt{C_0} R \, W_2(\mu, \nu) \, \sqrt{E_r}.
$$
This includes the case $E_r = 0$. Hence
$$
\left( \sum_i m_i |r_i|^2 \right)^{1/2} \leq \sqrt{C_0} R \, W_2(\mu, \nu).
$$

Define the interpolating potential
$$
u_t(x) = \max_i \{ x \cdot y_i(t) + c_i(t) \}.
$$
On $K$, the slopes $y_i(t)$ stay in $B(0, L)$ (a convex combination of two points in the ball), so $u_t$ is convex and $L$-Lipschitz. At a strict maximum index $i$, the time derivative is $\partial_t u_t(x) = x \cdot v_i + r_i$. Fubini identifies this derivative almost everywhere in $(t, x) \in (0,1) \times K$, since the tie set has zero $\rho$-measure for each fixed $t$. Consequently, for interior times,
$$
\|\partial_t u_t\|_{L^2(\rho)} \leq R \left( \sum_i m_i |v_i|^2 \right)^{1/2} + \left( \sum_i m_i |r_i|^2 \right)^{1/2} \leq (1 + \sqrt{C_0}) R \, W_2(\mu, \nu).
$$
The fundamental theorem of calculus together with Minkowski's inequality gives, for $0 < t < t' < 1$,
$$
\|u_t - u_{t'}\|_{L^2(\rho)} \leq (1 + \sqrt{C_0}) R \, W_2(\mu, \nu) |t - t'|.
$$

Extend to endpoints. For each $i$ the estimate $|c_i'(t)| \leq \sqrt{C_0} R W_2(\mu, \nu) / \sqrt{m_i}$ shows that $c_i(t)$ has finite limits as $t \downarrow 0$ and $t \uparrow 1$. Denote them again by $c_i(0)$ and $c_i(1)$, and set
$$
u_0(x) = \max_i \{ x \cdot y_i^0 + c_i(0) \}, \qquad u_1(x) = \max_i \{ x \cdot y_i^1 + c_i(1) \}.
$$
Because $y_i(t) \to y_i^0$ and $c_i(t) \to c_i(0)$ as $t \downarrow 0$, the potentials $u_t$ converge uniformly to $u_0$ on $K$; similarly $u_t \to u_1$ as $t \uparrow 1$. The stability inequality therefore holds at the endpoints:
$$
\|u_0 - u_1\|_{L^2(\rho)} \leq (1 + \sqrt{C_0}) R \, W_2(\mu, \nu).
$$

Finally, $u_0$ and $u_1$ are the finite-maximum potentials of $\mu$ and $\nu$. At each interior $t$, the slope map of $u_t$ produces exactly the law $\sum_i m_i \delta_{y_i(t)}$, because $M_i(t) = m_i$. Passing to the limit, the slope law of $u_0$ is $\mu$ and that of $u_1$ is $\nu$: for any continuous test $\varphi$,
$$
\int_K \varphi(\nabla u_t) \, d\rho = \sum_i m_i \varphi(y_i(t)) \xrightarrow[t \downarrow 0]{} \sum_i m_i \varphi(y_i^0) = \int_Y \varphi \, d\mu.
$$
By the optimality criterion, the slope selection of each finite maximum with a given slope law is that law's unique optimal map from $\rho$. Hence $T_\mu = \nabla u_0$ and $T_\nu = \nabla u_1$ almost everywhere.
pdfrefs_brenier/arxiv_2103.05934v2-delalande-merigot-2023.pdf
pdfrefs_brenier/arxiv_1603.05579v2-kitagawa-merigot-thibert-2019.pdf

Gradient interpolation. Let $u_0, u_1$ be convex and $L$-Lipschitz on $K$, and set $f = u_0 - u_1$. For each coordinate axis $e_j$, let $A_j$ be the $(d-1)$-dimensional volume of the projection of $K$ onto the hyperplane orthogonal to $e_j$, and set $P_K = \sum_{j=1}^d A_j$. Then for every $h > 0$,
$$
\|\nabla u_0 - \nabla u_1\|_{L^2(\rho)}^2 \leq \frac{12 d}{h^2} \|f\|_{L^2(\rho)}^2 + \frac{28 L^2 P_K}{|K|} h.
$$
Only the almost-everywhere coordinate derivatives of $u_0$ and $u_1$ are needed.

This is the coordinate-slice version of the convex-function interpolation inequality of Delalande and Merigot, with explicit constants.
pdfrefs_brenier/arxiv_2103.05934v2-delalande-merigot-2023.pdf

Fix $j$ and put
$$
K_{j,h} = \{ x \in K : x + h e_j \in K \}, \qquad D_{j,h} u(x) = \frac{u(x + h e_j) - u(x)}{h} \quad (x \in K_{j,h}).
$$
Almost every nonempty coordinate slice of $K$ (a line parallel to $e_j$) is an interval $[\alpha, \beta]$ on which the restriction of each $u_i$ is convex and $L$-Lipschitz. Its derivative exists almost everywhere, is monotone, and has absolute value at most $L$. Fubini then gives the existence of all coordinate derivatives almost everywhere in $K$.

At a point where both $D_{j,h} u$ and $\partial_j u$ exist, convexity gives
$$
0 \leq \delta_u(x) := D_{j,h} u(x) - \partial_j u(x) \leq \partial_j u(x + h e_j) - \partial_j u(x) \leq 2 L.
$$
The bounds come from the ordering of secant slopes between the one-sided derivatives at $x$ and at $x + h e_j$, and from the uniform $L$-Lipschitz bound on the derivative.

For a slice $[\alpha, \beta]$ of length at least $h$, cancellation of translated integrals gives
$$
\int_\alpha^{\beta - h} \bigl( \partial_j u(t + h) - \partial_j u(t) \bigr) \, dt
= \int_{\beta - h}^{\beta} \partial_j u(t) \, dt - \int_\alpha^{\alpha + h} \partial_j u(t) \, dt \leq 2 L h.
$$
Each boundary term is bounded by $L h$. Coordinates transverse to the slice are suppressed. For a slice of length less than $h$, its allowed portion is empty. Using $\delta_u^2 \leq 2 L \, \delta_u$ and integrating over slices, then over the orthogonal projection,
$$
\int_{K_{j,h}} \delta_u^2 \, d\rho \leq \frac{4 L^2 h A_j}{|K|}.
$$
The excluded portion of any slice has length at most $h$, so
$$
\rho(K \setminus K_{j,h}) \leq \frac{h A_j}{|K|}.
$$

On the allowed set where both $x$ and $x + h e_j$ lie in $K$, translation invariance of Lebesgue measure and $(a - b)^2 \leq 2 a^2 + 2 b^2$ give
$$
\int_{K_{j,h}} |D_{j,h} f|^2 \, d\rho \leq \frac{4}{h^2} \|f\|_{L^2(\rho)}^2.
$$
On the boundary strip $K \setminus K_{j,h}$, the Lipschitz bound gives $|\partial_j f| \leq 2 L$, so
$$
\int_{K \setminus K_{j,h}} |\partial_j f|^2 \, d\rho \leq 4 L^2 \rho(K \setminus K_{j,h}) \leq \frac{4 L^2 h A_j}{|K|}.
$$
On $K_{j,h}$ the pointwise identity $\partial_j f = D_{j,h} f - \delta_{u_0} + \delta_{u_1}$ together with $(a + b + c)^2 \leq 3 (a^2 + b^2 + c^2)$ yields
$$
\int_{K_{j,h}} |\partial_j f|^2 \, d\rho \leq \frac{12}{h^2} \|f\|_{L^2(\rho)}^2 + \frac{24 L^2 h A_j}{|K|}.
$$
Summing with the boundary strip contribution gives
$$
\int_K |\partial_j f|^2 \, d\rho \leq \frac{12}{h^2} \|f\|_{L^2(\rho)}^2 + \frac{28 L^2 h A_j}{|K|}.
$$
Summing over $j$ and using $\sum_j A_j = P_K$ proves the claimed inequality. All projection volumes $A_j$ are finite because $K$ is compact, and the argument applies to any $h > 0$, including values exceeding some or all slice lengths.

The coordinate derivatives form the full gradient almost everywhere. At an interior point $x$ where every coordinate partial of a convex function $u$ exists, set $p = (\partial_j u(x))_j$ and $q = \|v\|_1$ for a small nonzero vector $v$. The point $x + v$ is a convex combination of $x + \mathrm{sgn}(v_j) q e_j$ with weights $|v_j| / q$ (omitting zero weights). The one-sided derivative limits and convexity give $u(x + v) - u(x) \leq p \cdot v + o(q)$. Applying the same bound to $-v$ and using midpoint convexity gives the reverse, so $u$ is differentiable at $x$ with gradient $p$. Since $\partial K$ has zero volume, this applies almost everywhere to both $u_0$ and $u_1$.

The finite-case bound at exponent one-third. Set $C_0 = 162$ and
$$
C_* = \left( 12 d R^2 \left( 1 + \sqrt{C_0} \right)^2 + \frac{28 L^2 P_K}{|K|} \right)^{1/2}.
$$
Then, for any two finitely supported probabilities $\mu, \nu \in \mathrm{Prob}(Y)$,
$$
\|T_\mu - T_\nu\|_{L^2(\rho)} \leq C_* \, W_2(\mu, \nu)^{1/3}.
$$

The proof is a direct substitution. Write $w = W_2(\mu, \nu)$. The finite-potential estimate provides finite-maximum potentials $u_0, u_1$ with $\nabla u_0 = T_\mu$ and $\nabla u_1 = T_\nu$ almost everywhere and
$$
\|u_0 - u_1\|_{L^2(\rho)} \leq (1 + \sqrt{C_0}) R w.
$$
If $w > 0$, apply gradient interpolation with $h = w^{2/3}$:
$$
\|\nabla u_0 - \nabla u_1\|_{L^2(\rho)}^2
\leq \frac{12 d}{w^{4/3}} (1 + \sqrt{C_0})^2 R^2 w^2 + \frac{28 L^2 P_K}{|K|} w^{2/3}
= C_*^2 w^{2/3}.
$$
Taking the square root completes the finite case. If $w = 0$, then $\mu = \nu$ and $u_0 = u_1$ almost everywhere, so letting $h \downarrow 0$ in the interpolation inequality gives $\nabla u_0 = \nabla u_1$ almost everywhere and the bound $C_* w^{1/3} = 0$ holds.
pdfrefs_brenier/arxiv_2103.05934v2-delalande-merigot-2023.pdf

Extension to all Borel targets. For every $\lambda \in \mathrm{Prob}(Y)$ there is a unique quadratic optimal map $T_\lambda$ from $\rho$ to $\lambda$, up to a $\rho$-null set. If $F_n \subset Y$ are finite $\epsilon_n$-nets with $\epsilon_n \to 0$ and $Q_n : Y \to F_n$ are Borel nearest-point selections, then
$$
T_{(Q_n)_\#\lambda} \longrightarrow T_\lambda \quad \text{in } L^2(\rho).
$$

Compactness of $Y$ provides the finite nets. Enumerate each $F_n$ and choose the first closest point, obtaining a Borel selection $Q_n$ with $|Q_n(y) - y| \leq \epsilon_n$. Put $\lambda_n = (Q_n)_\#\lambda$, discard any zero-mass sites, and let $T_n = T_{\lambda_n}$. These are optimal finite-target maps by the cell calculus. The coupling induced by $y \mapsto (Q_n(y), Q_m(y))$ gives
$$
W_2(\lambda_n, \lambda_m) \leq \epsilon_n + \epsilon_m.
$$
By finite-target stability, $(T_n)$ is Cauchy in $L^2(\rho)$. Let $S$ be the $L^2(\rho)$ limit and choose a subsequence converging almost everywhere. Since $Y$ is closed, the limit lies in $Y$ almost everywhere; modifying $S$ on a null set makes it Borel with values in $Y$. For any continuous $\varphi : Y \to \mathbb{R}$, bounded convergence and uniform continuity give
$$
\int_K \varphi(S(x)) \, d\rho(x) = \lim_n \int_Y \varphi(Q_n(y)) \, d\lambda(y) = \int_Y \varphi(y) \, d\lambda(y),
$$
so $S_\#\rho = \lambda$.

Direct optimality. Choose finite maximum potentials $u_n$ for $T_n$ and subtract a constant so $u_n(x_*) = 0$ at a fixed $x_* \in K$. Slopes lie in $Y$, so each $u_n$ is $L$-Lipschitz and $|u_n| \leq L \, \mathrm{diam}(K)$ on $K$. Arzela-Ascoli gives a further subsequence converging uniformly on $K$ to a convex, $L$-Lipschitz function $u$. Intersect the full-measure set on which each of the countably many supporting inequalities holds with the set where the chosen subsequence of $T_n$ converges; this intersection has full $\rho$-measure. On the intersection, since each $u_n$ is a finite maximum with slopes in $Y$,
$$
u_n(z) \geq u_n(x) + T_n(x) \cdot (z - x) \qquad (z \in K).
$$
Passing to the limit gives, simultaneously for every $z \in K$,
$$
u(z) \geq u(x) + S(x) \cdot (z - x) \quad \text{for } \rho\text{-almost every } x.
$$
Define the restricted convex conjugate $u^*(y) = \max_{z \in K} (z \cdot y - u(z))$. It is finite and continuous, and bounded on $Y$. The same two displays give
$$
x \cdot y \leq u(x) + u^*(y) \quad (x \in K, \; y \in Y), \qquad
x \cdot S(x) = u(x) + u^*(S(x)) \quad \rho\text{-almost everywhere}.
$$

For any competing map $V$ with $V_\#\rho = \lambda$, the integral of $u(x) + u^*(V(x))$ d$\rho$ is fixed by the marginals. The first inequality of the display shows that $S$ maximizes $\int_K x \cdot V(x) \, d\rho(x)$. Combined with the cost identity
$$
\int_K |x - V(x)|^2 \, d\rho(x) = \int_K |x|^2 \, d\rho(x) + \int_Y |y|^2 \, d\lambda(y) - 2 \int_K x \cdot V(x) \, d\rho(x),
$$
where the first and second terms are fixed by the marginals, this identifies $S$ as an optimal map.

Uniqueness. If $V$ is another optimal map for $\lambda$, the nonnegative Fenchel gap $u(x) + u^*(V(x)) - x \cdot V(x)$ has zero integral, hence vanishes almost everywhere. Consequently $u(z) \geq u(x) + V(x) \cdot (z - x)$ for every $z \in K$ at almost every $x$. On almost every interior point of $K$ the coordinate partials of $u$ exist (Fubini and convexity of $u$ along each coordinate axis); at such $x$ a supporting slope $p$ must satisfy
$$
p_j \leq \lim_{t \downarrow 0} \frac{u(x + t e_j) - u(x)}{t}, \qquad p_j \geq \lim_{t \uparrow 0} \frac{u(x + t e_j) - u(x)}{t}.
$$
Both limits equal $\partial_j u(x)$, so $V(x) = \nabla u(x) = S(x)$. The boundary $\partial K$ has zero volume, so $V = S$ almost everywhere.

The same integrated inequality also handles arbitrary transport plans: equality in the Fenchel inequality forces the target coordinate to be the unique supporting slope over almost every source point, so the optimal plan is the one induced by $S$. Hence, for every $\lambda \in \mathrm{Prob}(Y)$, there is a unique quadratic optimal map $T_\lambda$, and $T_{(Q_n)_\#\lambda} \to T_\lambda$ in $L^2(\rho)$.

Assembly. Let $\mu, \nu \in \mathrm{Prob}(Y)$ be arbitrary, and use the same finite nets and selections as above. Put $\mu_n = (Q_n)_\#\mu$ and $\nu_n = (Q_n)_\#\nu$. For any coupling $\pi \in \Pi(\mu, \nu)$, its image under $(Q_n, Q_n)$ is a coupling of $\mu_n$ and $\nu_n$, and the triangle inequality in $L^2(\pi)$ gives
$$
\int |y - z|^2 \, d((Q_n, Q_n)_\#\pi)(y, z) \leq \left( \left( \int |y - z|^2 \, d\pi(y,z) \right)^{1/2} + 2 \epsilon_n \right)^2.
$$
Taking the infimum over $\pi$ yields $W_2(\mu_n, \nu_n) \leq W_2(\mu, \nu) + 2 \epsilon_n$. Applying finite-target stability to $\mu_n$ and $\nu_n$,
$$
\|T_{\mu_n} - T_{\nu_n}\|_{L^2(\rho)} \leq C_* \left( W_2(\mu, \nu) + 2 \epsilon_n \right)^{1/3}.
$$
Since $T_{\mu_n} \to T_\mu$ and $T_{\nu_n} \to T_\nu$ in $L^2(\rho)$, passing to the limit gives
$$
\|T_\mu - T_\nu\|_{L^2(\rho)} \leq C_* \, W_2(\mu, \nu)^{1/3}.
$$

The constant $C_*$ depends only on $K$ and the enclosing radius $L$ for $Y$. It is independent of atom counts, masses, site separations, and density bounds. No convexity of $Y$ is required: the site interpolation stays in $B(0, L)$ and the discretization takes values in $Y$. The construction covers the extreme case of a singleton $Y$, for which every target map coincides with the constant map.

The finite-target bound, the limiting extension, and the sharpness example together establish the two claims stated at the outset: the exponent $1/3$ is admissible with a uniform constant, and no larger exponent is. This closes the gap between the $W_1^{1/4}$ estimate that Merigot records for the nonuniform $L^1$ setting and the $W_2^{1/3}$ optimal rate in the quadratic setting.
pdfrefs_brenier/arxiv_2404.02855v1-divol-nilesweed-pooladian-2025.pdf
pdfrefs_brenier/merigot-delalande-chazal-2020-pmlr-v108.pdf
pdfrefs_brenier/merigot-2026-hal.pdf