I have the following identity,

$$\int_{-1}^{1} \! dz \, j_0\bigl(\sqrt{x^2 + y^2 - 2xy z}\bigr) \, P_n(z) = 2 \, j_n(x) \, j_n(y) \;,$$

where $x, y > 0$, $P_n$ is a Legendre polynomial, and $j_n$ is a spherical Bessel function. This follows from Gegenbauer's addition theorem (NIST 10.60.2),

$$ j_0\bigl(\sqrt{x^2 + y^2 - 2xy z}\bigr) = \sum_{l \ge 0} (2l+1) \, j_l(x) \, j_l(y) \, P_l(z) \;. $$

I am interested in the generalisation of the above to general Jacobi polynomials $P^{(\alpha, \beta)}_n$, $$ \int_{-1}^{1} \! dz \, j_0\bigl(\sqrt{x^2 + y^2 - 2xy z}\bigr) \, (1-z)^\alpha (1+z)^\beta \, P^{(\alpha,\beta)}_n(z) \;, $$ with $\alpha, \beta > -1$ and $x, y > 0$.

We have the following connection formula for Jacobi and Legendre polynomials, $$ P_l(z) = \sum_{k=0}^l \tfrac{\Gamma(k+\alpha+\beta+1) \, \Gamma(k+l+1)}{\Gamma(k+1) \, \Gamma(l-k+1) \, \Gamma(2k+\alpha+\beta+1)} \, _3F_2\bigl(\genfrac..{0pt}{1}{k-l, k+l+1, k+\alpha+1}{k+1, 2k+\alpha+\beta+2}; 1\bigr) \, P^{(\alpha,\beta)}_k(z) \;, $$ with $_3F_2$ the hypergeometric function. With this, we can rewrite the addition theorem for general Jacobi polynomials, $$ j_0\bigl(\sqrt{x^2 + y^2 - 2xyz}\bigr) = \sum_{l \ge 0} R^{(\alpha,\beta)}_l(x, y) \, P^{(\alpha,\beta)}_l(z) \;, $$ using the new coefficient functions $$ R_l^{(\alpha,\beta)}(x, y) = \sum_{k \ge l} (2k + 1) \, j_k(x) \, j_k(y) \, \tfrac{\Gamma(l+\alpha+\beta+1) \, \Gamma(k+l+1)}{\Gamma(l+1) \, \Gamma(k-l+1) \, \Gamma(2l+\alpha+\beta+1)} \, _3F_2\bigl(\genfrac..{0pt}{1}{l-k, l+k+1, l+\alpha+1}{l+1, 2l+\alpha+\beta+2}; 1\bigr) \;. $$

*Q:** Can the functions $R_l^{(\alpha,\beta)}(x, y)$ be simplified, even in special cases? Are there other cases besides $\alpha = \beta = 0$ where the result factors as $f(x)\,g(y)$?*

One special case might be $\beta = 0$, where the hypergeometric function can be replaced using Saalschütz's theorem.

*(This question was updated with results from a deleted answer.)*