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.