Deducing GJ electron density

Prove that plasma density surrounding Neutron Star follows:
\rho = \frac{\nabla\cdot \vec{E}}{4 \pi c}= \frac{-\vec{\Omega} \cdot \vec{B}}{2 \pi c }\frac{1}{1-(\Omega r/c)^2 \sin^2 \theta}
because the Lorentz force vanish here, we have:
\vec{E} = - \frac{\vec{\Omega}\times \vec{r} \times \vec{B}}{c}
So,
\frac{1}{4\pi} \nabla \cdot \vec{E} = - \frac{1}{4\pi} \nabla \cdot (\frac{\vec{\Omega}\times \vec{r} \times \vec{B}}{c}) = - \frac{1}{4\pi} \nabla \cdot [\frac{(\vec{\Omega}\times \vec{r}) \times \vec{B}}{c}]
make use of vector formulas, we obtain
\frac{1}{4\pi} \nabla \cdot \vec{E} = - \frac{1}{4\pi} \left\{ \vec{B} \cdot [ \nabla \times \frac{ \vec{\Omega} \times \vec{r}} {c}] - \frac{(\vec{\Omega}\times \vec{r})}{c} \cdot (\nabla \times \vec{B}) \right\}
due to the Maxwell’s equation \nabla \times \vec{B} = 4 \pi \vec{J} = 4 \pi q \vec{\beta_t} = \frac{\vec{\Omega} \times \vec{r}}{c} \cdot (\nabla \cdot \vec{E})
Thus,
\frac{1}{4\pi} \nabla \cdot \vec{E} = - \frac{1}{4\pi} \left\{ \vec{B} \cdot [ \nabla \times \frac{ \vec{\Omega} \times \vec{r}} {c}] - (\frac{\Omega r \sin \theta}{c}\vec{\phi}) (\frac{\Omega r \sin \theta}{c}(\nabla \cdot \vec{E})\vec{\phi}) \right\}

so,
[1-(\Omega r/c)^2 \sin^2 \theta] \frac{\nabla \cdot \vec{E} }{4 \pi} = - \frac{1}{4\pi} [\vec{B}\cdot(\nabla \times \frac{\Omega r \sin \theta}{c} \vec{\phi})]
where, \nabla \times \frac{\Omega r \sin \theta}{c} \vec{\phi} can be calculated by \nabla calculator as spherical coordinate form:

\nabla \times \frac{\Omega r \sin \theta}{c} \vec{\phi} = \frac{1}{r \sin \theta}[\frac{ \partial}{\partial \theta}(\sin \theta \frac{\Omega r \sin \theta}{c})- \frac{\partial 0 }{\partial \phi}] \vec{e_r} + \frac{1}{r}[\sin \theta \frac{\partial 0 }{\partial \phi} - \frac{\partial}{\partial r} ( r \frac{\Omega r \sin \theta}{c})]\vec{e_{\theta}} + \frac{1}{r}[\frac{\partial}{\partial r}(r 0) - \frac{\partial 0 }{\partial \theta}] \vec{e_{\phi}}
\nabla \times \frac{\Omega r \sin \theta}{c} \vec{\phi} = \frac{2\Omega \cos \theta}{c} \vec{e_r} - \frac{2 \Omega r \sin \theta}{c} \vec{e_{\theta}} 
\nabla \times \frac{\Omega r \sin \theta}{c} \vec{\phi} = \frac{2}{c} \vec{\Omega}
Therefore,
[1-(\Omega r/c)^2 \sin^2 \theta] \rho = - \frac{1}{4 \pi}(\vec{B} \cdot \frac{2}{c}\vec{\Omega})

finally, we obtain
\rho = \frac{-\vec{\Omega} \cdot \vec{B}}{2 \pi c }\frac{1}{1-(\Omega r/c)^2 \sin^2 \theta}