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!