Why does the gradient vector change with coordinate transformation? I know that a vector is invariant under coordinate transformation, that is:
$$\vec{v} = v^i e_i = \tilde{v}^i\tilde{e_i}$$
where here, I'm implying einstein's summation convention.
When I transform the gradient vector however, it does not follow that rule. For example:
$$\vec\nabla f = e_i g^{ij} \frac{\partial f}{\partial x^j} = e_r\frac{\partial f}{\partial r} + e_\theta\frac{1}{r}\frac{\partial f}{\partial \theta} \neq g^{rr}e_r\frac{\partial f}{\partial r} + g^{\theta \theta}e_\theta\frac{\partial f}{\partial \theta}$$
Why is this? Is this because I'm not using the covariant derivative? I'm not well versed in using it but I suspect the problem is because the partial derivatives do not transform like a vector. In other words, is this correct for all my choice of basis?:
$$\vec \nabla f = e_i g^{ij}D_jf$$
 A: $
\newcommand\R{\mathbb R}
\newcommand\PD[2]{\frac{\partial#1}{\partial#2}}
$You're using the wrong basis vectors; you have to use those that are directly associated with the coordinates.
We need to distinguish between points, vectors, and coordinates; let $M$ be our space of points, $V$ our space of vectors, and $C$ our space of coordinates. All of these happen to be modeled as $\R^n$, but that does not mean they are literally the same thing. If our manifold $M$ were not flat, then instead of $V$ we would need a tangent space $T_pM$ for each $p \in M$. A choice of coordinate space $C$ consists of two things: a point function $P : C \to M$ taking coordinates $(x^1,\dotsc,x^n)$ to the corresponding point $P(x^1,\dotsc, x^n)$, and coordinate functions $X^i : M \to C$ taking a point $p$ to its $i^{\text{th}}$ coordinate $X^i(p) = x^i$. These are inverses of each other:
$$
  P(X^1(p), \dotsc, X^n(p)) = p,\quad X^i(P(x^1,\dotsc, x^n)) = x^i.
$$
The coordinate basis of $V$ is then defined as
$$
  e_i = \PD{P(x^1,\dotsc,x^n)}{x^i}.
$$
Derivatives of $M$-valued functions are always valued in $V$. Note that $e_i$ is actually a function $C \to V$, but we will usually not write this coordinate dependence explicitly. Using the coordinate functions $X^i$, these can also be viewed as function $M \to V$.
If we also have a metric on $V$, in our case the usual dot product, then there is a reciprocal basis $e^i$ uniquely defined such that $e^i\cdot e_j = \delta^i_j$. Since $v\cdot\nabla$ ought to be the derivative in the $v$ direction, we see that
$$
  e_i\cdot\nabla X^j = \PD{X^j}{x^i} = \delta_i^j
$$
so we must have $\nabla X^j = e^j$. This gives us the coordinate expression for $\nabla$:
$$
  \nabla = \sum_{i=1}^ne^i\PD{}{x^i}.
$$
Now let's apply this to the specific case where $n = 2$, $M$ is modeled with Cartesian coordinates $C_{xy}$ so that $C_{xy} = M$, and we are interested in polar coordinates $C_{r\theta}$. The point function $P : C_{r\theta} \to M$ is
$$
  P(r, \theta) = (r\cos\theta, r\sin\theta),
$$
and the coordinate functions $R : M \to C_{r\theta}$ and $\Theta : M \to C_{r\theta}$ are
$$
  R(x, y) = \sqrt{x^2 + y^2},\quad \Theta(x, y) = \arctan\frac yx.
$$
(Or something like that, I am not being particularly careful.) The coordinate basis is
$$
  e_r = (\cos\theta)e_x + (\sin\theta)e_y,\quad e_\theta = -(r\sin\theta)e_x + (r\cos\theta)e_y.
$$
Notice that $e_\theta$ is not a unit vector. If $\hat e_\theta = e_\theta/|e_\theta|$ then
$$
  \hat e_\theta = e_\theta/r = -(\sin\theta)e_x + (\cos\theta)e_y.
$$
We can then compute the reciprocal basis to be
$$
  e^r = e_r,\quad e^\theta = e_\theta/r^2 = \hat e_\theta/r.
$$
Hence the gradient expressed in polar coordinates is
$$
  \nabla = e^r\PD{}r + e^\theta\PD{}\theta = e_r\PD{}r + e_\theta\frac1{r^2}\PD{}\theta = e_r\PD{}r + \hat e_\theta\frac1r\PD{}\theta.
$$

The short of it is: $\{e_r, \hat e_\theta\}$ is the wrong basis to transform; it needs to be $\{e_r, e_\theta\} = \{e_r, r\hat e_\theta\}$.
