Reference request: Recovering a Riemannian metric from the distance function 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.
 A: This is described in painstaking detail in the paper of Xavier Pennec (2017). (Hessian of he Riemannian Squared Distance).
A: For the requested reference: I believe it should follow from inequalities (5.6.6) in Jost (2011, p. 235) (plus user_1789’s polarization argument) because $r(x)\to0$ as $x\to p$.
A: I believe that the reason why you cannot find the result that you are asking about printed anywhere is that it is, after all, a mere exercise in Riemannian computation. First, it is easy to show that if $f$ is smooth around $p$, then $(Hf)_{ij} = \partial^2_{ij}f - \Gamma_{ij} ^k \partial_k f$ in any system of coordinates around $p$. Now, since your $f = d_p^2$ has radial symmetry, it is natural to continue the work in spherical normal coordinates, i.e. you go in $T_pM$ through $\exp_p ^{-1}$ and there you introduce spherical coordinates $r, \sigma_1, \dots, \sigma_n$, with $n = \dim M$. Since $\Gamma_{ij}^k (p) = 0$ as a consequence of your coordinates being normal, you will have $(Hf)_{ij} (p) = (\partial^2_{ij}f) (p) = (\partial ^2 _{rr} r^2) (p) = 2$ (all the other second-order partial derivatives vanish at $p$ because $f=r^2$ does not contain the variables $\sigma_1, \dots, \sigma_n$).
On the other hand, it is known that in normal spherical coordinates the expression of the metric tensor is $g_{ij} = \delta_{ij} + o(r)$, so that $g_{ij} (p) = \delta_{ij} (p)$ (the Kronecker symbol), whence it follows that $(Hf)(p) = 2g(p)$ (the metric evaluated at $p$). See p.114 of I. Chavel, "Riemannian Geometry - A Modern Introduction", 2006, or the more general  theorem 2.53 of Cartan on p.83 of S. Rosenberg, "The Laplacian on a Riemannian Manifold", 1997, or Petersen's book cited here.
A: 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).
A: If I am not mistaken, in order to define the Hessian you need to fix a connection. I suspect the Riemannian metric you get will depend on this connection, as Finsler metrics also have distance functions.
