[proofplan]
The transport polytope is a nonempty compact subset of the finite-dimensional space $\mathbb R^{m\times n}$, so the unregularized linear cost attains a minimum. We compare the entropic value with this minimum from above by testing the entropic problem on an unregularized minimizer. For the lower bound, we use the elementary uniform estimate $x(\log x-1)\ge -1$ on $[0,1]$, which bounds the entropy contribution from below independently of the plan. The two inequalities squeeze $\operatorname{OT}_\varepsilon(a,b;C)$ to the unregularized optimal value as $\varepsilon\downarrow0$.
[/proofplan]
[step:Verify compactness of the transport polytope and define the two cost functionals]
Define the matrix $P_{\mathrm{prod}}=(P_{\mathrm{prod},ij})\in[0,\infty)^{m\times n}$ by
\begin{align*}
P_{\mathrm{prod},ij}:=a_i b_j
\end{align*}
for $i\in\{1,\dots,m\}$ and $j\in\{1,\dots,n\}$. Then
\begin{align*}
\sum_{j=1}^n P_{\mathrm{prod},ij}=a_i\sum_{j=1}^n b_j=a_i
\end{align*}
for every $i\in\{1,\dots,m\}$, and
\begin{align*}
\sum_{i=1}^m P_{\mathrm{prod},ij}=b_j\sum_{i=1}^m a_i=b_j
\end{align*}
for every $j\in\{1,\dots,n\}$. Hence $P_{\mathrm{prod}}\in\Pi(a,b)$, so $\Pi(a,b)$ is nonempty.
Every $P\in\Pi(a,b)$ satisfies $0\le P_{ij}\le \sum_{k=1}^n P_{ik}=a_i\le 1$ for all $i,j$, so $\Pi(a,b)\subseteq[0,1]^{m\times n}$. Since $\Pi(a,b)$ is also cut out by finitely many affine hyperplanes in the finite-dimensional Euclidean space $\mathbb R^{m\times n}$, it is closed and bounded. By the Heine-Borel theorem in finite-dimensional Euclidean space, $\Pi(a,b)$ is compact.
Define the linear cost functional
\begin{align*}
L:\Pi(a,b)&\to\mathbb R,\qquad P\mapsto \sum_{i=1}^m\sum_{j=1}^n C_{ij}P_{ij}.
\end{align*}
Since $L$ is continuous on the compact set $\Pi(a,b)$, the extreme value theorem implies that $L$ attains its minimum. Define
\begin{align*}
M:=\min_{P\in\Pi(a,b)}L(P),
\end{align*}
and choose $P_\ast\in\Pi(a,b)$ such that $L(P_\ast)=M$.
Define $h:[0,1]\to\mathbb R$ by $h(x)=x(\log x-1)$ for $x\in(0,1]$ and $h(0)=0$. The limit $\lim_{x\downarrow0}x\log x=0$ shows that $h$ is continuous on $[0,1]$. Define the entropy correction functional
\begin{align*}
H:\Pi(a,b)&\to\mathbb R,\qquad P\mapsto \sum_{i=1}^m\sum_{j=1}^n h(P_{ij}).
\end{align*}
Since $H$ is a finite sum of continuous coordinate functions on $\Pi(a,b)$, the functional $P\mapsto L(P)+\varepsilon H(P)$ is continuous on the compact set $\Pi(a,b)$ for every $\varepsilon>0$. The extreme value theorem therefore implies that its minimum is attained. With this notation,
\begin{align*}
\operatorname{OT}_\varepsilon(a,b;C)=\min_{P\in\Pi(a,b)}\bigl(L(P)+\varepsilon H(P)\bigr)
\end{align*}
for every $\varepsilon>0$.
[/step]
[step:Obtain the upper bound by testing on an unregularized minimizer]
For every $\varepsilon>0$, the point $P_\ast$ is an admissible competitor in the minimization defining $\operatorname{OT}_\varepsilon(a,b;C)$. Therefore
\begin{align*}
\operatorname{OT}_\varepsilon(a,b;C)\le L(P_\ast)+\varepsilon H(P_\ast)=M+\varepsilon H(P_\ast).
\end{align*}
The number $H(P_\ast)$ is finite because each entry of $P_\ast$ lies in $[0,1]$ and $h$ is finite-valued on $[0,1]$. Taking the limit superior as $\varepsilon\downarrow0$ gives
\begin{align*}
\limsup_{\varepsilon\downarrow0}\operatorname{OT}_\varepsilon(a,b;C)\le M.
\end{align*}
[guided]
Recall that $L:\Pi(a,b)\to\mathbb R$ denotes the linear cost functional, $H:\Pi(a,b)\to\mathbb R$ denotes the entropy correction functional, $M=\min_{P\in\Pi(a,b)}L(P)$, and $P_\ast\in\Pi(a,b)$ is chosen so that $L(P_\ast)=M$. The upper bound uses the simplest possible comparison: instead of trying to identify the minimizer of the entropic problem, we plug in a plan that is already optimal for the unregularized problem.
The plan $P_\ast$ belongs to $\Pi(a,b)$, so it is allowed in the minimization defining $\operatorname{OT}_\varepsilon(a,b;C)$. Since a minimum is bounded above by the value of the objective at any admissible competitor, we obtain
\begin{align*}
\operatorname{OT}_\varepsilon(a,b;C)\le L(P_\ast)+\varepsilon H(P_\ast).
\end{align*}
The choice of $P_\ast$ gives $L(P_\ast)=M$, hence
\begin{align*}
\operatorname{OT}_\varepsilon(a,b;C)\le M+\varepsilon H(P_\ast).
\end{align*}
Now we must check that the coefficient $H(P_\ast)$ is a genuine finite real number. Every entry of a transport plan is nonnegative, and since all entries sum to $1$, each entry lies in $[0,1]$. The function $h:[0,1]\to\mathbb R$ is finite-valued under the convention $h(0)=0$, so the finite sum
\begin{align*}
H(P_\ast)=\sum_{i=1}^m\sum_{j=1}^n h((P_\ast)_{ij})
\end{align*}
is finite. Therefore the term $\varepsilon H(P_\ast)$ converges to $0$ as $\varepsilon\downarrow0$, and we get
\begin{align*}
\limsup_{\varepsilon\downarrow0}\operatorname{OT}_\varepsilon(a,b;C)\le M.
\end{align*}
[/guided]
[/step]
[step:Bound the entropy correction uniformly from below]
For $x\in(0,1]$, the derivative of the map $x\mapsto x(\log x-1)$ is $\log x$, which is nonpositive on $(0,1]$. Hence $h$ is decreasing on $(0,1]$. Since $h(1)=-1$ and $h(0)=0$, it follows that
\begin{align*}
h(x)\ge -1
\end{align*}
for every $x\in[0,1]$.
Every $P\in\Pi(a,b)$ satisfies $P_{ij}\in[0,1]$ for all $i,j$. Therefore
\begin{align*}
H(P)=\sum_{i=1}^m\sum_{j=1}^n h(P_{ij})\ge -mn
\end{align*}
for every $P\in\Pi(a,b)$.
[/step]
[step:Derive the lower bound from the unregularized minimum]
Let $\varepsilon>0$ and let $P\in\Pi(a,b)$. By definition of $M$ as the minimum of $L$ on $\Pi(a,b)$, we have
\begin{align*}
L(P)\ge M.
\end{align*}
By the uniform entropy lower bound from the previous step,
\begin{align*}
H(P)\ge -mn.
\end{align*}
Multiplying the second inequality by $\varepsilon>0$ and adding gives
\begin{align*}
L(P)+\varepsilon H(P)\ge M-\varepsilon mn.
\end{align*}
Taking the minimum over all $P\in\Pi(a,b)$ yields
\begin{align*}
\operatorname{OT}_\varepsilon(a,b;C)\ge M-\varepsilon mn.
\end{align*}
Consequently,
\begin{align*}
\liminf_{\varepsilon\downarrow0}\operatorname{OT}_\varepsilon(a,b;C)\ge M.
\end{align*}
[/step]
[step:Combine the two one-sided bounds]
The previous steps give
\begin{align*}
M\le \liminf_{\varepsilon\downarrow0}\operatorname{OT}_\varepsilon(a,b;C)\le \limsup_{\varepsilon\downarrow0}\operatorname{OT}_\varepsilon(a,b;C)\le M.
\end{align*}
Thus the lower and upper limits are equal to $M$, and therefore
\begin{align*}
\operatorname{OT}_\varepsilon(a,b;C)\longrightarrow M
\end{align*}
as $\varepsilon\downarrow0$. Recalling the definition of $M$, this is precisely
\begin{align*}
\operatorname{OT}_\varepsilon(a,b;C)\longrightarrow \min_{P\in\Pi(a,b)}\sum_{i=1}^m\sum_{j=1}^n C_{ij}P_{ij}.
\end{align*}
[/step]
