# Interesting special functions identity involving the inner product of real spherical harmonics with a cosecant weight function

In spherical coordinates $\Omega=(\theta,\phi)\in[0,\pi]\otimes[0,2\pi]$, define the inner product $$C_{L_1m_1}^{L_2m_2}:=\left\langle Y_{L_1m_1},\rho,Y_{L_2m_2}\right\rangle=\int Y_{L_1m_1}(\Omega)\rho(\Omega)Y_{L_2m_2}(\Omega)\,d\Omega$$ where $Y_{Lm}(\Omega)$ is the real spherical harmonic $$Y_{Lm}(\theta,\phi)=(-1)^m\sqrt{\frac{(2L+1)}{2\pi}\frac{(L-|m|)!}{(L+|m|)!}}P_L^{|m|}(\cos(\theta))\times\begin{cases} \sin(|m|\phi)&m<0\\ \sqrt{\frac{1}{2}}&m=0\\ \cos(|m|\phi)&m>0 \end{cases}$$ and $\rho(\theta,\phi)=\csc^2(\theta)$.

Testing in Mathematica has led me to conjecture that $$C_{L_1m_1}^{L_2m_2}=\frac{L_1+1/2}{|m_1|}\delta_{L_1L_2}\delta_{m_1m_2}$$ except when $m_1=m_2=0$, in which case the integral does not converge. So far, this appears to be valid for all integers up to $0\leq L_1,L_2\leq5$ and $0\leq |m_1|\leq L_1,0\leq |m_2|\leq L_2$ (ie, Mathematica was able to symbolically verify it for the first 1296 cases).

Assuming that this is true, is there a systematic way to prove it? I am vaguely aware of the usual identities involving the inner product of two real harmonics $\left\langle Y_{L_1m_1},Y_{L_2m_2}\right\rangle=\delta_{L_1,L_2}\delta_{m_1m_2}$ and the more complicated one involving the product of three (involving Clebsch-Gordan coefficients), but I don't see how they're directly applicable to the case where there's a $\frac{1}{\sin(\theta)^2}$ multiplying the two harmonics.

Also, is there a general formula where the weight function can be other trigonometric functions, such as $\rho(\theta,\phi)=\sin^{-1}(\theta),\cot(\theta)$, etc?

This falls under Sturm-Liouville Theory of second order symmetric differential operator $$LF=\frac{1}{w}\frac{d}{dx}\rho\frac{d}{dx}F,\;\;\; a < x < b,$$ where $\rho$ and $w$ are positive on $(a,b)$. Such operators have a useful so-called adjoint equation $$(LF)G-F(LG) = \frac{1}{w}\frac{d}{dx}\{\rho(F'G-FG')\}.$$ It's the adjoint equation that leads to orthogonality, provided conditions are carefully chosen at the endpoints. Specifically, suppose that $F$, $G$ are eigenfunctions of $L$ with eigenvalues $\lambda$, $\mu$, respectively; i.e., suppose $LF=\lambda F$ and $LG=\mu G$. Then $$(\lambda-\mu)wFG = \frac{d}{dx}\{\rho(F'G-FG')\},$$ which leads to $$\left.(\lambda-\mu)\int_{a}^{b}wFG\,dx = \rho(F'G-FG')\right|_{a}^{b}$$ If $F$ and $G$ are in $L^{2}_{w}(a,b)$, meaning that $$\int_{a}^{b}w|F|^{2}\,dx < \infty,\;\;\; \int_{a}^{b}w|G|^{2}\,dx < \infty,$$ then $(F,G)_{w}=\int_{a}^{b}wFG\,dx$ converges absolutely and, automatically, the evaluation terms must exist as $x$ approaches $a$ or $b$. By imposing conditions are the endpoints to make the evaluation terms vanish, you end up with $(\lambda -\mu)(F,G)=0$, which forces $(F,G)=0$ for eigenfunctions with different eigenvalues. In some cases you don't need conditions at the endpoints at all, which is your case. This can happen when $\rho$ vanishes at the endpoints of the interval (but it doesn't have to happen): you can see some reason how this might be related to the adjoint equation.
The equation you have for the the functions of $\theta$ is really unusual in that it may be viewed as an eigenfunction equation in two ways. That's how you end up with two sets of orthogonality conditions. You've just discovered the second orthogonality set of orthogonality conditions. Setting $x=\cos\theta$ transforms your equation to the form $$\left[\frac{d}{dx}(1-x^{2})\frac{d}{dx}-\frac{m^{2}}{1-x^{2}}+n(n+1)\right]P_{n}^{m}(x)=0,$$ For fixed $m$, you end up with weight $w=1$ and eigenvalues $n(n+1)$. For fixed $n$, you end up with weight $w=1/(1-x^{2})$ and eigenvalues $m^{2}$.
Except for the case where $m=0$, I know of no simple way to get the normalization constants. It's all a bit of tedious integration by parts applied to explicit forms of $P_{n}^{m}$ expressed using the Rodrigues formula.