[proofplan]
We use the stationarity of $u$ with respect to compactly supported domain variations. Testing the stationarity identity against radial vector fields centered at $x_0$ gives a distributional identity relating the derivative of the ball energy to the radial part of $|\nabla u|^2$. After rewriting this identity as an almost-everywhere derivative formula for $\Theta_u(x_0,t)$, integration in the radius and the polar integration formula give the stated annular identity. The monotonicity follows because the right-hand side is an integral of a nonnegative function.
[/proofplan]
[step:Write stationarity as a stress-energy identity]
Since $u \in W^{1,2}(U;N)$ and $N$ is isometrically embedded in $\mathbb{R}^q$, we view $u$ as an $\mathbb{R}^q$-valued Sobolev map with $u(x)\in N$ for $\mathcal{L}^m$-a.e. $x \in U$. Let
\begin{align*}
|\nabla u|^2 := \sum_{\alpha=1}^m |\partial_{x_\alpha}u|^2
\end{align*}
denote the Euclidean Hilbert-Schmidt norm of the weak differential.
The [stationarity identity for harmonic maps](/page/Stationary%20Harmonic%20Map), equivalently the vanishing first variation of the stress-energy tensor under domain variations, means that for every compactly supported $C^1$ vector field
\begin{align*}
X:U \to \mathbb{R}^m,
\end{align*}
one has
\begin{align*}
\int_U
\left(
|\nabla u|^2 \operatorname{div}X
-
2\sum_{\alpha,\beta=1}^m
\partial_{x_\alpha}u \cdot \partial_{x_\beta}u\,
\partial_{x_\alpha}X_\beta
\right)
\,d\mathcal{L}^m(x)
=0.
\end{align*}
Here $\cdot$ is the Euclidean [inner product](/page/Inner%20Product) in $\mathbb{R}^q$.
[/step]
[step:Test stationarity with a smooth radial vector field]
Fix a function $\eta \in C_c^1((0,R))$. Define the vector field $X_\eta: B(x_0,R) \to \mathbb{R}^m$ by $X_\eta(x)=\eta(|x-x_0|)(x-x_0)$, and extend it by $0$ outside $B(x_0,R)$ inside $U$. Since $\eta$ has compact support in $(0,R)$, this extension is compactly supported in $U$.
Write $s:=|x-x_0|$. Let $S^{m-1}:=\{w\in\mathbb{R}^m:|w|=1\}$ denote the Euclidean unit sphere, and let $\delta_{\alpha\beta}$ denote the Kronecker delta, equal to $1$ when $\alpha=\beta$ and equal to $0$ otherwise. Define the unit radial map $e_r:B(x_0,R)\setminus\{x_0\}\to S^{m-1}$ by
\begin{align*}
e_r(x)=\frac{x-x_0}{|x-x_0|}.
\end{align*}
For $x \ne x_0$,
\begin{align*}
\partial_{x_\alpha}(X_\eta)_\beta(x)
=
\eta(s)\delta_{\alpha\beta}
+
\eta'(s)\frac{(x_\alpha-(x_0)_\alpha)(x_\beta-(x_0)_\beta)}{s},
\end{align*}
and therefore
\begin{align*}
\operatorname{div}X_\eta(x)=m\eta(s)+s\eta'(s).
\end{align*}
Also,
\begin{align*}
\sum_{\alpha,\beta=1}^m
\partial_{x_\alpha}u\cdot \partial_{x_\beta}u\,
\frac{(x_\alpha-(x_0)_\alpha)(x_\beta-(x_0)_\beta)}{s}
=
s\left|\frac{\partial u}{\partial r}(x)\right|^2.
\end{align*}
Substituting these identities into the stationarity formula gives
\begin{align*}
\int_{B(x_0,R)}
\left((m-2)\eta(s)+s\eta'(s)\right)
|\nabla u|^2
\,d\mathcal{L}^m(x)
=
2\int_{B(x_0,R)}
s\eta'(s)
\left|\frac{\partial u}{\partial r}(x)\right|^2
\,d\mathcal{L}^m(x).
\end{align*}
[guided]
The purpose of the radial test field is to make the tangential and radial parts of the energy separate. We choose the map $X_\eta: B(x_0,R) \to \mathbb{R}^m$ defined by $X_\eta(x)=\eta(|x-x_0|)(x-x_0)$, where $\eta \in C_c^1((0,R))$. The compact support of $\eta$ away from $0$ and $R$ ensures that the extension of $X_\eta$ by $0$ is a compactly supported $C^1$ vector field in $U$, so it is an admissible domain variation in the [stationarity identity for harmonic maps](/page/Stationary%20Harmonic%20Map).
Set $s:=|x-x_0|$. Let $\delta_{\alpha\beta}$ denote the Kronecker delta: $\delta_{\alpha\beta}=1$ if $\alpha=\beta$ and $\delta_{\alpha\beta}=0$ otherwise. The codomain $S^{m-1}$ of the radial direction is the Euclidean unit sphere $\{w\in\mathbb{R}^m:|w|=1\}$. Differentiating the $\beta$-th component of $X_\eta$ gives
\begin{align*}
\partial_{x_\alpha}(X_\eta)_\beta(x)
=
\partial_{x_\alpha}\left(\eta(s)(x_\beta-(x_0)_\beta)\right)
=
\eta'(s)\frac{x_\alpha-(x_0)_\alpha}{s}(x_\beta-(x_0)_\beta)
+
\eta(s)\delta_{\alpha\beta}.
\end{align*}
Taking the trace over $\alpha=\beta$ gives
\begin{align*}
\operatorname{div}X_\eta(x)
=
\sum_{\alpha=1}^m
\left(
\eta(s)
+
\eta'(s)\frac{(x_\alpha-(x_0)_\alpha)^2}{s}
\right)
=
m\eta(s)+s\eta'(s).
\end{align*}
The second term in the stationarity identity contains the quadratic form of $\nabla u$ applied to the radial direction. Indeed, by the definition
\begin{align*}
\frac{\partial u}{\partial r}(x)
=
\sum_{\alpha=1}^m
\frac{x_\alpha-(x_0)_\alpha}{s}\,\partial_{x_\alpha}u(x),
\end{align*}
we have
\begin{align*}
s\left|\frac{\partial u}{\partial r}(x)\right|^2
=
s
\sum_{\alpha,\beta=1}^m
\frac{x_\alpha-(x_0)_\alpha}{s}
\frac{x_\beta-(x_0)_\beta}{s}
\partial_{x_\alpha}u(x)\cdot \partial_{x_\beta}u(x).
\end{align*}
Cancelling one factor of $s$ in the preceding expression gives
\begin{align*}
s\left|\frac{\partial u}{\partial r}(x)\right|^2
=
\sum_{\alpha,\beta=1}^m
\partial_{x_\alpha}u(x)\cdot \partial_{x_\beta}u(x)\,
\frac{(x_\alpha-(x_0)_\alpha)(x_\beta-(x_0)_\beta)}{s}.
\end{align*}
Substituting the computed derivative and divergence into stationarity yields
\begin{align*}
0
=
\int_{B(x_0,R)}
\left(
|\nabla u|^2(m\eta(s)+s\eta'(s))
-
2\eta(s)|\nabla u|^2
-
2s\eta'(s)\left|\frac{\partial u}{\partial r}\right|^2
\right)
\,d\mathcal{L}^m(x).
\end{align*}
Collecting the $|\nabla u|^2$ terms gives
\begin{align*}
0
=
\int_{B(x_0,R)}
\left((m-2)\eta(s)+s\eta'(s)\right)
|\nabla u|^2
\,d\mathcal{L}^m(x)
-
2\int_{B(x_0,R)}
s\eta'(s)
\left|\frac{\partial u}{\partial r}\right|^2
\,d\mathcal{L}^m(x).
\end{align*}
Rearranging gives the desired radial stationarity identity.
[/guided]
[/step]
[step:Convert the radial identity into a derivative formula for the scaled energy]
Define the energy distribution function $E:(0,R) \to [0,\infty)$ by
\begin{align*}
E(t)=\int_{B(x_0,t)} |\nabla u|^2\,d\mathcal{L}^m(x).
\end{align*}
Let $\mathcal{H}^{m-1}$ denote $(m-1)$-dimensional [Hausdorff measure](/page/Hausdorff%20Measure) on hypersurfaces in $\mathbb{R}^m$, and let $\mathcal{L}^1$ denote one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure). Define the radial boundary energy function $A:(0,R) \to [0,\infty]$, for those radii for which the boundary integral is defined, by
\begin{align*}
A(t)=\int_{\partial B(x_0,t)} \left|\frac{\partial u}{\partial r}(x)\right|^2\,d\mathcal{H}^{m-1}(x).
\end{align*}
The distance map $s:B(x_0,R)\to (0,R)$, $s(x)=|x-x_0|$, is Lipschitz and satisfies $|\nabla s|=1$ for $\mathcal{L}^m$-a.e. $x\ne x_0$. Since $|\nabla u|^2\in L^1(B(x_0,R))$, the [Coarea Formula](/theorems/23) applied to $s$ gives absolute continuity of $E$ and, for $\mathcal{L}^1$-a.e. $t\in(0,R)$,
\begin{align*}
E'(t)=\int_{\partial B(x_0,t)} |\nabla u|^2\,d\mathcal{H}^{m-1}(x).
\end{align*}
Define $B:(0,R)\to[0,\infty]$ for those radii for which the boundary integral is defined by
\begin{align*}
B(t)=\int_{\partial B(x_0,t)} |\nabla u|^2\,d\mathcal{H}^{m-1}(x).
\end{align*}
By the coarea formula, $B\in L^1((0,R),\mathcal{L}^1)$, $E'(t)=B(t)$ for $\mathcal{L}^1$-a.e. $t\in(0,R)$, and
\begin{align*}
\int_0^R A(t)\,d\mathcal{L}^1(t)
=
\int_{B(x_0,R)} \left|\frac{\partial u}{\partial r}(x)\right|^2\,d\mathcal{L}^m(x)
\leq
\int_{B(x_0,R)} |\nabla u|^2\,d\mathcal{L}^m(x)<\infty.
\end{align*}
Thus $A\in L^1((0,R),\mathcal{L}^1)$. Applying the coarea formula to the radial identity from the previous step gives, for every $\eta\in C_c^1((0,R))$,
\begin{align*}
\int_0^R \left((m-2)\eta(t)+t\eta'(t)\right)B(t)\,d\mathcal{L}^1(t)
=
2\int_0^R t\eta'(t)A(t)\,d\mathcal{L}^1(t).
\end{align*}
Define $G\in L^1_{\mathrm{loc}}((0,R),\mathcal{L}^1)$ by
\begin{align*}
G(t)=t(B(t)-2A(t))+(2-m)E(t)
\end{align*}
for those $t$ for which $A(t)$ and $B(t)$ are defined, and choose any measurable value elsewhere. Since $E$ is absolutely continuous and $E'=B$ for $\mathcal{L}^1$-a.e. $t$, [integration by parts](/theorems/2098) for compactly supported $C^1$ functions gives, for every $\eta\in C_c^1((0,R))$,
\begin{align*}
\int_0^R E(t)\eta'(t)\,d\mathcal{L}^1(t)= -\int_0^R \eta(t)B(t)\,d\mathcal{L}^1(t).
\end{align*}
By the definition of $G$,
\begin{align*}
\int_0^R G(t)\eta'(t)\,d\mathcal{L}^1(t)=\int_0^R t\eta'(t)(B(t)-2A(t))\,d\mathcal{L}^1(t)+(2-m)\int_0^R E(t)\eta'(t)\,d\mathcal{L}^1(t).
\end{align*}
Substituting the [integration by parts](/theorems/210) identity into this expression gives
\begin{align*}
\int_0^R G(t)\eta'(t)\,d\mathcal{L}^1(t)=\int_0^R t\eta'(t)(B(t)-2A(t))\,d\mathcal{L}^1(t)+(m-2)\int_0^R \eta(t)B(t)\,d\mathcal{L}^1(t).
\end{align*}
The radial identity above states exactly that the right-hand side is zero, so
\begin{align*}
\int_0^R G(t)\eta'(t)\,d\mathcal{L}^1(t)=0.
\end{align*}
Hence the [distributional derivative](/page/Distributional%20Derivative) of $G$ is zero on $(0,R)$, so $G$ is equal to a constant $c\in\mathbb{R}$ for $\mathcal{L}^1$-a.e. $t\in(0,R)$.
It remains to identify this constant. Since $|\nabla u|^2\in L^1(B(x_0,R),\mathcal{L}^m)$, the absolute continuity of the [Lebesgue integral](/page/Lebesgue%20Integral) gives $E(t)\to0$ as $t\downarrow0$. Also $A+B\in L^1((0,R),\mathcal{L}^1)$, so there exists a sequence $t_j\downarrow0$ such that $t_j(A(t_j)+B(t_j))\to0$ and such that $G(t_j)=c$ for every $j$. Passing to the limit in
\begin{align*}
c=G(t_j)=t_j(B(t_j)-2A(t_j))+(2-m)E(t_j)
\end{align*}
gives $c=0$. Therefore, for $\mathcal{L}^1$-a.e. $t\in(0,R)$,
\begin{align*}
tB(t)+(2-m)E(t)=2tA(t).
\end{align*}
Using $B(t)=E'(t)$ for $\mathcal{L}^1$-a.e. $t$, this becomes
\begin{align*}
tE'(t)+(2-m)E(t)=2tA(t).
\end{align*}
Therefore, for $\mathcal{L}^1$-a.e. $t\in(0,R)$,
\begin{align*}
\frac{d}{dt}\left(t^{2-m}E(t)\right)=2t^{2-m}A(t).
\end{align*}
[/step]
[step:Integrate the derivative formula over the annulus]
Let $0<\rho<r<R$. Since $t \mapsto t^{2-m}E(t)$ is absolutely continuous on $[\rho,r]$, integrating the derivative formula over $(\rho,r)$ gives
\begin{align*}
r^{2-m}E(r)-\rho^{2-m}E(\rho)=2\int_\rho^r t^{2-m}A(t)\,d\mathcal{L}^1(t).
\end{align*}
We apply the [Coarea Formula](/theorems/23) to the Lipschitz radius map $s:B(x_0,R)\to(0,R)$, $s(x)=|x-x_0|$. Since $|\nabla s|=1$ for $\mathcal{L}^m$-a.e. $x\ne x_0$, the measure decomposition is $d\mathcal{L}^m(x)=d\mathcal{H}^{m-1}(x)\,d\mathcal{L}^1(t)$ on the level sets $\partial B(x_0,t)$, and the domain $\{x:\rho<s(x)<r\}$ is exactly $B(x_0,r)\setminus \overline{B}(x_0,\rho)$ up to an $\mathcal{L}^m$-null boundary. Hence
\begin{align*}
2\int_\rho^r t^{2-m}A(t)\,d\mathcal{L}^1(t)=2\int_{B(x_0,r)\setminus B(x_0,\rho)} |x-x_0|^{2-m}\left|\frac{\partial u}{\partial r}(x)\right|^2\,d\mathcal{L}^m(x).
\end{align*}
Substituting the definition of $E$ gives
\begin{align*}
r^{2-m}\int_{B(x_0,r)} |\nabla u|^2\,d\mathcal{L}^m(x)-\rho^{2-m}\int_{B(x_0,\rho)} |\nabla u|^2\,d\mathcal{L}^m(x)=2\int_{B(x_0,r)\setminus B(x_0,\rho)} |x-x_0|^{2-m}\left|\frac{\partial u}{\partial r}(x)\right|^2\,d\mathcal{L}^m(x).
\end{align*}
[/step]
[step:Read monotonicity from the nonnegative annular term]
For $\mathcal{L}^m$-a.e. $x \in B(x_0,R)\setminus\{x_0\}$,
\begin{align*}
|x-x_0|^{2-m}\left|\frac{\partial u}{\partial r}(x)\right|^2 \geq 0.
\end{align*}
Hence the annular integral in the formula is nonnegative. Therefore, whenever $0<\rho<r<R$,
\begin{align*}
\Theta_u(x_0,r)-\Theta_u(x_0,\rho)\geq 0.
\end{align*}
Thus $t \mapsto \Theta_u(x_0,t)$ is nondecreasing on $(0,R)$.
[/step]
