I was revisiting an old paper by Leslie L Foldy where he states that the solution to $$\left(\nabla^2+k^2\right)\!\psi\!\left(\vec{r}\right)=-4\pi\,\delta\!\left(\vec{r}\right)$$ is $$\psi\!\left(\vec{r}\right)=\frac{e^{ikr}}{r}.$$ Here, $\delta\!\left(\vec{r}\right)$ is the (3D) Dirac delta and $r=|\vec{r}|$.
I can't find this identity anywhere. I've looked at a bunch of books on Green's functions without luck. Mathematical methods of physics by Mathews and Walker, for example, solves the very similar (in electrodynamics) $$\left(\nabla^2-\frac{1}{c}\,\partial_{tt}\right)\!\psi\!\left(\vec{r},t\right)=-4\pi\delta\!\left(\vec{r}\right)\delta(t)$$ ($\partial_{tt}$ is the second partial derivative with respect to $t$) and obtains $$\psi(\vec{r},t)=c\,\delta(r-ct),$$ which looks nothing like Foldy's outgoing spherical wave (an imaginary exponential divided by $r$).
I tried doing it myself but got stuck. I tried two methods:
If $\Psi\!\left(\vec{k}\right)$ is the (3D) Fourier transform of $\psi\!\left(\vec{r}\right)$, then we have $$\Psi\!\left(\vec{k}\right)=\frac{1}{(2\pi)^{3/2}}\int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\infty\psi\!\left(\vec{r}\right)e^{-i\vec{k}\cdot\vec{r}}\,\text{d}\vec{r},\\ \psi\!\left(\vec{r}\right)=\frac{1}{(2\pi)^{3/2}}\int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\infty\Psi\!\left(\vec{k}\right)e^{i\vec{k}\cdot\vec{r}}\,\text{d}\vec{k}.$$ Taking the Fourier transform of the equation we wish to solve, $$\left(-k_x\,\!^2-k_y\,\!^2-k_z\,\!^2+k^2\right)\!\Psi\!\left(\vec{k}\right)=-\frac{4\pi}{(2\pi)^{3/2}},$$ whereby $$\Psi\!\left(\vec{k}\right)=-\frac{4\pi}{(2\pi)^{3/2}\left(-k_x\,\!^2-k_y\,\!^2-k_z\,\!^2+k^2\right)}$$ and $$\psi\!\left(\vec{r}\right)=-\frac{1}{2\pi^2}\int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\infty\frac{e^{i\vec{k}\cdot\vec{r}}}{-k_x\,\!^2-k_y\,\!^2-k_z\,\!^2+k^2}\,\text{d}\vec{k}.$$ This is where I get stuck; I have no idea how to solve that.
The equation we wish to solve is an inhomogeneous differential equation. The most general homogeneous solution without poles is $$\psi^\text{h}\,\!_{\ell,m}\!\left((\vec{r}\right)=j_\ell(kr)Y_{\ell,m}(\theta,\phi),$$ where $j_\ell$ is the $\ell$th spherical Bessel function of the first kind and $Y_{\ell,m}$ is the spherical harmonic of order $(\ell,m)$; however, something closer to what we're looking for is $$\psi^\text{h}\!\left(\vec{r}\right)=\frac{e^{ikr}}{r},$$ which is the general solution when $\psi^\text{h}\!\left(\vec{r}\right)=\psi^\text{h}(r)$. But then I have no idea what sort of particular solution $\psi^\text{p}$ to propose in order to arrive at $$\psi\!\left(\vec{r}\right)=\psi^\text{h}\,\!\!\left(\vec{r}\right)+\psi^\text{p}\,\!\!\left(\vec{r}\right)=\frac{e^{ikr}}{r}.$$
There's an article$^1$ which follows the second procedure I tried but then uses generalised functions to find $\psi^\text{p}$, and I got completely lost trying to follow it (I'm not familiar with generalised functions; that article was the first time I ever heard of them), so that's no use to me. The authors do reach the desired $$\psi\!\left(\vec{r}\right)=\frac{e^{ikr}}{r},$$ though, which makes me think Foldy wasn't wrong, just terrible at being clear in his writing.
How can I obtain the desired result without using magic— er, I mean, generalised functions?
$^1$ Schmalz J E, Schmalz G, Gureyev T E & Pavlov K M (2010): On the derivation of the Green's function for the Helmholtz equation using generalized functions, American Journal of Physics 78, 181–186