[proofplan] The quadratic identity follows by expanding $|x-y|^2$ and observing that the marginal constraints make the two second-moment terms independent of the transport plan. For the structural statement, we use [Kantorovich duality](/theorems/6799) for the quadratic cost to obtain optimal dual potentials and then rewrite those potentials as a dual pair for maximizing the correlation $x \cdot y$. Convex conjugation replaces this correlation-dual pair by a Fenchel pair $(\varphi,\varphi^*)$ without increasing its integral value. [Complementary slackness](/theorems/2559) then forces equality in Fenchel's inequality on the optimal plan, and Fenchel-Young equality identifies the equality set with the graph of the convex subdifferential. [/proofplan]

[step:Expand the quadratic cost and separate the fixed marginal terms] For every $\pi \in \Pi(\mu,\nu)$, the function \begin{align*} (x,y) \mapsto x \cdot y \end{align*} is $\pi$-integrable because [Young's inequality](/theorems/244) gives \begin{align*} |x \cdot y| \leq \frac{1}{2}|x|^2 + \frac{1}{2}|y|^2 \end{align*} and $\mu,\nu \in \mathcal P_2(\mathbb R^n)$. Since $\pi$ has marginals $\mu$ and $\nu$, we have \begin{align*} \int_{\mathbb R^n \times \mathbb R^n} |x|^2 \, d\pi(x,y)=\int_{\mathbb R^n} |x|^2 \, d\mu(x) \end{align*} and \begin{align*} \int_{\mathbb R^n \times \mathbb R^n} |y|^2 \, d\pi(x,y)=\int_{\mathbb R^n} |y|^2 \, d\nu(y). \end{align*} Using the pointwise identity \begin{align*} |x-y|^2=|x|^2+|y|^2-2x \cdot y, \end{align*} we obtain, for every $\pi \in \Pi(\mu,\nu)$, \begin{align*} \int_{\mathbb R^n \times \mathbb R^n} |x-y|^2 \, d\pi(x,y)=A-2\int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi(x,y), \end{align*} where the finite constant $A$ is defined by \begin{align*} A:=\int_{\mathbb R^n} |x|^2 \, d\mu(x)+\int_{\mathbb R^n} |y|^2 \, d\nu(y). \end{align*} Taking the infimum over $\pi \in \Pi(\mu,\nu)$ on the left side is therefore the same as taking the supremum of the correlation term on the right side: \begin{align*} \inf_{\pi \in \Pi(\mu,\nu)} \int_{\mathbb R^n \times \mathbb R^n} |x-y|^2 \, d\pi(x,y)=A-2\sup_{\pi \in \Pi(\mu,\nu)}\int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi(x,y). \end{align*} [/step]

[step:Obtain quadratic dual potentials and complementary slackness]By the existence theorem for optimal transport plans with lower semicontinuous costs bounded from below by an integrable function, applied to the continuous nonnegative cost \begin{align*} c: \mathbb R^n \times \mathbb R^n \to [0,\infty), \qquad c(x,y):=|x-y|^2, \end{align*} there exists a plan $\pi_0 \in \Pi(\mu,\nu)$ attaining the quadratic-cost infimum. We now fix any such optimal plan $\pi_0$. We now use the Kantorovich duality hypothesis stated in the theorem statement. Applied to $\pi_0$, it gives Borel functions \begin{align*} f: \mathbb R^n \to \mathbb R, \qquad g: \mathbb R^n \to \mathbb R \end{align*} such that $\int_{\mathbb R^n} |f| \, d\mu < \infty$ and $\int_{\mathbb R^n} |g| \, d\nu < \infty$, \begin{align*} f(x)+g(y)\leq |x-y|^2 \quad \text{for all } x,y \in \mathbb R^n, \end{align*} and \begin{align*} \int_{\mathbb R^n} f \, d\mu+\int_{\mathbb R^n} g \, d\nu=\int_{\mathbb R^n \times \mathbb R^n} |x-y|^2 \, d\pi_0(x,y). \end{align*} Since the nonnegative function $(x,y)\mapsto |x-y|^2-f(x)-g(y)$ has integral zero with respect to $\pi_0$, it vanishes $\pi_0$-almost everywhere: \begin{align*} f(x)+g(y)=|x-y|^2 \quad \text{for } \pi_0\text{-almost every }(x,y). \end{align*}[/step]

[guided]The role of Kantorovich duality is to replace the global minimization over all couplings by pointwise inequalities involving two one-variable functions. The cost function is \begin{align*} c: \mathbb R^n \times \mathbb R^n \to [0,\infty), \qquad c(x,y):=|x-y|^2. \end{align*} It is continuous, hence lower semicontinuous, and nonnegative. Since $\mu,\nu \in \mathcal P_2(\mathbb R^n)$, the coupling set is nonempty: define the product coupling $\pi_{\mathrm{prod}} \in \Pi(\mu,\nu)$ by \begin{align*} \pi_{\mathrm{prod}}(A \times B):=\mu(A)\nu(B) \end{align*} for all Borel sets $A,B \subset \mathbb R^n$, and extend this prescription uniquely to the Borel product $\sigma$-algebra on $\mathbb R^n \times \mathbb R^n$. This coupling has finite quadratic cost because \begin{align*} |x-y|^2 \leq 2|x|^2+2|y|^2. \end{align*} Thus the usual existence theorem for lower semicontinuous transport costs applies and gives an optimal plan $\pi_0 \in \Pi(\mu,\nu)$. Next we use the Kantorovich duality hypothesis stated in the theorem statement. This is the strong finite-second-moment quadratic-cost version needed here: it gives real-valued Borel representatives of the dual potentials, integrability with respect to the two marginals, the pointwise dual inequality on all of $\mathbb R^n \times \mathbb R^n$, dual attainment, and complementary slackness for the fixed optimal plan $\pi_0$. Thus we obtain Borel functions \begin{align*} f: \mathbb R^n \to \mathbb R, \qquad g: \mathbb R^n \to \mathbb R \end{align*} such that $\int_{\mathbb R^n} |f| \, d\mu < \infty$ and $\int_{\mathbb R^n} |g| \, d\nu < \infty$, satisfying the dual feasibility inequality \begin{align*} f(x)+g(y)\leq |x-y|^2 \quad \text{for all } x,y \in \mathbb R^n. \end{align*} Dual optimality says that the dual value equals the primal value: \begin{align*} \int_{\mathbb R^n} f \, d\mu+\int_{\mathbb R^n} g \, d\nu=\int_{\mathbb R^n \times \mathbb R^n} |x-y|^2 \, d\pi_0(x,y). \end{align*} Now define the nonnegative slack function \begin{align*} s: \mathbb R^n \times \mathbb R^n \to [0,\infty], \qquad s(x,y):=|x-y|^2-f(x)-g(y). \end{align*} Dual feasibility gives $s \geq 0$. Integrating $s$ with respect to $\pi_0$ and using the marginal identities for $\pi_0$ gives \begin{align*} \int_{\mathbb R^n \times \mathbb R^n} s \, d\pi_0=\int_{\mathbb R^n \times \mathbb R^n} |x-y|^2 \, d\pi_0(x,y)-\int_{\mathbb R^n} f \, d\mu-\int_{\mathbb R^n} g \, d\nu=0. \end{align*} A nonnegative measurable function with integral zero is zero almost everywhere. Therefore \begin{align*} f(x)+g(y)=|x-y|^2 \quad \text{for } \pi_0\text{-almost every }(x,y). \end{align*} This is the complementary slackness relation that will later become equality in Fenchel's inequality.[/guided]

[step:Transform the quadratic dual pair into a correlation dual pair] Define the Borel functions \begin{align*} \alpha: \mathbb R^n \to \mathbb R \cup \{\infty\}, \qquad \alpha(x):=\frac{|x|^2-f(x)}{2} \end{align*} and \begin{align*} \beta: \mathbb R^n \to \mathbb R \cup \{\infty\}, \qquad \beta(y):=\frac{|y|^2-g(y)}{2}. \end{align*} The quadratic dual inequality is equivalent to \begin{align*} x \cdot y \leq \alpha(x)+\beta(y) \quad \text{for all } x,y \in \mathbb R^n. \end{align*} Also, \begin{align*} \int_{\mathbb R^n} \alpha \, d\mu+\int_{\mathbb R^n} \beta \, d\nu=\frac{1}{2}\left(A-\int_{\mathbb R^n} f \, d\mu-\int_{\mathbb R^n} g \, d\nu\right). \end{align*} Using the identity proved in the first step and the optimality of $\pi_0$, this becomes \begin{align*} \int_{\mathbb R^n} \alpha \, d\mu+\int_{\mathbb R^n} \beta \, d\nu=\sup_{\pi \in \Pi(\mu,\nu)}\int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi(x,y). \end{align*} Moreover, complementary slackness for $f,g$ is equivalent to \begin{align*} x \cdot y=\alpha(x)+\beta(y) \quad \text{for } \pi_0\text{-almost every }(x,y). \end{align*} [/step]

[step:Replace the correlation dual pair by a Fenchel dual pair]Since $f$ is finite $\mu$-almost everywhere with $\int_{\mathbb R^n} |f| \, d\mu < \infty$ and $g$ is finite $\nu$-almost everywhere with $\int_{\mathbb R^n} |g| \, d\nu < \infty$, the functions $\alpha$ and $\beta$ are finite $\mu$-almost everywhere and $\nu$-almost everywhere, respectively. Choose $x_0 \in \mathbb R^n$ with $\alpha(x_0)<\infty$ and choose $y_0 \in \mathbb R^n$ with $\beta(y_0)<\infty$. The correlation dual inequality gives the affine lower bound \begin{align*} \beta(y) \geq x_0 \cdot y-\alpha(x_0) \quad \text{for all } y \in \mathbb R^n. \end{align*} Define the convex conjugate of $\beta$ by the function \begin{align*} \varphi: \mathbb R^n \to (-\infty,\infty], \qquad \varphi(x):=\sup_{y \in \mathbb R^n}\{x \cdot y-\beta(y)\}. \end{align*} Because $x \cdot y \leq \alpha(x)+\beta(y)$ for all $x,y \in \mathbb R^n$, we have \begin{align*} \varphi(x)\leq \alpha(x) \quad \text{for all } x \in \mathbb R^n. \end{align*} The function $\varphi$ is convex and lower semicontinuous as a supremum of affine functions. The function $\beta$ is proper because it is finite at $y_0$, never takes the value $-\infty$, and has the affine minorant $y \mapsto x_0 \cdot y-\alpha(x_0)$. Hence its convex conjugate is a well-defined extended-real convex function. The function $\varphi$ is proper: the inequality above gives $\varphi(x_0)\leq \alpha(x_0)<\infty$, while $\beta(y_0)<\infty$ gives \begin{align*} \varphi(x)\geq x \cdot y_0-\beta(y_0)>-\infty \quad \text{for all } x \in \mathbb R^n. \end{align*} Moreover, $\varphi$ is bounded below by the affine function $x \mapsto x \cdot y_0-\beta(y_0)$, and $\varphi^*$ will be bounded below by the affine function $y \mapsto x_0 \cdot y-\varphi(x_0)$. Since $\mu,\nu \in \mathcal P_2(\mathbb R^n)$, both affine lower bounds are integrable with respect to $\mu$ and $\nu$. Thus the extended integrals of $\varphi$ and $\varphi^*$ are well-defined from below and cannot contain an ambiguous $\infty-\infty$ cancellation. Let $\varphi^*: \mathbb R^n \to (-\infty,\infty]$ denote the convex conjugate of $\varphi$: \begin{align*} \varphi^*(y):=\sup_{x \in \mathbb R^n}\{x \cdot y-\varphi(x)\}. \end{align*} Since $\beta$ is proper and bounded below by an affine function, the biconjugate $\beta^{**}$ is well-defined and the Fenchel-Moreau biconjugation inequality gives $\beta^{**}\leq \beta$ pointwise. Since $\varphi=\beta^*$, we have $\varphi^*=\beta^{**}$, and therefore $\varphi^*\leq \beta$ pointwise. Hence \begin{align*} \int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu\leq \int_{\mathbb R^n} \alpha \, d\mu+\int_{\mathbb R^n} \beta \, d\nu. \end{align*} On the other hand, Fenchel's inequality gives \begin{align*} x \cdot y \leq \varphi(x)+\varphi^*(y) \quad \text{for all } x,y \in \mathbb R^n. \end{align*} Integrating this inequality against any $\pi \in \Pi(\mu,\nu)$ gives \begin{align*} \int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi(x,y)\leq \int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu. \end{align*} Taking the supremum over $\pi \in \Pi(\mu,\nu)$ and using the optimality of $(\alpha,\beta)$ for the correlation dual problem yields \begin{align*} \int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu=\sup_{\pi \in \Pi(\mu,\nu)}\int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi(x,y). \end{align*}[/step]

[guided]The delicate point is that $\alpha$, $\beta$, $\varphi$, and $\varphi^*$ may be extended-valued, so we first put affine lower bounds in place before comparing integrals. Since $f$ is finite $\mu$-almost everywhere with $\int_{\mathbb R^n} |f| \, d\mu < \infty$ and $g$ is finite $\nu$-almost everywhere with $\int_{\mathbb R^n} |g| \, d\nu < \infty$, the functions $\alpha$ and $\beta$ are finite almost everywhere for their respective measures. Choose $x_0 \in \mathbb R^n$ with $\alpha(x_0)<\infty$ and $y_0 \in \mathbb R^n$ with $\beta(y_0)<\infty$. From \begin{align*} x \cdot y \leq \alpha(x)+\beta(y) \quad \text{for all } x,y \in \mathbb R^n \end{align*} we get, by fixing $x=x_0$, \begin{align*} \beta(y) \geq x_0 \cdot y-\alpha(x_0) \quad \text{for all } y \in \mathbb R^n. \end{align*} This affine lower bound is exactly what prevents the conjugate from degenerating at $x_0$. Define \begin{align*} \varphi: \mathbb R^n \to (-\infty,\infty], \qquad \varphi(x):=\sup_{y \in \mathbb R^n}\{x \cdot y-\beta(y)\}. \end{align*} For every $x,y \in \mathbb R^n$, the inequality $x \cdot y \leq \alpha(x)+\beta(y)$ gives $x \cdot y-\beta(y)\leq \alpha(x)$, and taking the supremum over $y$ gives \begin{align*} \varphi(x)\leq \alpha(x) \quad \text{for all } x \in \mathbb R^n. \end{align*} In particular, $\varphi(x_0)\leq \alpha(x_0)<\infty$. Also, because $\beta(y_0)<\infty$, the single competitor $y_0$ in the supremum gives \begin{align*} \varphi(x)\geq x \cdot y_0-\beta(y_0)>-\infty \quad \text{for all } x \in \mathbb R^n. \end{align*} Thus $\varphi$ is proper. It is convex and lower semicontinuous because it is the supremum of affine functions of $x$. Now define \begin{align*} \varphi^*: \mathbb R^n \to (-\infty,\infty], \qquad \varphi^*(y):=\sup_{x \in \mathbb R^n}\{x \cdot y-\varphi(x)\}. \end{align*} The same competitor $x_0$ gives the affine lower bound \begin{align*} \varphi^*(y)\geq x_0 \cdot y-\varphi(x_0) \quad \text{for all } y \in \mathbb R^n. \end{align*} Since $\mu$ and $\nu$ have finite second moments, every affine function on $\mathbb R^n$ is integrable with respect to $\mu$ and $\nu$. Therefore the extended integrals of $\varphi$ and $\varphi^*$ are defined from below; no expression of the form $\infty-\infty$ is being used. Because $\beta$ is proper and bounded below by an affine function, the biconjugate $\beta^{**}$ is well-defined and the Fenchel-Moreau biconjugation inequality gives $\beta^{**}\leq \beta$ pointwise. Since $\varphi=\beta^*$, we have $\varphi^*=\beta^{**}$, so $\varphi^*\leq \beta$ pointwise. Together with $\varphi\leq \alpha$, this yields \begin{align*} \int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu\leq \int_{\mathbb R^n} \alpha \, d\mu+\int_{\mathbb R^n} \beta \, d\nu. \end{align*} Fenchel's inequality for the proper lower semicontinuous convex function $\varphi$ gives \begin{align*} x \cdot y \leq \varphi(x)+\varphi^*(y) \quad \text{for all } x,y \in \mathbb R^n. \end{align*} Integrating against any coupling $\pi \in \Pi(\mu,\nu)$ gives \begin{align*} \int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi(x,y)\leq \int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu. \end{align*} Taking the supremum over all couplings gives \begin{align*} \sup_{\pi \in \Pi(\mu,\nu)}\int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi(x,y)\leq \int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu. \end{align*} On the other hand, the previous comparison with $\alpha$ and $\beta$ gives \begin{align*} \int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu\leq \int_{\mathbb R^n} \alpha \, d\mu+\int_{\mathbb R^n} \beta \, d\nu. \end{align*} It remains to identify the right-hand side. From the definitions of $\alpha$ and $\beta$, \begin{align*} \int_{\mathbb R^n} \alpha \, d\mu+\int_{\mathbb R^n} \beta \, d\nu=\frac{1}{2}\left(A-\int_{\mathbb R^n} f \, d\mu-\int_{\mathbb R^n} g \, d\nu\right), \end{align*} where \begin{align*} A:=\int_{\mathbb R^n} |x|^2 \, d\mu(x)+\int_{\mathbb R^n} |y|^2 \, d\nu(y). \end{align*} Dual optimality for $f,g$ says \begin{align*} \int_{\mathbb R^n} f \, d\mu+\int_{\mathbb R^n} g \, d\nu=\int_{\mathbb R^n \times \mathbb R^n} |x-y|^2 \, d\pi_0(x,y). \end{align*} Expanding $|x-y|^2=|x|^2+|y|^2-2x\cdot y$ and using that $\pi_0$ has marginals $\mu$ and $\nu$, we get \begin{align*} \int_{\mathbb R^n} \alpha \, d\mu+\int_{\mathbb R^n} \beta \, d\nu=\int_{\mathbb R^n \times \mathbb R^n} x\cdot y \, d\pi_0(x,y). \end{align*} Since $\pi_0$ is a coupling, this last quantity is at most the correlation supremum. Combining the three inequalities forces equality: \begin{align*} \int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu=\sup_{\pi \in \Pi(\mu,\nu)}\int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi(x,y). \end{align*}[/guided]

[step:Use complementary slackness for the Fenchel pair] Define the nonnegative measurable function \begin{align*} h: \mathbb R^n \times \mathbb R^n \to [0,\infty], \qquad h(x,y):=\varphi(x)+\varphi^*(y)-x \cdot y. \end{align*} Fenchel's inequality gives $h \geq 0$. The affine lower bounds obtained in the preceding step imply that the extended integrals of $\varphi$ and $\varphi^*$ are defined from below. Since their sum has already been proved equal to the finite correlation supremum, no undefined extended-real subtraction occurs below. Since $\pi_0$ maximizes the correlation functional by the identity in the first step, and since $(\varphi,\varphi^*)$ has the same dual value, we have \begin{align*} \int_{\mathbb R^n \times \mathbb R^n} h \, d\pi_0=\int_{\mathbb R^n} \varphi \, d\mu+\int_{\mathbb R^n} \varphi^* \, d\nu-\int_{\mathbb R^n \times \mathbb R^n} x \cdot y \, d\pi_0(x,y)=0. \end{align*} Therefore $h=0$ for $\pi_0$-almost every $(x,y)$, meaning \begin{align*} \varphi(x)+\varphi^*(y)=x \cdot y \quad \text{for } \pi_0\text{-almost every }(x,y). \end{align*} [/step]

[step:Identify the Fenchel equality set with the subdifferential graph] By the Fenchel-Young equality criterion for proper lower semicontinuous convex functions, for $x,y \in \mathbb R^n$ one has \begin{align*} \varphi(x)+\varphi^*(y)=x \cdot y \end{align*} if and only if $y \in \partial \varphi(x)$. The graph of $\partial \varphi$ is closed in $\mathbb R^n \times \mathbb R^n$, hence Borel measurable. Applying this equivalence to the $\pi_0$-almost everywhere equality obtained in the previous step gives \begin{align*} \pi_0\bigl(\{(x,y) \in \mathbb R^n \times \mathbb R^n : y \in \partial \varphi(x)\}\bigr)=1. \end{align*} Thus the optimal plan $\pi_0$ is concentrated on the graph of the convex subdifferential of the proper lower semicontinuous convex function $\varphi$, as claimed. [/step]