If $u$ is bounded and harmonic in $\mathbb{R}^n$, then $u$ is constant
For any twice differentiable function $u$ defined on an open subset $\Omega$ we have $$u(x)=\int\limits_\Omega G(x,y)\Delta u(y)dy+\int\limits_{\partial\Omega}\frac{\partial G}{\partial n}(x,y)u(y)dS_y$$ where $G$ is the corresponding Green's function.
We take the region to be $B_a(0)$ and take $G(x,y)=\Phi(x-y)-\Phi(\frac{|x|}{a}|x^*-y|)$ where $\Phi$ is the fundamental solution to the Laplace equation $\Delta u=0$ and $x^*$ is the inverse of the point $x$ w.r.t. the ball $B_a(0)$.It turns out then $$u(x)=\int\limits_{\partial\Omega}H(x,y)u(y)dS_y $$ where $$H(x,y)=\frac{a^2-|x|^2}{aw_n}\frac{1}{|x-y|^n}$$ $w_n$ being the surfacae integral of the unit sphere in $\mathbb{R}^n$. Partially differentiating $H$ we get $$|H_i(0,y)|\leq \frac{n}{w_na^n}$$ which gives $$|\frac{\partial u}{\partial x_i}(0)|\leq \frac{n}{a}\|u\|_\infty$$ Now I want something like $$|\frac{\partial u}{\partial x_i}(x)|\leq \frac{n}{a}\|u\|_\infty \forall x$$ to conclude $u$ is constant. I am not sure how to do so. Any help is appreciated.
In fact, Poisson integral formula implies the following lemma.
Proof. Without loss of a generality, we may assume $x_0=0$ and $u\in C(\overline{B}_R).$ By the Poisson integral formula, we have $$u(x)=\frac{R^2-|x|^2}{n\alpha(n)R}\int_{\partial B_R}\frac{u(y)}{|x-y|^n}\mathrm{~d}S$$ where $\alpha(n)$ is the volume of unit sphere.
Since $R-r\leq |x-y|\leq R+r$ with $|y|=R$, we have \begin{align*} \frac{1}{n\alpha(n)R}\cdot\frac{R-|x|}{R+|x|}&\left(\frac{1}{R+|x|}\right)^{n-2}\int_{\partial B_R} u(y)\mathrm{~d}S\leq u(x)\\&\leq \frac{1}{n\alpha(n)R}\cdot\frac{R+|x|}{R-|x|}\left(\frac{1}{R-|x|}\right)^{n-2} \int_{\partial B_R} u(y)\mathrm{~d}S. \end{align*} Then the Mean Value Property gives us that $$u(0)=\frac{1}{n\alpha(n)R^{n-1}}\int_{\partial B_R}u(y)\mathrm{~d}S,$$ which completes the proof.
Hence Liouville theorem is a corollary of the lemma.
In particular, We assume $u\geq 0$ in $\mathbb{R^n}$. Taking any point $x\in\mathbb{R^n}$ and applying the lemma to any ball $B_R(0)$ with $R > |x|$, we obtain $$\left(\frac{R}{R+r}\right)^{n-2}\frac{R-r}{R+r}u(0)\leq u(x)\leq \left(\frac{R}{R-r}\right)^{n-2}\frac{R+r}{R-r}u(0), $$ which yields $u(x)=u(0)$ by letting $R\to +\infty$.
If we don’t use Poisson's formula, we shall give another proof. The Mean Value Formula implies the following lemma.
Proof. Since u is smooth in $B_R$, we know $\Delta (D_{x_i}u)=0$, that is $D_{x_i}u$ is also harmonic in $B_R$. Hence $D_{x_i}u$ satisfies mean value formula. Then by divergence theorem and the nonnegativeness of $u$ we have \begin{align}|D_{x_i}u(x_0)|=\left|\frac{1}{\alpha(n)R^n}\int_{B_R}D_{x_i}u(y)\,dy\right|&=\left|\frac{1}{\alpha(n)R^n}\int_{\partial B_R} u(y) v_{x_i} \,dS\right|\\&\leq \frac{1}{\alpha(n)R^n}\int_{\partial B_R} u(y) \,dS\\&=\frac{n}{R}u(x_0) \end{align}
where in the last equality we used the mean value property.
Then letting $R\to +\infty$, we get $Du(x)=0$ for all $x\in\mathbb{R^n},$ which implies the Liouville theorem.
All above are from Qin Han and Fanghua Lin’s elliptic partial differential equations.