Evaluating $\int_0^{2 \pi} \int_0^\infty \int_0^\infty e^{-u} e^{-v} e^{i \alpha (\sin(\omega u + \phi) - \sin(\omega v + \phi))} du\ dv\ d\phi $ I recently came across the following definite integral while working on a model involving the overlap of quantum-mechanical photon wavefunctions under some periodic modulation:
$$\int_0^{2 \pi} \int_0^\infty \int_0^\infty \mathrm{e}^{-u} \mathrm{e}^{-v} \mathrm{e}^{i \alpha (\sin(\omega u + \phi) - \sin(\omega v + \phi))} \mathrm{d}u\  \mathrm{d}v\ \mathrm{d}\phi, $$
where $\alpha > 0$, $\omega > 0$.
A typical parameter regime I would be interested in for any approximate solutions would be $\omega \approx 2.5$, $\alpha \in [0, 3]$ (and in particular, $\alpha \ll 1$).
The result is real, and evaluated numerically across $\alpha \in [0, 10]$, looks roughly like a smeared-out Bessel function of the first kind, as you might expect from the sine in the exponent. Is there some clever contour integration trick for simplifying this? I didn't find anything using Mathematica/… (and a quick manual check, but Bessel functions, etc. were never my strength).
 A: Using a trigonometric identity, we can rewrite the integral as
$$ I (\alpha,\omega) = \int \limits_0^\infty \int \limits_0^\infty \int \limits_0^{2\pi} \mathrm{e}^{-(u+v)} \mathrm{e}^{2 \mathrm{i} \alpha \sin (\omega (u-v)/2) \cos(\phi + \omega (u+v)/2)} \, \mathrm{d} \phi \, \mathrm{d} u \, \mathrm{d} v \, . $$
The shift in the cosine can be ignored as the integration is over a whole period anyway, so the integral representation of the Bessel function of the first kind yields
$$ I(\alpha,\omega) = 2 \pi \int \limits_0^\infty \int \limits_0^\infty \mathrm{e}^{-(u+v)} \operatorname{J}_0 \left[2 \alpha \sin \left(\frac{\omega(u-v)}{2}\right)\right] \, \mathrm{d} u \, \mathrm{d} v \, . $$
Now we change variables to $x = u - v$ and $y = u + v$. The Jacobian is $1/2$ and the new integration region is $\{(x,y) \in \mathbb{R}^2 \colon y > \lvert x \rvert\}$. Just like the integrand, it is symmetric in $x$, so we can integrate over the part with $x>0$ twice to find
\begin{align}
I (\alpha, \omega) &= 2 \cdot 2 \pi \cdot \frac{1}{2} \int \limits_0^\infty \int \limits_x^\infty \mathrm{e}^{-y} \operatorname{J}_0 \left[2 \alpha \sin \left(\frac{\omega x}{2}\right)\right] \, \mathrm{d} y \, \mathrm{d} x \\
&= 2 \pi \int \limits_0^\infty \mathrm{e}^{-x} \operatorname{J}_0 \left[2 \alpha \sin \left(\frac{\omega x}{2}\right)\right] \, \mathrm{d} x \, .
\end{align}
While a closed-form expression for the remaining integral seems unlikely, it can easily be approximated. The series expansion of $\operatorname{J}_0$ leads to
$$ I(\alpha,\omega) = 2\pi \left[1 - \frac{\omega^2 \alpha^2}{2(1+\omega^2)} + \mathcal{O} \left(\alpha^4\right)\right] $$
for example and higher-order terms are readily computed.
