[proofplan]
Fix $\delta \in (0,b-a)$. We construct the local central field on the compact subarc $[a+\delta,b]$ by foliating a neighbourhood of that subarc by nearby extremals starting from the fixed point $(a,y_0(a))$ with varied initial velocity. The absence of conjugate points makes the parameter-to-position map locally invertible for every $x \in [a+\delta,b]$. The slope field is then obtained by assigning to each nearby point the velocity of the unique member of this extremal family passing through it, and uniqueness for the Euler--Lagrange initial value problem shows that these local definitions agree on overlaps.
[/proofplan]
[step:Convert the Euler--Lagrange equation into a first-order initial value problem]
For $x \in [a,b]$, $y \in U$, and $v \in \mathbb{R}^n$, define the Legendre matrix
\begin{align*}
A(x,y,v) := L_{vv}(x,y,v) \in \mathbb{R}^{n \times n}.
\end{align*}
By the strengthened Legendre condition along $y_0$ and continuity of $A$, after shrinking to an open neighbourhood $\mathcal{U}_0 \subset \mathcal{O}$ of the compact set
\begin{align*}
K_0 := \{(x,y_0(x),y_0'(x)) : x \in [a,b]\},
\end{align*}
we may assume that $A(x,y,v)$ is invertible for every $(x,y,v) \in \mathcal{U}_0$. Since $\mathcal{U}_0$ is open and contains $K_0$, the first-order formulation below is available on a slightly larger time interval near every point of $[a,b]$, including a right-sided extension past $b$ along nearby solutions.
The Euler--Lagrange equation can then be solved for $y''$. Indeed,
\begin{align*}\frac{d}{dx}\partial_v L(x,y,y') = L_{xv}(x,y,y') + L_{yv}(x,y,y')y' + L_{vv}(x,y,y')y'',\end{align*}
so the equation is equivalent on $\mathcal{U}_0$ to
\begin{align*}
y'' = F(x,y,y'),
\end{align*}
where the map
\begin{align*}
F: \mathcal{U}_0 &\to \mathbb{R}^n
\end{align*}
is defined by
\begin{align*}
F(x,y,v) := A(x,y,v)^{-1}\bigl(\partial_y L(x,y,v)-L_{xv}(x,y,v)-L_{yv}(x,y,v)v\bigr).
\end{align*}
Here $L_{xv}(x,y,v)$ denotes the $x$-derivative of the map $\partial_v L$, and $L_{yv}(x,y,v)v$ denotes the derivative of the map $y \mapsto \partial_v L(x,y,v)$ applied to the vector $v \in \mathbb{R}^n$. Since $L$ is $C^3$ and inversion is $C^1$ on the [open set](/page/Open%20Set) of invertible matrices, the map $F$ is $C^1$ on $\mathcal{U}_0$. Therefore the first-order system for $(y,w) \in U\times\mathbb{R}^n$,
\begin{align*}
y'=w, \qquad w'=F(x,y,w),
\end{align*}
has the standard local existence, uniqueness, and $C^1$ dependence on initial data, by the $C^1$ form of the [Picard-Lindelof Theorem](/theorems/69), because its phase-space vector field $(x,y,w)\mapsto (w,F(x,y,w))$ is $C^1$ on the open set $\mathcal{U}_0$. We use that local theory in the following form: for initial data sufficiently close to $(x_*,y_0(x_*),y_0'(x_*))$, there is a unique nearby phase-space solution $(Y,W)$ depending $C^1$ on the initial data, where $W=\partial_xY$ and the curve $x\mapsto Y(x)$ is an extremal.
[/step]
[step:Construct the fixed-initial-point family of extremals]
Apply the local existence part of the [Picard-Lindelof Theorem](/theorems/69) at the terminal phase-space point $(b,y_0(b),y_0'(b))$. Since $\mathcal{U}_0$ is open and the vector field is $C^1$ there, after decreasing $\eta_0>0$ if necessary, the reference solution extends uniquely from $[a,b]$ to $[a,b+\eta_0]$ while remaining in $\mathcal{U}_0$ near $b$. We keep the notation $Y(x,0)=y_0(x)$ and $W(x,0)=y_0'(x)$ for this extended reference solution.
For $s>0$, let $B(0,s)\subset\mathbb{R}^n$ denote the open Euclidean ball of radius $s$ centred at $0$. By $C^1$ dependence on initial data in the same theorem, compactness of $[a,b+\eta]$ for any fixed $\eta\in(0,\eta_0]$, and openness of $\mathcal{U}_0$, choose $\eta\in(0,\eta_0]$ and $\rho>0$ small enough that, for every $\alpha \in B(0,\rho)\subset\mathbb{R}^n$, the phase-space initial value problem
\begin{align*}Y(a,\alpha) = y_0(a), \qquad W(a,\alpha)=y_0'(a)+\alpha\end{align*}
has a unique solution on $[a,b+\eta]$ remaining in the neighbourhood where the first-order formulation is valid. Define the phase-space flow components
\begin{align*}
Y: [a,b+\eta]\times B(0,\rho) \to U
\end{align*}
and
\begin{align*}
W: [a,b+\eta]\times B(0,\rho) \to \mathbb{R}^n
\end{align*}
by assigning to $(x,\alpha)$ the position and velocity at time $x$ of this solution.
By construction, $Y(x,0)=y_0(x)$ and $W(x,0)=y_0'(x)$ for every $x \in [a,b]$. The initial-value dependence theorem gives $Y,W \in C^1$ on $[a,b+\eta]\times B(0,\rho)$, and for each fixed $\alpha$, the curve $x \mapsto Y(x,\alpha)$ is an extremal of $L$ with $\partial_xY(x,\alpha)=W(x,\alpha)$.
Let
\begin{align*}
J: [a,b] &\to \mathbb{R}^{n\times n}
\end{align*}
be the parameter derivative along the central extremal,
\begin{align*}
J(x) := \partial_\alpha Y(x,0).
\end{align*}
The standard linearized-flow equation for a $C^1$ ODE flow says that each column of $J$ satisfies the variational equation obtained by differentiating the Euler--Lagrange system with respect to the initial-velocity parameter; this is the Jacobi equation in the displayed form, consistent with the local [Existence and Uniqueness of Jacobi Fields](/theorems/3524). The hypotheses needed for this differentiation hold here: $L$ is $C^3$, $L_{vv}$ is invertible on $\mathcal{U}_0$, and the first-order flow depends $C^1$ on $\alpha$. The initial conditions are obtained by differentiating the defining initial conditions for $Y$ with respect to $\alpha$ at $\alpha=0$:
\begin{align*}
J(a)=0, \qquad J'(a)=I_n.
\end{align*}
Thus $J=Z$, the Jacobi matrix in the statement.
[/step]
[step:Invert the fixed-initial-point family away from $a$]
Fix $x_* \in (a,b]$. Since $J(x_*)=Z(x_*)$ and $\det Z(x_*)\neq 0$, the derivative
\begin{align*}
\partial_\alpha Y(x_*,0): \mathbb{R}^n \to \mathbb{R}^n
\end{align*}
is an isomorphism. By the [Inverse Function Theorem](/theorems/51) applied to the $C^1$ map $(x,\alpha)\mapsto (x,Y(x,\alpha))$ on the open time interval $(a,b+\eta)$, there exist an open interval $\widetilde I_{x_*}\subset (a,b+\eta)$ containing $x_*$, an open set $V_{x_*}\subset U$ containing $y_0(x_*)$, and an open parameter set $A_{x_*}\subset B(0,\rho)$ containing $0$ such that, with $I_{x_*}:=\widetilde I_{x_*}\cap(a,b]$,
\begin{align*}
\Phi_{x_*}: I_{x_*}\times A_{x_*} \to I_{x_*}\times V_{x_*}, \qquad \Phi_{x_*}(x,\alpha):=(x,Y(x,\alpha))
\end{align*}
is a $C^1$ diffeomorphism onto its image. Since this image is relatively open in $(a,b]\times U$ and contains $(x_*,y_0(x_*))$, choose a relatively open neighbourhood
\begin{align*}
\mathcal{N}_{x_*}\subset \Phi_{x_*}(I_{x_*}\times A_{x_*})
\end{align*}
of $(x_*,y_0(x_*))$.
Let
\begin{align*}
\alpha_{x_*}: \mathcal{N}_{x_*} &\to A_{x_*}
\end{align*}
denote the second component of $\Phi_{x_*}^{-1}$, so that
\begin{align*}
Y(x,\alpha_{x_*}(x,y))=y
\end{align*}
for every $(x,y)\in \mathcal{N}_{x_*}$. Define the local slope field
\begin{align*}
P_{x_*}: \mathcal{N}_{x_*} &\to \mathbb{R}^n
\end{align*}
by
\begin{align*}
P_{x_*}(x,y):=W(x,\alpha_{x_*}(x,y)).
\end{align*}
Then $P_{x_*}$ is $C^1$ because both $W$ and $\alpha_{x_*}$ are $C^1$. Since $\alpha_{x_*}(x,y_0(x))=0$ for $(x,y_0(x))$ in this chart,
\begin{align*}
P_{x_*}(x,y_0(x))=W(x,0)=y_0'(x).
\end{align*}
[guided]
The goal is to turn the family of extremals into a genuine field of directions. For that, a nearby point $(x,y)$ must lie on exactly one member of the family $Y(\cdot,\alpha)$, so that we can define the field value at $(x,y)$ as the velocity of that member.
Fix $x_* \in (a,b]$. The relevant map is
\begin{align*}
\alpha \mapsto Y(x_*,\alpha).
\end{align*}
Its derivative at $\alpha=0$ is
\begin{align*}
\partial_\alpha Y(x_*,0)=J(x_*)=Z(x_*).
\end{align*}
The no-conjugate-point hypothesis says $\det Z(x_*)\neq 0$, hence this derivative is an isomorphism of $\mathbb{R}^n$. Therefore the [Inverse Function Theorem](/theorems/51) applies to the map
\begin{align*}
\widetilde\Phi_{x_*}: \widetilde I_{x_*}\times A_{x_*} \to \widetilde I_{x_*}\times V_{x_*}, \qquad \widetilde\Phi_{x_*}(x,\alpha):=(x,Y(x,\alpha)),
\end{align*}
after choosing an open interval $\widetilde I_{x_*}\subset(a,b+\eta)$, an open parameter set $A_{x_*}$, and an open set $V_{x_*}$ sufficiently small. The derivative of this product map at $(x_*,0)$ is block triangular with identity in the $x$-coordinate and $\partial_\alpha Y(x_*,0)$ in the vertical block, so it is invertible. For the theorem we then restrict to $I_{x_*}:=\widetilde I_{x_*}\cap(a,b]$, which also covers the endpoint $x_*=b$ from the left.
The [Inverse Function Theorem](/theorems/51) gives a diffeomorphism onto an open image, not automatically a whole product $I_{x_*}\times V_{x_*}$. We therefore choose the actual domain of the local field to be a relatively open neighbourhood
\begin{align*}
\mathcal{N}_{x_*}\subset \Phi_{x_*}(I_{x_*}\times A_{x_*})
\end{align*}
of $(x_*,y_0(x_*))$. This choice guarantees that every point of $\mathcal{N}_{x_*}$ is represented by one and only one nearby parameter.
Let
\begin{align*}
\alpha_{x_*}: \mathcal{N}_{x_*} &\to A_{x_*}
\end{align*}
be the parameter selected by the inverse map. It is characterized by
\begin{align*}
Y(x,\alpha_{x_*}(x,y))=y.
\end{align*}
Now define
\begin{align*}
P_{x_*}: \mathcal{N}_{x_*} &\to \mathbb{R}^n
\end{align*}
by
\begin{align*}
P_{x_*}(x,y):=W(x,\alpha_{x_*}(x,y)).
\end{align*}
This is the only possible definition compatible with the extremal family: if $(x,y)$ lies on the member $Y(\cdot,\alpha)$, the field must assign the velocity $W(x,\alpha)=\partial_xY(x,\alpha)$.
The regularity follows from the phase-space construction. The velocity component $W$ is $C^1$ by $C^1$ dependence of the first-order flow, the inverse parameter map $\alpha_{x_*}$ is $C^1$ by the [Inverse Function Theorem](/theorems/51), and therefore $P_{x_*}$ is $C^1$. Finally, on the original curve $y=y_0(x)$, the corresponding parameter is $\alpha=0$ because $Y(x,0)=y_0(x)$. Hence
\begin{align*}
P_{x_*}(x,y_0(x))=W(x,0)=y_0'(x).
\end{align*}
[/guided]
[/step]
[step:Show the fixed-initial-point slope fields agree on overlaps]
Fix the compact subarc $K_\delta:=[a+\delta,b]$. For an invertible matrix $M\in\mathbb{R}^{n\times n}$, let $s_{\min}(M)$ denote its least singular value, equivalently
\begin{align*}
s_{\min}(M):=\inf\{|Mh|:h\in\mathbb{R}^n, |h|=1\}.
\end{align*}
Since $x\mapsto Z(x)$ is continuous and $\det Z(x)\neq 0$ on $K_\delta$, the least singular value function $x\mapsto s_{\min}(Z(x))$ has a positive minimum; define
\begin{align*}
\sigma_\delta := \min_{x\in K_\delta} s_{\min}(Z(x)) > 0.
\end{align*}
Shrink the neighbourhoods $\mathcal{N}_{x_*}$, if necessary, so that their inverse parameters take values in a common ball $B(0,r)\subset B(0,\rho)$. Here $\|\cdot\|_{\mathrm{op}}$ denotes the Euclidean operator norm. Shrink $\rho$ and the neighbourhoods once more, if necessary, so that
\begin{align*}
\|\partial_\alpha Y(x,\alpha)-Z(x)\|_{\mathrm{op}}\leq \frac{\sigma_\delta}{2}
\end{align*}
for every $x\in K_\delta$ and every $\alpha\in B(0,r)$. This is possible by continuity of $\partial_\alpha Y$ and compactness of $K_\delta$. Let $\mathcal{L}^1$ denote one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure). For $\alpha_1,\alpha_2\in B(0,r)$, the [fundamental theorem of calculus](/theorems/632) applied to the line segment $t\mapsto \alpha_2+t(\alpha_1-\alpha_2)$ gives
\begin{align*}
Y(x,\alpha_1)-Y(x,\alpha_2)=\int_0^1 \partial_\alpha Y(x,\alpha_2+t(\alpha_1-\alpha_2))(\alpha_1-\alpha_2)\,d\mathcal{L}^1(t).
\end{align*}
The preceding operator-norm bound and the definition of $\sigma_\delta$ imply that the average [linear map](/page/Linear%20Map) in this integral has lower bound at least $\sigma_\delta/2$ in the direction $\alpha_1-\alpha_2$. Hence $Y(x,\alpha_1)=Y(x,\alpha_2)$ forces $\alpha_1=\alpha_2$. Thus the maps $\alpha\mapsto Y(x,\alpha)$ are injective on $B(0,r)$ for every relevant $x\in K_\delta$.
Suppose $x_*,x_{**}\in K_\delta$ and $(x,y)$ belongs to the overlap $\mathcal{N}_{x_*}\cap\mathcal{N}_{x_{**}}$. Let $\alpha_1:=\alpha_{x_*}(x,y)$ and $\alpha_2:=\alpha_{x_{**}}(x,y)$. By the defining property of the inverse charts,
\begin{align*}
Y(x,\alpha_1)=y=Y(x,\alpha_2).
\end{align*}
Both parameters lie in the injectivity ball $B(0,r)$, so $\alpha_1=\alpha_2$. Therefore the selected extremal is the same curve, and the velocity component of the phase-space flow gives
\begin{align*}
P_{x_*}(x,y)=W(x,\alpha_1)=W(x,\alpha_2)=P_{x_{**}}(x,y).
\end{align*}
[/step]
[step:Assemble the central field on a neighbourhood of the separated subarc]
Fix $\delta \in (0,b-a)$. The compact set
\begin{align*}
\Gamma_{y_0,\delta}=\{(x,y_0(x)):x\in[a+\delta,b]\}
\end{align*}
is covered by finitely many neighbourhoods $\mathcal{N}_{x_1},\dots,\mathcal{N}_{x_m}$ from the away-from-endpoint construction. Define the relatively open set
\begin{align*}
\mathcal{N}_\delta:=\left(\mathcal{N}_{x_1}\cup\cdots\cup\mathcal{N}_{x_m}\right)\cap\left([a+\delta,b]\times U\right).
\end{align*}
This set is a relatively open neighbourhood of $\Gamma_{y_0,\delta}$ in $[a+\delta,b]\times U$. By the overlap agreement proved above, the formula
\begin{align*}
P_\delta(x,y)=P_{x_j}(x,y)
\end{align*}
on $\mathcal{N}_\delta\cap\mathcal{N}_{x_j}$ is well-defined. Since each local field is $C^1$ and the definitions agree on overlaps, the restricted map $P_\delta:\mathcal{N}_\delta\to\mathbb{R}^n$ is $C^1$ in the relative sense.
For every $x\in[a+\delta,b]$, the point $(x,y_0(x))$ lies in one of the local charts, and in each chart the local field satisfies
\begin{align*}
P_\delta(x,y_0(x))=y_0'(x).
\end{align*}
It remains to verify that an integral curve of $P_\delta$ is one of the constructed extremals locally. Let $y:I\to U$ be a $C^1$ integral curve of $y'=P_\delta(x,y)$, and restrict to a subinterval $I_0\subset I$ on which $(x,y(x))\in\mathcal{N}_{x_*}$. Define $\alpha:I_0\to A_{x_*}$ by $\alpha(x):=\alpha_{x_*}(x,y(x))$. Then $Y(x,\alpha(x))=y(x)$ for $x\in I_0$. Differentiating this identity gives
\begin{align*}
\partial_xY(x,\alpha(x))+\partial_\alpha Y(x,\alpha(x))\alpha'(x)=y'(x).
\end{align*}
Since $\partial_xY(x,\alpha(x))=W(x,\alpha(x))=P_{x_*}(x,y(x))=y'(x)$, we obtain
\begin{align*}
\partial_\alpha Y(x,\alpha(x))\alpha'(x)=0.
\end{align*}
The chart was chosen so that $\partial_\alpha Y(x,\alpha(x))$ is invertible, hence $\alpha'(x)=0$ on $I_0$. Thus $\alpha$ is constant on $I_0$, and $y(x)=Y(x,\alpha_0)$ there for a fixed parameter $\alpha_0$. Therefore $y$ locally satisfies the Euler--Lagrange equation, and being an extremal is local in $x$. Every integral curve of $P_\delta$ is consequently an extremal of $L$. Thus $P_\delta$ is a central field on the neighbourhood $\mathcal{N}_\delta$ of the separated graph segment $\Gamma_{y_0,\delta}$, completing the proof.
[/step]