force as well as gravitational force because both of them are represented by the inverse-square law: $$\boldsymbol{E}_i = - \sum_{j \neq i}^N q_j \frac{\boldsymbol{r}_j - \boldsymbol{r}_i}{( |\boldsymbol{r}_j - \boldsymbol{r}_i|^2 + \epsilon^2)^{3/2}},$$ where the subscripts $i$, $j$ represent $i/j$th particle; and $N$,$\boldsymbol{r}$, $\epsilon$, $q$, and $\boldsymbol{E}$ are the total particle number, position vector, softening parameter, charge, and electric field, respectively.

In our simulation, we assumed that a star is an aligned rotator and a perfect conductor. To specify the charge sign, we chose aligned ($\boldsymbol{\mu} \cdot
\boldsymbol{\Omega} > 0$), not counter aligned ($\boldsymbol{\mu} \cdot
\boldsymbol{\Omega} < 0$), ie, the poles are negatively charged, while the equator is positively charged. The numerical method follows WS. Magnetospheric plasmas are represented as super-particles. The electric field on the $i$th particle is calculated by the superposition of the vacuum solution $\boldsymbol{E}_{\mathrm v}$ and the space-charge solution $\boldsymbol{E}_{\mathrm q}$, where $$\begin{aligned}
\boldsymbol{E}_{\mathrm v} (\boldsymbol{r}_i) = - \left[ \frac{\Omega \mu R^2 (3\cos^2 \theta -1)}{c r^4_i} - \frac{Q_{\mathrm {sys}}}{r^2_i} \right] \boldsymbol{e}_{\mathrm r} \nonumber
\\ - \left( \frac{2 \Omega \mu R^2 \sin \theta \cos \theta }{c r^4_i}  \right) \boldsymbol{e}_{\mathrm \theta},\end{aligned}$$ $$\begin{aligned}
\boldsymbol{E}_{\mathrm q}  (\boldsymbol{r}_i) = \sum_{j \neq i}^N  \biggl[
\frac{\boldsymbol{r}_i - \boldsymbol{r}_j }{|\boldsymbol{r}_i - \boldsymbol{r}_j|^3}
- \frac{R}{r_j} \frac{\boldsymbol{r}_i - 