Let $M$ be a Riemannian manifold of dimension $n$, $d$ be the natural distance function on it, $\Delta_y$ be the Laplace-Beltrami operator acting in the variable $y \in M$, and $\operatorname{Ric}$ be its Ricci tensor. I have found in some notes of mine from a few years ago the formula $$\Delta_y [d(x, y)^2] = 2n - \frac 2 3 \operatorname{Ric}_y (\exp_y ^{-1} x, \exp_y ^{-1} x) + O(d(x,y)^3)$$ whwre the last term is written in the big-O notation.
I have completely forgotten how I have got to this formula, and even whether it is correct at all. In particular, even if I worked in normal coordinates around $x$, I would still not know how to come up with $\operatorname{Ric}_y$ (notice that the tensor is at $y$, not at $x$). Could anyone please help me with either a reference, or with the sketch of a proof (if the main lines of the proof were clearly given I could fill the details in myself)? Thank you.