I was trying as exercise to derive the Laplacian in spherical coordinates using riemannian geometry tools.
Definition If $\mathcal{M}$ is a smooth manifold of dimension $n$ and $\left\{E_1,\ldots,E_n \right\}$ is an orthonormal frame then we define for $f : \mathcal{M} \to \mathbb{R}$ $$ \text{grad} f = \sum_{i=1}^n E_i(f)E_i $$ and $$ \Delta f = \sum_{i=1}^n E_i(E_i(f)). $$
Let $$ \mathbb{S}^2(r) = \left\{ (x,y,z) \in \mathbb{R}^3 : \left\{ \begin{array}{l} x = r\sin \phi \cos \theta \\ y = r\sin \phi \sin \theta \\ z = r\cos\phi \end{array} \right. \;\;,0\leq \theta \leq 2\pi, 0 \leq \phi \leq \pi \right\} $$
So there's single chart for the sphere, name this $\varphi$ . If I did the calculation right I have
$$ \varphi_{*,(\theta,\phi)} = \begin{pmatrix} -r\sin\phi \sin\theta & r\cos\phi \cos\theta \\ r\sin\phi \cos\theta & r\cos\phi \sin\theta \\ 0 & -r\sin\phi \end{pmatrix} $$
To define an o.n. reference frame I did the following calculations
$$ \begin{array}{l} \left\langle \frac{\partial}{\partial \theta}, \frac{\partial}{\partial \phi} \right\rangle_{\mathbb{S}^2(r),(\theta,\phi)} = \left\langle \varphi_{*,(\theta,\phi)}\frac{\partial}{\partial \theta}, \varphi_{*,(\theta,\phi)} \frac{\partial}{\partial \phi} \right\rangle_{\mathbb{R}^3,\varphi(\theta,\phi)} = 0\\ \left\langle \frac{\partial}{\partial \theta}, \frac{\partial}{\partial \theta} \right\rangle_{\mathbb{S}^2(r),(\theta,\phi)} = \left\langle \varphi_{*,(\theta,\phi)} \frac{\partial}{\partial \theta},\varphi_{*,(\theta,\phi)} \frac{\partial}{\partial \theta} \right\rangle_{\mathbb{R}^3,\varphi(\theta,\phi)} = r^2 \sin \phi \\ \left\langle \frac{\partial}{\partial \phi}, \frac{\partial}{\partial \phi} \right\rangle_{\mathbb{S}^2(r),(\theta,\phi)} = \left\langle \varphi_{*,(\theta,\phi)} \frac{\partial}{\partial \phi}, \varphi_{*,(\theta,\phi)} \frac{\partial}{\partial \phi} \right\rangle_{\mathbb{R}^3,\varphi(\theta,\phi)} = r^2 \end{array} $$
We have then
$$ T_{(\theta,\phi)}\mathbb{S}^2(r) = \text{span} \left\{ \frac{1}{r\sin\phi} \frac{\partial}{\partial \theta}, \frac{1}{r} \frac{\partial}{\partial \phi} \right\} $$ Which is an orthonormal frame for every $(\theta,\phi)$. This gives me the following expression of the laplacian $$ \Delta_{\mathbb{S}^2(r)}f = \left( \frac{1}{r^2 \sin^2\phi} \frac{\partial^2}{\partial \theta^2} + \frac{1}{r^2} \frac{\partial^2}{\partial \phi^2} \right)f $$
According to few references I've seen around the expression is wrong. I did the calculation few times but I cannot spot the error.
Can you maybe tell me what I am doing wrong? I am assuming the metric on $\mathbb{S}^2(r)$ is the one induced by $\mathbb{R}^3$.
Your definition for $\text{grad}(f)$ is fine, however, the expression you gave for $\Delta f$ only works if you fix a point $p\in M$ and consider a 'geodesic frame at $p$', i.e an orthonormal frame $\{E_1,\dots, E_n\}$ defined in an open neighborhood of $p$ such that for all $i,j\in\{1,\dots, n\}$, $(\nabla_{E_i}E_j)(p)=0$. The orthonormal frame you gave doesn't satisfy the 'geodesic' condition, which is why you're not getting the right result.
The 'correct' definition for the Laplacian is as the divergence of the gradient of $f$, where the divergence of a vector field $X$ is defined as $\text{div}(X)(p):=\text{trace}\left(v_p\mapsto\nabla_{v_p}X\right)$, the trace of a linear mapping $T_pM\to T_pM$. Let us now derive the formula the divergence in terms of a geodesic frame at $p$. Say $X=\sum_{j=1}^nX^jE_j$, then \begin{align} (\text{div }X)(p)&=\text{trace}(v_p\mapsto\nabla_{v_p}X)\\ &=\sum_{i=1}^n\left\langle \nabla_{E_i(p)}X, E_i(p)\right\rangle\\ &=\sum_{i,j=1}^n \left\langle \nabla_{E_i(p)}(X^jE_j), E_i(p)\right\rangle\\ &=\sum_{i,j=1}^n \left\langle (E_i(X^j))(p)\cdot E_j(p)+ X^j(p)\nabla_{E_i(p)}E_j, E_i(p)\right\rangle\\ &=\sum_{i=1}^n(E_i(X^i))(p), \end{align} where we first used the product rule on $\nabla_{E_i(p)}(X^jE_j)$, and then used that we have an orthonormal and geodesic frame at $p$. This is the general formula for the divergence of any vector field. If you specialize to $X=\text{grad}(f)$, then $X^i=E_i(f)$, so you recover your stated formula for the Laplacian.
Now, since you're not actually using a geodesic frame in your computation, you won't get the right answer. An alternative approach to computing the Laplacian is to use the Voss-Weyl formula for the divergence of any vector field applied to the gradient: \begin{align} \Delta f&=\frac{1}{\sqrt{|g|}}\frac{\partial}{\partial x^{\mu}}\left(\sqrt{|g|}g^{\mu\nu}\frac{\partial f}{\partial x^{\nu}}\right). \end{align} You already worked out the metric components $g_{\mu\nu}$ in your post, so carrying out this summation should be trivial.