This isn't a full answer but hopefully the beginning of one. I'm assuming it's a matrix group but it might work for all Lie groups anyway.
Let $G$ be a matrix group with Lie algebra $\mathfrak{g}$. Suppose $\gamma(t)$ is a geodesic in $G$ and $x(t)$ is the geodesic in the exponential coordinates, ie. $\gamma(t) = \exp(x(t))$. Then $x(t)$ satisfies an ODE $x''(t) + Q_{x(t)}(x'(t),x'(t)) = 0$, where, for given $x$, the $Q_x$ is some quadratic map $\mathfrak{g} \to \mathfrak{g}$. If we can find $Q$ then we just need to polarise it to a bilinear map $\mathfrak{g} \times \mathfrak{g} \to \mathfrak{g}$ and this bilinear map will contain all the Christoffel symbols.
There is a power series formula for the differential of the exponential map. See G. M. Tuynman, The Derivation of the Exponential Map of Matrices http://www.jstor.org/pss/2974511 - it implies:
$\frac{d}{d t}(\exp(x(t)) = \exp(x(t)) \left( \frac{1 - \exp(-\mathrm{ad}\ x(t))}{\mathrm{ad}\ x(t)}\right)x'(t)$
where $\mathrm{ad}\ x$ is the map which takes $y$ to $[x,y]$ and the fraction $\left( \frac{1 - \exp(-\mathrm{ad}\ x(t))}{\mathrm{ad}\ x(t)}\right)$ means $1 - \frac{1}{2!} \mathrm{ad}\ x(t) + \frac{1}{3!} \mathrm{ad}\ x(t)^2 - \frac{1}{4!} \mathrm{ad}\ x(t)^3 + \cdots,$ ie. it is the power series for the fraction as if $x(t)$ was a real number.
Because the covariant derivative is bi-invariant, the geodesic $\gamma$ must be of the form $\gamma(t) = \gamma(0) \exp(A t)$ for some $A \in \mathfrak{g}$. Differentiating $\exp(x(t)) = \gamma(0) \exp(A t):$
$\exp(x(t)) \left( \frac{1 - \exp(-\mathrm{ad}\ x(t))}{\mathrm{ad}\ x(t)}\right)x'(t) = \gamma(0) \exp(A t) A = \exp(x(t)) A.$
So $\left( \frac{1 - \exp(-\mathrm{ad}\ x(t))}{\mathrm{ad}\ x(t)}\right)x'(t) = A.$ Differentiating this,
$\frac{d}{dt}\left( \frac{1 - \exp(-\mathrm{ad}\ x(t))}{\mathrm{ad}\ x(t)}\right)x'(t) + \left( \frac{1 - \exp(-\mathrm{ad}\ x(t))}{\mathrm{ad}\ x(t)}\right)x''(t) = 0.$ So
$x''(t) + \left( \frac{1 - \exp(-\mathrm{ad}\ x(t))}{\mathrm{ad}\ x(t)}\right)^{-1} \frac{d}{dt}\left( \frac{1 - \exp(-\mathrm{ad}\ x(t))}{\mathrm{ad}\ x(t)}\right)x'(t)=0$, ie.
$\begin{align}x''(t) &+ \left( 1 - \frac{1}{2!} \mathrm{ad}\ x(t) + \frac{1}{3!} \mathrm{ad}\ x(t)^2 - \frac{1}{4!} \mathrm{ad}\ x(t)^3 + \cdots \right)\\\\&\quad\cdot \left( - \frac{1}{2!} \mathrm{ad} x'(t) + \frac{1}{3!} (\mathrm{ad}\ x'(t) \mathrm{ad}\ x(t) + \mathrm{ad}\ x(t) \mathrm{ad}\ x'(t)) - \cdots \right) x'(t) = 0\\\\\end{align}$
I'm not sure what to do now. We could compute a few terms of it to get some terms of $Q_x$, but it seems like there should be a nice formula, like the one above for the differential of $\exp$, or at least there should be a recursive formula for the coefficients of the series...? It's even possible that the formula simplifies WAY down, note eg that all the terms ending with $\mathrm{ad}(x'(t)) x'(t)$ die instantly.
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$.
Best Answer
Using the reference https://arxiv.org/pdf/0903.2087.pdf, which agrees with https://arxiv.org/pdf/hep-th/0001078v1.pdf, which agrees with the reference U. Müller, C. Schubert and Anton M. E. van de Ven, J. Gen. Rel. Grav. 31 (1999) 1759-1768 [arXiv], one sees that the expansion in normal coordinates about a point $p\in M$ of a metric $g$ is given by $$ g_{ij} = g_{ij}(0) + \sum_{n\ge 2} \frac1{n!} C_{ijk_1\cdots k_n}\,x^{k_1}\cdots x^{k_n} $$ where $C_{ijk_1\cdots k_n}$ is the value at $p$ of a tensor of the form $$ -\frac{2(n{-}1)}{n{+}1}\,\nabla_{k_3}\cdots\nabla_{k_n} R_{ik_1jk_2} + LOT_n $$ where $LOT_n$ is a polynomial in the curvature tensor and its derivatives of order strictly less than $n{-}2$. (Note that, in particular, $LOT_2 = 0$.)
Note: The first two references only verify the crucial coefficient $-2(n{-}1)/(n{+}1)$ as above for $n=2,3,4,5$, but, in any case, the coefficients for $n=2$ and $n=3$ are well-known.
Meanwhile, if one does not want to rely on these sources, one can note that these coefficients clearly cannot depend on the dimension of the underlying manifold, so it suffices to verify them in the simplest nontrivial case, i.e., when the dimension of the manifold is $2$. In this case, using the Gauss Lemma, a metric $g$ in geodesic normal coordinates $(x,y)$ centered on $p$ takes the form $$ g = \mathrm{d}x^2 + \mathrm{d}y^2 + h(x,y)\bigl(x\,\mathrm{d}y-y\,\mathrm{d}x)^2, $$ where the function $h$ is arbitrary, subject to the condition that $(x^2{+}y^2)h(x,y)+1>0$.
Letting $r^2 = x^2 + y^2$ and letting $T$ be the radial vector field $x\,\partial_x + y\,\partial_y$, one computes the formula for the Gauss curvature of $g$ to be $$ K = -\frac{2(1+r^2h)(TTh) - r^2(Th)^2+2(5+3r^2h)(Th) + 8r^2h^2+12h}{4(1+r^2h)^2}. $$ Of course, one has the formula for the Riemann curvature tensor in the form $$ R = K\,(1{+}r^2h)\,(\mathrm{d}x\wedge\mathrm{d}y)\otimes(\mathrm{d}x\wedge\mathrm{d}y), $$ while the tensor $(1{+}r^2h)\,(\mathrm{d}x\wedge\mathrm{d}y)\otimes(\mathrm{d}x\wedge\mathrm{d}y)$ is parallel with respect to the Levi-Civita connection. Then, by assuming that $h$ vanishes to order $n{-}2\ge0$, one can easily verify that the above coefficient is correct.