Peskin and Schroeder, Eq 2.51 with a spacelike interval. From the book "An Introduction to Quantum Field Theory" by Peskin and Schroeder, in page 27, they give the following derivation:
We have the amplitude:
$$
D(x - y) = \frac{1}{(2\pi)^3} \int \text{d}^{3}p \frac{1}{2 E_{p}} \exp \left(-i p \cdot (x- y) \right)
$$
And if we consider the case where $x - y$ is spacelike: $x^{0} - y^{0} = 0$, $\vec{x} - \vec{y} = \vec{r}$, then the amplitude becomes:
$$
D(x - y) = \frac{1}{(2\pi)^3} \int \text{d}^{3}p \frac{1}{2 E_{p}} \exp \left(i \vec{p} \cdot \vec{r} \right) \tag{1}
$$
Changing to spherical coordinate, we get:
$$
= \frac{2 \pi}{(2\pi)^3} \int_{0}^{\infty} \text{d}p \frac{p^{2}}{2 E_{p}} \frac{\exp\left(ipr\right) - \exp\left(-ipr\right)}{i p r} \tag{2}
$$
This part is where I get lost. I don't know how they got from the exponential function, $\exp \left(i \vec{p} \cdot \vec{r} \right)$, to something of the form:
$$\frac{\sin(pr)}{pr}$$
My initial guess was to used identities with Bessel functions. Knowing that $\vec{p} \cdot \vec{r} = p r \cos\left(\alpha\right)$, for some $\alpha$, I tried with the identities mentioned in NIST DLMF and wikipedia, but to be honest, I don't know how to continue.
If someone could give me a clue, I will appreciate. Thank you!.
 A: The solution is very simple, and I feel dumb for not catch it before. Thank you @Sal.
We consider a changing of coordinates that make the vector $\vec{r}$ rest over the z-axis, with this the dot product becomes $\vec{p} \cdot \vec{r} = p r \cos(p_{\theta})$ where $p_{\theta}$ is the zenith angle of $\vec{p}$. Plugging this in (1), we get:
$$
D(x - y) = \frac{2 \pi}{(2\pi)^3} \int_{0}^{\infty} \text{d}p \frac{p^{2}}{2 E_{p}} \left[ \int_{0}^{\pi} \text{d}p_{\theta} \sin \left(p_{\theta}\right) \exp \left(i p r \cos \left(p_{\theta}\right) \right) \tag{3} \right]
$$
From here the integral is straightforward.
The reason why this works is because $\vec{r}$ is fixed, the integral (1) only depends on the magnitude of $\vec{p}$ and $\vec{r}$ and the relative position between both vectors, then the expression is invariant under a rotation, that give us the freedom to set $\vec{r}$ in wherever place we want. One solution would be to make the dependency of the projection angle only of one of the spherical angle (in the original expression, $\alpha$ depends on both) because the dependency of the azimuthal angle is more hard to work, then we chose the polar angle, which give us the integral from before.
A: Peskin and Schroeder give an argument that this is essentially:
\begin{equation}\tag{1}
D(x - y) = \frac{1}{4\pi r^{2}}\int^{\infty}_{m}\sqrt{E^{2}-m^{2}}\exp(-\mathrm{i}Et)\,\mathrm{d}E.
\end{equation}
If we just blindly compute this integral (or plug it into Mathematica), we get a modified Bessel function of the second kind:
\begin{equation}\tag{2}
D(x - y) = \frac{-\mathrm{i}}{4\pi^{2}t}K_{1}(\mathrm{i}mt).
\end{equation}
This is a bit handwavy, since we need $m\gt0$ and the imaginary part of time be negative --- i.e., $\mathrm{Im}(t)\lt0$. But let us soldier on, and see what happens.
Using the asymptotic expansion of $K_{\alpha}(z)$ as found in Abramowitz and Stegun 9.7.2:
\begin{equation}\tag{3}
K_{\alpha}(z) \sim \sqrt{\frac{\pi}{2z}} e^{-z} \left(1 + \frac{4 \alpha^2 - 1}{8z}+ \cdots \right).
\end{equation}
Taking the leading term, then plugging it into Eq (2), we get
\begin{equation}\tag{4}
D(x - y) \sim \frac{-\mathrm{i}}{4\pi^{2}t}
\sqrt{\frac{\pi}{2\mathrm{i}mt}} e^{-\mathrm{i}mt} \left(1 + \frac{3}{8\mathrm{i}mt}+ \cdots \right).
\end{equation}
This requires a bit of cavalier calculation, though.
