The Hessian matrix $\{\partial_i \partial_j f \}$ of a function $f:\mathbb{R}^n \to \mathbb{R}$ depends on the coordinate system you choose. If $x_1,\cdots,x_n$ and $y_1,\cdots,y_n$ are two sets of coordinates (say, in some open neighborhood of a manifold), then $\frac{\partial f(y(x))}{\partial x_i} = \sum_{k} \frac{\partial f}{\partial y_k} \frac{\partial y_k}{\partial x_i}$. Differentiating again, this time with respect to $x_j$, we get $\frac{\partial^2 f(y(x))}{\partial x_i \partial x_j} = \sum_{k} \sum_{l} \frac{\partial^2 f}{\partial y_k \partial y_l} \frac{\partial y_l}{\partial x_j} \frac{\partial y_k}{\partial x_i}+\frac{\partial f(y(x))}{\partial y_k}\frac{\partial^2y}{\partial x_i \partial x_j}$. At a critical point, the second term goes away, so we will consider such a case.

In other words, if the derivative is a differential $1$-form, i.e. $\sum_{i} \frac{\partial f}{\partial x_i} dx_i$, a section of the cotangent bundle, then the second derivative should be $\sum_{k,l} \frac{\partial^2 f(y(x))}{\partial y_k \partial x_l} dy_k \otimes dy_l$. This makes sense since $dy_k=\sum_{i} \frac{\partial y_k}{\partial x_i} dx_i$, and $dy_l=\sum_{j} \frac{\partial y_l}{\partial x_j} dx_i$, meaning that $\sum_{k,l} \frac{\partial^2 f(y(x))}{\partial y_k \partial x_l} dy_k \otimes dy_l = \sum_{k,l} \frac{\partial^2 f(y(x))}{\partial y_k \partial x_l} (\sum_{i} \frac{\partial y_k}{\partial x_i} dx_i) \otimes (\sum_{j} \frac{\partial y_l}{\partial x_j} dx_j) = \sum_{i,j,k,l} \frac{\partial^2 f(y(x))}{\partial y_k \partial x_l} \frac{\partial y_k}{\partial x_i} \frac{\partial y_l}{\partial x_j} dx_i dx_j = \sum_{i,j} \frac{\partial^2 f}{\partial x_i \partial x_j}$, making it coordinate independent. Note that I did not use exterior powers, I used tensor powers, since I wanted to actually find a way to make sense of second derivatives, rather than having $d^2=0$. This means the Hessian should be a rank $2$ tensor ((2,0) or (0,2), I can't remember which, but definitely not (1,1)).

Does this make sense? Can we then express the third, etc, derivative as a tensor? More interestingly, how can this help us make sense of Taylor's formula? Can we come up with a coordinate-free Taylor series of a function at a point on a manifold?

EDIT: An in general, if the first $n$ derivatives vanish, then the $n+1$ derivative should be a rank $n+1$ tensor, right?

Natural Operations in Differential Geometryby Kolar, Michor, and Slovak. (My apologies for not typing in the proper diacritics on the letters.) – Willie Wong Jul 15 '10 at 10:32