Klein-Gordon field commutator integral? Consider a Klein-Gordon field $\phi$, which satisfies $$(\Box+ \omega_0^2)\phi=0$$ on points $x \equiv \{x_0,\vec{x}\},y\equiv \{y_0,\vec{y} \}$ of 4D Minkowski-spacetime. The field commutator is
$$
[\phi(x),\phi(y)]=c \int_{\mathbb{R}^3} \frac{1}{(2\pi)^3 2 \sqrt{\lvert \vec{k} \rvert^2 + \omega_0^2}}\left ( e^{-i\sqrt{\lvert \vec{k} \rvert^2 + \omega_0^2} \ (x_0-y_0)} - e^{i\sqrt{\lvert \vec{k} \rvert^2 + \omega_0^2} \ (x_0-y_0)} \right )e^{i \vec{k}\cdot (\vec{x}-\vec{y})} \ dk.
$$
I want to verify that this integral is equal to $$
[\phi(x),\phi(y)]=c \  \text{sgn}(x_0-y_0) \left ( i \omega_0 \theta (\tau^2) \frac{J_1(\omega_0 \tau)}{4\pi \tau} - \frac i {2\pi} \delta (\tau^2) \right )
$$
where $\tau\equiv\sqrt{(x_0-y_0)^2 - \lvert \vec{x}-\vec{y} \rvert^2}$, $\theta$ is the Heavyside function, and $J_1$ is the Bessel function. This identity looks intractable to me.
Writing the integral in polar coordinates and doing the integral over the angles simplifies the expression somewhat but I still cannot derive the result using either theorems regarding the Fourier transform or substitution techniques.
 A: ($c\equiv 1$) The commutator can be written as the sum of the positive- and negative-frequency part of the Pauli-Jordan function:
$$
i[\phi(x),\phi(y)]=D^+(x-y)+D^-(x-y),
$$
where
\begin{align}
D^\pm(x) &= \pm\frac{1}{i(2\pi)^3}\int\frac{d^3k}{2\sqrt{\mathbf k^2+m^2}}
e^{\pm i\sqrt{\mathbf k^2+m^2}\,x^0-i\mathbf k\cdot\mathbf x}.
\end{align}
Performing the $3$-momentum integration and letting $r\equiv|\mathbf x|$, one has
$$
D^+(x) = -\frac{1}{8\pi^3r}\int_{-\infty}^{+\infty}\frac{e^{i\sqrt{\mathbf \rho^2+m^2}x^0+i\rho r}}{\sqrt{\rho^2+m^2}}\rho d\rho=\frac{1}{4\pi r}\frac{\partial}{\partial r}\,f(r),
$$
where
$$
f(r) = \frac{i}{2\pi}\int_{-\infty}^{+\infty}\frac{e^{i\sqrt{\mathbf \rho^2+m^2}x^0+i\rho r}}{\sqrt{\rho^2+m^2}}d\rho
=\frac{i}{2\pi}\int_{-\infty}^{+\infty}e^{im(x^0\cosh\varphi+r\sinh\varphi)}d\varphi,
$$
using $\rho = m \sinh\varphi$. Now: let $\lambda\equiv(x^0)^2-r^2$,
1) if $x^0,\lambda>0$, let 
$x^0=\sqrt{\lambda}\cosh\varphi_0$ and
$r=\sqrt{\lambda}\sinh\varphi_0$, so that
\begin{align}
-f(r)=\frac{1}{i2\pi}\int_{-\infty}^{+\infty}e^{im\sqrt{\lambda}\cosh(\varphi+\varphi_0)}d\varphi=\frac{1}{2}H_0^1(m\sqrt{\lambda})=\frac{1}{2}J_0(m\sqrt{\lambda})+\frac{i}{2}Y_0(m\sqrt{\lambda})
\end{align}
(see here for reference on integral representations of Bessel functions);
2) if $x^0>0,\lambda<0$, let 
$x^0=\sqrt{-\lambda}\sinh\varphi_0$ and
$r=\sqrt{-\lambda}\cosh\varphi_0$, so that
\begin{align}
f(r)=\frac{i}{2\pi}\int_{-\infty}^{+\infty}e^{im\sqrt{-\lambda}\sinh(\varphi+\varphi_0)}d\varphi=\frac{i}{\pi}K_0(m\sqrt{-\lambda});
\end{align}
3) if $x^0<0,\lambda>0$, let 
$x^0=-\sqrt{\lambda}\cosh\varphi_0$ and
$r=\sqrt{\lambda}\sinh\varphi_0$, so that
\begin{align}
f(r)=\frac{i}{2\pi}\int_{-\infty}^{+\infty}e^{-im\sqrt{\lambda}\cosh(\varphi-\varphi_0)}d\varphi=\frac{1}{2}H_0^2(m\sqrt{\lambda})=\frac{1}{2}J_0(m\sqrt{\lambda})-\frac{i}{2}Y_0(m\sqrt{\lambda});
\end{align}
4) if $x^0,\lambda<0$, let 
$x^0=-\sqrt{-\lambda}\sinh\varphi_0$ and
$r=\sqrt{-\lambda}\cosh\varphi_0$, so that
\begin{align}
f(r)=\frac{i}{2\pi}\int_{-\infty}^{+\infty}e^{-im\sqrt{\lambda}\sinh(\varphi-\varphi_0)}d\varphi=\frac{i}{\pi}K_0(m\sqrt{-\lambda}).
\end{align}
From all this: let $\text{sign}(x^0)\equiv \varepsilon(x^0)$
$$
f(r) = \begin{cases}
\frac{1}{2i}Y_0(m\sqrt{\lambda})-\frac{\varepsilon(x^0)}{2}J_0(m\sqrt{\lambda})\text{ for }\lambda>0\\
\frac{i}{\pi}K_0(m\sqrt{-\lambda})\text{ for }\lambda<0.
\end{cases}
$$
Now, to take the derivative with respect to $r$, use
$$
\frac{1}{4\pi r}\frac{\partial}{\partial r}=-\frac{1}{2\pi}\frac{\partial}{\partial \lambda}
$$
together with the fact that the function $f$ is discontinuous at $\lambda=0$ due to of $J_0(0^+)=1$ showing up for for positive $\lambda$, then:
$$
D^+(x)=\frac{1}{4\pi}\varepsilon(x^0)\delta(\lambda)-\frac{im}{8\pi\sqrt{\lambda}}\theta(\lambda)\left[Y_1(m\sqrt{\lambda})-i\varepsilon(x^0) J_1(m\sqrt\lambda) \right]-\frac{im}{4\pi^2\sqrt{-\lambda}}\theta(-\lambda)K_1(m\sqrt{-\lambda}).
$$
Working in a similar fashion,
$$
D^-(x)=\frac{1}{4\pi}\varepsilon(x^0)\delta(\lambda)-\frac{im}{8\pi\sqrt{\lambda}}\theta(\lambda)\left[-Y_1(m\sqrt{\lambda})-i\varepsilon(x^0) J_1(m\sqrt\lambda) \right]+\frac{im}{4\pi^2\sqrt{-\lambda}}\theta(-\lambda)K_1(m\sqrt{-\lambda}).
$$
Finally
$$\boxed{
i[\phi(x), \phi(y)]=\frac{1}{2\pi}\varepsilon(x^0)\delta(\lambda)-\frac{m}{8\pi\sqrt{\lambda}}\theta(\lambda)\varepsilon(x^0) J_1(m\sqrt\lambda).}
$$
(Ref: This closely follows the derivation in Bogoliubov's book.)
