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][1], 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).$$


  [1]: https://mathoverflow.net/questions/37651/riemannian-surfaces-with-an-explicit-distance-function/360046#360046