Xrfgtjtk

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)$$

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^inftyfrac{e^{-beta^2left( 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)},dtasympsqrt{frac{2pi}{Mleft|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}}}{betasqrt{4alpha^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.