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.
I found on the final of the page $14$ of this thesis how to prove that $\nabla_i X_j + \nabla_j X_i = (\mathscr{L}_X g)_{ij}$, which is more clear than in Wald's book. I will put here the development:
$\textit{Proof.}$ Let $\omega$ be the $1$-form dual to the vector field $X$, $\omega(Y) = \langle X, Y \rangle$. Using the product rule (from Lemma 1.6) and the metric compatibility and torsion-free conditions on the Levi-Civita connection we have
\begin{align*}
\mathscr{L}_X g(Y,Z) &= X(g(Y,Z)) - g(\mathscr{L}_X Y, Z) - g(Y, \mathscr{L}_X Z)\\
&= \langle \nabla_X Y, Z \rangle + \langle Y, \nabla_X Z \rangle - \langle [X, Y], Z \rangle - \langle Y, [X, Z] \rangle\\
&= \langle \nabla_X Y - [X,Y], Z \rangle + \langle Y, \nabla_X Z - [X,Z] \rangle\\
&= \langle \nabla_Y X, Z \rangle + \langle Y, \nabla_Z X \rangle\\
&= Y \langle X, Z \rangle - \langle X, \nabla_Y Z \rangle + Z \langle Y, X \rangle - \langle \nabla_Z Y, X \rangle\\
&= Y(\omega(Z)) - \omega(\nabla_Y Z) + Z(\omega(Y)) - \omega(\nabla_Z Y)\\
&= (\nabla_Y \omega) (Z) + (\nabla_Z \omega) (Y),
\end{align*}
which is the coordinate-free way of expressing the identity we wanted. Note that we use the product rule again to get the last line. $\square$
Best Answer
In the paper they include also the following:
This is indeed how one proves the Poincare lemma for 1 one form. Indeed you are always dealing with the one form $X = X_i dx^i$. Since
$$ (dX)_{ij} = \partial _i X_j - \partial _j X_i = \nabla _j X_i - \nabla _i X_j = 0,$$
$X$ is a closed one form and thus $X = df$ for some smooth function $f$.