[proofplan]
Flatness plus torsion-freeness of the Levi-Civita connection produces a local parallel frame, and the affine-structure argument converts this frame into local coordinates $\tilde x_1, \dots, \tilde x_n$ whose coordinate vector fields are parallel. In these coordinates the metric components $g_{ij}$ are constants, because parallel frames plus metric compatibility force $d g_{ij} = 0$. A constant symmetric positive-definite matrix is diagonalisable by an orthogonal linear change of coordinates, which puts $g$ in Euclidean form.
[/proofplan]
[step:Produce local parallel coordinates from flatness and zero torsion]
Fix $p \in M$. The [Levi-Civita connection](/theorems/1552) $A_{LC}$ is torsion-free by construction, and by hypothesis $F_{LC} \equiv 0$. By the [Flatness Characterises Parallel Frames](/theorems/???) theorem, a flat connection on $TM$ admits, in a neighbourhood of any point, a local frame of covariant-constant sections:
\begin{align*}
s_1, \ldots, s_n : U &\to TM, & A_{LC}\, s_i &= 0, \quad 1 \le i \le n,
\end{align*}
where $U \subset M$ is an open neighbourhood of $p$. By the [affine-structure / parallel-frame-to-coordinates argument](/theorems/???) (torsion-freeness is the input: $[s_i, s_j] = A_{LC}_{s_i} s_j - A_{LC}_{s_j} s_i - \tau_{LC}(s_i, s_j) = 0 - 0 - 0 = 0$, so the $s_i$ pairwise commute and integrate to coordinates), shrinking $U$ if needed there exist coordinates
\begin{align*}
\tilde x_1, \ldots, \tilde x_n : U &\to \mathbb R
\end{align*}
with $\partial_{\tilde x_i} = s_i$ for each $i$. Thus the coordinate fields $\partial_{\tilde x_1}, \dots, \partial_{\tilde x_n}$ are $A_{LC}$-parallel on $U$.
[guided]
Why do flatness and torsion-freeness together yield coordinates whose coordinate fields are parallel?
Flatness produces a parallel *frame* $s_1, \dots, s_n$ by [Flatness Characterises Parallel Frames](/theorems/???) — this is a standard ODE/integrability argument on the bundle $TM$.
To upgrade a frame to coordinates, we need the frame fields to commute pairwise. Commutativity is read off from the fundamental identity relating Lie bracket, connection, and torsion:
\begin{align*}
[X, Y] = A_{LC}_X Y - A_{LC}_Y X - \tau_{LC}(X, Y).
\end{align*}
Plugging $X = s_i, Y = s_j$, both $A_{LC}$-parallel, gives $A_{LC}_{s_i} s_j = 0$ and $A_{LC}_{s_j} s_i = 0$. Torsion-freeness of $A_{LC}$ gives $\tau_{LC} = 0$. So $[s_i, s_j] = 0$.
A commuting frame on $U$ integrates to local coordinates by the classical Frobenius/commuting-flows result: there exist $\tilde x_i$ with $\partial_{\tilde x_i} = s_i$ on a possibly smaller neighbourhood. These coordinate fields are the original parallel frame, so are themselves $A_{LC}$-parallel.
[/guided]
[/step]
[step:Deduce that the metric components in these coordinates are constants]
In the coordinates $(\tilde x_1, \dots, \tilde x_n)$ on $U$, write the metric as
\begin{align*}
g = \sum_{i, j = 1}^n g_{ij}\, d\tilde x_i \otimes d\tilde x_j, \qquad g_{ij} := g(\partial_{\tilde x_i}, \partial_{\tilde x_j}) \in C^\infty(U).
\end{align*}
Fix any $1 \le i, j \le n$ and any $Y \in \Gamma(TU)$. By metric compatibility of the Levi-Civita connection,
\begin{align*}
Y \cdot g_{ij} = Y\cdot g(\partial_{\tilde x_i}, \partial_{\tilde x_j}) = g(A_{LC}_Y \partial_{\tilde x_i}, \partial_{\tilde x_j}) + g(\partial_{\tilde x_i}, A_{LC}_Y \partial_{\tilde x_j}).
\end{align*}
Since $\partial_{\tilde x_i}$ and $\partial_{\tilde x_j}$ are both $A_{LC}$-parallel on $U$ by the previous step, both $A_{LC}_Y \partial_{\tilde x_i}$ and $A_{LC}_Y \partial_{\tilde x_j}$ vanish. Thus $Y \cdot g_{ij} = 0$ for every $Y \in \Gamma(TU)$, which means $d g_{ij} = 0$. Since $U$ is (after shrinking) connected, $g_{ij}$ is a real constant $c_{ij} \in \mathbb R$.
The matrix $(c_{ij}) \in \mathbb R^{n \times n}$ inherits symmetry ($c_{ij} = c_{ji}$ because $g$ is symmetric) and positive-definiteness (because $g_p$ is an inner product at each $p$, including at $p$ itself).
[/step]
[step:Orthogonally diagonalise the constant metric matrix]
Let $C := (c_{ij}) \in \mathbb R^{n \times n}$ be the constant matrix of metric components. By the [Spectral Theorem for symmetric matrices](/theorems/???), there exists an orthogonal matrix $P \in O(n)$ (i.e.\ $P^\top P = I$) such that $P C P^\top$ is diagonal with positive entries $\lambda_1, \dots, \lambda_n > 0$. Composing with a further diagonal rescaling by $\mathrm{diag}(\lambda_1^{-1/2}, \dots, \lambda_n^{-1/2})$ (which preserves the "constant matrix of $g$" property and just changes the diagonal from $(\lambda_i)$ to $(1, \dots, 1)$), we may assume $P C P^\top = I_n$.
Define a new coordinate system
\begin{align*}
y_j : U &\to \mathbb R, & y_j &:= \sum_{i = 1}^n P_{ji}\, \tilde x_i, \quad 1 \le j \le n.
\end{align*}
The Jacobian matrix of the change of coordinates is the constant matrix $P$, which is orthogonal and in particular invertible, so $(y_1, \dots, y_n)$ is a valid coordinate system on $U$ (possibly after shrinking to ensure the coordinate patch condition). The differentials transform by
\begin{align*}
dy_j = \sum_i P_{ji}\, d\tilde x_i.
\end{align*}
[/step]
[step:Verify the metric is Euclidean in the new coordinates]
Substituting,
\begin{align*}
g = \sum_{i, k} c_{ik}\, d\tilde x_i \otimes d\tilde x_k.
\end{align*}
We invert the relation $dy_j = \sum_i P_{ji} d\tilde x_i$ using $P^{-1} = P^\top$: $d\tilde x_i = \sum_j P_{ji}\, dy_j$ (using $P^\top P = I$ so $\sum_j P_{ji} P_{jk} = \delta_{ik}$). Plugging in,
\begin{align*}
g &= \sum_{i, k} c_{ik} \sum_{j, \ell} P_{ji} P_{\ell k}\, dy_j \otimes dy_\ell = \sum_{j, \ell}\Bigl(\sum_{i, k} P_{ji}\, c_{ik}\, P_{\ell k}\Bigr) dy_j \otimes dy_\ell \\
&= \sum_{j, \ell} (P C P^\top)_{j\ell}\, dy_j \otimes dy_\ell = \sum_{j, \ell} \delta_{j\ell}\, dy_j \otimes dy_\ell = \sum_{j = 1}^n dy_j \otimes dy_j.
\end{align*}
Evaluating at the point $p$ gives $g_p = \sum_{j=1}^n dy_j|_p \otimes dy_j|_p$. Renaming $y_j$ as $x_j$ gives local coordinates $(x_1, \dots, x_n)$ near $p$ with $g_p = \sum_{i=1}^n dx_i \otimes dx_i$, establishing the theorem.
[guided]
The final step is a straightforward bilinear change of basis, but it is worth unpacking the algebra.
We have $g = \sum_{i, k} c_{ik}\, d\tilde x_i \otimes d\tilde x_k$ with constant matrix $(c_{ik})$. A linear change of coordinates $y_j = \sum_i P_{ji}\tilde x_i$ with constant orthogonal $P$ produces $dy_j = \sum_i P_{ji}\, d\tilde x_i$, inverted as $d\tilde x_i = \sum_j P_{ji}\, dy_j$ via $P^\top P = I$.
Substituting and collecting,
\begin{align*}
g = \sum_{j,\ell}(PCP^\top)_{j\ell}\, dy_j \otimes dy_\ell.
\end{align*}
We chose $P$ so that $PCP^\top = I_n$ (after the diagonal rescaling that made the eigenvalues all $1$). So the coefficient matrix becomes $\delta_{j\ell}$ and the metric reduces to $\sum_j dy_j \otimes dy_j$, the Euclidean metric.
Note: the result is local and concerns $g_p$ *as a tensor at* $p$ — the statement is actually stronger because our construction gives $g \equiv \sum dy_j \otimes dy_j$ on the whole neighbourhood $U$, not merely at $p$. Flatness is the hypothesis that trivialises the metric, not just diagonalises it at a single point.
[/guided]
[/step]