creation and loss of the particles in the simulation balance with each other. Thus, the total charge of the system is determined automatically.

As for normalization, the position vector, magnetic field, electric field, and velocity are normalized as follows: $\bar{ \boldsymbol{r} } = \boldsymbol{r} / R$, $\bar{ \boldsymbol{B} } = \boldsymbol{B} / B_*$, $\bar{ \boldsymbol{E} } = \boldsymbol{E} / B_*$, and $\boldsymbol{\beta} = \boldsymbol{v} / c$, where $B_*$ is the magnetic field strength on the poles; the variables with bar are non-dimensional in this paper. In our simulation, the strength of the magnetic field is normalized by $B_*$, so that our simulation can treat pulsars which have various magnetic field strengths basically. The other major simulation parameters are $\bar{ \Omega } = R/R_{\mathrm L} =
0.2$, $\bar{m} = 10^{-10}$, and $\bar{q} = 10^{-5}$, where $R_{\mathrm L}$ is the light cylinder radius; $\bar{m}$ and $\bar{q}$ are the mass and charge of a super-particle, respectively, in non-dimensional form.

The simulation parameters should be chosen so that what takes place in the simulation reproduces the actual phenomena in the pulsar magnetosphere. At the same time, we have a restriction of the computational power. Our parameters include the mass and charge of the super-particles 