Neumann point source on a half-space boundary ($\Phi (x,y,z)=1/[2\pi \sqrt{x^2+y^2+z^2}]$)

In the upper half-space $z>0$, the displayed decaying harmonic function has outward normal derivative ${\partial }_n\Phi =-{\partial }_z\Phi =z/[2\pi (x^2+y^2+z^2{)}^{3/2}]$. This is the Poisson kernel for the upper half-space, whose integral on a horizontal plane is one and whose limit as $z\downarrow 0$ is the two-dimensional Dirac delta. Thus ${\partial }_n\Phi =\delta (x)\delta (y)$ in the distributional boundary sense. The compensating flux through a large hemisphere is minus one. This illustrates why a nonzero total boundary source is possible on an unbounded domain even when the interior source is zero.