Removable singularity for harmonic functions There seems to be a removable singularity theorem for harmonic functions that somewhat goes like:

Let $\Omega\subset \Bbb R^n$ be an open region and $x_0\in\Omega$, $u$ is harmonic in $\Omega-\{x_0\}$. Now if we have $u(x)=o(\Phi(x-x_0))$ around $x_0$ then there exists $\tilde u$ harmonic in all of $\Omega$ and $u=\tilde u$ in $\Omega-\{x_0\}$. 

Here $\Phi$ is the fundamental solution to $\Delta u=0$ in $\Bbb R^n$, and, since dropping the normalising constant doesn't affect our statement above, you may just take $\Phi=-\log|x-x_0|$ in 2D and $=1/|x-x_0|^{n-2}$ in 3D or higher. 
I have two questions concerning this theorem:
1). I can't find a satisfactory proof. So far I only found this one, which however doesn't look okay to me because there seems to be an abuse of $\epsilon-\delta$ language: when the author concluded the proof by letting $\epsilon\to 0$, he seemed to forget the previous part of his proof depended on $\delta$ which depended on $\epsilon$. So is there any other available proof?
2). I think the statement of the theorem is unnecessarily complicated. If $x_0$ is really removable for $u$, then I think it's necessary $u$ is bounded around $x_0$. And then boundedness clearly implies the condition that $u(x)=o(\Phi(x-x_0))$. So why not just say "if $u$ is bounded around $x_0$ then it's a removable singularity"? Why all the fuss about comparing $u$ against $\Phi$?
 A: (1) Let $B(x_0, r) =  \{x: |x-x_0|< r\}$. Fix $r$ such that $\overline{B(x_0, r)}\subset \Omega$. Let $v$ be the harmonic function in the $B(x_0, r)$ that agrees with $u$ on $\partial B(x_0, r)$, such a function exists as it's given by the Poisson integral. For $\epsilon>0$, consider 
$$
w(x) = u(x)-v(x)-\epsilon \Phi(x, x_0) 
$$
which is harmonic in $B(x_0, r)\setminus \{x_0\}$ and tends to $-\infty$ as $x\to x_0$. By the maximum principle (you may want to remove a tiny neighborhood of $x_0$ before applying it), 
$$
\sup \{w(x): x\in B(x_0, r)\setminus \{x_0\}\}) = -\epsilon \sup_{|x-x_0|}\Phi(x, x_0)
$$ 
Letting $\epsilon\to 0$, we obtain that $u(x)-v(x) \le 0$ for every $x\in B(x_0, r)\setminus \{x_0\}$ (think of such $x$ as fixed while $\epsilon\to 0$).
The same argument applies to $w(x) = v(x)-u(x)-\epsilon \Phi(x, x_0)$ and shows that  $v(x)-u(x) \le 0$ for every $x\in B(x_0, r)\setminus \{x_0\}$. Conclusion: $v$ provides the required harmonic extension of $u$.
(2) There is a benefit in having weaker assumptions in a theorem, because a weaker statement may be easier to verify. You noted that replacing "bounded" by "$o(\Phi)$" does not expand the set of functions to which the theorem applies; but it can make it easier to prove that a particular function falls in this set. The various forms of Phragmén–Lindelöf principle are another example of this. E.g., if a  function $f$ is holomorphic in the the strip  $|\operatorname{Im} z|<\pi/2$, its boundary values are bounded by $1$, and $$f(z)< \exp(A\exp(c|\operatorname{Re} z|))$$ for some finite $A$ and $c<1$, then $f$ is actually bounded by $1$ in the strip. Why such a fussy condition, when the function is in fact bounded? It's so that we don't have to prove that it's bounded in order to apply the theorem.
