[proofplan]
We cut off the source term by multiplying $f$ with $\zeta$ and extending the product by zero to obtain a compactly supported $C^2$ function $g$ on $\mathbb{R}^n$. The Newtonian potential $v=Ng$ then solves the whole-space equation $-\Delta v=g$. Since $\zeta=1$ on $B(x_0,r)$, the equation for $v$ agrees with the equation for $u$ on $B(x_0,r)$, so their difference is harmonic there.
[/proofplan]
[step:Construct the compactly supported source on $\mathbb{R}^n$]
Let $K=\operatorname{supp}\zeta$ denote the closed support of $\zeta$ inside $B(x_0,R)$. Since $\zeta \in C_c^\infty(B(x_0,R))$, the set $K$ is compact and satisfies $K \subset B(x_0,R)$. Define
\begin{align*}
g: \mathbb{R}^n \to \mathbb{R}
\end{align*}
by
\begin{align*}
g(x)=\zeta(x)f(x) \quad \text{for } x \in B(x_0,R),
\end{align*}
and
\begin{align*}
g(x)=0 \quad \text{for } x \in \mathbb{R}^n \setminus B(x_0,R).
\end{align*}
Because $\zeta$ vanishes on a neighbourhood of $\mathbb{R}^n \setminus B(x_0,R)$ and $f \in C^2(B(x_0,R))$, this zero extension satisfies $g \in C_c^2(\mathbb{R}^n)$. Moreover, for every $x \in B(x_0,r)$, the hypothesis $\zeta=1$ gives
\begin{align*}
g(x)=f(x).
\end{align*}
[/step]
[step:Apply the Newtonian potential to solve the whole-space equation]
Define
\begin{align*}
v: \mathbb{R}^n \to \mathbb{R}
\end{align*}
by
\begin{align*}
v(x)=Ng(x)=\int_{\mathbb{R}^n}\Gamma(x-y)g(y)\,d\mathcal{L}^n(y).
\end{align*}
Since $g \in C_c^2(\mathbb{R}^n)$, the standard Newtonian potential theorem for the whole-space Poisson equation applies: the potential $v=Ng$ belongs to $C^2(\mathbb{R}^n)$ and satisfies
\begin{align*}
-\Delta v = g
\end{align*}
pointwise on $\mathbb{R}^n$.
(citing a result not yet in the wiki: Newtonian potential solves the whole-space Poisson equation for compactly supported $C^2$ data)
Restricting this identity to $B(x_0,r)$ and using $g=f$ there, we obtain
\begin{align*}
-\Delta v = f
\end{align*}
pointwise on $B(x_0,r)$.
[guided]
The purpose of the cutoff is to create a source term on all of $\mathbb{R}^n$, because the Newtonian potential is a whole-space operator. We have already defined
\begin{align*}
g: \mathbb{R}^n \to \mathbb{R}
\end{align*}
as the zero extension of $\zeta f$. The support of $\zeta$ is compactly contained in $B(x_0,R)$, so the extension introduces no loss of regularity across the boundary of $B(x_0,R)$: near that boundary, $\zeta f$ is identically zero. Hence $g \in C_c^2(\mathbb{R}^n)$.
Now define
\begin{align*}
v: \mathbb{R}^n \to \mathbb{R}
\end{align*}
by the Newtonian potential formula
\begin{align*}
v(x)=Ng(x)=\int_{\mathbb{R}^n}\Gamma(x-y)g(y)\,d\mathcal{L}^n(y).
\end{align*}
The relevant whole-space result says that if the source $g$ is compactly supported and $C^2$, then its Newtonian potential is a $C^2$ solution of the Poisson equation
\begin{align*}
-\Delta v = g
\end{align*}
on $\mathbb{R}^n$.
(citing a result not yet in the wiki: Newtonian potential solves the whole-space Poisson equation for compactly supported $C^2$ data)
This is exactly why we introduced $g$: on the smaller ball $B(x_0,r)$, the cutoff satisfies $\zeta=1$, so for every $x \in B(x_0,r)$ we have
\begin{align*}
g(x)=\zeta(x)f(x)=f(x).
\end{align*}
Therefore the whole-space equation for $v$ restricts to
\begin{align*}
-\Delta v = f
\end{align*}
pointwise on $B(x_0,r)$.
[/guided]
[/step]
[step:Subtract the Newtonian potential to obtain the harmonic remainder]
Define
\begin{align*}
h: B(x_0,r) \to \mathbb{R}
\end{align*}
by
\begin{align*}
h(x)=u(x)-v(x).
\end{align*}
Since $u \in C^2(B(x_0,R))$ and $v \in C^2(\mathbb{R}^n)$, the restriction of $h$ to $B(x_0,r)$ belongs to $C^2(B(x_0,r))$. On $B(x_0,r)$, using $-\Delta u=f$ and $-\Delta v=f$, we compute
\begin{align*}
-\Delta h = -\Delta(u-v)= -\Delta u+\Delta v = f-f=0.
\end{align*}
Thus $\Delta h=0$ on $B(x_0,r)$, so $h$ is harmonic there. By the definition of $h$, for every $x \in B(x_0,r)$,
\begin{align*}
u(x)=v(x)+h(x)=Ng(x)+h(x).
\end{align*}
This is the required local [representation formula](/theorems/39).
[/step]