How are Fourier Transforms, Green's Functions and the Poisson Equation related? So I have a question. It starts, as many good things in life do, with Fourier Transforms.
I have been playing around with these things ever since I read somewhere that:
$$
\mathcal{F} \left\{ \frac{d u}{dx} \right\}  = \omega \hat u
$$
Where the notation $\hat u$ is the Fourier Transform of $u(\vec r)$ and $\omega = i \hat x$ where $\hat x$ is the reciprocal coordinate such that $\hat u = \hat u(\hat x)$. I managed to, using the definition of Fourier Transform found in Wikipedia, convince myself that this identity is true. So I ran with it.
If that thing is true, than this following expression must also be true:
$$
\frac{d u}{dx} = \mathcal{F}^{-1} \{\omega \hat u\}
$$
Which is a pretty weird way of differentiating something when you know how to do that symbolically, but by this point I had used discrete Fourier Transforms in a lot of grid-computations and I knew this to be quite a handy way of doing calculus with discretized functions for which I had no analytical form. So I got excited like a musician. You know what musicians do when they get excited and moved by the beauty in front of them? They compose. So I did, and here is what I got:
$$
\frac{d }{dx} \left[ \frac{d u}{dx} \right] = \mathcal{F}^{-1}\{ \omega \mathcal{F}\{ \mathcal{F}^{-1} \{\omega \hat u\}\}\}
$$
Which really is just:
$$
\frac{d }{dx} \left[ \frac{d u}{dx} \right] = \mathcal{F}^{-1}\{ \omega^2 \hat u\}
$$
This started looking a lot like this:
$$
\frac{d^n u}{dx^n} = \mathcal{F}^{-1}\{ \omega^n \hat u\}
$$
Which look like a general expression for the $n$-th derivative of any function $u$. Even one that has no analytical form. I googled for about ten minutes and figured out that, indeed. That is true. And not only that but when $n$ is negative, it turns out that it stop being a derivative and starts being an antiderivative, a definite integral. I was mindblown.
Fourier, it turns out, had unified calculus under a single computation framework. It is no wonder that this operator on $u$ is called the Fourier Differintegral. Charming.
Now, what else can we do with derivatives and integrals? In 1D, it turns out, not much. But things get juicy when you go full planar, adding a whole new dimension to play with. With it, comes a coordinate, $y$, a reciprocal coordinate $\hat y$  and an imaginary reciprocal coordinate $\sigma = i \hat y$. Ah, and also a position vector $\vec x = (x,y)$, a reciprocal position vector $\hat {\vec x} = (\hat x,\hat y)$ and its imaginary counterpart ${\vec \omega} = i \hat {\vec x} = (\omega, \sigma)$
Now we arrive at things I pulled out of my ass. For partial derivatives, I'm doing it like this:
$$
\mathcal{F}\left\{ \frac{\partial u}{\partial x} \right\} = \omega \hat u
$$
Since $u$ is now a function of $\vec x$ I understand that $\hat u$ is a multidimensional Fourier Transform in 2D, and a function of $\vec \omega$.
Now, I guess that, equivalently,
$$
\mathcal{F}\left\{ \frac{\partial u}{\partial y} \right\} = \sigma \hat u
$$
Notice how the imaginary reciprocal variable that multiplies the transform is actually related to the variable with respect to which the derivative is being taken.
I don't know for sure that this is the case. I just assumed this was right and then started computing things. One thing, in particular that I calculated was this:
$$
\frac{\partial }{\partial y} \left[ \frac{\partial u}{\partial x} \right] = \mathcal{F}^{-1}\{ \sigma \mathcal{F}\{ \mathcal{F}^{-1} \{\omega \hat u\}\}\}
$$
Which is really:
$$
\frac{\partial }{\partial y} \left[ \frac{\partial u}{\partial x} \right] = \mathcal{F}^{-1}\{ \sigma \omega \hat u\}
$$
I was then very confident to say that:
$$
\frac{\partial^2 u}{\partial x^2} = \mathcal{F}^{-1}\{ \omega^2 \hat u\}
$$
And likewise for
$$
\frac{\partial^2 u}{\partial y^2} = \mathcal{F}^{-1}\{ \sigma^2 \hat u\}
$$
And here I noticed that:
$$
\nabla^2 u = \mathcal{F}^{-1} \{[\omega^2 + \sigma^2] \hat u\} = \mathcal{F}^{-1} \{ |\vec \omega|^2 \hat u\}
$$
By this point I was really confident, so I decided to go up a notch and add a third dimension with coordinate $z$, reciprocal $\hat z$ and imaginary reciprocal coordinate $\kappa$.
So in 3D I got:
$$
\nabla^2 u = \mathcal{F}^{-1} \{[\omega^2 + \sigma^2 + \kappa^2] \hat u\} = |\vec \omega|^2 \hat u
$$
And here I have grown arrogant (and potentially sloppy) and I thought I could solve the Poisson Equation just using this Fourier magic calculus I had uncovered, so the gods of mathematics punished me for my hubris. It went like this:
Here is the Poisson Equation:
$$
\nabla^2 \varphi = \rho
$$
Now if I take the Fourier Transform on both sides, here is what I get:
$$
|\vec \omega|^2 \hat \varphi = \hat \rho
$$
Which suggests that I can just...
$$
\hat \varphi = |\vec \omega|^{-2} \hat \rho
$$
And taking the inverse Fourier Transform on both sides we get that:
$$
\varphi = \mathcal{F}^{-1} \{ |\vec \omega|^{-2} \hat \rho \}
$$
But then I found out that this type of equation is actually solved using Green's Functions, like so:
$$
\varphi = -\int \int \int \frac{\rho({\vec r}')}{4\pi |\vec r - {\vec r}'|} d^3 \vec r'
$$
Which looks rather... strange. Suddenly there is a 6D function ($|\vec r - {\vec r}'|^{-1}$) multiplying $\rho$ and the entire thing looks like it was solved using volume integrals, which should be like this:
$$
\int \int \int u d^3 \vec r = \mathcal{F}^{-1}\left\{\frac{1}{\omega \sigma \kappa}\hat u \right\}
$$
And not this:
$$
\varphi = \mathcal{F}^{-1} \left\{ \frac{1}{\omega^2 + \sigma^2 + \kappa^2} \hat \rho \right\}
$$
... which is the expression I had found using the Fourier Differintegral.
So I can't help but ask:
Where did I drop the ball?
Where in my reasoning did I assume wrong things? I already know I dropped the ball somewhere, since my method did not arrive at something resembling Green's Functions, I just don't know where and I feel the answer will teach me a valuable lesson about the relationship between Fourier Transforms, Green's Functions, convolutions and Poisson's Equations.
Maybe that was the lesson the gods were trying to teach me. May wisdom make humble people of us all.
Cheers.
EDIT 1
This edit is intended as a response to @martin 's answer.
Martin, you said in your answer that:
$$
\varphi = \mathcal{F}^-1\{|\vec \omega|^{-2}\} \ast \rho
$$
So I went around trying to figure out what \mathcal{F}^-1{|\vec \omega|^{-2}} is and how to solve a convolution in 3D. Here's what I got:
So formula 502 in this table shows me that:
$$
\mathcal{F} \{ |\vec x|^{-a} \} = \frac{(2 \pi)^a}{c_{n,a}} |\vec \omega|^{-(n-a)} 
$$
if $0 < \Re a < n$
where:
$$
c_{n,a} = \pi^\frac{n}{2} 2^a \frac{\Gamma(a/2)}{\Gamma \left( \frac{n-a}{2} \right) }
$$
I noticed that the pair $a=1$ and $n=3$ satisfies the constraint that $0 < \Re a < n$ and gives formula 502 a promising look:
$$
\mathcal{F} \{ |\vec x|^{-1} \} = \frac{2 \pi}{c_{3,1}} |\vec \omega|^{-2} 
$$
also, for this particular pair,  the complicated $c_{n,a}$ formula acquires a much simpler form, check it out:
$$
c_{3,1} = \pi^\frac{3}{2} 2 \frac{\Gamma(1/2)}{\Gamma(1)}
$$
since $\Gamma(1/2) = \sqrt{\pi}$  and $\Gamma(1) = 1$ we have that:
$$
c_{3,1} = \pi^\frac{3}{2} 2 \sqrt{\pi}$
$$
And we can use the relationships between exponentiation and square roots to simplify it even further like so:
$$
\pi^\frac{3}{2} \sqrt{\pi} = \pi^\frac{3}{2} \pi^\frac{1}{2} = \pi^\frac{3+1}{2} = \pi^\frac{4}{2} = \pi^2
$$
Which means that $c_{3,1}=2\pi^2$! How neat is that?
-ahem- Anyway, we can plug that back into the formula as:
$$
\mathcal{F} \{ |\vec x|^{-1} \} = \frac{2 \pi}{2\pi^2} |\vec \omega|^{-2} 
$$
That warrants some simplification, which gives us
$$
\mathcal{F} \{ |\vec x|^{-1} \} = \frac{1}{\pi} |\vec \omega|^{-2} 
$$
Because the Fourier Transform and its inverse are integrals and $\frac{1}{\pi}$ is a constant, I assume I can do this:
$$
 \pi |\vec x|^{-1} = \mathcal{F}^{-1} \{|\vec \omega|^{-2} \}  
$$
I can then define a function $f(\vec x) = \pi |\vec x|^{-1}$ and, going back to our Poisson equation $\nabla \varphi = \rho$, we can state that
$$
\varphi = f \ast \rho
$$
Now, I didn't find any compact, understandable formula for a 3D integral, but I found this one:
$$(f \ast g)(t)=\int_{-\infty}^{\infty} f(\tau) g(t-\tau) d \tau$$
And since convolutions are commutative then:
$$(f \ast g)(t)=(g \ast f)(t)=\int_{-\infty}^{\infty} f(t-\tau)  g(\tau) d \tau$$
so now I'm praying that this thing, when blown to 3D dimensionality, looks like this:
$$
(f \ast \rho)(\vec x)=\int\int\int_{-\infty}^{\infty} f({\vec x}-{\vec x}')  g(\vec x ') d^3 \vec x '
$$
When substituting $r=|{\vec x} - {\vec x}'|$ and $\rho(\vec x ') = \rho'$ we get that
$$
\varphi = (f \ast \rho)(\vec x)=\int\int\int_{-\infty}^{\infty}\frac{\pi \rho'}{r} d^3 \vec x '
$$
Which looks different from the correct one:
$$
\varphi = (f \ast \rho)(\vec x)=-\int\int\int_{-\infty}^{\infty}\frac{\rho'}{4 \pi r} d^3 \vec x '
$$
I'm thinking I dropped the ball somewhere around the transition from 1D convolutions to 3D convolutions, but IDK. Could be an incorrect usage of formula 502. Once again I have to ask Where did I drop the ball? but at least now I know where the 6D pairwise potential comes from.
 A: It's great that you discovered those identities.
One thing that should help you make sense of the 3D part (and might also lead you to revisit your 1D and 2D exploration), and something that you might know already from your experience with discrete Fourier transform, is the relation between multiplication and convolution in Fourier analysis. Namely, given two sufficiently "nice" functions $u$ and $v$, defined on $\mathbb{R}^3$ (or $\mathbb{R}^n$, for any $n \in \mathbb{N}$...) one has the following classical relations:
$$\mathcal{F}(u\,v) = \mathcal{Fu}*\mathcal{Fv}\,, \quad \mathcal{F}(u*v) = \mathcal{Fu}\,\mathcal{Fv}$$
where $*$ is the convolution operation (depending on the constants chosen in the definition of the Fourier transform, you  might have to amend the previous identities with multiplicative constants, but let's forget about that here). The same relations are also satisfied when replacing $\mathcal{F}$ by $\mathcal{F}^{-1}$, again up to some multiplicative constants. By definition, the convolution of two functions $f,g:\mathbb{R}^n \to \mathbb{C}$, denoted by $f*g:\mathbb{R}^n \to \mathbb{C}$, is
$$\forall x \in \mathbb{R}^n\,, \quad (f*g)(x) := \int_{\mathbb{R}^n} f(x - y)g(y) dy\,.$$
Again, you need $f$ and $g$ to be quite "nice" for this to make sense (try the constant functions $f = g = 1$ to see how it could go very, very wrong).
When you correctly (up to technical details) write that $-\Delta u = \varphi$ implies
$$u = \mathcal{F}^{-1}(|\vec{\omega}|^{-2} \hat{\varphi})$$
you may then deduce
$$u = \mathcal{F}^{-1}(|\vec{\omega}|^{-2})* \mathcal{F}^{-1}(\hat{\varphi})\,.$$
There is some simplification going on in the rhs, since of course $\mathcal{F}^{-1}(\hat{\varphi}) = \varphi$, and thus you get $u = G * \varphi $ where
$$G = \mathcal{F}^{-1}(|\vec \omega|^{-2})$$
I wonder what this last inverse Fourier transform could be... But I trust the mathematical adventurer in you to answer that.
N.B.: Everything above is very "formal", i.e. it is a bit challenging to make all of this precise. But the formal calculations are always illuminating, nevertheless. For example, when I said that you need "nice" functions, you could understand "functions that are infinitely many times differentiable, with a fast decay at infinity" (Schwartz functions). But in this regard, is $|\vec \omega|^{-2}$ that nice? What about the singularity at $0$? Those are standard issues in the theory of Partial Differential Equations (of which the Poisson equation is the most prominent example).
Edit: answer to your edit
Your computations are all exact, and you almost got to the correct conclusion up to a constant factor, which is furthermore a power of $2\pi$. This clearly indicates that the issue merely comes from consistency of conventions in the Fourier transform definition. After looking carefully for the faulty constant, I finally tracked it down. If you are to use formula 502 in the wiki link, then it means you are using the convention
$$\mathcal{F}(f)(\omega) := \int_{\mathbb{R}^n} f(x) e^{-2i\pi x\cdot \omega}\,d\omega$$
(I replaced the $\xi$s by $\omega$s, which looks strange, but this is nothing else than a change of notation).
Under that convention, one has (you can easily prove this using integration by parts)
$$\mathcal{F}\left(\frac{\partial f}{\partial x_k}\right)(\omega) := {\color{red}{2i\pi}}\, \omega_k \hat{f}(\omega)$$
so that
$$\mathcal{F}\left(\Delta f\right)(\omega) := -{\color{red}{4\pi^2}} |\omega|^2 \hat{f}(\omega)\,.$$
There you have the missing factor.
P.S.: I now realize you have called $\hat{x}$ your Fourier variable and defined $\omega = i \hat{x}$ at the beginning of your answer, so my answer could confuse you a bit. If that's the case, check precisely what is the definition of the Fourier transform that you want to use. If its (in 1D)
$$ \mathcal{F}(f)(\hat{x}) = \int_{\mathbb{R}} f(x) e^{-2i \pi x \hat{x}}d\hat{x}$$
then the differentiation theorem reads
$$ \mathcal{F}(\frac{df}{dx})(\hat{x}) = 2i\pi\hat{x} \hat{f}(\hat{x}) = 2\pi \omega \hat{f}(\hat{x})$$
(you miss the $2\pi$ there)
Later on, you write $|\vec \omega|^2 = \omega^2 + \sigma^2 + \kappa^2$. But written as a function of $\hat{x},\hat{y},\hat{z}$, this is really
$|\vec \omega|^2 = -\hat{x}^2 - \hat{y}^2 - \hat{z}^2$. In fact, your `` vector" $\vec \omega$ has imaginary coordinates, and $|\vec \omega|^2$ is negative... This is very misleading. Maybe try to rewrite everything just with $\hat{x}, \hat{y}, \hat{z}$. You could call $\hat{X} = (\hat{x},\hat{y},\hat{z})$, and then the Fourier inversion you're looking for should be
$$\hat{G}(\hat{X}) = \mathcal{F}^{-1}\left(\frac{-1}{4\pi^2|\hat{X}|^2}\right)$$
