Is there a closed form for the integral $\int_0^\infty \frac{e^{-x^2} I_0 \left(\beta x \right) d x}{\sqrt{ \alpha^2+x^2}}$ I encountered this integral in my work, and it would be really convenient if it had a closed form in terms of any known special functions (which Mathematica could handle):
$$J(\alpha,\beta)=\int_0^\infty \frac{e^{-x^2} I_0 \left(\beta x \right) d x}{\sqrt{ \alpha^2+x^2}}$$
It's easy to see that:
$$J(\alpha,0)=\frac12 e^{\alpha^2/2} K_0 \left( \frac{\alpha^2}{2} \right)$$
Using series expansion of the Bessel function and a little help from Mathematica, I was able to get a series expression:
$$J(\alpha,\beta)=\frac12 e^{\alpha^2/2} K_0 \left( \frac{\alpha^2}{2} \right)+\frac{\sqrt{\pi}}{16} \beta^2 \sum_{k=0}^\infty \frac{(2k+1)!}{k!^3 (k+1)^2} U \left(\frac{1}{2},-k,\alpha^2 \right) \frac{\beta^{2k}}{4^{2k}}$$
Where $U$ is the confluent hypergeometric function, which in this case (according to Mathematica) is a linear combination of Bessel functions $K_0$ and $K_1$ though I wasn't able to get the general expression for the coefficients yet.
This series is good for small values of $\beta$, but for large $\beta$ I need too many terms, and the computation time is slower than I'd like.
I know that it's simple to get a quadrature formula for the integral, but still, a closed form would be better.

The integral can be thought as a Hankel transform with imaginary argument, but I wasn't able to find it in the tables of Hankel transforms either.

Clarification re: R. Burton's comment:
The final expression has the form:
$$\alpha e^{-\beta^2/2} J(\alpha,\beta)$$
Which should make the expression finite for every choice of variables.

Now that I think about it, I should probably use asymptotics for the Bessel function for large $\beta$:
$$I_0 \left(\beta x \right) \asymp \frac{1}{\sqrt{2 \pi \beta x}}e^{\beta x}$$
Then I get:
$$J(\alpha,\beta) \asymp \frac{1}{\sqrt{2 \pi \beta}} \int_0^\infty \frac{e^{-x^2+\beta x}  d x}{\sqrt{x}\sqrt{ \alpha^2+x^2}}, \quad \beta \to \infty$$
Which I still don't know the closed form for.

A more simple integral which does have a closed form:
$$\int_0^\infty e^{-x^2} I_0 \left(\beta x \right) d x= \frac{\sqrt{\pi}}{2} e^{\beta^2/8} I_0  \left(\frac{\beta^2}{8} \right)$$
 A: Changing $x=\beta t$, one can express the asymptotic expression as
  \begin{align}
   J(\alpha,\beta) &\asymp \frac{1}{\sqrt{2 \pi \beta}} \int_0^\infty \frac{e^{-x^2+\beta x}  d x}{\sqrt{x}\sqrt{ \alpha^2+x^2}}, \quad \beta \to \infty\\
    &\asymp \frac{1}{\sqrt{2 \pi }} \int_0^\infty\frac{e^{-\beta^2\left( t^2-t \right)}}{\sqrt{\alpha^2+\beta^2t^2}}\frac{dt}{\sqrt{t}}
  \end{align}
The exponential term reaches its maximum when $t=1/2$, while the remaining part of the integrand behaves smoothly. One may use the Laplace method:
\begin{equation}
 \int_a^b h(t)e^{Mg(t)}\,dt\asymp\sqrt{\frac{2\pi}{M\left|g''(t_0)\right|}}h(t_0)e^{Mg(t_0)}, \quad M \to \infty
\end{equation} 
here $M=\beta^2,g(t)=-t^2+t,h(t)=t^{-1/2}\left( \alpha^2+\beta^2t^2 \right)^{-1/2},t_0=1/2,g(t_0)=1/4,g''(t_0)=-2$
\begin{equation}
   J(\alpha,\beta) \simeq \frac{2e^{\frac{\beta^2}{4}}}{\beta\sqrt{4\alpha^2+\beta^2}}
\end{equation} 
For $\alpha=1.2345,\beta=10$, we find numerically $\alpha e^{-\beta^2/2}J(\alpha,\beta)=3.4603.10^{-13}$, while the approximation gives $3.3289^{-13}$. Additional terms, both in the $I_0$ expansion and (more laboriously) in the Laplace method, may be added in this way.
