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?

(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. \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*}