Tal vez no sea realmente una respuesta, pero puedes captar la idea: esto es cierto más o menos por construcción. Los armónicos esféricos (como ejemplo) se obtienen como funciones propias de la parte angular del operador de Laplace, es decir, satisfacen $$\Delta_{S^2} Y_{lm}(\vartheta,\phi) = \lambda Y_{lm} (\vartheta,\phi)$$ (En realidad resulta que esto implica $\lambda = -l(l+1)$ con un número entero $l$ ) Si se tienen tales funciones propias para diferentes valores propios, es cuestión de álgebra lineal demostrar que son ortogonales, observando $$\int_{S^2}\langle \nabla_{S^2}Y_{lm}, \nabla_{S^2}Y_{l'm'}\rangle d\mu_{S^2}= -\int_{S^2}\langle Y_{lm}, \Delta_{S^2}Y_{l'm'}\rangle d\mu_{S^2}$$ Esto implica que las funciones son ortogonales si $l\neq l'$ ya que de lo contrario se podría derivar $l(l+1) = l'(l'+1)$ de esto. Para los fijos $l$ resulta que se puede resolver la ecuación por una aproximación de separación que lleva a una EDO que se sabe que es resoluble por polinomios ortogonales por la teoría de las EDO.
También puede anotar el $Y_{ml}$ de forma bastante explícita, véase, por ejemplo, la página de la wikipedia alemana sobre "Kugelflächenfunktionen" http://de.wikipedia.org/wiki/Kugelfl%C3%A4chenfunktionen . Si se observan con más detenimiento y se tiene algo de experiencia en ODE, se puede observar que el $\vartheta$ son polinomios ortogonales bien conocidos en $\cos(\vartheta)$ (Polinomios de Legendre), mientras que el $\phi$ la parte es más o menos $e^{im\phi}$ que se conoce a un sistema de funciones ortogonales. Demostrar la ortogonalidad es todavía un poco de trabajo, pero te muestra la dirección que debes tomar.
Para hacerlos realmente ortonormal hay que normarlos, por supuesto, de ahí vienen los factores de aspecto complicado.