This thread is meant to record a question that I feel interesting during my self-study. I'm very happy to receive your suggestion and comments. See: SE blog: Answer own Question and MSE meta: Answer own Question.
Let $X,Y$ be topological spaces. Let $\mathcal P(X), \mathcal P(Y)$ be the spaces of all Borel probability measures on $X,Y$ respectively. Let $c: X \times Y \to [0, +\infty]$. Fix $\mu \in \mathcal P(X)$ and $\nu \in \mathcal P(Y)$.
$\Pi(\mu, \nu)$ is the set of $\pi \in \mathcal P(X \times Y)$ such that for all measurable subsets $A \subset X$ and $B \subset Y$, $$ \pi (A \times Y) = \mu (A) \quad \text{and} \quad \pi (X \times B) = \nu (B). $$
$\Phi_{c}$ is the set of all $(\varphi, \psi) \in L_1 (\mu) \times L_1 (\nu)$ satisfying $$ \varphi(x)+\psi(y) \leq c(x, y) $$ for $\mu$-a.e. $x \in X$ and $\nu$-a.e. $y \in Y$.
For $\pi \in \mathcal P(X \times Y)$ and $(\varphi, \psi) \in L_1 (\mu) \times L_1 (\nu)$, let $$ \mathbb K (\pi) := \int_{X \times Y} c d \pi \quad \text{and} \quad \mathbb J(\varphi, \psi) := \int_{X} \varphi d \mu+\int_{Y} \psi d \nu . $$
Theorem: Let $X = Y = \mathbb R^d$ and assume that $\mu, \nu$ both have finite second moment. Then $\pi$ is a minimizer of $\mathbb K$ over $\Pi(\mu, \nu)$ with cost $c(x,y) := \frac{1}{2} |x-y|^2$ if and only if there is $\varphi \in L_1 (\mu)$ convex l.s.c. such that $y \in \partial f (x)$ for $\pi$-a.e. $(x, y) \in X \times Y$. Moreover, the pair $(\varphi, \varphi^*)$ is a minimizer of $\mathbb J$ over $\tilde \Phi$ where $\varphi^*$ is the convex conjugate and $$ \tilde \Phi := \{(\tilde \varphi, \tilde \psi) \in L_1(\mu) \times L_1 (\nu) \mid \tilde \varphi (x) + \tilde \psi (y) \ge \langle x, y \rangle\}. $$