Prove Liouville theorem without using mean value property 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}|\nabla 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.

From the comments: is it rigorous to do it like this: $-\Delta u = 0 \implies \int_{\mathbb R^n} |\nabla u|^2 = 0$ (integrating by parts, hence $u$ is constant? It seems to easy, probably I'm missing something.
 A: Here is a more thorough write up of my comment.
Fix a non-negative smooth function $\phi$ which is identically $1$ on $B_1$ and vanishes identically outside $B_2$.  Pick $M$ so $|\Delta \phi| \leq M$.   Set $\phi_R(x)=\phi(x/R)$.  We have $|\Delta \phi_R(x)|\leq M/R^2$.
By the Bochner identity
$$
\Delta \frac{1}{2} |\nabla u|^2= \nabla u \cdot \nabla \Delta u + |\nabla^2 u|^2=|\nabla^2 u|^2
$$
We have
$$
\int_{B_R} |\nabla^2 u|^2\leq \int_{\mathbb{R}^n} \phi_R |\nabla^2 u|^2=\frac{1}{2} \int_{\mathbb{R}^n} \phi_R \Delta |\nabla u|^2\leq \frac{M}{2R^2} \int_{B_{2R}} |\nabla u|^2 \leq \frac{CM}{2 R^2}.
$$
Sending $R\to \infty$ implies $\int_{\mathbb{R}^n}|\nabla^2 u|^2=0$ so $\nabla^2 u$ vanishes identically.
This means $u(x)=\mathbf{a}\cdot x +b$ but finite energy forces $\mathbf{a}=0$.
A: Below we will prove, by purely energy methods, the following sharper statement: Let $u$ be a harmonic function which satisfies $$\liminf_{r \to \infty} \frac1{r^2} \frac{1}{|B_r|}\int_{B_r} |u|^2 = 0.$$ Then $u$ is constant.
All Liouville theorems of this sort are soft/qualitative versions of harder/more quantitative regularity statements. In this case, the Liouville statement follows from the interior gradient $L^\infty$ bound for harmonic functions (as we will see below). So another version of your question could be: Is there an energy methods proof of the interior $L^\infty$ bound for harmonic functions which in particular does not use the mean value property?
The classical energy methods proof of the pointwise estimates for harmonic functions uses only two ingredients: the Caccioppoli inequality (the most basic energy estimate), and the Sobolev inequality (because you have to get pointwise bounds from $L^2$ bounds somehow!).
The Caccioppoli inequality says that
\begin{equation*}
\int_{B_r}
|\nabla u|^2
\leq 
\frac{C}{r^2} \int_{B_{2r}}  u^2 \,.
\end{equation*}
You get this by testing the equation with $\varphi^2u$ where $\varphi$ is an appropriate cutoff function.
Iterating the Caccioppoli inequality many times, we get
\begin{equation*}
\int_{B_r} |\nabla^k u|^2
\leq 
\frac{C}{r^{2k}} \int_{B_{2^kr}} u^2 \,.
\end{equation*}
If $k > 1+\frac d2$, then we can the Sobolev inequality to get
\begin{equation*}
\| \nabla u \|_{L^\infty(B_r)}^2 
\leq 
\frac{C}{r^2} 
\frac{1}{|B_r|}
\int_{B_{2^kr}} u^2\,.
\end{equation*}
Now send $r\to \infty$ along a good subsequence for the liminf, and you discover that $\nabla u=0$, which means $u$ is constant.
