Integral over Modified Bessel Functions I've stumbled upon this integral
\begin{equation}
\mathscr{I}_\nu^k(a,b)=\int_0^\infty dr\frac{r^k}{r^2+m^2}e^{-ar^2}I_\nu(br)
\end{equation}
(where $I_\nu(x)$ is the modified Bessel function) during some QFT research, but I cannot seem to crack it. I'm specifically interested in the integrals $\mathscr{I}_1^2$, $\mathscr{I}_1^4$, and $\mathscr{I}_2^3$, each of which have an odd integrand, so that the typical semicircular contours are seemingly inutile. I've tried Weierstrass transforms, Hankel transforms, Laplace transforms, and the like, but each made the evaluation even more complicated. Additionally, I cannot find any similar integrals in the literature. Does anyone know a solution or a plan of attack? Any input would be greatly appreciated.
 A: Let's do the $(1,2)$ case. First let's denote
$$J(a,b) = \int_0^\infty \frac{r^2}{r^2+m^2}e^{-ar^2}I_1(br)dr = \int_0^\infty e^{-ar^2}I_1(br)dr - \int_0^\infty \frac{m^2}{r^2+m^2}e^{-ar^2}I_1(br)dr$$
Which means we have that
$$\partial_a J(a,b) = \int_0^\infty -r^2e^{-ar^2}I_1(br)dr + \int_0^\infty \frac{m^2r^2}{r^2+m^2}e^{-ar^2}I_1(br)dr $$
$$\implies \partial_a J(a,b) = m^2 J(a,b) - \int_0^\infty r^2e^{-ar^2}I_1(br)dr$$
which gives us a nice ODE with homogeneous solution $J_h(a,b) = B(b)e^{m^2 a}$. Using variation of parameters (i.e. assuming the particular solution is of the form $v\cdot J_h$) gives
$$\partial_a v = - e^{-m^2a}\int_0^\infty r^2e^{-ar^2}I_1(br)dr$$
We'll focus on the integral inside first. Using the series representation for the Modified Bessel function of the first kind, we have that 
$$\int_0^\infty r^2e^{-ar^2}I_1(br)dr = \frac{1}{2}\sum_{k=0}^\infty \frac{4^{-k}b^{2k+1}}{k!\Gamma(k+2)} \int_0^\infty r^{2k+3}e^{-ar^2}dr$$
Then using the substitution $s=\sqrt{a}r$ we get
$$\frac{1}{2a^{\frac{3}{2}}}\sum_{k=0}^\infty \frac{4^{-k}\left(\frac{b}{\sqrt{a}}\right)^{2k+1}}{k!\Gamma(k+2)} \int_0^\infty s^{2k+3}e^{-s^2}ds$$
$$ = \frac{1}{2a^{\frac{3}{2}}}\sum_{k=0}^\infty \frac{4^{-k}\left(\frac{b}{\sqrt{a}}\right)^{2k+1}}{k!\Gamma(k+2)} \cdot \frac{\Gamma(k+2)}{2} = \frac{b}{4a^2}e^{\frac{b^2}{4a}}$$
So we have 
$$ v = \int -e^{-m^2a}\cdot\left(\frac{b}{4a^2}e^{\frac{b^2}{4a}}\right)da = \frac{1}{b}e^{-m^2a}e^{\frac{b^2}{4a}} + \int \frac{m^2}{b}e^{-m^2a}e^{\frac{b^2}{4a}}da$$
Thus our (almost) final answer is
$$ \mathscr{I}_1^2(a,b) = B(b)e^{m^2a}+\frac{1}{b}e^{\frac{b^2}{4a}} + e^{m^2a}\int_1^a \frac{m^2}{b}e^{-m^2 \alpha}e^{\frac{b^2}{4\alpha}}d\alpha$$
where $B(b)$ can be determined by plugging in $a=1$ for the original integral and solving for the series.
What's nice is that once you have the $(1,2)$ case with whatever method you choose, you have the $(1,4)$ case for free because
$$\mathscr{I}_1^4(a,b) = -\partial_a \mathscr{I}_1^2(a,b)$$
A: We can transform
\begin{align}
\mathscr{I}_\nu^k(a,b)&=\int_0^\infty dr\frac{r^k}{r^2+m^2}e^{-ar^2}I_\nu(br)\\
&=e^{am^2}\int_0^\infty dr\frac{r^k}{r^2+m^2}e^{-a(r^2+m^2)}I_\nu(br)\\
&=e^{am^2}J_\nu^k(a,b)
\end{align} 
We have 
\begin{equation}
\frac{\partial J_\nu^k(a,b)}{\partial a}=-e^{-am^2}\int_0^\infty r^ke^{-ar^2}I_\nu(br)\,dr
\end{equation} 
with $\lim_{a\to\infty} J_\nu^k(a,b)=0$. When $k=\nu+1$, this integral is tabulated (DLMF):
\begin{equation}
\int_{0}^{\infty}t^{\nu+1}I_{\nu}\left(bt\right)\exp\left(-at^{2}\right)\,dt=\frac{b^\nu}{(2a)^{\nu+1}}\exp\left(\frac{b^{2}}{4a}\right)
\end{equation} 
then
\begin{equation}
\frac{\partial J_\nu^{\nu+1}(a,b)}{\partial a}=-\frac{b^\nu}{(2a)^{\nu+1}}\exp\left(\frac{b^{2}}{4a}-am^2\right)
\end{equation} 
We deduce
\begin{equation}
J_\nu^{\nu+1}(a,b)=\frac{b^{\nu}}{2^{\nu+1}}\int_a^\infty\exp\left(\frac{b^{2}}{4\alpha}-\alpha m^2\right)\frac{d\alpha}{\alpha^{\nu+1}}
\end{equation} 
Thus
\begin{align}
\mathcal{I}_\nu^{\nu+1}(a,b)&=\frac{b^{\nu}e^{am^2}}{2^{\nu+1}}\int_a^\infty\exp\left(\frac{b^{2}}{4\alpha}-\alpha m^2\right)\frac{d\alpha}{\alpha^{\nu+1}}\\
&=\frac{b^{\nu}e^{am^2}}{2^{\nu+1}}\int_0^{1/a}\exp\left(\frac{b^{2}x}{4}-\frac{m^2}{x}\right)x^{\nu-1}\,dx
\end{align}
For the cases of interest,
\begin{align}
\mathscr{I}_1^{2}(a,b)&=\frac{be^{am^2}}{4}\int_0^{1/a}\exp\left(\frac{b^{2}x}{4}-\frac{m^2}{x}\right)\,dx\\
\mathscr{I}_2^{3}(a,b)&=\frac{b^{2}e^{am^2}}{8}\int_0^{1/a}\exp\left(\frac{b^{2}x}{4}-\frac{m^2}{x}\right)x\,dx
\end{align}
and, as noted in @NinadMunshi answer,
\begin{equation}
\mathscr{I}_1^{4}(a,b)=-\partial_a\mathscr{I}_1^{2}(a,b)
\end{equation} 
Using derivations and another tabulated integral (DLMF), it seems possible to obtain  integral expressions for other values of $\nu$ and $k$. It should also be noticed that obtained integrals could be expressed in terms of the  incomplete Bessel functions
