Proof. For each \(j\), define the Borel set \[S_j=\{x\in B:T(x)=H_j(x)\}.\] Its image \(H_j(S_j)\) is Borel. Pairwise disjointness of the component images and injectivity of \(H_j\) give \[\mu(S_j)=\frac13\nu(S_j),\qquad \mu(H_j(S_j))=\frac16\nu(S_j).\] Since \(S_j\subset T^{-1}(H_j(S_j))\), preservation of \(\mu\) would imply \[\frac13\nu(S_j) \le \mu(T^{-1}(H_j(S_j))) =\mu(H_j(S_j)) =\frac16\nu(S_j).\] Thus \(\nu(S_j)=0\) for all four indices, contradicting (3). ◻

Proof. For a parameter \(K>0\), use the potentials in the statement and set \[F_k(x,y,z)=c(x,y,z)-u_0(x)-v_k(z).\] By (4), its \((x,z)\)-gradient vanishes at the \(k\)-th center triple. We choose \(K\) to obtain a unique minimizing pair on a small product neighborhood for each fixed \(y\), and to make each coordinate of the resulting minimizing pair an invertible function of \(y\).

The Coulomb cost is smooth near both center triples. At the \(k\)-th center, write \[M_k=D^2_{(x,z)(x,z)}c,\qquad N_k=D_y\nabla_{(x,z)}c.\] The \((x,z)\)-Hessian of \(F_k\) there is \(M_k+KI_6\), while its mixed derivative in \(y\) is \(N_k\). For \(h(b)=|b|^{-1}\), \(b\ne0\), differentiation gives \[ D^2h(b)=|b|^{-3} \left(3\frac{bb^{\mathsf T}}{|b|^2}-I_3\right). \tag{6}\] This matrix has eigenvalues \(2|b|^{-3},-|b|^{-3},-|b|^{-3}\) and is therefore invertible. The two \(3\times3\) blocks of \(N_k\) are \[-D^2h(-a_k),\qquad -D^2h(-2a_k),\] and hence are both invertible. The matrices \(M_k\) are fixed symmetric \(6\times6\) matrices. For all sufficiently large \(K\), both \(M_k+KI_6\) are positive definite. Moreover, \[K(M_k+KI_6)^{-1}N_k =(I_6+K^{-1}M_k)^{-1}N_k\longrightarrow N_k.\] Invertibility is an open condition on \(3\times3\) matrices, so we may choose a single \(K\) for which both blocks of \((M_k+KI_6)^{-1}N_k\) are invertible for each \(k\).

Fix such a \(K\). The implicit function theorem applied to \(\nabla_{(x,z)}F_k=0\) gives a smooth stationary branch \[\nabla_{(x,z)}F_k(X_k(y),y,Z_k(y))=0,\] with the prescribed center values. Differentiating this equation at \(y=a_k\) gives \[D(X_k,Z_k)(a_k)=-(M_k+KI_6)^{-1}N_k.\] The two block invertibilities allow the inverse function theorem to be applied separately to \(X_k\) and \(Z_k\). After restricting to one sufficiently small neighborhood of \(a_k\), both maps are diffeomorphisms onto their open images.

By continuity, the \((x,z)\)-Hessian of \(F_k\) remains positive definite on a product of sufficiently small open balls around the three centers. Shrink the \(y\)-ball so that the stationary pair \((X_k(y),Z_k(y))\) belongs to the product of the other two balls. For fixed \(y\), this latter product is convex. The positive definite Hessian makes \(F_k(\cdot,y,\cdot)\) strictly convex on it, so its stationary pair is its unique minimum there.

Define its minimum value by \[w_k(y)=F_k(X_k(y),y,Z_k(y)).\] The preceding minimum property gives (5) with exactly the stated equality set. The center values in (4) give \(w_k(a_k)=0\). All neighborhoods can be restricted further as needed to retain the separate diffeomorphism properties. ◻

Proof. Use the functions and neighborhoods of Lemma 3. Choose small balls contained in these neighborhoods, taking \(U_0\) inside both central neighborhoods. On the respective components define \[u=u_0\ \text{on }U_0,\qquad u=w_k\ \text{on }U_k^+,\qquad u=v_k\ \text{on }U_k^-.\] The balls will be shrunk below. Since the functions are smooth near their centers, we can first restrict to fixed preliminary balls on which their absolute values are bounded by one constant \(M>0\). Take the final balls sufficiently small that their closures are pairwise disjoint.

At the five centers the values of \(u\) are \(5/2,0,0,0,0\). For a triple of distinct centers other than a permutation of \((0,a_k,-a_k)\), there are two possibilities. If zero is absent, its potential sum is zero and its cost is positive. If zero is present, the other two centers come from different opposite pairs. Their cost is \[2+\frac1{\sqrt2}>\frac52,\] so again there is a strict gap. There are finitely many such triples, and the cost and local potentials are continuous near each one. All these strict inequalities therefore persist on the corresponding products after sufficiently small restrictions.

For a repeated-center type, two coordinates belong to the same ball. If all radii are at most \(r\), the distance between those two coordinates is at most \(2r\), and their pair cost is at least \(1/(2r)\), with the usual interpretation at coincidence. Taking \(r\) small enough that \(1/(2r)>3M\) makes the cost strictly larger than the potential sum on every such product.

The inequalities of Lemma 3 remain valid when their domains are restricted, with the same equality characterization. Permuting a triple leaves both \(c\) and the sum of the three values of \(u\) unchanged. Thus the two prescribed component types and their permutations satisfy (7), with equality exactly as in (8); every other type has strict inequality. Finally, keep \(U_0,U_k^-\) fixed and shrink each \(U_k^+\) once more. Since \(X_k(a_k)=0\) and \(Z_k(a_k)=-a_k\), continuity ensures \[X_k(U_k^+)\subset U_0,\qquad Z_k(U_k^+)\subset U_k^-.\] All inequalities and diffeomorphism properties survive this restriction. The function \(u\) is bounded and continuous on the final disjoint union \(U\). ◻

Proof. Both \(X_k(U_k^+)\) are open neighborhoods of zero. Choose a positive-radius ball \(B\) about zero inside their intersection and \(U_0\). The inverse maps then take \(B\) into \(U_k^+\), and Proposition 4 gives \(Z_k(Y_k(B))\subset U_k^-\). The inverse and composition of the relevant diffeomorphisms are diffeomorphisms onto open images. The identities \(X_k(Y_k(x))=x\) and \(W_k=Z_k\circ Y_k\) give equality by Proposition 4. Disjointness follows from the five component balls. ◻

Proof of Proposition 6. Write \(K_0=\operatorname{supp}g\), and let \(K_1,\ldots,K_4\) be its images under the four maps in \(\mathcal H\). These are nonempty compact sets in the five disjoint components, and \(\mu(K)=1\) for \(K=\bigcup_{a=0}^4K_a\). In particular, \[\eta=\min_{a\ne b}\operatorname{dist}(K_a,K_b)>0.\] The optimal coupling \(\pi_0\) constructed in Section 5 is supported on products of three distinct component sets.

Fix \(\delta>0\). Partition each \(K_a\) into finitely many Borel cells of diameter at most \(\delta\), for example by a finite grid of half-open cubes. Discard the cells of zero \(\mu\)-mass and denote the remaining cells by \(A_1,\ldots,A_N\). They partition a set of full \(\mu\)-measure and each belongs to a single \(K_a\). Put \[q_{ijk}=\pi_0(A_i\times A_j\times A_k).\] The marginal constraints on \(\pi_0\) read \[\sum_{j,k}q_{ijk}=\mu(A_i),\qquad \sum_{i,k}q_{ijk}=\mu(A_j),\qquad \sum_{i,j}q_{ijk}=\mu(A_k).\] For each fixed \(i\), split \(A_i\) into Borel sets \(B_{ijk}\) of masses \(q_{ijk}\), ignoring zero-mass pieces. Such a split can be made by successive first-coordinate slabs: the first-coordinate marginal of \(\mu|_{A_i}\) is absolutely continuous, so its cumulative distribution is continuous and takes every value between zero and \(\mu(A_i)\). Endpoints of the slabs have zero mass; any unused null set is included in a fixed piece.

For each \(q_{ijk}>0\), Lemma 7 supplies maps on \(B_{ijk}\) sending the probability \(q_{ijk}^{-1}\mu|_{B_{ijk}}\), respectively, to \[\frac{\mu|_{A_j}}{\mu(A_j)} \quad\hbox{and}\quad \frac{\mu|_{A_k}}{\mu(A_k)}.\] Reset any exceptional values outside the indicated target cell to a fixed point in that cell. These changes are Borel and occur on null sets. Glue the maps over the finitely many pieces to obtain \(T_2,T_3\), and give them a fixed value on the remaining Borel null set in \(\mathbb R^3\). Different source pieces may use the same normalized target law. Their weights \(q_{ijk}\) add to the full mass of the target cell, so no target partition is needed. Indeed, for any Borel set \(D\), \[\mu(T_2^{-1}(D)) =\sum_{i,j,k:q_{ijk}>0}q_{ijk} \frac{\mu(D\cap A_j)}{\mu(A_j)} =\sum_j\mu(D\cap A_j)=\mu(D),\] and the same calculation proves \((T_3)_\#\mu=\mu\). The graph coupling \(\pi_\delta=(\operatorname{Id},T_2,T_3)_\#\mu\) gives each product \(A_i\times A_j\times A_k\) exactly the mass \(q_{ijk}\).

Whenever \(q_{ijk}>0\), the three cells lie in three distinct component sets. Thus all pairwise distances in their product are at least \(\eta\), including after the redistribution within the cells. If \(x\) and \(y\) are two triples in the same such product, then \(|x_\ell-y_\ell|\le\delta\) for each coordinate. Consequently \[\left|\frac1{|x_\ell-x_r|}-\frac1{|y_\ell-y_r|}\right| \le\frac{2\delta}{\eta^2}, \qquad |c(x)-c(y)|\le\frac{6\delta}{\eta^2}.\] Both costs are finite, bounded by \(3/\eta\) on all active products. Since \(\pi_\delta\) and \(\pi_0\) have identical masses on these products, summing their oscillation bounds gives \[\left|\int c\,d\pi_\delta-\int c\,d\pi_0\right| \le\frac{6\delta}{\eta^2}.\] Letting \(\delta\downarrow0\) proves the required upper bound on the Monge infimum. The reverse bound follows because every admissible graph plan is a coupling. Strict nonattainment is Theorem 1. ◻

Proof of Lemma 7. We construct a Borel map from each such probability to uniform measure on \((0,1)^d\), and a Borel map in the reverse direction. Let \(f\) be a nonnegative Borel probability density. The forward map is Rosenblatt’s successive conditional-distribution transform (Rosenblatt 1952); we also give the inverse quantile construction and make the Borel choices explicit in both directions.

The forward map. We successively transform each coordinate by its conditional distribution function given the preceding coordinates. Write \(x_{1:j}=(x_1,\ldots,x_j)\), with \(x_{1:0}\) the empty prefix, and set \[f_{1:j}(x_{1:j})=\int_{\mathbb R^{d-j}}f(x_{1:j},z)\,dz \quad(1\le j<d),\qquad f_{1:d}=f,\qquad f_{1:0}=1.\] The prefix densities are Borel, have integral one, and are finite almost everywhere. Tonelli’s theorem gives \(\int_\mathbb Rf_{1:j}(a,r)\,dr=f_{1:j-1}(a)\) for every prefix \(a\), including the empty prefix when \(j=1\). For \(j=1,\ldots,d\), define \[F_j(a,t)=\frac{\int_{-\infty}^t f_{1:j}(a,r)\,dr} {f_{1:j-1}(a)},\qquad a\in\mathbb R^{j-1}.\] This formula is used where its denominator belongs to \((0,\infty)\); on the remaining fibers use instead the fixed distribution function \(t\mapsto(1+e^{-t})^{-1}\). Thus every fiber is a continuous distribution function with limits zero and one, and all \(F_j\) are jointly Borel. For \(j\ge2\), the exceptional prefix fibers have zero probability under \(f_{1:j-1}(a)\,da\); there are no exceptional fibers when \(j=1\).

For a random variable with a continuous distribution function \(F\), the variable \(F(X)\) is uniform on \((0,1)\). Indeed, continuity and the limits of \(F\) allow a cut at each probability level \(u\in(0,1)\); any interval on which \(F=u\) has probability zero, and hence \(\mathbb P(F(X)\le u)=u\). At stage \(j\), this identity and Fubini show that, conditionally on the preceding original prefix, the new transformed coordinate is uniform. It is therefore independent of that prefix, and hence of all preceding transformed coordinates. A finite induction over \(j=1,\ldots,d\) shows that the Borel map \[R_f(x)=\bigl(F_j(x_{1:j-1},x_j)\bigr)_{j=1}^d\] sends \(f\,dx\) to product uniform measure on \((0,1)^d\). Its values lie in \([0,1]^d\), including on exceptional fibers. This argument does not require a density to be positive everywhere or a distribution function to be strictly increasing.

The reverse map. Given independent uniform coordinates, we successively recover the original conditional laws by quantiles. For each \(j\) and \(0<u<1\), define \[q_j(a,u)=\inf\{r\in\mathbb Q:F_j(a,r)\ge u\}.\] This value is finite, and \[\{(a,u):q_j(a,u)<t\} =\bigcup_{\substack{r\in\mathbb Q\\r<t}} \{(a,u):F_j(a,r)\ge u\}\] proves joint Borel measurability on \(0<u<1\). At \(u=0,1\), define \(q_j(a,u)=0\). For \(0<u<1\), continuity of \(F_j(a,\cdot)\) gives \[q_j(a,u)\le t\quad\Longleftrightarrow\quad u\le F_j(a,t).\] Consequently \(q_j(a,U)\) has distribution function \(F_j(a,\cdot)\) for uniform \(U\). For \(u\in[0,1]^d\), recursively set \[z_j=q_j(z_{1:j-1},u_j)\quad(j=1,\ldots,d), \qquad S_f(u)=(z_1,\ldots,z_d),\] starting with the empty prefix. This finite recursion defines a Borel map. For independent uniform inputs, induction shows that after step \(j\) the output prefix has density \(f_{1:j}\): the preceding prefix has density \(f_{1:j-1}\), its exceptional fibers have probability zero, and on every other fiber the quantile has the conditional distribution just defined. Thus \(S_f\) has output law \(f\,dx\). If \(f\) and \(g\) are densities of \(\alpha\) and \(\beta\), respectively, then \(S_g\circ R_f\) has the required property. ◻

Proof. Fix \(d\ge2\) and \(s>0\), and put \(a_1=e_1\) and \(a_2=e_2\) in \(\mathbb R^d\). The measure argument already applies in this dimension. We verify the three properties of the cost that it needs: local equality families given by diffeomorphisms, strict supporting inequalities on all other component products, and a uniform cost-oscillation bound on separated products.