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.