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.
($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.)