 Boundary flux from a singular harmonic potential

Although $1/r$ is harmonic away from the origin in three dimensions, it is not a regular harmonic function on a ball containing the origin. Its gradient is $-\hat{r}/r^2$, so the displayed improper volume integral follows by canceling the $r^2$ in the spherical volume element. On the annulus $\epsilon <r<a$, Green's first identity gives outer boundary contribution $-4\pi$ and inner contribution $4\pi \epsilon /a$, recovering the same limit. Its global Laplacian contains a Dirac delta distribution, so the regular-harmonic gradient-pairing theorem is inapplicable.