Solution (source code)

= Solution

Let $\Phi=\varphi_\varepsilon$ and let $\Phi^*$ denote its <convex conjugate>. Its <Fenchel–Young gap>
$$
H(x,y)=\Phi(x)+\Phi^*(y)-x\cdot y
$$
is nonnegative by the <Fenchel–Young inequality>, and the hypothesis says $\int H\,d\pi_\varepsilon\leq\varepsilon$.

First justify the integrability needed for the certificate. Finite second <moments> and the <Cauchy-Schwarz inequality> imply, for every <transport plan> $\pi$,
$$
\int|x\cdot y|\,d\pi\leq
\left(\int|x|^2\,d\mu\right)^{1/2}
\left(\int|y|^2\,d\nu\right)^{1/2}<\infty.
$$
Since $\Phi\in L^1(\mu)$ and $H\in L^1(\pi_\varepsilon)$, the identity $\Phi^*(y)=H(x,y)-\Phi(x)+x\cdot y$ proves $\Phi^*\in L^1(\nu)$ by the <marginal distribution> property. In particular, no subtraction of infinite integrals is being used.

Set
$$
u(x)=\tfrac12|x|^2-\Phi(x),\qquad
v(y)=\tfrac12|y|^2-\Phi^*(y),\qquad
D=\int u\,d\mu+\int v\,d\nu.
$$
These are integrable <Kantorovich potentials>. The <Fenchel–Young inequality> gives $u(x)+v(y)\leq\tfrac12|x-y|^2$. More precisely, for every <transport plan> $\pi$,
$$
\mathbb K(\pi)=D+\int H\,d\pi\geq D,
\qquad\mathbb K(\pi_\varepsilon)=D+\int H\,d\pi_\varepsilon\leq D+\varepsilon.
$$
Taking the infimum in the first inequality therefore proves
$$
\boxed{\mathbb K(\pi_\varepsilon)\leq\inf_{\pi\in\Pi(\mu,\nu)}\mathbb K(\pi)+\varepsilon.}
$$
The factor $1/2$ in the quadratic cost is essential for this exact gap identity.