Looking for a reference of integral involving product of four spherical harmonics We know
$$\int d \Omega Y_{l_1m_1}(\theta,\phi) Y_{l_2 m_2}(\theta,\phi)  Y_{l_3 m_3 } (\theta,\phi) = \sqrt{ \frac{ (2l_1 + 1)(2 l_2+1)(2l_3+1)}{4\pi} } \pmatrix{ l_1 l_2 l_3 \\ m_1 m_2 m_3 } \pmatrix{ l_1 l_2 l_3 \\ 0 0 0 } $$
Is there any reference for integral involving four spherical harmonics? I could try to work it out myself, but a cross check with literature maybe useful
 A: We can always reduce the product of two spherical harmonics to a sum over one spherical harmonic (DLMF 34.3.20),
$$
Y_{l_1m_1} \, Y_{l_2m_2}
= \sum_{lm} \sqrt{\tfrac{(2l_1+1)(2l_2+1)(2l+1)}{4\pi}} \, \bigl(\begin{smallmatrix}l_1&l_2&l\\0&0&0\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_1&l_2&l\\m_1&m_2&m\end{smallmatrix}\bigr) \, \overline{Y_{lm}} \\
= \sum_{lm} (-1)^m \, \sqrt{\tfrac{(2l_1+1)(2l_2+1)(2l+1)}{4\pi}} \, \bigl(\begin{smallmatrix}l_1&l_2&l\\0&0&0\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_1&l_2&l\\m_1&m_2&-m\end{smallmatrix}\bigr) \, Y_{lm}
$$
Hence for the integral over $n+1$ spherical harmonics:
$$
\int Y_{l_1m_1} \cdots Y_{l_nm_n} \, Y_{l_{n+1}m_{n+1}} \\
= \sum_{lm} (-1)^m \, \sqrt{\tfrac{(2l_n+1)(2l_{n+1}+1)(2l+1)}{4\pi}} \, \bigl(\begin{smallmatrix}l_n&l_{n+1}&l\\0&0&0\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_n&l_{n+1}&l\\m_n&m_{n+1}&-m\end{smallmatrix}\bigr) \int Y_{l_1m_1} \cdots Y_{l_{n-1}m_{n-1}} \, Y_{lm}
$$
We can build expressions recursively from there. For example, the next formula after $n = 3$:
$$
\int Y_{l_1m_1} \, Y_{l_2m_2} \, Y_{l_3m_3} \, Y_{l_4m_4} \\
= \sum_{lm} (-1)^m \, \sqrt{\tfrac{(2l_3+1)(2l_4+1)(2l+1)}{4\pi}} \, \bigl(\begin{smallmatrix}l_3&l_4&l\\0&0&0\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_3&l_4&l\\m_3&m_4&-m\end{smallmatrix}\bigr) \int Y_{l_1m_1} \, Y_{l_2m_2} \, Y_{lm} \\
= \tfrac{\sqrt{(2l_1+1)(2l_2+1)(2l_3+1)(2l_4+1)}}{4\pi} \sum_{lm} (-1)^m \, (2l+1) \, \bigl(\begin{smallmatrix}l_1&l_2&l\\0&0&0\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_3&l_4&l\\0&0&0\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_1&l_2&l\\m_1&m_2&m\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_3&l_4&l\\m_3&m_4&-m\end{smallmatrix}\bigr)
$$
The sum over $m$ picks out the condition that $m_1 + m_2 + m_3 + m_4 = 0$, and the integral vanishes otherwise. Assuming the condition is fulfilled:
$$
\int Y_{l_1m_1} \, Y_{l_2m_2} \, Y_{l_3m_3} \, Y_{l_4m_4} \\
= (-1)^{m_1+m_2} \, \tfrac{\sqrt{(2l_1+1)(2l_2+1)(2l_3+1)(2l_4+1)}}{4\pi} \sum_{l} (2l+1) \, \bigl(\begin{smallmatrix}l_1&l_2&l\\0&0&0\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_3&l_4&l\\0&0&0\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_1&l_2&l\\m_1&m_2&m_3+m_4\end{smallmatrix}\bigr) \, \bigl(\begin{smallmatrix}l_3&l_4&l\\m_3&m_4&m_1+m_2\end{smallmatrix}\bigr)
$$
I don't immediately see if the remaining sum reduces to a known expression, but the $3j$-symbols with $000$ at least can be reduced further, so perhaps there is something more to be done. I am marking this answer community wiki in the hope that someone might be able to improve it.
