Why is $\Gamma^i_{\ \ jk}(x) \approx \frac{1}{2}R^i_{\ j\ell k}(x_0) (x^\ell-x_0^\ell)$ in normal coordinates of a curved space? This comes from an earlier post. Why is this approximation true?

We can always choose coordinates such that the Christoffel symbols vanish at an arbitrary point $x_0$ of our choosing, but the presence of spacetime curvature means that they will not vanish in a nonzero neighborhood of $x_0$.  The best we can do is choose Riemann normal coordinates, in which case
$$\Gamma^i_{\ \ jk}(x) \approx \frac{1}{2}R^i_{\ j\ell k}(x_0) (x^\ell-x_0^\ell)$$

-J.Murray. https://physics.stackexchange.com/questions/614243/how-do-you-describe-a-geometry-where-the-christoffel-symbols-vanish
I believe $\Gamma$ and ${R^a}_{bcd}$ should be of the same order, due to how ${R^a}_{bcd}$ is defined in terms of $\Gamma$ by parallel transport. This approximation seems related to Taylor's theorem.
 A: I think that the correct expressions for Riemann normal coordinates about an arbitrarily chosen origin $x^\mu=0$ are
$$
g_{\mu\nu}(x)= \delta_{\mu\nu}- \frac 13 R_{\mu\sigma \nu\tau}(0) x^\sigma x^\tau + O(|x|^3).
$$
and
$$
{\Gamma^{\lambda}}_{\mu\nu}(x)= -\frac 13 (R_{\lambda\nu\mu\tau}+R_{\lambda\mu\nu\tau})x^\tau+ O(|x|^2).
$$
Suppose that we have constructed  coordinates at O such that
$$
g_{\mu\nu}(x)= \delta_{\mu\nu}+\frac 12 A_{\mu \nu\sigma \tau} x^\sigma x^\tau + O(|x|^3),
$$
where $A_{\mu\nu\sigma\tau}$  is symmetric under $\mu\leftrightarrow \nu$ and under $\sigma\leftrightarrow \tau$.
Then
$$
{\Gamma^{\lambda}}_{\mu\nu}(x)=\frac 12(A_{\lambda\mu\nu\tau}+A_{\lambda\nu\mu\tau}-A_{\mu\nu\lambda\tau})x^\tau+O(|x^2|),
$$
and
$$
R_{\rho\sigma\mu\nu}(0)= \frac 12(A_{\rho\nu\sigma\mu}-A_{\nu\sigma\mu\rho}+A_{\mu\sigma\nu\rho}-A_{\rho\mu\sigma\nu}).
$$
We can verify that this curvature tensor  satisfies the pair exchange symmetry $R_{\mu\nu\rho\sigma}= R_{\rho\sigma\mu\nu}$, the two antisymmetries
$$
R_{\rho\sigma\mu\nu}=- R_{\sigma\rho\mu\nu}= -R_{\rho\sigma\nu\mu}
$$
and  the first Bianchi identity:
$$
R_{\rho\sigma\mu\nu}+R_{\rho\mu\nu\sigma}+R_{\rho\nu\sigma\mu}=0.
$$
Now in d=4, for example, the  array of numbers  $A_{\mu\nu\sigma\tau}$ has 10$\times$10=100 independent entries  while its  symmetries lead $R_{\rho\sigma\mu\nu}$ to have  only 20. We have,  however, at our disposal  80
degrees of freedom in the  coefficients ${b^\mu}_{\rho\sigma\tau}$
$$
x^\mu \to x^\mu + {b^\mu}_{\rho\sigma\tau} x^\rho x^\sigma x^\tau+\ldots
$$
of a  local change-of-coordinates expansion that keeps the metric euclidean up to quadratic corrections.
We should therefore be able to find a co-ordinate system in which $A_{\nu\sigma\mu\rho}$  is expressed in terms of $R_{\rho\sigma\mu\nu}$.  Indeed, counting degrees of freedom shows that  in any number of dimensions   we can use a  cubic change of co-ordinate to  reduce $A_{\mu\nu\sigma\tau}$ to the form
$$
A_{\mu\nu\sigma\tau}= \frac{c}{2}(R_{\mu\sigma\nu\tau}+R_{\nu\sigma\mu\tau}).
$$
Note that this expression  has the correct $\mu\leftrightarrow \nu$ and $\sigma\leftrightarrow \tau$ symmetries.
If we plug this form for $A_{\mu\nu\sigma\tau}$ into the formula for $R_{\rho\sigma \mu\nu}$ and the use the symmetries (including the first Bianchi identity)  we find
$$
R_{\rho\sigma \mu\nu} =\frac{c}{4}(-6 R_{\rho\sigma \mu\nu}),
$$
and so $c=-2/3$.
Thus we can  find coordinates in which
$$
g_{\mu\nu}(x)= \delta_{\mu\nu}- \frac 13 R_{\mu\sigma \nu\tau}(0) x^\sigma x^\tau + O(|x|^3).
$$
and
$$
{\Gamma^{\lambda}}_{\mu\nu}(x)= -\frac 13 (R_{\lambda\nu\mu\tau}+R_{\lambda\mu\nu\tau})x^\tau+ O(|x|^2).
$$
These  new coordinates  are precisely the geodesic normal coordinates. In other words,  straight lines through the origin in ${\mathbb R}^d$ are geodesics to $O(|x|^2)$.
We can do something similar for Vielbeins and the assocated "spin connection". We seek an orthonormal  vielbein co-frame ${\bf e}^{*a}$ such that
$$
e^{*a}_\mu e^{*a}_\nu=g_{\mu\nu}= \delta_{\mu\nu}- \frac 13 R_{\mu\sigma \nu\tau}(0) x^\sigma x^\tau + O(|x|^3).
$$
We can take
$$
e^{*a}_\mu= \delta_{a \mu}- \frac 16  R_{a \sigma \mu\tau}(0) x^\sigma x^\tau +\ldots
$$
and confirm that
$$
g^{\mu\nu} e^{*a}_\mu e^{*b}_\nu= \delta^{ab}+ O(|x|^3).
$$
Knowing the ${\bf e}^{*a}$, we  compute the spin connection from the definition
$$
\nabla_\mu {\bf e}^{*a}=-{\omega^{a}}_{b\mu} {\bf e^{*b}}
$$
as
$$
{\omega^a}_{b\mu}(x)= -e^\nu_b(\partial_\mu e^{*a}_\nu-{\Gamma^\sigma}_{\nu\mu}e^{*a}_\sigma)\\
= \frac 16 ({R^a}_{\mu b\tau}+{R^a}_{\tau b\mu})x^\tau-\frac 13({R^a}_{\mu b\tau}+{R^a}_{b\mu\tau})x^\tau+\ldots\\
= \left(-\frac 16 {R^a}_{\mu b\tau}-\frac 16 {R^a}_{\tau\mu b}-\frac 13 {R^a}_{b\mu\tau}\right)x^\tau+\ldots\\
=\left(+\frac 16 {R^a}_{b\tau\mu}-\frac 13 {R^a}_{b\mu\tau}\right)x^\tau+\ldots\quad\hbox{(Bianchi)} \\
= - \frac 12 {R^a}_{b\mu\tau}(0)x^\tau+O(|x|^2).
$$
This last equation looks like the one you cited, so it is for the orthonormal frame-field "spin" connection and not for the coordinate frame Christofflel symbols
A: This is not a general expression, it applies in a very special set of coordinates $x^{\hat{\mu}}$ known as Riemann normal coordinates. These are defined by the derivatives of the metric being zero at the origin of the coordinates $g_{\hat{\mu}\hat{\nu},\hat{\kappa}}|_{x^{\hat{\lambda}} = 0} = 0$. This leads to non-trivial coordinate conditions that generic coordinates simply do not fulfill!
As a result, you have vanishing Christoffel symbols at the origin, and
$$g_{\hat{\mu}\hat{\nu}} = g_{\hat{\mu}\hat{\nu}}|_0 + g_{\hat{\mu} \hat{\nu},\hat{\kappa}\hat{\lambda}}|_{0}x^{\hat{\kappa}} x^{\hat{\lambda}} +... $$
From that you can derive that
$$g_{\hat{\mu} \hat{\nu},\hat{\kappa}\hat{\lambda}}|_{0} = -\frac{1}{3} R_{\hat{\mu}\hat{\kappa}\hat{\nu}\hat{\lambda}}$$
Then it is easy to derive the identity in the OP.
A: I think I made a mistake in that answer - the standard expression for the Christoffel symbols in normal coordinates is
$$\Gamma^i_{\ \ jk} = \frac{1}{3}(R^i_{\ \ j\ell k} + R^i_{\ \ k \ell j})(x^\ell-x_0^\ell)$$
I think my expression is for the spin connection $\omega^i_{\ \ jk}$.  I'll fix that, but it doesn't affect the motivation for your question.

In Riemann normal coordinates, the metric goes like
$$g_{\mu\nu} \approx \eta_{\mu\nu} -\frac{1}{3} R_{\mu\alpha \nu\beta} (0) x^\alpha x^\beta +\mathcal O(x^3)$$  Since the $\Gamma$'s are constructed from the derivatives of the metric, you can construct them immediately and see that they are linear in $x$.  Note that the $R$ which appears in this expansion is the Riemann-Christoffel tensor evaluated at the coordinate origin.
