[proofplan]
The strategy is to show that locally, on each small subinterval, $\gamma$ coincides with a smooth radial geodesic from a normal-coordinate chart centred at a point on $\gamma$. Without loss of generality, assume $\gamma$ has unit speed (rescale time). Fix any $t_0 \in [0, 1]$ and let $p_0 := \gamma(t_0)$. Choose $\varepsilon_0 > 0$ such that $\exp_{p_0}\big|_{B(0, \varepsilon_0)}$ is a diffeomorphism onto an open neighbourhood $U_0$ of $p_0$, by [Exponential Map as a Local Diffeomorphism](/theorems/2712). Choose $\delta > 0$ small enough that $\gamma([t_0 - \delta, t_0 + \delta] \cap [0, 1]) \subseteq U_0$ and that the length of $\gamma$ on this subinterval is less than $\varepsilon_0$. The minimality and unit speed of $\gamma$ guarantee $\operatorname{length}(\gamma|_{[t_0, t_0 + \delta]}) = \delta = d(\gamma(t_0), \gamma(t_0 + \delta))$, so $\gamma|_{[t_0, t_0 + \delta]}$ is a length minimizer to a point in $U_0$ at distance $\delta < \varepsilon_0$. Apply [Geodesics Minimize Length Locally](/theorems/2720) part (3): the unique constant-speed length minimizer of length $< \varepsilon_0$ from $p_0$ to its image is the radial geodesic $t \mapsto \exp_{p_0}((t - t_0)\cdot v_0)$ for some $v_0 \in T_{p_0}M$ with $|v_0| = 1$. Hence $\gamma$ coincides with this smooth radial geodesic on $[t_0, t_0 + \delta]$. Since $t_0$ was arbitrary, $\gamma$ is locally a smooth geodesic, hence smooth and a geodesic on $[0, 1]$.
[/proofplan]
[step:Reduce to the case of unit-speed parametrisation]
By hypothesis, $\gamma$ is a piecewise $C^1$ minimal geodesic traversed at constant speed $c \ge 0$. (Here "minimal geodesic" means: $\operatorname{length}(\gamma) = d(\gamma(0), \gamma(1))$, and $\gamma$ minimises length among piecewise $C^1$ curves from $\gamma(0)$ to $\gamma(1)$.)
If $c = 0$, then $\gamma$ is the constant curve at some point, which is $C^\infty$ and a (degenerate) geodesic. Assume henceforth $c > 0$.
Reparametrize: define
\begin{align*}
\bar\gamma: [0, c] &\to M \\
s &\mapsto \gamma(s/c).
\end{align*}
Then $|\dot{\bar\gamma}(s)|_{g_{\bar\gamma(s)}} = (1/c)\, |\dot\gamma(s/c)|_{g_{\gamma(s/c)}} = (1/c) \cdot c = 1$, so $\bar\gamma$ has unit speed. The length and minimality are preserved (length is invariant under reparametrization; the image and minimality are unchanged). It suffices to prove $\bar\gamma$ is $C^\infty$ and a geodesic, since then $\gamma(t) = \bar\gamma(ct)$ is the affine reparametrization, which is also $C^\infty$ and a geodesic (the geodesic equation $\nabla_{\dot\gamma} \dot\gamma = 0$ is preserved under affine reparametrization).
Henceforth, assume $\gamma: [0, L] \to M$ has unit speed, where $L = \operatorname{length}(\gamma) = d(\gamma(0), \gamma(L))$. The interval $[0, 1]$ in the theorem statement is replaced by $[0, L]$ to ease notation; the conclusion is unaffected.
[/step]
[step:For each $t_0 \in [0, L]$, choose a normal-coordinate chart and a small subinterval]
Fix $t_0 \in [0, L)$. We will show $\gamma$ is smooth on a closed interval $[t_0, t_0 + \delta]$ for some $\delta > 0$. Set $p_0 := \gamma(t_0) \in M$.
By [Exponential Map as a Local Diffeomorphism](/theorems/2712), there exists $\varepsilon_0 > 0$ such that
\begin{align*}
\exp_{p_0}\big|_{B(0, \varepsilon_0)}: B(0, \varepsilon_0) \subseteq T_{p_0}M \to U_0 \subseteq M
\end{align*}
is a diffeomorphism onto an open neighbourhood $U_0$ of $p_0$.
By continuity of $\gamma$ and openness of $U_0$, there exists $\delta_1 > 0$ such that $\gamma([t_0, t_0 + \delta_1] \cap [0, L]) \subseteq U_0$. Set
\begin{align*}
\delta := \min\bigl(\delta_1, \, \tfrac{1}{2}\varepsilon_0, \, L - t_0\bigr) > 0.
\end{align*}
Then $[t_0, t_0 + \delta] \subseteq [0, L]$, $\gamma([t_0, t_0 + \delta]) \subseteq U_0$, and the length of $\gamma|_{[t_0, t_0 + \delta]}$ — which equals $\delta$ since $\gamma$ has unit speed — satisfies $\delta \le \tfrac{1}{2}\varepsilon_0 < \varepsilon_0$.
[/step]
[step:Show that $\gamma|_{[t_0, t_0 + \delta]}$ is a length minimizer of length $< \varepsilon_0$ from $p_0$ to $\gamma(t_0 + \delta)$]
Set $q_0 := \gamma(t_0 + \delta) \in U_0$. We show $\operatorname{length}(\gamma|_{[t_0, t_0 + \delta]}) = d(p_0, q_0) = \delta$.
Since $\gamma$ has unit speed,
\begin{align*}
\operatorname{length}(\gamma|_{[t_0, t_0 + \delta]}) = \int_{t_0}^{t_0 + \delta} |\dot\gamma(t)|\, d\mathcal{L}^1(t) = \int_{t_0}^{t_0 + \delta} 1\, d\mathcal{L}^1(t) = \delta.
\end{align*}
Now we show $d(p_0, q_0) = \delta$. The inequality $d(p_0, q_0) \le \delta$ holds because $\gamma|_{[t_0, t_0 + \delta]}$ is a piecewise $C^1$ curve from $p_0$ to $q_0$ of length $\delta$, and $d(p_0, q_0)$ is by definition the infimum of lengths of such curves. For the reverse inequality, suppose for contradiction that $d(p_0, q_0) < \delta$. Then there exists $\tilde\gamma_0 \in \Omega(p_0, q_0)$ (piecewise $C^1$) with $\operatorname{length}(\tilde\gamma_0) < \delta$. Construct a competitor curve from $\gamma(0)$ to $\gamma(L)$ by concatenating:
\begin{align*}
\tilde\gamma := \gamma|_{[0, t_0]} \cdot \tilde\gamma_0 \cdot \gamma|_{[t_0 + \delta, L]},
\end{align*}
where $\cdot$ denotes the concatenation of piecewise $C^1$ curves (with parameter rescaling so the result is defined on $[0, L]$). This is piecewise $C^1$ and goes from $\gamma(0)$ to $\gamma(L)$, with total length
\begin{align*}
\operatorname{length}(\tilde\gamma) = t_0 + \operatorname{length}(\tilde\gamma_0) + (L - t_0 - \delta) < t_0 + \delta + (L - t_0 - \delta) = L.
\end{align*}
But this contradicts $L = d(\gamma(0), \gamma(L))$ (the minimality hypothesis on $\gamma$). Hence $d(p_0, q_0) \ge \delta$, and combined with $\le$, $d(p_0, q_0) = \delta$.
We have established
\begin{align*}
\operatorname{length}(\gamma|_{[t_0, t_0 + \delta]}) = \delta = d(p_0, q_0) < \varepsilon_0.
\end{align*}
So $\gamma|_{[t_0, t_0 + \delta]}$ is a length-minimizing piecewise $C^1$ curve from $p_0$ to $q_0$, with length less than $\varepsilon_0$.
[/step]
[step:Apply the local minimization theorem to identify $\gamma|_{[t_0, t_0 + \delta]}$ as a smooth radial geodesic]
By [Geodesics Minimize Length Locally](/theorems/2720) part (1), applied with $p = p_0$ and $\varepsilon = \varepsilon_0$: there is a unique geodesic $\eta \in \Omega(p_0, q_0)$ with $\operatorname{length}(\eta) < \varepsilon_0$, and $\eta$ is the unique constant-speed length-minimizing curve from $p_0$ to $q_0$.
We must verify that $\gamma|_{[t_0, t_0 + \delta]}$ has constant speed. By construction, $\gamma|_{[t_0, t_0 + \delta]}$ is a unit-speed reparametrisation of itself; since unit speed is constant, it has constant speed.
By the uniqueness clause of part (1) of [Geodesics Minimize Length Locally](/theorems/2720), the constant-speed length-minimising curve from $p_0$ to $q_0$ on $[t_0, t_0 + \delta]$ (after affine reparametrisation to $[0, 1]$) is uniquely determined as the radial geodesic
\begin{align*}
s \mapsto \exp_{p_0}(s \cdot a_0), \quad s \in [0, 1],
\end{align*}
where $a_0 := (\exp_{p_0}\big|_{B(0, \varepsilon_0)})^{-1}(q_0) \in B(0, \varepsilon_0)$. Translating back to the original parametrisation on $[t_0, t_0 + \delta]$,
\begin{align*}
\gamma(t) = \exp_{p_0}\bigl((t - t_0) \cdot v_0\bigr) \quad \text{for } t \in [t_0, t_0 + \delta],
\end{align*}
where $v_0 := a_0/\delta \in T_{p_0}M$ (so that the affine rescaling $s = (t - t_0)/\delta$ converts $s \cdot a_0$ to $(t - t_0)\, v_0$). Note $|v_0|_{g_{p_0}} = |a_0|_{g_{p_0}}/\delta = \delta/\delta = 1$, so $v_0$ is a unit vector.
The map $t \mapsto \exp_{p_0}((t - t_0)v_0)$ is $C^\infty$ on its domain — the exponential map is smooth on its domain of definition by the standard theory of smooth dependence of geodesics on initial data. Hence $\gamma$ is $C^\infty$ on $[t_0, t_0 + \delta]$, and being expressed as a radial geodesic, it satisfies $\nabla_{\dot\gamma}\dot\gamma = 0$ on $(t_0, t_0 + \delta)$. So $\gamma$ is a smooth geodesic on $[t_0, t_0 + \delta]$.
[guided]
This is the heart of the proof. We have built a curve $\gamma|_{[t_0, t_0 + \delta]}$ with three crucial properties: it is piecewise $C^1$, it has length $\delta = d(p_0, q_0)$ (it minimises distance), and that length is strictly less than $\varepsilon_0$ (it lives inside the normal chart). Why does this combination matter? Because the local minimization theorem characterises *uniquely* the constant-speed length-minimisers of length $< \varepsilon_0$ inside a normal chart — they must be radial geodesics. So once we feed our curve into the theorem, its identity is forced.
We invoke [Geodesics Minimize Length Locally](/theorems/2720) part (1) with $p = p_0$ and $\varepsilon = \varepsilon_0$. The theorem requires (i) that $\exp_{p_0}|_{B(0, \varepsilon_0)}$ is a diffeomorphism — guaranteed by our choice of $\varepsilon_0$ in Step 2; and (ii) a target point $q_0$ in the image with $d(p_0, q_0) < \varepsilon_0$ — guaranteed by Step 3, which gave $d(p_0, q_0) = \delta < \varepsilon_0$. The conclusion of part (1): there is a unique geodesic $\eta \in \Omega(p_0, q_0)$ with $\operatorname{length}(\eta) < \varepsilon_0$, and $\eta$ is the unique constant-speed length-minimising curve from $p_0$ to $q_0$ (up to reparametrisation).
Before we can identify $\gamma|_{[t_0, t_0 + \delta]}$ with $\eta$, we must verify that $\gamma|_{[t_0, t_0 + \delta]}$ has constant speed. This is immediate: $\gamma$ has unit speed by the Step 1 reduction, so any restriction of $\gamma$ has unit speed, and unit speed is constant. Why does it minimise length among piecewise $C^1$ curves from $p_0$ to $q_0$? Because Step 3 showed $\operatorname{length}(\gamma|_{[t_0, t_0 + \delta]}) = d(p_0, q_0)$, and $d(p_0, q_0)$ is by definition the infimum of lengths of piecewise $C^1$ curves from $p_0$ to $q_0$ — a curve realising the infimum is a minimiser.
Now we apply the uniqueness clause. The theorem describes the unique minimiser explicitly: in the normal chart $\exp_{p_0}|_{B(0, \varepsilon_0)}$, the constant-speed length-minimising curve from $p_0$ to $q_0$, parametrised on $[0, 1]$, is the radial line in the tangent space pushed forward by the exponential map:
\begin{align*}
s \mapsto \exp_{p_0}(s \cdot a_0), \quad s \in [0, 1],
\end{align*}
where $a_0 := (\exp_{p_0}|_{B(0, \varepsilon_0)})^{-1}(q_0) \in B(0, \varepsilon_0)$ is the unique tangent vector that the exponential map sends to $q_0$. Geometrically, $a_0$ is the "initial velocity, integrated over unit time": following it for time $1$ lands at $q_0$.
Our curve $\gamma$ is parametrised on $[t_0, t_0 + \delta]$, not $[0, 1]$, so we must translate the conclusion back. The affine reparametrisation $s = (t - t_0)/\delta$ converts the displayed formula. Setting $v_0 := a_0 / \delta$, the substitution $s \cdot a_0 = ((t-t_0)/\delta) \cdot a_0 = (t - t_0) v_0$ gives
\begin{align*}
\gamma(t) = \exp_{p_0}\bigl((t - t_0) \cdot v_0\bigr) \quad \text{for } t \in [t_0, t_0 + \delta],
\end{align*}
where $v_0 \in T_{p_0}M$. A sanity check on the speed: $|v_0|_{g_{p_0}} = |a_0|_{g_{p_0}} / \delta$. By the Gauss lemma / radial isometry property baked into the normal chart, $|a_0|_{g_{p_0}} = d(p_0, q_0) = \delta$, so $|v_0|_{g_{p_0}} = 1$. This matches the unit-speed assumption on $\gamma$ — the formula is internally consistent.
Smoothness of $\gamma$ on $[t_0, t_0 + \delta]$ follows because the exponential map $\exp_{p_0}: T_{p_0}M \supseteq \text{(domain)} \to M$ is $C^\infty$ on its domain (this is the standard smooth-dependence-on-initial-data result for the geodesic ODE). The map $t \mapsto (t - t_0) v_0$ from $[t_0, t_0 + \delta]$ to $T_{p_0}M$ is affine, hence $C^\infty$. Composition of $C^\infty$ maps is $C^\infty$, so $\gamma$ is $C^\infty$ on $[t_0, t_0 + \delta]$.
Finally, the geodesic equation $\nabla_{\dot\gamma}\dot\gamma = 0$ on $(t_0, t_0 + \delta)$: any radial curve $t \mapsto \exp_{p_0}((t - t_0) v_0)$ is, by the very definition of the exponential map, the geodesic with initial point $p_0$ and initial velocity $v_0$. So it satisfies the geodesic equation. We conclude that $\gamma$ is a smooth geodesic on $[t_0, t_0 + \delta]$.
[/guided]
[/step]
[step:Patch the local conclusions to obtain $C^\infty$ and the geodesic equation on all of $[0, L]$]
For each $t_0 \in [0, L)$, the previous step gives an open interval $(t_0 - \delta, t_0 + \delta) \cap [0, L]$ (using the same construction with $\delta$ extended to a two-sided interval; the same argument applies for $t_0 \in (0, L]$ on $[t_0 - \delta, t_0]$) on which $\gamma$ is $C^\infty$ and satisfies $\nabla_{\dot\gamma}\dot\gamma = 0$. The collection of these open intervals covers $[0, L]$. Hence $\gamma$ is $C^\infty$ on $[0, L]$ and satisfies the geodesic equation everywhere.
To make the patching explicit: $\gamma$ being $C^\infty$ and a geodesic on each open subinterval of a cover means $\gamma$ is $C^\infty$ on the union $[0, L]$ (smoothness is local) and satisfies the geodesic equation pointwise everywhere on $[0, L]$ (each point lies in some open subinterval where the equation holds, and the equation is pointwise).
The endpoint behaviour: at $t = 0$ and $t = L$, smoothness in the one-sided sense follows from smoothness on a one-sided open neighbourhood by the same argument applied to the one-sided radial geodesic chart. The geodesic equation is satisfied at the endpoints by continuity of $\nabla_{\dot\gamma}\dot\gamma$ (which is continuous because $\gamma$ is $C^\infty$, so $\dot\gamma$ is continuous, and $\nabla$ acts smoothly).
We have proved: the unit-speed reparametrisation $\bar\gamma$ of the original $\gamma$ is $C^\infty$ on its full parameter interval and is a geodesic. Reparametrising back via the affine map $t \mapsto t/c$ (with $c$ the original constant speed), we conclude that the original $\gamma: [0, 1] \to M$ is $C^\infty$ and a geodesic.
[guided]
The patching argument is local-to-global. We have shown that for every parameter $t_0$, $\gamma$ is smooth on a neighbourhood of $t_0$ and satisfies the geodesic equation there. Two facts then complete the proof:
1. **Smoothness is a local property.** A function is $C^\infty$ on a set if and only if it is $C^\infty$ on a neighbourhood of each point in the set. So smoothness on each neighbourhood implies smoothness on the union.
2. **The geodesic equation is pointwise.** At each $t \in [0, L]$, the equation $\nabla_{\dot\gamma}\dot\gamma|_t = 0$ refers to the value of a vector field at a single point. We have established this value is zero on a neighbourhood of $t$, hence at $t$.
Both arguments use that the open neighbourhoods constructed in Step 4 form an open cover of the compact $[0, L]$, but in fact the local-to-global passage does not need compactness — pointwise validity on a neighbourhood of every point is exactly equivalent to pointwise validity everywhere.
A subtle issue is whether the smoothness improvement is automatic given the geodesic equation. The answer is yes: any solution of the geodesic equation, even a priori only $C^2$, is automatically $C^\infty$ by the smoothness of the connection coefficients (the Christoffel symbols are smooth, so a $C^2$ solution of $\ddot x_i + \Gamma_{jk}^i(\gamma)\dot x_j \dot x_k = 0$ inherits regularity from the right-hand side). This is a standard result on smooth dependence in ODE theory.
[/guided]
[/step]