For question 2, indeed, you cannot expect the covariant derivatives to vanish on all of $U$. That would imply the full Riemann-curvature tensor vanishes. You’d have to jump through some hoops depending where you start from/how you defined things/what you already know etc, but for instance, using Cartan’s moving frames, this is pretty straight forward since the connection 1-forms by definition vanish, and the curvature is determined by the connection 1-forms. Alternatively, if $\nabla_{E_i}E_j=0$ for all $i,j$ and everywhere on $U$, then in particular, the Lie brackets vanish (since the Levi-Civita connection is torsion-free): $[E_i,E_j]=\nabla_{E_i}E_j-\nabla_{E_j}E_i=0$. So, by Frobenius, there is a local coordinate system $(x^1,\dots, x^n)$ which has $E_i=\frac{\partial}{\partial x^i}$, and hence $g_{ij}=\delta_{ij}$ in these coordinates (expressing orthonormality). From here, we see that we’re locally isometric to $\Bbb{R}^n$, i.e flat.
For the first question, it all boils down to ODE theory. Let’s start with somewhat simpler questions, which you haven’t asked: how do we know geodesics exist? How do we know that the exponential map exists and is smooth? How do we know that we can make sense of parallel transport (i.e why does such a linear map exist between tangent spaces of different points joined by a curve)? The answer to all these questions is that we have existence and uniqueness theorems for ODEs, and smooth dependence on parameters. These results can be found in any sufficiently general ODE textbook/chapters (e.g if I remember correctly, Hirsch and Smale’s ODE book has info about this, Henri Cartan’s Differential Calculus has an entire chapter (i.e the second half of his book) dealing with ODEs. You can also try Dieudonne’s Treatise on Analysis Vol I).
- Geodesics exist because locally in the manifold, you choose coordinates $(x^1,\dots, x^n)$, and write down the geodesic equation $\ddot{x}^i+\Gamma^{i}_{jk}(x)\dot{x}^j\dot{x}^k=0$. Smoothness of $\Gamma$’s (on the chart domain) implies smoothness of $x$ as functions of $t$ (ODE theory).
- You could think of the exponential map as just taking (sufficiently small) tangent vectors and producing points in the manifold by shooting out geodesics for time $1$. Well, why is this mapping smooth? Well, you could consider first the geodesic vector field, $G$, on $TM$, which gives rise to a smooth flow map (ODE theory, with smooth dependence on initial conditions), in this case aptly called the geodesic flow, and then you’re projecting that information from $TM\to M$. This is why $\exp$ is a smooth mapping.
- Parallel transport along a curve $t\mapsto \gamma(t)$ has the same story as above: ODE smoothness theory tells us if we start with a vector $v_0$ in one tangent space, then we can parallel-transport this to produce a mapping $t\mapsto v(t)$ along the curve $\gamma$ and such that $v$ is a smooth function of $t$.
For your question, we’re trying to generalize the third bullet point: rather than just having one fixed curve $t\mapsto \gamma(t)$, I have for each $q\in U$ a curve $t\mapsto \gamma(t,q)$ (a geodesic aimed from $p$ to $q$). In fact, the mapping $(t,q)\mapsto \gamma(t,q)$ is smooth (consequence of second bullet point, which itself is a consequence of ODE theory concerning smoothness with respect to initial conditions). As a result, if I start with a vector $\xi\in T_pM$, then parallel-transporting this along the family of curves $(t,q)\mapsto \gamma(t,q)$ means I have to solve the following ODE (where $(x^1,\dots, x^n,v^1,\dots, v^n)$ are coordinate functions on $TM$)
\begin{align}
\frac{dv^i}{dt}+\Gamma^i_{jk}(\gamma(t,q))\frac{\partial (x^j\circ \gamma)}{\partial t}\bigg|_{(t,q)}v^k&=0
\end{align}
For each point $q\in U$, this becomes a standard ODE, with which we’re satisfied that it has unique smooth solutions. Well, since all the coefficients are smooth functions of $(t,q)$, our ODE theory tells us that the solution $v^i$ will be smooth functions of $(t,q)$ whose initial condition is the original vector $\xi$ we started with (you’ll have to do this $n$ times of course, once for each of the basis vectors of $T_pM$).