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))=vw^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?

$\begingroup$ See chapter 5 in DoCarmo's Riemannian Geometry book. Specifically, this one is Proposition 2.7 there $\endgroup$– David LehaviJun 6, 2021 at 10:02
2 Answers
The "standard" proof using Jacobi fields can be found in section 1.3 here https://www2.math.upenn.edu/~wziller/math660/TopogonovTheoremMyer.pdf

$\begingroup$ hi Otis, I am studying this expansion provided by Myer's paper you suggested. There is an equality in display (24) on page 8 of the paper that I do not understand. It is $\nabla_{D_{\epsilon}} R(E,T)E _{0,t} = R(E, \nabla_{D_{t}}E)E_{0,t}$, where $\nabla_{D_{t}}E = \nabla_{D_{\epsilon}}T$. Do you know how this is derived? Thanks. $\endgroup$– CheeOct 13, 2021 at 2:23

$\begingroup$ @Chee, can you be more specific about your issue with this? I think you should use that T vanishes at eps =0 so the given term is the only one that can appear. $\endgroup$ Oct 13, 2021 at 3:05

$\begingroup$ hi Otis, thank you for your reply. I checked every thing in the proof but the identity was the only thing I did not understand. I do not see how $T$ vanishes at $\epsilon=0$ leads to the identity. The lefthand side is the covariant derivative of $R(E,T)E$ whereas the righthand side is $R(E,\cdot)E$ with $\cdot$ as the covariant derivative of $T$. In a sense, for this specific case, $\nabla_{\epsilon}$ "commutes" with curvature endomorphism. Strange ... $\endgroup$– CheeOct 13, 2021 at 3:24

$\begingroup$ Hi Otis, in the display (24) of Myer's paper, the first equality was obtained by the definition of curvature via 3thorder mixed directional derivatives, and to prove this identity, we should not expand the curvature via its definition and then covariant differentiate each term in the expansion, since this will lead us back to the first equality there in (24). i have tried $(\nabla_{D_{\epsilon}} R) (E,T)E$, i.e., the covariant derivative of the curvature tensor, which equivalently should be $0$ at $(0,t)$ if the identity holds. I was not able to prove that it is $0$ at $(0,t)$ $\endgroup$– CheeOct 13, 2021 at 3:33

1$\begingroup$ Morning Otis, thank you for guiding me through it. I got it. $\nabla{_{\cdot}}$ preserves the type of the tensor it is differentiating, and since $T_{0,t}=0$, the covariant derivative i wrote and you pointed to is $0$ at $(0,t)$ $\endgroup$– CheeOct 13, 2021 at 15:33
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^iy^i\bigr)\bigl(x^jy^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).$$