Show harmonic function is constant on $\mathbb{R}^n$

I'm trying to solve the following question (this is just for practice):

If $u$ is harmonic within $\mathbb{R}^n$ with $\int_{\mathbb{R}^n}|Du|^2 dx \leq C$ for some $C > 0$, then show that u is a constant in $\mathbb{R}^n$.

I guess the idea is to somehow show that $Du = 0$ which implies $u$ is constant, or otherwise show that $u$ is bounded and thus constant by Liouville's theorem. I can't quite see how to do this though. Of course if it were on a bounded domain $U$ I know I could use the integration by parts formula $$0 = - \int_U u \Delta u dx = \int_U |Du|^2 dx - \int_{\partial U} u^2 dS$$ which would imply that $u$ is bounded ... but then Liouville's theorem doesn't apply because it's not defined on all of $\mathbb{R}^n$ (I think).

Can anyone point me in the right direction? (Also, out of curiosity, is there some kind of analogue of integration by parts for unbounded domains?)

If I understand your problem, you are trying to prove that an harmonic function on $\mathbb R^N$ whose gradient is $L^2(\mathbb R^N)$ is necessarily constant.

First of all, you should point out what definition of harmonic function you use; I will use the fact that a continuous function $u$ which satisfies the mean property is harmonic (it is equivalent to $u \in C^2$ and $\Delta u =0$). I remark also the fact that an harmonic function is necessarily $C^{\infty}$ (indeed, analytic) - you can easily see this by using mollifiers.

If $u$ is harmonic, then $\partial_{x_i}u$ is harmonic for every $i=1,\ldots , N$. It follows that $\nabla u$ has the mean property: using Cauchy-Schwarz, $$ \vert \nabla u(x) \vert =\left\vert \frac{1}{\omega_N R^N}\int_{B(x,R)}\nabla u(y)dy \right\vert \leq \frac{1}{\omega_N R^N}\int_{B(x,R)} \vert \nabla u(y) \vert dy\leq \frac{\sqrt{\omega_N R^N}}{\omega_N R^N}\left(\int_{B(x,R)}|\nabla u(y)|^2 dy\right)^{1/2}\leq\frac{||\nabla u||_2}{\sqrt{\omega_N R^N}} $$

Let $R \to + \infty$ and you get $\nabla u(x)=0$, which is the claim. Hope this helps (but I am afraid of misunderstanding the question).

Remark: if this is correct, we can substitute the hypothesis $\nabla u \in L^2(\mathbb R^N)$ with the more general one $\nabla u \in L^p(\mathbb R^N)$, for some $p \in [1,+\infty]$.

This problem is really interesting and after a long search on internet, I have found the solution of the problem: I will post here the crucial steps to prove it

I - Weitzenbock–Bochner Formula

For every harmonic function $u$ we have that $$\tag{1}\frac{1}{2}\Delta (|\nabla u|^2)=|\operatorname{Hess}(u)|^2$$

II - Another Fomurla (this is the formula for $\Delta(uv)$.

$$\tag{2}|\nabla u|\Delta (|\nabla u|)=-|\nabla (|\nabla u|)|^2+\frac{1}{2}\Delta (|\nabla u|^2)$$

III - Kato Refined Inequality

$$\tag{3}|\operatorname{Hess}(u)|^2-|\nabla (|\nabla u|)|^2\geq\frac{1}{n-1}|\nabla (|\nabla u|)|^2$$

for $x$ a.e. on the open dense subset of $\mathbb{R}^n$: $\Omega=\{x\in \mathbb{R}^n:\ |\nabla u|\neq 0\}$.

IV - Ambrosio and Xavier Proposition 2.1

Now, we can prove the statement.

By joining $(1)-(3)$, we have that $$\tag{4}|\nabla u|\Delta(|\nabla u|)\geq\frac{1}{N-1}|\nabla(|\nabla u|)|^2$$

in the sense of distributions. Because the $L^2(\mathbb{R}^n)$ energy of $u$ is bounded, we conclude by IV that $|\nabla u|$ is constant and then $\nabla u(x)=0$.

Remark: All I have done here can be generalized for Manifolds. (see Proposition 6.1). Also, maybe there is a more straightforward prove of this fact in the $\mathbb{R}^n$ case!