# Simplification of integration of product of Bessel functions

I have been trying to evaluate the following integral involving modified Bessel functions:

\begin{align} \int r I_1 (r) K_1(r) dr \end{align}

This integral has an explicit expression given in the Wolfram documentation in terms of other Bessel functions. Granted, this expression has divergent components for integer-order Bessel functions but still has a well-defined limit at least for the above case. After performing these manually, one gets that the integral evaluates to

1/2 (-1 + r BesselI[0, r] (r BesselK[0, r] + BesselK[1, r]) +
r BesselI[1, r] (-BesselK[0, r] + r BesselK[1, r]))


However, when evaluating in Mathematica, this returns

In[1]:= Integrate[r BesselI[1, r] BesselK[1, r], r]

Out[1]= MeijerG[{{1, 3/2}, {}}, {{1, 2}, {0, 0}}, r, 1/2]/(4 Sqrt[\[Pi]])


Furthermore, this Meijer-G function refuses to simplify with FunctionExpand.

Is there a way get Mathematica to simplify the output of the integration? As far as I can see, the Mathematica internals should already know how to perform this simplification (or at least return the indefinite integral not in Meijer-G form).

Edit:

Furthermore, if one tries to evaluate

Integrate[z BesselI[\[Nu], a z] BesselK[\[Nu], a z], z]


as in the docs, Mathematica doesn't actually apply the rule discussed above, instead expressing in terms of hypergeometricPFQ functions. Why isn't this integration rule ever being used?

• Moreover, Series[MeijerG[{{1, 3/2}, {}}, {{1, 2}, {0, 0}}, r, 1/2]/( 4 Sqrt[\[Pi]]), {r, 0, 2}] produces two "Power::infy: Infinite expression 1/0 encountered." and $$\frac{G_{2,4}^{2,2}\left(0,\frac{1}{2}| \begin{array}{c} 1,\frac{3}{2} \\ 1,2,0,0 \\ \end{array} \right)}{4 \sqrt{\pi }}+\text{ComplexInfinity} r+\text{ComplexInfinity} r^2+O\left(r^3\right).$$ Aug 23, 2023 at 5:45

You could try

D[Series[
Integrate[r BesselI[n,a r] BesselK[n,a r],r] //FunctionExpand, {n,1,0}]
//Normal//FullSimplify,r]//FullSimplify
(* r BesselI[1, a r] BesselK[1, a r]  *)


Edit:

By trying different values of n the general case seems to be

Integrate[r*BesselI[n,a*r]*BesselK[n,a*r],r]=
-(n/(2*a^2))+(1/2)*r^2*(BesselI[n-1,a*r]*BesselK[n-1,a*r]+
BesselI[n,a*r]*BesselK[n,a*r])+((n*r)/(2*a))*(BesselI[n-1,a*r]*BesselK[n,a*r]-
BesselI[n,a*r]*BesselK[n-1,a*r])

• Perhaps you can shed some light on how you obtained the bottom expression? At least, that's not what I get from naive FunctionExpand and FullSimplify routines on the generic form of the integral. Aug 24, 2023 at 3:31
• @Aaron if the limit n-> integer evaluates to indeterminate one can try the series approach as I did above. Then I tried for different n, integer and non integer and checked the derivative symbolically and numerically. Aug 24, 2023 at 8:44
• the term -(n/(2*a^2)) may be omitted as it is just an integration constant Aug 24, 2023 at 12:04