[proofplan] We compute the Euclidean Laplacian in the spherical coordinate chart by expressing the gradient and divergence in the orthonormal moving frame determined by the coordinate map. The metric scale factors are $1$, $r$, and $r\sin\theta$, so the radial contribution becomes the divergence term $r^{-2}\partial_r(r^2\partial_r\tilde{u})$. The remaining two terms involve only $\theta$ and $\phi$, and they combine exactly into $r^{-2}\Delta_{S^2}\tilde{u}$. [/proofplan] [step:Introduce the spherical coordinate map and its scale factors] Define the spherical coordinate map \begin{align*} \Phi: (0,\infty)\times(0,\pi)\times\mathbb{R} &\to \mathbb{R}^3\setminus\{0\} \end{align*} by \begin{align*} \Phi(r,\theta,\phi) = (r\sin\theta\cos\phi, r\sin\theta\sin\phi, r\cos\theta). \end{align*} The coordinate derivatives are \begin{align*} \partial_r\Phi = (\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta), \end{align*} \begin{align*} \partial_\theta\Phi = (r\cos\theta\cos\phi,r\cos\theta\sin\phi,-r\sin\theta), \end{align*} and \begin{align*} \partial_\phi\Phi = (-r\sin\theta\sin\phi,r\sin\theta\cos\phi,0). \end{align*} Their Euclidean lengths are \begin{align*} |\partial_r\Phi|=1,\qquad |\partial_\theta\Phi|=r,\qquad |\partial_\phi\Phi|=r\sin\theta, \end{align*} and their Euclidean pairwise inner products are zero. Thus the spherical coordinates are orthogonal with scale factors \begin{align*} h_r=1,\qquad h_\theta=r,\qquad h_\phi=r\sin\theta. \end{align*} [/step] [step:Write the gradient in the spherical orthonormal frame] Define the orthonormal frame fields \begin{align*} e_r(r,\theta,\phi)&:=\partial_r\Phi(r,\theta,\phi), \end{align*} \begin{align*} e_\theta(r,\theta,\phi)&:=\frac{1}{r}\partial_\theta\Phi(r,\theta,\phi), \end{align*} and \begin{align*} e_\phi(r,\theta,\phi)&:=\frac{1}{r\sin\theta}\partial_\phi\Phi(r,\theta,\phi). \end{align*} Let $v := \tilde{u}=u\circ\Phi$. Since $u\in C^2(\mathbb{R}^3\setminus\{0\})$ and $\Phi$ is smooth on $(0,\infty)\times(0,\pi)\times\mathbb{R}$, the function $v$ is $C^2$ on this coordinate domain. For a coordinate curve, the chain rule gives \begin{align*} \partial_r v=(\nabla u\circ\Phi)\cdot \partial_r\Phi, \end{align*} \begin{align*} \partial_\theta v=(\nabla u\circ\Phi)\cdot \partial_\theta\Phi, \end{align*} and \begin{align*} \partial_\phi v=(\nabla u\circ\Phi)\cdot \partial_\phi\Phi. \end{align*} Using the scale factors and the orthonormal frame, these identities imply \begin{align*} \nabla u\circ\Phi = (\partial_r v)e_r + \frac{1}{r}(\partial_\theta v)e_\theta + \frac{1}{r\sin\theta}(\partial_\phi v)e_\phi. \end{align*} [guided] The goal is to express $\nabla u$ in directions adapted to spheres centered at the origin. The coordinate map \begin{align*} \Phi: (0,\infty)\times(0,\pi)\times\mathbb{R} &\to \mathbb{R}^3\setminus\{0\} \end{align*} sends $(r,\theta,\phi)$ to \begin{align*} (r\sin\theta\cos\phi, r\sin\theta\sin\phi, r\cos\theta). \end{align*} The function being differentiated in coordinates is \begin{align*} v: (0,\infty)\times(0,\pi)\times\mathbb{R} &\to \mathbb{R}, \end{align*} defined by $v=u\circ\Phi$. Since $u$ is $C^2$ away from the origin and $\Phi$ is smooth on its coordinate domain, $v$ is $C^2$. The coordinate vector $\partial_r\Phi$ has length $1$, while $\partial_\theta\Phi$ has length $r$ and $\partial_\phi\Phi$ has length $r\sin\theta$. Therefore the corresponding unit vectors are \begin{align*} e_r=\partial_r\Phi,\qquad e_\theta=\frac{1}{r}\partial_\theta\Phi,\qquad e_\phi=\frac{1}{r\sin\theta}\partial_\phi\Phi. \end{align*} These vectors are pairwise orthonormal. The chain rule says that differentiating $v=u\circ\Phi$ in a coordinate direction computes the directional derivative of $u$ in the corresponding coordinate vector: \begin{align*} \partial_r v=(\nabla u\circ\Phi)\cdot \partial_r\Phi, \end{align*} \begin{align*} \partial_\theta v=(\nabla u\circ\Phi)\cdot \partial_\theta\Phi, \end{align*} and \begin{align*} \partial_\phi v=(\nabla u\circ\Phi)\cdot \partial_\phi\Phi. \end{align*} Because $\partial_\theta\Phi=r e_\theta$ and $\partial_\phi\Phi=r\sin\theta\, e_\phi$, the components of $\nabla u\circ\Phi$ in this orthonormal frame are \begin{align*} (\nabla u\circ\Phi)\cdot e_r=\partial_r v, \end{align*} \begin{align*} (\nabla u\circ\Phi)\cdot e_\theta=\frac{1}{r}\partial_\theta v, \end{align*} and \begin{align*} (\nabla u\circ\Phi)\cdot e_\phi=\frac{1}{r\sin\theta}\partial_\phi v. \end{align*} Hence \begin{align*} \nabla u\circ\Phi = (\partial_r v)e_r + \frac{1}{r}(\partial_\theta v)e_\theta + \frac{1}{r\sin\theta}(\partial_\phi v)e_\phi. \end{align*} [/guided] [/step] [step:Apply the divergence formula in the same orthogonal coordinates] Let \begin{align*} A: \mathbb{R}^3\setminus\{0\} &\to \mathbb{R}^3 \end{align*} be a $C^1$ vector field whose spherical-coordinate components in the frame $(e_r,e_\theta,e_\phi)$ are \begin{align*} A\circ\Phi=A_r e_r+A_\theta e_\theta+A_\phi e_\phi, \end{align*} where \begin{align*} A_r,A_\theta,A_\phi:(0,\infty)\times(0,\pi)\times\mathbb{R}\to\mathbb{R} \end{align*} are $C^1$. For orthogonal coordinates with scale factors $h_r=1$, $h_\theta=r$, and $h_\phi=r\sin\theta$, the Euclidean divergence is \begin{align*} (\operatorname{div}A)\circ\Phi = \frac{1}{r^2\sin\theta} \frac{\partial}{\partial r}(r^2\sin\theta\,A_r) + \frac{1}{r^2\sin\theta} \frac{\partial}{\partial\theta}(r\sin\theta\,A_\theta) + \frac{1}{r^2\sin\theta} \frac{\partial}{\partial\phi}(rA_\phi). \end{align*} Apply this to $A=\nabla u$. From the preceding step, \begin{align*} A_r=\partial_r v,\qquad A_\theta=\frac{1}{r}\partial_\theta v,\qquad A_\phi=\frac{1}{r\sin\theta}\partial_\phi v. \end{align*} Therefore \begin{align*} (\Delta u)\circ\Phi = \frac{1}{r^2\sin\theta} \frac{\partial}{\partial r}(r^2\sin\theta\,\partial_r v) + \frac{1}{r^2\sin\theta} \frac{\partial}{\partial\theta}(\sin\theta\,\partial_\theta v) + \frac{1}{r^2\sin\theta} \frac{\partial}{\partial\phi}\left(\frac{1}{\sin\theta}\partial_\phi v\right). \end{align*} Since $\sin\theta$ is independent of $r$ and $\phi$, this becomes \begin{align*} (\Delta u)\circ\Phi = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2\partial_r v) + \frac{1}{r^2\sin\theta} \frac{\partial}{\partial\theta}(\sin\theta\,\partial_\theta v) + \frac{1}{r^2\sin^2\theta}\partial_{\phi\phi}v. \end{align*} [/step] [step:Identify the angular part with the spherical Laplacian] Define the angular differential operator \begin{align*} \Delta_{S^2}v = \frac{1}{\sin\theta} \frac{\partial}{\partial\theta}(\sin\theta\,\partial_\theta v) + \frac{1}{\sin^2\theta}\partial_{\phi\phi}v. \end{align*} Substituting this definition into the formula obtained above gives \begin{align*} (\Delta u)\circ\Phi = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2\partial_r v) + \frac{1}{r^2}\Delta_{S^2}v. \end{align*} Since $v=\tilde{u}=u\circ\Phi$, this is exactly the claimed decomposition at every $(r,\theta,\phi)\in(0,\infty)\times(0,\pi)\times\mathbb{R}$. This completes the proof. [/step]