Definition of the Hessian on Riemannian manifolds I can't understand the following definition given by Petersen in his Riemannian Geometry book.

Let $M$ be a Riemannian manifold and let $f \colon M \to \mathbb{R}$ be a smooth function. Let $\nabla f$ denotes its gradient. Then the Hessian Hess$f$ is defined as the symmetric $(0,2)$-tensor $\frac{1}{2}L_{\nabla f}g$, where $g$ is the metric on $M$, and $L$ is the Lie derivative.

How can I prove that in $\mathbb{R}^n$ this does coincides with the usual definition?
Moreover the book says also that we can define the Hessian as a self-adjoint $(1,1)$-tensor by $S(X) = \nabla_X \nabla f$ and the two tensors are naturally related by:
$$
\text{Hess}f(X,Y) = g(S(X),Y) \qquad \text{for every smooth vector fields $X$, $Y$.}$$
How can I check that?
 A: (1) Note that $L_X Y = [X,Y]= \nabla_X Y - \nabla_Y X$. Moreover, on $\mathbb{R}^n$, the covariant derivative is given by $\nabla_X Y = (X(Y_1), \ldots, X(Y_n))$. If $\partial_i$ is the $i$-th coordinate vector field on $\mathbb{R}^n$, then
$$
 \begin{align*} 
 L_{\nabla f}g(\partial_i,\partial_j) &= (\nabla f)(g(\partial_i,\partial_j))-g([\nabla f,\partial_i],\partial_j) - g(\partial_i, [\nabla f, \partial_j]) \\
&= (\nabla f)(\delta_{ij}) + g(\partial_i \nabla f,\partial_j) + g(\partial_i, \partial_j \nabla f) \\
&= 0 + \partial_i\partial_j f + \partial_j \partial_i f = 2\partial_i\partial_j f.
 \end{align*}
$$
And this is equal to the $(i,j)$-th entry of the Hessian matrix of $f$ (apart from the factor 2 of course).
(2) For all smooth vector fields $X$, $Y$ we have
$$
 \begin{align*}
   L_{\nabla f}g(X,Y) &= (\nabla f)(g(X,Y)) - g([\nabla f,X], Y)- g(X,[\nabla f, Y])\\
&= g(\nabla_{\nabla f} X, Y) +g(X,\nabla_{\nabla f}Y)  \\
  &\qquad- g(\nabla_{\nabla f}X-\nabla_X \nabla f, Y)-g(X,\nabla_{\nabla f}Y- \nabla_Y\nabla f)\\
&=  g(\nabla_X \nabla f, Y) + g(X,\nabla_Y\nabla f).
 \end{align*}
$$
If we show that the two last terms are equal, then we are done. This follows from the definition of a gradient and the fact that the covariant derivative is torsion free (because we are working with the Levi-Civita connection).
$$
  \begin{align*}
   g(\nabla_X \nabla f,Y) -g(X,\nabla_Y\nabla f) &= X g(\nabla f,Y) - g(\nabla f, \nabla_X Y) \\
& \qquad- Y g(\nabla f,X) + g(\nabla f,\nabla_Y X) \\
&= X(Y(f)) - Y(X(f)) - [X,Y](f) = 0
  \end{align*}
$$
