Let $U: \mathbb C \to \mathbb R$ be a harmonic function such that $U(z)\ge 0 \ \forall z\in \mathbb C$ ; then $U$ is constant.
I don't know about simply connected domain. Can we do it without using that?
Edit: The more elementary the proof, the better. I'm working my way through a complex analysis workbook, and by this exercise, the substantial things covered are Lioville's theorem, Maximum modulus theorem.
On
Another, similar, way to go:
Again, since $U: \Bbb C \to \Bbb R$ is defined everywhere, it has a global harmonic conjugate $V: \Bbb C \to \Bbb R$ and then
$F(z) = U(z) + iV(z) \tag 1$
is entire. Thus
$e^{F(z)} = e^{U(z) + iV(z)}, \tag 2$
is also entire, as is
$e^{-F(z)} = e^{-U(z) - iV(z)}; \tag 3$
we have
$\vert e^{-F(z)} \vert = \vert e^{-U(z) - iV(z)} \vert = \vert e^{-U(z)}e^{-iV(z)} \vert = \vert e^{-U(z)} \vert \vert e^{-iV(z)} \vert; \tag 4$
now
$\vert e^{-iV(z)} \vert = 1, \tag 5$
since $V(z)$ is real, thus
$\vert e^{-F(z)} \vert = \vert e^{-U(z)} \vert \le 1, \tag 6$
since $U(z) \ge 0$ implies $-U(z) \le 0$. Thus $e^{-F(z)}$ is a bounded entire function, hence constant by Liouville's theorem. Thus $-F(z)$, and hence $F(z)$, must also be constant. Thus $U(z)$ is constant as well.
On
Liouville Type (Theorem's) Suppose $u:\mathbb R^d\to \mathbb R$ be a nonegaitve harmonic function Then, $u$ is constant.
Proof $u\ge 0$ then $\inf_{\mathbb R^d} u<\infty$. Hence we set $$v= u-\inf_{\mathbb R^d} u$$ and
\begin{split} \begin{cases}\Delta v=0\\ v\ge 0\\ \inf_{\mathbb R^d} v =0\end{cases}\end{split}
Let $\varepsilon>0,$ and $y\in\mathbb R^d$ then there exists $x_\varepsilon$ such that,
$$ v(x_\varepsilon)\le \inf_{\mathbb R^d} v+\varepsilon$$
Let $R>\max(|x_\varepsilon|,|y|)+1$ therefore, $x_\varepsilon,y\in B_R(0)$ and
$$ B_R(y)\subset B_{3R}(x_\varepsilon) $$
By Mean value property (The mean value is true in any dimension for harmonics functions),
\begin{split} v(y) &=& \frac{1}{|B_R(y)|}\int_{B_R(y)} v(z) dz\\ &=& \frac{3^d}{|B_{3R}(x_\varepsilon)|}\int_{B_R(y)} v(z) dz \\&\le&\frac{3^d}{|B_{3R}(x_\varepsilon)|}\int_{B_{3R}(x_\varepsilon)} v(z) dz \\&= & 3^d v(x_\varepsilon) \end{split} This leads to,
$$ v(y)\le 3^d v(x_\varepsilon)\le 3^d(\inf_{\mathbb R^d} v+\varepsilon)~~\forall y\in \mathbb R^d$$
That is $$ \sup_{\mathbb R^d} v < 3^d \varepsilon $$
Since $\inf_{\mathbb R^d} v=0$ ~~ but $\varepsilon>0$ was arbitrarily chosen, letting $\varepsilon \to 0$ we get
$$ 0\le \sup_{\mathbb R^d} v \le 0$$
i.e $v = 0$ or $u= \inf_{\mathbb R^d} u $
The proof given HERE relies on the fact that a harmonic function on a simply connected domain is the real part of a holomorphic function on that simply connected domain (For a proof, SEE HERE). Hence, we have $u(z)=\text{Re}(f(z))$, for an enire function $f$.
Next, since we assume that $u(z)=\text{Re}(f(z))>0$ then $f(z)+1=\left(u(z)+1\right)+i\,\text{Im}(f(z))$ can never be zero. Hence, the function $g$ given by
$$g(z)=\frac{f(z)-1}{f(z)+1}$$
is entire.
Finally, since
$$\left|\frac{f(z)-1}{f(z)+1}\right|=\sqrt{\frac{(u(z)-1)^2+\left(\text{Im}(f(z))\right)^2}{(u(z)+1)^2+\left(\text{Im}(f(z))\right)^2}}<1$$
is bounded, then Liouville's Theorem guarantees that $g$ is a constant. In turn, this implies that $f$ is a constant, which implies that $u$ is a constant.
And we are done!