Divergence of a radial vector field I am reading Modern Electrodynamics by Zangwill and cannot verify equation (1.61) [page 7]:
\begin{equation}
\nabla \cdot  \textbf{g}(r)=\textbf{g}^{\prime}\cdot \mathbf{\hat{r}},
\end{equation}
where the vector field $\textbf{g}(r)$ is only nonzero in the radial direction.
By using the divergence formula in Spherical coordinates, I get:
\begin{align}
\nabla \cdot  \textbf{g}(r)&=\frac{1}{r^2} \partial_r (r^2 g_r) + \frac{1}{r \sin \theta} \partial_{\theta} (g_{\theta} \sin \theta) + \frac{1}{r \sin \theta} \partial_{\phi} g_{\phi}\\
&=\frac{2}{r}g_r + \frac{d}{dr}g_r\\
&= \frac{2}{r}\textbf{g}\cdot \mathbf{\hat{r}}+\textbf{g}^{\prime}\cdot \mathbf{\hat{r}}
\end{align}
What is going wrong?
 A: The function $\mathbf g(r)$ is not a radial vector field; it is a vector field which depends only on $r\equiv\sqrt{x^2+y^2+z^2}$.  Indeed unless the vector field is the trivial $\mathbf g(r) = \vec 0$, then it cannot be radial; you can see this by noting that if $\mathbf g$ is radial, then the vectors at the north and south poles must point in opposite directions, but if $\mathbf g=\mathbf g(r)$ then the vectors at the north and south poles must be the same.
A: In $\S\texttt{1.3 Derivatives }$ of the textbook $''$Modern Electrodynamics$''$ by A.Zangwill the function $\,\mathbf g\left(r\right)\,$ is considered as a vector function of the position radius $\,r$. Following the notation therein this means that
\begin{equation}
\mathbf g\left(r\right)\boldsymbol= \mathrm g_x\left(r\right)\mathbf{\hat{x}}\boldsymbol+\mathrm g_y\left(r\right)\mathbf{\hat{y}}\boldsymbol+\mathrm g_z\left(r\right)\mathbf{\hat{z}} 
\boldsymbol=
\begin{bmatrix}
\mathrm g_x\left(r\right) & \mathrm g_y\left(r\right)  & \mathrm g_z\left(r\right)
\end{bmatrix}
\begin{bmatrix}
\mathbf{\hat{x}}\vphantom{\dfrac{a}{b}}\\
\mathbf{\hat{y}}\vphantom{\dfrac{a}{b}}\\ 
\mathbf{\hat{z}}\vphantom{\dfrac{a}{b}}
\end{bmatrix}
\tag{01}\label{01}  
\end{equation}
where $\,\mathbf{\hat{x}},\mathbf{\hat{y}}\,$ and $\,\mathbf{\hat{z}}\,$  the unit vectors along the cartesian coordinates axes $\,x,y\,$ and $\,z\,$ respectively.
Now, for the analysis of this vector function by components in the spherical coordinate system we make use of the relation between the local spherical coordinate system $\,\{\mathbf{\hat{r}},\boldsymbol{\hat{\theta}},\boldsymbol{\hat{\phi}}\}\,$ and the global cartesian $\,\{\mathbf{\hat{x}},\mathbf{\hat{y}},\mathbf{\hat{z}}\}\,$
\begin{align}
\mathbf{\hat{x}} & \boldsymbol=\left(\sin\theta\cos\phi\right)\mathbf{\hat{r}}\boldsymbol+\left(\cos\theta\cos\phi\right)\boldsymbol{\hat{\theta}}\boldsymbol+\left(\boldsymbol-\sin\phi\right)\boldsymbol{\hat{\phi}}  
\tag{02a}\label{02a}\\
\mathbf{\hat{y}} & \boldsymbol=\left(\sin\theta\sin\phi\right)\mathbf{\hat{r}}\boldsymbol+\left(\cos\theta\sin\phi\right)\boldsymbol{\hat{\theta}}\boldsymbol+\left(\cos\phi\right)\boldsymbol{\hat{\phi}} 
\tag{02b}\label{02b}\\
\mathbf{\hat{z}} & \boldsymbol=\left(\cos\theta\right)\mathbf{\hat{r}}\boldsymbol+\left(\boldsymbol-\sin\theta\right)\boldsymbol{\hat{\theta}}
\tag{02c}\label{02c}
\end{align}
or in matrix form
\begin{equation}
\begin{bmatrix}
\mathbf{\hat{x}}\vphantom{\dfrac{a}{b}}\\
\mathbf{\hat{y}}\vphantom{\dfrac{a}{b}}\\ 
\mathbf{\hat{z}}\vphantom{\dfrac{a}{b}}
\end{bmatrix}
\boldsymbol=
\begin{bmatrix}
\sin\theta\cos\phi & \cos\theta\cos\phi & \boldsymbol-\sin\phi \vphantom{\dfrac{a}{b}}\\
\sin\theta\sin\phi & \cos\theta\sin\phi & \hphantom{\boldsymbol-}\cos\phi  \vphantom{\dfrac{a}{b}}\\ 
\cos\theta & \boldsymbol-\sin\theta & \hphantom{\boldsymbol-}0 \vphantom{\dfrac{a}{b}}
\end{bmatrix}
\begin{bmatrix}
\mathbf{\hat{r}}\vphantom{\dfrac{a}{b}}\\
\boldsymbol{\hat{\theta}}\vphantom{\dfrac{a}{b}}\\ 
\boldsymbol{\hat{\phi}}\vphantom{\dfrac{a}{b}}
\end{bmatrix}
\tag{03}\label{03}  
\end{equation}
Inserting this expression in equation \eqref{01} we have
\begin{equation}
\mathbf g\left(r\right)\boldsymbol= \mathrm g_r\left(r,\theta,\phi\right)\mathbf{\hat{r}}\boldsymbol+ \mathrm g_\theta\left(r,\theta,\phi\right)\boldsymbol{\hat{\theta}}\boldsymbol+\mathrm g_\phi\left(r,\theta,\phi\right)\boldsymbol{\hat{\phi}} 
\tag{04}\label{04}  
\end{equation}
where
\begin{align}
\mathrm g_r\left(r,\theta,\phi\right) & \boldsymbol=\sin\theta\cos\phi\cdot\mathrm g_x\left(r\right)\boldsymbol+\sin\theta\sin\phi\cdot\mathrm g_y\left(r\right)\boldsymbol+\cos\theta\cdot\mathrm g_z\left(r\right)
\tag{05a}\label{05a}\\
\mathrm g_\theta\left(r,\theta,\phi\right) & \boldsymbol=\cos\theta\cos\phi\cdot\mathrm g_x\left(r\right)\boldsymbol+\cos\theta\sin\phi\cdot\mathrm g_y\left(r\right)\boldsymbol-\sin\theta\cdot\mathrm g_z\left(r\right)
\tag{05b}\label{05b}\\
\mathrm g_\phi\left(r,\theta,\phi\right) & \boldsymbol=\boldsymbol-\sin\phi\cdot\mathrm g_x\left(r\right)\boldsymbol+\cos\phi\cdot\mathrm g_y\left(r\right)
\tag{05c}\label{05c}
\end{align}
So all components $\,\mathrm g_r,\mathrm g_\theta,\mathrm g_\phi\,$ are functions of $\,r,\theta,\phi\,$ and especially $\,\mathrm g_\theta,\mathrm g_\phi\,$ are not zero in general.
It's not convenient to find the divergence of such a vector field using spherical coordinates. To the contrary, by cartesian coordinates
\begin{equation}
\begin{split}
\boldsymbol\nabla\boldsymbol\cdot\mathbf g\left(r\right) & \boldsymbol= \dfrac{\partial\mathrm g_x\left(r\right)}{\partial x}\boldsymbol+\dfrac{\partial\mathrm g_y\left(r\right)}{\partial y}\boldsymbol+\dfrac{\partial\mathrm g_z\left(r\right)}{\partial z}\\ 
&\boldsymbol=
\dfrac{\mathrm d\mathrm g_x\left(r\right)}{\mathrm dr}\dfrac{\partial r}{\partial x}\boldsymbol+\dfrac{\mathrm d\mathrm g_y\left(r\right)}{\mathrm dr}\dfrac{\partial r}{\partial y}\boldsymbol+\dfrac{\mathrm d\mathrm g_z\left(r\right)}{\mathrm dr}\dfrac{\partial r}{\partial z}\\
&\boldsymbol=
\dfrac{\mathrm d\mathrm g_x\left(r\right)}{\mathrm dr}\dfrac{x}{r}\boldsymbol+\dfrac{\mathrm d\mathrm g_y\left(r\right)}{\mathrm dr}\dfrac{y}{r}\boldsymbol+\dfrac{\mathrm d\mathrm g_z\left(r\right)}{\mathrm dr}\dfrac{z}{r}\\
&\boldsymbol=
\dfrac{\mathrm d\mathbf g\left(r\right)}{\mathrm dr}\boldsymbol\cdot\dfrac{\mathbf r}{r}\boldsymbol=\mathbf g'\boldsymbol\cdot\mathbf{\hat{r}}\\
\end{split}
\tag{06}\label{06}  
\end{equation}
$=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!=\!$
ADDENDUM : $\texttt{ Gymnastics with Spherical Coordinates}$
If we insist to follow the difficult path of finding the divergence using spherical coordinates then
\begin{equation}
\boldsymbol\nabla\boldsymbol\cdot\mathbf g\left(r\right) 
\boldsymbol=\underbrace{\dfrac{1}{r^2}\dfrac{\partial\left(r^2\mathrm g_r\right)}{\partial r}}_{\boxed{\boldsymbol 1}}\boldsymbol+\underbrace{\dfrac{1}{r\sin\theta}\dfrac{\partial\left(\sin\theta\,\mathrm g_\theta\right)}{\partial \theta}}_{\boxed{\boldsymbol 2}}\boldsymbol+\underbrace{\dfrac{1}{r\sin\theta}\dfrac{\partial\mathrm g_\phi}{\partial \phi}}_{\boxed{\boldsymbol 3}}
\tag{07}\label{07}  
\end{equation}
Using the expressions \eqref{05a},\eqref{05b} and \eqref{05c} for the components $\,\mathrm g_r,\mathrm g_\theta\,$ and $\,\mathrm g_\phi\,$ respectively we have in sequence
\begin{equation}
\begin{split}
\boxed{\boldsymbol 1} & \boldsymbol= \dfrac{1}{r^2}\dfrac{\partial\left(r^2\mathrm g_r\right)}{\partial r}\boldsymbol=\dfrac{2}{\,r\,}\mathrm g_r\boldsymbol+\dfrac{\partial\mathrm g_r}{\partial r}\quad \boldsymbol\implies\\
\boxed{\boldsymbol 1} & \boldsymbol= \dfrac{2\sin\theta\cos\phi}{r}\mathrm g_x\boldsymbol+\dfrac{2\sin\theta\sin\phi}{r}\mathrm g_y\boldsymbol+\dfrac{2\cos\theta}{r}\mathrm g_z\\ 
& \hphantom{\boldsymbol=}\boldsymbol+\sin\theta\cos\phi\dfrac{\mathrm d\mathrm g_x}{\mathrm dr}\boldsymbol+\sin\theta\sin\phi\dfrac{\mathrm d\mathrm g_y}{\mathrm dr}\boldsymbol+\cos\theta\dfrac{\mathrm d\mathrm g_z}{\mathrm dr}\\
\end{split}
\tag{08}\label{08}  
\end{equation}
and
\begin{equation}
\begin{split}
\boxed{\boldsymbol 2} & \boldsymbol= \dfrac{1}{r\sin\theta}\dfrac{\partial\left(\sin\theta\,\mathrm g_\theta\right)}{\partial \theta}\boldsymbol=\dfrac{\cot\theta }{r}\mathrm g_\theta\boldsymbol+\dfrac{1}{r}\dfrac{\partial\mathrm g_\theta}{\partial \theta}\quad \boldsymbol\implies\\
\boxed{\boldsymbol 2} & \boldsymbol=\dfrac{\cos^2\theta\cos\phi }{r\sin\theta}\mathrm g_x\boldsymbol+\dfrac{\cos^2\theta\sin\phi }{r\sin\theta}\mathrm g_y\boldsymbol-\dfrac{\cos\theta }{r}\mathrm g_z\\ 
 & \hphantom{\boldsymbol=}\boldsymbol-\dfrac{\sin\theta\cos\phi}{r}\mathrm g_x\boldsymbol-\dfrac{\sin\theta\sin\phi}{r}\mathrm g_y\boldsymbol-\dfrac{\cos\theta}{r}\mathrm g_z\\
\end{split}
\tag{09}\label{09}  
\end{equation}
and
\begin{equation}
\begin{split}
\boxed{\boldsymbol 3} & \boldsymbol= \dfrac{1}{r\sin\theta}\dfrac{\partial\mathrm g_\phi}{\partial \phi}\quad \boldsymbol\implies\\
\boxed{\boldsymbol 3} & \boldsymbol=\boldsymbol-\dfrac{\cos\phi}{r\sin\theta}\mathrm g_x\boldsymbol-\dfrac{\sin\phi}{r\sin\theta}\mathrm g_y\\
\end{split}
\tag{10}\label{10}   
\end{equation}
So
\begin{equation}
\begin{split}
\boldsymbol\nabla\boldsymbol\cdot\mathbf g\left(r\right) 
&\boldsymbol=\boxed{\boldsymbol 1} \boldsymbol+\boxed{\boldsymbol 2} \boldsymbol+\boxed{\boldsymbol 3}\\  &\boldsymbol=\overbrace{\left(\dfrac{2\sin\theta\cos\phi}{r}\boldsymbol+\dfrac{\cos^2\theta\cos\phi }{r\sin\theta}\boldsymbol-\dfrac{\sin\theta\cos\phi}{r}\boldsymbol-\dfrac{\cos\phi}{r\sin\theta}\right)}^{0}\mathrm g_x\\ 
& \hphantom{\boldsymbol=}\boldsymbol+\overbrace{\left(\dfrac{2\sin\theta\sin\phi}{r}\boldsymbol+\dfrac{\cos^2\theta\sin\phi }{r\sin\theta}\boldsymbol-\dfrac{\sin\theta\sin\phi}{r}\boldsymbol-\dfrac{\sin\phi}{r\sin\theta}\right)}^{0}\mathrm g_y\\
& \hphantom{\boldsymbol=}\boldsymbol+\overbrace{\left(\dfrac{2\cos\theta}{r}\boldsymbol-\dfrac{\cos\theta}{r}\boldsymbol-\dfrac{\cos\theta}{r}\right)}^{0}\mathrm g_z\\
& \hphantom{\boldsymbol=}\boldsymbol+\underbrace{\sin\theta\cos\phi}_{x/r}\dfrac{\mathrm d\mathrm g_x}{\mathrm dr}\boldsymbol+\underbrace{\sin\theta\sin\phi}_{y/r}\dfrac{\mathrm d\mathrm g_y}{\mathrm dr}\boldsymbol+\underbrace{\cos\theta}_{z/r}\dfrac{\mathrm d\mathrm g_z}{\mathrm dr}\\
&\boldsymbol=\dfrac{\mathrm d\mathbf g\left(r\right)}{\mathrm dr}\boldsymbol\cdot\dfrac{\mathbf r}{r}\boldsymbol=\mathbf g'\boldsymbol\cdot\mathbf{\hat{r}}\\
\end{split}
\tag{11}\label{11}  
\end{equation}
qed.
A: The first equation is presumptive written without rigorous thinking, and it is not right.  For a radial vector field, $\vec g(r) = g_r(r) \hat r $, its divergence is:
\begin{align}
\nabla \cdot  \vec g(r) &= \vec \nabla \cdot   g_r(r) \hat r \\
&= \vec \nabla g_r(r) \cdot \hat r + g_r(r) \vec \nabla \cdot \hat r \\
&= g'_r(r) +  g_r(r) \vec \nabla \cdot \frac{\vec r}{r}\\
&= g'_r(r) +  g_r(r) \left(\frac{\vec \nabla \cdot \vec r}{r} - \frac{\vec r \cdot \vec \nabla r}{r^2}  \right)\\
&= g'_r(r) + g_r(r) \left(\frac{3}{r} - \frac{r}{r^2}  \right) \\
&= g'_r(r) + g_r(r) \frac{2}{r}\\
& = \vec g'(r) \cdot \hat r + \frac{2}{r} \vec g \cdot \vec r.
\end{align}
The result is same as your second equation.
A: $
\renewcommand\vec\mathbf
\newcommand\R{\mathbb R}
\newcommand\adj\overline
\newcommand\diff\underline
$
I want to add that the desired result is a straightforward application of an extremely general chain rule for $\nabla$. For any $f : \R^n \to R$ where $R = \R$ or $R = \R^n$
$$
  \nabla = \adj f_{\vec x}(\partial_f).
$$
This requires some unpacking.

*

*This rewriting of $\nabla$ can be applied in any expression where $\nabla$ appears in a linear slot, i.e. where if you replaced $\nabla$ with a vector $\vec v$ then the resulting expression is linear in $\vec v$.

*$\nabla = \nabla_{\vec x}$ is implicitly differentiating with respect to $\vec x$.

*$\partial_f = \partial/\partial f$ when $R = \R$ and $\partial_f = \nabla_f$ when $R = \R^n$; in either case, what this notation means is that for an appropriate function $g$
$$
   \partial_fg(f(x)) = [\partial_yg(y)]_{y=f(x)},
 $$
or similar expressions for divergence, curl, etc.

*$\adj f_{\vec x} : R \to \R^n$ is the adjoint of the total differential $\diff{f_{\vec x}} : \R^n \to \R$ at $\vec x$ under a suitable inner product, either multiplication for $R = \R$ or the standard inner product for $R = \R^n$. The total differential $\diff{f_{\vec x}}$ at $\vec x$ is the best linear approximation of $f$ at that point; the matrix of $v \mapsto \diff{f_{\vec x}}(v)$ in the standard basis is the Jacobian. When $R = \R$, we can write
$$
 \adj f_{\vec x}(v) = \nabla(vf(\vec x)),
   $$
and when $R = \R^n$ we can write
$$
 \adj f_{\vec x}(\vec v) = \nabla(\vec v\cdot f(\vec x)).
   $$
In your specific case, $f(\vec x) = r = |\vec x|$, and the chain rule becomes
$$\begin{aligned}
  \nabla\cdot\vec g(r)
&= \nabla\cdot\vec g(f(\vec r))
= \adj f_{\vec x}\left(\frac\partial{\partial f}\right)\cdot\vec g(f)
= (\nabla|\vec x|\frac\partial{\partial f})\cdot\vec g(f)
\\
&= (\nabla|\vec x|)\cdot\frac{\partial\vec g}{\partial f}
= \frac{\vec x}{|\vec x|}\cdot\vec g'(f(\vec x))
\\
&= \hat{\vec r}\cdot\vec g'(r).
\end{aligned}$$
