Calculating the divergence of the Gravitational field $\nabla \cdot \vec{F}$ I want to calculate the divergence of the Gravitational field: $$\nabla\cdot \vec{F}=\nabla\cdot\left( -\frac{GMm}{\lvert \vec{r} \rvert^2} \hat{e}_r\right )$$ in spherical coordinates.
I know that in spherical coordinates: $$\begin{aligned} & x=r \sin\theta \cos \phi \\&y=r\sin\theta \sin \phi \\& z=r\cos\theta \end{aligned}$$
and the unit vector are:
$$\begin{aligned} & e_r=\begin{pmatrix}\sin\theta\cos\phi\\\sin\theta \sin\phi\\\cos\theta \end{pmatrix} \\ & e_{\theta}=\begin{pmatrix}\cos\theta\cos\phi\\\cos\theta \sin\phi\\-\sin\theta \end{pmatrix}\\&e_{\phi}=\begin{pmatrix}-\sin\phi\\\cos\phi\\0\end{pmatrix}\end{aligned}$$
Now I need to convert my original vector field into spherical coordinates (this is the part I am not really sure about):
$$\vec{F}=-\frac{GMm}{x^2+y^2+z^2} \hat{e}_x-\frac{GMm}{x^2+y^2+z^2} \hat{e}_y-\frac{GMm}{x^2+y^2+z^2} \hat{e}_z $$
transforming the coordinates: $x^2+y^2+z^2=(r\sin\theta\cos\phi)^2+(r\sin\theta\sin\phi)^2+(r\cos\theta)^2=r^2$
$$\implies\vec{F}=\frac{-GMm}{r^2}\left(\hat{e}_x +\hat{e}_y +\hat{e}_z \right )$$
How can I transform the unit vectors now? Do I just replace them by the spherical unit verctors?
Is there a short really cool way to calculate the divergence of this vector field? I know that the answer should be zero except at $r=0$ the divergence should be undefined.
 A: Since there's only $r$ dependence,
\begin{align*}
  \nabla \cdot \mathbf{F} &=
  \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} F_{r}) \\
  &= \frac{1}{r^{2}} \frac{\partial}{\partial r} (-GMm) \\
  &= 0 
\end{align*}
for $\mathbf{r}\in \mathbb{R}^{3} \backslash \{ \mathbf{0} \}$.
May refer to this
A: $\newcommand{\angles}[1]{\left\langle\,{#1}\,\right\rangle}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\Li}[1]{\,\mathrm{Li}}
 \newcommand{\ol}[1]{\overline{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
\color{#f00}{-\nabla\cdot\pars{{1 \over r^{2}}\,{\vec{r} \over r}}} & =
-\nabla\cdot\nabla\pars{-\,{1 \over r}} = \nabla^{2}\pars{1 \over r} = \color{#f00}{-4\pi\, \delta\pars{\vec{r}}}
\\[5mm]
\mbox{because}\
\color{#f00}{\int_{V}
\bracks{-\nabla\cdot\pars{{1 \over r^{2}}\,{\vec{r} \over r}}}\dd^{3}\vec{r}}
& =
-\int_{S}{1 \over r^{2}}\,{\vec{r} \over r}\cdot\dd\vec{\mathsf{S}} =
-\int\dd\Omega_{\vec{r}} = \color{#f00}{-4\pi}\,,\qquad\vec{0} \in V
\end{align}
