I am interested in the following integral $$\int_0^{2\pi}\int_0^\pi\mathop{\mathrm{d}\theta}\mathop{\mathrm{d}\phi} \sin\theta Y_l^{m*}(\theta,\phi)Y_{l'}^{m'}(\theta, \phi)\cos^2\theta\cos^2\phi,$$ where $Y_l^{m*}$ is the complex conjugate of spherical harmonics $Y_{l}^m$.
Any hints on how to proceed? Thanks.
A common approach to solve these types of problems is to expand the integrand until we have a sum of products of two spherical harmonics and then using the orthogonality relation
$$\int \overline{Y}^{m}_\ell Y^{m'}_{\ell'}{\rm d}\Omega = \delta_{\ell\ell'}\delta_{mm'}\tag{1}$$
to evaluate the integrals in the sum. I will here give the general outline for how to evaluate the integral
$$\int \overline{Y}^m_\ell Y^{m'}_{\ell'} f(\theta,\phi)\,{\rm d}\Omega$$
where $f(\theta,\phi) = \cos^2\theta\cos^2\phi$ for this particular question. Here and below I use the convention ${\rm d}\Omega = \sin\theta\,{\rm d}\theta\,{\rm d}\phi$, $\sum\limits_{\ell,m} \equiv \sum\limits_{\ell=0}^\infty\sum\limits_{m=-\ell}^\ell$ and an overbar denotes complex conjugation.
We start with a well known, and very useful, result (see e.g this page). We can expand a product $Y^{m_1}_{\ell_1}Y^{m_2}_{\ell_2}$ in a series of spherical harmonics as follows
$$Y^{m_1}_{\ell_1}Y^{m_2}_{\ell_2} = \sum_{\ell,m}\sqrt{\frac{(2\ell_1+1)(2\ell_2+1)(2\ell+1)}{4\pi}}\begin{pmatrix} \ell_1 & \ell_2 & \ell \\[8pt] 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_1 & \ell_2 & \ell\\ m_1 & m_2 & m \end{pmatrix}\overline{Y}^{m}_\ell$$
where $(\cdot)$ is the so-called Wigner $3j$ symbol. Since the spherical harmonics form a complete basis of the Hilbert space of square-integrable functions we can expand the function $f(\theta,\phi)$ in a series
$$f(\theta,\phi) = \sum_{\ell,m} a_{\ell m}Y^m_\ell~~~~\text{where}~~~a_{\ell m} = \int f(\theta,\phi) Y^m_\ell{\rm d}\Omega\tag{2}$$
Using the statements above we can write
$$\overline{Y}^m_\ell Y^{m'}_{\ell'}f(\theta,\phi) = \overline{Y}^m_\ell Y^{m'}_{\ell'} \sum_{L,M} a_{LM}Y^M_L = \overline{Y}^m_\ell \sum_{L,M} a_{LM}Y^{m'}_{\ell'}Y^M_L \\= \sum_{L,M}\sum_{L',M'}a_{LM}\sqrt{\frac{(2\ell'+1)(2L+1)(2L'+1)}{4\pi}}\begin{pmatrix} \ell' & L & L' \\[8pt] 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell' & L & L'\\ m' & M & M' \end{pmatrix}\overline{Y}^m_\ell\overline{Y}^{M'}_{L'}$$
and by using the orthogonality relation $(1)$ and $\overline{Y}^{m}_{\ell} = (-1)^m Y^{-m}_{\ell}$ get ridd of one of the conjugates we get that the integral evaluates to
$$\int \overline{Y}^m_\ell Y^{m'}_{\ell'} f(\theta,\phi)\,{\rm d}\Omega \\= \sum_{L,M}a_{LM}(-1)^m\sqrt{\frac{(2\ell'+1)(2L+1)(2\ell+1)}{4\pi}}\begin{pmatrix} \ell' & L & \ell \\[8pt] 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell' & L & \ell\\ m' & M & -m \end{pmatrix}$$
To solve your particular problem you'll need to compute $a_{\ell m}$ for $f(\theta,\phi) = \cos^2\theta\cos^2\phi$ using $(2)$. In this case it's only when $\ell$ is even and $m \in \{-2,2\}$ plus the case $m=0$ for $\ell\in \{0,2\}$ that will give rise to non-zero terms. Then in the end you can use some of the properties of the Wigner $3j$ symbol to simplify the final result (for example the last term in the last equation above is zero unless $m' + M = m$).