[guided]The purpose of this step is to prove that the nonlinear part is genuinely second order in the error, uniformly in time.
Fix $v\in B(0,\rho)$ and $t\in I$. Introduce the line segment map
\begin{align*}
\gamma_{v,t}:[0,1]\to V
\end{align*}
by
\begin{align*}
\gamma_{v,t}(s):=\hat{x}(t)+sv.
\end{align*}
For every $s\in[0,1]$ we have $|sv|\le |v|<\rho$, so $\gamma_{v,t}(s)\in\mathcal{T}_{\rho}\subset V$. This verifies the domain condition for applying Taylor's formula along the whole segment.
For the state equation, consider the function
\begin{align*}
\phi_{v,t}:[0,1]\to\mathbb{R}^n
\end{align*}
defined by
\begin{align*}
\phi_{v,t}(s):=f(\gamma_{v,t}(s),u(t)).
\end{align*}
The parameter $u(t)$ is fixed while $s$ varies. Since $f$ is $C^2$, the map $\phi_{v,t}$ is $C^2$. Its first derivative is
\begin{align*}
\phi_{v,t}'(0)=J_x f_{(\hat{x}(t),u(t))}v=A(t)v.
\end{align*}
Its second derivative is obtained by differentiating again in the state direction $v$:
\begin{align*}
\phi_{v,t}''(s)=D_x^2 f_{(\gamma_{v,t}(s),u(t))}[v,v].
\end{align*}
By the definition of $L_f$ and by the operator norm bound for bilinear maps,
\begin{align*}
|\phi_{v,t}''(s)|\le L_f|v|^2
\end{align*}
for all $s\in[0,1]$. Taylor's formula with quadratic remainder for $\phi_{v,t}$ at $s=0$ gives
\begin{align*}
R_f(v,t)=\phi_{v,t}(1)-\phi_{v,t}(0)-\phi_{v,t}'(0).
\end{align*}
The remainder estimate is therefore
\begin{align*}
|R_f(v,t)|\le \frac{1}{2}L_f|v|^2.
\end{align*}
The same argument applies to the output map. Define
\begin{align*}
\psi_{v,t}:[0,1]\to\mathbb{R}^p
\end{align*}
by
\begin{align*}
\psi_{v,t}(s):=h(\gamma_{v,t}(s)).
\end{align*}
Then $\psi_{v,t}$ is $C^2$,
\begin{align*}
\psi_{v,t}'(0)=Jh_{\hat{x}(t)}v=C(t)v,
\end{align*}
and
\begin{align*}
\psi_{v,t}''(s)=D^2h_{\gamma_{v,t}(s)}[v,v].
\end{align*}
Since $\gamma_{v,t}(s)\in\mathcal{T}_{\rho}$, the definition of $L_h$ gives
\begin{align*}
|\psi_{v,t}''(s)|\le L_h|v|^2.
\end{align*}
Taylor's formula gives
\begin{align*}
|R_h(v,t)|\le \frac{1}{2}L_h|v|^2.
\end{align*}
Now combine the two bounds. Since
\begin{align*}
r(v,t)=R_f(v,t)-K(t)R_h(v,t),
\end{align*}
the triangle inequality and the induced operator norm give
\begin{align*}
|r(v,t)|\le |R_f(v,t)|+\|K(t)\|_{\mathrm{op}}|R_h(v,t)|.
\end{align*}
Using $\|K(t)\|_{\mathrm{op}}\le K_0$, we obtain
\begin{align*}
|r(v,t)|\le \frac{1}{2}\bigl(L_f+K_0L_h\bigr)|v|^2.
\end{align*}
Thus, with
\begin{align*}
M:=\frac{1}{2}\bigl(L_f+K_0L_h\bigr),
\end{align*}
or with $M:=1$ in the degenerate case where this number is $0$, the desired uniform quadratic estimate holds for every $v\in B(0,\rho)$ and every $t\in I$.[/guided]