Let $M = (M, g)$ be a Riemannian manifold, and let $p \in M$.
Writing $d$ for the geodesic distance in $M$, there is a function
$$
d(-, p)^2 : M \to \mathbb{R}.
$$
This function is smooth near $p$. Hence for each point $x \in M$ sufficiently close to $p$, we have the Hessian
$$
\text{Hess}_x(d(-, p)^2)
$$
(defined using the Levi-Civita connection), which is a bilinear form on $T_x M$. In particular, we can take $x$ to be equal to $p$ itself, giving a bilinear form
$$
\text{Hess}_p(d(-, p)^2)
$$
on $T_p M$. But of course, we already have another bilinear form on $T_p M$, namely, the Riemannian metric $g_p$ itself. And the fact is that up to a constant factor, these two forms are equal:
$$
g_p = \frac{1}{2} \text{Hess}_p(d(-, p)^2).
$$
I'm looking for a reference for this fact. For the purposes of what I'm writing, it would ideally be a reference that states this fact in the same simple direct terms as above, without involving any other differential-geometric concepts (e.g. normal coordinates).
I understand that this is a basic fact of Riemannian geometry, so I've already looked for it in various introductions to the subject, including those by do Carmo, Jost, Lee, and Petersen. But I haven't found it stated in any of those sources (which isn't to say it's not there). I have found more sophisticated stuff about $\text{Hess}_x(d(-, p)^2)$ for points $x$ different from $p$, but not the simple fact I'm looking for.
Requests for references often result in people giving their favourite proofs rather than a reference. While that doesn't do any harm (and can be quite interesting), I emphasize that it's a reference I'm looking for, not a proof.
Best Answer
While it does not answer your question, the following direct argument may clarify certain things:
Since the Hessian is a symmetric bilinear form, it suffices to show $\frac{1}{2}Hess_p(d^2(\cdot,p))(v,v)=|v|^2$.
If $p$ is a critical point of a smooth function $f$ on $\mathbb{R}^n$, then $Hess_p(f)(v,v)=\frac{d^2}{dt^2}\vert_{t=0}f(\gamma(t)) $, where $\gamma$ is any smooth path with $\gamma(0)=p$ and $\gamma'(0)=v$. This formula continues to hold, if $p$ is a critical point of a function $f$ on a manifold (in this case the definition of the Hessian does not rely on the choice of a Riemannian metric).
If $\gamma$ is the geodesic through $p$ with $\gamma'(0)=v$, then, since $\gamma$ is locally distance minimizing, $d(\gamma(t),p)=|tv|$ for $t$ near $0$. Combined with the above this gives the result.
(If we use a Riemannian metric $g$ and its associated Levi-Civita connection $\nabla$ to define $Hess_p^g(f)$ at a noncritical point $p$ of $f$, then the formula $\frac{d^2}{dt^2}\vert_{t=0}f(\gamma(t)) =Hess^g_p(f)(v,v)$ still holds, if $\nabla_t\gamma'(0)=0$. This is however not used above).