Computing the integral $\int_{-1}^1 dx \, |x| J_0(\alpha \sqrt{1 - x^2}) P_\ell(x)$ I would like to compute the following integral:
$$
I_\ell(\alpha) := \int_{-1}^1 dx \, |x| J_0(\alpha \sqrt{1 - x^2}) P_\ell(x)
\tag{1}
\label{1}
$$
where $\alpha \geq 0$, $J_0$ is the zeroth-order Bessel function of the first kind, $P_\ell(x)$ is the Legendre polynomial of order $\ell$, and $\ell$ is an arbitrary positive integer or zero.
Since the integrand is odd if $\ell$ is odd, we have that $I_\ell(\alpha) = 0, \ell \text{ odd}$, so we just need to care about even $\ell$s.
Mathematica reports remarkably simple results for some actual values of $\ell$:
$$
I_0(\alpha) = \frac{2 J_1(\alpha )}{\alpha }\\
I_2(\alpha) = \frac{6 J_2(\alpha )-\alpha  J_1(\alpha )}{\alpha ^2}\\
I_4(\alpha) = \frac{3 \alpha ^2 J_1(\alpha )-60 \alpha  J_2(\alpha )+280 J_3(\alpha )}{4 \alpha ^3}
$$
This seems to suggest that we have something along the lines of (purely heuristically, not necessarily true):
$$
I_\ell(\alpha)
= \sum\limits_k a_k \alpha^{b_k} J_k(\alpha)
$$
where $b_k$ seem to be integers.
Now, one idea I had in mind was to use the expansion (DLMF 10.60, written out in a more suitable form):
$$
J_0\left(\alpha\sqrt{1 - x^2}\right)=\sum_{n=0}^\infty (4n+1) \frac{(2n)!}{2^{2n}(n!)^2} j_{2n}(\alpha) P_{2n} (x)
$$
along with the following identities (see here and here):
$$
P_k P_\ell
= \sum\limits_{m=|k - \ell|}^{k + \ell}
\begin{pmatrix}k & \ell & m\\ 0 & 0 & 0\end{pmatrix}^2
(2m + 1) P_m \\
|x| = \begin{cases} -P_1(x),\quad x \leq 0\\ P_1(x),\quad x > 0 \end{cases}
\\
\int_0^1 dx\; P_m P_n
= \begin{cases}
\frac{1}{2n + 1}, & m=n\\
0, & m \neq n,m,n \text{ both even or odd}\\
f_{m,n}, & m \text{ even},n\text{ odd}\\
f_{n,m} ,& m \text{ odd},n\text{ even}
\end{cases}
$$
where I'll call $g(m,n) \equiv \int_0^1 dx\; P_m P_n$ for brevity, and:
$$
f_{m,n}
\equiv
\frac{(-1)^{(m+n+1)/2}m!n!}{2^{m+n-1} (m - n) (m + n + 1)
\big[\big(\frac{1}{2}m\big)!\big]^2 \big\{\big[\frac{1}{2}(n - 1)\big]!\big\}^2 }
$$
We can rewrite:
$$
\int_{-1}^1
d\mu\;
|\mu|
P_{2n} P_\ell
=
[(-1)^\ell + 1]
\int_0^1
d\mu\;
P_1 P_{2n} P_\ell
$$
and likewise:
\begin{align*}
\int_0^1
d\mu\;
P_1 P_{2n} P_\ell
&=
\sum\limits_{m=|2n - \ell|}^{2n + \ell}
\begin{pmatrix}
2n & \ell & m\\
0 & 0 & 0
\end{pmatrix}^2
(2m + 1)
\int_0^1
d\mu\;
P_1 P_m\\
&=
\sum\limits_{m=|2n - \ell|}^{2n + \ell}
\begin{pmatrix}
2n & \ell & m\\
0 & 0 & 0
\end{pmatrix}^2
(2m + 1)
g(1, m)
\end{align*}
so that we have:
\begin{align}
I_\ell(\alpha)
&=
\sum_{n=0}^\infty
\sum\limits_{m=|2n - \ell|}^{2n + \ell}
(4n+1) \frac{(2n)!}{2^{2n}(n!)^2} j_{2n}\left(\alpha\right)
[(-1)^\ell + 1]
\begin{pmatrix}
2n & \ell & m\\
0 & 0 & 0
\end{pmatrix}^2
(2m + 1) g(1, m)
\tag{2}
\label{2}
\end{align}
This is where I kind of got stuck, as I have no idea how to evaluate the double sum.
An alternate method would be to use the expansion from here:
$$
J_0(\alpha \sqrt{1 - x^2})
= e^{-\alpha x} \sum\limits_{n=0}^\infty
\frac{P_n(x)}{n!}\alpha^n
$$
but then I end up with integrals of the form:
$$
\int_{-1}^1 dx\; |x| P_\ell (x) P_n(x) e^{-\alpha x}
$$
which seems even more challenging to evaluate.
One idea for this one would be to expand $e^{-\alpha x} = \sum_k \frac{1}{k!} (-1)^k \alpha^k x^k$, and then rewrite $x^k$ as a linear combination of Legendre polynomials, but this again yields an integral over three Legendre polynomials, so I'd probably just obtain eq. \ref{2} in a more roundabout way.
Any hints would be appreciated!
 A: Thanks to the comment by Johannes, the solution can indeed be obtained by using the following identities:
\begin{equation}
P_\ell(z)
=
\frac{1}{2^\ell}
\sum\limits_{k=0}^{\left\lfloor \frac{\ell}{2}\right\rfloor}
(-1)^k
\begin{pmatrix} \ell \\ k \end{pmatrix}
\begin{pmatrix} 2\ell - 2k \\ \ell \end{pmatrix}
z^{\ell - 2k}
\tag{3}
\label{3}
\end{equation}
and:
\begin{equation}
\int_{0}^{\frac{1}{2}\pi}J_{\mu}\left(z\sin\theta\right)(\sin\theta)^{\mu+1}(%
\cos\theta)^{2\nu+1}\mathrm{d}\theta=2^{\nu}\Gamma\left(\nu+1\right)z^{-\nu-1}%
J_{\mu+\nu+1}\left(z\right)
\tag{4}
\label{4}
\end{equation}
Transforming the integral yields:
\begin{align}
I_\ell(\alpha)
&=
\int_{-1}^1 dx\, |x|\, J_0(\alpha\sqrt{1 - x^2}) P_\ell(x)\\
&=
[(-1)^\ell + 1]
\int_0^1 dx\, x\, J_0(\alpha \sqrt{1 - x^2}) P_\ell(x)\\
&=
|\mathrm{substitution}\;x = \cos \phi|\\
&=
[(-1)^\ell + 1]
\int_0^\frac{\pi}{2} d\phi\, \sin \phi\, \cos \phi\, P_\ell(\cos \phi)\, J_0 (\alpha \sin \phi)\\
&=
|\mathrm{expansion\;of}\;P_\ell|\\
&=
[(-1)^\ell + 1]
\sum\limits_{k=0}^{\left\lfloor \frac{\ell}{2}\right\rfloor}
(-1)^k
\begin{pmatrix} \ell \\ k \end{pmatrix}
\begin{pmatrix} 2\ell - 2k \\ \ell \end{pmatrix}
\int_0^\frac{\pi}{2} d\phi\, \sin \phi\, \cos \phi\, (\cos \phi)^{\ell - 2k} J_0 (\alpha \sin \phi)\\
&=
[(-1)^\ell + 1]
\sum\limits_{k=0}^{\left\lfloor \frac{\ell}{2}\right\rfloor}
(-1)^k
\begin{pmatrix} \ell \\ k \end{pmatrix}
\begin{pmatrix} 2\ell - 2k \\ \ell \end{pmatrix}
\int_0^\frac{\pi}{2} d\phi\, \sin \phi\, (\cos \phi)^{\ell - 2k + 1} J_0 (\alpha \sin \phi)
\end{align}
The integral in the above sum has the same form as the Bessel identity \ref{4}, with $\mu = 0$ and $\nu = \ell / 2 - k$, so that the final result is:
\begin{equation}
\boxed{
I_\ell(\alpha)
=
\frac{
[(-1)^\ell + 1]
}
{
2^\frac{\ell}{2}
}
\sum\limits_{k = 0}^{\left \lfloor \frac{\ell}{2} \right \rfloor}
\frac{(-1)^k}{2^k}
\begin{pmatrix} \ell \\ k \end{pmatrix}
\begin{pmatrix} 2\ell - 2k \\ \ell \end{pmatrix}
\Gamma\left[\frac{\ell}{2} - k + 1\right]
\frac{
J_{\frac{\ell}{2} - k + 1} (\alpha)
}
{
\alpha^{\frac{\ell}{2} - k + 1}
}
}
\tag{5}
\label{5}
\end{equation}
