Riemann's formula for the metric in a normal neighborhood I would love to understand the famous formula $g_{ij}(x) = \delta_{ij} + \frac{1}{3}R_{kijl}x^kx^l +O(\|x\|^3)$, which is valid in Riemannian normal coordinates and possibly more general situations.
I'm aware of 2 proofs: One using Jacobi fields [cf. e.g. S.Sternberg's "Curvature in Mathematics and Physics" from which the question title and formula is stolen :-) or cf. S.Lang's "Differential and Riemannian Manifolds"]. The other proof involves computing that $\partial_k\partial_lg_{ij}(x)$ shares some symmetries of curvature [cf. M.Spivak's "A Comprehensive Introduction to Differential Geometry, Vol. 2" where it is a several page "hairy computation" or cf. H.Weyl's 1923 edition of Riemann's Habilitationsvortrag (reprinted in a recent German book by Jürgen Jost) which I find uncomprehensible.]
Are you aware of any other proof? Are normal coordinates necessary?
While the Jacobi fields proof is short and elegant enough, it irks me that it requires "higher technology" not involved in the endproduct. Somehow the formula should be provable by pure calculus. Indeed, it is stated as an exercise in P.Petersen's "Riemannian Geometry": From the context I guess he thinks it should follow from the expression of $\partial_lg_{ij}$ as a sum of 2 Christoffel symbols and the simplified expression for curvature at $x=0$ where the Christoffel symbols vanish. Alas my attempts at this go in circles...
I find the situation quite amazing: Not many textbooks treat this fundamental and historic formula. (Estimating from the sample on my shelf it is $3/17.$ E.g. it seems it's not even in Levi-Civita's classic.)
Update/Scholium:
In classical language: The knackpoint seems to be a "differential Bianchi formula" for the Christoffel symbols at $0$. This follows from the geodesic equation. I see no other way yet.
A more modern approach minimizing (but not eliminating) the role of geodesics is in A.Gray's Tubes book. (Noted in comments. I'm waiting for www.amazon.de to deliver this treasure.)
$\bullet$ While geodesics are very geometric and normal coordinates are very practical, methinks the formula is a tad ungeometric. What I'm hoping/asking for is a coordinate-independent formula for the second derivative of $g$ in terms of a suitable "reference connection".
 A: I came to this post many years later, since I too was concerned about the absence of Riemann's formula in most texts, lengthy treatment in others, or reliance on more advanced techniques like Jacobi fields. I include here a direct concise proof which I think would be well suited for beginning students. We want to show that
$$
g_{ij}(x)=\delta_{ij}+\frac{1}{3} R_{k ij\ell}(o)x_kx_\ell+\mathcal{O}(|x|^3),
$$
where $x=(x_1,\dots, x_n)$ are normal coordinates centered at a point $o$ in a Riemannian manifold $M$.  By Taylor's theorem, we need to check that
\begin{eqnarray}
\tag{1}  g_{ij}(o)&=&\delta_{ij}\label{eq:1},\\ 
\tag{2} g_{ij,k}(o)&=&0\label{eq:2},\\
\tag{3} g_{ij,k\ell}(o)&=&\frac{1}{3}\big(R_{kij\ell }(o)+R_{\ell ij k}(o)\big)\label{eq:3},
\end{eqnarray}
where we use the notation $f_{,i}:=\partial_i f$, and $f_{,ij}:=\partial^2_{ij} f$.
\eqref{eq:1} follows immediately from the construction of normal coordinates, and \eqref{eq:2} is not difficult to establish either. \eqref{eq:3}, which is the heart of the matter, requires a bit more work.
In most sources, like Spivak, vol II or the  relatively recent  book by Jost, which gives a very comprehensive treatment of Riemann's lecture, \eqref{eq:3} is established via symmetry properties of $g_{ij,k\ell}$ which involve long computations; although 2011 Lecture Notes of John Douglas Moore gives a very nice and efficient proof of them. Instead I will demonstrate \eqref{eq:3} via a Bianchi type cyclic identity for derivatives of Christoffel symbols, which follows quickly from the geodesic equation (as the OP had mentioned in the "Update/Scholium" above). The only place I have seen this approach is in 2013 Lecture Notes of Christian Bär.
Proofs of (1) and (2)
All indices here range from $1$ to $n$, and any term which involves repeated indices stands for a sum over that index.
Let $\exp_o\colon T_o M\to M$ be the exponential map, and $U\subset M$ be a ball centered at $o$ such that $\exp_o\colon \exp_o^{-1}(U)\to U$ is a diffeomorphism. Let $e_i$ be an orthonormal basis for $T_o M$, i.e.,
$
g(e_i, e_j)=\delta_{ij}.
$
Then the normal coordinates $x_i\colon U\to \mathbf{R}$ (with respect to $e_i$) are given by
$$
\exp_o(x_i(p) e_i)=p.
$$
The mapping $x\colon U\to\mathbf{R}^n$, given by $x:=(x_1,\dots, x_n)$ identifies $U$ with a ball centered at the origin in $\mathbf{R}^n$, which  we again denote by $U$. Let $E_i(x):=\partial_i|_x$ be the coordinate vector fields on $U$. Then $g_{ij}\colon U\to\mathbf{R}$ are given by
$$
g_{ij}(x):=g(E_i(x),E_j(x)).
$$
Since $E_i(o)=e_i$ we immediately obtain \eqref{eq:1}. To see \eqref{eq:2} note that, by the definition of normal coordinates $x_i$, the geodesics in $U$ passing through $o$ and another point $x$ of $U$ are given by $\gamma(t):=tx$. Since $\gamma$ is a geodesic,
$$
\gamma_k''(t)+\Gamma_{ij}^k(\gamma(t))\gamma_i'(t)\gamma_j'(t)=0,
$$
which yields
\begin{equation}\label{eq:4}
\tag{4}\Gamma_{ij}^k(t x)x_i x_j=0.
\end{equation}
Setting  $t=0$, observing that $x_i$, $x_j$ may assume any values,  and recalling that $\Gamma_{ij}^k=\Gamma_{ji}^k$, we obtain
\begin{equation*}\label{eq:Gamma0}
\Gamma_{ij}^k(o)=0.
\end{equation*}
Since
$
\Gamma_{ij}^k=\frac{1}{2} g^{k\ell}(g_{\ell i,j}+g_{\ell j,i}-g_{ij,\ell}),
$
and $g^{ij}(o)=\delta_{ij}$, we have
\begin{equation*}\label{eq:8}
0=\Gamma_{ij}^k(o)=\frac{1}{2} \big(g_{k i,j}(o)+g_{k j,i}(o)-g_{ij,k}(o)\big).
\end{equation*}
Adding the above equation to itself, after a cyclic permutation of indices yields \eqref{eq:2}.
Proof of (3)
Since
$
\nabla_{E_i}E_j=\Gamma_{ij}^kE_k,
$
$$
g_{ij,k}=g(\nabla_{E_k}E_i, E_j)+g(E_i, \nabla_{E_k}E_j)=\Gamma_{k i}^\ell g_{\ell j}+\Gamma_{k j}^\ell g_{i\ell}.
$$
Differentiating again, and using \eqref{eq:2}, yields
\begin{equation}\label{eq:5}
\tag{5}g_{ij,k\ell}(o)=\Gamma_{k i,\ell}^j(o)+\Gamma_{k j,\ell}^i(o).
\end{equation}
Next we differentiate \eqref{eq:4} at $t=0$ to obtain
$$
\Gamma_{ij,\ell}^k(o)x_i x_jx_\ell=0,
$$
a homogeneous polynomial of degree $3$ which vanishes identically. The coefficient of each term $x_i x_jx_\ell$ is the sum of all $6$ permutations of lower indices of $\Gamma_{ij,\ell}^k(o)$. Since $\Gamma_{ij,\ell}^k=\Gamma_{ji,\ell}^k$, we obtain
\begin{equation*}\label{eq:cyclic}
\Gamma_{ij,\ell}^k(o)+\Gamma_{j\ell,i}^k(o)+\Gamma_{\ell i,j}^k(o)=0.
\end{equation*}
Now note that, since
$
R_{ijk}^\ell=\Gamma_{ik,j}^\ell-\Gamma_{jk,i}^\ell+\Gamma_{ik}^p\Gamma_{pj}^\ell-\Gamma_{jk}^p\Gamma_{ip}^\ell,
$
\eqref{eq:2} yields that
$$
R_{ijk\ell}(o)=\Gamma_{ik,j}^\ell(o)-\Gamma_{jk,i}^\ell(o).
$$
Here we have also used the fact that $R_{ijk}^\ell(o)=R_{ijk\ell}(o)$ due to \eqref{eq:1}. The last two displayed equations yield
\begin{eqnarray*}
R_{ik\ell j}(o)+R_{i\ell kj}(o)
= \Gamma_{i\ell,k}^j(o)+\Gamma_{ik,\ell}^j(o)-2\Gamma_{\ell k,i}^j(o)
=-3\Gamma_{\ell k,i}^j(o).
\end{eqnarray*}
The last equality together with \eqref{eq:5} and symmetries of $R$ now yields
\begin{eqnarray*}
g_{ij,k\ell}(o)
&=&-\frac{1}{3}\big(R_{\ell ik j}(o)+\require{cancel}\cancel{R_{\ell kij}(o)}+R_{\ell jk i}(o)+\require{cancel}\cancel{R_{\ell kji}(o)\big)}\\
&=&\frac{1}{3}\big(R_{\ell i jk}(o)+R_{kij\ell}(o)\big),
\end{eqnarray*}
as desired.
A: Still another approach to the Riemann normal coordinates expansion formula can be found in http://arxiv.org/abs/gr-qc/9712092 (A Closed Formula for the Riemann Normal Coordinate Expansion, by U. Mueller, C. Schubert and A. van de Ven). This approach indicates that the Riemann normal coordinates are the gravity analogue of the Fock-Schwinger gauge in gauge theory. Fock-Schwinger gauge (centered at the origin) is defined by the condition $$x^\mu A_\mu(x)=0. \tag{1}$$
In the local neighbourhood of the origin, the condition (1) can be solved in terms of the following integral representation
$$A_\mu(x)=x^\nu\int\limits_0^1F_{\nu\mu}(s x)s ds, \tag{2}$$
which connects the gauge potential $A_\mu$ and the field strength tensor $F_{\mu\nu}$. As a result, the Taylor expansion coefficients of $A_\mu$ at the origin is expressed through the covariant derivatives of $F_{\mu\nu}$.
In analogy, Riemann normal coordinates centered at the origin can be defined by the conditions
$$g_{\mu\nu}(0)=\delta_{\mu\nu},\;\;\;x^\mu g_{\mu\nu}(x)=x^\mu g_{\mu\nu}(0). \tag{3}$$
(the second condition is equivalent to the following condition on the Chrisoffel symbol $x^\mu x^\nu \Gamma_{\mu\nu}^\lambda(x)=0$, which determines the coordinate system locally up to a rigid rotation). 
Clearly, (3) is the analog of (1). While the analog of (2), proved in the Mueller, Schubert and van de Ven paper, is
\begin{eqnarray} &&
g(x)=\sum\limits_{k=0}^\infty\int\limits_0^1ds_1\,(1-s_1)\int\limits_0^1ds_2\,(1-s_2)\cdots \int\limits_0^1ds_k\,(1-s_k) \\ &&
\times \sum\limits_{l=0}^ks_1s_2^3\cdots s_l^{2l-1}s_{l+1}^{2k-2l-1}
s_{l+2}^{2k-2l-3}\cdots s_k \\ && \times {\cal R}(s_1 s_2 \cdots s_l x,x){\cal R}(s_2 s_3 \cdots s_l x,x)\cdots {\cal R}(s_l x,x)
\\ && \times {\cal R}(s_{l+1} x,x){\cal R}(s_{l+1}s_{l+2} x,x)\cdots
{\cal R}(s_{l+1}s_{l+2}\cdots s_k x,x), \tag{4}
\end{eqnarray}
where $${\cal R}^\mu_{\;\nu}(x,y)=R^\mu_{\;\alpha\beta \nu}(x)y^\alpha y^\beta.$$
The generalization of (4) to the Fermi normal coordinates in tubular geometry is considered in http://arxiv.org/abs/1203.1151 (All order covariant tubular expansion, by P. Mukhopadhyay).
A: Perhaps the simplest way to understand this formula is to think about how you would go about deriving it:  Try to find the 'best' coordinates you can centered on a given point and see what doesn't change in such coordinates.
Suppose $g$ is a Riemannian metric on $M$ and $p\in M$ is fixed.  Start by choosing a $p$-centered local coordinate system $x = (x^i)$ on $U\subset M$ and write
$$
g = g_{ij}(x)\,\mathrm{d}x^i\mathrm{d}x^j.
$$
Since $\bigl(g_{ij}(0)\bigr)$ is a positive definite matrix, you can make a linear change of coordinates in $x$ so that $g_{ij}(0) = \delta_{ij}$.  Call such a $p$-centered coordinate system $0$-adapted to $g$ at $p$.
Now, ask what would be the effect of expressing $g$ in the coordinates $y=(y^i)$ that are related to the coordinate $x$ by $x^i = y^i + \tfrac12a^i_{jk} y^jy^k$ for some $a^i_{jk} = a^i_{kj}$.  It is easy to see by Taylor series expansion that you can uniquely choose the $a^i_{jk}$ so that, when we write 
$$
g = \bar g_{ij}(y)\,\mathrm{d}y^i\mathrm{d}y^j,
$$
we have, for all $i$, $j$, and $k$,
$$
\frac{\partial\bar g_{ij}}{\partial y^k}(0) = 0.
$$
(It's clear that this is the same number of equations as unknowns for the $a^i_{jk}$, one just has to check that the inhomogeneous system of equations has only the zero solution when the inhomogeneous part is set to zero.)  Call such a system of $p$-centered coordinates $1$-adapted to $g$ at $p$.  Thus, for a system of coordinates $y$ that is $1$-adapted to $g$ at $p$, one has
$$
g = \left(\delta_{ij} 
+ \tfrac12 \frac{\partial^2g_{ij}}{\partial y^k\partial y^l}(0)\, y^ky^l
+ R^3_{ij}(y)\right)
\ \mathrm{d}y^i\mathrm{d}y^j,
$$
where $R^3_{ij}(y)$ vanishes to order $3$ at $y=0$.
Finally, consider what such a metric would look like in the coordinates $z = (z^i)$
that are defined by $y^i = z^i + \tfrac16 b^i_{jkl} z^jz^kz^l$ for some constants $b^i_{jkl} = b^i_{kjl} = b^i_{jlk}$.  Now, there are $n^2(n{+}1)(n{+}2)/6$ unknowns $b^i_{jkl}$, but there are $n^2(n{+}1)^2/4$ quantities in the second-order Taylor expansion of 
$g = {\bar g}_{ij}(z)\mathrm{d}z^i\mathrm{d}z^j$, i.e.,
$$
g = \left(\delta_{ij} 
+ \tfrac12 \frac{\partial^2{\bar g}_{ij}}{\partial z^k\partial z^l}(0)\, z^kz^l
+ {\bar R}^3_{ij}(z)\right)
\ \mathrm{d}z^i\mathrm{d}z^j.
$$
Thus, the equations $\frac{\partial^2{\bar g}_{ij}}{\partial z^k\partial z^l}(0)=0$,
as linear equations for the $b^i_{jkl}$, are overdetermined by 
$$
n^2(n{+}1)^2/4 - n^2(n{+}1)(n{+}2)/6 = n^2(n^2{-}1)/12
$$
equations.  
It is not hard to see that the corresponding homogeneous equations in the $b^i_{jkl}$ have only the solution $b^i_{jkl}=0$.  In fact, the $b^i_{jkl}$ are uniquely determined by requiring that, when we compute the Taylor expansion about $z=0$ we get
$$
g = \left(\delta_{ij} 
+ \tfrac12 h_{ij,kl}\, z^kz^l
+ R^3_{ij}(z)\right)
\ \mathrm{d}z^i\mathrm{d}z^j,
$$
with $h_{ij,kl}+h_{ik,lj}+h_{il,jk}=0$ (which is $n^2(n{+}1)(n{+}2)/6$ independent equations on the $b^i_{jkl}$).  Say that a system of coordinates $z = (z^i)$ centered at $p$ for which $g$ has its Taylor expansion at $p$ of the above form is $2$-adapted to $g$ at $p$.  Two such coordinate systems at $p$ are related in the form $z^i = a^i_j\,\bar z^j + O(|{\bar z}|^4)$, where $a = (a^i_j)$ is an orthogonal matrix.
Thus, the $2$-adapted condition forces the $h_{ij,kl}$ to lie in a vector space of dimension $n^2(n^2{-}1)/12$, as explained above. 
It's now a matter of linear algebra to show, as Riemann did, that these conditions imply that the $h_{ij,kl}$ can be written uniquely in the form
$$
h_{ij,kl} = \tfrac13(R_{kijl}+R_{lijk})
$$
where $R_{ijkl}=-R_{jikl}=-R_{ijlk}$ and $R_{ijkl}+R_{iklj}+R_{iljk}=0$.
A: My question has been answered in comments by Liviu Nicolaescu:
The (almost) ultimate proof (for my taste) is via A. Gray's formula(e) for (symmetric) higher covariant derivative(s) of normal coordinate vector fields. Any exposition of normal coordinates lacking this formula is severely lacking. (I'd prefer symmetrized c.d. of differentials of normal coordinates...). It is valid for any symmetric connection and gives Riemann's formula by application to some parallel bilinear form (e.g. the Riemann metric). (Riemann's really original Habilitationsvortrag formula needs a little additional algebra, cf. Dedekind/Weber [below] or Spivak.)
The second best approach is using Jacobi fields (cf. e.g. Le Spectre (LNM 194) for higher order terms without detail. (For hardcore syntacticists: Strook, Intro to An. of Paths on a Riem. Mf.)).
Third best: Repere mobile, but this seems Riemannian. (Cf. e.g. Heat Kernels and Dirac Operators, or Atiyah/Bott/Patodi appendix, or Cartan's Geometrie des Espaces de Riemann.)
The historically first proof is due to Dedekind/Weber (Anmerkung in "Bernhard Riemann's Gesammelte Mathematische Werke und Wissenschaftlicher Nachlaß", 1876/1892) via proving Riemann's Commentatio formula (which is also in Spivak) for curvature. (I could streamline this quite amazingly by doing the Levi-Civita in cotangent space.)
H.Weyl (1919/1923) has the historically second proof in his German edition of Riemann's Habilitationsvortrag, based on the extrinsic Gauss lemma (i.e. e.g. Besse, Einstein Mfs, Thrm. 1.45), which is what Spivak has worked out, and possibly has inspired the repere mobile proof (explicit in Heat Kernels [above]) - I don't (and am unwilling to) comprehend Cartan's proof.)
The nicest classical Riemannian proof is in L.P.Eisenhart, Riemannian Geometry, 1925/1949
I still don't have an opinion on R.Bryant's answer. It is possibly exactly what Riemann did for his inaugural lecture 1854. But then, Riemann had at least in later work (Commentatio, 1861) introduced Christoffel symbols (written $p_{\iota,\iota^\prime,\iota^{\prime\prime}}$) and other formulas for curvature.
A: There is the following way of doing it. Let $\Gamma_{ijk}\equiv g_{is}\Gamma^{s}_{jk}$. Then (see here)
$$
g_{jk}(x)=\delta_{jk}+\int_0^1 \alpha ~x^i\Gamma_{[ji]k}(\alpha x)~d\alpha
$$
If we expand $\Gamma_{jik}(x)\equiv\sum_{n=1}^{\infty}\Gamma_{jika_1...a_n}x^{a_1}...x^{a_n}$ (where $\Gamma_{jika_1...a_n}$ are coefficients) then the formula above says
$$
g_{jk}(x)=\delta_{jk}+\frac{1}{3}x^ix^a\Gamma_{[ji]ka}+O(|x|^3).
$$
Using the same expansion in the formula defining the Riemann curvature in terms of the connection, shows that
$$
\Gamma_{ab[dc]}=R_{abcd}(0)
$$
Next we use the fact that $x^ix^j\Gamma^a_{ij}(x)=0$ in normal coordinates, which implies $\Gamma_{j(ika)}=0$. This together with the equation above implies $\Gamma_{ijkl}=-\frac{1}{3}R_{i(kj)l}(0)$ which in turn means
$$
g_{jk}=\delta_{jk}-\frac{1}{3}x^ix^aR_{jika}+O(|x|^3).
$$
I'm using $(...)$ and $[...]$ here without the symmetrisation factors.
