The definition of the green's function for the Klein-Gordon equation reads: $$ (\partial_t^2-\nabla^2+m^2)G(\vec{x},t)=-\delta(t)\delta(\vec{x}) $$ According to these resources:
- Green's function for the inhomogenous Klein-Gordon equation , the green's function looks like this:
$$ G(\vec{x},t) = \frac {\theta(t)} {2 \pi} \delta\Big( t^2 - |\vec{x}|^2 \Big) - \frac {m} {2 \pi} \theta(t - |\vec{x}|) \frac {J_{1}\left(m \sqrt{t^2 - |\mathbf{x}|^2}\right)} {m\sqrt{t^2 - |\mathbf{x}|^2}} $$
- From the wikipedia, there are several kinds of propagator, which depends on the choice of contour. They looks very like the above one, but still not identical to me. For example, the advanced and retarded green's function looks very similar to the above expression.
I want to derive the explicit form of the above listed green's function, here is what I tried:
Assume: $$ G(\vec{x},t)=\int~\mathrm{d}^3p~\mathrm{d}\omega \,e^{i(\vec{p}\cdot\vec{x}-\omega t)}G(\vec{p},\omega) $$
Substitute to Klein-Gordon equation one gets:
$$ G(\vec{p},\omega)=\frac{1}{(2\pi)^4}\frac{1}{\omega^2-p^2-m^2} $$
Then: $$ G(\vec{x},t)=\frac{1}{(2\pi)^4}\int~\mathrm{d}^3p~\mathrm{d}\omega \,e^{i(\vec{p}\cdot\vec{x}-\omega t)}\frac{1}{\omega^2-p^2-m^2} $$
In the above formula, the pole are on real axis, to get a finite answer, one need manipulate the pole slightly away from the real axis. In Feynman's choice, one has the left pole slightly above and right pole slightly down, like this:
Then we have: \begin{align} G(\vec{x},t)&=\frac{1}{(2\pi)^4}\int~\mathrm{d}^3p~\mathrm{d}\omega \,e^{i(\vec{p}\cdot\vec{x}-\omega t)}\frac{1}{\omega^2-p^2-m^2+i\epsilon} \\ &\text{integrate over $\omega$} \\ &=\theta(t)\frac{-i}{(2\pi)^3}\int~\mathrm{d}^3\vec{p}\frac{e^{i\vec{p}\cdot{x}}e^{-i\sqrt{p^2+m^2}t}}{2\sqrt{p^2+m^2}}+\theta(-t)\frac{-i}{(2\pi)^3}\int~\mathrm{d}^3\vec{p}\frac{e^{i\vec{p}\cdot{x}}e^{i\sqrt{p^2+m^2}t}}{2\sqrt{p^2+m^2}} \\ &\text{integrate over $\phi$ and $\theta$} \\ &=\frac{-i}{(2\pi)^2}\frac{\theta(t)}{|\vec{x}|}\int_0^\infty\frac{p\sin(p|\vec{x}|)e^{-i\sqrt{p^2+m^2}t}}{\sqrt{p^2+m^2}} + \frac{-i}{(2\pi)^2}\frac{\theta(-t)}{|\vec{x}|}\int_0^\infty\frac{p\sin(p|\vec{x}|)e^{i\sqrt{p^2+m^2}t}}{\sqrt{p^2+m^2}} \end{align}
The above derivation seems has no flaw and I don't know how to proceed the $p$ integral, and I can't see the resemblance of the current formula to the closed formula given in resources 1 and 2.
To integrate $p$, I found this integral from the book might help: $$ \int_0^\infty e^{-\beta\sqrt{\gamma^2+x^2}} \cos bx =\frac{\beta\gamma}{\sqrt{\gamma^2+\beta^2}}K_1\left(\gamma\sqrt{\beta^2+b^2}\right), \text{ with: $\mathrm{Re}\beta>0, \mathrm{Re}\gamma>0$} $$
Edit2
---to address the updated answer by @Solenodon Paradoxus
The answer suggested to rotate the contour counter-clockwise ($\omega=i\omega'$), therefore: \begin{align} G_F(\vec{p},\omega)&=\frac{1}{-\omega'^2-|\vec{p}|^2-m^2+i\epsilon} \\ &=-\int_0^\infty \mathrm{d} L\,e^{-L(\omega'^2+|\vec{p}|^2+m^2-i\epsilon)} \end{align}
Plug the above $G_F(\vec{p},\omega)$ into the integral of four momentum: \begin{align} G_F(\vec{x},t)&=\frac{-i}{(2\pi)^4}\int_0^\infty~\mathrm{d}L\,e^{-(m^2-i\epsilon)L}\int\mathrm{d}\omega'\mathrm{d}^3\vec{p}\,e^{-L\omega'^2+t\omega'}e^{-Lp^2+i\vec{x}\cdot\vec{p}} \\ &=\frac{-i}{16\pi^2}\int_0^\infty~\mathrm{d}L\,e^{-(m^2-i\epsilon)L}\frac{1}{L^2}e^{\frac{\tau^2}{4L}}\quad\text{ with $\tau^2=t^2-|\vec{x}|^2$} \end{align}
Since now $\epsilon$ is unimportant to conduct the above integral, just omit it. \begin{align} G_F(\vec{x},t)&=\frac{-i}{16\pi^2}\int_0^\infty\mathrm{d}L\,e^{-m^2L}\frac{1}{L^2}e^{\frac{\tau^2}{4L}} \\ &=\frac{-im^2}{16\pi^2}\int_0^\infty\mathrm{d}\xi\,e^{-\xi}e^{\frac{m^2\tau^2}{4\xi}}\frac{1}{\xi^2} \end{align} When $\tau^2>0$, the integrand diverges at $\xi\to0$, the integral can't be conducted.
When $\tau^2<0$, we have: $$ -\frac{im}{4\pi^2\sqrt{-\tau^2}}K_1\left(m\sqrt{-\tau^2}\right) $$ where $K_1$ is modified Bessel function, this is exactly form of the Feynmann propagator when $\tau^2<0$.
Question: Although our result triumphs on one side, what about the other side( $\tau^2>0$)? How can we get it from the above procedure? In which step did we exclude this possibility?