Preliminaries: In their QFT text, Peskin and Schroeder give the KG propagator (eq. 2.50)
$$ D(x-y)\equiv\langle0|\phi(x)\phi(y)|0\rangle = \int\frac{d^3p}{(2\pi)^3}\frac{1}{2\omega_\vec{p}}e^{-ip\cdot(x-y)}, $$
where $\omega_\vec{p}\equiv\sqrt{|\vec{p}|^2+m^2}$. For lightlike separations, we can choose a frame where $x-y$ is purely in the time-direction and the propagator can be put into the form (2.51)
$$ D(x-y)=\frac{1}{4\pi^2}\int^\infty_m d\omega\sqrt{\omega^2-m^2}e^{-i\omega (y^0-x^0)} \tag{1}\label{timelike_prop}, $$
where I use the $\text{diag }\eta=(-,+,+,+)$ convention.
Now, one has the following integral representation of the modified Bessel function (http://dlmf.nist.gov/10.32.8)
\begin{align} K_1(z) &= z\int^\infty_1 dt \sqrt{t^2-1} e^{-zt} \\ &= \frac{z}{m^2} \int^\infty_m dt \sqrt{t^2-m^2} e^{-zt/m}, \tag{2}\label{int_rep} \end{align}
where we go to the second line by rescaling the integration variable $t \to t/m$. Comparing \eqref{timelike_prop} with \eqref{int_rep} suggests
$$ D(x-y)=\frac{m}{(2\pi)^2|y-x|}K_1(m|y-x|), $$
where we have written the time separation in terms of the Lorentz invariant $i (y^0-x^0)=|y-x|$. (Note: there is an issue in what I've written here in that the integral representation \eqref{int_rep} is only valid for $|arg z|<\pi/2$ and $|y-x|$ is on the imaginary axis ($|arg z|=\pi$), but I think one could infinitesimally displace $z$ off of the imaginary axis to get a convergent integral. Check me on that.)
Anyway, for spacelike separations, we can choose a frame where $y-x=\vec{y}-\vec{x}\equiv\vec{r}$. Performing the polar integrations yield
$$ D(x-y)=\frac{-i}{2(2\pi)^2 r}\int^\infty_{-\infty}dp\frac{p e^{ipr}}{\sqrt{p^2+m^2}}. $$
Finally, PS claim that taking a contour integral in the upper half plane (making sure to avoid the branch cut at +im) will give
$$ D(x-y)= \frac{1}{(2\pi)^2r}\int^\infty_m d\rho \frac{\rho e^{-\rho r}}{\sqrt{\rho^2-m^2}}, \tag{3}\label{spacelike_prop} $$ where $\rho\equiv-ip$.
Question: I know from plugging into Mathematica that the spacelike propagator \eqref{spacelike_prop} can also be expressed as a modified Bessel function $K_1$. Moreover, the integration bounds of \eqref{spacelike_prop} and \eqref{int_rep} are even the same. However, I don't see how to transform the spacelike propagator integral \eqref{spacelike_prop} into the form of \eqref{int_rep}. Any ideas?
(I'd prefer, if at all possible, to use the integral representation that I've quoted \eqref{int_rep} and used for the timelike case rather than some other representation of the modified Bessel function.)