I think the answer is yes, but you probably have to express your question in a more generalzied form. i.e. Let $M$ be a Riemannian manifold and $f\in C^{\infty}(M)$. Consider a smooth level set of $f$, i.e. say $\Sigma:=f^{-1}(0)$. Let $\nabla f$ be the gradient vector of $f$. Then the hypersurface $\Sigma$ has fundamental form
$$II(v,w)=\langle \nabla_v \vec{n}, w \rangle$$ where $\vec{n}$ is the normal vector of $\Sigma$ which is exatctly $\nabla f$. And Hessian can be written as $$Hess(f)(v, w)=\langle \nabla_v \nabla f, w \rangle$$
where $v, w\in T_p(\Sigma)$ and $p\in \Sigma$. If you are not familer with the definition of Hessian or Second fundamental form, you can find them in any basic Riemannian Geometry book, for example Kobayashi-Nomizu's Fundamental of Diff. Geo. II
(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*}
$$
Best Answer
Regarding question (1): When $M$ is compact and $f$ is a positive $C^2$ function, such an $h$ exists if and only if $f$ is constant on each component of $M$.
Proof: Suppose $M$ is compact and $f\colon M\to \mathbb R$ is a positive $C^2$ function. If $f$ is constant on each component, clearly $h\equiv 0$ works. Conversely, suppose $h\colon M\to \mathbb R$ is a $C^2$ function such that $$H^h(X,Y)=\frac 1fH^f(X,Y).$$ Taking the trace of both sides, we conclude that $\Delta h = \Delta f/f$. Now compute $$ \Delta ( h - \log f)= \frac{\Delta f}{f} - \left(\frac{\Delta f}{f} - \frac{|\text{grad}\, f|^2}{f^2}\right) = \frac{|\text{grad}\, f|^2}{f^2}. $$ Integrating both sides over $M$, we conclude that $|\text{grad}\, f|^2/f^2\equiv 0$ (because the integral of a Laplacian on a compact manifold is zero). This implies that $f$ is constant on each component of $M$. $\square$
When $M$ is noncompact, it might be possible to find such an $h$. For example, if $M=\mathbb R$ (with the Euclidean metric) and $f(x)=e^x$, we can take $h(x) = \tfrac12 x^2$. At the moment, I can't think of any necessary or sufficient conditions on $M$ or $f$ in that case.