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.