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.