Characterizing Hessians among symmetric bilinear tensors I apologize in advance if this is somewhat elementary, but:

Let $(M,g)$ be a compact Riemannian manifold. Is there a "characterization" of which symmetric bilinear tensors $B\in Sym^2(M)$ are Hessians of functions $f\colon M\to \mathbb R$?
In other words, for which $B\in Sym^2(M)$ there exists $f\colon M\to\mathbb R$ such that $B(X,Y)=Hess \ f(X,Y)$, where $Hess\ f(X,Y)=g(\nabla_X \nabla f,Y)=X(Y(f))-\mathrm df(\nabla_X Y)$?

Obviously, since $M$ is assumed compact, any continuous function will attain a min and a max, so a symmetric bilinear form that is a candidate to Hessian cannot be definite. This maybe pretty naive, but I cannot think of any other property that distinguishes Hessians among bilinear symmetric tensors in this context... For example, it doesn't seem to be the case that $B$ has to be parallel, or satisfy any other sort of similar "nice" conditions.
Note that even if we impose the extra condition that $B$ is the Hessian of a Morse function, I don't quite see if the topology of $M$ will impose any properties on $B$ via Morse theory, since $B$ does not know if it is at a (future) critical point or not, right?
 A: This is more of a comment/thinking aloud than an answer, but perhaps the symplectic viewpoint helps:
Think of the Hessian with respect to a Riemannian metric as follows: through the Legendre transfom, the Levi-Civita connection can be seen as a splitting of the tangent bundle of the cotangent bundle of $M$ into vertical and horizontal bundles: $T(T^*M) = V(T^*M) \oplus H(T^*M)$. At each point $p_x \in T^*M$, the subspace $V_{p_x}(T^*M) = V_{p_x}$ and $H_{p_x}(T^*M) = H_{p_x}$ are Lagrangian subspaces.
If $f$ is a function on $M$, the graph of its differential is a Lagrangian submanifold $L \subset T^*M$. The tangent space of $L$ at the point $df(x)$ is a Lagrangian subspace of $T_{df(x)} T^*M$ which is transverse to the vertical space $V_{df(x)}$ and, therefore, we can think of $T_{df(x)}L$ as the graph of a linear transformation from $H_{df(x)}$ to $V_{df(x)}$. Modulo some basic identifications, this linear map is the Hessian of $f$ at $x$.
A nice thing about this description is that it works on Finsler spaces since the splitting $T(T^*M) = V(T^*M) \oplus H(T^*M)$ just comes from the linearlization of the geodesic flow and the convexity of the Hamiltonian. 
If you start with a Hessian, this description suggests the following procedure to construct a function $f$: given a point $p_x \in T^*M$, use the Hessian of $f$ at $x$, the identification
of $V_{p_x}$ and $H_{p_x}$ with $T_xM$ and the decomposition $T_{p_x}T^*M$ to construct a Lagrangian subspace in  $T_{p_x}T^*M$. One ends up with a field of Lagrangian subspaces in
$TT^*M$. We have to find an integral manifold $L$. Such a manifold will be Lagrangian by construction and, given that the  constructed Lagrangian subspaces are transverse to the vertical subspaces, it will also project in a locally-diffeomorphic way onto $M$. It's a good candidate for the graph of $df$.
A: There are local conditions, but they typically involve the curvature tensor of the underlying metric.  For example, if the metric is flat, so that one can choose orthonormal coordinates $x_i$ in which the covariant derivatives are the ordinary derivatives, then the condition that a quadratic form $H = h_{ij}\ dx_idx_j$ be a Hessian is just that
$$
\frac{\partial h_{ij}}{\partial x_k} = \frac{\partial h_{ik}}{\partial x_j} \ .
$$
This is an overdetermined system, and, if it is satisfied, there are solutions $f$ of 
$$
\frac{\partial^2 f}{\partial x_i\partial x_j} = h_{ij}
$$
and they are uniquely determined (locally) up to the addition of a function linear in the $x_i$.
In the more general case, one has local conditions in terms of the curvature of $g$, so the answer is more complicated (and more interesting).  I calculate using moving frames (see below at the end for the 'Nomizu-style' interpretation), so I would explain it this way:  
Let $g = {\omega_1}^2 +\cdots + {\omega_n}^2$ be the expression of $g$ in a local orthonormal coframing $\omega$.  There will exist unique $1$-forms $\theta_{ij}=-\theta_{ji}$ so that (using the summation convention) $d\omega_i = -\theta_{ij}\wedge\omega_j$.  (These are the first structure equations.)
If $f$ is a function defined in the domain of the coframing, one will have 
$$
df = f_i\ \omega_i
\qquad\text{and}\qquad
df_i = -\theta_{ij}\ f_j + f_{ij}\ \omega_j
$$
and, by definition, $\mathrm{Hess}_g(f) = f_{ij}\ \omega_i{\circ}\omega_j$.  
To determine whether a given symmetric quadratic form $H = h_{ij}\ \omega_i{\circ}\omega_j$ can be written as $H = \mathrm{Hess}_g(f)$ for some $f$, one computes the exterior derivative of the equations $df_i = -\theta_{ij}\ f_j + h_{ij}\ \omega_j$ 
and finds that one must have
$$
h_{ikl}-h_{ilk} = -R_{ijkl}\ f_j\ , \tag{1}
$$
where $dh_{ij} = -\theta_{ik}\ h_{kj} + \theta_{kj}\ h_{ik} + h_{ijk}\ \omega_k$ and the $R_{ijkl}=-R_{ijlk}$ are the components of the Riemann curvature tensor in this coframing, as defined by the second structure equations
$$
d\theta_{ij} = -\theta_{ik}\wedge\theta_{kj} + \tfrac12\ R_{ijkl}\ \omega_k\wedge\omega_l\ .
$$
The system (1) gives necessary conditions for $f$ to exist, since, in most cases, it completely determines the only possible functions $f_i$ that could be its derivatives in this coframing.  
It may be worth pausing to interpret (1) as a global equation.  In Nomizu-style notation, one has $\mathrm{Hess}_g(f) = (\nabla\nabla f)^{\flat\flat}$, but I am going to ignore the distinction between $T=TM$ and $T^\ast=T^*M$ since we have a metric $g$ that is $\nabla$-parallel and just write $\mathrm{Hess}_g(f) = \nabla\nabla f$.  Applying the Bianchi identities, one has 
$$
\sigma\bigl(\nabla\bigl(\mathrm{Hess}_g(f)\bigr)\bigr) = \mathsf{R}(\nabla f),
$$
where $\sigma:\mathsf{S}^2(T^\ast)\otimes T^\ast\to T^\ast\otimes\Lambda^2(T^\ast)$ is the canonical skew-symmetrization operator, and $\mathsf{R}:T\to T^\ast\otimes\Lambda^2(T^\ast)$ is the usual mapping induced by the curvature operator.  In particular, if $H = \mathrm{Hess}_g(f)$, then one must have 
$$
\sigma\bigl(\nabla H\bigr) = \mathsf{R}(\nabla f)\tag{1'}
$$
which, in my moving frames notation, is the system (1).  
Now, the image of $\sigma$ is, as usual, the kernel of the further skew-symmetrization map $\sigma: T^\ast\otimes\Lambda^2(T^\ast)\to \Lambda^3(T^\ast)$ and, fortunately, the first Bianchi identity implies that $\mathsf{R}$ also takes values in this kernel, which has dimension $n\cdot \tfrac12 n(n{-}1) - \tfrac16 n(n{-}1)(n{-}2)= \tfrac13 n(n^2{-}1)$.  
We already know what happens when $\mathsf{R}\equiv0$, namely, the first order system of equations $\sigma\bigl(\nabla H\bigr)=0$ are necessary and sufficient that $H$ be locally expressible as a Hessian.  
However, the 'generic' situation is that $\mathsf{R}$ is injective, so let's consider that case.  (I'll say a bit more about what happens in the intermediate cases at the end.)  When $\mathsf{R}$ is injective, let $\mathsf{R}(T)\subset T^\ast\otimes\Lambda^2(T^\ast)$ be the image subbundle of rank $n$.  The equations $(1')$ then say that $H$ must satisfy the system of $\tfrac13 n(n^2{-}1) - n = \tfrac13 n(n^2{-}4)$ first order equations
$$
\sigma(\nabla H) \equiv 0\ \text{modulo}\ \mathsf{R}(T).\tag{1''}
$$
Assuming that these equations do hold, then there is a unique vector field $F(H)$ on $M$ such that
$$
\sigma(\nabla H) = \mathsf{R}\bigl(F(H)\bigr),
$$
and this $F(H)$, which is a linear, first order differential expression in $H$ (whose coefficients involve the curvature of $g$), is the only possible candidate for $\nabla f$.  However, in order for this to solve our problem, it must satisfy
$$
\nabla \bigl(F(H)\bigr) - H = 0,\tag{2}
$$ 
and this is a system of $n^2$ first-order equations on $F(H)$, which is, of course, a system of $n^2$ second-order equations on $H$.  If $H$ does satisfy this system, then, because $\nabla\bigl(F(H)\bigr)=H$ is symmetric, it will follow that $F(H)$ is (locally) a gradient vector field of some function $f$, uniquely determined up to an additive constant.
For example, when $n=2$ and the Gauss curvature of $g$ is nowhere zero, the equations $(1'')$ are trivial since $\mathsf{R}(T) = T^*\otimes\Lambda^2(T^\ast)$.  Thus, the result of this analysis is that there exists a linear, second order differential operator $\mathsf{S}_g$ from symmetric quadratic forms to quadratic forms, namely $\mathsf{S}_g(H) = \nabla\bigl(F(H)\bigr)-H$, such that a symmetric quadratic form $H$ is (locally) a Hessian with respect to $g$ if and only if $\mathsf{S}_g(H)=0$.  Moreover, when the first deRham cohomology group of $M$ vanishes and $\mathsf{S}_g(H)=0$, there will exist a global function $f$ such that $\mathrm{Hess}_g(f)=H$, and this $f$ is unique up to an additive constant.
When $n>2$, the equations $(1'')$ are never trivial, so the condition that $H$ be a Hessian with respect to $g$ involves first order conditions and second order conditions.  What is interesting is that, when $n$ is sufficiently large and $\mathsf{R}$ is sufficiently 'generic', it appears (I haven't checked for sure) that $(1'')$ can actually imply $(2)$, so that the conditions become first order again.  (This is not true when $n=3$, though.)
Finally, if $\mathsf{R}$ is not injective, one may have to go to higher order derivatives of $H$ to determine whether it is a Hessian. This is an interesting case, but I don't have time to go into it right now. 
Remark:  Of course, the metric $g$ is not really essential in this story, since it really depends more on $\nabla$, than on $g$.  The equation $\nabla(df) = H$ makes sense without any metric and is a reasonable equation, as long as $\nabla$ is a torsion-free connection on $M$.  Thus, one could ask for characterizations of those symmetric quadratic differentials $H$ that can be written in the form $H = \nabla(df)$ where $\nabla$ is a given torsion-free connection on $M$.  The answer in this case would have essentially the same form as the answer above since, if one carries out the computations carefully, one sees that the metric $g$ plays essentially no rôle in the problem.
A: Given $f : N \to \mathbb R$, there is its derivative $Df : TN \to \mathbb R$, which has its dual gradient $\nabla f: N \to TN$, which you can covariantly differentiate $c \circ D\nabla f : TN \to TN$, and this map is dual to the Hessian $Hf : TN \oplus TN \to \mathbb R$.  From this perspective there's two things to determine.  Here $c : T^2 N \to TN$ is the connection, in the Ehresmann formalism. 
1) Your Hessian is adjoint to a bundle map $TN \to TN$, is this the covariant derivative of a vector field on $N$? 
2) If the answer to (1) is yes, among the solution vector fields from (1) is there a vector field which is the gradient of a real-valued function on $N$?  
(2) Has the traditional cohomological answer so I'll focus on (1). 
A bundle map $TN \to TN$ is the covariant derivative of a vector field $N \to TN$ means that you could write the vector field as a type of holonomy integral, by adding to your parallel transport an integral of the bundle map $TN \to TN$. So there will be a local triviality condition on this holonomy integral, as well as a global triviality condition. 
