I recently saw a question here about bounded/unbounded preimages of a set under a harmonic function. The question asked did not seem to make sense as it was talking about harmonic functions on $\mathbb{R}$, but got me thinking about a corresponding problem for a harmonic function on $\mathbb{R}^n$, and I would appreciate it if anyone could tell me if my argument below is correct, or if there is a neater way to do it.
Suppose that $u:\mathbb{R}^n\to\mathbb{R}$ is harmonic and suppose that the set $C = u^{-1}(\{c\})$ is bounded. Since adding a constant to $u$ leaves it harmonic, we may assume that $c=0$. Take any point $x\in\mathbb{R}^n$ and consider a ball $B(x,r)$ with $r$ sufficiently large that all points of $C$ are in the interior of $B$. Now either $u>0$ or $u<0$ on $\partial B$ since $u$ is continuous and non-zero on $\partial B$. Take the case $u>0$ on $\partial B$ (the other case is similar), and recall that $u(x)$ is the average of the values of $u$ on the surface $\partial B$, so that $u(x)>0$. Since $x$ is arbitrary, it follows that $u$ is bounded below on $\mathbb{R}$, and harmonic, and hence constant.
Rephrasing: "If $u$ is a non-constant harmonic function on $\mathbb{R}$, then $u^{−1}(c)$ is unbounded."
Any comments or corrections would be gratefully received.