Greens function of Laplace operator

I am having trouble deriving the Greens function for the three dimensional Laplacian $\nabla^2$. I wish to solve the equation

$$ -\nabla_r^2 G(r,r') = \delta(r-r') \quad\quad\quad r,r'\in\mathbb{R}^3$$

I know that $\frac{1}{4\pi|r-r'|}$ is the solution, but I would like to derive this for myself. I believe that the correct approach is to apply a Fourier transform to go from $r$ space to $k$ space:

$$ -\int dr \; e^{-ikr}\nabla_r^2 G(r,r') = \int dr\; e^{-ikr}\delta(r-r') $$

$$ k^2 \tilde{G}(k,r') = e^{-ikr'} $$

Where $\tilde{G}(k,r') \equiv \int dr \; e^{-ikr} G(r,r')$ is the Fourier transform of $G$. We can solve for $\tilde{G}$ easily:

$$ \tilde{G}(k,r') = \frac{e^{-ikr'}}{k^2} $$

And now we transform back into $r$ space:

$$ G(r,r') = \int dk\; e^{+ikr}\tilde{G}(k,r') = \int dk\;\frac{e^{ik(r-r')}}{k^2}$$

But unfortunately this integral does not converge! We can see this by integrating by parts:

$$ G(r,r') = -\int dk\;e^{ik(r-r')}\frac{d}{dk}\frac{1}{k} $$

$$ G(r,r') = -\frac{e^{ik(r-r')}}{k} + i(r-r')\int dk\;\frac{e^{ik(r-r')}}{k}$$

$$ G(r,r') = -\frac{e^{ik(r-r')}}{k} - (r-r')^2 Ei\left(ik(r-r')\right)$$

Neither the boundary term nor the error function converge when we evaluate for $k\in\mathbb{R}^3$, and it is not at all obvious that they should cancel in any meaningful way. I don't see how to arrive at the correct solution $G(r,r')=\frac{1}{4\pi|r-r'|}$

Let us define the Green's function by the equation,

$$\nabla^2 G(\mathbf{r},\mathbf{r}_0) = \delta(\mathbf{r}-\mathbf{r}_0).$$

Now let us define $S_\epsilon = \{\mathbf{r} : |\mathbf{r}-\mathbf{r}_0| \leq \epsilon \}$, from which we thus have that,

$$\int_{S_\epsilon} \nabla^2 G(\mathbf{r},\mathbf{r}_0) \, d\mathbf{r} = \int_{S_\epsilon} \delta(\mathbf{r}-\mathbf{r}_0) \, d\mathbf{r} = 1$$

by properties of the delta function. By Green's theorem (or analogously, Stokes') we have that,

$$\int_{\partial S_\epsilon} \nabla G(\mathbf{r},\mathbf{r}_0) \cdot \mathbf{n} ds = 1.$$ In two dimensions, $S_\epsilon$ is a circle, and we find,

$$\lim_{\epsilon \to 0} \int_0^{2\pi} \frac{\partial G}{\partial r} r \bigg\rvert_{r=\epsilon} d\theta = 1$$

from which one can deduce the leading order behaviour of the Green's function in two dimensions,

$$G \sim \frac{1}{2\pi} \ln r = \frac{1}{2\pi} \ln |\mathbf{r}-\mathbf{r}_0|.$$

Repeating this analysis for the case of a sphere instead, in three dimensions, you will find that,

$$G \sim -\frac{1}{4\pi r} = -\frac{1}{4\pi |\mathbf{r}-\mathbf{r}_0|}.$$

Note these are the behaviours of the Green's function near $\mathbf{r}_0$. However, if there are no boundary conditions and we consider all of $\mathbb R^3$, then this corresponds to the free space Green's function, as we only require it to have the correct behaviour near the source.

If you choose the Laplace operator as $-\nabla^2$ so that it is elliptic, then this accounts for the difference in sign, compared to your cited result.

I prefer to start with the Green's function for a massive field (defined by: $\left[-\nabla^2 + m^2\right] G = \delta(\mathbf{x}-\mathbf{x}')$) and take the limit as $m\rightarrow 0$. That Green's function is: $$G(r=|\mathbf{x}-\mathbf{x}'|) = \frac{1}{(2\pi)^{d/2}} \left(\frac{m}{r}\right)^{d/2-1} K_{d/2-1}(mr),$$ where $K_\nu(x)$ is the modified Bessel function of the second kind, $d$ is the number of dimensions of the space. Applying the limiting forms of the modified Bessel functions to take the limit gives: $$\begin{align} \lim_{m\rightarrow 0} G(r) &= \frac{1}{(2\pi)^{d/2}} \lim_{m\rightarrow 0}\left(\frac{m}{r}\right)^{d/2-1} \frac{\Gamma(d/2-1)}{2} \left(\frac{2}{mr}\right)^{d/2-1}\\ & = \frac{\Gamma\left(\frac{d}{2} - 1\right)}{4 \pi^{d/2} r^{d-2} } \end{align}$$ when $d > 2$. For $d=2$ we get: $$\lim_{m\rightarrow0} G(r) = \frac{1}{2\pi}\lim_{m\rightarrow 0} -\ln(mr).$$ The limit in this case diverges because the massive case requires that $\lim_{r\rightarrow \infty} G(r) = 0$, and that's not required for the massless case. In the massless case we can add a constant with respect to $r$, though, so the answer is: $$G_2(r) = -\frac{\ln r}{2\pi}.$$

The 1-d case has a similar divergence as the 2-d case, but the pattern fits with: $$G_1(x) = -\frac{|x|}{2}.$$

The pattern is: $$-\frac{\partial G(r)}{\partial r} = \frac{1}{S_{d-1}(r)},$$ where $S_n(r)$ is the surface area of the $n$-sphere. The $1$-sphere is the circle, and the $0$ sphere is two isolated points (the measure of isolated points is a count of them).