$\Delta \mathbf n = -2 \mathbf n$ on the Euclidean sphere Let us consider the Euclidean two-sphere, defined by the embedding in the three dimensional Euclidean space as
$$
\mathbf n \cdot \mathbf n = 1\,,
$$
where $\cdot$ denotes the standard scalar product. The metric on the sphere, in some coordinates $x^i$, is expressed as 
$$
\gamma_{ij}=\mathbf e_i \cdot \mathbf e_j\,,
$$
where $\mathbf e_i=\partial_i\mathbf n$. For instance, in the standard spherical coordinates 
$$
\mathbf n=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)
$$
and
$$
\gamma_{\theta\theta}=1\,,\qquad \gamma_{\phi\phi}=\sin^2\theta\,,\qquad
\gamma_{\theta\phi}=0\,.
$$
We define the Laplace-Beltrami operator on the sphere by
$
\Delta = \gamma^{ij} D_iD_i
$,
where $\gamma^{ij}$ is its inverse and $D_i$ is the associated Levi-Civita connection.
I would like to prove that
$$\Delta \mathbf n = -2 \mathbf n$$
and that in higher dimensions the same holds with $2$ replaced by the dimension of the sphere. I came to believe that this is true by an explicit check in spherical coordinates in dimensions $3$, $4$ and $6$.
Considering that parallel transport of a given tangent vector $\mathbf v$ defined at the point $x+dx$ to the point $x$ is defined by keeping constant its embedding components and then projecting it on the sphere at the point $x$, we have
$$
\mathbf v_{\parallel}(x+dx,x)=\mathbf v(x+dx)-\mathbf v(x+dx)\cdot \mathbf n (x)\, \mathbf n(x)
$$
hence
$$
D_i\mathbf v\, dx^i = \mathbf v_{\parallel}(x+dx,x)- \mathbf v (x)= (\partial_i\mathbf v+\partial_i\mathbf n \cdot \mathbf v\, \mathbf n)dx^i
$$
where we have used $\mathbf n \cdot \partial_i\mathbf v+\partial_i\mathbf n \cdot \mathbf v=0$ and
$$
D_i\mathbf v = \partial_i\mathbf v+\partial_i\mathbf n \cdot \mathbf v\, \mathbf n\,.
$$
Applying this to the basis vectors $\mathbf e_j =\partial_j\mathbf n$ affords
$$
D_i\mathbf e_j = \partial_i \mathbf e_j+\gamma_{ij}\mathbf n\,.
$$
But unfortunataly I am not able to go further.
 A: There is a more general formula. For every immersion $\mathbf x \colon M^n \to \mathbb{E}^m$ of an $n$-dimensional manifold $M$ into $\mathbb{E}^m$ that
$$
  \Delta \mathbf{x} = n H
$$
where $H$ is the mean curvature vector of $\mathbf{x}$. This is sometimes called the formula of Beltrami. For $M=S^n$ the mean curvature is $1$, so one obtains
$$
 \Delta \mathbf{x} = n \mathbf{n}.
$$
Note that this formula differs from yours by a minus sign. I can think of two reasons. 


*

*Some authors put a minus sign in their definition of the Laplacian, some don't. Your definition doesn't have a minus sign, so in the calculation below there doesn't appear a minus sign. 

*Replacing the normal vector $\mathbf{n}$ by $-\mathbf{n}$ also gives a change of sign.


Proof. Let $v$ be an arbitrary vector in $\mathbb{E}^m$ and $p\in M$. If $\{e_1,\ldots, e_n\}$ is an orthonormal basis of $T_p M$, we can extend $e_1,\ldots, e_n$ to an orthonormal frame $E_1,\ldots, E_n$ such that 
$$
   D_{E_i} E_j = 0 \quad \text{at $p$ for $i,j=1,\ldots,n$,}
$$ 
where $D$ is the Levi-Civita connection of $M$. Then at $p$ we have
$$
 \begin{align*}
  (\Delta \langle \mathbf{x},v\rangle)_p &= 
     \sum_{i=1}^n e_i\langle  E_i,v\rangle = \sum_{i=1}^n \langle \bar D_{e_i}E_i,v\rangle \\
 &= \sum_{i=1}^n \langle h(e_i,e_i),v\rangle = n \langle H,v\rangle(p).
 \end{align*}
$$
Here $\bar D$ stands for the Levi-Civita connection on $\mathbb{E}^n$.
Since both $\Delta x$ and $H$ are independent of the choice of local basis, we have $\langle\Delta x, v\rangle = n \langle H,v \rangle$ for any $v$. Since $v$ was arbitrary and the inner product is non-degenerate, the formula of Beltrami follows.
Reference: Pseudo-Riemannian Geometry and Delta Invariants by B.-Y. Chen.
As Yuri Vyatkin pointed out in the comments, your question is related to this question, which was also recently asked.
A: A way to prove the above is the following (this is probably a special case of the more general answer given by @Ernie060, but I still need to fill in a few details). 
In three-dimensional Euclidean space $\mathbb R^3$, the metric in Cartesian coordinates is $\delta_{IJ}=\mathrm{diag}(1,1,1)$ and in spherical coordinates, defined by $\mathbf x = r\,\mathbf n(x^i)$, reads $g_{rr}=1$ and $g_{ij}=r^2\gamma_{ij}$, with $\gamma_{ij}=\partial_i\mathbf n\cdot\partial_j\mathbf n$.  Comparing the two expressions for the Laplacian in the given coordinate systems, we have
$$
0=\Delta_{\mathbb R^{3}}\mathbf x=\frac{1}{r^2}\partial_r(r^2 \mathbf n)+\frac{1}{r}\Delta_{S^2}\mathbf n
$$
which yields precisely
$$
\Delta_{S^2}\mathbf n = -2\mathbf n\,.
$$
This is actually a consequence of the fact that any second covariant derivative of $\mathbf x$ vanishes (since it vanishes in the Cartesian coordiante frame), hence in particular
$$
0=\nabla_i \nabla_j \mathbf x =r D_i D_j\mathbf n-\Gamma_{ij}^r\mathbf n\,,
$$
but explicit calculation affords $\Gamma_{ij}^r=-r\gamma_{ij}$ and hence
$$
D_i D_j \mathbf n = - \gamma_{ij} \mathbf n\,.
$$
A: You can also develop in spherical harmonics. This works equally well on spheres of arbitrary dimension. In Cartesian coordinates $x_1, x_2, \ldots, x_d$, the vector $\boldsymbol{n}$ on $\mathbb{S}^{d-1}$ is
$$
\boldsymbol n = (x_1, x_2, \ldots, x_d), $$
and since each entry is, manifestly, a homogeneous harmonic polynomial of degree $1$, $\boldsymbol n$ is a spherical harmonic of degree $1$. Now, it is well-known that a spherical harmonic of degree $\ell$ is an eigenvector of $\Delta_{\mathbb{S}^{d-1}}$ with eigenvalue $-\ell(\ell+d-2)$. Specializing to $\ell=1$ completes the proof that
$$
\Delta_{\mathbb{S}^{d-1}} \boldsymbol n= -(d-1)\boldsymbol n.$$
REMARK. I will prove that the eigenfunctions of the Laplace-Beltrami are precisely the homogeneous harmonic polynomials, as stated above. This illustrates a computational technique that can be often used in practice for the Laplace-Beltrami operator.
Let $H_\ell=H_\ell(x_1, x_2, \ldots, x_d)$ be a homogeneous harmonic polynomial of degree $\ell$. Then, by harmonicity, $\Delta_{\mathbb R^d}H_\ell = 0$. Expanding the Laplacian in polar coordinates yields
$$\tag{*}\Delta_{\mathbb R^n} H_{\ell}=
\frac{1}{r^{d-1}}\partial_r(r^{d-1}\partial_r H_\ell) + \frac{1}{r^2}\Delta_{\mathbb S^{d-1}}H_\ell=0.$$
Since $H_\ell$ is homogeneous, by the Euler theorem on homogeneous functions we have that $x\cdot \nabla H_\ell=\ell H_\ell$, that is
$$
r\partial_r H_\ell = \ell H_\ell.$$
Plugging this into (*), we arrive at
$$
\frac{\ell(\ell+d-2)}{r^2}H_\ell + \frac{\Delta_{\mathbb{S}^{d-1}}H_\ell}{r^2}=0, $$
from which our claim immediately follows.
