Solution (source code)

= Solution

Define the candidate potential by convolving the decaying <Green function> with the source. Acting with the <Laplacian> in the observation variable gives, in the sense of <distributions>,
$$
\begin{aligned}
\nabla_{\mathbf r}^2\int_{\mathbb R^3}G(\mathbf r,\mathbf r')\rho(\mathbf r')\,d^3r'
&=-\int_{\mathbb R^3}\delta^{(3)}(\mathbf r-\mathbf r')\rho(\mathbf r')\,d^3r'\\
&=-\rho(\mathbf r).
\end{aligned}
$$
The final equality is the sifting property of the <Dirac delta>. The kernel found in part (ii) is locally <integrable>, and its convolution with a sufficiently regular rapidly decaying source also vanishes at spatial infinity. Thus it satisfies the <Poisson equation> and the required condition at infinity.

The difference from any other regular decaying solution is a <harmonic function> on all of space vanishing at infinity, hence is zero by the <maximum principle>. Therefore
$$
\boxed{\phi(\mathbf r)=\int_{\mathbb R^3}G(\mathbf r,\mathbf r')\rho(\mathbf r')\,d^3r'.}
$$
One can also obtain this expression from <Green's second identity> on larger and larger balls, with the singularity of the kernel treated distributionally. The surface term at infinity vanishes for the usual $1/r$ potential decay and $1/r^2$ gradient decay.