[proofplan]
We construct the normal-coordinate chart by composing $\exp_p^{-1}$ with a linear isometry $T_pM \cong \mathbb{R}^n$ chosen so that the inner product $g_p$ becomes the standard Euclidean inner product. Property 1 ($p$ at the origin) holds because $\exp_p(0) = p$. Property 2 ($g_{ij}(p) = \delta_{ij}$) holds because $(d\exp_p)_0 = \operatorname{id}$ pushes the Euclidean inner product on $T_0(T_pM)$ to $g_p$ on $T_pM$, and the linear identification was chosen to make $g_p$ Euclidean. Property 3 ($\Gamma_{jk}^i(p) = 0$) is the substantive geometric content: the geodesics through $p$ in normal coordinates are the straight lines $t \mapsto ta$, which substituted into the geodesic equation force a quadratic-form vanishing condition on the symmetric part of the Christoffel symbols. Symmetry of the Levi-Civita connection then upgrades this to vanishing of all $\Gamma_{jk}^i(p)$.
[/proofplan]
[step:Construct the normal-coordinate chart at $p$]
Fix $p \in M$. By the [Exponential Map as a Local Diffeomorphism](/theorems/2712), there exists $\delta > 0$ such that
\begin{align*}
\exp_p|_{B(0,\delta)}: B(0, \delta) \subseteq T_pM &\to U \subseteq M
\end{align*}
is a diffeomorphism onto an open neighbourhood $U$ of $p$ in $M$, where $B(0, \delta)$ is the open ball with respect to the inner product $g_p$ on $T_pM$.
Choose a $g_p$-orthonormal basis $\{e_1, \ldots, e_n\}$ of $T_p M$. This basis induces a linear isometry
\begin{align*}
L: \mathbb{R}^n &\to T_p M \\
(x_1, \ldots, x_n) &\mapsto \sum_{i=1}^n x_i\, e_i,
\end{align*}
where $\mathbb{R}^n$ carries the standard Euclidean inner product $\delta_{ij}$ and $T_pM$ carries $g_p$. Orthonormality of $\{e_i\}$ gives $g_p(L\xi, L\eta) = \sum_i \xi_i \eta_i$ for all $\xi, \eta \in \mathbb{R}^n$, i.e., $L$ pulls $g_p$ back to the Euclidean inner product.
Define the **normal-coordinate chart** $(U, \varphi)$ at $p$ by
\begin{align*}
\varphi: U &\to L^{-1}(B(0, \delta)) \subseteq \mathbb{R}^n \\
q &\mapsto L^{-1}\bigl(\exp_p^{-1}(q)\bigr).
\end{align*}
Since both $\exp_p^{-1}: U \to B(0, \delta)$ and $L^{-1}: T_pM \to \mathbb{R}^n$ are diffeomorphisms (the latter is a linear isomorphism, hence smooth in both directions), $\varphi$ is a diffeomorphism from $U$ onto the open set $L^{-1}(B(0, \delta)) \subseteq \mathbb{R}^n$. Write $(x_1, \ldots, x_n) = \varphi(q)$ for the coordinates of a point $q \in U$.
[/step]
[step:Verify that $p$ has coordinates $(0, \ldots, 0)$]
By construction, $\exp_p(0) = p$, hence $\exp_p^{-1}(p) = 0 \in T_pM$. Then
\begin{align*}
\varphi(p) &= L^{-1}(\exp_p^{-1}(p)) = L^{-1}(0) = (0, \ldots, 0).
\end{align*}
This is property (1).
[/step]
[step:Verify $g_{ij}(p) = \delta_{ij}$ via $(d\exp_p)_0 = \operatorname{id}$]
The coordinate basis vectors $\partial_{x_i}|_p \in T_p M$ in the chart $\varphi$ are, by definition, the images of the standard basis vectors of $\mathbb{R}^n$ under the inverse of the chart's differential:
\begin{align*}
\partial_{x_i}\big|_p &= (d\varphi^{-1})_0(\hat e_i),
\end{align*}
where $\hat e_i$ is the $i$-th standard basis vector of $\mathbb{R}^n$ at the origin (under the identification $T_0 \mathbb{R}^n \cong \mathbb{R}^n$).
We compute $(d\varphi^{-1})_0$. Since $\varphi^{-1} = \exp_p \circ L$, the chain rule and the [Differential of the Exponential Map at the Origin](/theorems/2711) give
\begin{align*}
(d\varphi^{-1})_0 = (d\exp_p)_0 \circ (dL)_0 = \operatorname{id}_{T_pM} \circ L = L,
\end{align*}
where we have used: $(d\exp_p)_0 = \operatorname{id}_{T_pM}$ under the identification $T_0(T_pM) \cong T_pM$; and $(dL)_0 = L$ because $L$ is linear and the differential of a linear map at any point is the linear map itself (under the canonical identifications $T_0 \mathbb{R}^n \cong \mathbb{R}^n$ and $T_0(T_pM) \cong T_pM$).
Hence $\partial_{x_i}|_p = L(\hat e_i) = e_i$, the original orthonormal basis of $T_pM$. The metric in coordinates at $p$ is therefore
\begin{align*}
g_{ij}(p) &= g_p(\partial_{x_i}\big|_p, \partial_{x_j}\big|_p) = g_p(e_i, e_j) = \delta_{ij},
\end{align*}
where the last equality is $g_p$-orthonormality of $\{e_i\}$. This is property (2).
[/step]
[step:Identify radial geodesics with straight lines in normal coordinates]
Let $a \in T_pM$ be such that the maximal geodesic $\gamma_p(\cdot, a)$ is defined on an interval containing $[0, 1]$ and satisfies $a \in B(0, \delta)$. By the [definition of the exponential map](/page/Exponential%20Map), $\exp_p(ta) = \gamma_p(1, ta)$, and by [Geodesic Rescaling](/theorems/2710), $\gamma_p(1, ta) = \gamma_p(t, a)$. So
\begin{align*}
\exp_p(ta) &= \gamma_p(t, a) \qquad \text{for all sufficiently small } t.
\end{align*}
In normal coordinates, the geodesic $\gamma_p(\cdot, a)$ has coordinate representation
\begin{align*}
x(t) &= \varphi(\gamma_p(t, a)) = L^{-1}(\exp_p^{-1}(\exp_p(ta))) = L^{-1}(ta) = t \cdot L^{-1}(a).
\end{align*}
Setting $\xi := L^{-1}(a) \in \mathbb{R}^n$, the geodesic $\gamma_p(\cdot, a)$ is given in normal coordinates by the straight line
\begin{align*}
x(t) &= t \xi, \qquad t \in (-\eta, \eta)
\end{align*}
for $\eta > 0$ small enough that $|t\xi| < \delta$. Every $\xi \in \mathbb{R}^n$ in a neighbourhood of $0$ is realisable in this way (take $a = L\xi$).
[/step]
[step:Substitute straight lines into the geodesic equation to deduce $\Gamma_{jk}^i(p) = 0$]
In the chart $(U, \varphi)$, the geodesic equation reads
\begin{align*}
\ddot x_i(t) + \sum_{j, k = 1}^n \Gamma_{jk}^i(x(t))\, \dot x_j(t)\, \dot x_k(t) &= 0, \qquad i = 1, \ldots, n,
\end{align*}
where $\Gamma_{jk}^i$ are the Christoffel symbols of the Levi-Civita connection in this chart.
For the straight-line geodesic $x(t) = t\xi$, we have $\dot x_j(t) = \xi_j$ and $\ddot x_i(t) = 0$. Substituting and evaluating at $t = 0$ (where $x(0) = 0 = \varphi(p)$):
\begin{align*}
0 + \sum_{j,k=1}^n \Gamma_{jk}^i(0)\, \xi_j\, \xi_k &= 0, \qquad i = 1, \ldots, n.
\end{align*}
That is,
\begin{align*}
\sum_{j,k=1}^n \Gamma_{jk}^i(p)\, \xi_j\, \xi_k &= 0 \qquad \text{for every } \xi \in \mathbb{R}^n \text{ in a neighbourhood of } 0. \tag{$\star$}
\end{align*}
The left-hand side of $(\star)$ is a homogeneous quadratic polynomial in $\xi$ — concretely, $\xi^\top \Gamma^i_{(\cdot,\cdot)}(p)\, \xi$ where $\Gamma^i_{(\cdot,\cdot)}(p)$ is the $n \times n$ matrix with $(j,k)$-entry $\Gamma_{jk}^i(p)$. A homogeneous polynomial that vanishes on an open subset of $\mathbb{R}^n$ vanishes identically (it is determined by its values on any open set), so $(\star)$ holds for **all** $\xi \in \mathbb{R}^n$.
A quadratic form $\xi \mapsto \xi^\top A \xi$ vanishes for all $\xi$ if and only if the symmetric part $\frac{1}{2}(A + A^\top)$ of $A$ is zero. Hence the symmetric part of the matrix $(\Gamma_{jk}^i(p))_{j, k}$ vanishes:
\begin{align*}
\Gamma_{jk}^i(p) + \Gamma_{kj}^i(p) &= 0 \qquad \text{for all } i, j, k \in \{1, \ldots, n\}.
\end{align*}
The Levi-Civita connection is **torsion-free**, which in coordinates means $\Gamma_{jk}^i = \Gamma_{kj}^i$ identically. Combining with the above:
\begin{align*}
2 \Gamma_{jk}^i(p) = \Gamma_{jk}^i(p) + \Gamma_{kj}^i(p) = 0,
\end{align*}
hence $\Gamma_{jk}^i(p) = 0$ for all $i, j, k$. This is property (3), completing the proof.
[guided]
We have, from the previous step, that the radial geodesic $\gamma_p(\cdot, a)$ in normal coordinates is the straight line $x(t) = t\xi$ where $\xi = L^{-1}(a)$. We now extract the vanishing of $\Gamma_{jk}^i(p)$ from this geometric fact. The plan: (i) substitute the straight line into the geodesic equation, (ii) extend the resulting identity from a neighbourhood of $0$ to all of $\mathbb{R}^n$ via polynomial identity, (iii) deduce the symmetric part of $(\Gamma_{jk}^i(p))_{j,k}$ vanishes, and (iv) use torsion-freeness to upgrade this to vanishing of every $\Gamma_{jk}^i(p)$.
In any chart, the geodesic equation for a curve $x: I \to \mathbb{R}^n$ is the second-order ODE system
\begin{align*}
\ddot x_i(t) + \sum_{j, k = 1}^n \Gamma_{jk}^i(x(t))\, \dot x_j(t)\, \dot x_k(t) &= 0, \qquad i = 1, \ldots, n,
\end{align*}
where $\Gamma_{jk}^i$ are the Christoffel symbols of the Levi-Civita connection in the chart $(U, \varphi)$. The bilinear form $\Gamma_{jk}^i(x)\, \dot x_j\, \dot x_k$ is the "covariant correction" to coordinate acceleration: it measures how much the chart's coordinate axes themselves bend, so that a curve with zero coordinate acceleration is geodesic only when the connection coefficients act trivially on its velocity.
Why does this give us information at $p$? Because we already know — from the previous step — a family of geodesics whose coordinate acceleration is identically zero. For the straight line $x(t) = t\xi$ we compute the velocity and acceleration directly:
\begin{align*}
\dot x_j(t) = \xi_j, \qquad \ddot x_i(t) = 0.
\end{align*}
Substituting into the geodesic equation and evaluating at $t = 0$ (where $x(0) = 0 = \varphi(p)$, so $\Gamma_{jk}^i(x(0)) = \Gamma_{jk}^i(p)$):
\begin{align*}
0 + \sum_{j,k=1}^n \Gamma_{jk}^i(p)\, \xi_j\, \xi_k &= 0, \qquad i = 1, \ldots, n.
\end{align*}
That is,
\begin{align*}
\sum_{j,k=1}^n \Gamma_{jk}^i(p)\, \xi_j\, \xi_k &= 0 \qquad \text{for every } \xi \in \mathbb{R}^n \text{ in a neighbourhood of } 0. \tag{$\star$}
\end{align*}
Why only "in a neighbourhood of $0$" and not "for all $\xi \in \mathbb{R}^n$"? Because we only know straight lines $t \mapsto t\xi$ are geodesics for $\xi$ small enough that the line stays inside $L^{-1}(B(0, \delta))$ — i.e., for $\xi$ in some open neighbourhood of $0$. Concretely, as $a$ varies over $B(0, \delta) \subseteq T_pM$, the substitution $\xi = L^{-1}(a)$ produces $\xi$ in the open ball $L^{-1}(B(0, \delta)) \subseteq \mathbb{R}^n$.
We now extend $(\star)$ to all of $\mathbb{R}^n$. The left-hand side, viewed as a function of $\xi$, is a homogeneous quadratic polynomial $Q^i(\xi) = \xi^\top A^i \xi$ where $A^i$ is the $n \times n$ matrix with $(j,k)$-entry $\Gamma_{jk}^i(p)$. A polynomial that vanishes on a non-empty open subset of $\mathbb{R}^n$ vanishes identically: its Taylor coefficients at any interior point of the open set must all be zero, and a polynomial is determined by its Taylor coefficients. Hence $(\star)$ holds for **all** $\xi \in \mathbb{R}^n$:
\begin{align*}
\sum_{j,k=1}^n \Gamma_{jk}^i(p)\, \xi_j\, \xi_k &= 0 \qquad \text{for every } \xi \in \mathbb{R}^n.
\end{align*}
Next, we extract the vanishing of the symmetric part of $A^i = (\Gamma_{jk}^i(p))_{j,k}$. The standard polarisation argument: substituting $\xi + \eta$ for $\xi$ and subtracting,
\begin{align*}
0 = (\xi+\eta)^\top A^i (\xi+\eta) - \xi^\top A^i \xi - \eta^\top A^i \eta = \xi^\top (A^i + (A^i)^\top) \eta \qquad \text{for all } \xi, \eta \in \mathbb{R}^n.
\end{align*}
A bilinear form vanishing on $\mathbb{R}^n \times \mathbb{R}^n$ has zero matrix, so $A^i + (A^i)^\top = 0$, i.e.,
\begin{align*}
\Gamma_{jk}^i(p) + \Gamma_{kj}^i(p) &= 0 \qquad \text{for all } i, j, k \in \{1, \ldots, n\}.
\end{align*}
Why does the quadratic form not detect the antisymmetric part? Because $\xi^\top (A^i - (A^i)^\top) \xi = 0$ for every $\xi$ — antisymmetric matrices are invisible to quadratic forms.
Finally, we close the gap with torsion-freeness. The Levi-Civita connection is **torsion-free**, which in coordinates is the symmetry $\Gamma_{jk}^i = \Gamma_{kj}^i$, valid identically as a function on $U$, and in particular at $p$. Adding the two relations:
\begin{align*}
2 \Gamma_{jk}^i(p) = \Gamma_{jk}^i(p) + \Gamma_{kj}^i(p) = 0,
\end{align*}
hence $\Gamma_{jk}^i(p) = 0$ for all $i, j, k$. This is property (3), completing the proof.
Two remarks on the structure of the argument. First, torsion-freeness was essential: without it, the antisymmetric part of $(\Gamma_{jk}^i(p))_{j,k}$ could be non-zero, and on manifolds with torsion normal coordinates do not in general kill all Christoffel symbols at $p$. Second, $\Gamma_{jk}^i(p) = 0$ does **not** mean $\Gamma_{jk}^i$ vanishes in a neighbourhood of $p$ — it vanishes only at the centre point. The first derivatives $\partial_l \Gamma_{jk}^i(p)$ encode the curvature tensor, so normal coordinates "linearise" the metric at one point but preserve all second-order curvature information.
[/guided]
[/step]