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)$$
This seems to be the Lambertian reflectance for parallelly polarised light passing from medium 0 to medium 1, so let us introduce $n = \frac{n_1}{n_0} > 1$ and write your integral as
\begin{align}
r_\text{p} (n) &= \int \limits_0^{\pi/2} \left(\frac{n^2 \cos(f) - \sqrt{n^2 - \sin^2(f)}}{n^2 \cos(f) + \sqrt{n^2-\sin^2(f)}}\right)^2 2 \sin(f) \cos(f) \, \mathrm{d} f \\
&\!\!\!\!\!\!\stackrel{\sin^2(f) = t}{=} \int \limits_0^1 \left(\frac{1 - \frac{n^2 \sqrt{1-t}}{\sqrt{n^2-t}}}{1 + \frac{n^2 \sqrt{1-t}}{\sqrt{n^2-t}}}\right)^2 \, \mathrm{d} t \stackrel{\frac{n \sqrt{1-t}}{\sqrt{n^2-t}} = u}{=} 2 n^2 (n^2 - 1) \int \limits_0^1 \frac{u (1-nu)^2}{(1+nu)^2(n+u)^2(n-u)^2} \, \mathrm{d} u \\
&= 2 n^2 (n^2 - 1) \int \limits_0^1 \left[\frac{8 n^3(n^4+1)}{(n^4-1)^3(1+nu)} - \frac{4n^3}{(n^4-1)^2 (1+nu)^2} \right. \\
&\phantom{= 2 n^2 (n^2 - 1) \int \limits_0^1 \left[\vphantom{\frac{8 n^3(n^4+1)}{(n^4-1)^3(1+nu)}}\right.} \left. - \, \frac{n^2-1}{(n^2+1)^3 (n-u)} + \frac{(n^2-1)^2}{4(n^2+1)^2n(n-u)^2} \right. \\
&\phantom{= 2 n^2 (n^2 - 1) \int \limits_0^1 \left[\vphantom{\frac{8 n^3(n^4+1)}{(n^4-1)^3(1+nu)}}\right.} \!\left. - \, \frac{n^2+1}{(n^2-1)^3(n+u)} - \frac{(n^2+1)^2}{4(n^2-1)^2 n (n+u)^2}\right] \mathrm{d} u \, .
\end{align}
The remaining integrals are elementary and after some simplification we end up with
$$ r_\text{p} (n) = 1 - \frac{4 n^3 (n^2+2n-1)}{(n^2-1)(n^2+1)^2} + \frac{16 n^4 (n^4+1)}{(n^2-1)^2(n^2+1)^3} \ln(n) - \frac{4 n^2 (n^2-1)^2}{(n^2+1)^3} \operatorname{arcoth} (n) \, .$$
Best Answer
Probably not the most direct way to obtain the result, but one possible method: From the integral representation DLMF \begin{equation} \int_{0}^{\infty}\frac{e^{-a^2t}}{\sqrt{t+\cos^{2}\theta}}\mathrm{d}t=\frac{\sqrt{\pi}}{a}e^{a^2\cos^{2}\theta}\operatorname{erfc}\left(a\cos\theta\right) \tag{1}\label{1} \end{equation} the integral can be transformed into \begin{equation} I=\frac{a}{\sqrt{\pi}}\int_{0}^{\pi/2}\cos \theta \,d\theta \int_{0}^{\infty}\frac{e^{-a^2t}}{\sqrt{t+\cos^{2}\theta}}\,\mathrm{d}t \end{equation} By changing the integration order, \begin{align} I&=\frac{a}{\sqrt{\pi}}\int_{0}^{\infty}e^{-a^2t}\,dt\int_{0}^{\pi/2}\frac{\cos \theta }{\sqrt{t+1-\sin^2\theta}}\,d\theta\\ &=\frac{a}{\sqrt{\pi}}\int_{0}^{\infty}e^{-a^2t}\arcsin \frac{1}{\sqrt{t+1}}\,dt \end{align} Now, integrating by parts, ($u'=e^{-a^2t},v=\arcsin \frac{1}{\sqrt{t+1}}$), we get \begin{align} I&=\frac{a}{\sqrt{\pi}}\left[\frac{\pi}{2a^2}-\frac{1}{2a^2}\int_{0}^{\infty}\frac{e^{-a^2t}}{\sqrt{t}(t+1)} \right]\\ &=\frac{a}{\sqrt{\pi}}\left[\frac{\pi}{2a^2}-\frac{1}{a^2}\int_{0}^{\infty}\frac{e^{-a^2u^2}}{u^2+1} \right] \end{align} (by changing $t=u^2$). Using the integral representation DLMF \begin{equation} \operatorname{erfc}(a)=\frac{2}{\pi}e^{-a^{2}}\int_{0}^{\infty}\frac{e^{-a^{2}u^% {2}}}{u^{2}+1}\mathrm{d}u \tag{2}\label{2} \end{equation} the final result is obtained: \begin{equation} I= \frac{\sqrt{\pi } \left(1-e^{a^2} \operatorname{erfc}\left(a\right)\right)}{2 a} \end{equation} \eqref{1} This identity is retrieved by changing $t=u^2-\cos^2\theta$ in the integral.
\eqref{2} With $A=a^2$, defining \begin{align} J(A)&=\frac{2}{\pi}\int_0^\infty\frac{e^{-A\left( u^2+1 \right)}}{u^2+1}\\ \frac{dJ(A)}{dA}&=-\frac{2}{\pi}\int_0^\infty e^{-A\left( u^2+1 \right)}\,du\\ &=-\frac{1}{\sqrt{\pi A}}e^{-A} \end{align} and remarking that $J(0)=1$, \begin{align} J(A)&=1-\frac{1}{\sqrt{\pi}}\int_0^A\frac{e^{-s}}{\sqrt{s}}\,ds\\ &=1-\operatorname{erf}(a) \end{align}