I am trying to prove the Liouville Theorem using the Mean Value Property for Harmonic Functions. The question is:
Let $u$ be a Harmonic Function and $|u|$ is bounded. Prove that $u$ is constant. Interpret this when $u$ represents the temperature at a point $\textbf x$ in the plane.
Here is my $\textit{attempt} $ at a solution.
We begin by showing: $$|\partial_{x_i} u(\textbf x_0)|\leq \frac{3|n_i|M}{r}$$ Noting that $\partial_{x_i} u(\textbf x_0)$ is Harmonic. Applying the MVP: $$|\partial_{x_i} u(\textbf x_0)| = \bigg|\frac{3}{4\pi r^3}\iiint_{B(\textbf x_0,r)}\partial_{x_i} u(\vec x) d\vec x\bigg|$$ Using the Divergence Theorem: $$\bigg|\frac{3}{4\pi r^3}\iiint_{B(\textbf x_0,r)}\partial_{x_i} u(\vec x) d\vec x\bigg| = \frac{3}{4\pi r^3}\bigg|\iint_{\partial B(\textbf x_0, r)} u(\vec x)\cdot n_i dS_{\vec x}\bigg|$$ Where $n_i$ denotes the $i$-th component of the unit normal vector. Now using the Cauchy Inequality: $$\frac{3}{4\pi r^3}\bigg|\iint_{\partial B(\textbf x_0, r)} u(\vec x)\cdot n_i dS_{\vec x}\bigg| \leq \frac{3}{4\pi r^3}\iint_{\partial B(\textbf x_0,r)}| u(\vec x)|\cdot|n_i|dS_{\vec x}$$ The above yields: $$\frac{3}{4\pi r^3}\iint_{\partial B(\textbf x_0,r)}| u(\vec x)|\cdot|n_i|dS_{\vec x}\leq \frac{3}{4\pi r^3}\cdot 4\pi r^2|n_i|\cdot M = \frac{3|n_i|\cdot M}{r}$$ Now, applying this to: $$\bigg|\sum_{i=1}^3\partial x_i u(\textbf x_0)\bigg|\leq \frac{C M}{r}$$ Where $C$ is some constant. Now, since $u$ and $\partial_{x_i} u$ are both Harmonic on $\mathbb R^3$, taking $r\rightarrow\infty$ yields that: $$|\nabla u(\textbf x_0)| =0$$ Thus we conclude that $u$ is a constant.
Am I missing anything? Also if someone could explain how to interpret this as a temperature it would be appreciated. Thanks!
Since the argument seems all right to me, I'll just add a few words on two generalizations (as partially requested in the comments).
First, the argument actually shows that if $u$ has polynomial growth, i.e. $|u(x)| \lesssim (1+|x|)^k$ for some $k$, then $|\nabla u| \lesssim (1+|x|)^{k-1}$. Iterating this, we get $$ |u(x)| \lesssim (1+|x|)^k \quad \Rightarrow \quad D^{k+1} u \equiv 0 \quad \Rightarrow \quad u \text{ is a polynomial of degree } \le k. $$
Second, it is enough to assume that $u$ is bounded from below (or from above).
Without loss of generality, assume that $u \ge 0$ in $\mathbb R^n$ (otherwise consider $\pm u + c$). Take any two points $x,y \in \mathbb R^n$ and denote $d = |x-y|$. For any $r>0$ we have $B(x,r) \subseteq B(y,r+d)$ and so $$ \omega_n r^n u(x) = \int_{B(x,r)} u \le \int_{B(y,r+d)}u = \omega_n(r+d)^n u(y) $$ by mean value property. Thus, we obtained the inequality $$ u(x) \le \left( 1 + \frac d r \right)^n u(y) $$ valid for any $r>0$, in consequence $u(x) \le u(y)$. Since $x,y$ are arbitrary (and in particular the order doesn't matter), $u$ is constant.