A problematic integral in calculating the entanglement entropy in 1+1 D free massive bosonic field theory I encountered a curious integration identity when I was reading the paper by Pasquale Calabrese and John Cardy on the entanglement entropy of 1+1D quantum field theory (arXiv). The identity is given below:
$$
\int^\infty_0 x I_{\alpha}(x)K_{\alpha}(x) dx \,\, ``=" \frac{\alpha}{2}.
$$
Here $I_\alpha(x)$ and $K_{\alpha}(x)$ are the modified Bessel functions of the fist and second kind. The above identity is implicitly used in Eq.4.19 and Eq.4.22 of the said paper. However, the above identity looks false to me as the left hand side is apparently divergent. Recall that we have the following asymptotic expansion:
$$
I_\alpha (x) \sim \frac{e^{x}}{\sqrt{2\pi x}};\\
K_\alpha (x) \sim \sqrt{\frac{\pi}{2 x}}e^{-x},
$$
from which we obtain
$$
xI_\alpha(x)K_\alpha(x) \sim \frac{1}{2}.
$$
I believe that some sort of regularization has been done on the integral or I simply miss something. I hope someone here could help me understanding this identity.
EDIT: I think the following might be a useful hint.
The above integral rises from the expression of $\int^\infty_0 rG(r,r)dr$, where $r$ is the radial coordinate and $G$ is the usual Green's function (see Eq.4.18). Now let us consider $\int^\infty_0 G((1+\epsilon)r,r)rdr$, motiviated by the short distance expansion. The corresponding integral is:
$$
\int^\infty_0 x I_{\alpha}((1+\epsilon)x)K_{\alpha}(x) dx = \frac{(1+\epsilon)^{-\alpha}} {(1+\epsilon)^2-1} = \frac{1}{2\epsilon}-\frac{1}{2}(\alpha+\frac{1}{2})+O(\epsilon).
$$
The above identity can be found in Gradshteyn & Ryzhik's table, 7th Edition, entries (6.521.2), (6.521.8), and (6.521.9).
Now we take the limit $\epsilon\to 0$. The first term of the above is divergent, which could be interpreted as the usual UV divergence part. The second term, however, is finite. If somehow we can argue that it is meaningful to drop the first term and take the $\alpha$-dependent as the "physical" answer of the integral, then we are done.
 A: Even though I haven't been able to come up with a satisfactory interpretation for the integral identity I asked about, I did find a way to circumvent the problem and derived the final result in the said paper. After waiting for a week without receiving any definitive solution, I post my calculation below for it is not entirely irrelevant to my original question and it might be useful for fellow students.
Define
$$
I_n = \int^\infty_0 r G(r,r) dr = \sum_{k\ge0} d_k \int^\infty_0 rI_{k/n}(mr) K_{k/n}(mr)dr.
$$
Here $d_0=1/2$ and $d_k=1$ for $k>1$. The goal is to compute $I_n-nI_1$.
As I mentioned earlier, the above integral over radius $r$ is divergent, not mentioning the divergent summation over $k$. Here we consider the following modified sum,
$$
I_n(\epsilon) = \int^\infty_0 r G[(1+\epsilon)r,r] dr = \sum_{k\ge0} d_k \int^\infty_0 rI_{k/n}(mr) K_{k/n}[(1+\epsilon)mr]dr.
$$
We will keep $\epsilon>0$ and take the limit $\epsilon\to 0$ later. The integral over $r$ now can be done and the result is
$$
I_n(\epsilon)=\frac{1}{m^2[(1+\epsilon)^2-1]}\sum_{k\ge0}d_k(1+\epsilon)^{-k/n}=\frac{1}{m^2(2\epsilon+\epsilon^2)}\left(\frac{1}{1-(1+\epsilon)^{-1/n}}-\frac{1}{2}\right),
$$
from which we find
$$
I_n(\epsilon) = \frac{1}{m^2(2\epsilon+\epsilon^2)}\left(\frac{n}{\epsilon}+\frac{n}{2}+\frac{n\epsilon}{12}(\frac{1}{n^2}-1)+O(\epsilon^2)\right).
$$
Likewise,
$$
nI_1(\epsilon) = \frac{1}{m^2(2\epsilon+\epsilon^2)}(\frac{n}{\epsilon}+\frac{n}{2}).
$$
Thus,
$$
\lim_{\epsilon\to0}[I_n(\epsilon)-nI_1(\epsilon)]=
\lim_{\epsilon\to0}\frac{1}{2m^2\epsilon}\times\frac{n\epsilon}{12}(\frac{1}{n^2}-1)=\frac{1}{24m^2}(\frac{1}{n}-n),
$$
which is the desired result. Comparing to the original derivation in the paper, the above avoids the dubious integral identity and the Euler-Maclaurin summation formula altogether.
A: Have you tried with Mathematica?
I plugged in the original expression
$$\int\limits_0^\infty I_{k/n}(mr)K_{k/n}(mr)rdr$$
into Mathematica and got (with assumtions that $k,\ n\in\mathbb{Z}$)
 Integrate[BesselI[k/n, m*r]*BesselK[k/n, m*r]*r, {r, 0, \[Infinity]}, 
  Assumptions -> 
   k \[Element] Integers && n \[Element] Integers] // FullSimplify

ConditionalExpression[-(k/(2 m^2 n)), n (k + n) > 0 && Re[m] > 0]

$$\text{ConditionalExpression}\left[-\frac{k}{2 m^2 n},n (k+n)>0\land
   \Re(m)>0\right]$$
And with $\alpha=k/n$ you get the correct result (up to the minus sign, but that could be 'eaten away' by the minus in front of the derivative of the logarithm of partition function).
Look at comments for more explanations about $m$ and $r$.
