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.