[proofplan]
Choose a smooth cutoff function $\zeta$ supported in $B(x_0,R)$, equal to $1$ on $B(x_0,r)$, and with gradient bounded by a constant times $(R-r)^{-1}$. Since $u$ is only locally $H^1$, first extend the weak harmonic identity from smooth compactly supported test functions to the admissible Sobolev [test function](/page/Test%20Function) $\zeta^2 u$. Expanding $\nabla(\zeta^2 u)$ produces the energy term $\zeta^2|\nabla u|^2$ and a cross term. The cross term is controlled by Cauchy-Schwarz and [Young's inequality](/theorems/244), and the resulting estimate is restricted to $B(x_0,r)$ because $\zeta = 1$ there.
[/proofplan]
[step:Choose a cutoff adapted to the two concentric balls]
Write $B_R := B(x_0,R)$ and $B_r := B(x_0,r)$. Choose a cutoff function $\zeta \in C_c^\infty(B_R)$ such that $0 \leq \zeta \leq 1$ on $B_R$, $\zeta = 1$ on $B_r$, and
\begin{align*}
|\nabla \zeta(x)| \leq \frac{A_n}{R-r}
\end{align*}
for every $x \in B_R$, where $A_n > 0$ depends only on $n$. Since $B_R \subset U$ and $u \in H^1_{\mathrm{loc}}(U)$, the restrictions $u|_{B_R}$ and $\nabla u|_{B_R}$ belong to $L^2(B_R)$.
[/step]
[step:Extend the weak harmonic identity to the Sobolev test function $\zeta^2u$]
Define the linear functional
\begin{align*}
L: H^1_0(B_R) &\to \mathbb{R}
\end{align*}
by
\begin{align*}
L(v) = \int_{B_R} \nabla u(x) \cdot \nabla v(x)\, d\mathcal{L}^n(x).
\end{align*}
This is well-defined and continuous because $\nabla u \in L^2(B_R)$ and Cauchy-Schwarz gives
\begin{align*}
|L(v)| \leq \|\nabla u\|_{L^2(B_R)}\|\nabla v\|_{L^2(B_R)}
\end{align*}
for every $v \in H^1_0(B_R)$. If $\psi \in C_c^\infty(B_R)$, then $\psi \in C_c^\infty(U)$, so the weak harmonicity hypothesis gives $L(\psi)=0$. Since $H^1_0(B_R)$ is the closure of $C_c^\infty(B_R)$ in the $H^1(B_R)$ norm, continuity implies $L(v)=0$ for every $v \in H^1_0(B_R)$.
The product $\zeta^2u$ belongs to $H^1_0(B_R)$ because $\zeta \in C_c^\infty(B_R)$ and $u \in H^1(B_R)$. Therefore
\begin{align*}
\int_{B_R} \nabla u(x) \cdot \nabla(\zeta^2u)(x)\, d\mathcal{L}^n(x) = 0.
\end{align*}
[guided]
The weak formulation is stated only for smooth compactly supported test functions, while the natural test function for the energy estimate is $\zeta^2u$, which is generally not smooth. We therefore extend the identity by continuity.
Define
\begin{align*}
L: H^1_0(B_R) &\to \mathbb{R}
\end{align*}
by
\begin{align*}
L(v) = \int_{B_R} \nabla u(x) \cdot \nabla v(x)\, d\mathcal{L}^n(x).
\end{align*}
This integral is finite because $\nabla u \in L^2(B_R)$ and $\nabla v \in L^2(B_R)$. More precisely, Cauchy-Schwarz gives
\begin{align*}
|L(v)| \leq \|\nabla u\|_{L^2(B_R)}\|\nabla v\|_{L^2(B_R)}.
\end{align*}
Thus $L$ is a continuous linear functional on $H^1_0(B_R)$.
Now take any $\psi \in C_c^\infty(B_R)$. Since $B_R \subset U$, the same function $\psi$ is also an element of $C_c^\infty(U)$ after extending it by zero outside $B_R$. Hence the weak harmonicity assumption gives
\begin{align*}
L(\psi) = \int_{B_R} \nabla u(x) \cdot \nabla \psi(x)\, d\mathcal{L}^n(x) = \int_U \nabla u(x) \cdot \nabla \psi(x)\, d\mathcal{L}^n(x) = 0.
\end{align*}
Because $H^1_0(B_R)$ is defined as the closure of $C_c^\infty(B_R)$ in $H^1(B_R)$, and because $L$ is continuous, it follows that $L(v)=0$ for every $v \in H^1_0(B_R)$.
Finally, $\zeta^2u \in H^1_0(B_R)$: multiplication by the smooth compactly supported function $\zeta^2$ preserves $H^1(B_R)$ and gives compact support in $B_R$. Therefore this particular Sobolev function is an admissible test function, so
\begin{align*}
\int_{B_R} \nabla u(x) \cdot \nabla(\zeta^2u)(x)\, d\mathcal{L}^n(x) = 0.
\end{align*}
[/guided]
[/step]
[step:Expand the Sobolev product rule and isolate the energy term]
The Sobolev product rule gives
\begin{align*}
\nabla(\zeta^2u) = \zeta^2\nabla u + 2\zeta u\nabla\zeta
\end{align*}
as an identity in $L^2(B_R;\mathbb{R}^n)$. Substituting this into the previous identity yields
\begin{align*}
0 = \int_{B_R} \zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) + 2\int_{B_R}\zeta(x)u(x)\nabla u(x)\cdot \nabla\zeta(x)\, d\mathcal{L}^n(x).
\end{align*}
Therefore
\begin{align*}
\int_{B_R} \zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) = -2\int_{B_R}\zeta(x)u(x)\nabla u(x)\cdot \nabla\zeta(x)\, d\mathcal{L}^n(x).
\end{align*}
Taking absolute values gives
\begin{align*}
\int_{B_R} \zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq 2\int_{B_R}\zeta(x)|u(x)|\,|\nabla u(x)|\,|\nabla\zeta(x)|\, d\mathcal{L}^n(x).
\end{align*}
[/step]
[step:Control the cross term by Young's inequality and absorb half of the energy]
For each $x \in B_R$, apply the elementary Young inequality $2ab \leq \frac{1}{2}a^2 + 2b^2$ with
\begin{align*}
a = \zeta(x)|\nabla u(x)|
\end{align*}
and
\begin{align*}
b = |u(x)|\,|\nabla\zeta(x)|.
\end{align*}
This gives
\begin{align*}
2\zeta(x)|u(x)|\,|\nabla u(x)|\,|\nabla\zeta(x)| \leq \frac{1}{2}\zeta(x)^2|\nabla u(x)|^2 + 2|u(x)|^2|\nabla\zeta(x)|^2.
\end{align*}
Integrating over $B_R$ with respect to $\mathcal{L}^n$ and combining with the previous estimate, we obtain
\begin{align*}
\int_{B_R} \zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq \frac{1}{2}\int_{B_R}\zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) + 2\int_{B_R}|u(x)|^2|\nabla\zeta(x)|^2\, d\mathcal{L}^n(x).
\end{align*}
Subtracting the first term on the right-hand side from both sides yields
\begin{align*}
\int_{B_R} \zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq 4\int_{B_R}|u(x)|^2|\nabla\zeta(x)|^2\, d\mathcal{L}^n(x).
\end{align*}
[guided]
The term
\begin{align*}
2\int_{B_R}\zeta(x)|u(x)|\,|\nabla u(x)|\,|\nabla\zeta(x)|\, d\mathcal{L}^n(x)
\end{align*}
contains one factor of the quantity we want to estimate, namely $\zeta|\nabla u|$, and one error factor, namely $|u||\nabla\zeta|$. The standard way to separate such a product is Young's inequality.
For each fixed $x \in B_R$, set
\begin{align*}
a = \zeta(x)|\nabla u(x)|
\end{align*}
and
\begin{align*}
b = |u(x)|\,|\nabla\zeta(x)|.
\end{align*}
Young's inequality in the form $2ab \leq \frac{1}{2}a^2 + 2b^2$ gives
\begin{align*}
2\zeta(x)|u(x)|\,|\nabla u(x)|\,|\nabla\zeta(x)| \leq \frac{1}{2}\zeta(x)^2|\nabla u(x)|^2 + 2|u(x)|^2|\nabla\zeta(x)|^2.
\end{align*}
Integrating this pointwise inequality over $B_R$ with respect to $\mathcal{L}^n$ gives
\begin{align*}
2\int_{B_R}\zeta(x)|u(x)|\,|\nabla u(x)|\,|\nabla\zeta(x)|\, d\mathcal{L}^n(x) \leq \frac{1}{2}\int_{B_R}\zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) + 2\int_{B_R}|u(x)|^2|\nabla\zeta(x)|^2\, d\mathcal{L}^n(x).
\end{align*}
Combining this with the previous cross-term estimate gives
\begin{align*}
\int_{B_R} \zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq \frac{1}{2}\int_{B_R}\zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) + 2\int_{B_R}|u(x)|^2|\nabla\zeta(x)|^2\, d\mathcal{L}^n(x).
\end{align*}
The coefficient $\frac{1}{2}$ is the reason for choosing this version of Young's inequality: it lets us absorb half of the desired energy term into the left-hand side. Subtracting that term from both sides yields
\begin{align*}
\int_{B_R} \zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq 4\int_{B_R}|u(x)|^2|\nabla\zeta(x)|^2\, d\mathcal{L}^n(x).
\end{align*}
[/guided]
[/step]
[step:Use the cutoff bounds to obtain the Caccioppoli estimate on $B(x_0,r)$]
Because $\zeta = 1$ on $B_r$ and $0 \leq \zeta \leq 1$ on $B_R$,
\begin{align*}
\int_{B_r}|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq \int_{B_R}\zeta(x)^2|\nabla u(x)|^2\, d\mathcal{L}^n(x).
\end{align*}
Using the estimate from the previous step and the cutoff gradient bound,
\begin{align*}
\int_{B_r}|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq 4\int_{B_R}|u(x)|^2|\nabla\zeta(x)|^2\, d\mathcal{L}^n(x).
\end{align*}
Since $|\nabla\zeta(x)|^2 \leq A_n^2(R-r)^{-2}$ for every $x \in B_R$, we obtain
\begin{align*}
\int_{B_r}|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq \frac{4A_n^2}{(R-r)^2}\int_{B_R}|u(x)|^2\, d\mathcal{L}^n(x).
\end{align*}
Defining $C_n := 4A_n^2$ gives
\begin{align*}
\int_{B(x_0,r)}|\nabla u(x)|^2\, d\mathcal{L}^n(x) \leq \frac{C_n}{(R-r)^2}\int_{B(x_0,R)}|u(x)|^2\, d\mathcal{L}^n(x),
\end{align*}
which is the desired estimate.
[/step]