From numerical experiments in Mathematica, I have found the following expression for the integral: $$ \int_{-1}^{1}h_{n}^{(1)}\left(\sqrt{a^{2}+b^{2}+2ab\tau}\right)P_{n}^{m}\left(\frac{a\tau+b}{\sqrt{a^{2}+b^{2}+2ab\tau}}\right)P_{n'}^{m}(\tau)d\tau=\sum_{k=0}^{\min(n,n')-m}\frac{(-1)^{k+m+n'}(n+m)!(n'+m)!(2k+2m)!}{2^{k+m-1}k!(k+m)!(k+2m)!(n-m-k)!(n'-m-k)!b^{k+m}}j_{n'}(a)h_{n+n'-m-k}^{(1)}(b) $$ where $j_{n}$ are the spherical Bessel functions of the first kind, $h_{n}^{(1)}$ are the spherical Hankel functions of the first kind, $P_{n}^{m}$ are the associated Legendre polynomials, $b>a>0$ and $m,n,n'$ are non-negative integers such that $m\leq \min(n,n')$. This can be seen by running the code:
b = 5;
nmax = 5;
Flatten[Table[Table[{m, n, n1,
ListPlot[
Table[{a,
Log10@Abs[NIntegrate[
SphericalHankelH1[n,
Sqrt[a^2 + b^2 + 2 a b \[Tau]]] LegendreP[n,
m, (a \[Tau] + b)/
Sqrt[a^2 + b^2 + 2 a b \[Tau]]] LegendreP[n1,
m, \[Tau]], {\[Tau], -1,
1}] - (-1)^(m + n1) (n + m)! (n1 + m)! SphericalBesselJ[n1,
a] Sum[(-1)^
k (2 k + 2 m)! SphericalHankelH1[n + n1 - m - k,
b]/(2^(k + m - 1) k! (k + m)! (k + 2 m)! (n - m -
k)! (n1 - m - k)! b^(k + m)), {k, 0,
Min[n, n1] - m}]]}, {a, 0, 4, 0.5}], Joined -> True]}, {n,
m, nmax}, {n1, m, nmax}], {m, 0, nmax}],2] //Quiet
nmax =.
b =.
where it is seen that the differences between the left and right hand sides are numerically very small.
The integral arises as the result of considering the scattering from multiple spheres. Essentially modes on a sphere can be represented as 'vector spherical harmonics'. These can be written in terms of spherical Bessel functions and Legendre polynomials. What the integral represents is the interaction between charge distributions on different spheres described by these modes.
I would like to be able to formally show this but (unsurprisingly) I'm having a few difficulties and I haven't been able to find what I'm looking for in places like DLMF and Gradshteyn and Ryzhik. I've tried using various definitions and expansions for the Legendre and spherical Hankel functions but it all just gets too messy. Also, I'm keen to know if the right hand side can be simplified to something a bit more manageable.
Progress update: I feel the substitution $\sigma=\sqrt{a^{2}+b^{2}+2ab\tau}$ may help. This converts the integral to: $$ \frac{1}{ab}\int_{b-a}^{b+a}\sigma h_{n}^{(1)}(\sigma)P_{n}^{m}\left(\frac{\sigma^{2}-a^{2}+b^{2}}{2b\sigma}\right)P_{n'}^{m}\left(\frac{\sigma^{2}-a^{2}-b^{2}}{2ab}\right)d\sigma $$ At least this gets rid of the square root. I'm not sure how to proceed next, though.
For full disclosure, I have already posted this question on StackExchange at https://math.stackexchange.com/questions/4257140/how-to-calculate-the-integral-of-a-product-of-a-spherical-hankel-function-with-a, but without success (despite multiple bounties placed and, as of writing, 17 upvotes and nearly 600 views).
