Gegenbauer's addition theorem for Jacobi polynomials 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.)
 A: Too long for a comment: Since the Legendre polynomials define the zonal harmonics on the sphere, a similar formula should hold at least for the Jacobi polynomials which correspond to harmonics for the other two-point homogeneous spaces. The spherical Bessel function may not appear in this case, but instead some other function. If you are particulary interested in an equation which is exactly in the form as the spherical $d=3$ case, the below will likely not be very relevant.
For compact two-point homogeneous spaces generally, there exists a decomposition of $L^2(\Omega)$, the space of complex-valued square-integrable functions on $ \Omega $, into irreducible representations:
$L^2(\Omega)=\bigoplus\limits_{n\geq 0}V_n ,$ 
 where spaces $ V_n $ are finite-dimensional and invariant under the isometries of $ \Omega $. Moreover, they can  be chosen as the eigenspaces  of the Laplace--Beltrami operator on $\Omega$  corresponding to the $n$-th eigenvalue in the increasing order. 
Let $Y_{n,k}$, $ k =1,\ldots,\dim V_n $, be an orthonormal basis in $V_n$. Because of the invariance of $ V_n $ and due to the two-point homogeneity of $\Omega$, the reproducing kernel for $ V_n $ only depends on the distance $\vartheta(x,y)$  between points . Furthermore, as a function of $\tau(x,y) := \cos \vartheta(x,y)$  the reproducing kernel is a polynomial $C_n$ of degree $n$, which satisfies
$$C_n(\tau(x,y))= \frac1{\dim V_n} \sum\limits_{k=1}^{\dim V_n} Y_{n,k}(x) \overline{Y_{n,k}(y)}.$$
The formula above is known as the addition formula on these spaces.
The Jacobi polynomials $C_n$ given by the addition formula satisfy $ C_n(1) = 1 $ and are orthogonal with respect to the probability measure
$$d\nu^{(\alpha, \beta)}  =  \frac1{\gamma_{\alpha,\beta}} (1-t)^{\alpha} (1+t)^{\beta} dt,$$
where $\alpha = (d-1) \dim_{\mathbb R} (\mathbb F)/{2}-1$ and
$$  \beta=\alpha, \text{ if } \Omega = \mathbb{S}^{d-1} \\ $$
$$\beta={\dim_{\mathbb{R}}(\mathbb{F})}/{2} -1, \text{ if } \Omega = \mathbb{FP}^{d-1}, $$
 and the normalization factor is given by
$ \gamma_{\alpha,\beta} = 2^{\alpha +\beta +1}B(\alpha+1, \beta+1),$
 where $B$ is the beta function.
