Gradient of a real valued function defined on a sphere 
Points on a sphere of radius $R$ is expressed in spherical coordinates as $\left(\varphi,\theta\right)$.
For a real valued, continuous and differentiable function $f:\mathbb R\times\left[0,\pi\right] \to \mathbb R$ evaluated on said sphere centered at $\left(0,0,0\right)$ with radius $R$, what is its gradient written in the form of angles (radians) in the azimuthal and zenith directions? 
Additional notes and assumptions:
  
  
*
  
*$f\left(\varphi+2\pi,\theta\right)=f\left(\varphi,\theta\right)$;
  
*$f\left(\varphi,0\right) = f\left(0,0\right)$, $f\left(\varphi,\pi\right) = f\left(0,\pi\right)$;
  
*$\varphi\in\mathbb R$ is the azimuthal angle and
  $\theta\in\left[0,\pi\right]$ is the zenith angle.
  

By design, the direction of the gradient has to be parallel to the tangent plane at any position, and points in the direction of greatest increase of $f$. Meanwhile, the magnitude of the gradient is the "slope" of $f$ in said direction.
I'm aware that the gradient operator of $\mathbb R^3$ in the spherical coordinate system can be written as a linear combination of the basis vectors, namely
$$
  \nabla =
 \mathbf{e}_{r}                             \frac{\partial}{\partial r}
   + \mathbf{e}_{\varphi} \frac{1}{r\sin\theta} \frac{\partial}{\partial \varphi}
   + \mathbf{e}_{\theta}  \frac{1}{r}           \frac{\partial}{\partial \theta}.
$$
Since $f$ is only mapped from points on the sphere, I removed the term with $\mathbf{e}_{r}$ to arrive at
$$
  \nabla f =
 \mathbf{e}_{\varphi} \frac{1}{R\sin\theta}\frac{\partial f}{\partial\varphi}
   + \mathbf{e}_{\theta}  \frac{1}{R}          \frac{\partial f}{\partial\theta},
$$
from where I naively translate the "arc lengths" into their respective "angles" by
\begin{align}
  \frac{1}{R  \sin  \theta} \frac{\partial f}{\partial\varphi}
\to
  \frac{1}{R^2\sin^2\theta} \frac{\partial f}{\partial\varphi},
\quad &
  \text{because the circle of latitude has radius of }R\sin\theta,
\text{and} \\
  \frac{1}{R}   \frac{\partial f}{\partial\theta}
\to
  \frac{1}{R^2} \frac{\partial f}{\partial\theta},
\quad &
  \text{because the circle of longitude has radius of }R.
\end{align}
Therefore, I believe that the gradient described in azimuthal and zenith angles is the following pair:
$$
  \frac{1}{R^2} \left(
\frac{1}{\sin^2\theta} \frac{\partial f}{\partial\varphi},
                       \frac{\partial f}{\partial\theta}
  \right).
$$

My questions would be:
  
  
*
  
*Is there a rigorous formulation for this kind of gradient operator that I'm looking for?
  
*What about gradient at the poles? I'm asking because $\sin\theta$ is obviously $0$ when $\theta=0,\pi$. 

 A: Rigorous formulation
This involves the definition of the surface gradient operator, which is defined as
$$
  \nabla_{\Gamma} = \nabla -
       {\mathbf e}_{r} \left({\mathbf e}_{r} \cdot \nabla\right).
$$
The projection of the gradient along the unit normal ${\mathbf e}_{r}$ is evaluated by
$$
  {\mathbf e}_{r} \left({\mathbf e}_{r} \cdot \nabla f\right)
    = {\mathbf e}_{r} \left(\frac{\partial f}{\partial r} {\mathbf e}_{r}^2\right)
    = {\mathbf e}_{r}       \frac{\partial f}{\partial r}
$$
which when subtracted from $\nabla f$ gives (at $r=R$)
$$
  {\nabla}_{\Gamma} f =
      \mathbf{e}_{\varphi} \frac{1}{R\sin\theta}\frac{\partial f}{\partial\varphi}
    + \mathbf{e}_{\theta}  \frac{1}{R}          \frac{\partial f}{\partial\theta}.
$$
Angular gradient field
Because of $\mathbf{e}_{\theta}\cdot\mathbf{e}_{\varphi} = 0$ and that the longitudinal rotation and the latitudinal rotations are independent of each other, the component of ${\nabla_{\Gamma}}f$ in the direction of $\mathbf{e}_{\theta}$ and $\mathbf{e}_{\varphi}$ can be individually resolved into angular components.
Polar regions
The original extent for asking the case of $\sin\theta=0$ is for numerical purposes. However, it should be noted that the gradient operator cannot be defined on $\sin\theta=0$ due to the choice of coordinate system and boundary conditions.
