Your definition of a Killing vector field $X$ implies that
\begin{equation}
g_t(Y,Z)=g_0(Y,Z)+O(t^2)\\\lim_{t \to 0}\frac{g_t(Y,Z)-g_0(Y,Z)}{t}=0
\end{equation}
the LHS is the definition of the Lie derivative i.e.
\begin{equation}
(L_Xg)(Y,Z)=0
\end{equation}
for any vector fields $Y$ and $Z$. So
\begin{equation}
(L_Xg)(X,Y)=0\\
Xg(X,Y)-g(L_XX,Y)-g(X,L_XY)=0\\
g(\nabla _XX,Y)+g(X,\nabla _XY)-g([X,X],Y)-g(X,[X,Y])=0\\
g(\nabla _XX,Y)+g(X,\nabla _YX)=0
\end{equation}
If the integral curve of $X$ is a geodesic, then $\nabla _XX=0$ and so
\begin{equation}
g(X,\nabla _YX)=0\\
0.5Yg(X,X)=0
\end{equation}
Then $X$ has constant length. Conversely, if $X$ has constant length then $Yg(X,X)=0$ and so $g(\nabla _XX,Y)=0$ for any vector field $Y$ and so $\nabla _XX=0$. the result.
To show $\Psi_t$ of $V$ is an isometry for each $t$, that is,
\begin{equation*}
g_{\Psi_t(x)}(D_x\Psi_t X, D_x\Psi_t Y)=g_x(X,Y)\quad\forall\,\,X,Y\in T_p M
\end{equation*}
Instead, we shall show
\begin{equation*}
\frac{\partial}{\partial t}( X_i, X_j) =0,
\end{equation*}
which is the same as showing
\begin{equation*}
\frac{\partial}{\partial t}(D_x\Psi_t X, D_x\Psi_t Y)=\text{const}\quad\dagger
\end{equation*}
Proof:
$V$ is a killing field, then we have for every $p\in M$ and $X,Y\in T_p M$.
\begin{equation*}
g(\nabla_X V,Y)+g(X,\nabla_Y V)=0
\end{equation*}
Write $Y=V_i=\frac{\partial}{\partial x_i}$,\quad $X= V_j=\frac{\partial}{\partial x_j}$. Since the partial derivative commute in $\mathbb{R}^n$, that is,
\begin{equation*}
[V_i,V_j]=0\quad\forall\,\,i,j\,\,\text{and}\,\,[V, V_i]=[V_n,V_i]=0
\end{equation*}
It suffices to show that $V$ is a Killing field if and only $L_V g(V_i,V_j)=0$ for all $i,j$. First note that local flows for any given time $t$ always give local diffeomorphisms, since their inverse is provided by the local flow of the vector field $-V$. Thus, for $V$ to be a Killing field near $x$ is equivalent to have the property $(u,v)=((D_x\Psi_t u, D_x\Psi_t v))$ for all $u,v\in T_p M$, and all $p$ near $x$. We claim moreover that for fixed $x_1,\cdots,x_{n-1}$, as $t$ varies the coefficients of $D_x\Psi_t u$ and $D_x\Psi_t v$ in term of the $V_i$ are constant.
Then
\begin{align*}
&\ g(\nabla_{V_j}V,V_i)+g(\nabla_{V_i}V,V_j)\\
= &\ g([V_j,V]-\nabla_V V_j, V_i)+g([V_i,V]-\nabla_V V_i,V_j)\\
= &\ ([V_j,V],V_i)-(\nabla_V V_j,V_i)+([V_i,V],V_j)-(\nabla_V V_i,V_j)\\
= &\ -g([V,V_j],V_i)-g([V,V_i],V_j)-{\color{blue}g((\nabla_V,V_j,V_i)+(\nabla_V V_i,V_j))}\\
= &\ -([V,V_j],V_i)-([V,V_i],V_j)-{\color{blue}g(\nabla_V V_j\cdot V_i+ V_j\nabla_V V_i)}\\
= &\
-([V,V_j],V_i)-([V,V_i],V_j)-V(V_j,V_i)\\
= &\
-\frac{\partial}{\partial t}g(V_j,V_i)\\
= &\
-\frac{\partial}{\partial t}g(D_x\Psi_t V_j, D_x\Psi_t V_i)
\end{align*}
The last two lines indicate $\Psi_t$ of $V$ is an isometry of $(M,g)$ for each $t$.
I always feel the minus sigh looks suspicious…
Any correction is appreciated.
Best Answer
I think that the main point of the proof i.e. the reason that $\Phi_s^* g = g$ is missing in your proof. We can fill that gap by observing that it follows from the uniqueness of solutions of the equation:
$\left.\frac{d}{ds}\right|_{s = 0} \Phi_s^* g = 0 $
with initial condition $\Phi_0^* g = g$.