 Solution

Write $q=|r-r^{\prime }|$. Away from the source, the radial Laplacian gives ${\nabla }^2(-1/(4\pi q))=q^{-2}{\partial }_q[q^2/(4\pi q^2)]=0$. Through a small source-centred sphere, the outward flux of its gradient isThis identifies a unit Dirac delta at the source. More explicitly, for a smooth compactly supported test function $\psi$, apply Green's second identity outside the radius-$\epsilon$ ball. The inner outward normal of that punctured domain is $-e_q$. The term $H{\partial }_n\psi$ has integral $O(\epsilon )$, while $-\psi {\partial }_{n}H$ tends to $\psi (r^{\prime })$. Thus $\int H{\nabla }^2\psi dV=\psi (r^{\prime })$, proving the distributional identityWhenever the integral converges, convolution with this Green function of the Laplacian solves the Poisson equation for Newtonian gravity:The harmonic term represents the homogeneous freedom and is fixed by boundary conditions. For an isolated finite mass distribution, regularity away from the prescribed sources and the normalization $\Phi \to 0$ at infinity select $h=0$.

For a point mass at the origin, the delta integral gives the Newtonian potential of a point mass,