Existence of a local geodesic frame

Let $(M,g)$ be a Riemannian manifold of dimension $n$ with Riemannian connection $\nabla,$ and let $p \in M.$ Show that there exists a neighborhood $U \subset M$ of $p$ and $n$ (smooth) vector fields $E_1,...,E_n \in \chi (U),$ orthonormal at each point of $U,$ s.t. at $p,$ $\nabla_{E_i}E_j(p)=0.$

What I've found so far;

Let $V$ be a normal neighborhood at $p.$ Let $(W,x(x_1,...,x_n))$ be a local coordinate system at $p$ s.t. $x(W) \subset V.$ Let $\{X_i:=\frac{\partial}{\partial x_i}\}$ be the standard basis of $T_pM.$ By Gram-Schmidt, there exists an orthonormal basis $\{E_i\}$ of $T_pM.$ Consider the orthonormal frame $E_i(t),$ by parallel transporting $\{E_i(0)=E_i\}$ along a a geodesic $\gamma :I \to M$ starting at $\gamma(0)=p$ and ending in $\gamma(1)=q,$ where $q \in V.$ This gives us an orthonormal frame at each point of $V.$ To satisfy the second condition, let $E_i(t)=\sum_{l}a_{li}X_l$ and $E_j(t)=\sum_{s}b_{sj}(t)X_s.$ Then

$$\nabla_{E_i}E_j=\sum_{k}(\sum_{l,s} a_{li}(t)b_{sj}(t)\Gamma^k_{ls}+E_i(b_{kj}))X_k,$$

where $\Gamma^k_{is}$ are Christoffel symbols.

To have $\nabla_{E_i}E_j(p)=0,$ we must have $\sum_{l,s}a_{li}(0)b_{sj}(0)\Gamma^k_{ls}+E_i(b_{kj})(p)=0$ for each $k.$

How can I conclude the argument?

Let $U=B_r(p)\subset M^n$ be a normal neighborhood. For each $q\in U$, there is a normalized geodesic $\gamma_q$ joining $p$ with $q$ (radial geodesic). Let $\{v_1,\ldots,v_n\}$ be an orthonormal basis of $T_pM$ and let $\{V_1,\ldots,V_n\}$ be their respective parallel transports along $\gamma_q$. For each $j=1,\ldots,n$, define the field $E_j$ by $$E_j(q)=V_j(d(p,q)),$$ where $d$ is the Riemannian distance. One has that $E_j$ is a $C^{\infty}$ field, because the curves $\gamma_q$ vary $C^{\infty}$ with $q$, in the sense that the EDO's of the geodesics $\gamma_q$ have their coefficients depending $C^{\infty}$ on $q$.

Now, consider $\sigma_i(s)$ the normalized geodesic such that $\sigma_i(0)=p$ and $$\sigma_i'(0)=v_i=V_i(0)=E_i(p).$$ One has, $$\nabla_{E_i}E_j(p)=\nabla_{E_i(p)}E_j=\nabla_{\sigma_i'(0)}E_j=\frac{D(E_j\circ\sigma_i)}{ds}(s)\Big|_{s=0}$$ Since $(E_j\circ\sigma_i)(s)=V_j(d(p,\sigma_i(s)))=V_j(s)$ is a parallel field along $\gamma_{\sigma_i(s)}=\sigma_i\big|_{[0,s]}$, we have that $$\nabla_{E_i}E_j(p)=\frac{D(E_j\circ\sigma_i)}{ds}(s)\Big|_{s=0}=\frac{DV_j}{ds}(0)=0.$$

Let $(U,\phi)$ be a coordinate neighborhood of $p$ and let $g_{ij}$ and $\Gamma_{ij}^{k}$ denote the Riemannian metric tensors and the Christoffel symbols, respectively. If we recall that $g_{ij}=(E_i,E_j)$ and $\Gamma_{ij}^{k}=\nabla_{E_i}E_j(E_k)$ for all $1\leq i,j,k\leq n$ (where $M$ is a smooth $n$-manifold and $E_1,\dots,E_n$ are the coordinate frames on $(U,\phi)$), then we need only have $g_{ij}=\delta_{ij}$ on $U$ and $\Gamma_{ij}^{k}(p)=0$.

You are right in choosing a normal coordinate system at $p$, that is, choosing a normal coordinate neighborhood $(U,\phi)$ at $p$. Let me recall that this means choosing an orthonormal basis $F_1,\dots,F_n$ of the tangent space $T_p(M)$, choosing a star-shaped neighborhood of the origin in $T_p(M)$ which is diffeomorphically mapped onto $U$ under the exponential mapping $\text{exp}_p:D_p\to M$ (where $D_p$ is an open subset of $T_p(M)$ containing $0$) and defining $\phi=\exp_{p}^{-1}$ on $U$ (here we identify $T_p(M)$ with $\mathbb{R}^n$ by the linear isomorphism mapping $F_i$ onto $e_i$, $1\leq i\leq n$, where $e_1,\dots,e_n$ is the standard basis of $\mathbb{R}^n$).

In normal coordinates, the geodesics are of the form $y^i=a^it$ where $a_i$ is a constant for all $1\leq i\leq n$. If you prove this, then you should easily be able to see that $\Gamma_{ij}^{k}(p)=0$ for all $1\leq i,j,k\leq n$ by looking at the second-order ordinary differential equation of geodesics. You can also check that $g_{ij}(p)=\delta_{ij}$ for all $1\leq i,j\leq n$.

I hope this helps!

Choose a normal neighborhood $U\subset M$ containing $p$. Let $\varphi:U\to\mathbb{R}^n$ be a diffeomorphism and $\mathbb{i}:T_pM\to\mathbb{R}^n$ be an isomorphism. Suppose $\{e_i\}_{i=1}^n$ is an orthonormal base in $\mathbb{R}^n$. Then $\widetilde E_i:=\mathbb{i}^{-1}(e_i)$ is a base in $T_pM$. By using the Gram-Schmidt method, we can find an orthonormal base $\{F_i\}_{i=1}^n$ in $T_pM$, such that $g(F_i,F_j)=\delta_{ij},$ and $F_1=\widetilde E_1$.

Now for any $q\in U_p:=U\setminus\{p\}$, there exists a unique $v_q\in T_pM$, such that $\exp_p(v_q)=q$. Let $\gamma_{p,q}(t) :=\exp_p(tv_q)$ and $\widetilde F_i$ be the $g$-parallel vector field along $\gamma_{p,q}$ such that $\widetilde F_i(0)=F_i$. Let $E_i(q):=\widetilde F_i(1)$ and $E_i(p):=F_i$. Then, $\{E_i\}_{i=1}^n$ forms an orthonormal base at each $q\in U$. It is sufficient to prove that $\nabla_{E_i}E_j(p)=0$.

For any $v=\sum_{i=1}^nv_iE_i(p)\in T_pM$. Consider $\gamma_v(t):=\exp_p(tv)$, which is a geodesic with $\dot\gamma_v(t)=\sum_{i=1}^nv_iE_i(\gamma_{p,q}(t))$ since $E_i$ is parallel along each geodesic. Then, we have $\ddot\gamma_v(t)=0$ for any $t\in[0,1]$ by that $v_i$ is a constant and the continuity of $\ddot\gamma_v$ at the endpoints. Then, by the geodesic equation, we have \begin{align*} \sum_{i,j=1}^n\Gamma_{ij}^k(p)v_iv_j= -(\ddot\gamma_v)_k(0)=0,\quad\text{ for any }k\in\{1,\ldots,n\}. \end{align*} Due to the arbitrary of $v\in T_pM$, we can deduce that $\Gamma_{ij}^k(p)=0$ for any $i,j,k\in\{1,\ldots,n\}$, which follows, by the definition of $\Gamma_{ij}^k$, that $g(\nabla_{E_i}E_j(p),E_k(p))=0$ for any $i,j,k\in\{1,\ldots,n\}$, which further implies that $\nabla_{E_i}E_j(p)=0$ for any $i,j\in\{1,\ldots,n\}$.$\square$

In your proof, it is not necessary to use $\{X_i\}$ again in the last paragraph, since $\{X_i\}$ will not give any more information on the Christoffel coefficients.