I try to solve $\int_0^T e^{-x}\frac{I_n(\alpha x)}{x}dx$ where $I_n(x)$ is the modified Bessel function of the first kind and $0<\alpha<1$.

My first approach was to turn this integral into an infinite sum to fit a hypergeometric series:

• Using the infinite series representation of the Bessel function, I got incomplete gamma functions in the sum, which does not sound promising.
• The multiplication theorem yields an infinite series of integrals where we get rid of the $\alpha$: $$\int_0^T e^{-x}\frac{I_n(\alpha x)}{x}dx=\alpha^n\sum_{m=0}^{\infty}\frac {\big(\frac {\alpha ^{2}-1}{2}\big)^m}{m!}\int_0^T e^{-x}x^{m-1}I_{n+m}(x)dx$$

According to a table of integrals, the new integrals are: \begin{align} \int_0^T e^{-x}x^{m-1}I_{n+m}(x)dx&=\Big(\frac{T}{2}\Big)^{n+m}\frac{\Gamma(2m+n)}{\Gamma(m+n+1)\Gamma(2m+n+1)}\\ &\times{}_2F_2[\{m+n+\frac{1}{2},2m+n\};\{2m+2n,2m+n+1\};-2T] \end{align}

Expanding ${}_2F_2$ (let's call $k$ the summation index), we get a double infinite series which might fit the definition of a hypergeometric function of 2 variables. However, I get several Pochhammer symbols with coupled summations indices like this: $\frac{(n)_{2m+k}}{(n+1)_{2m+k}(2n)_{2m+k}}$ which does not fit any hypergeometric function definition.

Another approach could be to get inspiration from the limit $T\rightarrow \infty$ which is the Laplace transform of $\frac{I_n(x)}{x}$ (up to a constant) and it has a closed form (according to a table): $$\int_0^\infty e^{-x}\frac{I_n(\alpha x)}{x}dx=\frac{\big(\frac{\alpha}{1+\sqrt{1-\alpha^2}}\big)^n}{n}$$

However, I don't find any reference on the way to compute this.

Do you have any suggestion about how to compute it?