There are several important and useful formulas for the integration of spherical harmonic. We simplify the notation,

$\displaystyle \int_0^{\pi} \sin(\theta)d\theta\int_0^{2\pi}d\phi = \int d\Omega$

The first one is the average of spherical harmonic.

$\displaystyle \int Y_{lm} d\Omega = \sqrt{4\pi} \delta_{l0}\delta_{m0}$

The 2nd one is the orthonormal  condition.

$\displaystyle \int Y^{*}_{l'm'}Y_{lm} d\Omega = \delta_{l'l}\delta_{m'm}$

The 3rd one is triplet integral, we use the product of spherical harmonic,

$\displaystyle \int Y_{l_1m_1}Y_{l_2m_2} Y^*_{l_3m_3} d\Omega \\ = \int \sum_{lm} \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi(2l+1)}} C_{l_10l_20}^{l0} C_{l_1m_1l_2m_2}^{lm} Y_{lm} Y^*_{l_3m_3} d\Omega \\= \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi(2l_3+1)}} C_{l_10l_20}^{l_30} C_{l_1m_1l_2m_2}^{l_3m_3}$

The 4th one is another triple integral,

$\displaystyle \int Y_{l_1m_1}Y_{l_2m_2} Y_{l_3m_3} d\Omega \\ = \int \sum_{lm} \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi(2l+1)}} C_{l_10l_20}^{l0} C_{l_1m_1l_2m_2}^{lm} Y_{lm} Y_{l_3m_3} d\Omega \\ = \int \sum_{lmLM} \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi(2l+1)}} C_{l_10l_20}^{l0} C_{l_1m_1l_2m_2}^{lm} \\ \sqrt{\frac{(2l+1)(2l_3+1)}{4\pi(2L+1)}} C_{l0l_30}^{L0} C_{lml_3m_3}^{LM}Y_{LM}d\Omega$

$\displaystyle = \sum_{lm} \sqrt{\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi}} C_{l_10l_20}^{l0} C_{l_1m_1l_2m_2}^{lm} C_{l0l_30}^{00} C_{lml_3m_3}^{00}$

Notice that

$C_{lmLM}^{00} = (-1)^{L+M} \sqrt{\frac{1}{2L+1}} \delta_{Ll}\delta_{-m,M}$

$\displaystyle \int Y_{l_1m_1}Y_{l_2m_2} Y_{l_3m_3} d\Omega = \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi (2l_3+1)}} C_{l_10l_20}^{l_30} C_{l_1m_1l_2m_2}^{l_3,-m_3} (-1)^{m_3}$

using Wigner 3-j symbol,

$C_{l_1m_1l_2m_2}^{l_3m_3} = (-1)^{l_1-l_2+m_3} \sqrt{2l_3+1} \begin{pmatrix} l_1 & l_2 & l_3 \\ m_1 & m_2 & -m_3 \end{pmatrix}$

$\displaystyle \int Y_{l_1m_1}Y_{l_2m_2} Y_{l_3m_3} d\Omega \\= \sqrt{\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi}} \begin{pmatrix} l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{pmatrix}$

For other integral, we can use

$Y^*_{lm}(\theta, \phi) = (-1)^{m}Y_{l(-m)}(\theta,\phi) = Y_{lm}(\theta, -\phi)$