Let $\Omega$ be an open subset of $\mathbb{R}^N$ and $u:\Omega\to \mathbb{R}$ be a harmonic function. I need to prove that, for every $x_0\in\Omega$, the function
$$
\rho\mapsto \frac{1}{\rho^N}\int_{|x-x_0|<\rho} |\nabla u|^2 \; dx
$$
is non-decreasing on $(0,+\infty)$.
I proved that for any $u\in C^2(\Omega,\mathbb{R})$ (not necessarily harmonic) we have that
$$
\frac{d}{d\rho} \left( \frac{1}{\rho^{N-1}}\int_{|x-x_0|=\rho} u \; d\sigma \right) = \frac{1}{\rho^{N-1}}\int_{|x-x_0|<\rho} \Delta u\; dx,
$$
and from this, if $u$ is harmonic, that the function
$$
\rho \mapsto \frac{1}{\rho^N}\int_{|x-x_0|<\rho} u \; dx
$$
is constant on $(0,+\infty)$.
In the book Partial Differential Equations by Evans, Section 8.6, he proves that the function
$$
\rho \mapsto \frac{1}{\rho^{N-2}}\int_{|x-x_0|<\rho} |\nabla u|^2 \; dx
$$
is non-decreasing, but I can't see (if possible) how to apply the ideas given there to this case.
Best Answer
Forgive me for writing down part of this in Einstein summation notation, i'm just not used enough to vector calculus to be able to write this out in reasonable length otherwise.
If I'm not mistaken, $\left|\nabla u \right|^2$ is a subharmonic function, i.e. $\Delta\left|\nabla u \right|^2>0$ on $\Omega$. Indeed, \begin{align*} \Delta\left|\nabla u \right|^2 = \nabla \cdot \nabla (\nabla u \cdot \nabla u) = \partial^k\partial_k \left(\partial^ju \partial_j u\right) = 2\partial^k\left(\left(\partial_k\partial^j u \right)\partial_j u\right) = \\ 2\left(\left(\partial^j \partial^k\partial_ku \right)\partial_j u + \left(\partial_k\partial^j u \right)\left(\partial^k\partial_j u \right)\right) = 2\left(\left(\partial^j \Delta u \right)\partial_j u + tr(H_uH_u^T\right)) = 2 tr(H_uH_u^T) > 0, \end{align*} where $H_u$ denotes the Hessian of $u$.
Then $F(\rho) := \frac{1}{\rho^{N-1}} \int_{|x-x_0| = \rho} |\nabla u|^2dx$ is a non-decreasing function by the OP's remark. A change of variables yields
$$\frac{1}{\rho^{N}}\int_{|x-x_0| < \rho} |\nabla u|^2dx = \frac{1}{\rho^{N}} \int_0^\rho r^{n-1}\frac{1}{r^{n-1}}\int_{|y-x_0| = r}|\nabla u|^2dy dr = \\ \frac{1}{\rho^{N}} \int_0^\rho r^{n-1} F(r) dr = \int_0^1 s^n F(s\rho) ds.$$
We can now differentiate inside the integral to obtain
$$\frac{d}{d\rho} \int_0^1 s^n F(s\rho) ds = \int_0^1 s^{n+1} F'(s\rho) ds > 0,$$
finishing the proof.