Let $(M,g)$ be a Lorentzian manifold. Denote $g_{ij}$ to be the metric components and $(g^{-1})^{ij}$ to be the metric components of its inverse. Let $u$ and $\underline{u}$ satisfy eikonal equations, i.e. $$(g^{-1})^{\mu\nu}\partial_{\mu}u\partial_{\nu}u=0,\quad(g^{-1})^{\mu\nu}\partial_{\mu}\underline{u}\partial_{\nu}\underline{u}=0.$$
Define the vector fields $L'$ and $\underline{L'}$ by
$$ L'^{\mu}=-2(g^{-1})^{\mu\nu}\partial_{\nu}u,\quad \underline{L}'^{\mu}=-2(g^{-1})^{\mu\nu}\partial_{\nu}\underline{u}.$$
Then,
$$\nabla_{L'}L'=0, \quad\nabla_{\underline{L}'}\underline{L}'=0.$$
Then, there is a computation that I don know how to varify:
in arbitrary coordinates,
$$g_{\lambda\mu}L'^{\nu}\nabla_{\nu}L'^{\mu}=-2L'^{\nu}\nabla_{\nu}\partial_{\lambda}u=-2L'^{\nu}\nabla_{\lambda}\partial_{\nu}u=4g^{\nu k}\partial_{k}u\nabla_{\lambda}\partial_{\nu}u=2\partial_{\lambda}(g^{\nu k}\partial_{\nu}u\partial_{k}u)=0.$$
Can anyone please explain how to show this computation? Thanks a lot.
The first equality follows from the fact that: $$ \nabla_\nu L'^{\mu} = -2 g^{\mu\alpha}\nabla_\nu (d u)_\alpha = -2g^{\mu \alpha} \text{Hess}(u)_{\nu\alpha}. $$ Some notation writes this as $\text{Hess}(u)_{\nu\alpha} = \nabla_\nu \partial_\alpha u$. Then, contraction with the metric factor $g_{\lambda \mu}$ gives the first equality you have: $$ g_{\lambda \mu}\nabla_\nu L'^{\mu} = -2 g_{\lambda\mu}g^{\mu\alpha}\nabla_\nu \partial_\alpha u = -2\nabla_\nu\partial_\lambda u. $$ The next equality follows from the symmetry of the Hessian, so the indices $\nu$ and $\lambda$ can be interchanged.
The third equality uses the definition of $L'$ once more, the fourth uses the Leibniz rule, and the last the fact that $g^{\lambda k}\partial_\lambda u \partial_k u$ is (in particular) constant.
Edit: The fact that $L'^{\mu} = -2g^{\mu\nu}\partial_\nu u$ is just the coordinate-notation way to write that $L' = -2\text{grad}(u)$. Here, grad means the gradient defined by the Lorentz metric $g$, i.e. it is the metric dual to the one-form $du$. In general the gradient of a function $f$ is given by $g^{\mu\nu}\partial_\nu f$. Really this is the expression for the $\mu$th coordinate of $\text{grad}(f)$ in the coordinate basis $\partial_\mu$.
In general for a vector field $X$, $\nabla_\mu X^\nu$ means the $\nu$th coordinate of $\nabla_\mu X$ (which is a vector field, and is shorthand for $\nabla_{\partial_\mu} X$). More accurately it would be written as $(\nabla_\mu X)^\nu$. This also happens to be the (lower) $\mu$th and (upper) $\nu$th coordinate of $\nabla X$, which is written as $(\nabla X)_\mu^{\ \nu}$.
For a covector field $\omega$, $\nabla_\mu\omega_\alpha$ is the $\alpha$th component of $\nabla_\mu \omega$; it would more accurately be written $(\nabla_\mu \omega)_\alpha$. It is also the $\mu, \alpha$th component of $\nabla\omega$. But here one has to be careful in the use of indices; does $(\nabla\omega)_{\alpha\beta}$ equal $(\nabla_{\partial_\alpha} \omega)(\partial_\beta)$, or does it equal $(\nabla_{\partial_\beta}\omega)(\partial_\alpha)$? This is a convention and depends on who is writing. In the case of this question, $\omega = du$ is the differential of a function $u$, and $\omega_\alpha = \partial_\alpha u$. I wrote $\nabla_\nu (du)_\alpha$ to emphasize that it is the $\alpha$th component of $\nabla_\nu du$, as you mentioned; when $\nabla_\nu \partial_\alpha u$ is written, it usually means the same thing, but I get confused because it looks like it could also be $\nabla_\nu(\partial_\alpha u)$, which is just the second partial derivative $\partial^2_{\nu\alpha} u$.