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?
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.