I'm trying to make a numerical simulation of pulsar magnetosphere using FDTD on a log-spherical Yee lattice for fields and PIC for plasma particles. Field part is working like charm but issue arises when adding the plasma.
To simulate the dynamics of charged plasma particles I use Newton-Lorenz equation:
$\frac{d\gamma\vec{v}}{dt} = \frac{e}{m}\left( \vec{E} + \frac{\vec{v}}{c} \times \vec{B} \right)$
Adopting CGS system of units, characteristic values for $\vec{E}$ and $\vec{B}$ are of order of magnitule $10^{12}$, charge to mass ratio for electron is of order of magnitude $10^{17}$, so the value for acceleration is ~$10^{29}$.
In the code itself I normalized all variables in the following manner:
$ \vec{u} = \frac{\gamma\vec{v}}{c}\\ \vec{E}^{*} = \vec{E}/B_0\\ \vec{B}^{*} = \vec{B}/B_0\\ \tau = \frac{ct}{R}\\ $
Where $B_0=10^{12}$, $R=10^6$ (radius of the pulsar) and $c$ is the speed of light. As a result, my equation looks as follows:
$\frac{d\vec{u}}{dt} = \frac{eRB_0}{mc^2}\left(\vec{E}^{*}+\frac{\vec{u}\times\vec{B}^{*}}{\gamma}\right)$
Normalization factor $\frac{eRB_0}{mc^2}\sim 10^{15}$, acceleration is still huge.
I may choose different normalization for time as follows:
$t = t_0\tau\\ \frac{eB_0}{mc}t_0 = 1$
But with $\tau$ chosen in that way it takes literally forever for fields to get to stationary configuration.
What can I do?