We are solving the integral:
$$I=\int_{0}^{\pi} \int_{0}^{2 \pi} (\cos{\phi})^{n_1}(\sin{\phi})^{n_2}(\cos{\theta})^{n_3}(\sin{\theta})^{n_4}Y^u_p Y^v_q Y^w_r \sin{\theta}\, d\phi\, d\theta.$$
Due to the nature of the spherical harmonics, which can be written as a product of associated Legendre polynomials and an exponential :
$$Y_l^m(\theta,\phi)=\sqrt{\frac{2l+1}{4\pi}\cdot\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta)e^{im\phi},$$
you can split the integral in two parts:
$$\Phi(n_1,n_2,u+v+w)=\int_0^{2\pi}(\cos\phi)^{n_1} (\sin\phi)^{n_2} e^{i(u+v+w)\phi}d\phi,$$
$$\Theta(n_3,n_4,p,q,r,u,v,w)=\int^1_{-1} x^{n_3}(1-x^2)^\frac{n_4}{2} P_p^u(x)P_q^v(x)P_r^w(x)dx.$$
The first integral $\Phi$, can be solved by converting the cosine and sine into their exponential form and using the Binomial theorem and taking into account that
$$\int_0^{2\pi} e^{im\phi}d\phi=2\pi\delta_{m,0}.$$
Thus,
$$\Phi(n_1,n_2,m)=\frac{2\pi}{2^{n_1+n_2}i^{n_2}}\sum_{k_1=0}^{n_1}\sum_{k_2=0}^{n_2}(-1)^{n_2-k_2}\binom{n_1}{k_1}\binom{n_2}{k_2}\delta_{n_1+n_2-m,2(k_1+k_2)}.$$
This gives you your first selection rules:
- if $n_1+n_2-m$ is odd, then $\Phi=0$,
- if $n_1+n_2-m<0$, then $\Phi=0$,
- if $-m > n_1+n_2$, then $\Phi=0$.
To solve the integral of $\Theta$, you should try to reduce it to the form :
$$\int^1_{-1}P_p^u(x)P_q^v(x)P_r^w(x)dx,$$
which has a closed-form expression. It is related to 3j-symbols/Clebsch-Gordan coefficients and is known as Gaunt's Formula (See Dong S.H., Lemus R., (2002), "The overlap integral of three associated Legendre polynomials", Appl. Math. Lett. 15, 541-546.).
The integral of $\Theta$ can be reduced to this form using the recurrence relations of the associated Legendre polynomials:
$$ P_{l-1}^m(x) - P_{l+1}^m(x) = (2l+1)\sqrt{1-x^2}P_l^{m-1}(x)$$
$$ (2l+1)xP_l^m(x) = (l-m+1)P^m_{l+1}(x)+(l+m)P^m_{l-1}(x)$$