Schwarzschild geometry in Schwarzschild coordinates $(t,r,\theta,\phi)$ is time-symmetric
\begin{equation}
ds^2=-\left(1-\frac{2GM}{c^2 r}\right)c^2dt^2+\left(1-\frac{2GM}{c^2 r}\right)^{-1}dr^2+r^2\left(d \theta^2 +\sin^2 \theta d \phi^2\right) \;.
\end{equation}
Novikov coordinate system is defined by a set of geodesic clocks. The coordinate clocks are freely falling from some maximal radius $r_m$ towards $r=0$, where $r_m$ is different for each clock. All clocks start falling at the same Schwarzschild time $t_0$ and they are synchronized in such a manner that each clock shows $0$ at $r_m$. Novikov coordinate is defined to stay constant along the trajectory of each clock, while for time coordinate proper time is taken.
From now on the angular part metric will be omitted, since is stays the same. We also take $r_s=2M$ and $G=c=1$:
\begin{equation}\label{eq:sch-met2}
ds^2=-\left(1-\frac{r_s}{r}\right)dt^2+\left(1-\frac{r_s}{r}\right)^{-1}dr^2 \;.
\end{equation}
Geodesics in Schwarzschild gometry
To get the equation of geodesics in Schwarzschild geometry we have to solve equations of motion of a free particle:
\begin{equation}\label{eq:lagrangian}
\mathcal{L}=\frac{1}{2}mg_{\mu\nu}\dot{x}^\mu\dot{x}^\nu \;,
\end{equation}
\begin{equation}\label{eq:dot}
\dot{x}^\mu=\frac{dx^\mu}{d\tau}=u^\mu \;.
\end{equation}
\begin{equation}\label{eq:lagrangian2}
\mathcal{L}=-\frac{m}{2}\left(1-\frac{r_s}{r}\right)\dot{t}^2+\left(1-\frac{r_s}{r}\right)^{-1}\dot{r}^2 \;,
\end{equation}
\begin{equation}\label{eq:EL}
\frac{d}{d\tau}\frac{\partial\mathcal{L}}{\partial \dot{x}^\mu}-\frac{\partial\mathcal{L}}{\partial x^\mu}=0 \;,
\end{equation}
For $\mu=0$ we get a constant of motion
\begin{equation}\label{eq:ConstOfMotion}
\frac{\partial}{\partial\tau}\left[\left(1-\frac{r_s}{r}\right)\dot{t}\right]=0 \qquad \Rightarrow \qquad \left(1-\frac{r_s}{r}\right)\dot{t}=a \;,
\end{equation}
For timelike geodesics: $ds^2=-d\tau^2$ the radial geodesic equation becomes
\begin{equation}\label{eq:orbit}
\left(\frac{d\tau}{dr}\right)^2=\frac{1}{a^2-\left(1-\frac{r_s}{r}\right)} \;.
\end{equation}
Maximal radius is ($dr/d\tau=0$)
\begin{equation}\label{eq:maximal}
r_m=\frac{r_s}{1-a^2} \;.
\end{equation}
We use $\frac{dt}{dr}=\frac{dt}{d\tau}\frac{d\tau}{dr}$ and obtain the following relations:
\begin{eqnarray}
\frac{d\tau}{dr} &=& \frac{\varepsilon}{\sqrt{\frac{r_s}{r}-\frac{r_s}{r_m}}} \;,\label{eq:orbit1} \
\frac{dt}{dr} &=& \frac{\varepsilon\sqrt{1-\frac{r_s}{r_m}}}{\left(1-\frac{r_s}{r}\right)\sqrt{\frac{r_s}{r}-\frac{r_s}{r_m}}} \;, \label{eq:orbit2}
\end{eqnarray}
where $\varepsilon$ is $+1$ or $-1$. For falling particles we choose $\varepsilon=-1$.
Novikov time coordinate
We first transform from $(r,t)$ to $(r,\tau)$.
From last two equations we obtain for $d\tau(dt,dr)$
\begin{equation}\label{eq:bit}
d\tau=\left(1-\frac{r_s}{r_m}\right)^{1/2}dt+\frac{\left(\frac{r_s}{r}-\frac{r_s}{r_m}\right)^{1/2}}{1-\frac{r_s}{r}}dr \;.
\end{equation}
where we assumed $t$ are $r$ known.
This can be integratet from $r$ to $r_m$, where we take into account that all clocks reach their maximum radius at $\tau_{0i}=0$. It follows
\begin{equation}\label{eq:integral}
\tau=\left(1-\frac{r_s}{r_m}\right)^{1/2}(t-t_0)+\int_{r_m}^{r}\frac{\left(\frac{r_s}{y}-\frac{r_s}{r_m}\right)^{1/2}}{1-\frac{r_s}{y}}dy \;.
\end{equation}
maximal radius $r_m$ is here a function of $r$ and \tau$. Their implicit relationship is
\begin{equation}\label{eq:implicit1}
\tau=-f(r,r_m)\;,
\end{equation}
where
\begin{equation}
f(r,r_m) = \int_{r_m}^{r}\frac{dy}{\sqrt{\frac{r_s}{y}-\frac{r_s}{r_m}}} \label{eq:integral3}
= -\left[\frac{rr_m}{r_s}(r_m-r)\right]^{1/2}-\frac{r_m^{3/2}}{\sqrt{r_s}}\arccos\left[\left(\frac{r}{r_m}\right)^{1/2}\right] \;.\label{eq:f}
\end{equation}
We can now eliminate coordinate $t$ from the line element
\begin{equation}\label{eq:sch-met3}
ds^2=-d\tau^2+\frac{1}{1-\frac{r_s}{r_m}}\left[- dr-\left(\frac{r_s}{r}-\frac{r_s}{r_m}\right)^{1/2}d\tau\right]^2 \;.
\end{equation}
Novikov radial coordinate
For radial coordinate we take the maximal Schwarzschild radius $r_m$, which remains constant along the worldline of a geodesic clock.
\begin{equation}\label{eq:relation2}
- dr-\left(\frac{r_s}{r}-\frac{r_s}{r_m}\right)^{1/2}d\tau=\left(\frac{r_s}{r}-\frac{r_s}{r_m}\right)^{1/2}\frac{\partial f}{\partial r_m}dr_m \;.
\end{equation}
With this we can eliminate the other Schwrazschild coordinate $r$:
\begin{equation}\label{eq:sch-met4}
ds^2=-d\tau^2+\frac{\left[g(r,r_m)\right]^2}{1-\frac{r_s}{r_m}}dr_m^2 \;.
\end{equation}
Here we $g(r,r_m)$ is the following
\begin{eqnarray}
g(r,r_m)&=&-\left(\frac{r_s}{r}-\frac{r_s}{r_m}\right)^{1/2}\frac{\partial f}{\partial r_m} \label{eq:g}\
&=&1+\frac{1}{2}\left(1-\frac{r}{r_m}\right)-\frac{3}{4}\left(\frac{r_m}{r}-1\right)^{1/2}\left[\sin^{-1}\left(\frac{2r}{r_m}-1\right)-\frac{\pi}{2}\right] \;. \nonumber
\end{eqnarray}
$r$ is not a radial coordinate anymore, but a metric function of coordinates $r_m$ and $\tau$, which is given implicitly by equation ().
Novikov metric
By introducting $r_m$ the metric became diagonal as in Schwarzschild coordinates. It also stays diagonal by introducig a new radial coordinate, that is only functionally related to the old one. Novikov-s choice is $r^*$ with the following monotonic relation to $r_m$:
\begin{equation}\label{eq:r*}
r^*=\left(\frac{r_m}{r_s}-1\right)^{1/2}\;.
\end{equation}
The metric now becomes
\begin{equation}\label{eq:novikov-met}
ds^2=-d\tau^2+4r_s^2\left({r^*}^2+1\right)\left[g(r,r^*)\right]^2d{r^*}^2 \;.
\end{equation}
We can show that the following also holds
\begin{equation}\label{eq:relation}
4Mg(r,r^*)=\frac{1}{r^*}\frac{\partial r}{\partial r^*}\;.
\end{equation}
With this the metric gets the form that is standard in literature [MTW,p. 826]:
\begin{equation}\label{eq:novikov}
ds^2=-d\tau^2+\left(\frac{{r^*}^2+1}{{r^*}^2}\right)\left(\frac{\partial r}{\partial r^*}\right)^2d{r^*}^2+r^2\left(d \theta^2 +\sin^2 \theta d \phi^2\right) \;,
\end{equation}
where we also included the angular part.
Relations among coordinates
We now give the relations between Schwarzschild coordinates $(t,r)$ and Novikov coordinates $(\tau,r^*)$. The first one, $r=(\tau,r^*)$, is obtained from equations () and ()
\begin{equation}\label{eq:CoordRela1}
\tau=r_s\left({r^*}^2+1\right)\left[\frac{r}{r_s}-\frac{(r/r_s)^2}{{r^*}^2+1}\right]^{1/2}+r_s\left({r^*}^2+1\right)^{3/2}\arccos\left[\left(\frac{r/r_s}{{r^*}^2+1}\right)^{1/2}\right] \;.
\end{equation}
The second one, $t=(\tau,r^*)$, is obtained by integration from ()
\begin{equation}
t=r_s\ln\left|\frac{r^*+\left(\frac{r_s}{r}\left({r^}^2+1\right)-1\right)^{1/2}}{r^-\left(\frac{r_s}{r}\left({r^*}^2+1\right)-1\right)^{1/2}}\right|+r_sr^\left[\left({r^}^2+3\right)\arctan\left(\frac{r_s}{r}\left({r^}^2+1\right)-1\right)^{1/2}+\left({r^}^2+1\right)\frac{\sqrt{\frac{r_s}{r}\left({r^*}^2+1\right)-1}}{\frac{r_s}{r}\left({r^*}^2+1\right)}\right] \;.
\end{equation}
The problem here is with assuming that because $n_{\alpha}\sim \delta_{t\alpha}$, then the $n_{\alpha;\beta}$ will also vanish for $\alpha \neq t$. This is not true. Even though a field might have only some component, it does not mean that its covariant derivative does not have any other components.
Let us see, what hides behind the $n_{\alpha;\beta}$.
In semicolon notation, we have to be careful and remember that we first take the covariant derivative of the tensor field in the $\frac{\partial}{\partial x^{\beta}}$ direction, and only then evaluate its $\alpha$-th component.
Therefore:
$$n_{\alpha;\beta}\equiv(\nabla_{\beta}n)_{\alpha}$$
where the parentheses serve to indicate the order of operations. It's analogous to how we first differentiate a function and then take its value at a point. In the opposite order, we'd always get zero, like you're getting here!
I encourage you to try and calculate the covariant derivative for
$\beta=(t,r,\theta,\varphi)$ like this:
$$\nabla_{\beta}n = \Big{[} \frac{\partial n_{\alpha}}{\partial x^{\beta}} - \Gamma_{\beta \alpha}^{\lambda}n_{\lambda} \Big{]}dx^{\alpha}.$$
And then read off the relevant entries $(\nabla_{\beta}n)_{\alpha}$ of the one-form field $\nabla_{\beta}n$ .
Alternatively, since you already have $n$, we can do it like this:
$$ \nabla_{\beta}n = \nabla_{\beta}(n_{t}dt) = \nabla_{\beta}(n_{t})dt+ n_{t}\nabla_{\beta}dt = \partial_{\beta}n_{t}dt - n_{t}\Gamma_{\beta \lambda}^{t}dx^{\lambda} $$
where we used, accordingly, the Leibniz rule and the fact that on functions, the covariant derivative reduces to a partial derivative, and we know its actions on coordinate one-forms and coordinate basis vector fields.
Second problem I see is that in the question's body, you normalized incorrectly the one-form $n$. You should use the $g^{tt}$ that's given below.
For $n_{\alpha}$, I am getting:
$$ n = -\frac{\sqrt{2m+r}}{(2mr+r^{2})\sqrt{r}}\sqrt{4j^{2}\sin(\theta)^{2} +r^{4} -4m^{2}r^{2}} \;dt$$
After an edit, the answer below can be now treated as illustrating an alternative way of calculating the extrinsic curvature tensor:
I allowed myself the liberty to run calculations partly through SageMath.
$$g_{\mu\nu} = \begin{bmatrix}
-(1-\frac{2M}{r}) & 0 & 0 & \frac{-2j\sin^{2}\theta}{r} \\
0 & 1+\frac{2M}{r} & 0 & 0 \\
0 & 0 & r^{2}(1+\frac{2M}{r}) & 0\\
\frac{-2j\sin^{2}\theta}{r} & 0& 0& r^{2}(1+\frac{2M}{r})\sin^{2}\theta
\end{bmatrix}
$$
Inverting:
$$g^{\mu\nu} = \begin{bmatrix}
\frac{r^{3}(2M+r)}{4M^{2}r^{2}-r^{4}-4j^{2}\sin^{2}\theta} & 0 & 0 & \frac{2jr}{4M^{2}r^{2}-r^{4}-4j^{2}\sin^{2}\theta} \\
0 & \frac{r}{2M+r} & 0 & 0 \\
0 & 0 & \frac{1}{2Mr+r^{2}} & 0\\
\frac{2jr}{4M^{2}r^{2}-r^{4}-4j^{2}\sin^{2}\theta} & 0& 0& - \frac{r(2M-r)\,\sin\theta^{-2}}{4M^{2}r^{2}-r^{4} -4j^{2}\sin^{2}\theta}
\end{bmatrix}
$$
You can calculate the inverse easily if you treat the $t-\varphi$ block of the metric as 2-by-2 and calculate its inverse separately, and the other block to be $r-\theta$ and obtain the corresponding components of the inverse metric simply by raising the entries of that block to the power of $-1$.
Now you should calculate $n_{\mu}= At_{,\mu}$, normalize it: $g^{\mu\nu}n_{\mu}n_{\nu}=-1$ to get the scaling factor $A$ and finally, raise the index to get the normal vector field instead of the normal form.
You should get that the normal vector field has two nonvanishing components - in the $t$ and $\varphi$ coordinates.
$$ n = n^{t}\partial_{t} + n^{\varphi}\partial_{\varphi}$$
Then you can use the formula:
$$ K_{ij}= -\frac{1}{2}(\mathcal{L}_{n}\gamma)_{ij},$$
where $\mathcal{L}$ is the Lie derivative and $\gamma_{ij}$ is the induced metric on a $t=$const slice.
Best Answer
I don't think the way you take the covariant derivatives is correct. I will show that, in fact, there is no time derivative of the lapse in the trace of extrinsic curvature.
We'll be using the four-dimensional point of view, in which fields living on $\mathbb{R} \times \Sigma$ are treated as defined on the whole spacetime manifold $\mathcal{M}$.
The extensions of tensor fields on $\Sigma$ to $\mathcal{M}$ we call spatial, as their contractions with the normal vector field $\mathbf{n}$ vanish, and the action of the orthogonal projector $\gamma^{\mu}_{\;\;\nu}$ on these fields is an identity. For example, $\beta$ is defined to be a spatial vector field that quantifies the difference between normal observers and coordinate time observers:
$$ \partial_{t} = \alpha n + \beta, \\[0.5em] g(\beta, n) = 0.$$
From the latter, knowing that $n_{\mu}=(-\alpha, 0,0,0)$, we can derive:
$$ 0= \beta^{\mu}n_{\mu} = \beta^{t}n_{t} + \beta^{i}n_{i}= \beta^{t}n_{t}.$$
And therefore the temporal component of the shift vector is identically zero $$ \beta^{t}=0.$$
CONTENT OF THE EDIT:
Contract $K_{\mu\nu}=-\nabla_{\nu}n_{\mu} -a_{\mu}n_{\nu}$ with the inverse 4-metric:
$$ K=K_{\mu\nu}g^{\mu\nu} = -g^{\mu\nu}\nabla_{\nu}n_{\mu} = \nabla_{\mu}n^{\mu} = -\partial_{\mu}n^{\mu} - \Gamma^{\mu}_{\mu\lambda}n^{\lambda}.$$
Here $a_{\mu}n_{\nu}g^{\mu\nu}=0$, since $a_{\mu}=n^{\lambda}\nabla_{\lambda}n^{\mu}.$
Now let us introduce the notation: $g$, $\gamma$ are the determinants of the 4 and 3-metrics, respectively, satisfying the relation $\sqrt{-g}=\alpha \sqrt{\gamma}$ (ref 1, eq. 4.55) and we will use the fact that:
$$\Gamma^{\mu}_{\mu\lambda}=\frac{1}{2}g^{\mu\nu}g_{\mu\nu,\lambda}=\frac{1}{\sqrt{-g}}\frac{\partial \sqrt{-g}}{\partial x^{\lambda}}=\frac{1}{\alpha\sqrt{\gamma}}\frac{\partial (\alpha\sqrt{\gamma})}{\partial x^{\lambda}}, $$
as well as:
$$ n^{\mu}=(\frac{1}{\alpha}, \frac{-\beta^{i}}{\alpha})$$
Expanding the expression of interest:
$$-\partial_{\mu}n^{\mu} - \Gamma^{\mu}_{\mu\lambda}n^{\lambda} = -\partial_{t}n^{t} - \partial_{i}n^{i} - \frac{1}{\alpha\sqrt{\gamma}}\partial_{\lambda}(\alpha\sqrt{\gamma}) n^{\lambda} =\\ \frac{1}{\alpha^{2}}\partial_{t}\alpha - \partial_{i}n^{i} -\frac{n^{\lambda}}{\sqrt{\gamma}}\partial_{\lambda}\sqrt{\gamma} -\frac{1}{\alpha}n^{\lambda}\partial_{\lambda}\alpha = \\ \frac{1}{\alpha^{2}}\partial_{t}\alpha - \partial_{i}n^{i} - \frac{n^{\lambda}}{\sqrt{\gamma}}\partial_{\lambda}\sqrt{\gamma} - \frac{1}{\alpha}n^{t}\partial_{t}\alpha - \frac{1}{\alpha}n^{i}\partial_{i}\alpha=\\ \frac{1}{\alpha^{2}}\partial_{t}\alpha - \partial_{i}n^{i} - \frac{n^{\lambda}}{\sqrt{\gamma}}\partial_{\lambda}\sqrt{\gamma} - \frac{1}{\alpha^{2}}\partial_{t}\alpha - \frac{1}{\alpha}n^{i}\partial_{i}\alpha$$ As you see, the time derivatives of the lapse cancel out. Going further: $$ -\partial_{i}n^{i} - \frac{n^{\lambda}}{\sqrt{\gamma}}\partial_{\lambda}\sqrt{\gamma} - \frac{1}{\alpha}n^{i}\partial_{i}\alpha =\\ \partial_{i}(\frac{\beta^{i}}{\alpha}) - \frac{n^{t}}{\sqrt{\gamma}}\partial_{t}\sqrt{\gamma} - \frac{n^{i}}{\sqrt{\gamma}}\partial_{i}\sqrt{\gamma} + \frac{1}{\alpha^{2}}\beta^{i}\partial_{i}\alpha=\\ \frac{1}{\alpha}\partial_{i}\beta^{i} - \frac{1}{\alpha^{2}}\beta^{i}\partial_{i}\alpha - \frac{1}{\alpha\sqrt{\gamma}}\partial_{t}\sqrt{\gamma} + \frac{\beta^{i}}{\alpha\sqrt{\gamma}}\partial_{i}\sqrt{\gamma} + \frac{1}{\alpha^{2}}\beta^{i}\partial_{i}\alpha=\\ \frac{-1}{\alpha}\Big{[} \frac{1}{\sqrt{\gamma}}\partial_{t}\sqrt{\gamma} - \frac{1}{\sqrt{\gamma}}\partial_{i}(\sqrt{\gamma}\beta^{i})\Big{]}=\\ \frac{-1}{\alpha}\Big{[} \frac{1}{\sqrt{\gamma}}\partial_{t}\sqrt{\gamma} - D_{i}\beta^{i} \Big{]}$$
Where the last expression involving the shift vector is just an alternative formula for the divergence of a vector field.
We can get to the expression we have just derived by using an an alternative (but equivalent) formula for the extrinsic curvature tensor.
$$ K_{ij} = \frac{-1}{2\alpha}(\mathcal{L}_{m}\gamma)_{ij}=\frac{-1}{2\alpha}\big{[}\partial_{t}\gamma_{ij} - (\mathcal{L}_{\beta}\gamma)_{ij} \big{]}.$$
Upon contraction and using $(\mathcal{L}_{\beta}\gamma)_{ij}= D_{i}\beta_{j} + D_{j}\beta_{j}$:
$$ K=K_{ij}\gamma^{ij} = \frac{-1}{\alpha}\Big{(} \frac{1}{\sqrt{\gamma}} \frac{\partial \sqrt{\gamma}}{\partial t} - D_{i}\beta^{i}\Big{)} .$$
An easier derivation could be made from the start, based on the fact that $\mathbf{K}$ is a spatial tensor, and that $g^{\mu\nu}= \gamma^{\mu\nu} + n^{\mu}n^{\nu}$.
Contracting as before:
$$ K=K_{\mu\nu}g^{\mu\nu} = K_{\mu\nu}\gamma^{\mu\nu}=K_{ij}\gamma^{ij} = -\gamma^{ij}\nabla_{j}n_{i}$$
and once again, no time derivative of the lapse is present in the end.
In the above we use the fact that $\gamma^{t\mu}=0$ is zero, as it is a spatial tensor with an upper temporal index $\gamma^{\mu\nu}n_{\nu}=0$. Additionally, we use the fact that $n_{i}=0$ (but not the lowered-index version covariant derivative of its covariant derivative, $(\nabla_{j}n)_{i}$, mind you!
This $\gamma$ is, by the way, the extension of the 3-metric $\gamma_{ij}$ to a spatial tensor field in 4-dimensional spacetime.
As for the derivations of different slicings:
to obtain maximal slicing prescription for the lapse, you should contract the evolution equation for $K_{\mu\nu}$ and obtain an evolution equation for its trace; then, you enforce for all times $t$ the vanishing of $K$, i.e. $K=0$ and $\partial_{t}K=0$, thereby obtaining an elliptic equation for the lapse (ref 1, eq. 4.63, 4.64, 6.90, 9.15); the maximal slicing prescription requires that at each time step, the elliptic equation for the lapse on the hypersurface is solved for before we evolve further
to obtain the harmonic slicing evolution equation, we demand that the time coordinate of our foliation is a harmonic function: $\Box_{g}t=0$. (ref 1, section 9.2.3); in there, the time derivative of lapse appears in a natural way and the derivation answers your question regarding the minus sign, I think; you can generalize this prescription and obtain a family of slicings based on different choices of the $f(\alpha)$ function
The arbitrariness in either of these choices is now evident:
for maximal slicing, you choose that your foliation has $K=0$ everywhere, and you use a no-longer-an-evolution-equation to fix your lapse; inverting the reasoning, you pick such a lapse that will guarantee $K=0$ for each hypersurface
for harmonic slicing, the choice that the time coordinate of the slicing is a harmonic function is totally arbitrary
For details I refer you to the 3+1 Formalism and Bases of Numerical Relativity, specifically Chapter 9 for details on the lapse conditions.
DESCRIPTION OF THE EDIT: I owe you an apology - the answer was not correct. In fact, both parts of the full expression (the one with the Christoffel symbol as well) contribute to the absence of the temporal derivative of the lapse. From one part we get $\sim \partial_{t}\alpha$ and from the other $\sim - \partial_{t}\alpha$, so these cancel out.
In particular, I was wrong stating that $g^{t \mu} = \frac{\beta^{\mu}}{\alpha^{2}}$. This concerns the $g^{ti}$, but not $g^{tt}$. In the quoted part, you can even see that $g^{tt}=\frac{-1}{\alpha^{2}}$.
The fact that $\beta^{t}=0$ is of course true, but does not come into play in this calculation in the end. Also, one must be careful. The orthogonality of a tensor field does not mean that its temporal component vanishes. When the one of the indices of a spatial tensor field component is $t$, the component will be zero by virtue of this contraction: $T_{\;\;\;\mu\nu}^{\alpha\beta}n_{\alpha}=0$ (since $n_{i}=0$). However, for a covariant index, this is not longer the case, easily visible even in the case of $\beta$.
$$ 0 = \beta_{\mu}n^{\mu} = \beta_{t}n^{t} + \beta_{i}n^{i}$$ which leads to: $$\beta_{t}=g_{t\mu}\beta^{\mu}=g_{ti}\beta^{i} = \beta_{i}\beta^{i}$$