Lets start from the definition of the modified Bessel Function. Recall
$$I_{0}(z)=\sum_{n=0}^{\infty}\frac{1}{4^{n}}\frac{z^{2n}}{n!n!}.$$ So our integral is
$$\int_{0}^{\infty}e^{-\frac{1}{2}\left(x^{2}+b^{2}\right)}x\sum_{n=0}^{\infty}\frac{1}{4^{n}}\frac{b^{2n}x^{2n}}{n!n!}dx.$$
Switching the order of sumation and integration we get
$$e^{-\frac{1}{2}b^{2}}\sum_{n=0}^{\infty}\frac{b^{2n}}{4^{n}n!n!}\int_{0}^{\infty}e^{-\frac{1}{2}x^{2}}x^{2n+1}dx.$$
Now, let $u=\frac{1}{2}x^{2}$ and $du=xdx$, to see that
$$\int_{0}^{\infty}e^{-\frac{1}{2}x^{2}}x^{2n+1}dx=2^{n}\int_{0}^{\infty}e^{-u}u^{n}du=2^{n}\Gamma(n+1)=2^{n}n!.$$
Hence
$$e^{-\frac{1}{2}b^{2}}\sum_{n=0}^{\infty}\frac{b^{2n}}{4^{n}n!n!}\int_{0}^{\infty}e^{-\frac{1}{2}x^{2}}x^{2n+1}dx=e^{-\frac{1}{2}b^{2}}\sum_{n=0}^{\infty}\frac{b^{2n}}{4^{n}n!n!}\left(2^nn!\right)$$
$$=e^{-\frac{1}{2}b^{2}}\sum_{n=0}^{\infty}\frac{\left(\frac{b^{2}}{2}\right)^{n}}{n!}=e^{-\frac{1}{2}b^{2}}e^{\frac{1}{2}b^{2}}=1$$
as desired.
Notice that the exact same argument shows
$$\int_0^\infty e^{\left(\frac{1}{2}b^2-\frac{1}{2}x^2\right)}J_0(bx)xdx=1$$ where $J_0(x)$ is the Bessel Function.
Hope that helps,
One can try to do the following derivations.
$$\mathrm{Int}=\frac{1}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty} \, x \, e^{(-\frac{x}{\alpha}-\frac{1}{2w^2}(x-\hat{x})^2)} \, I_0\left(\frac{x}{\beta}\right)\,dx$$
You can simplify the power of the exponent:
$$-\left(\frac{x}{\alpha}+\frac{1}{2w^2}(x-\hat{x})^2\right)=-\left(\frac{(x-\mu)^2}{2w^2}+\gamma \right),$$
where $\gamma=\left(\frac{\hat{x}}{\alpha}-\frac{w^2}{2\alpha^2}\right)$, $\mu=\hat{x}-\frac{w^2}{\alpha^2}$. So:
$$\mathrm{Int}=\frac{1}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty} \, x \, e^{(-\frac{x}{\alpha}-\frac{1}{2w^2}(x-\hat{x})^2)} \, I_0\left(\frac{x}{\beta}\right)\,dx=\frac{e^{-\gamma}}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty} \, x \, e^{-\frac{(x-\mu)^2}{2w^2}} \, I_0\left(\frac{x}{\beta}\right)\,dx$$
Then you can change the variable of integration $y=x-\mu$:
$$\mathrm{Int}=\frac{e^{-\gamma}}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty} \, (y+\mu) \, e^{-\frac{y^2}{2w^2}} \, I_0\left(\frac{y+\mu}{\beta}\right)\,dx$$
and make use of the Neumann’s addition theorem:
$$\mathop{I_{{\nu}}}\nolimits\!\left(u\pm v\right)=\sum _{{k=-\infty}}^{\infty}(\pm 1)^{k}\mathop{I_{{\nu+k}}}\nolimits\!\left(u\right)\mathop{I_{{k}}}\nolimits\!\left(v\right)$$
setting $\nu=0$, assuming that $y,\mu,\beta\in\mathbb{R}$ and using connection formulas one will get:
$$\mathop{I_{0}}\nolimits\!\left(\frac{y+\mu}{\beta}\right)=\sum _{{k=-\infty}}^{\infty}\mathop{I_{k}}\nolimits\!\left(\frac{y}{\beta}\right)\mathop{I_{{k}}}\nolimits\!\left(\frac{\mu}{\beta}\right)=\mathop{I_{0}}\nolimits\!\left(\frac{y}{\beta}\right)\mathop{I_{{0}}}\nolimits\!\left(\frac{\mu}{\beta}\right)+2\sum _{{k=1}}^{\infty}\mathop{I_{k}}\nolimits\!\left(\frac{y}{\beta}\right)\mathop{I_{{k}}}\nolimits\!\left(\frac{\mu}{\beta}\right)$$
$$\mathrm{Int}=\!
\frac{\mathop{I_0}\nolimits\!\left(\frac{\mu}{\beta}\right)e^{-\gamma}}{\sqrt{2\pi w^2}}\!\int_{-\infty}^{+\infty} \, (y+\mu)\! e^{-\frac{y^2}{2w^2}}\! \mathop{I_{0}}\nolimits\!\left(\frac{y}{\beta}\right)\!dx\!
+\!\frac{2 e^{-\gamma}}{\sqrt{2\pi w^2}}\int_{-\infty}^{+\infty}\! (y+\mu) \, e^{-\frac{y^2}{2w^2}} \! \sum _{k=1}^{\infty}\!\mathop{I_k}\!\!\left(\frac{y}{\beta}\right)\!\!\mathop{I_k}\nolimits\!\left(\frac{\mu}{\beta}\right) \,dx$$
Or reorganising terms:
$\mathop{I_0}\nolimits\!\left(\frac{\mu}{\beta}\right)$
$$\mathrm{Int}=\!
\frac{e^{-\gamma}}{\sqrt{2\pi w^2}}\!\left(\mathop{I_0}\nolimits\!\left(\frac{\mu}{\beta}\right)\int_{-\infty}^{+\infty} \, (y\!+\!\mu)\! e^{-\frac{y^2}{2w^2}}\! \mathop{I_{0}}\nolimits\!\left(\frac{y}{\beta}\right)\!dx\!
+\!2\!\sum _{k=1}^{\infty}\!\mathop{I_k}\!\left(\frac{\mu}{\beta}\right)\!\int_{-\infty}^{+\infty}\! (y\!+\!\mu) \, e^{-\frac{y^2}{2w^2}} \! \!\mathop{I_k}\!\left(\frac{y}{\beta}\right)\!\!\,dx \right)$$
Let's set
$$\mathrm{Int}_k(\mu,\alpha,\beta)\!=\!\int_{-\infty}^{+\infty} \, (y\!+\!\mu) e^{-\frac{y^2}{2w^2}}\! \mathop{I_{k}}\nolimits\!\left(\frac{y}{\beta}\right)\!dx$$
$$\mathrm{Int}_0(\mu,\alpha,\beta)\!=\!\int_{-\infty}^{+\infty} \, (y\!+\!\mu) e^{-\frac{y^2}{2w^2}}\! \mathop{I_{0}}\nolimits\!\left(\frac{y}{\beta}\right)\!dx$$
Here one can play with the evenness/unevenness of the functions under integrals (as I did at first), or just use CAS and get:
$$\mathrm{Int}_k(\mu,\alpha,\beta)\!=\!\sqrt{\frac{\pi \alpha ^2}{8 \beta ^2}}\! e^{\frac{\alpha ^2}{4 \beta ^2}} \left(\!\!2 \beta \left(\!\!(-1)^k+1\!\right)\! \mu I_{\frac{k}{2}}\!\left(\!\!\frac{\alpha ^2}{4 \beta ^2}\!\!\right)+\alpha ^2 \left(\!1-(-1)^k\!\right)\!\! \left(\!\!I_{\frac{k-1}{2}}\!\!\left(\!\!\frac{\alpha ^2}{4 \beta ^2}\!\!\right)+I_{\frac{k+1}{2}}\!\left(\!\!\frac{\alpha ^2}{4 \beta ^2}\!\!\right)\!\!\right)\!\!\right)$$
$$\mathrm{Int}_0(\mu,\alpha,\beta)\!=\! \sqrt{2 \pi } \alpha \mu e^{\frac{\alpha ^2}{4 \beta ^2}} I_0\left(\frac{\alpha ^2}{4 \beta ^2}\right)$$
So the original integral can be represented in the following way:
$$\mathrm{Int}(\mu,\alpha,\beta)\!=\! \frac{e^{-\gamma}}{\sqrt{2\pi w^2}}\left(\mathrm{Int}_0(\mu,\alpha,\beta)+
2\!\sum _{k=1}^{\infty}\mathop{I_k}\left(\frac{\mu}{\beta}\right)
\mathrm{Int}_k(\mu,\alpha,\beta)
\right)$$
And I am not shure that there is any chance to simplify it futher (meaning finding the sum of the series).
Best Answer
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).