Proof. The first-index choice is Borel because its level sets are given by finitely many strict and non-strict affine inequalities. Merge repeated slopes, retaining the largest intercept for each slope, and write the resulting intercept as \(c(y)\). Ties between distinct slopes lie in a finite union of hyperplanes. Outside this null set, the maximizing affine function remains strictly maximal in a neighborhood, so its slope is the gradient of \(u\).

For every transport plan \(\pi\) from \(\rho\) to \(\lambda\), the inequality \(x\cdot y+c(y)\le u(x)\) gives \[\int x\cdot y\,d\pi(x,y) \le \int_K u\,d\rho-\int c(y)\,d\lambda(y).\] The graph plan of \(T\) attains equality. Since the two marginals fix the integrals of \(|x|^2\) and \(|y|^2\), maximizing the displayed scalar product is equivalent to minimizing quadratic cost. Equality for any optimal plan requires it to be concentrated on pairs \((x,y)\) with \(x\cdot y+c(y)=u(x)\). The maximizing slope is unique for \(\rho\)-almost every \(x\), so the plan is concentrated on the graph of \(T\). This proves uniqueness of the plan. ◻

Proof. For \(0<a<1/2\), put \(b=a/2\) and define \[ u_a(x)=\max\{-x_1,x_1,a\},\qquad \widetilde u_a(x)=\max\{-x_1,x_1,a+bx_2\}. \tag{1}\] Let \(\mu_a=(\nabla u_a)_\#\rho\) and \(\nu_a=(\nabla\widetilde u_a)_\#\rho\), with any fixed first-index selection on the null tie sets. Their slopes are, respectively, \(-e_1,e_1,0\) and \(-e_1,e_1,be_2\), all in the fixed set \(Y\). By 2, these gradients are the unique optimal maps to their laws.

Since \(a/2\le a+bx_2\le3a/2<1\), the central cells are, up to null boundaries, \[A_a=\{x\in K:|x_1|<a\},\qquad B_a=\{x\in K:|x_1|<a+bx_2\}.\] Tilting the central strip preserves its mass but changes which points receive a central or an outer slope. 1 shows the symmetric difference where this reassignment occurs.

On \(A_a\), the first map is zero; on its complement it is \(\operatorname{sign}(x_1)e_1\). The second map equals \(be_2\) on \(B_a\) and \(\operatorname{sign}(x_1)e_1\) on its complement. Under \(\rho\), the coordinates are independent and uniform on \([-1,1]\). On the slice \(x_2=t\), the second central interval has probability \(a+bt\), and each outer interval has probability \((1-a-bt)/2\). Averaging over \(t\) gives \[ \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_{be_2}. \tag{2}\] Both targets have three distinct atoms of positive mass.

The coupling that fixes the outer atoms and moves the central atom from \(0\) to \(be_2\) has squared cost \(ab^2\). Conversely, every point in the support of \(\mu_a\) has second coordinate zero. Every coupling \(\pi\) of the two targets therefore satisfies \[\int|y-z|^2\,d\pi(y,z) \ge\int|y_2-z_2|^2\,d\pi(y,z)=ab^2.\] Consequently, \[ W_2(\mu_a,\nu_a)^2=ab^2. \tag{3}\] The same matching coupling costs \(ab\) for \(W_1\), while every coupling has cost at least \(\int|y_2-z_2|\,d\pi=ab\). Thus \[ W_1(\mu_a,\nu_a)=ab=\frac{a^2}{2}. \tag{4}\]

The two central intervals on the slice \(x_2=t\) are nested, and their symmetric difference has conditional probability \(b|t|\). Thus \[\rho(A_a\mathbin{\triangle}B_a) =\frac12\int_{-1}^1b|t|\,dt=\frac b2.\] The first coordinate of the map difference contributes the indicator of that symmetric difference, and its second coordinate contributes \(b^2\) on \(B_a\). Orthogonality gives the pointwise identity \[|T_{\mu_a}-T_{\nu_a}|^2 =1_{A_a\mathbin{\triangle}B_a}+b^2 1_{B_a} \quad\rho\hbox{-almost everywhere}.\] Since \(\rho(B_a)=a\), this proves \[ \|T_{\mu_a}-T_{\nu_a}\|_{L^2(\rho)}^2=\frac b2+ab^2. \tag{5}\] The remaining \(d-2\) coordinates integrate out, so all computations hold in every fixed dimension \(d\ge2\).

Substituting \(b=a/2\) into (3) and (5) gives \[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 every \(\alpha>1/3\), their ratio is \[\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} \longrightarrow+\infty\qquad(a\downarrow0).\] The source, target container, and number of atoms stay fixed, proving the proposition. ◻

Proof. We first record the geometric facts that permit differentiation for an arbitrary convex body. The boundary of a full-dimensional compact convex set has zero volume: contractions toward an interior point lie in the interior and their volumes tend to the full volume. The same argument applies within an affine subspace. Moreover, if a hyperplane \(H\) meets \(\operatorname{int}K\), then \[\operatorname{relint}(K\cap H)=\operatorname{int}K\cap H.\] Indeed, choose \(p\in\operatorname{int}K\cap H\). A relative interior point \(z\ne p\) of the section lies strictly between \(p\) and a point of the section beyond \(z\). A strict convex combination with the ambient interior point \(p\) lies in \(\operatorname{int}K\). The reverse inclusion is immediate. It follows that \(\partial K\cap H\) has zero \((d-1)\)-dimensional measure.

Since \(M_i>0\), cell \(i\) contains an interior point of \(K\) at which its affine function is strictly maximal: remove the boundary of \(K\) and the finite union of tie planes from the positive-volume cell. Strictly active points for \(i\) and \(j\) give opposite signs of their affine difference, so \(H_{ij}\) meets \(\operatorname{int}K\). Also, \(H_{ij}\) and \(H_{ik}\) cannot coincide for three distinct indices. Coincidence would make their affine differences proportional. The three affine functions would then lie on one affine line, with one a strict convex combination of the other two; that function could never be strictly maximal. Hence the intersections \(H_{ij}\cap H_{ik}\) are empty or have dimension at most \(d-2\).

We next prove the differential formulas. For a continuous function \(f:\mathbb R^d\to\mathbb R\), write \[I_i^f=\frac1{|K|}\int_K f(x) \prod_{j\ne i}\mathbf1_{\{\ell_{ij}(x)\ge0\}}\,dx, \qquad \ell_{ij}(x)=x\cdot(y_i-y_j)+c_i-c_j.\] Only \(f=1\) and \(f(x)=x_k\) will be needed. Perturb \((c,y)\) by \(h(r,v)\), with \(h\downarrow0\), and telescope the difference of the products of indicators one factor at a time. On bounded \(K\), the term in which the \(j\)-th factor changes is supported in \(|\ell_{ij}(x)|\le Ch\), with \(C\) fixed for the configuration and direction. Put \(n_{ij}=(y_i-y_j)/|y_i-y_j|\) and use coordinates \[x=z+htn_{ij},\qquad z\in H_{ij}.\] Both \(z\) and \(t\) can be restricted to fixed bounded sets. The factor \(h\) in the volume element cancels the difference quotient. Outside \[(\partial K\cap H_{ij}) \ \cup\!\bigcup_{k\ne i,j}(H_{ik}\cap H_{ij}),\] a surface-null set by the preceding geometric observations, all remaining indicators converge to their values at \(z\). The changed factor converges for almost every \((z,t)\) to \[\mathbf1_{\{|y_i-y_j|t+r_i-r_j+z\cdot(v_i-v_j)\ge0\}} -\mathbf1_{\{t\ge0\}}.\] Its integral in \(t\) is \((r_i-r_j+z\cdot(v_i-v_j))/|y_i-y_j|\). Dominated convergence therefore gives \[ D I_i^f[r,v] =\frac1{|K|}\sum_{j\ne i}\frac1{|y_i-y_j|} \int_{F_{ij}} f(z) \bigl(r_i-r_j+z\cdot(v_i-v_j)\bigr)\,dS(z). \tag{10}\] Applying the same argument in the opposite direction supplies the two-sided derivative. Taking \(f=1,x_1,\ldots,x_d\) gives (7) and (8).

For completeness, these derivatives are continuous under joint changes of sites and intercepts. Identify a nearby equality plane with the fixed \(H_{ij}\) as an affine graph in the normal direction \(n_{ij}\). The graph maps and their surface Jacobians converge uniformly on a common bounded integration region. The face-defining indicators converge outside the same surface-null sets used above. The remaining integrands and the positive denominators \(|y_i-y_j|\) converge as well. Dominated convergence proves continuity of each coefficient of the differential, including those containing \(z_kz_l\). Thus all coordinate partial derivatives exist and are continuous, which yields joint \(C^1\) regularity. The admissible parameter set is open: the masses themselves are continuous by dominated convergence away from the original tie hyperplanes.

To prove connectivity, fix strictly active interior points in any two cells. They may be chosen so that the first point \(p\) avoids all equality planes and the second point \(q\) avoids all equality planes and the affine spans of \(p\) with every triple equality set. These are finitely many proper affine subsets, so they cannot cover a strictly active open region. The segment \([p,q]\) lies in \(\operatorname{int}K\) and meets no triple equality set. At any change of active index along it, precisely two affine functions are maximal. The other inequalities are strict, so a small relatively open disk in their equality plane belongs to the common face. That face has positive area and gives an edge of the graph. The successive active indices therefore give a graph path between the chosen cells.

Finally, symmetry of the face coefficients gives \[q\cdot Gq=\sum_{i<j}g_{ij}(q_i-q_j)^2.\] Connectivity identifies the kernel as the constant vectors. Symmetry then identifies the range as their orthogonal complement, as claimed. ◻

Proof. The case \(N=1\) has \(c_1=0\), so assume \(N\ge2\). On the normalization space \(H_m=\{c\in\mathbb R^N:\sum_i m_i c_i=0\}\), minimize the convex function \[\mathcal F_y(c)=\int_K\max_i(x\cdot y_i+c_i)\,d\rho(x) -\sum_i m_i c_i.\] Choose \(R\) with \(|x|\le R\) on \(K\), and set \(S=\max_i|y_i|\). On \(H_m\), \[\mathcal F_y(c)\ge\max_i c_i-RS.\] Positivity of all the weights implies that \(\max_i c_i\) tends to infinity as \(|c|\to\infty\) within \(H_m\): an upper bound on all coordinates and the weighted mean-zero identity also bound every coordinate below. Hence \(\mathcal F_y\) is coercive on \(H_m\) and attains its minimum.

Distinct sites give an almost-everywhere unique maximizing index for every \(c\). Differentiation of the maximum under the integral, justified by bounded difference quotients, shows \[\nabla_c\mathcal F_y=(M_i-m_i)_{i=1}^N.\] At a minimum on \(H_m\), this gradient lies in \(\operatorname{span}\{m\}\). Its coordinate sum is zero because the masses and weights each sum to one. Since \(\sum_i m_i=1\), the gradient must vanish. Thus a normalized solution of (11) exists. Conversely, every such solution is a minimizer by convexity.

At any solution all cells have positive mass, so 4 applies. The intercept derivative of the mass residual is the restriction \[G:H_m\longrightarrow H_0, \qquad H_0=\{z\in\mathbb R^N:\textstyle\sum_i z_i=0\}.\] It is an isomorphism: both spaces have dimension \(N-1\), and the only constant vector in \(H_m\) is zero. The finite-dimensional implicit function theorem therefore gives a locally unique \(C^1\) solution as the sites vary. In particular, each solution is locally isolated for its fixed sites. If two normalized solutions existed, their segment would consist of minimizers of the convex objective, contradicting this local isolation. Global uniqueness follows, and the local \(C^1\) solution branches agree wherever they overlap. This proves the claimed \(C^1\) dependence on the sites. Only the invertibility of the restricted Laplacian is used here; no quantitative bound on its inverse is needed. ◻

Proof. We parameterize the two uniform laws by successive one-dimensional quantiles on the same cube. Their diagonal derivatives have constant products. Averaging the two maps and substituting coordinates successively will give the density bound.

We need the boundary observation from the proof of 4 for arbitrary positive-dimensional affine slices, not only hyperplanes. The same segment-extension argument shows that a relative interior point of a slice of \(U\) meeting \(\operatorname{int}U\) is an ambient interior point: unless it equals a chosen ambient interior point in the slice, it lies strictly between that point and a farther point of the section. The contraction argument for convex boundaries, applied within the slice, then shows that its intersection with \(\partial U\) has zero slice volume.

For \(0\le k<d\), let \[p_k(\xi)=\int_{\mathbb R^{d-k}} 1_{\operatorname{int}U}(\xi,\eta)\,d\eta, \qquad p_d=1_{\operatorname{int}U}.\] Thus \(p_0=|U|\). Write \(D_k\) for the projection of \(\operatorname{int}U\) onto its first \(k\) coordinates. Each \(p_k\) is positive and continuous on \(D_k\). For \(k<d\), continuity follows by dominated convergence: at a prefix in \(D_k\), the preceding boundary property makes the indicators of the fibers converge almost everywhere as the prefix varies, and boundedness of \(U\) provides a common integrable bound. For \(k=d\) the assertion is immediate. The fiber integrals are Borel functions on their full domains.

Given \(s=(s_1,\ldots,s_d)\in(0,1)^d\), construct \(X(s)\) recursively. For a previously constructed 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\longmapsto \frac{p_k(\xi,t)}{p_{k-1}(\xi)}.\] Fubini’s theorem shows that this density integrates to one. It is positive and continuous on the bounded open interval \(\{t:(\xi,t)\in D_k\}\), and vanishes outside that interval. Its cumulative distribution function is therefore strictly increasing from zero to one there. This defines the quantile uniquely, keeps the new prefix in \(D_k\), and gives \[\partial_{s_k}X_k =\frac{p_{k-1}(X_1,\ldots,X_{k-1})} {p_k(X_1,\ldots,X_k)}>0.\] For each fixed earlier parameter prefix, \(X_k\) is \(C^1\) in \(s_k\). The conditional cumulative distribution functions are Borel in their arguments, so comparison at rational points shows that their quantiles are Borel. Thus the whole construction is measurable. With \(s\) uniform on the cube, successive conditioning gives the uniform law on \(U\). The diagonal derivatives telescope to \[\prod_{k=1}^d\partial_{s_k}X_k=|U|.\] Construct \(Z(s)\) in the same way for \(V\), using the same parameters. For \(H=(X+Z)/2\), the arithmetic–geometric mean inequality yields \[\prod_{k=1}^d\partial_{s_k}H_k \ge \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|}.\]

We justify substitution using only the stated regularity in each coordinate. For a nonnegative Borel function \(g\) on \(\mathbb R^d\), fix \(a=(s_1,\ldots,s_{d-1})\). One-dimensional substitution in the strictly increasing \(C^1\) function \(s_d\mapsto H_d(a,s_d)\), followed by enlargement of its image interval to \(\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 \le \int_{\mathbb R}g\bigl(H_1(a),\ldots,H_{d-1}(a),z\bigr)\,dz.\] Here each \(H_k(a)\) uses only the first \(k\) entries of \(a\). The integral of \(g\) in its last variable is again a nonnegative Borel function. The earlier derivative factors depend only on earlier parameters. Tonelli’s theorem therefore permits multiplication by their product and integration in \(a\). Repeating this argument backwards gives \[\int_{(0,1)^d}g(H(s)) \prod_{k=1}^d\partial_{s_k}H_k(s)\,ds \le \int_{\mathbb R^d}g(x)\,dx.\] All derivative factors are Borel, as also follows from their explicit formulas above. No differentiability in earlier parameters or boundedness of the derivatives near the boundary of the cube is required. Combining the last two inequalities proves \[\mathbb Eg\!\left(\frac{X+Z}{2}\right) \le \frac{1}{\sqrt{|U||V|}}\int_{\mathbb R^d}g(x)\,dx,\] which is the asserted measure domination. ◻

Proof of 6. Fix the intercept direction \(q\) and hold the sites fixed. For sufficiently small \(s>0\), let \(P_i^-\) and \(P_i^+\) be the cells at \(c-sq\) and \(c+sq\), respectively, and write \[\alpha_i=\rho(P_i^-),\qquad \beta_i=\rho(P_i^+),\qquad \gamma_i=\sqrt{\alpha_i\beta_i}.\] These masses are positive by continuity from 4. The perturbed cells are compact convex sets of positive volume. Couple their uniform probability measures by 7, writing \(X_i\in P_i^-\), \(Z_i\in P_i^+\), and \(D_i=\mathbb E|X_i-Z_i|^2\). For each competing index \(j\), the defining inequalities give \[\begin{align*} X_i\cdot(y_i-y_j)+c_i-c_j-s(q_i-q_j)&\ge0,\\ Z_i\cdot(y_i-y_j)+c_i-c_j+s(q_i-q_j)&\ge0. \end{align*}\] Averaging shows that \((X_i+Z_i)/2\) satisfies the \(i\)-cell inequalities at \(c\); convexity also places it in \(K\). Its law is supported on \(P_i\). Since \(|P_i^-|=|K|\alpha_i\) and \(|P_i^+|=|K|\beta_i\), the midpoint bound has exactly the normalization \[\gamma_i\operatorname{Law}\!\left(\frac{X_i+Z_i}{2}\right) \le \frac{1_{P_i}\,dx}{|K|}=\rho|_{P_i}.\]

Use the nonnegative function \(\Phi(x)=R^2-|x|^2\) on \(K\), and set \(A_i^-=\mathbb E\Phi(X_i)\) and \(A_i^+=\mathbb E\Phi(Z_i)\). The 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}\] and the cellwise measure domination imply \[\sum_i\gamma_i\left(\frac{A_i^-+A_i^+}{2}+\frac{D_i}{4}\right) \le\int_K\Phi\,d\rho =\frac12\sum_i(\alpha_iA_i^-+\beta_iA_i^+).\] The last equality uses the partitions at the two perturbed intercepts; all three partitions have null overlaps because the sites are distinct.

Put \(\delta_i=\alpha_i-\beta_i\) and \[\mathcal A=\sum_i\gamma_iD_i, \qquad \mathcal B=\sum_i\frac{\delta_i^2}{\gamma_i}.\] Subtracting the averaged endpoint terms in the preceding inequality gives \[\frac{\mathcal A}{4} \le\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].\] Here \(0\le A_i^\pm\le R^2\), and \[\frac{\alpha_i+\beta_i}{2}-\gamma_i =\frac{\delta_i^2}{2(\sqrt{\alpha_i}+\sqrt{\beta_i})^2} \le\frac{\delta_i^2}{8\gamma_i}.\] Moreover, since \(|X_i|,|Z_i|\le R\), \[|A_i^--A_i^+| \le 2R\mathbb E|X_i-Z_i|\le2R\sqrt{D_i}.\] Cauchy–Schwarz therefore yields \[\frac{\mathcal A}{4} \le\frac{R^2\mathcal B}{8} +\frac R2\sqrt{\mathcal A\mathcal B}.\] Using \(\frac R2\sqrt{\mathcal A\mathcal B} \le\mathcal A/8+R^2\mathcal B/2\) and rearranging, we obtain \[\mathcal A\le5R^2\mathcal B.\] This argument includes the cases where either quantity vanishes.

The cell-moment difference is \[\Delta B_i:=B_i(c-sq)-B_i(c+sq) =\alpha_i\mathbb EX_i-\beta_i\mathbb EZ_i,\] where the notation suppresses the fixed sites. Consequently, \[|\Delta B_i| \le R|\delta_i|+\alpha_i\sqrt{D_i}.\] For this fixed finite configuration and fixed direction \(q\), continuity allows \(s\) to be chosen small enough that every \(\alpha_i,\beta_i,\gamma_i\) lies between \(M_i/2\) and \(2M_i\). Set \(\mathcal E=\sum_i\delta_i^2/M_i\). These comparisons give \(\alpha_i^2/M_i\le8\gamma_i\) and \(\mathcal B\le2\mathcal E\). It follows that \[\begin{align*} \sum_i\frac{|\Delta B_i|^2}{M_i} &\le2R^2\mathcal E+2\sum_i\frac{\alpha_i^2}{M_i}D_i\\ &\le2R^2\mathcal E+16\mathcal A\\ &\le162R^2\mathcal E. \end{align*}\] Divide by \((2s)^2\) and let \(s\) decrease to zero. By 4, \(\Delta B_i/(2s)\to-D B_i[q,0]\) and \(\delta_i/(2s)\to-D M_i[q,0]\), which proves the claim.

The permitted threshold for \(s\) may depend on the direction and the finite configuration, including its smallest mass. Only the uniform constant \(162\) remains after taking the derivative limit. Thus the estimate has no dependence on the number of cells, their masses, or the separation of the sites. ◻

Proof. Choose an optimal finite coupling of the targets. After deleting zero weights and combining identical pairs, write it as \[\pi=\sum_{i=1}^N m_i\delta_{(y_i^0,y_i^1)},\qquad m_i>0,\quad \sum_i m_i=1.\] Finite-dimensional compactness gives existence of such a coupling. Set \[v_i=y_i^1-y_i^0,\qquad y_i(t)=(1-t)y_i^0+ty_i^1,\qquad w=W_2(\mu,\nu).\] Then \[ \sum_i m_i|v_i|^2=w^2. \tag{13}\]

Optimality implies \[ (y_i^0-y_j^0)\cdot(y_i^1-y_j^1)\ge0. \tag{14}\] Indeed, a negative scalar product would permit swapping equal positive portions of the masses at these two coupling pairs, decreasing their quadratic cost. For \(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+2t(1-t)a\cdot b>0 \qquad(0<t<1).\] The last inequality holds because the endpoint pairs are distinct. Thus the sites are distinct at every interior time, even if some endpoint sites coincide.

By 5, for \(0<t<1\) there are unique intercepts \(c_i(t)\), normalized by \[\sum_i m_ic_i(t)=0,\] whose cells have masses \(M_i(t)=m_i\). These intercepts are \(C^1\) in \(t\). Let \(r_i=c_i'(t)\). In particular, \[ \sum_i m_ir_i=0. \tag{15}\] All derivatives in the next calculation are taken at this fixed interior configuration.

Differentiate the fixed-mass constraints using 4. For every intercept direction \(q\), symmetry of the mass and moment face coefficients gives \[ \sum_i r_i\,DM_i[q,0] =-\sum_i v_i\cdot DB_i[q,0]. \tag{16}\] To see the symmetry directly, the mass part of \(\sum_iq_iDM_i[r,v]\) is \(\sum_{i<j}g_{ij}(q_i-q_j)(r_i-r_j)\). Its site part is \(\sum_{i<j}(q_i-q_j)b_{ij}\cdot(v_i-v_j)\), which equals \(\sum_iv_i\cdot DB_i[q,0]\).

The range statement in 4, together with (15), permits a choice of \(q\) satisfying \[DM_i[q,0]=m_ir_i\qquad(1\le i\le N).\] No bound on this \(q\) is needed. Put \(E_r=\sum_i m_ir_i^2\). Using (16), weighted Cauchy–Schwarz, (13), and 6, we obtain \[\begin{split} E_r &\le \left(\sum_i m_i|v_i|^2\right)^{1/2} \left(\sum_i\frac{|DB_i[q,0]|^2}{m_i}\right)^{1/2}\\ &\le \sqrt{C_0}Rw \left(\sum_i\frac{|m_ir_i|^2}{m_i}\right)^{1/2} =\sqrt{C_0}Rw\sqrt{E_r}. \end{split}\] Consequently \[ \left(\sum_i m_i|r_i|^2\right)^{1/2} \le\sqrt{C_0}Rw. \tag{17}\] This includes \(E_r=0\). If \(N=1\), the normalization gives \(r_1=0\), and the same conclusion is immediate.

Define \[u_t(x)=\max_i\{x\cdot y_i(t)+c_i(t)\}.\] On compact time intervals inside \((0,1)\), this function is uniformly Lipschitz in time for \(x\in K\). At a strict maximum its time derivative is \(x\cdot v_i+r_i\). For each fixed time, the exceptional tie set has zero source measure. Fubini therefore identifies this derivative for almost every \((t,x)\). Since its active cells have masses \(m_i\), \[\|\partial_tu_t\|_{L^2(\rho)} \le R\left(\sum_i m_i|v_i|^2\right)^{1/2} +\left(\sum_i m_i|r_i|^2\right)^{1/2} \le(1+\sqrt{C_0})Rw.\] The fundamental theorem of calculus and Minkowski’s inequality imply \[ \|u_t-u_{t'}\|_{L^2(\rho)} \le(1+\sqrt{C_0})Rw\,|t-t'| \qquad(0<t,t'<1). \tag{18}\]

It remains to justify the endpoint maps without imposing separation of their sites. For each fixed \(i\), (17) gives \[|c_i'(t)|\le\frac{\sqrt{C_0}Rw}{\sqrt{m_i}}.\] Hence \(c_i(t)\) has finite limits at \(0\) and \(1\). This individual estimate is used only for endpoint existence; its mass dependence does not enter (18). The finite maxima converge uniformly on \(K\) to their endpoint maxima \(u_0,u_1\), so that inequality extends to the endpoints.

At either endpoint, merge repeated slopes by retaining their largest intercept. Outside a finite union of hyperplanes there is a unique maximizing slope. For any such \(x\), all maximizing slopes at approaching interior times converge to this endpoint slope. Otherwise a subsequence with a fixed maximizing index would contradict the strict gap at the endpoint. The interior slope law is \(\sum_i m_i\delta_{y_i(t)}\), and all slopes are bounded by \(L\). Bounded convergence against continuous tests shows that the endpoint slope laws are \(\mu\) and \(\nu\). By 2, these slopes are the unique optimal maps to their laws. Taking \(t=0,t'=1\) in (18) proves the proposition. ◻

Proof. Fix \(j\) and write \[K_{j,h}=\{x\in K:x+he_j\in K\},\qquad D_{j,h}u(x)=\frac{u(x+he_j)-u(x)}{h} \quad(x\in K_{j,h}).\] Almost every nonempty coordinate slice of \(K\) 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 therefore gives all coordinate derivatives almost everywhere in \(K\). For either function \(u\), convexity gives, at points where both derivatives exist, \[0\le \delta_u(x):=D_{j,h}u(x)-\partial_j u(x) \le \partial_j u(x+he_j)-\partial_j u(x)\le 2L.\] If \(\beta-\alpha\ge h\), cancellation of translated integrals gives \[\begin{align*} &\int_\alpha^{\beta-h} \bigl(\partial_j u(t+h)-\partial_j u(t)\bigr)\,dt\\ &\hspace{15mm}= \int_{\beta-h}^{\beta}\partial_j u(t)\,dt -\int_\alpha^{\alpha+h}\partial_j u(t)\,dt \le 2Lh. \end{align*}\] The coordinates transverse to the slice have been suppressed here. If the slice length is less than \(h\), its allowed portion is empty. Thus, in all cases, integrating \(\delta_u^2\le 2L\delta_u\) over slices and then over the orthogonal projection yields \[ \int_{K_{j,h}}\delta_u^2\,d\rho \le \frac{4L^2hA_j}{|K|}. \tag{20}\] The excluded portion of any slice has length at most \(h\), so \[\rho(K\setminus K_{j,h})\le \frac{hA_j}{|K|}.\] There we use \(|\partial_j f|\le 2L\).

On the allowed set, both \(x\) and \(x+he_j\) belong to \(K\). Translation invariance of Lebesgue measure and \((a-b)^2\le 2a^2+2b^2\) give \[\int_{K_{j,h}}|D_{j,h}f|^2\,d\rho \le \frac{4}{h^2}\|f\|_{L^2(\rho)}^2.\] Since \(\partial_jf=D_{j,h}f-\delta_{u_0}+\delta_{u_1}\), the inequality \((a+b+c)^2\le 3(a^2+b^2+c^2)\) and (20) imply \[\int_K|\partial_jf|^2\,d\rho \le \frac{12}{h^2}\|f\|_{L^2(\rho)}^2 +\frac{(24+4)L^2hA_j}{|K|}.\] Summing over \(j\) proves (19). All projection measures are finite because \(K\) is compact. The slice argument applies to every \(h>0\), including values larger than some or all slice lengths.

The coordinate derivatives form the full gradient almost everywhere. Indeed, at an interior point \(x\) where all partials of a convex function \(u\) exist, set \(p_j=\partial_j u(x)\) and \(q=\|v\|_1\) for a small nonzero vector \(v\). The point \(x+v\) is the convex combination of \(x+\operatorname{sign}(v_j)q e_j\) with weights \(|v_j|/q\), omitting zero weights. The finitely many coordinate derivative limits and convexity give \(u(x+v)-u(x)\le p\cdot v+o(q)\). Apply this bound to \(-v\) and use midpoint convexity to obtain the reverse bound. Since \(q\le\sqrt d\,|v|\), the function is differentiable at \(x\) with gradient \(p\). The boundary of \(K\) has zero volume, as proved in 4, so this applies almost everywhere to both \(u_i\). ◻

Proof. Write \(w=W_2(\mu,\nu)\). By 8, the maps have convex, \(L\)-Lipschitz maximum potentials \(u_0,u_1\) satisfying \[\|u_0-u_1\|_{L^2(\rho)} \le (1+\sqrt{C_0})Rw.\] If \(w>0\), apply 10 with \(h=w^{2/3}\). The resulting bound on the squared map norm is \(C_*^2w^{2/3}\). Taking its square root proves the claim. If \(w=0\), the potential difference has zero \(L^2(\rho)\) norm, and sending \(h\) to zero in (19) gives the same conclusion. ◻

Proof. Compactness of \(Y\) gives the finite nets. Enumerating each net and choosing the first closest point gives a Borel selection with \(|Q_n(y)-y|\le\epsilon_n\). Put \(\lambda_n=(Q_n)_\#\lambda\), discard any zero-mass sites, and let \(T_n=T_{\lambda_n}\). These finite maps exist by [cells:weights,cells:optimality]. The coupling induced by \(y\mapsto(Q_n(y),Q_m(y))\) gives \[W_2(\lambda_n,\lambda_m)\le\epsilon_n+\epsilon_m.\] By 11, \((T_n)\) is Cauchy in \(L^2(\rho)\). Let \(S\) be its limit and choose a subsequence converging almost everywhere. Its limit lies in \(Y\) almost everywhere because \(Y\) is closed. Modifying \(S\) on a null set gives a Borel map with values in \(Y\). For any continuous \(\phi:Y\to\mathbb R\), bounded convergence and uniform continuity give \[\int_K\phi(S(x))\,d\rho(x) =\lim_n\int_Y\phi(Q_n(y))\,d\lambda(y) =\int_Y\phi(y)\,d\lambda(y).\] Here and below a limit along the chosen subsequence is understood where needed. It follows that \(S_\#\rho=\lambda\).

To prove optimality directly, choose finite maximum potentials \(u_n\) for \(T_n\) and subtract a constant so that \(u_n(x_*)=0\) at a fixed \(x_*\in K\). Their slopes lie in \(Y\), so they are \(L\)-Lipschitz and \(|u_n|\le L\mathop{\mathrm{diam}}(K)\) on \(K\). Arzelà–Ascoli gives a further subsequence converging uniformly on \(K\) to a convex, \(L\)-Lipschitz function \(u\). Intersect the full-measure sets on which each of the countably many finite-map supporting inequalities holds with the set on which the selected subsequence of \(T_n\) converges. On this common set, \[u_n(z)\ge 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)\ge u(x)+S(x)\cdot(z-x) \quad\hbox{for $\rho$-almost every }x. \tag{22}\] Define the restricted convex conjugate \[u^*(y)=\max_{z\in K}\bigl(z\cdot y-u(z)\bigr).\] It is finite and continuous, and hence bounded on \(Y\). The defining inequality and (22) give \[ x\cdot y\le u(x)+u^*(y)\quad(x\in K,\ y\in Y), \qquad x\cdot S(x)=u(x)+u^*(S(x))\quad\rho\hbox{-almost everywhere}. \tag{23}\] For any competing map \(V\) with \(V_\#\rho=\lambda\), the integral of \(u(x)+u^*(V(x))\) is fixed by the marginals. Thus (23) shows that \(S\) maximizes \(\int_Kx\cdot V(x)\,d\rho(x)\). The 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_Kx\cdot V(x)\,d\rho(x)\] proves quadratic optimality.

If \(V\) has the same optimal cost, its nonnegative Fenchel gap \(u(x)+u^*(V(x))-x\cdot V(x)\) has integral zero. Equality and the definition of \(u^*\) imply \(u(z)\ge u(x)+V(x)\cdot(z-x)\) for every \(z\in K\) at almost every \(x\). On almost every interior point of \(K\), all coordinate partials of \(u\) exist: on each coordinate interval the one-sided derivatives of a convex function are monotone, and Fubini applies. At such a point \(x\), a supporting slope \(p\) must satisfy \[p_j\le\lim_{t\downarrow0}\frac{u(x+te_j)-u(x)}{t}, \qquad p_j\ge\lim_{t\uparrow0}\frac{u(x+te_j)-u(x)}{t}.\] Both limits equal \(\partial_j u(x)\), so the supporting slope is unique. The differentiability argument in the proof of 10 identifies it with \(\nabla u(x)\). The boundary of \(K\) has zero volume by 4; hence \(V=S\) almost everywhere. The same integrated inequality also proves optimality among arbitrary transport plans. Equality forces their target coordinate to be this unique supporting slope over almost every source point, so the optimal plan is the one induced by \(S\).

We have identified the \(L^2\) limit as the unique map \(T_\lambda\). The original sequence was Cauchy, so the asserted convergence holds for the full sequence. ◻

Proof of 1. Use the same finite nets and selections for \(\mu\) and \(\nu\), and set \(\mu_n=(Q_n)_\#\mu\), \(\nu_n=(Q_n)_\#\nu\). For any coupling \(\pi\) of \(\mu,\nu\), its image under \((Q_n,Q_n)\) is a coupling of \(\mu_n,\nu_n\). The triangle inequality in \(L^2(\pi)\) gives \[W_2(\mu_n,\nu_n) \le\left(\int_{Y\times Y}|y-z|^2\,d\pi(y,z)\right)^{1/2} +2\epsilon_n.\] Taking the infimum over \(\pi\) yields \(W_2(\mu_n,\nu_n)\le W_2(\mu,\nu)+2\epsilon_n\). Apply 11 and then use the strong \(L^2(\rho)\) convergence from 12 to conclude \[\|T_\mu-T_\nu\|_{L^2(\rho)} \le C_*W_2(\mu,\nu)^{1/3}.\] The constant in (21) depends only on \(K\) and an enclosing radius for \(Y\). In particular it is independent of atom counts, masses, site separations, or target densities. No convexity of \(Y\) is required: the site interpolation used earlier stays in the ball of radius \(L\), and the discretizations take their values in \(Y\). The argument includes a singleton \(Y\), for which all target maps coincide. The sharpness assertion was proved in 3. ◻