In my studies I've come across the hessian in the context of Riemannian geometry. I use the following definition of the hessian
$$ H^f(X,Y)=XYf-(\nabla_XY)f=\langle \nabla_X(\operatorname{grad} f),T\rangle. $$
I here want to show that at a critical point $p$ of $f$ we get the following
$$ H^f(X_p,Y_p)=X_p(Yf)=Y_p(Xf),\quad \text{for all } X,Y. $$
I'm not quite sure how to start.
I'll record a quick answer here: if $M$ is a smooth manifold and $\nabla$ is a linear connection in $M$, then the covariant hessians according to $\nabla$ are symmetric tensors if and only if $\nabla$ is torsion free. This happens, for instance, if $(M,\langle\cdot,\cdot\rangle)$ is a pseudo-Riemannian manifold and $\nabla$ is its Levi-Civita connection.
If given a tensor field $A \in \mathscr{T}^r_{\;\;s}(M)$, we define $\nabla A \in \mathscr{T}^r_{\;\;s+1}(M)$ by $$\begin{align}\nabla A(&\omega^1,\ldots, \omega^r, X_1,\ldots,X_s;X) = \nabla_XA(\omega^1,\ldots,\omega^r,X_1,\ldots,X_s)\\ &= X(A(\omega^1,\ldots, \omega^r,X_1,\ldots, X_s)) - \sum_{i=1}^r A(\omega^1,\ldots, \nabla_X\omega^i, \ldots, \omega^r,X_1,\ldots,X_r) \\ & \qquad\qquad - \sum_{i=1}^r A(\omega^1,\ldots, \omega^r, X_1,\ldots \nabla_XX_i,\ldots, X_s),\end{align}$$then by definition $${\rm Hess}(f)(X,Y) \doteq \nabla(\nabla f)(X;Y) = Y(X(f)) - {\rm d}f(\nabla_YX).$$If $\tau(X,Y) = \nabla_XY - \nabla_YX - [X,Y]$ is the torsion of the connection, we have $${\rm Hess}(f)(X,Y) = {\rm Hess}(f)(Y,X) + {\rm d}f(\tau(X,Y)).$$If we have no torsion, and $p$ is a critical point of $f$, then ${\rm d}f_p = 0$, and so $${\rm Hess}(f)_p(X_p,Y_p) = {\rm Hess}(f)_p(Y_p,X_p) = X_p(Y(f)) = Y_p(X(f)).$$