In this question @barto explains the definition of $Df|_p$ in differential geometry.
This question is about the formal definition of $D^2f|_p$, given a function $f:M\to \mathbb R$ for some manifold $M$.
Given that $Df_p$ gives us a mapping from the tangent space of $f$ at $p\in M$ to $\mathbb R$, we can state that $D$ is a map $(M\to \mathbb R)\times M \to (T(M)\to \mathbb R)$, where $T(M)$ is the union of tangent spaces on $M$ at all points.
(I am not sure if the way I defined this is correct now.)
My question is, is $D^2$ defined rigorously as an operator? By $D^2$ I mean the operator that takes $f:M\to\mathbb R$, $p\in M$, such that $D^2f|_p=D(Df|_p)$. Basically, $D^2$ should map $f$ and $p$ to a linear operator that is represented by the Hessian matrix.
This answer became overly long in an attempt to give some motivation for the constructions done, generalizing the familiar case in $\mathbb{R}^n$, but for a TL;DR see the last paragraph. You might refer, for instance, to Lee's "Riemannian Manifolds" for more info on these things.
I think that the appropriate notion to generalize is the Hessian matrix as a bilinear form. We work first in $\mathbb{R}^n$. For vectors $X$ and $Y$, we can take $H_f(X,Y) = Y^TH_f X = Y^iX^j\frac{\partial^2f}{\partial x^i\partial x^j}$, where $H_f$ is the Hessian matrix for some function $f$. We want some sort of interpretation of what this means, in order to generalize it properly. It looks something like a second derivative along the directions of $X$ and $Y$, but the problem with this interpretation is that $Df_p(X)$ is just a number, not a function which we can differentiate again. However, if we extend $X$ to a vector field $X(p)$ on $\mathbb{R}^n$ having constant components (the same as $X$), $Xf(p)=Df_p(X(p))=X^i\frac{\partial f}{\partial x^i}(p)$ is a function on $\mathbb{R}^n$, which we can again differentiate in the direction $Y$. We then have $D(Xf)_p(Y) = H_f(X,Y)$. By symmetry we might equally well extend $Y$ to a global vector field, and see $D(Yf)_p(X)=H_f(Y,X)=H_f(X,Y)$.
Now what is the generalization of this construction to any manifold $M$? A vector field being constant (having constant components) does not make sense in this setting, since this is coordinate-dependent. However, it makes sense for vector field to be covariantly constant - or rather to have covariant derivative zero. By this covariant derivative we mean that coming from the Levi-Civita connection $\nabla$ of some riemannian metric $g$ on $M$. We then say that $X$ has covariant derivative zero at $p$ if $\nabla_YX(p)=0$ for any $Y\in T_pM$.
What we want to do, then, is to take any two tangent vectors $X$ and $Y$ at a point $p\in M$, i.e. in $T_pM$, extend them to vector fields on some neighborhood of $p$ with covariant derivatives zero at $p$, and define $H_f(X,Y) = Y(Xf)$. Now we might ask if this construction makes $H_f$ well-defined, i.e. if this depends on the extensions chosen. First we show that this is in fact symmetric, i.e. that $Y(Xf)=X(Yf)$ whenever $Y$ and $X$ have covariant derivative zero at $p$. We use the fact that the Levi-Civita connection is torsion-free to write \begin{equation} Y(Xf)(p)-X(Yf)(p) = [Y,X]f(p) = (\nabla_YX - \nabla_XY)f(p) = ((0-0)f)(p)=0, \end{equation} since the vector fields have covariant derivative zero at $p$. Then if $X_1$ and $X_2$ are different extensions, we have \begin{equation} H_f(X_1,Y)-H_f(X_2,Y) = X_1(Yf) - X_2(Yf) = (X_1-X_2)(Yf) \end{equation} and this is zero at $p$ since $X_1$ and $X_2$ take the same value there. Symmetrically, $H_f(X,Y)$ does not depend on the extension of $Y$.
The object we seek in order to define your $D^2$ rigorously is the covariant derivative $\nabla$. Aside from acting on vector fields, this also acts functions simply as a differential: $\nabla f=df$. It also acts on a one-form $\omega$, mapping it to a two-tensor as $\nabla\omega(X,Y) = Y(\omega(X)) - \omega(\nabla_YX)$. Applying it twice to a function $f$ gives the covariant Hessian $\nabla^2f$, with \begin{equation} \nabla^2f(X,Y)=Y(\nabla f(X)) - \nabla f(\nabla_YX) = Y(Xf) - (\nabla_YX)f. \end{equation} Again, since the Levi-Civita connection is torsion-free, we have $Y(Xf) - (\nabla_YX)f = X(Yf) - (\nabla_XY)f = \nabla^2f(Y,X)$, so the covariant Hessian is symmetric. With these two expressions we can show in the same way as above that $\nabla^2f(X,Y)(p)$ depends only on the values of $X$ and $Y$ at $p$. If $X$ and $Y$ are the extensions used in defining $H_f$, we get that $\nabla^2f(X,Y) = Y(Xf)= X(Yf) = H_f(X,Y)$.