Taylor expansion of the square of the distance function on a Riemannian manifold I have recently read the problem named "Square of the distance function on a Riemannian manifold"(enter link description here) and I am interested in the formula
$ d^2(exp_{x_0}(tv),exp_{x_0}(tw))=|v-w|^2t^2-\frac{1}{3}R(v,w,w,v)t^4+O(t^5) $.
However I can not solve it by myself. Can one give me a complete solution?
 A: The "standard" proof using Jacobi fields can be found in section 1.3 here https://www2.math.upenn.edu/~wziller/math660/TopogonovTheorem-Myer.pdf
A: This is just a remark about an alternative, more `low tech', derivation of this famous formula by using Taylor series.
It relies on this property of a Riemannian metric:  If $(M,g)$ is a Riemannian manifold and $y\in M$ is a given point, let $\delta_y(x)$ be the length of the shortest path from $y$ to $x$.  Then $\delta_y$ is not differentiable at $y$, but is smoothly differentiable on a punctured neighborhood of $y$, and satisfies $|\nabla \delta_y|=1$ there, i.e., its gradient has length $1$ and, in fact, $(\nabla \delta_y) (x)$ for $x$ in this punctured neighborhood is equal to the velocity at $x$ of the unit speed geodesic that starts at $y$ and passes through $x$.  Consequently, if $\sigma_y = (\delta_y)^2$ then one finds that this function is differentiable on a neighborhood of $y$ and it satisfies the smooth differential equation
$$
|\nabla \sigma_y|^2 = |\nabla d_y^2|^2 = 4d_y^2\,|\nabla d_y|^2 =  4\sigma_y.\tag1
$$
Thus, when $g$ is expressed in local coordinates $(x^1,\ldots, x^n)$ centered on $p\in U\subset M$, one can write $\sigma = d(x,y)^2$ on $U\times U$ (at least near the diagonal) as a smooth function of $(x,y)\in U\times U$ that satisfies $\sigma(x,x) = 0$, $\sigma(x,y) = \sigma(y,x)$, and the first order PDE
$$
g^{ij}(x)\frac{\partial\sigma}{\partial x^i}\frac{\partial\sigma}{\partial x^j} - 4\sigma = g^{ij}(y)\frac{\partial\sigma}{\partial y^i}\frac{\partial\sigma}{\partial y^j} - 4\sigma = 0.
$$
The function $\sigma$ is determined by these conditions plus the 'initial condition' $\sigma(x,0) = g_{ij}(0)\,x^ix^j+O(3)$.
Then, in local coordinates, expanding the above equation out in Taylor series and using the 'initial conditions' determines the Taylor series for $\sigma$.  In particular, in geodesic normal coordinates centered on $p$, where
$$
g_{ij}(x) = \delta_{ij} -\tfrac13\,R_{ikjl}\,x^kx^l + O(3),
$$
and $R_{ijkl}=-R_{jikl}=-R_{ijlk}=-R_{iklj}-R_{iljk}$,
examining the first three terms of the above Taylor series expansion of the PDE yields
$$
\sigma = \delta_{ij}\,\bigl(x^i-y^i\bigr)\bigl(x^j-y^j\bigr) -\tfrac1{12}\,R_{ikjl}\,\bigl(x^iy^k{-}x^ky^i\bigr)\bigl(x^jy^l{-}x^ly^j\bigr) + O(5),
$$
which is equivalent to the desired formula.
Remark:  To see another application of the formula (1), one might consult this answer of mine, where it is used to compute the explicit distance function for the complete metric of negative curvature on $\mathbb{R}^2$ given by $$g = (x^2+y^2+2)\bigl(\mathrm{d}x^2+\mathrm{d}y^2\bigr).$$
