Derive vector gradient in spherical coordinates from first principles

It follows from the general definition of the gradient as $$ langlenabla f(p)|vrangle=d_pf(v)=sum_ileft.frac{partial f}{partial x^i}right|_pdx^i(v) $$ where p is a point in space and v a vector in the tangent space. The summation is over the basis vectors of the tangent space. You can try to expand this expression to get the final result for the component $i$ $$ (nabla f)_i=frac{1}{h_i}frac{partial f}{partial x^i} $$ This is the most useful formula. The quantity $h_i$ is the modulus of the $i$th tangent vector.

Example: you want to compute the gradient in spherical coordinates. The basis of the tangent space is ${frac{partial}{partial r},frac{partial}{partial theta},frac{partial}{partialphi}}$. Since $$ begin{split} left|frac{partial}{partial theta}right|^2&= left| frac{partial x}{partialtheta}frac{partial}{partial x}+ frac{partial y}{partialtheta}frac{partial}{partial y}+ frac{partial z}{partialtheta}frac{partial}{partial z} right|^2 &= r^2cos^2thetacos^2phiunderbrace{left|frac{partial}{partial x}right|^2}_{=1}+ r^2cos^2thetasin^2phiunderbrace{left|frac{partial}{partial y}right|^2}_{=1}+ r^2sin^2thetaunderbrace{left|frac{partial}{partial z}right|^2}_{=1} &=r^2 end{split} $$ Thus we get $$h_theta=left|frac{partial}{partial theta}right|=r$$ In the same spirit you can calculate that $$ h_r=1quadtext{and}quad h_phi=rsintheta $$ giving us the gradient in spherical coordinates $$ nabla f=frac{partial f}{partial r}hat e_r+frac{1}{r}frac{partial f}{partialtheta}hat e_theta+frac{1}{rsintheta}frac{partial f}{partialphi}hat e_phi $$

Proof for the first step

Expand the vector $|nabla frangle$ in terms of basis vectors $$ |nabla frangle =sum_i(nabla f)_i|e_irangle =sum_i(nabla f)_ifrac{1}{h_i}|frac{partial}{partial x^i}rangle $$ This is basically where the factor $h_i$ comes from. Now take $v=|frac{partial}{partial x^k}rangle$ and insert it in the first expression given above. Note that by definition of a dual vector we get $dx^i(|frac{partial}{partial x^k}rangle)=delta_k^i$. The left-hand side is $$ begin{split} langle f|frac{partial}{partial x^k}rangle &=sum_i(nabla f)_ifrac{1}{h_i}langlefrac{partial}{partial x^i}|frac{partial}{partial x^k}rangle &=sum_i(nabla f)_ifrac{1}{h_i}h_i^2delta_{ik} &=(nabla f)_kh_k end{split} $$ Whreas the right-hand side $$ sum_ileft.frac{partial f}{partial x^i}right|_pdx^ileft(|frac{partial}{partial x^k}rangleright) =sum_ileft.frac{partial f}{partial x^i}right|_pdelta^i_k =frac{partial f}{partial x^k} $$ By comparing both expressions you obtain the claim.

We take: $$x=rsinthetacosphi$$ $$y=rsinthetasinphi$$ $$z=rcostheta$$

Now, you know the definition of the gradient in Cartesian coordinates: $vec{nabla}=frac{partial}{partial x}hat{x}+frac{partial}{partial y}hat{y}+frac{partial}{partial z}hat{z}$

Now, we use the chain rule or each component. For instance, $$frac{partial}{partial x}=frac{partial r}{partial x}frac{partial}{partial r}+frac{partial theta}{partial x}frac{partial}{partial theta}+frac{partial phi}{partial x}frac{partial}{partial phi}$$

After lots of cumbersome algebra, this will give you the correct form.

You asked for a proof from "first principles". So let's do it. I'll highlight the most common sources of errors and I'll show an alternative proof later that doesn't require any knowledge of tensor calculus or Einstein notation.

The hard way

First, the coordinates convention:

$$(r,theta,phi) rightarrow (x,y,z) = (rsinthetacosphi,;rsinthetasinphi,;rcostheta)$$

The same way we can express $(x,y,z)$ as $x,mathbf{hat e}_x + y,mathbf{hat e}_y + z,mathbf{hat e}_z$, we can also express $(r,theta,phi)$ as $r',mathbf{hat e}_r + theta',mathbf{hat e}_theta + phi',mathbf{hat e}_phi$, but now the coefficients are not the same: $(r',theta',phi') neq (r,theta,phi)$, in general. This is because spherical coordinates are curvilinear, so the basis vectors are not the same at all points. For small variations, however, they are very similar. More precisely, relative to a point $vec{mathbf p}_0 = (x,y,z)$, a neighbor point $vec{mathbf p}_1 = (x+Delta x,;y+Delta y,;z+Delta z)$ can be described by $Delta vec{mathbf p} = (Delta x,Delta y,Delta z)$ and, in spherical coordinates, if this variation is "infinitesimal", then $dvec{mathbf p} = (dr, dtheta, dphi) = dr,mathbf{hat e}_r + dtheta,mathbf{hat e}_theta + dphi,mathbf{hat e}_phi$. This is basically the motivation for defining the (unnormalized) basis as:

$$ vec{mathbf e}_r = frac{partialvec{mathbf p}}{partial r}, quad vec{mathbf e}_theta = frac{partialvec{mathbf p}}{partial theta}, quad vec{mathbf e}_phi = frac{partialvec{mathbf p}}{partial phi} $$

But this is not normalized yet. Coincidentally, $||partialvec{mathbf p}/partial r||$ turns out to be $1$, but $||partialvec{mathbf p}/partial theta|| = r$, as we'll see. So the actual basis should be defined as:

$$ hat{mathbf e}_r = frac{vec{mathbf e}_r}{||vec{mathbf e}_r||}, quad hat{mathbf e}_theta = frac{vec{mathbf e}_theta}{||vec{mathbf e}_theta||}, quad hat{mathbf e}_phi = frac{vec{mathbf e}_phi}{||vec{mathbf e}_phi||} $$

Explicitly:

$$begin{align} vec{mathbf e}_r &=;( & sinthetacosphi&, & sinthetasinphi&, & costheta) vec{mathbf e}_theta &=;( & rcosthetacosphi&, & rcosthetasinphi&, & -rsintheta) vec{mathbf e}_phi &=;( & -rsinthetasinphi&, & rsinthetacosphi&, & 0) end{align}$$

$$begin{align} ||vec{mathbf e}_r||^2 &= sin^2theta(cos^2phi + sin^2phi) + cos^2theta & &= 1 ||vec{mathbf e}_theta||^2 &= r^2cos^2theta(cos^2phi + sin^2phi) + r^2sintheta & &= r^2 ||vec{mathbf e}_phi||^2 &= r^2sin^2theta(sin^2phi + cos^2phi) & &= r^2sin^2theta end{align}$$

$$begin{align} hat{mathbf e}_r &= vec{mathbf e}_r & &= & &( & sinthetacosphi&, & sinthetasinphi&, & costheta) hat{mathbf e}_theta &= vec{mathbf e}_theta/r & &= & &( & costhetacosphi&, & costhetasinphi&, & -sintheta) hat{mathbf e}_phi &= vec{mathbf e}_phi/(rsintheta) & &= & &( & -sinphi&, & cosphi&, & 0) end{align}$$

You may verify that this also forms an orthogonal basis (hence orthonormal). For example:

$$begin{align} hat{mathbf e}_r cdot hat{mathbf e}_theta &= sinthetacosthetacos^2phi + sinthetacosthetasin^2phi - sinthetacostheta &= 0 end{align}$$

That does not need to happen in general.

To go from one set of coordinates to the other using the basis vectors, we solve:

$$ begin{bmatrix}hat{mathbf e}_r hat{mathbf e}_theta hat{mathbf e}_phiend{bmatrix} = begin{bmatrix} sinthetacosphi & sinthetasinphi & costheta costhetacosphi & costhetasinphi & -sintheta -sinphi & cosphi & 0 end{bmatrix} begin{bmatrix}hat{mathbf e}_x hat{mathbf e}_y hat{mathbf e}_zend{bmatrix} $$

for $hat{mathbf e}_x$, $hat{mathbf e}_y$, and $hat{mathbf e}_z$ in terms of $hat{mathbf e}_r$, $hat{mathbf e}_theta$, and $hat{mathbf e}_phi$. Then any vector $vec{mathbf p} = x,mathbf{hat e}_x + y,mathbf{hat e}_y + z,mathbf{hat e}_z$ can be written in the form $r',mathbf{hat e}_r + theta',mathbf{hat e}_theta + phi',mathbf{hat e}_phi$ by simple substitution. Since this particular basis is orthonormal, there's an alternative way: simply use the dot product. For example, to get $r'$:

$$begin{align} vec{mathbf p}cdotmathbf{hat e}_r &= r',mathbf{hat e}_rcdotmathbf{hat e}_r + theta',mathbf{hat e}_thetacdotmathbf{hat e}_r + phi',mathbf{hat e}_phicdotmathbf{hat e}_r &= r' end{align}$$

Now to the gradient. Using matrix notation, we can write the gradient as a row vector and the formula for the chain rule becomes:

$$begin{align} vecnabla f &= begin{bmatrix}frac{partial f}{partial x} & frac{partial f}{partial y} & frac{partial f}{partial z}end{bmatrix} &= begin{bmatrix}frac{partial f}{partial r} & frac{partial f}{partial theta} & frac{partial f}{partial phi}end{bmatrix} begin{bmatrix} frac{partial r}{partial x} & frac{partial r}{partial y} & frac{partial r}{partial z} frac{partial theta}{partial x} & frac{partial theta}{partial y} & frac{partial theta}{partial z} frac{partial phi}{partial x} & frac{partial phi}{partial y} & frac{partial phi}{partial z} end{bmatrix} end{align}$$

Call the matrix on the right $J$ (it's the Jacobian matrix). Note that this also works the other way around too:

$$ begin{bmatrix}frac{partial f}{partial r} & frac{partial f}{partial theta} & frac{partial f}{partial phi}end{bmatrix} = begin{bmatrix}frac{partial f}{partial x} & frac{partial f}{partial y} & frac{partial f}{partial z}end{bmatrix} begin{bmatrix} frac{partial x}{partial r} & frac{partial y}{partial r} & frac{partial z}{partial r} frac{partial x}{partial theta} & frac{partial y}{partial theta} & frac{partial z}{partial theta} frac{partial x}{partial phi} & frac{partial y}{partial phi} & frac{partial z}{partial phi} end{bmatrix} $$

And call this other matrix $J'$. We can invert the first equation to get $vecnabla f ,J^{-1} = vecnabla f ,J'$ $Rightarrow$ $vecnabla f ,left(J^{-1}-J'right) = 0$. Since this works for an arbitrary $f$, we have $J^{-1} - J' = 0$ $Rightarrow$ $J' = J^{-1}$. An important consequence is that, in general:

$$ frac{partial a}{partial b} neq left(frac{partial b}{partial a}right)^{-1} $$

It appears that the OP made this mistake in a comment, confusing $partial r/partial x$ with $(partial x/partial r)^{-1} = 1/(sinthetacosphi)$, as would be the case if we were using regular (instead of partial) derivatives.

Now we have two ways of calculating the matrix $J$. Directly or by calculating $J'$ first and then inverting it. Let's do it directly. We're gonna need the expressions for $r$, $theta$, and $phi$ in terms of $x$, $y$, and $z$ (for other coordinate systems this might be very difficult to obtain):

$$begin{align} r &= sqrt{x^2+y^2+z^2} theta &= arctanleft(frac{sqrt{x^2+y^2}}{z}right) phi &= arctanleft(frac{y}{x}right) end{align}$$

The partial derivatives are:

$$begin{align} frac{partial r}{partial x} &= frac{x}{sqrt{x^2+y^2+z^2}} & &= sinthetacosphi frac{partial r}{partial y} &= frac{y}{sqrt{x^2+y^2+z^2}} & &= sinthetasinphi frac{partial r}{partial z} &= frac{z}{sqrt{x^2+y^2+z^2}} & &= costheta end{align}$$

$$begin{align} frac{partial theta}{partial x} &= frac{zx}{sqrt{x^2+y^2}left(x^2+y^2+z^2right)} & &= frac{costhetacosphi}{r} frac{partial theta}{partial y} &= frac{zy}{sqrt{x^2+y^2}left(x^2+y^2+z^2right)} & &= frac{costhetasinphi}{r} frac{partial theta}{partial z} &= -frac{sqrt{x^2+y^2}}{x^2+y^2+z^2} & &= -frac{sinphi}{r} end{align}$$

$$begin{align} frac{partial phi}{partial x} &= -frac{y}{x^2+y^2} & &= -frac{sinphi}{rsintheta} frac{partial phi}{partial y} &= frac{x}{x^2+y^2} & &= frac{cosphi}{rsintheta} frac{partial phi}{partial z} &= 0 & &= 0 end{align}$$

Our Jacobian is then:

$$ J = begin{bmatrix} sinthetacosphi & sinthetasinphi & costheta frac{costhetacosphi}{r} & frac{costhetasinphi}{r} & -frac{sinphi}{r} -frac{sinphi}{rsintheta} & frac{cosphi}{rsintheta} & 0 end{bmatrix} $$

Alternatively, we could have calculated the inverse Jacobian (which is straightforward) and then inverted it (which is a nightmare). We can use Wolfram Alpha to confirm that it gives the same result:

Finally, we use the dot product to find the coefficients $r'$, $theta'$, and $phi'$:

$$ r' = vecnabla fcdotmathbf{hat e}_r = begin{bmatrix}frac{partial f}{partial r} & frac{partial f}{partial theta} & frac{partial f}{partial phi}end{bmatrix} J begin{bmatrix}sinthetacosphi sinthetasinphi costhetaend{bmatrix} = begin{bmatrix}frac{partial f}{partial r} & frac{partial f}{partial theta} & frac{partial f}{partial phi}end{bmatrix} begin{bmatrix}1 0 0end{bmatrix} = frac{partial f}{partial r} $$

$$ theta' = vecnabla fcdotmathbf{hat e}_theta = begin{bmatrix}frac{partial f}{partial r} & frac{partial f}{partial theta} & frac{partial f}{partial phi}end{bmatrix} J begin{bmatrix}costhetacosphi costhetasinphi -sinthetaend{bmatrix} = begin{bmatrix}frac{partial f}{partial r} & frac{partial f}{partial theta} & frac{partial f}{partial phi}end{bmatrix} begin{bmatrix}0 1/r 0end{bmatrix} = frac{1}{r}frac{partial f}{partial theta} $$

$$ phi' = vecnabla fcdotmathbf{hat e}_phi = begin{bmatrix}frac{partial f}{partial r} & frac{partial f}{partial theta} & frac{partial f}{partial phi}end{bmatrix} J begin{bmatrix}-sinphi cosphi 0end{bmatrix} = begin{bmatrix}frac{partial f}{partial r} & frac{partial f}{partial theta} & frac{partial f}{partial phi}end{bmatrix} begin{bmatrix}0 0 1/(rsintheta)end{bmatrix} = frac{1}{rsintheta}frac{partial f}{partial phi} $$

Therefore:

$$ vecnabla f = frac{partial f}{partial r}mathbf{hat e}_r + frac{1}{r}frac{partial f}{partial theta}mathbf{hat e}_theta + frac{1}{rsintheta}frac{partial f}{partial phi}mathbf{hat e}_phi $$

---

A much better way

We're gonna need a new notation to avoid having to use different letters for $x$, $y$, and $z$, for example. Instead let's use indices from $1$ to $3$. For Cartesian coordinates we'll use the letter $x$, and for spherical coordinates we'll use the letter $r$. The following should be self explanatory:

$$ vec{mathbf p} = sum_i x_imathbf{hat x}^i = sum_k r_kmathbf{hat r}^k $$

From the definition of the basis vectors:

$$ mathbf{vec r}^k = frac{partial vec{mathbf p}}{partial r_k}, quad mathbf{hat r}^k = frac{mathbf{vec r}^k}{||mathbf{vec r}^k||} = frac{1}{h_k}frac{partial vec{mathbf p}}{partial r_k} $$

Where $h_k triangleq ||mathbf{vec r}^k||$. Expanding in the $x$ basis:

$$ mathbf{hat r}^k = sum_jfrac{1}{h_k}frac{partial x_j}{partial r_k}mathbf{hat x}^j $$

Now the gradient is just:

$$ vecnabla f = mathbf{vec F} = sum_i F_imathbf{hat x}^i = sum_i frac{partial f}{partial x_i}mathbf{hat x}^i $$

To get the $k$'th component in spherical coordinates ($F'_k$), use the dot product:

$$begin{align} F'_k &= mathbf{vec F} cdot mathbf{hat r}^k &= left(sum_i frac{partial f}{partial x_i}mathbf{hat x}^iright) cdot left(sum_jfrac{1}{h_k}frac{partial x_j}{partial r_k}mathbf{hat x}^jright) &= frac{1}{h_k}sum_i frac{partial f}{partial x_i}frac{partial x_i}{partial r_k} &= frac{1}{h_k}frac{partial f}{partial r_k} end{align}$$

and we're done.