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.

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