1)
On first glance it seems to be correct. If your hypersurface is given by $r=R$, then it is the level set of the coordinate function $r$, thus $\mathrm{d}r\sim n$.
2)
The covariant derivative is calculated with respect to the original metric. The connection induced on $\Sigma$ by the induced metric is technically a connection that only works on fields that are intrinsically defined on $\Sigma$. The normal vector is not such a a field.
3)
First of all, let us note that in component-index notation, you can define the extrinsic curvature with spacetime indices, or with 3-indices.
The (2) definition seems awfully confused, with this regard. It is clear that $a$ has to be a 3-index, because $a$ numbers the coordinate-frame vector fields on $\Sigma$, and you have 3 of those, however the $b$ index labels $\nabla$, so it must be a 4-index.
The formula is either wrong, or by what they mean $\nabla_b$ is actually $e^\mu_{\ b}\nabla_\mu$, the restriction of the spacetime connection to directions that are purely tangential to $\Sigma$.
If the latter is the case, then the formula is (aside from a sign) consistent with your (1) definition, since writing the inner product out in index notation nets $$ e_a\cdot\nabla_bn=e^\mu_{\ a}\nabla_bn_\mu=e^\mu_{\ a}e^\nu_{\ b}\nabla_\nu n_\mu, $$ and since $e^\mu_{\ a}=\partial x^\mu/\partial y^a$, this is exactly your formula.
Aside from a sign that is, but some authors define the extrinsic curvature with a different sign.
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 think you're pretty much almost there. I'd just expand out the symmetrization of $\nabla_{(a}n_{b)}$ and collect terms into your $P$'s. One trick to realize is that, since you have a foliation of 3-surfaces that can be labeled by some function $\tau$, you have:
$$n_{a} = \alpha \nabla_{a}\tau$$
where $\alpha = \frac{1}{\sqrt{|\nabla_{a}\tau\nabla^{a}\tau|}}$, which will let you work out identities involving $P^{b}_{c}P^{a}_{d}\nabla_{[a}n_{b]}$ and $P^{b}_{c}P^{a}_{d}\nabla_{(a}n_{b)}$