How can I prove the following Liouville theorem without using the mean value property?
If $u$ is harmonic on $\mathbb{R}^n$ and $\int_{\mathbb{R}^n}|D u|^2 dx \leq C$ for some $C > 0$, then $u$ is constant.
The proof that I know indeed uses the mean value property for harmonic functions.
Let me first say that the statement $$ \int_{\mathbb{R}^n} \vert Du \vert^2 d x = -\int_{\mathbb{R}^n} u \Delta u d x$$ which is in the comments is false - integration by parts only works on bounded sets. As a simple counter example to the above claim take $u(x) = x_1$.
One way to prove this for $n>2$ without the Mean Value Formula makes use of a monotonicity formula. This is kind of similar to the Mean Value Formula in spirit (you can think of the MVF as a monotonicity formula where the monotone quantity is in fact constant).
Let $B_r$ be a ball of radius $r$. Then $$\phi : r\mapsto \frac{1}{r^{n-2}} \int_{B_r} \vert Du \vert^2 d x $$ is nondecreasing.
Observe that \begin{align*} \frac{d}{dr} \bigg ( \frac{1}{r^{n-2}} \int_{B_r} \vert Du \vert^2 d x\bigg ) &= - \frac{n-2}{r^{n-1}} \int_{B_r} \vert Du \vert^2 d x + \frac 1 {r^{n-2}} \int_{\partial B_r} \vert Du \vert^2 d \mathcal{H}^{n-1}. \end{align*} It is possible to show that for harmonic functions \begin{equation}(n-2)\int_{B_r} \vert Du \vert^2 d x = r \int_{\partial B_r} \vert Du \vert^2 - 2 \bigg ( \frac{\partial u}{\partial \nu }\bigg )^2d \mathcal{H}^{n-1}. \label{eq:2} \end{equation} (To prove this calculate $ \int_{B_r} \mathrm{div} (x \vert Du \vert^2) dx $ in two different ways: first by using the divergence theorem and secondly by expanding $\mathrm{div} (x \vert Du \vert^2)$ then applying integration by parts to one of the terms. Another way to prove this is using Neother's Theorem, see Chapter 8.6.2 in Partial Differnetial Equations by Lawrence C. Evans.) Hence, it follows $$\frac{d}{dr} \bigg ( \frac{1}{r^{n-2}} \int_{B_r} \vert Du \vert^2 d x\bigg ) = \frac 2 {r^{n-2}} \int_{\partial B_r} \bigg ( \frac{\partial u}{\partial \nu }\bigg )^2d \mathcal{H}^{n-1} \geqslant 0.$$ However, by our assumption $$ 0 \leqslant \phi(r) \leqslant \frac{C}{r^{n-2}} \to 0$$ as $r \to \infty$ provided $n>2$. Hence, $\phi(r) = 0$ for all $r >0$, so $u$ is a constant. The only issue with this is it doesn't apply to the case $n=2$.