Let $(M^n,g)$ be a smooth Riemannian manifold. Consider the square of the distance function $$dist^2\colon M\times M\to \mathbb{R}$$ given by $(x,y)\mapsto dist^2(x,y)$. It is easy to see that this function is infinitely smooth near the diagonal.

Now fix a point $a\in M$. Consider the Taylor series of $dist^2$ at $(a,a)$ in some coordinate system (say normal). Can one compute explicitly its coefficients up to the third order?