By experimenting in Mathematica, I have found the following expression for the integral: $$ \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=\\ aj_{n'}(a)\sum_{k=0}^{\min(n,n')-m}c_{n,n'}^{m,k}\frac{h_{n+n'-m-k}^{(1)}(b)}{b^{k+m-1}} $$
where $$ c_{n,n'}^{m,k}=\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)!} $$ and $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')$.
An equivalent, potentially tidier, formulation is: $$ \lim_{\delta\to 0}\frac{1}{\delta^{n'}}\int_{-1}^{1}h_{n}^{(1)}(b(1+\delta\sigma))P_{n}^{m}\left(1-\frac{\delta^{2}(1-\sigma^{2})}{2(1+\delta\sigma)}\right)P_{n'}^{m}\left(\sigma-\frac{\delta(1-\sigma^{2})}{2}\right)(1+\delta\sigma)d\sigma=\frac{2^{n'}(n')!}{(2n'+1)!}\sum_{k=0}^{\min(n,n')-m}c_{n,n'}^{m,k}\frac{h_{n+n'-m-k}^{(1)}(b)}{b^{k+m-n'}} $$ The result can be seen by running the code:
j[n_?IntegerQ, x_] :=
Simplify[(-x)^n Nest[D[#, z]/z &, Sin[z]/z, n] /. z -> x];
h[n_?IntegerQ, x_] :=
Simplify[-I (-x)^n Nest[D[#, z]/z &, Exp[I z]/z, n] /. z -> x];
m = 0;
n = 2;
np = 1;
int = Assuming[b > a > 0,
FullSimplify[
Integrate[\[Sigma] SphericalHankelH1[n, \[Sigma]] LegendreP[n,
m, (\[Sigma]^2 - a^2 + b^2)/(2 b \[Sigma])] LegendreP[np,
m, (\[Sigma]^2 - a^2 - b^2)/(2 a b)], {\[Sigma], b - a, b + a}]]]
SolveAlways[
int Exp[-I b]/(a j[np, a]) ==
Sum[a1[k] h[n + np - k - m, b] Exp[-I b]/b^(k + m - 1), {k, 0,
Min[n, np] - m}], b][[1, All, 2]]
where any combination of $m,n,n'$ such that $m\leq \min(n,n')$ can be chosen. The first output gives the full expression for the integral; the second output shows the coefficients, $c_{n,n'}^{m,k}$.
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. In the second formulation, I've tried expanding out the Legendre functions in powers of $\delta$, but, again, it gets very messy. An inductive proof would be elegant, but I don't see how to do this either.
