In most texts on Riemannian geometry that I've seen, (Riemannian) geodesics are first defined as smooth curves $x: [0, 1] \to M$ on a Riemannian manifold $(M, g)$ satisfying $D_t \dot{x} = 0$, where $D_t$ denotes the Levi-Civita connection corresponding to $g$ taken along $x$. It is then shown via variational analysis that critical points of the length/energy functional (let's say energy functional for convenience) $E(x) = \int_0^1 g(\dot{x}, \dot{x})dt$ are geodesics.
The domain of $E$ is often chosen as piecewise smooth curves satisfying some boundary conditions in position. However, one could reasonably define $E$ over curves with less regularity (for example $H^1$ regularity). To that end, let $\Omega \subset H^1([0, 1], M)$ be the space of Sobolev class $H^1$ curves on $M$ satisfying some fixed boundary conditions in position, and reinterpret $E$ as a functional $E: \Omega \to \mathbb{R}$. $\Omega$ is a smooth Hilbert manifold, and the tangent space $T_x \Omega$ consists of Sobolev class $H^1$ vector fields along $x$ which vanish at the endpoints. Hence, given $x \in \Omega$ and $V \in T_x \Omega$, we can calculate $E'(x)V$ by considering a family of curves $x_s: (-\epsilon, \epsilon) \times [0, 1] \to M$ such that $x_s \in \Omega$ for all $s \in (-\epsilon, \epsilon)$, $x_0 = x$, and $\partial_s x_s \vert_{s=0} = V$, and then calculating $\partial_s E(x_s) \vert_{s=0}$. In particular, we find:
\begin{align*} E'(x)V &= \int_0^1 g(D_s \dot{x}_s, \dot{x}_s)\big\vert_{s=0}dt \\ &= \int_0^1 g(D_t V, \dot{x})dt \end{align*}
Which vanishes for all $V \in T_x \Omega$ if $x$ is a critical point of $E$. An article I read (dealing with a problem that is of higher order, but similar in spirit) seems to suggest that one could prove that any critical point of $E$ is smooth via a "bootstrap method," and from there the proof continues as it normally would (by integrating by parts and then choosing $V = D_t \dot{x}$). I'm not familiar with bootstrap methods in this context, could anyone help show me how such an argument may go in this case?
More concretely, how could I show that the $H^1$ regularity of $x$ together with the condition $E'(x)V = \int_a^b g(\dot{x}, D_t V)dt = 0$ for all $V \in T_x \Omega$ implies that $x$ is $H^2$, and inductively that $x$ is smooth?
ADDED: There two key observations.
One is that if $u, u' \in L^1([0,1])$, then $u$ is continuous and satisfies the fundamental theorem of calculus. There are (at least) two ways to show this. One is to apply the fundamental theorem of calculus to a smooth approximation of $f$ and take a limit. The other is to use, for each $\epsilon > 0$, a nonnegative test function $\chi_\epsilon$ that is equal to $1$ on $[a+\epsilon, b-\epsilon]$ and zero at $a$ and $b$. You then show that $$ \lim_{\epsilon\rightarrow 0} \int_0^1 u'\chi_\epsilon\, dt = u(b) - u(a). $$
The second is that, if $V$ is assumed to be smooth and zero at $t= 0, 1$, then in local coordinates, integration by parts implies \begin{align*} 0 &= \int_0^1 g(D_tV,\dot{x})\,dt\\ &= \int_0^1 g_{ij}(\dot{V}^i + \Gamma^i_{pq}V^p\dot{x}^q)\dot{x}^j\,dt\\ &= \int_0^1 -g_{ij}V^j(\ddot{x}^i + \Gamma^i_{pq}\dot{x}^p\dot{x}^q)\,dt \end{align*} Since this holds for any smooth $V$ supported in $(0,1)$, it follows that $x$ is a weak solution to $$ \ddot{x}^i = \Gamma^i_{jk}\dot{x}^j\dot{x}^k. $$ Since the right side is in $L^1$, so is the left. This implies that $$ \dot{x}(t) = \dot{x}(0) + \int_0^t \Gamma^k_{ij}\dot{x}^i\dot{x}^j\,dt. $$ is a continuous function of $t$. This implies that $x$ is a strong solution to the differential equation. Now the argument below applies.
This is covered in most textbooks on ODEs.
The geodesic equation can always be written as a first order system of ODEs, $$ \dot{u} = \Phi(u), $$ by setting $u = (x,\dot{x})$. Integrating this, $u$ satisfies $$ u(t) = u(0) + \int_{s=0}^{s=t} \Phi(u(s))\,ds. $$ Since $u$ is known to be $C^0$, this equation implies that $u$ is $C^1$.
Bootstrapping means applying the above argument to the derivatives of $u$, which themselves are solutions to a system of ODEs obtained by differentiating the original system. For example, $v = (u,\dot{u})$ satisfies $$\dot{v} = \Psi(v),$$ where $\Psi$ is obtained by differentiating the right side of the original system. The argument above shows that $v$ is $C^1$ and therefore $u$ is $C^2$. If $\Phi$ is a smooth function of $u$, then this argument can be continued indefinitely, implying that $u$ is smooth.