Laplace operator in spherical coordinates, abstract approach I'd like to show the well-known formula of the Laplacian operator for euclidean $\mathbb{R}^3$ in spherical coordinates:
$$
\Delta U = \frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial U}{\partial r}\right) + \frac{1}{r^2 \sin\vartheta}\frac{\partial }{\partial \vartheta}\left(\sin \vartheta \frac{\partial U}{\partial \vartheta}\right)+\frac{1}{r^2 \sin^2\vartheta}\frac{\partial U}{\partial\varphi}.
$$
Where as usual
$$
x=r\sin\vartheta \cos\varphi\\
y=r\sin\vartheta \sin\varphi\\
z=r\cos\vartheta.
$$
To gain something from this proof, I would like to avoid the "mechanical" use of the chain rule, and rather follow a more abstract path.
 A: Since we have a metric, the usual euclidean metric, which is expressed as $g=g_{ij}\tilde{d}x^i\otimes \tilde{d}x^j=\delta_{ij}\tilde{d}x^i\otimes \tilde{d}x^j$ in cartesian coordinates $x,y,z$, we can establish the euclidean preferred volume form to be:
$$
\tilde{\omega }= \tilde{d}x\wedge \tilde{d}y\wedge \tilde{d}z.
$$
Then, in spherical coordinates we have: $g= g_{i'j'}\tilde{d}x^{i'}\otimes \tilde{d}x^{j'}$
$$
g_{i'j'} 
= \left(
      \begin{array}{ccc}
        1 & 0 &   0 \\
        0 & r^2 & 0 \\
        0 & 0 & r^2\sin^2\vartheta
      \end{array} \right)
$$
and therefore
$$
\tilde{\omega} = \big|\text{det} (g_{i'j'})\big|^{1/2}\ \tilde{d}r\wedge \tilde{d}\vartheta \wedge \tilde{d}\varphi = r^2\sin\vartheta \ \tilde{d}r\wedge \tilde{d}\vartheta \wedge \tilde{d}\varphi .
$$
Now we can think of the coordinate-free definition of the laplacian, which is namely:
$$
\Delta U = \text{div}_{\tilde{\omega}} \bar{d}U
$$
where the divergece of a vector field $\bar{X}$ with respect to the given volume form is
$$
\left(\text{div}_{\tilde{\omega}}\bar{X} \right)\tilde{\omega} = \tilde{d}\left[\tilde{\omega}(\bar{X},\ ,\ )\right],
$$
whereas the gradient of a smooth function $U$ is:
$$
\bar{d}U = g^{-1}(\tilde{d}U,\ ).
$$
This tells us we have to do the following calculations: let $\bar{X} = X^r \partial/{\partial r} + X^\vartheta \partial/\partial\vartheta + X^\varphi \partial/\partial\varphi,$ then
$$
\tilde{d}\left[\tilde{\omega}(\bar{X},\ ,\ )\right]=
\tilde{d}\left[r^2\sin\vartheta\left(X^r \tilde{d}\vartheta\wedge \tilde{d}\varphi + X^\vartheta \tilde{d}\varphi \wedge \tilde{d}r+ X^\varphi \tilde{d}r \wedge \tilde{d}\vartheta \right)\right]=\\
\left[\sin\vartheta \frac{\partial}{\partial r}\left(r^2 X^r\right) + r^2 \frac{\partial}{\partial\vartheta}\left( \sin\vartheta X^\vartheta\right) + r^2\sin\vartheta \frac{\partial X^\varphi}{\partial \varphi}\right]\ \tilde{d}r\wedge \tilde{d}\vartheta \wedge \tilde{d}\varphi  = 
\left[\frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 X^r\right) + \frac{1}{\sin\vartheta} \frac{\partial}{\partial\vartheta}\left( \sin\vartheta X^\vartheta\right) + \frac{\partial X^\varphi}{\partial \varphi}\right]\ \tilde{\omega}.
$$
We have shown:
$$
\text{div}_{\tilde{\omega}}\bar{X} = \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 X^r\right) + \frac{1}{\sin\vartheta} \frac{\partial}{\partial\vartheta}\left( \sin\vartheta X^\vartheta\right) + \frac{\partial X^\varphi}{\partial \varphi}.
$$
Now, 
$$
\tilde{d}U = \frac{\partial U}{\partial r}\tilde{d} r+\frac{\partial U}{\partial \vartheta}\tilde{d} \vartheta +\frac{\partial U}{\partial \varphi}\tilde{d} \varphi,
$$
but since the inverse of the metric is
$$
g^{i\ 'j'} =
\left(
      \begin{array}{ccc}
        1 & 0 &   0 \\
        0 & \frac{1}{r^2} & 0 \\
        0 & 0 & \frac{1}{r^2\sin^2\vartheta}
      \end{array} \right)
$$
the gradient is just
$$
\bar{d}U=  \frac{\partial U}{\partial r} \frac{\partial}{\partial r} + \frac{1}{r^2} \frac{\partial U}{\partial \vartheta}\frac{\partial}{\partial \vartheta} + \frac{1}{r^2\sin^2\theta}\frac{\partial U}{\partial \varphi}\frac{\partial}{\partial \varphi}.
$$
Finally we have:
$$
\Delta U = \text{div}_{\tilde{\omega}}\bar{d} U = \frac{1}{r^2}\frac{\partial U}{\partial r}\left(r^2\frac{\partial}{\partial r}\right) + \frac{1}{r^2 \sin\vartheta}\frac{\partial }{\partial \vartheta}\left(\sin \vartheta \frac{\partial U}{\partial \vartheta}\right)+\frac{1}{r^2 \sin^2\vartheta}\frac{\partial U}{\partial\varphi}.
$$
Notice that, while the laplacian is a scalar operator acting on functions and therefore the form of the result depends only on the coordinates we choose,  divergence and gradient vary greatly depending also on what vector basis we select; in our derivation we used the unnormalized coordinate basis associated to spherical coordinates, while the common form of these operators is expressed with respect to normalized non-coordinate basis.
If however we need the expressions of the latter with respect to unit vectors $\hat{r},\ \hat{\vartheta},\ \hat{\varphi}$, we have:
$$
\hat{r} = \frac{\partial}{\partial r}\\
\hat{\vartheta} = \frac{1}{r}\frac{\partial}{\partial \vartheta}\\
\hat{\varphi} = \frac{1}{r\sin\vartheta}\frac{\partial}{\partial \varphi}.
$$
This implies that the components of our vector $\bar{X}$ with respect to this new basis are $X^{\hat{r}} = X^r$, $X^{\hat{\vartheta}}=rX^\vartheta$,  $X^{\hat{\varphi}} = r \sin\theta X^\varphi$, yielding:
$$
\text{div}_{\tilde{\omega}}\bar{X} = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2 X^{\hat{r}}\right) + \frac{1}{r \sin\vartheta}\frac{\partial}{\partial\vartheta}\left(\sin\vartheta X^{\hat{\vartheta}}\right) + \frac{1}{r \sin\vartheta}\frac{\partial X^{\hat{\varphi}}}{\partial \varphi}
$$
$$
\bar{d}U =  \frac{\partial U}{\partial r}\hat{r} +\frac{1}{r}\frac{\partial U}{\partial\vartheta}\hat{\vartheta}+\frac{1}{r \sin\vartheta}\frac{\partial U}{\partial \varphi}\hat{\varphi}.
$$
