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)
= sumlimits_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_0left(alphasqrt{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
= sumlimits_{m=|k – ell|}^{k + ell}
begin{pmatrix}k & ell & m\ 0 & 0 & 0end{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},ntext{ odd}\
f_{n,m} ,& m text{ odd},ntext{ 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}mbig)!big]^2 big{big[frac{1}{2}(n – 1)big]!big}^2 }
$$

We can rewrite:

$$
int_{-1}^1
dmu;
|mu|
P_{2n} P_ell
=
[(-1)^ell + 1]
int_0^1
dmu;
P_1 P_{2n} P_ell
$$

and likewise:

begin{align*}
int_0^1
dmu;
P_1 P_{2n} P_ell
&=
sumlimits_{m=|2n – ell|}^{2n + ell}
begin{pmatrix}
2n & ell & m\
0 & 0 & 0
end{pmatrix}^2
(2m + 1)
int_0^1
dmu;
P_1 P_m\
&=
sumlimits_{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
sumlimits_{m=|2n – ell|}^{2n + ell}
(4n+1) frac{(2n)!}{2^{2n}(n!)^2} j_{2n}left(alpharight)
[(-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} sumlimits_{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!