Proof for covariant vector transformation law (I'm asking this on the physics exchange not on the math one because i don't need an extremely rigorous explanation)
I understand the derivation for the contravariant vector transformation law is given by
$$a^i = \frac{\partial \xi^i}{\partial x^K}a^K,$$
but I can't seem to find any proof for why the covariant vector transformation law is given by
$$a_i = \frac{\partial x^K}{\partial \xi^i}a_K.$$
The explanation for this that I keep seeing is just "the covariant vector transformation is similar to the contravariant one, but you flip the derivative", which doesn't help. Why is it flipped?
 A: In the more general language of vector fields and differential geometry, this benefits from some careful manipulations.
So the different coordinates you have $x^{1,2,\dots D}$ and $\xi^{1,2,\dots D}$ are scalar fields, maps from your underlying space of points $\mathcal P$ (imagine, say, the 2D surface of Earth) to real numbers $\mathbb R$, usually satisfying some smoothness properties. Then a given other smooth scalar field $f:\mathcal P\to\mathbb R$, you can find out that it is represented by two different smooth functions $\mathbb R^D\to \mathbb R,$ different for each coordinate system:
$$ f(p) = f_x\left(x^1(p),\dots x^D(p)\right)=f_\xi\left(\xi^1(p),\dots \xi^D(p)\right).$$You can then also do this trick for the scalar fields $x^i$ and $\xi^i$, coming up with functions like $x^i_\xi(\xi^1,\dots \xi^D)$ and $\xi^i_x(x^1,\dots x^D)$ which allow you to convert between them.
Maybe the most obvious thing you can say here is that they are inverses in the sense that
$$x^k_\xi\left(\xi^1_x(x^1,\dots x^D),\dots \xi^D_x(x^1,\dots x^D)\right)=x^k.$$As a result, taking a partial with respect to $x^i,$ $$
\tag{1}
\delta^k_i = \sum_{n=1}^D\frac{\partial x^k_\xi}{\partial \xi^n}~\frac{\partial \xi^n_x}{\partial x^i}.$$
The above should not seem technically difficult to you, but it is a different perspective from how some undergraduate courses approach this subject. It may be helpful to re-read it in a quiet corner until you are sure you understand the perspective.
So in this context a vector field $v^i$ is usually identified with a directional derivative, they are sometimes called “derivations.” That they are Leibniz-linear means that this plays really well with the above coordinates-are-just-other-scalar-fields approaches. So you've got a directional derivative $V$, its action on a scalar field $V(f)$ is given by the Leibniz property as
$$
V(f)=\sum_{n=1}^D\frac{\partial f_x}{\partial x^n}~V(x^n).$$(If you want to you can just take this as a definition of what it means for something to be a directional derivative.) So you see that we can write this as $V=\sum_n v_x^n~(\partial/\partial x^n)$ where $v_x^k=V(x^k)$ is just the directional derivative of the underlying scalar field we are using as our $k$th coordinate. And of course to find the expansion of $V$ in the $\xi$-coordinates we have $V=\sum_n v_\xi^n~(\partial/\partial \xi^n)$ and its action on a scalar field in those coordinates is
$$
\begin{align}
V(f)&=\sum_{n=1}^D\frac{\partial f^n_\xi}{\partial \xi^n}~V(\xi^n)\\
&=\sum_{m,n}^D\frac{\partial f^k_\xi}{\partial \xi^n}~\frac{\partial \xi^n_x}{\partial x^m}~V(x^m),\tag 2
\end{align}
$$
and there is your transformation law that you already have, in this new perspective.
How this answers your question
Once you have equations (1) and (2) you technically already have everything. A covector is a linear map from a vector field to a scalar field, physically we sometimes visualize them as a sort of repeating set of 2D planes across a 3D space, they allow you to count how many planes a given vector crosses, from its base to its tip. (This is also why in solid state physics you will be using the dual lattice to understand crystalline planes.) You want these numbers to to be coordinate independent, you have some $u_i~v^i$ expression giving a scalar, you know that the vector side needs to transform to contain $\partial \xi/\partial x$ due to (2), but you know that the scalar field should be invariant. The Kronecker $\delta$ from (1) is kind of the only obvious place to go to get this invariance! In equations,
$$
\sum_i u_i v^i = \sum_{ik} u_k\delta^k_i v^i= \sum_{ikn}u_k~\frac{\partial x^k_\xi}{\partial \xi^n}~\frac{\partial \xi^n_x}{\partial x^i}~v^i,
$$and if we identify the right hand part of this with $v_\xi$ then the left hand side must be $u_\xi$.
But that is a little unprincipled and a principled answer is only about one line longer... A covector $\mathbf f(V)$, being a linear operation on the tangent space, must have some linear expansion in those $x$-coordinates for that tangent space,
$$\mathbf f(V)
= \sum_n f_n(x^1\dots)~V(x^n))\\
=\sum_{m,n} f_n(x^1\dots)~\frac{\partial x_\xi^n}{\partial \xi^m}~V(\xi^m)\\
$$and we are, perhaps surprisingly, immediately done.
So the intuition is that these two things are inverses and you need the inverse to keep the dot product invariant, so if $\tilde V = q(V)$ then you need this inverse aspect $\tilde u(\dots) =u(q^{-1}(\dots)$ to get $u(V) = \tilde u(\tilde V).$ But you can also just view this as saying that $f$ naturally acts in $x$-space so we need the coordinates of $V$ in $x$-space, if we don't have those because we are in $\xi$-space then we need to convert by multiplying by  $\partial x/\partial \xi$.
