Without loss of generality, we can create spherical coordinates $(p, \theta, \phi)$ such that $\mathbf{r}$ is along the $\theta=0$ axis. This may seem a little backwards, since we usually have $\mathbf{r}$ as the position vector, which naturally varies. But here we leave $\mathbf{r}$ constant, and vary in $\mathbf{p}$-space, since that's what we're integrating over. In particular, variations in $\theta$ and $\phi$ correspond to variations in $\mathbf{p}$ rather than $\mathbf{r}$. Then, we can transform the integral into those coordinates, perform the trivial integral over $\phi$, and then change coordinates again using $\sin\theta\, d\theta = -d\cos\theta$, and finally perform that integral too, leaving only the integral over $p$. \begin{align} D(x-y) &= \int_0^\infty \int_0^\pi \int_0^{2\pi} \frac{1}{(2\pi)^3}\, \frac{1}{2E_{\mathbf{p}}} e^{i\, p\, r\, \cos\theta}\, p^2\, \sin\theta\, d\phi\, d\theta\, dp \\ &= \frac{1}{(2\pi)^2}\, \int_0^\infty \frac{p^2}{2E_{\mathbf{p}}} \left[\int_0^\pi e^{i\, p\, r\, \cos\theta}\, \sin\theta\, d\theta\right]\, dp \\ &= \frac{1}{(2\pi)^2}\, \int_0^\infty \frac{p^2}{2E_{\mathbf{p}}} \left[\int_1^{-1} -e^{i\, p\, r\, x}\, dx\right]\, dp \\ &= \frac{1}{(2\pi)^2}\, \int_0^\infty \frac{p^2}{2E_{\mathbf{p}}} \frac{-e^{-i\, p\, r} - -e^{i\, p\, r}}{i\, p\, r}\, dp. \end{align} Note that this assumes $E_{\mathbf{p}}$ is independent of $\theta$ and $\phi$, but that's fine because we usually have $E_{\mathbf{p}} = \sqrt{ \lvert \mathbf{p} \rvert^2 + m^2} = \sqrt{ p^2 + m^2}$.