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$ that satisfies $\sigma(x,x) = 0$, $\sigma(x,y) = \sigma(y,x)$, and the first order PDE (regarding $y$ as a parameter) $$ g^{ij}(x)\frac{\partial\sigma}{\partial x^i}\frac{\partial\sigma}{\partial x^j} - 4\sigma = 0. $$ The function $\sigma$ is completely determined by these conditions.

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, where $$ g_{ij}(x) = \delta_{ij} -\tfrac13\,R_{ikjl}\,x^kx^l + O(3), $$ 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).$$

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).$$

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. $$

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$ that satisfies $\sigma(x,x) = 0$, $\sigma(x,y) = \sigma(y,x)$, and the first order PDE (regarding $y$ as a parameter) $$ g^{ij}(x)\frac{\partial\sigma}{\partial x^i}\frac{\partial\sigma}{\partial x^j} - 4\sigma = 0. $$ The function $\sigma$ is completely determined by these conditions.

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, where $$ g_{ij}(x) = \delta_{ij} -\tfrac13\,R_{ikjl}\,x^kx^l + O(3), $$ 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.

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$ that satisfies $\sigma(x,x) = 0$, $\sigma(x,y) = \sigma(y,x)$, and the first order PDE (regarding $y$ as a parameter) $$ g^{ij}(x)\frac{\partial\sigma}{\partial x^i}\frac{\partial\sigma}{\partial x^j} - 4\sigma = 0. $$ The function $\sigma$ is completely determined by these conditions.

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, where $$ g_{ij}(x) = \delta_{ij} -\tfrac13\,R_{ikjl}\,x^kx^l + O(3), $$ 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).$$