I'm looking for the proof of very nice identity about the Laplace-Beltrami operator of the Gauss map $N$ of a regular surface in $\mathbb{R}^3$ given by a patch $X$. I want to show that
$$\Delta N = -2g^{ij} \partial_i H \partial_j X - 2(2H^2-K)N,$$
where $\Delta$ denotes the Laplace-Beltrami operator on $M$ (applied componentwise on the Gauss map $N$). Unfortunately I did not find a proof for this in any book on differential geometry. It would be awesome if somebody could give a hint or a reference!
Assuming that we are calculating at a normal coordinate $(x_1, x_2)$ at $x$. Then
$$\begin{split}\nabla_i N &= \langle \nabla_i N, e_1\rangle e_1 + \langle \nabla_i N, e_1\rangle e_1 \\ &= -h_{1i} e_1 -h_{2i} e_2 \end{split}$$
Differentiate again,
$$\nabla_i \nabla _i N = -(\nabla_i h_{1i}) e_1 - (\nabla_i h_{2i}) e_2 - (h_{1i})^2N - (h_{2i})^2 N . $$
Using the Codazzi equation, $\nabla_i h_{jk} = \nabla_jh_{ik}$, we have \
$$\nabla_i \nabla _i N = -( \nabla_1 h_{ii} e_1 - \nabla_2 h_{ii}) e_2 -( h_{1i}^2 + h_{2i}^2) N . $$
If we sum over $i =1, 2$, using $H = \frac 12 (h_{11} + h_{22})$ and $K = h_{11} h_{22} - h_{12}^2$, we have
$$\Delta N = \sum \nabla_ i \nabla_i N = -2\nabla_i H e_i -2(2H^2 -K) N,$$
which is what you want.