Hessian matrix in spherical coordinates This seems like a straightfoward question but I cannot find the answer anywhere.
I need to implement the Hessian matrix of a real scalar function f (an Hamiltonian, to be specific) in spherical coordinates (1,$\theta$,$\phi$) on the unit sphere and evaluate it at a turning point where the first derivatives vanish ($\frac{\partial f}{\partial \theta} =\frac{\partial f}{\partial \phi}=0 $). 
My first try was $$
H=
\begin{pmatrix}
    \frac{\partial^2 f}{\partial \theta^2} &  \frac{\partial^2 f}{\partial \theta \partial \phi} \\
    \frac{\partial^2 f}{\partial \theta \partial \phi} &  \frac{\partial^2 f}{\partial \phi^2} 
  \end{pmatrix},
$$
which is what I found in many (physics) papers where $\theta$ is conveniently set to $\frac{\pi}{2}$.
However, what I observed numerically is that because the elementary displacement on the sphere $\delta \vec{m}$ caused by a variation d$\phi$ depends on $\theta$, this is valid at the equator and then gradually breaks down closer to the poles because the  $\frac{\partial^2 f}{\partial \phi^2}$ and $\frac{\partial^2 f}{\partial \theta \partial \phi}$ terms tend to $0$ when they shouldn't.
This problem is corrected in the spherical gradient by the $\frac{1}{\sin \theta}$ factor on the $\hat{e}_{\phi}$ component, ie:
$$
\nabla f = \frac{\partial f}{\partial \theta}\hat{e}_{\theta}+\frac{1}{sin \theta} \frac{\partial f}{\partial \phi}\hat{e}_{\phi}.$$
From various sources (incl. wikipedia), the Hessian is sometimes defined as the Jacobian of the gradient,  which in spherical coordinates would be $H=J(\nabla f)= \big[ \frac{\partial \nabla f}{\partial \theta} ,\frac{\partial \nabla f}{\partial \phi} \big]$.
If I apply this at a turning point where the first derivatives vanish, I end up with
$$ H=
\begin{pmatrix}
    \frac{\partial^2 f}{\partial \theta^2} &  \frac{\partial^2 f}{\partial \theta \partial \phi} \\
   \frac{1}{\sin \theta}  \frac{\partial^2 f}{\partial \theta \partial \phi} &  \frac{1}{\sin \theta} \frac{\partial^2 f}{\partial \phi^2} 
  \end{pmatrix}.
$$
This is not symmetric anymore, which bothers me because in every book I looked, it says the Hessian is symmetric as long as the mix derivatives commute.
Does this mean this definition does not apply in spherical coordinates? Is there a more general definition I can use? Additionally, I found this post which defines $H=\nabla \nabla^T$. This yields
$$ H=
\begin{pmatrix}
    \frac{\partial^2 f}{\partial \theta^2} & \frac{1}{\sin \theta} \frac{\partial^2 f}{\partial \theta \partial \phi} \\
   \frac{1}{\sin \theta}  \frac{\partial^2 f}{\partial \theta \partial \phi} &  \frac{1}{\sin^2 \theta} \frac{\partial^2 f}{\partial \phi^2} 
  \end{pmatrix}
$$
which is a lot more intuitive, but it clashes with the previous Jacobian of the gradient definition
Any help/clarifications or suggestions greatly appreciated.
 A: It sounds like you want the Hessian of a Riemannian manifold.
(If you just want the Hessian in spherical coordinates in $\mathbb{R}^3$, then you should just use the chain rule on the regular Hessian).
In this case, the Hessian is written: $$ \mathcal{H}[f] = (\partial_{ij}f-\Gamma^k_{ij}\partial_kf)[dx^i\otimes dx^j] $$
In local coordinates, the unit sphere is given by: $$ \Psi(\theta,\phi)=(\cos(\theta)\sin(\phi),\sin(\theta)\sin(\phi),\cos(\phi)) $$
So the metric tensor is written: $$ g=\begin{bmatrix}\sin^2(\phi) & 0 \\ 0 & 1\end{bmatrix} $$
with Christoffel symbols given by (e.g. see here):
\begin{align} 
\Gamma_{11}^1 &= \Gamma_{22}^1 = \Gamma_{12}^2 = \Gamma_{22}^2=0\\[2mm]
\Gamma_{12}^1 &= \frac{\cos(\phi)}{\sin(\phi)}\\[1mm]
\Gamma_{11}^2 &= -\sin(\phi)\cos(\phi)
\end{align}
which gives me the result:
$$
\mathcal{H}[f] = \begin{bmatrix}
\partial_{\theta\theta}f + \cos(\phi)\sin(\phi)\partial_\phi f & \;\partial_{\theta\phi}f-\dfrac{\cos(\phi)}{\sin(\phi)}\partial_\theta f\\
\partial_{\theta\phi}f-\dfrac{\cos(\phi)}{\sin(\phi)}\partial_\theta f & \;\partial_{\phi\phi}f
\end{bmatrix}
$$
Apologies if there are any computation errors.
See also: [1], [2].
One check is that the Laplace-Beltrami operator is the trace of the Hessian, which is written as
$$ \text{tr}(\mathcal{H}[f]) = g^{ij}\mathcal{H}[f]_{ij} = 
\frac{1}{\sin^2(\phi)}\left[ \partial_{\theta\theta}f+\sin(\phi)\cos(\phi)\partial_\phi f \right] + \partial_{\phi\phi}f $$
From here, it is given by:
$$
\Delta_g f = (\sin(\phi))^{-1}{\partial}_{\phi}\left(
\sin(\phi)\partial_\phi f
\right) + (\sin(\phi))^{-2}\partial_{\theta\theta}f
$$
which corroborates it.
