The covariant derivative of the metric is $0$ since we’re using the Levi-Civita connection, so
\begin{align}
k_{\mu\nu}X^{\mu}Y^{\nu}&=(h^{\rho}_{\mu}X^{\mu})(\nabla_{\rho}u_{\sigma})(h^{\sigma}_{\nu}Y^{\nu})= (h^{\rho}_{\mu}X^{\mu})(g_{\sigma\alpha}\nabla_{\rho}u^{\alpha})(h^{\sigma}_{\nu}Y^{\nu}).
\end{align}
In index-free notation, let us write $h(X)=X+g(u,X)u$. Then, the components $h^{\mu}_{\nu}$ match exactly what’s written in equation (6.16) and the above equation becomes
\begin{align}
k(X,Y)&=[h(X)]^{\rho}(g_{\sigma\alpha}\nabla_{\rho}u^{\alpha})[h(Y)]^{\sigma}=g\left(\nabla_{h(X)}u,h(Y)\right).
\end{align}
In particular, if $g(u,X)=0$, i.e $X$ is orthogonal to $u$, then the projection $h(X)$ equals $X$, and likewise for $Y$, so $k(X,Y)=g(\nabla_Xu,Y)$. Hence, the symmetry of $k$ is equivalent to $(6.24)$.
Edit:
In your course, you said you’re using equation (6.21) as the definition for $k$. Let me just warn you that this is not a true ‘definition’, but rather only holds when $u$ is a geodesic vector field (i.e $\nabla_uu=0$, i.e $u$ has vanishing acceleration). Let me thus temporarily define $r_{\mu\nu}=\nabla_{\mu}u_{\nu}$, i.e more explicitly, $r$ is the $(0,2)$-tensor field given by $r:=\nabla(g^{\flat}(u))$. Our goal is to show that $r$ is a symmetric tensor field. We start off with a simple computation (i.e it’s essentially what I did above, but without indices)
\begin{align}
r(X,Y)&:=[\nabla(g^{\flat}(u))](X,Y):=[\nabla_X(g^{\flat}u)](Y)=[g^{\flat}(\nabla_Xu)](Y):=g(\nabla_Xu,Y),\tag{$1$}
\end{align}
where we have used metric-compatibility in the third equal sign. Now, since we assumed $u$ is a unit-timelike geodesic vector field, we have that
- $\nabla_uu=0$ by definition
- $g(u,u)=-1$ is constant, hence for all vector fields $X$, we have $g(\nabla_Xu,u)=0$, or more succinctly, $g(\nabla u,u)=0$.
So, decomposing $X=h(X)+X^{\perp}=h(X)+\xi u$ for some function $\xi$ (namely $\xi=-g(X,u)$), and $Y=h(Y)+\eta u$ for some function $\eta$, we have
\begin{align}
r(X,Y)&=g(\nabla_Xu,Y)\\
&= g(\nabla_Xu,h(Y)+\eta u)\\
&=g(\nabla_Xu,h(Y))+0\tag{since $g(\nabla u,u)=0$}\\
&=g(\nabla_{h(X)}u+0,h(Y))\tag{since $\nabla_{\xi u}u=\xi\nabla_uu=0$}\\
&=k(X,Y).
\end{align}
Hence, now we’re back to the above situation, and the symmetry of $r=k$ follows from (6.24).
Another remark: you may wonder why after obtaining equation $(1)$, $r(X,Y)=g(\nabla_Xu,Y)$, I did not simply apply equation (6.24) and claim that $r$ is a symmetric tensor field. Well, the reason is that equation (6.24) only holds for vector fields $X,Y$ which are orthogonal to $u$. So, this would only prove symmetry of $r$ on the orthogonal complement of $u$. In order to prove the full symmetry of $r$, I had to first show that $r$ is a ‘transverse’ tensor field, i.e $r(X,Y)=r(h(X),h(Y))$ (which equals $k(X,Y)$), i.e depends only on the projections of the vector fields; for this, we needed to use the extra assumption that $u$ is a geodesic vector field, which allowed us to show that $r=k$. This is why I don’t like using equation (6.21) as a definition; it really is a theorem to be proved (the more general definitions are of course given below).
Upon closer reading, your derivation of (**) holds for all $X,Y$ orthogonal to $u$, and so it implies symmetry of $r_{\mu\nu}$ on the orthogonal complement of $u$. If you now use the fact that $r$ is transverse relative to $u$ (e.g as I proved above, or as you quote from Natario’s book… though as you said my $r_{\mu\nu}$ is Natario’s $B_{\nu\mu}$), then in fact we get that $r_{\mu\nu}X^{\mu}Y^{\nu}=r_{\nu\mu}X^{\mu}Y^{\nu}$, for all $X,Y$. In particular, take $X,Y$ to be vector fields with components $X^{\mu}=\delta^{\mu}_{\rho},Y^{\nu}=\delta^{\nu}_{\sigma}$. Then, we get $r_{\rho\sigma}=r_{\sigma\rho}$, which proves the symmetry in index notation… but really knowing $r_{\mu\nu}X^{\mu}Y^{\nu}=r_{\nu\mu}X^{\mu}Y^{\nu}$ for all $X,Y$ means $r(X,Y)=r(Y,X)$, so it directly proves symmetry. This intermediate step of taking special vector fields $X,Y$ is simply to bridge the gap between the ‘symmetry of arguments’ as a $(0,2)$ tensor field, and symmetry of the components, and perhaps this is what you were getting at.
One can generalize the definitions and results much more. The general theme here is that we don’t care so much so about particular vector fields, but rather we care about properties of various subbundles of a given vector bundle (in this case the tangent bundle). For instance, we care about how connections in the big bundle induce connections in (pursuing this line of thought, e.g as I explain here and here gives us the notion of shape tensor or second-fundamental form); in the case of subbundles of the tangent bundle, we are also interested in when the subbundle is Frobenius integrable, etc.
Here are the definitions (written out merely for the sake of completeness):
Definition (Shape tensor, Vorticity).
Let $(M,g)$ be a semi-Riemannian manifold, and $L$ a smooth non-degenerate subbundle of $TM$ (i.e we have an orthogonal direct sum decomposition $TM=L\oplus L^{\perp}$).
Recall that relative to this direct sum decomposition, we can define the shape tensor $\alpha$ on $L$ of the connection $\nabla$ relative to the decomposition $L\oplus L^{\perp}$; so $\alpha:TM\oplus L\to L^{\perp}$ is the bilinear bundle morphism characterized by the fact that for all local vector fields $X$ on $M$ and all local sections $\psi$ of $L$,
$\alpha(X,\psi):= (\nabla_X\psi)_{\perp}$. This is equivalent to saying that for all local vector fields $X$ on $M$, local sections $\psi$ of $L$, and local sections $\zeta$ of $L^{\perp}$, we have $g(\alpha(X,\psi),\zeta):=g(\nabla_X\psi,\zeta)=-g(\psi,\nabla_X\zeta)$.
Let us define the enlarged shape tensor (a made up term… I just didn’t want to abuse notation by calling this $\alpha$ as well) by enlarging the domain of $\alpha$, i.e define $\widetilde{\alpha}:TM\oplus TM\to L^{\perp}$ as $\widetilde{\alpha}(X,Y):=\alpha(X_{\parallel},Y_{\parallel})$. For each $\zeta\in L^{\perp}$, we define the enlarged second-fundamental form along $\zeta$ to be $\mathrm{II}_{\zeta}:TM\oplus TM\to L^{\perp}$ as $\mathrm{II}_{\zeta}(X,Y):=g(\widetilde{\alpha}(X,Y),\zeta)$.
Next, for any local section $\zeta$ of $L^{\perp}$ we define the vorticity of $\zeta$ relative to the decomposition $L\oplus L^{\perp}$ to be the 2-form $\omega_{\zeta}$ on (an open subset of) $M$ defined by
\begin{align}
\omega_{\zeta}(X,Y):=[d(g^{\flat}\zeta)](X_{\parallel}, Y_{\parallel}),
\end{align}
where $X_{\parallel},Y_{\parallel}$ are the projections of $X,Y$ to $L$.
In other words, we convert the vector field into a 1-form, then take its exterior derivative to get a 2-form, and then we only look at $L$-projected inputs.
Now, keep in mind that if $v$ is a 1-form and $X,Y$ are vector fields on $M$, and $\nabla$ is the connection on $TM$, then the coordinate-free definition of the exterior derivative gives
\begin{align}
(dv)(X,Y)&=\mathscr{L}_X(v(Y))-\mathscr{L}_Y(v(X))-v([X,Y])\\
&=\nabla_X(v(Y))-\nabla_Y(v(X))-v([X,Y])\\
&=(\nabla_Xv)(Y)+v(\nabla_XY)-(\nabla_Yv)(X)-v(\nabla_YX)-v([X,Y])\\
&=(\nabla_Xv)(Y)-(\nabla_Yv)(X)+v(\text{Tor}_{\nabla}(X,Y))\\
&= (\nabla_Xv)(Y)-(\nabla_Yv)(X),
\end{align}
where in the second equal sign, we used that the Lie and covariant derivatives are the same when acting on smooth real-valued functions, and in the end we used the fact that $\nabla$ is torsion-free. Now, if the 1-form $v$ is obtained from a vector field $\zeta$ by the metric, i.e $v=g^{\flat}(\zeta)$, then metric compatibility tells us $\nabla_X(g^{\flat}\zeta)=g^{\flat}(\nabla_X\zeta)$, and so the above equation becomes
\begin{align}
d(g^{\flat}(\zeta))(X,Y)&=g^{\flat}(\nabla_X\zeta)(Y)-g^{\flat}(\nabla_Y\zeta)(X)\\
&=g(\nabla_X\zeta,Y)-g(\nabla_Y\zeta,X).
\end{align}
Hence, another way of expressing the vorticity is
\begin{align}
\omega_{\zeta}(X,Y)&:=
d(g^{\flat}(\zeta))(X_{\parallel},Y_{\parallel})\\
&=g\left(\nabla_{X_{\parallel}}\zeta,Y_{\parallel}\right)-g\left(\nabla_{Y_{\parallel}}\zeta,X_{\parallel}\right)\\
&=-g(\alpha(X_{\parallel},Y_{\parallel}),\zeta)+ g(\alpha(Y_{\parallel},X_{\parallel}),\zeta)\\
&=-g(\widetilde{\alpha}(X,Y),\zeta)+g(\widetilde{\alpha}(Y,X),\zeta).\tag{$*$}
\end{align}
Now, here comes the theorem:
Theorem.
Let $(M,g)$ be a semi-Riemannian manifold and $L$ a smooth non-degenerate subbundle. Then, the following statements are equivalent:
- $L$ is Frobenius-integrable.
- for every local section $\zeta$ of $L^{\perp}$, its vorticity $\omega_{\zeta}$ vanishes.
- for every local frame $\{\zeta_i\}_{i=1}^{\text{rank}(L^{\perp})}$ of $L^{\perp}$, the vorticities $\omega_{\zeta_i}$ all vanish.
- the shape tensor $\alpha:TM\oplus L\to L^{\perp}$ when restricted to $L\oplus L$ is symmetric, or equivalently, $\widetilde{\alpha}:TM\oplus TM\to L^{\perp}$ is symmetric.
- for every $\zeta\in L^{\perp}$, the map $\mathrm{II}_{\zeta}:TM\oplus TM\to L^{\perp}$ is symmetric (or equivalently, its restriction to $L\oplus L$ is symmetric).
For the equivalence of $(1)$ and $(2)$, you mimic the proof given in the book. The equivalence of $(2)$ and $(3)$ is pretty simple. The equivalence of $(2)$ and $(4)$ is by formula $(*)$ above, The equivalence of $(4)$ and $(5)$ is also simple.