Bug in Integrate involving hypergeometric functions?

Question

I am working with v. 13.0.0.0 for Mac. I am attempting to compute the following integral with Mathematica for integer values of L>=2 and with a>0 and r>0:

I1[L_,a_,r_]=Integrate[(1 - x^2)^
  L (1/(r + I a x)^(L + 5) (2 I a x - r (L + 2)) + (-1)^(L + 1)
      1/(r - I a x)^(L + 5) (2 I a x + r (L + 2))), {x, -1, 1}, 
 Assumptions -> r > 0 && a > 0 && L >= 2 && L \[Element] Integers];

Mathematica does give an expression written in terms of Hypergeometric functions. However, if I for instance set L=2 first, and then compute the same integral, I am not finding compatible results. For instance, if I execute the following command:

L=2;
I2=Integrate[(1 - x^2)^
  L (1/(r + I a x)^(L + 5) (2 I a x - r (L + 2)) + (-1)^(L + 1)
      1/(r - I a x)^(L + 5) (2 I a x + r (L + 2))), {x, -1, 1}, 
 Assumptions -> r > 0 && a > 0 && L >= 2 && L \[Element] Integers];

I am getting

I2=(64 (a^2 - 2 r^2))/(15 (a^2 + r^2)^4).

which is the correct result. However, from I1 above I am getting

I1[2,a,r]=-((8 (5 a^3 + 20 a^2 r + 29 a r^2 + 16 r^3))/(15 r^5 (a + r)^4))

What is going on here?

Answer 1

Reproducing OP's result

The notebook can be found on my github. My version is "13.1.0 for Mac OS X ARM (64-bit) (June 16, 2022)" and I can reproduce the OP's result by

Clear[f,I1]
f[x_,{r_,a_,L_}]:=(1-x^2)^L ((r+I a x)^(-5-L) (-((2+L) r)+2 I a x)+(-1)^(1+L) (r-I a x)^(-5-L) ((2+L) r+2 I a x))

I1[r_,a_,L_]:=Integrate[f[x,{r,a,L}],{x,-1,1},Assumptions->L>-1&&r>0&&a>0]

I1[r,a,2]
(*(64 (a^2-2 r^2))/(15 (a^2+r^2)^4)*)

I1[r,a,L]/.{L->2}//Simplify
(*-((8 (5 a^3+20 a^2 r+29 a r^2+16 r^3))/(15 r^5 (a+r)^4))*)

Actually if we modify the assumptions of I1 without $a>0$

I11[r_,a_,L_]:=Integrate[f[x,{r,a,L}],{x,-1,1},Assumptions->L>-1&&r>0]

Mathematica only returns results with $\Im a \neq 0$

I11[r,a,L]/.{L->2}//Simplify

Indeed there is an inconsistency of Mathematica.

---

A workaround

We need to absorb the factor I into a by

Clear[g,I2];
g[x_,{r_,a1_,L_}]:=Evaluate[f[x,{r,a,L}]/.a->-I a1]

I2[r_,a1_,L_]:=Integrate[g[x,{r,a1,L}],{x,-1,1},Assumptions->L>-1&&r>0]

Now the order of integration and limit is commutative:

I2[r,a1,2]/.a1->I a

I2[r,a1,L]
I2[r,a1,L]/.{L->2,a1->I a}//Simplify

---

Some explanations

As a multi-variable complex function, the hypergeometric ${}_p F_{q}(\{a\},\{b\},z)$ has branch cut at $z>1$ for generic values of $\{a\},\{b\}$, and has simple poles with respect to ${b}$. It's often hard to implement the analytic structures in symbolic computation.

The branch cut can be detected in the integral expression of the hypergeometric function as follows. For simplicity we only consider the ${}_2 F_{1}$ from dlmf

\begin{equation} \mathbf{F}(a, b ; c ; z)=\frac{1}{\Gamma(b) \Gamma(c-b)} \int_0^1 \frac{t^{b-1}(1-t)^{c-b-1}}{(1-z t)^a} d t, \end{equation}

The singularities at $t=0,1$ are okay for suitable choices of $b,c$. The only divergence happens at $t=1/z$. To see the branch cut we take $a=1$ and consider the distribution

\begin{equation} \int^{1}_{0}dt f(t)\frac{1}{t-1/z} \end{equation}

If $1/z\in (0,1)$ this is divergent except for test functions satisfying $f(1/z)=0$. But if $1/z$ has a small nonvanishing imaginary part $1/z -I \epsilon \in (0,1)$ it is convergent due to

\begin{equation} \lim _{\epsilon \rightarrow 0^{+}} \frac{1}{x \pm i \epsilon}=\mp i \pi \delta(x)+\mathcal{P}\left(\frac{1}{x}\right) \end{equation}

(This is called the regularization of distributions.)

The ambiguity of the sign of $\epsilon$ corresponds to the direction of $z$ approaching to the real interval, either from above or below, which means we have a branch cut of $f(z)$ at $z>1$.

Back to this question, I believe that starting from real axis and then analytic continue it onto the imaginary axis is more safe.