For context: there are Schwartz distributions called Pauli-Jordan (or Schwinger) functions which come up in quantum field theory. Given $m>0$, $x_0\in\mathbb{R}$ and $\vec{x}\in \mathbb{R}^3$, $D_m^{\pm}(x_0,\vec{x})\in \mathcal{S}(\mathbb{R}^4)'$ is defined by \begin{align*} D_m^\pm(x_0,\vec{x}):=f&\mapsto \frac{1}{(2\pi)^3}\int_{\mathbb{R}^4}e^{\pm i px}f(p)\delta(p^2-m^2)\theta(x_0)d^4p \\ &=\frac{1}{(2\pi)^3}\int_{\mathbb{R}^3}\frac{e^{\pm ix_0\sqrt{|\vec{p}|^2+m^2}\mp i\vec{x}\cdot\vec{p}}f\left(\sqrt{|\vec{p}|^2+m^2},\vec{p}\right)}{2\sqrt{|\vec{p}|^2+m^2}}~d^3\vec{p} \end{align*} where $\delta$ is the Dirac delta, $\theta$ is the Heavyside step function, and $p^2:=p_0^2-p_1^2-p_2^2-p_3^2$. As usual in distribution theory, the $f\mapsto$ is left implicit, and with the shorthand $\omega(\vec{p}):=\sqrt{|\vec{p}|^2+m^2}$ this becomes \begin{align} D_m^\pm(x_0,\vec{x}):=\frac{1}{(2\pi)^3}\int_{\mathbb{R}^3}\frac{e^{\pm ix_0\omega(\vec{p})\mp i\vec{x}\cdot\vec{p}}}{2\omega(\vec{p})}~d^3\vec{p} \end{align}
This integral is quite difficult, but it can be evaluated to \begin{equation} \frac{\epsilon(x_0)\delta(x^2)}{4\pi}\mp \frac{im\theta(x^2)\big(Y_1(m\sqrt{x^2})\mp i\epsilon(x_0)J_1(m\sqrt{x^2})\big)}{8\pi \theta(\sqrt{x^2}) }\mp\frac{im\theta(-x^2)\big(J_1(iz)+iY_1(iz)\big)}{8\pi \theta(\sqrt{-x^2})}\end{equation} where $\epsilon(x_0)$ is the sign function, and $J_1(z), Y_1(z)$ are the Bessel functions of the first and second kind. My question concerns representing $D^{\pm}(x)$ as the boundary-value of a holomorphic function. More specifically, let $V^+:=\{(x_0,\vec{x})\in \mathbb{R}^4\mid x_0>|\vec{x}|\}$ be the cone of "forward directed time-like vectors". It can be shown that if $z=(z_0,\vec{z})\in \mathbb{R}^4\pm iV^+$, then the integral \begin{align} F_{\pm}(z):=\frac{1}{(2\pi)^3}\int_{\mathbb{R}^3}\frac{e^{\pm iz_0\omega(\vec{p})\mp i\vec{z}\cdot\vec{p}}}{2\omega(\vec{p})}~d^3\vec{p} \end{align} obtained by replacing $x$ with $z$ converges in the ordinary sense, and the function $F_{\pm}:\mathbb{R}^4\pm iV^+\to \mathbb{C}$ is holomorphic. I am trying to explicitly evaluate the above integral in terms of Bessel functions. So far, I tried converting to spherical coordinates:
\begin{align*} F_-(z)&=\frac{1}{(2\pi)^3}\int_{\mathbb{R}^3}\frac{e^{- iz_0\omega(\vec{p})+ i\vec{z}\cdot\vec{p}}}{2\omega(\vec{p})}~d^3\vec{p} \\ &=\frac{1}{2\pi}^3\int_0^\pi\int_0^{2\pi}\int_0^\infty\frac{\rho^2\sin\phi}{2\sqrt{\rho^2+m^2}}\mathrm{exp}\big\{-iz_0 \sqrt{\rho^2+m^2}+iz_1\rho\cos\theta\sin\phi+iz_2\rho\sin\theta\sin\phi+iz_3\cos\phi\big\}~d\rho d\theta d\phi \\ &= \frac{1}{(2\pi)^3}\int_0^\pi\int_0^\infty\frac{\rho^2\sin\phi}{2\sqrt{\rho^2+m^2}}e^{-iz_0\sqrt{\rho^2+m^2}+iz_3\rho\cos\phi}\left(\int_0^{2\pi}e^{iz_1\rho\cos\theta\sin\phi+iz_2\rho\sin\theta\sin\phi}d\theta\right) d\rho d\phi \end{align*} Using the Bessel function identities $e^{iz\cos\theta}=\sum_{n\in\mathbb{Z}}i^nJ_n(z)e^{in\theta}$, $e^{izsin\theta}=\sum_{n\in\mathbb{Z}}J_n(z)e^{in\theta}$, and $n\in\mathbb{Z}\Longrightarrow J_{-n}(z)=(-1)^nJ_n(z)$, I was able to evaluate the $\theta$ integral \begin{align*} e^{ia\cos\theta}e^{ib\sin\theta}&=\left(\sum_{n\in\mathbb{Z}}i^nJ_n(a)e^{in\theta}\right)\left(\sum_{n\in\mathbb{Z}}J_n(b)e^{in\theta}\right) \\ &=\sum_{n\in\mathbb{Z}}\left(\sum_{k\in\mathbb{Z}}i^kJ_k(a)J_{n-k}(b)\right)e^{in\theta} \\ \int_0^{2\pi} e^{ia\cos\theta}e^{ib\sin\theta}d\theta &=\int_0^{2\pi}\sum_{n\in\mathbb{Z}}\left(\sum_{k\in\mathbb{Z}}i^kJ_k(a)J_{n-k}(b)\right)e^{in\theta}d\theta \\ &=2\pi \sum_{n\in\mathbb{Z}}i^nJ_n(a)J_{-n}(b) \\ &=2\pi J_0(a)J_0(b) \end{align*} Substituting $a=z_1\rho\sin\phi$ and $b=z_2\rho\sin\phi$ gives the overall integral to be \begin{equation} F_-(z)=\frac{1}{(2\pi)^2}\int_0^\pi\int_0^\infty \frac{\rho^2\sin\phi J_0(z_1\rho\sin\phi)J_0(z_2\rho\sin\phi)e^{-iz_0\sqrt{\rho^2+m^2}+iz_3\rho \cos\phi}}{2\sqrt{\rho^2+m^2}}d\rho d\phi \end{equation} I was optimistic after all the cancellation using Bessel function identities, but after that I got stuck. Using the same identities again doesn't seem to help this time because the integrand has other terms involving $\phi$. Is there some other way to go about this? I tried to look up the derivation of the expansion of $D_m^\pm(x)$ into Bessel functions, but in both references I found they change to spherical coordinates and then say "integrating out the angle coordinates, this becomes..." and proceed to write an integral involving only $\rho$ as if integrating out the angle coordinates is so easy that you can do it in your head in one step. For example, see "Introduction to the theory of quantized fields (1980) pg. 147 or Calculating Pauli–Jordan function (2020) pg. 3. Is there some easier way to integrate out the angle coordinates that I don't know about? I was able to find an identity for the $n$-dimensional Dirac delta as an integral over angle coordinates that seems relevant \begin{equation} \delta^n(x)=\frac{(n-1)!}{(-2\pi i)^n}\int_{S^{n-1}}\frac{d\xi}{(x\xi +i0)^n} \end{equation} Any help would be greatly appreciated!