Liouville Theorem for Harmonic Functions 
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.
 A: In fact, Poisson integral formula implies the following lemma.

(Harnack’s Inequality) Suppose $u$ is harmonic in $B_R(x_0)$ and $u\geq0$.Then there holds
  $$\left(\frac{R}{R+r}\right)^{n-2}\frac{R-r}{R+r}u(x_0)\leq u(x)\leq \left(\frac{R}{R-r}\right)^{n-2}\frac{R+r}{R-r}u(x_0) $$
  where $r=|x-x_0|<R.$

Proof. Without loss of a generality, we may assume $x_0=0$ and $u\in C(\overline{B}_R).$ By Poisson integral formula, we get
$$u(x)=\frac{R^2-|x|^2}{n\alpha(n)R}\int_{\partial B_R}\frac{u(y)}{|x-y|^n}\,dS$$
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)\, dS\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)\, dS.
\end{align} 
Mean value property implies
$$u(0)=\frac{1}{n\alpha(n)R^{n-1}}\int_{\partial B_R}u(y)\, dS.$$
That finish the proof.
Hence Liouville theorem is a corollary of the lemma.
In particular, We assume $u\geq 0$in $\mathbb{R^n}$.  Take any point $x\in\mathbb{R^n}$ and apply 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 implies $u(x)=u(0)$ by letting $R\to +\infty. $

If we don’t use Poisson formula, we will give another proof. Mean value formula implies the following lemma.

Suppose $u\in C(\overline{B}_R(x_0))$ is a nonnegative harmonic function in $B_R=B_R(x_0)$. Then there holds
  $$|Du(x_0)|\leq \frac{n}{R}u(x_0).$$

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.
A: I hope the E. Nelson proof of the Liouville Theorem for Harmonic Functions will be useful here, see references [1] and [2]. 
Mean-Value Property: If $u$ is a harmonic function on the ball $B(P, r)$ with the ball volume $V$, then $u(P)$ equals average of $u$ over ball $B(P, r)$
$$
  u(P) = \frac{1}{V} \int_B u \, dV
$$
Liouville Theoreme: A bounded harmonic function on $R^n$ is constant
Suppose $u$ is a harmonic function on $R^n$ bounded by a constant $C$. Consider two points $P$, $Q$ and consider two balls with the given points as centers and of equal radius $r$: $B_1(P, r)$ and $B_2(Q, r)$. Denote by $V$ the volume of each ball.
By the mean-value property we get:
$$
  | u(P) - u(Q) | = \frac{1}{V} \left| \int_{B_1} u \, dV- \int_{B_2} u \, dV \right| \leq C \frac{V_{B_1 \setminus B_2}}{V},
$$
where $V_{B_1 \setminus B_2}$ is the symmetric difference between two balls.
If $r \to \infty$, then the last expression $\to 0$ because the value of $V_{B_1 \setminus B_2}$ comparing to the full volume $V$ becomes really small. So $u(P) = u(Q)$ and $u$ is a constant.
References
[1] Nelson, Edward. "A proof of Liouville’s theorem." Proc. AMS. Vol. 12. 1961.
[2] Bounded harmonic function is constant
A: You're basically there already! Two hints:


*

*Note that $a>0$ in your argument was arbitrary. What does that mean for the inequality
$$
\left|\frac{\partial u(0)}{\partial x_i}\right| \le \frac{n}{a}\Vert u \Vert_{L^\infty}?
$$

*Note that you can replace the origin $x_0=0$ in your argument by arbitrary $x_0 \in \mathbb{R}^n$.

