integral involving modified bessel function I would like to put in closed form the integral:
$\int{e^{-k x} I_0(x) dx } $
where $I_\alpha(x)$ is the modified Bessel function of the first kind.
The closest I have found in tables is for k=1
$\int{e^{-x} I_0(x) dx } =  x e^{-x} (I_0(x)+I_1(x))$
It would be interesting to see if it is solvable,  at least for a numerable set of k values.  
The integral arises in a rather fundamental problem: the probability of being inside a ball of defined radious  for a bivariate normal. So if it is not solvable, perhaps it worth to define a ad hoc function for it.
Thanks for your interest
 A: For any $\alpha\in\mathbb{N}$ we have:
$$ I_\alpha(x) = \frac{1}{\pi}\int_{0}^{\pi}e^{x\cos\theta}\cos(\alpha\theta)\,d\theta,$$
$$ I_\alpha'(x) = \frac{1}{\pi}\int_{0}^{\pi}e^{x\cos\theta}\cos(\alpha\theta)\cos\theta\,d\theta=\left\{\begin{array}{rcl}\frac{1}{2}\left(I_{\alpha-1}(x)+I_{\alpha+1}(x)\right)&\text{if}&\alpha\geq 1,\\I_1(x)&\text{if}&\alpha=0.\end{array}\right.$$
$$I_{\alpha+1}(x)=\frac{x}{2\alpha+2}\left(I_{\alpha}(x)-I_{\alpha+2}(x)\right).$$
With these identities, it is straightforward to check that the primitive of $I_0(x)e^{-x}$ is $xe^{-x}(I_0(x)+I_1(x))$. Differentiating both sides, we have that the primitive of $(I_1(x)-I_0(x))e^{-x}$ is $I_0 e^{-x}$, hence the primitive of $I_1(x)e^{-x}$ is $I_0(x) e^{-x}+xe^{-x}(I_0(x)+I_1(x))$. Along the same lines, we can compute the primitive of $I_{\alpha}(x)e^{-x}$ with a recursive argument.
Since for any $k>1$ we have:
$$\int f(x)e^{-kx}\,dx = -e^{-kx}\sum_{j=0}^{+\infty}\frac{f^{(j)}(x)}{k^{j+1}}$$
it follows that:
$$\mathcal{J}_k(x)=\int I_0(x)e^{-kx}=-\frac{e^{-kx}}{\pi}\int_{0}^{\pi}e^{-x\cos\theta}\sum_{j=0}^{+\infty}\frac{\cos^j\theta}{k^{j+1}}\,d\theta=-\frac{e^{-kx}}{\pi}\int_{0}^{\pi}\frac{e^{-x\cos\theta}}{k-\cos\theta}\,d\theta,$$
or, by exploiting the Fourier cosine series of $e^{-x\cos\theta}$:
$$-\pi e^{kx}\mathcal{J}_k(x)=\int_{0}^{\pi}\frac{I_0(x)}{k-\cos\theta}d\theta+2\sum_{m=1}^{+\infty}(-1)^m\int_{0}^{\pi}\frac{I_m(x)\cos(m\theta)}{k-\cos\theta}\,d\theta,$$
$$-e^{kx}\mathcal{J}_k(x)=\frac{I_0(x)}{\sqrt{k^2-1}}+2\sum_{j=1}^{+\infty}(-1)^j I_j(x)\int_{-\infty}^{+\infty}\frac{\cos(2j\arctan u)}{k(1+u^2)-(1-u^2)}\,du$$
where the last integrals can be computed through the residue theorem and they always equal
$$\left(A+\frac{B}{\sqrt{k^2-1}}\right)$$
for some $A,B\in\mathbb{Z}$.
A: Start with
$$
\int e^{-k x}\;dx = -\frac{e^{-k x}}{k}
$$
by inspection
$$
\frac{d^n}{dk^n}\left[-\frac{e^{-k x}}{k}\right]=\sum_{i=1}^{n+1}\frac{(-1)^{n+1}n!}{(n-i+1)!}\frac{x^{n+1-i}e^{-kx}}{k^i}
$$
we also have
$$
\frac{d^n}{dk^n}\left[e^{-k x}\right]=(-x)^n e^{-kx}
$$
generate the operator 
$$
\hat{O}=\sum_{n=0}^\infty \frac{1}{4^nn!^2}\frac{d^{2n}}{dk^{2n}}
$$
apply to both sides of the first equation
$$
\hat{O}\left[\int e^{-k x}\;dx\right] = \hat{O}\left[-\frac{e^{-k x}}{k}\right]
$$
$$
\int \sum_{n=0}^\infty \frac{1}{4^nn!^2}\frac{d^{2n}}{dk^{2n}}\left[e^{-k x}\right]\;dx = \sum_{n=0}^\infty \frac{1}{4^nn!^2}\frac{d^{2n}}{dk^{2n}}\left[-\frac{e^{-k x}}{k}\right]
$$
$$
\int \sum_{n=0}^\infty \frac{1}{4^nn!^2}(-x)^{2n}e^{-kx}\;dx = \sum_{n=0}^\infty \frac{1}{4^nn!^2}\sum_{i=1}^{2n+1}\frac{(-1)^{2n+1}(2n)!}{(2n-i+1)!}\frac{x^{2n+1-i}e^{-kx}}{k^i}
$$
we have 
$$
I_0(x)=\sum_{n=0}^\infty \frac{x^{2n}}{4^nn!^2}
$$
$$
\int e^{-kx}I_0(x)\;dx = -e^{-kx}\sum_{n=0}^\infty \frac{(2n)!x^{2n}}{4^nn!^2}\sum_{i=1}^{2n+1}\frac{x^{1-i}}{k^i(2n-i+1)!}
$$
$$
\int e^{-kx}I_0(x)\;dx = -e^{-kx}\sum_{n=0}^\infty \frac{(2n)!x^{2n}}{4^nn!^2}\left(-\frac{x}{(1+2n)!}+e^{kx}\frac{\Gamma(2+2n,k x)}{k(kx)^{2n}(1+2n)!} \right)
$$
$$
\int e^{-kx}I_0(x)\;dx = xe^{-kx}\,_1F_2\left(\frac{1}{2};1,\frac{3}{2},\frac{x^2}{4}\right) - \sum_{n=0}^\infty \frac{\Gamma(2+2n,k x)}{4^nn!^2 (1+2n)k^{2n+1}}
$$
The incomplete gamma function does not help reduce the last sum to a hypergeometric function.
