Finding some orthogonality in a convolution-like integral over Legendre polynomials I encountered the following integral in my (physics) research, and I've yet to find an analytic solution:
$$I(n_1,n_2,n_3) = \int_{-1}^{1} d(\cos\theta_1) \int_{-1}^{1} d(\cos\theta_2) P_{n_1}(\cos\theta_1) P_{n_2}[\cos(\theta_1-\theta_2)] P_{n_3}(\cos\theta_2)$$
where $P_n(x)$ is the nth Legendre polynomial. For reference, the orthogonal polynomial relation for Legendre polynomials is:
$$\int_{-1}^{1} d(\cos\theta) P_{n_1}(\cos\theta) P_{n_2}(\cos\theta) = \int_{0}^{\pi} d\theta (\sin\theta) P_{n_1}(\cos\theta) P_{n_2}(\cos\theta) = \frac{2}{2n_1+1}\delta_{n_1,n_2}$$
Mathematica does not have much trouble in evaluating the integral exactly (Integrate, not NIntegrate) for a given $\{n_1, n_2, n_3\}$.
k[n1_,n2_,n3_]:=Integrate[Sin[x]Sin[y]LegendreP[n1,Cos[x]]LegendreP[n2,Cos[x-y]]LegendreP[n3,Cos[y]],{x,0,Pi},{y,0,Pi}]

I can see already that $I(n_1, n_2, n_3)$ is not diagonal. Even the diagonal terms contain factors which I've yet to figure out, here's what I see so far:
$$I(n,n,n) = 2^n \left[\frac{2}{2n+1}\right]^2 \times \left\{1,\frac{1}{2},\frac{1}{3},\frac{1}{5},\frac {4}{35},\frac{4}{63},\frac{8}{231},\frac{8}{429}, \frac{64}{6435}, ... \right\}$$
where the bracketed terms correspond to $n=0,1,2,...$
Many of the off-diagonal terms appear to contain a $\pi^2$ term and some power of 2.
Do you have any suggestions?  Useful variable transformations, pattern recognition?  Thanks in advance.
 A: I guess this is one of these cases where the spherical harmonics addition theorem
$$ P_l(\cos \gamma)  = \frac{4\pi}{2l+1} \sum_{m=-l}^l Y^*_{lm}(\theta_1,\phi_1) Y_{lm}(\theta_2,\phi_2)$$
 comes in handy; here, $\gamma$ is the angle between two vectors having spherical coordinates $(\theta_1,\phi_1)$ and $(\theta_2,\phi_2)$ and $Y_{lm}$ are the spherical harmonics.
In your case, we choose $\phi_1=\phi_2=0$ and the fact that
$$Y_{lm}(\theta,0)= \sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}} P_l^m(\cos\theta) $$
with $P_l^m$ the associated Legendre polynomials. The addition theorem thus reads
$$P_l [\cos(\theta_1 - \theta_2)]= \sum_{m=-l}^l \frac{(l-m)!}{(l+m)!}P_l^m(\cos\theta_1)P_l^m(\cos\theta_2).$$
Plugging the addition theorem into your integral yields,
$$I(n_1,n_2,n_3) = \int_{-1}^{1}\! dx_1 \int_{-1}^{1}\! dx_2\, P_{n_1}(x_1) \left[\sum_{m=-n_2}^{n_2} \frac{(n_2-m)!}{(n_2+m)!}P_{n_2}^m(x_1)P_{n_2}^m(x_2) \right] P_{n_3}(x_2).$$
In a last step, we need to know the integral
$$p(k,l,m)=\int_{-1}^1 \!dx \,P_k(x) P_l^m(x)$$ in terms of which we can find the expression
$$I(n_1,n_2,n_3) = \sum_{m=-n_2}^{n_2} \frac{(n_2-m)!}{(n_2+m)!} p(n_1,n_2,m) p(n_3,n_2,m).$$
It turns out that the general expression for $p(k,l,m)$ is rather complicated.
However, for $k=l$, we have the simple expression $$p(l,l,m)= \begin{cases}
\frac{2(-1)^{m/2}  l! (2l)!}{(1+2l)! (l-m)!} & m\text{ even},\\
0& m\text{ odd}.\end{cases} $$
For the diagonal elements, we obtain
$$I(n,n,n)= \sum_{\substack{m=-n;\\m\text{ even}}}^n \frac{1}{(n+m)!}\frac{4 \,n!^2\, (2 n)!^2}{(2 n+1)!^2\, (n-m)!} .
$$
As the expression for the diagonal element is already rather complicated, and I am afraid the off diagonal elements will be even more cumbersome. However, using  the paper linked above, you should be able to obtain an explicit expression.
A: This solution is a continuation of Fabian's answer.  Thanks for pointing me in the right direction!
Starting from the original definition of $I(n_1,n_2,n_3)$:
$$I(n_1,n_2,n_3) = \int_{-1}^{1} d(\cos\theta_1) \int_{-1}^{1} d(\cos\theta_2) P_{n_1}(\cos\theta_1) P_{n_2}[\cos(\theta_1-\theta_2)] P_{n_3}(\cos\theta_2)$$
we can expand the middle Legendre polynomial in terms of associated legendre polynomials as follows:
$$P_l [\cos(\theta_1 - \theta_2)]= \sum_{m=-l}^l \frac{(l-m)!}{(l+m)!}P_l^m(\cos\theta_1)P_l^m(\cos\theta_2)$$
Then $I(n_1,n_2,n_3)$ is written:
$$I(n_1,n_2,n_3) = \int_{-1}^{1} dx_1 \int_{-1}^{1} dx_2\, P_{n_1}(x_1) \left[\sum_{m=-n_2}^{n_2} \frac{(n_2-m)!}{(n_2+m)!}P_{n_2}^m(x_1)P_{n_2}^m(x_2) \right] P_{n_3}(x_2)$$
Defining the overlap of a Legendre polynomial with an associated Legendre polynomial as
$$X(n,n_2,m_2)=\int_{-1}^1 dx\, P_{n}(x) P_{n_2}^{m_2}(x)$$
so that
$$I(n_1,n_2,n_3) = \sum_{m_2=-n_2}^{n_2} \frac{(n_2-m)!}{(n_2+m)!} X(n_1,n_2,m_2) X(n_3,n_2,m_2)$$
A single sum expression has been derived for the overlap integral of two associated Legendre polynomials in J. Phys. A: Math. Gen. 32 (1999) 2601-2603 (along with a small correction in a comment paper: J. Phys. A: Math. Gen. 35 (2002) 4187-4188). The solution in that paper is for
$$ Z(n_1,m_1,n_2,m_2) = \int_{-1}^1 dx\, P_{n_1}^{m_1}(x) P_{n_2}^{m_2}(x) $$
but we only need $X(n,n_2,m_2)=Z(n,0,n_2,m_2)$. This simplifies the resulting expressions slightly. When $m_2=0$ we have the overlap of two Legendre polynomials, so
$$X(n,n_2,0) = \frac{2}{2n+1} \delta_{n,n_2}$$
When $m_2 \neq 0$,
$$X(n,n_2,m_2 \neq 0) = A(n_2,m_2) \sum_{l=l_{\text{min}}}^{l_{\text{max}}} D(l,m_2)(2l+1) 
\begin{pmatrix}
n & n_2 & l \\
0 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
n & n_2 & l \\
0 & m_2 & -m_2
\end{pmatrix}
$$
where
$
\begin{pmatrix}
j_1 & j_2 & j_3 \\
m_1 & m_2 & m_3
\end{pmatrix} \text{ is the Wigner 3-j symbol}
$
and
\begin{aligned}
 A(n,m) =& (-1)^{\tau(m)} |m|\, 2^{|m|-2} \sqrt{\frac{(n+m)!}{(n-m)!}} \\
 \tau(m) =& \begin{cases}0 & \text{if } m \geq 0 \\ m & \text{if } m < 0\end{cases} \\
 D(l,m) =& \begin{cases}
  \left(1+(-1)^{l+|m|}\right) \sqrt{\frac{(l-|m|)!}{(l+|m|)!}} \frac{\Gamma(l/2)\Gamma[(l+|m|+1)/2]}{\Gamma[(l+3)/2][(l-|m|)/2]!} & \text{if } l+m \text{ is even} \\
  0 & \text{if } l+m \text{ is odd}
 \end{cases} \\
 l_{min} =& \begin{cases}|n-n_2| & \text{if } |n-n_2| \geq |m_2| \\ |m_2| & \text{if } |n-n_2| < |m_2|\end{cases} \\
 l_{max} =& n + n_2
\end{aligned}
This completes the general solution for $I(n_1,n_2,n_3)$. As pointed out by Fabian, the diagonal elements have a considerably simpler form:
$$X(n,n,m) = \begin{cases}
\frac{2(-1)^{m/2}  n! (2n)!}{(2n+1)! (n-m)!} & m\text{ even}\\
0& m\text{ odd}\end{cases}$$
so that $I(n,n,n)$ can be quickly calculated from
$$I(n,n,n) = \sum_{\substack{m=-n\\m\text{ even}}}^{n} \frac{1}{(n+m)!}\frac{4[n!]^2[(2n)!]^2}{[(2n+1)!]^2(n-m)!}$$
