Given a surface defined by spherical harmonic terms
$f(\theta,\phi) = \sum_{l=0}^\infty\sum_{m=-l}^lf_l^mY_l^m(\theta,\phi)$
How can I evaluate the normal vector at each location on the surface? i.e. $\mathbf{N}(\theta,\phi)=N_r(\theta,\phi)\hat{\mathbf{r}}+N_\theta(\theta,\phi)\hat{\mathbf{\theta}}+N_\phi(\theta,\phi)\hat{\mathbf{\phi}}$
I imagine this answer lies somewhere in the vector spherical harmonics.. it feels like it should be true that the vector field produced by converting our scalar harmonics to vector terms as follows (using the notation from https://en.wikipedia.org/wiki/Vector_spherical_harmonics) would produce the desired normal field but I have no idea how to verify this. Or how one would go about taking the gradient of spherical harmonics.
$\mathbf{Y}_{lm}=Y_{lm}\hat{\mathbf{r}}$
$\mathbf{\Psi}_{lm}=r\nabla Y_{lm}$
$\mathbf{\Phi}_{lm}=r\times\nabla Y_{lm}$
Take the derivative of $f$ along $\frac{\partial}{\partial \theta},\frac{\partial}{\partial \phi}$. Now take the cross product of \begin{equation} \begin{pmatrix} 1\\0\\\frac{\partial f}{\partial \phi} \end{pmatrix}\times\begin{pmatrix} 0\\1\\\frac{\partial f}{\partial \theta} \end{pmatrix}=N \end{equation} Normalize and you are done.