Let $\mu,\nu \in \mathcal{P}_2(\mathbb{R}^n)$, and assume $\mu \ll \mathcal{L}^n$. Let $\Pi(\mu,\nu)$ denote the set of Borel probability measures on $\mathbb{R}^n \times \mathbb{R}^n$ with first marginal $\mu$ and second marginal $\nu$. Then the quadratic Kantorovich problem

\begin{align*} \inf_{\pi \in \Pi(\mu,\nu)} \int_{\mathbb{R}^n \times \mathbb{R}^n} |x-y|^2 \, d\pi(x,y) \end{align*}

has a unique minimizer $\pi_*$. Moreover, there exist a proper lower semicontinuous convex function $\varphi:\mathbb{R}^n \to (-\infty,\infty]$, finite and differentiable $\mu$-a.e., and a Borel map $T:\mathbb{R}^n \to \mathbb{R}^n$ such that $T=\nabla\varphi$ $\mu$-a.e. and

\begin{align*} \pi_* = (\operatorname{id}_{\mathbb{R}^n},T)_\#\mu. \end{align*}

Equivalently, $\pi_*=(\operatorname{id}_{\mathbb{R}^n},\nabla\varphi)_\#\mu$ after choosing any Borel representative of $\nabla\varphi$ on the $\mu$-null nondifferentiability set. If $\psi:\mathbb{R}^n \to (-\infty,\infty]$ is another convex Brenier potential whose gradient has a Borel representative inducing an optimal plan, then $\nabla\psi=\nabla\varphi$ $\mu$-a.e.