How are we able to identify the differential with a covector field? I have been going through Dr. Frederic Schuller's Lectures on the Geometric Anatomy of Theoretical Physics, in these lectures, he states that the gradient of a function is the covector field
$$df: M\to T^*M$$ that sends a point $p\in M$ to the covector $$d_pf: T_pM\to\mathcal R$$ (where $\mathcal R$ is the real numbers) defined by $$d_pf(X):= X(f).$$  In Lee's book Introduction to Smooth Manifolds, he states that the differential (or push-forward) of a smooth map $F:M\to N$ is $$d_pF: T_pM\to T_{F(p)}N$$ that sends $X\in T_pM$ to $d_pF(X) \in T_{F(p)}N$ where the last function acts on $f \in C^{\infty}(N)$ to produce $$d_pF(X)(f) := X(f \circ F).$$ Lee says that we can use the differential of a smooth map $f:M\to \mathcal R$ and the fact that $$T_p\mathcal R \cong \mathcal R$$ in order to define the gradient as a covector field (Chapter 11: The differential of a function). How is it that we can do this? Why do we call the gradient a covector if it is really just a special case of the differential?
 A: The differential of a real-valued function is a covector by definition. I'll work with real-valued functions on the plane for simplicity of notation but this works with functions on manifolds. If $f:\Bbb R^2 \to \Bbb R$ is smooth, then its partial derivatives exist and at any point $p$ the differential $df_p:\Bbb R^2 \to \Bbb R $ is a linear map which best approximates the real variation of $f$ near $p$. (Here technically the domain and codomain of the differential are tangent spaces at $p$ and $f(p)$, which we can harmlessly identify with Euclidean spaces.) So if $v = (h,k)$ is a tangent vector at $p$ (picture a displacement away from the base point $p$), the differential $df_p(v) \approx f(p + v) - f(p)$ where the error goes to $0$ faster than $||v||$ does. For calculations using partial derivatives, one can easily show that $$df_p(v) = df_p(h,k) = f_x(p)h + f_y(p)k = \begin{pmatrix} f_x(p) & f_y(p)\end{pmatrix}\begin{pmatrix} h \\ k\end{pmatrix},$$ i.e., the matrix associated to the linear map or covector $df_p$ is the row of partial derivatives at $p$. It acts on displacements or column vectors via matrix multiplication. (More generally, the matrix of $df_p$ will be the Jacobian matrix of the given $f$ between two manifolds.) One can also write $df = f_x dx + f_y dy$ globally, where $x,y$ are the canonical projection maps.
The differential $df$, a covector field, is thus a natural object obtained directly from the smooth function $f$. On the other hand, the gradient $\nabla f$, a vector field, requires additional structure to define, namely one needs an inner product (Riemannian metric on manifolds). In Euclidean space, it is obvious that the formula above for computing $df_p(h,k)$ might as well have been formulated using vectors instead:
$$df_p(h,k) = f_x(p)h + f_y(p)k  = \begin{pmatrix} f_x(p) \\ f_y(p)\end{pmatrix} \cdot\begin{pmatrix} h \\ k\end{pmatrix} = \nabla f_p \cdot (h,k).$$
This is the formula that defines the gradient vector on a general manifold: if you have a smooth inner product $g_p$ on each tangent space, you define $\nabla f$ via the formula $$df_p(v) = g_p(\nabla f_p,v) \quad \forall v.$$
Abstractly, the non-degenerate bilinear form $g_p$ gives you a natural "musical" isomorphism $T_p M \to T_p^* M$ mapping $v$ to $g_p(v,-)$, so by "raising or lowering the indices using the metric tensor" as physicists would say we can treat vectors and covectors alike. One needs to be careful when the metric is non-Euclidean though: the formula for $\nabla f$ will generally not be so simple as a vector of partial derivatives.
For a simple example, consider $\Bbb R^2 - \{ 0 \}$ with polar coordinates $(r,\theta)$; then for any smooth function $f:\Bbb R^2 - \{ 0 \} \to \Bbb R$ we have the usual natural formula $$df = f_r dr + f_\theta d\theta,$$ but
$$\nabla f = f_r \partial_r + \frac1{r^2} f_\theta \partial_\theta$$
(note the extra $\frac 1{r^2}$ factor coming from the Euclidean metric in polar coordinates).
