Orthogonality of spherical harmonics under a rotation Consider the orthonormalized spherical harmonics defined with respect to Cartesian axes $(x,y,z)$ which obey the relation,
$$
\int d\Omega Y_l^m Y_{l'}^{m'*} = \int_0^{2\pi} d\phi \int_0^{\pi} d\theta \sin{\theta} Y_l^m(\theta,\phi) Y_{l'}^{m'*}(\theta,\phi) = \delta_{ll'} \delta_{mm'}
$$
Suppose now I define new Cartesian axes $(x',y',z')$ which are a rotation of the original $(x,y,z)$. Now I want to integrate products of spherical harmonics defined on the new coordinates $(\theta',\phi')$ but the integral should be done with respect to the original coordinates. In other words, I want to calculate:
$$
\int_0^{2\pi} d\phi \int_0^{\pi} d\theta \sin{\theta} Y_l^m(\theta',\phi') Y_{l'}^{m'*}(\theta',\phi')
$$
where $\theta' = \theta'(\theta,\phi)$ and $\phi' = \phi'(\theta,\phi)$ are the new spherical coordinates defined wrt $(x',y',z')$.
I have done numerical experiments which indicate that the new integral is also equal to
$$
\delta_{ll'} \delta_{mm'}
$$
but I don't see an obvious reason why this should be true. For an arbitrary rotation, the equations for $\theta',\phi'$ in terms of $\theta,\phi$ get somewhat complicated.
Is there an easy way to see why the spherical harmonics should be orthogonal with respect to an integral over rotated coordinates?
 A: According to Steinborn and Ruedenberg 1973, Eq. 189, under a rigid rotation with Euler angles $\alpha,\beta,\gamma$, a spherical harmonic of degree $l$ transforms as,
$$
Y_l^m(\theta',\phi') = \sum_{m'=-l}^l D_{m'm}^{(l)}(\alpha,\beta,\gamma) Y_l^{m'}(\theta,\phi)
$$
where the $D^{(l)}$ matrices denote the $(2l+1)$ dimensional irreducible represenation of the rotation group. Explicit expressions for the elements $D_{m'm}^{(l)}$ are given in Eqs. 185 and 201 of the paper.
Working from this expression, we would find
\begin{align}
\int d\Omega Y_l^m(\theta',\phi') Y_k^{n*}(\theta',\phi') &= \sum_{m'=-l}^l \sum_{n'=-k}^k D_{m'm}^{(l)}(\alpha,\beta,\gamma) D_{n'n}^{(k)}(\alpha,\beta,\gamma) \int d\Omega Y_l^{m'}(\theta,\phi) Y_k^{n'*}(\theta,\phi) \\
&= \delta_{lk} \sum_{m'=-l}^l D_{m'm}^{(l)}(\alpha,\beta,\gamma) D_{m'n}^{(l)}(\alpha,\beta,\gamma)
\end{align}
This shows that the integral vanishes for $l \neq k$.
Above Eq. 193, the authors state that the matrices $D^{(l)}$ are unitary. This means,
\begin{align}
\sum_{m'=-l}^l D_{m'm}^{(l)} D_{m'n}^{(l)} &= \sum_{m'=-l}^l (D^{(l)T})_{mm'} D^{(l)}_{m'n} \\
&= (D^{(l)T} D^{(l)})_{mn} \\
&= \delta_{mn}
\end{align}
which proves the result
$$
\int d\Omega Y_l^m(\theta',\phi') Y_k^{n*}(\theta',\phi') = \delta_{lk} \delta_{mn}
$$
for any rigid rotation.
