[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*}[/step]

[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: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]