= Solution
A <fundamental solution> of a constant-coefficient operator $P(D)$ is a <distribution> $E$ satisfying $P(D)E=\delta_0$. The proposed $N$ is <locally integrable> in three dimensions and harmonic away from zero. For a smooth <compactly supported> <test function> $\varphi$, apply Green's identity outside the ball of radius $\eta$. On the inner boundary the outward normal is $-\hat r$, so $\partial_\nu N=1/(4\pi\eta^2)$. The outer boundary contributes nothing. Thus
$$
\int_{|x|>\eta}N(-\Delta\varphi)\,dx=\int_{|x|=\eta}\left(\frac{\varphi}{4\pi\eta^2}+\frac{\partial_r\varphi}{4\pi\eta}\right)dS.
$$
The first term tends to $\varphi(0)$ and the second is $O(\eta)$. <Local integrability> permits taking the limit on the left. Therefore \b[$-\Delta N=\delta_0$ in <distributions>].
For a <harmonic function> on a ball, define its <spherical mean> $M(r)=(4\pi)^{-1}\int_{S^2}u(x_0+r\omega)\,d\omega$. Differentiation and the <divergence theorem> give $M'(r)=(4\pi r^2)^{-1}\int_{B_r(x_0)}\Delta u=0$. Its limit at zero is $u(x_0)$, so
$$
\boxed{u(x_0)=\frac1{4\pi r^2}\int_{\partial B_r(x_0)}u\,dS=\frac1{|B_r|}\int_{B_r(x_0)}u\,dx.}
$$
The ball formula follows by integrating <spherical means> in the radius. If two solutions of the prescribed <Poisson equation> tend to zero at infinity, their difference is globally harmonic and also tends to zero. Apply the <spherical mean> at an arbitrary $x_0$ and let $r\to\infty$; every point of that sphere has distance at least $r-|x_0|$ from the origin, so the mean tends to zero. Hence the difference is zero everywhere, proving \b[uniqueness].
Back to article page
