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.