<s>[INST] Integral Involving Bessel function and trigonometric function Consider the following integral
\begin{equation}
\int_0^{2 \pi} \int_0^{\pi} e^{A \cos \phi \cos \theta} \, I_0(B \sin \phi \sin \theta) \sin \theta \, d\theta \, d\phi
\label{1}
\end{equation}
where $A$ and $B$ are two constants and $I_0(.)$ is the modified Bessel function of the first kind.
To evaluate the integral, I tried to write the expansion of $I_0$ which (by Methemtica) results in
$$
\sum_{n=0}^{\infty} \frac{1}{4^n} \frac{1}{(n!)^2} \int_0^{2 \pi} \int_0^{\pi} e^{A \cos \phi \cos \theta} \, (B \sin \phi \sin \theta)^{2n} \sin \theta \, d\theta \, d\phi
$$
$$
= \sum_{n=0}^{\infty} \frac{1}{4^n} \frac{1}{(n!)^2} \sqrt{\pi} \, \Gamma(n+1) \int_0^{2 \pi} \operatorname{HGR}_{01}\left(n+\frac{3}{2}, \frac{1}{4} A^2 \cos^2 \phi\right) \, (B^2 \sin^2 \phi)^{n} \sin \theta  \, d\phi
$$
where $\operatorname{HGR}_{01}$ is the regularized confluent hypergeometric function.
Again, by expanding this function I found (up to some constants)
$$
\sum_{n=0}^{\infty}\sum_{m=0}^{\infty} A'^m B'^n \frac{\Gamma(n+\frac{1}{2}) \Gamma(m+\frac{1}{2})\Gamma(m + n + \frac{3}{2})}{\Gamma(n+1) \Gamma(m+1)\Gamma(m+n+1)}
$$
where $A',B'$ are (properly) rescaled version of $A,B$.
I couldn't go further, and stucked with this double sum.
I appreciate if anyone can provide me some insights on this sum, or any hints on the original integral .
 [/INST] In this paper by VNP Anghel, a result obtained by Glasser is generalized (eq. 32):
\begin{equation}
 I_{\frac{m+1}{2}}(b)=\left( \frac{b}{2\pi \sin^m\alpha}\right)^{1/2} \int_0^\pi \exp(b\cos\alpha\cos\theta)I_{\frac m2}(b\sin\alpha\sin\theta)(\sin\theta)^{\frac{m+2}{2}}\,d\theta \tag{1}\label{eq1}
\end{equation}
Then, in the case $A=B=b$, by choosing $m=0$, we have directly
\begin{align}
 J(b,b)&=\int_0^{2 \pi} \int_0^{\pi} e^{b \cos \phi \cos \theta} \, I_0(b \sin \phi \sin \theta) \sin \theta \, d\theta \, d\phi\\
 &=\int_0^{2 \pi} \, d\phi \left( \frac{2\pi}{b} \right)^{1/2}I_{1/2}(b)\\
 &=4\pi\frac{\sinh b}{b}
\end{align}
The result \eqref{eq1} may also be used in the case $A\ne B$ to express the integral as a series. We remark first that $A$ can be chosen to be positive (parity of the integral wrt $A$ is clear from the substitution $\phi\to \pi-\phi$ in the integral). By the multiplication theorem,
\begin{equation}
 I_{\nu}\left(\lambda z\right)=\lambda^{\nu}\sum_{k=0}^{\infty}
\frac{(\lambda^{2}-1)^{k}(\frac{1}{2}z)^{k}}{k!}I_{\nu+ k}\left(z
\right)
\end{equation}
with $\nu=0,\lambda=B/A,z=A\sin\phi\sin\theta$, one can express
\begin{equation}
 I_0\left(B\sin\phi\sin\theta\right)=\sum_{k=0}^{\infty}
\frac{\left( B^2-A^2 \right)^k}{2^kk!A^k}\sin^k\phi\sin^k\theta \,I_{ k}\left(A\sin\phi\sin\theta\right)
\end{equation}
and thus, by interverting integral and summation,
\begin{align}
 J(A,B)&=\int_0^{2 \pi} \int_0^{\pi} e^{A \cos \phi \cos \theta} \, I_0(B \sin \phi \sin \theta) \sin \theta \, d\theta \, d\phi\\
    &=\sum_{k=0}^{\infty}\frac{\left( B^2-A^2 \right)^k}{2^kk!A^k}\int_0^{2 \pi}\sin^k\phi\,d\phi
    \int_0^{\pi} e^{A \cos \phi \cos \theta}I_{ k}\left(A\sin\phi\sin\theta\right) \sin^{k+1}\theta  \, d\theta
\end{align}
With $m=2k$ in eq. \eqref{eq1}, it may be expressed as
\begin{align}
 J(A,B)&=\sum_{k=0}^{\infty}\frac{\left( B^2-A^2 \right)^k}{2^kk!A^k}\int_0^{2 \pi}\left( \frac{2\pi\sin^{2k}\phi}{A} \right)^{1/2}\sin^k\phi \,I_{k+1/2}(A)\,d\phi\\
 &=\sqrt{\frac{2\pi}{A}}\sum_{k=0}^{\infty}\frac{\left( B^2-A^2 \right)^k}{2^kk!A^k}I_{k+1/2}(A)\int_0^{2 \pi}\sin^{2k}\phi\,d\phi\\
 &=\frac{(2\pi)^{3/2}}{\sqrt{A}}\sum_{k=0}^{\infty}\frac{(2k)!}{2^{3k}(k!)^3}\frac{\left( B^2-A^2 \right)^k}{A^k}I_{k+1/2}(A)
\end{align}
Numerical experiments show that this series converges very quickly, which is related to the asymptotic expansion of the modified Bessel function for large orders (DLMF).
</s>