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?



Characterizing Hessians among symmetric bilinear tensors
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.


