[proofplan]
Choose a smooth cutoff $\zeta$ which equals $1$ on $B(x_0,r)$, vanishes outside $B(x_0,R)$, and has gradient bounded by a constant times $(R-r)^{-1}$. Use $\phi=\zeta^2u$ as a [test function](/page/Test%20Function) in the weak formulation, justified by approximating in $H^1_0(B(x_0,R))$. Expanding $\nabla(\zeta^2u)$ gives the energy term $\int \zeta^2|\nabla u|^2$, while the cutoff-error term and forcing term are controlled by pointwise Young inequalities. Absorbing the smaller energy contribution yields the desired estimate on $B(x_0,r)$ because $\zeta=1$ there.
[/proofplan]
[step:Choose a cutoff adapted to the two balls]
Set $\delta := R-r > 0$ and $B_R := B(x_0,R)$. Choose a one-dimensional cutoff $\theta: \mathbb{R} \to [0,1]$ with $\theta \in C^\infty(\mathbb{R})$, $\theta(t)=1$ for $t \le 0$, $\theta(t)=0$ for $t \ge 1$, and $|\theta'(t)| \le C_0$ for every $t \in \mathbb{R}$, where $C_0>0$ is an absolute constant. Define the function $\zeta: U \to [0,1]$ by setting, for each $x \in U$,
\begin{align*}
\zeta(x) := \theta\left(\frac{|x-x_0|-r}{\delta}\right).
\end{align*}
Because $\theta$ is constant on $(-\infty,0]$, this radial composition is smooth at $x_0$ as well as away from $x_0$. It satisfies $\zeta \in C_c^\infty(B_R)$, $\zeta=1$ on $B(x_0,r)$, and the chain rule gives
\begin{align*}
|\nabla \zeta(x)| \le \frac{C_0}{\delta}
\end{align*}
for every $x \in U$. Since $u \in H^1_{\mathrm{loc}}(U)$ and $\zeta \in C_c^\infty(B_R)$, the product $\zeta^2u$ belongs to $H^1_0(B_R)$ and has compact support in $U$.
Let $K := \operatorname{supp}\zeta$, so $K \subset B_R$ is compact and hence $K \Subset U$. Choose an [open set](/page/Open%20Set) $W$ with $K \Subset W \Subset U$. Since $u \in H^1_{\mathrm{loc}}(U)$ and $f \in L^2_{\mathrm{loc}}(U)$, we have $u \in H^1(W)$ and $f \in L^2(W)$. The product $\zeta^2u$ belongs to $H^1_0(W)$ because $\zeta^2u \in H^1(W)$ and its support is contained in $K \Subset W$. By the definition of [Sobolev spaces with zero trace](/page/Sobolev%20Space), choose a sequence $(\phi_k)_{k=1}^{\infty}$ in $C_c^\infty(W)$ such that $\phi_k \to \zeta^2u$ in $H^1(W)$, and extend each $\phi_k$ by zero to an element of $C_c^\infty(U)$. Applying the weak formulation to each extended $\phi_k$ and passing to the limit by the $L^2(W)$ [Cauchy-Schwarz Inequality](/theorems/432) gives
\begin{align*} \int_{B_R} \nabla u \cdot \nabla(\zeta^2u) \, d\mathcal{L}^n = \int_{B_R} f \zeta^2u \, d\mathcal{L}^n. \end{align*}
[/step]
[step:Expand the energy identity and isolate the main gradient term]
Using the [Product Rule for Weak Derivatives](/theorems/3098),
\begin{align*} \nabla(\zeta^2u) = \zeta^2\nabla u + 2\zeta u\nabla \zeta. \end{align*}
as an $L^2(B_R;\mathbb{R}^n)$ function. Substituting this into the weak formulation gives
\begin{align*}
\int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n
+
2\int_{B_R} \zeta u \nabla u \cdot \nabla \zeta \, d\mathcal{L}^n
=
\int_{B_R} f\zeta^2u \, d\mathcal{L}^n.
\end{align*}
Taking absolute values of the two non-main terms and applying the Euclidean [Cauchy-Schwarz Inequality](/theorems/432) pointwise to $\nabla u \cdot \nabla\zeta$ yields
\begin{align*}
\int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n
\le
2\int_{B_R} \zeta |u| |\nabla u| |\nabla \zeta| \, d\mathcal{L}^n
+
\int_{B_R} |f| \zeta^2 |u| \, d\mathcal{L}^n.
\end{align*}
[guided]
We choose the cutoff so that the weak equation measures the gradient only where we need it. Set $\delta := R-r > 0$ and $B_R := B(x_0,R)$. Choose $\zeta: U \to [0,1]$ with $\zeta \in C_c^\infty(B_R)$, $\zeta=1$ on $B(x_0,r)$, and $|\nabla \zeta(x)| \le C_0\delta^{-1}$ for every $x \in U$, where $C_0>0$ depends only on the dimension.
Let $K := \operatorname{supp}\zeta$ and choose an open set $W$ with $K \Subset W \Subset U$. Since $u \in H^1_{\mathrm{loc}}(U)$ and $f \in L^2_{\mathrm{loc}}(U)$, we have $u \in H^1(W)$ and $f \in L^2(W)$. The function $\zeta^2u$ belongs to $H^1_0(W)$, so by the definition of [Sobolev spaces with zero trace](/page/Sobolev%20Space) there are functions $\phi_k \in C_c^\infty(W)$ with $\phi_k \to \zeta^2u$ in $H^1(W)$. Extending $\phi_k$ by zero to $U$, applying the weak formulation, and passing to the limit using the $L^2(W)$ [Cauchy-Schwarz Inequality](/theorems/432) gives
\begin{align*} \int_{B_R} \nabla u \cdot \nabla(\zeta^2u) \, d\mathcal{L}^n = \int_{B_R} f\zeta^2u \, d\mathcal{L}^n. \end{align*}
By the [Product Rule for Weak Derivatives](/theorems/3098),
\begin{align*} \nabla(\zeta^2u) = \zeta^2\nabla u + 2\zeta u\nabla \zeta. \end{align*}
Substitution and the Euclidean [Cauchy-Schwarz Inequality](/theorems/432), applied pointwise to $\nabla u \cdot \nabla\zeta$, give
\begin{align*} \int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n \le 2\int_{B_R} \zeta |u| |\nabla u| |\nabla \zeta| \, d\mathcal{L}^n + \int_{B_R} |f|\zeta^2|u| \, d\mathcal{L}^n. \end{align*}
Now apply [Young's Inequality](/theorems/244) in the pointwise form
\begin{align*}
2ab \le \frac{1}{2}a^2+2b^2
\end{align*}
with $a:=\zeta|\nabla u|$ and $b:=|u||\nabla\zeta|$. This yields
\begin{align*} \frac{1}{2}\int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n \le 2\int_{B_R} |u|^2|\nabla \zeta|^2 \, d\mathcal{L}^n + \int_{B_R} |f|\zeta^2|u| \, d\mathcal{L}^n. \end{align*}
Using $|\nabla\zeta| \le C_0\delta^{-1}$ and multiplying by $2$ gives
\begin{align*} \int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n \le 4C_0^2\delta^{-2}\int_{B_R} |u|^2 \, d\mathcal{L}^n + 2\int_{B_R} |f|\zeta^2|u| \, d\mathcal{L}^n. \end{align*}
Apply [Young's Inequality](/theorems/244) again, now in the form $2ab \le a^2+b^2$ with $a:=\delta |f|\zeta$ and $b:=\delta^{-1}|u|\zeta$. Since $0\le \zeta\le 1$,
\begin{align*} \int_{B_R} |f|\zeta^2|u| \, d\mathcal{L}^n \le \frac{\delta^2}{2}\int_{B_R} |f|^2 \, d\mathcal{L}^n + \frac{1}{2\delta^2}\int_{B_R} |u|^2 \, d\mathcal{L}^n. \end{align*}
Combining the preceding two estimates gives
\begin{align*} \int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n \le (4C_0^2+1)\delta^{-2}\int_{B_R} |u|^2 \, d\mathcal{L}^n + \delta^2\int_{B_R} |f|^2 \, d\mathcal{L}^n. \end{align*}
Finally, $\zeta=1$ on $B(x_0,r)$ and the integrand is non-negative, so restricting from $B_R$ to $B(x_0,r)$ and recalling $\delta=R-r$ gives the claimed Caccioppoli inequality with $C:=\max\{4C_0^2+1,1\}$.
[/guided]
[/step]
[step:Absorb the cutoff-gradient error into the energy]
Apply [Young's Inequality](/theorems/244) in the pointwise form
\begin{align*}
2ab \le \frac{1}{2}a^2 + 2b^2
\end{align*}
with $a := \zeta|\nabla u|$ and $b := |u||\nabla \zeta|$.
Then
\begin{align*}
2\zeta |u| |\nabla u| |\nabla \zeta|
\le
\frac{1}{2}\zeta^2|\nabla u|^2
+
2|u|^2|\nabla \zeta|^2.
\end{align*}
Hence
\begin{align*}
\frac{1}{2}
\int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n
\le
2\int_{B_R} |u|^2|\nabla \zeta|^2 \, d\mathcal{L}^n
+
\int_{B_R} |f| \zeta^2 |u| \, d\mathcal{L}^n.
\end{align*}
Using $|\nabla\zeta| \le C_0\delta^{-1}$, we obtain
\begin{align*}
\int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n
\le
4C_0^2\delta^{-2}
\int_{B_R} |u|^2 \, d\mathcal{L}^n
+
2\int_{B_R} |f| \zeta^2 |u| \, d\mathcal{L}^n.
\end{align*}
[/step]
[step:Estimate the forcing term at the correct scale]
Apply [Young's Inequality](/theorems/244) in the pointwise form $2ab \le a^2+b^2$ with $a := \delta |f|\zeta$ and $b := \delta^{-1}|u|\zeta$.
This gives
\begin{align*}
|f|\zeta^2|u|
\le
\frac{\delta^2}{2}|f|^2\zeta^2
+
\frac{1}{2\delta^2}|u|^2\zeta^2.
\end{align*}
Since $0 \le \zeta \le 1$,
\begin{align*}
\int_{B_R} |f|\zeta^2|u| \, d\mathcal{L}^n
\le
\frac{\delta^2}{2}
\int_{B_R} |f|^2 \, d\mathcal{L}^n
+
\frac{1}{2\delta^2}
\int_{B_R} |u|^2 \, d\mathcal{L}^n.
\end{align*}
Substituting this estimate into the previous bound yields
\begin{align*}
\int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n
\le
(4C_0^2+1)\delta^{-2}
\int_{B_R} |u|^2 \, d\mathcal{L}^n
+
\delta^2
\int_{B_R} |f|^2 \, d\mathcal{L}^n.
\end{align*}
[/step]
[step:Restrict the cutoff estimate to the smaller ball]
Because $\zeta=1$ on $B(x_0,r)$ and $\zeta^2|\nabla u|^2 \ge 0$ on $B_R$,
\begin{align*}
\int_{B(x_0,r)} |\nabla u|^2 \, d\mathcal{L}^n
\le
\int_{B_R} \zeta^2|\nabla u|^2 \, d\mathcal{L}^n.
\end{align*}
Combining this with the preceding estimate and recalling that $\delta=R-r$, we get
\begin{align*} \int_{B(x_0,r)} |\nabla u|^2 \, d\mathcal{L}^n \le C\left(\frac{1}{(R-r)^2}\int_{B(x_0,R)} |u|^2 \, d\mathcal{L}^n +(R-r)^2\int_{B(x_0,R)} |f|^2 \, d\mathcal{L}^n\right), \end{align*}
where $C := \max\{4C_0^2+1,1\}$ depends only on the cutoff construction and hence only on the dimension. This is the desired Caccioppoli inequality.
[/step]