Integrate Product of Modified Bessel Functions of First and Second Kind I need the following integration to solve a much larger problem:
$$\int rI_0(kr)K_0(kr)dr$$
I've looked through the typical integration tables but haven't found this particular combination. Maybe I've somehow missed it? Seems like this would be well known.
I've tried integrating by parts but it doesn't seem to get anywhere:
$$\int I_0(kr)rK_0(kr)dr=\frac{1}{k}I_0(kr)K_1(kr)+\int I_1(kr)rK_1(kr)$$
Thanks!
 A: Coming from a CAS,
$$\int r\,I_0(kr)\, K_0(kr)\,dr=\frac{r^2}{4 \sqrt{\pi }}\,\, G_{1,3}^{2,1}\left(k r,\frac{1}{2}|
\begin{array}{c}
 \frac{1}{2} \\
 0,0,-1
\end{array}
\right)$$ where appears Meijer G-function (it is very close to linearity).
A: From the integral representation
\begin{equation}
I_{\mu}\left(x\right)K_{\nu}\left(x\right)=\int_{0}^{\infty}J_{\mu\pm\nu}\left(2x\sinh t\right)e^{(-\mu\pm\nu)t}\mathrm{d}t
\end{equation} 
we can express, by changing the order of integration
\begin{align}
\int r\,I_0(kr)\, K_0(kr)\,dr&=\frac 1 {k^2}\int x\,I_0(x)\, K_0(x)\,dx\\
&=\frac1{k^2}\int x\,dx\int_{0}^{\infty}J_0\left(2x\sinh t\right)\,dt\\
&=\frac1{k^2}\int_0^\infty\,dt\int_0^\infty xJ_0\left(2x\sinh t\right)\,dx
\end{align}
Now, the $x$ integration is simple:
\begin{equation}
\int r\,I_0(kr)\, K_0(kr)\,dr=\frac1{k^2}\int_0^\infty \frac{xJ_1\left(2x\sinh t\right)}{2\sinh t}\,dt
\end{equation} 
We use the recurrence relation for the Bessel functions
\begin{equation}
2\frac{J_1(X)}{X}=J_0(X)+J_1(X)
\end{equation}
to write
\begin{equation}
\frac{J_1\left(2x\sinh t\right)}{2\sinh t}=\frac{1}{2x}\left(  J_0\left(2x\sinh t\right)+J_2\left(2x\sinh t\right)\right)
\end{equation} 
and by using twice the integral representation given above with the positive sign, $\mu=\nu=0$ and $\mu=\nu=1$:
\begin{align}
\int r\,I_0(kr)\, K_0(kr)\,dr&=\frac{x^2}{2k^2}\left[
\int_0^\infty J_0\left(2x\sinh t\right)\,dt+\int_0^\infty J_2\left(2x\sinh t\right)\,dt
\right]\\
&=\frac{x^2}{2k^2}\left[I_0(x)K_0(x)+I_1(x)K_1(x)\right]
\end{align}
Reintroducing the original variables:
$$
\int r\,I_0(kr)\, K_0(kr)\,dr=\frac{r^2}2\left[I_0(kr)K_0(kr)+I_1(kr)K_1(kr)\right]
$$
which can be checked by direct derivation. However I could not find the correspondance with the Meijer function proposed by @ClaudeLeibovici
A: We can actually do this with a fairly simple manipulation of the differential equations satisfied by $u_a(x) = I_0(ax)$ and $v_a(x) = K_0(ax)$. These are both solutions to
$$ (xy')' = a^2 x y . $$
Hence
$$ (xu_a')'v_b - (xv_b')' v_a = (a^2-b^2) x u_a v_b , $$
and the left-hand side is
$$ xu_a'' v_b - x u_a v_b'' + u_a''v_b - u_a v_b'' = (x(u_a' v_b - u_a v_b'))' $$
Integrating thus gives
$$ \int x u_a(x) v_b(x) \, dx = \frac{x(u_a' v_b - u_a v_b')}{a^2-b^2}  $$
But we want $a=b=1$. Applying L'Hopital's rule to $a \to b$ gives
$$ \int x u_b(x) v_b(x) \, dx = \frac{x}{2b} \left. \partial_a (u_a' v_b - u_a v_b') \right|_{a=b} . $$
We calculate
$$ \begin{align}
u_a(x) &= I_0(ax) & v_a'(x) &= a I_0'(ax) \\
\left. \partial_a u_a(x) \right|_{a=b}  &= x I_0'(bx) & \left. \partial_a u_a'(x) \right|_{a=b} &= I_0'(bx) + b I_0''(bx) ,
\end{align} $$
so
$$ \int x I_0(bx)K_0(bx) \, dx = \frac{x}{2b} ( (I_0'(bx) + b I_0''(bx)) K_0(bx) - x I_0'(bx) K_0'(bx)) , $$
and of course $ I'_0(bx) = b I_1(bx) $, $K_0'(bx) = -K_1(bx)$ and $I_0'(bx) + b I_0''(x) = bx I_0(bx) $ using the differential equation, so
$$ \int x I_0(bx)K_0(bx) \, dx = \frac{x^2}{2} (I_0(bx)K_0(bx)+ I_1(bx) K_1(bx)) . $$
Of course the same will work with other combinations of Bessel functions of the same order.
