3.1. Frame where $\boldsymbol {E}$ and $\boldsymbol {B}$ are parallel

A useful way to follow particle trajectories in any prescribed electromagnetic field given in an inertial frame $K$ consists of Lorentz transforming the electromagnetic field into a special frame $K'$ in which the electric field is parallel to the magnetic field. In the general case, there is an infinite number of frames for which the electric field is parallel to the magnetic field. We can, however, choose the particular velocity given by $\boldsymbol {V}_\parallel = c\boldsymbol {\beta }_\parallel = \alpha \boldsymbol {E} \wedge \boldsymbol {B}$ where $\alpha$ is a constant to be determined. In the frame $K'$, the motion is along the common direction of $\boldsymbol {E}'$ and $\boldsymbol {B}'$ that is the electromagnetic field as measured in this frame $K'$. The constant $\alpha$ is the solution given by

(3.1)\begin{equation} \alpha = \frac{E^2+c^2B^2 - \sqrt{{\mathcal{I}}_1^2 + 4c^2{\mathcal{I}}_2^2}}{2 (\boldsymbol{E} \wedge \boldsymbol{B})^2} . \end{equation}

The minus sign in front of the square root enforces a speed less than that of light. The electric and magnetic fields in the frame moving at speed $\boldsymbol {V}_\parallel$ are found by a special relativistic Lorentz boost of the electromagnetic field that gives

(3.2a)\begin{gather} \boldsymbol{E}' = \varGamma_\parallel [ ( 1 - \alpha B^2 ) \boldsymbol{E} + \alpha ( \boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{B}) \boldsymbol{B} ], \end{gather}

(3.2b)\begin{gather}\boldsymbol{B}' = \varGamma_\parallel [ ( 1 - \alpha E^2/c^2 ) \boldsymbol{B} + \alpha ( \boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{B}) \boldsymbol{E}/c^2 ] . \end{gather}

The Lorentz factor of the frame in which $\boldsymbol {E}'$ and $\boldsymbol {B}'$ are parallel is defined by $\varGamma _\parallel = (1-\beta _\parallel ^2)^{-1/2}$. In this frame, the particle trajectory is decomposed into a motion along the common direction of $\boldsymbol {E}'$ and $\boldsymbol {B}'$ and a gyration around the magnetic field $\boldsymbol {B}'$. Thus, the local tangent to the trajectory becomes $\boldsymbol {t}'_\parallel = \pm \boldsymbol {E}'/E' = \pm \boldsymbol {B}'/B'$, the sign being chosen such that particles flow outwards. Our expression for the particle velocity resembles the Aristotelian expression given by Gruzinov (Reference Gruzinov2013). Our velocity prescription is, however, more general because we do not assume that particles travel exactly at the speed of light. The speed along the common $\boldsymbol {E}$ and $\boldsymbol {B}$ direction is constrained by the electric field acceleration along $\boldsymbol {B}$, contrary to Aristotelian electrodynamics. Gruzinov (Reference Gruzinov2013) introduced two new quantities $E_0>0$ and $B_0$ according to the following invariants (but see also Mestel (Reference Mestel1999)):

(3.3a)\begin{gather} {\mathcal{I}}_1 = E^2 - c^2 B^2 = E_0^2 - c^2 B_0^2, \end{gather}

(3.3b)\begin{gather}{\mathcal{I}}_2 = \boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{B} = E_0 B_0 . \end{gather}

Solving for the magnetic field strength $B_0$ and keeping only the real solution with a positive sign $+$ we get

(3.4)\begin{equation} B_0^2 = \frac{-{\mathcal{I}}_1 + \sqrt{{\mathcal{I}}_1^2+4c^2{\mathcal{I}}_2^2}}{2c^2}. \end{equation}

In such a way, plugging this expression into the electromagnetic field transform given in (3.2) we find

(3.5a)\begin{gather} \boldsymbol{E}' = \frac{\varGamma E_0}{E_0^2/c^2 + B^2} \left[ \frac{E_0}{c^2} \boldsymbol{E} + B_0 \boldsymbol{B} \right], \end{gather}

(3.5b)\begin{gather}\boldsymbol{B}' = \frac{\varGamma B_0}{E_0^2/c^2 + B^2} \left[ B_0 \boldsymbol{B} + \frac{E_0}{c^2} \boldsymbol{E} \right]. \end{gather}

They are therefore colinear because $E_0 \boldsymbol {B}' = B_0 \boldsymbol {E}'$. The frame velocity consequently simplifies into

(3.6)\begin{equation} \boldsymbol V = \frac{\boldsymbol{E} \wedge \boldsymbol{B}}{E_0^2/c^2 + B^2} . \end{equation}

Note in this expression the mixing between field strengths in both frames, one with subscript $0$ and the other without any subscript. There exists, however, a symmetry in the sense that $E_0^2 + c^2 B^2 = E^2 + c^2 B_0^2$ so we can use either $E_0$ or $B_0$ but not both simultaneously. This velocity is always less than the speed of light if ${\mathcal {I}}_1\neq 0$ and ${\mathcal {I}}_2\neq 0$. Therefore, there always exist a frame where $\boldsymbol {E}$ and $\boldsymbol {B}$ are parallel, whatever the strength of $\boldsymbol {E}$ compared to $\boldsymbol {B}$. The vanishing magnetic fields $\boldsymbol {B}=0$ or electric fields $\boldsymbol {E}=0$ are special cases of the general treatment presented here, the second one reducing to the electric drift motion.

The special case of a null field for which ${\mathcal {I}}_1 = {\mathcal {I}}_2 = 0$ is treated separately because then $E=cB$ and the speed of the frame is exactly equal to $c$ and thus is not a physical frame.

3.2. Motion in the frame where $\boldsymbol {E}$ and $\boldsymbol {B}$ are parallel

In the previous section, we showed that it is always possible to reduce the problem of particle motion into a configuration where $\boldsymbol {E}$ and $\boldsymbol {B}$ are parallel except for light-like fields. In this frame, integration of the trajectory is particularly simple when expressed in terms of the proper time $\tau$ of the particle. Any electromagnetic field configuration can always be reduced to parallel electric and magnetic fields by an appropriate change of reference. In order not to overload the notations, in this and the following subsections we omit the primes to designate quantities expressed in the frame $K'$ where $\boldsymbol {E}$ and $\boldsymbol {B}$ are parallel. The integration of the equation of motion of a charged particle in such a field is relatively simple and straightforward. Indeed, let us consider an electromagnetic field such that $\boldsymbol {E}$ and $\boldsymbol {B}$ are directed along the $\boldsymbol {e}_z$ axis in a Cartesian coordinate system. The initial position of the particle is $(x_0, y_0, z_0)$ and its initial velocity is $\boldsymbol {v} = (v_x^0, v_y^0, v_z^0)$. The equation of motion in covariant form is

(3.7)\begin{equation} \frac{\textrm{d} p^i}{\textrm{d}\tau} = q F^{ik} u_k \end{equation}

or in terms of the 4-velocity only

(3.8)\begin{equation} \frac{\textrm{d} u^i}{\textrm{d}\tau} = \frac{q}{m} F^{ik} u_k, \end{equation}

where $u^i=(\gamma c,\gamma \boldsymbol {v})$ is the 4-velocity and $p^i=mu^i$ the 4-momentum. In the Cartesian coordinate system, the electromagnetic field tensor is anti-diagonal and given by

(3.9)\begin{equation} F^{ik} =\begin{pmatrix} 0 & 0 & 0 & -E/c \\ 0 & 0 & -B & 0 \\ 0 & B & 0 & 0 \\ E/c & 0 & 0 & 0 \end{pmatrix}. \end{equation}

Introducing $\omega _E = {qE}/{m c}$ and $\omega _B = {qB}/{m}$, the equation of motion reduces to

(3.10a)\begin{gather} \frac{\textrm{d} u^0}{\textrm{d}\tau} = \frac{q}{m} F^{03} u_3 = - \omega_E u_3, \end{gather}

(3.10b)\begin{gather}\frac{\textrm{d} u^1}{\textrm{d}\tau} = \frac{q}{m} F^{12} u_2 = - \omega_B u_2, \end{gather}

(3.10c)\begin{gather}\frac{\textrm{d} u^2}{\textrm{d}\tau} = \frac{q}{m} F^{21} u_1 = + \omega_B u_1, \end{gather}

(3.10d)\begin{gather}\frac{\textrm{d} u^3}{\textrm{d}\tau} = \frac{q}{m} F^{30} u_0 = + \omega_E u_0 . \end{gather}

All the components of the 4-velocity are brought back to their contravariant expressions so that by index elevation $u^i = \eta ^{ik} u_k$. This implies that $u^0 = u_0$ for the temporal index and $u^a = - u_a$ for the spatial indices. The system decouples into two size 2 subsystems each so that

(3.11a)\begin{gather} \frac{\textrm{d} u^0}{\textrm{d}\tau} = \omega_E u^3, \end{gather}

(3.11b)\begin{gather}\frac{\textrm{d} u^1}{\textrm{d}\tau} = \omega_B u^2, \end{gather}

(3.11c)\begin{gather}\frac{\textrm{d} u^2}{\textrm{d}\tau} = - \omega_B u^1, \end{gather}

(3.11d)\begin{gather}\frac{\textrm{d} u^3}{\textrm{d}\tau} = \omega_E u^0 . \end{gather}

Two variables are eliminated to reduce the system to the velocity components $u^0$ and $u^1$ only

(3.12a)\begin{gather} \frac{\textrm{d}^2 u^0}{\textrm{d}\tau^2} = \omega_E^2 u^0, \end{gather}

(3.12b)\begin{gather}\frac{\textrm{d}^2 u^1}{\textrm{d}\tau^2} = - \omega_B^2 u^1. \end{gather}

The general solutions are given by

(3.13a)\begin{gather} u^0 = A \, \textrm{e}^{\omega_E\tau} + B \, \textrm{e}^{-\omega_E\tau}, \end{gather}

(3.13b)\begin{gather}u^1 = C \cos(\omega_B\tau) + D \sin(\omega_B\tau) . \end{gather}

At the initial time, we have $t = t_0$ corresponding to $\tau =0$ and $\boldsymbol {v} = v_0 \boldsymbol {t}$ thus $\boldsymbol {u}^0 = \gamma _0 (c , \boldsymbol {v}_0)$ with $\gamma _0=(1-(\beta _0)^2)^{-1/2}$. These initial conditions enforce

(3.14a)\begin{gather} A = \gamma_0 \frac{c+v_0^z}{2}, \end{gather}

(3.14b)\begin{gather}B = \gamma_0 \frac{c-v_0^z}{2}, \end{gather}

(3.14c)\begin{gather}C = \gamma_0 v_0^x, \end{gather}

(3.14d)\begin{gather}D = \gamma_0 v_0^y, \end{gather}

thus the 4-velocity evolution is given in terms of the proper time according to

(3.15a)\begin{gather} u^0 = \gamma_0 c \left[ \textrm{ch} (\omega_E\tau) + \beta_0^z \,\textrm{sh} (\omega_E\tau) \right], \end{gather}

(3.15b)\begin{gather}u^3 = \gamma_0 c \left[ \textrm{sh} (\omega_E\tau) + \beta_0^z \,\textrm{ch} (\omega_E\tau) \right], \end{gather}

(3.15c)\begin{gather}u^1 = \gamma_0 c \left[ \beta_0^x \cos (\omega_B\tau) + \beta_0^y \sin (\omega_B\tau) \right], \end{gather}

(3.15d)\begin{gather}u^2 = \gamma_0 c \left[ -\beta_0^x \sin (\omega_B\tau) + \beta_0^y \cos (\omega_B\tau) \right] . \end{gather}

All what remains is to integrate with respect to the proper time to find the trajectory of the particle

(3.16a)\begin{gather} c ( t-t_0) = \frac{\gamma_0 c }{\omega_E} \left[ \textrm{sh} (\omega_E \tau) + \beta_0^z ( \textrm{ch} (\omega_E \tau) - 1 ) \right], \end{gather}

(3.16b)\begin{gather}x - x_0 = \frac{\gamma_0 c }{\omega_B} \left[ \beta_0^x \sin (\omega_B\tau) - \beta_0^y ( \cos (\omega_B\tau) - 1 ) \right], \end{gather}

(3.16c)\begin{gather}y - y_0 = \frac{\gamma_0 c }{\omega_B} \left[ \beta_0^x ( \cos (\omega_B\tau) - 1 ) + \beta_0^y \sin (\omega_B \tau) \right], \end{gather}

(3.16d)\begin{gather}z - z_0 = \frac{\gamma_0 c }{\omega_E} \left[ \textrm{ch} (\omega_E \tau) - 1 + \beta_0^z \,\textrm{sh} (\omega_E \tau) \right] . \end{gather}

The trajectory is thus entirely determined as a function of proper time $\tau$ in an analytical way with simple expressions including trigonometric and hyperbolic functions. These equations for 4-velocity and 4-position are implemented in the code.

Note that the equation of motion could be resolved immediately in the observer's frame of reference by diagonalizing the tensor of the electromagnetic field $F^{ik}$. This would be the same as the change of reference frame made above (Vandervoort Reference Vandervoort1960).

For testing the numerical algorithm and the limiting Lorentz factor due to numerical round off error, we check the implementation on a purely electric and a purely magnetic field.

3.3. Vanishing magnetic field

In the case of a vanishing magnetic field $\boldsymbol {B}=0$, the frames $K$ and $K'$ are identical and there is no need to boost from one frame to the other. The 4-velocity reduces to

(3.17a)\begin{gather} u^0 = \gamma_0 c \left[ \textrm{ch} (\omega_E \tau) + \beta_0^z \,\textrm{sh} (\omega_E \tau) \right], \end{gather}

(3.17b)\begin{gather}u^1 = \gamma_0 c \beta_0^x, \end{gather}

(3.17c)\begin{gather}u^2 = \gamma_0 c \beta_0^y, \end{gather}

(3.17d)\begin{gather}u^3 \gamma_0 c \left[ \textrm{sh} (\omega_E \tau) + \beta_0^z \,\textrm{ch} (\omega_E \tau) \right] \end{gather}

and the trajectory simply into

(3.18a)\begin{gather} c ( t-t_0) = \frac{\gamma_0 c }{\omega_E} \left[ \textrm{sh} (\omega_E \tau) + \beta_0^z ( \textrm{ch} (\omega_E \tau) - 1 ) \right], \end{gather}

(3.18b)\begin{gather}x - x_0 = \gamma_0 c \beta_0^x \tau, \end{gather}

(3.18c)\begin{gather}y - y_0 = \gamma_0 c \beta_0^y \tau, \end{gather}

(3.18d)\begin{gather}z - z_0 = \frac{\gamma_0 c }{\omega_E} \left[ \textrm{ch} (\omega_E \tau) - 1 + \beta_0^z \,\textrm{sh} (\omega_E \tau) \right] . \end{gather}

3.4. Vanishing electric field

In the case of a vanishing electric field $\boldsymbol {E}=0$, the frames $K$ and $K'$ are again identical and there is no need to boost from one frame to the other. The 4-velocity reduces to

(3.19a)\begin{gather} u^0 = \gamma_0 c , \end{gather}

(3.19b)\begin{gather}u^1 = \gamma_0 c \left[ \beta_0^x \cos (\omega_B \tau) + \beta_0^y \sin (\omega_B \tau) \right], \end{gather}

(3.19c)\begin{gather}u^2 = \gamma_0 c \left[ -\beta_0^x \sin (\omega_B \tau) + \beta_0^y \cos (\omega_B \tau) \right], \end{gather}

(3.19d)\begin{gather}u^3 = \gamma_0 c \beta_0^z. \end{gather}

All that remains is to integrate with respect to the proper time to find the trajectory of the particle

(3.20a)\begin{gather} (t-t_0) = \gamma_0 \tau, \end{gather}

(3.20b)\begin{gather}x - x_0 = \frac{\gamma_0 c }{\omega_B} \left[ \beta_0^x \sin (\omega_B \tau) - \beta_0^y ( \cos (\omega_B \tau) - 1 ) \right], \end{gather}

(3.20c)\begin{gather}y - y_0 = \frac{\gamma_0 c }{\omega_B} \left[ \beta_0^x ( \cos (\omega_B \tau) - 1 ) + \beta_0^y \sin (\omega_B \tau) \right] , \end{gather}

(3.20d)\begin{gather}z - z_0 = \gamma_0 c \beta_0^z \tau = v_0^z ( t-t_0) . \end{gather}

Testing our algorithm in a constant and uniform electromagnetic field is meaningless because the solution is known analytically in any reference frame. Even the force-free field test introduced by Ripperda et al. (Reference Ripperda, Bacchini, Teunissen, Xia, Porth, Sironi, Lapenta and Keppens2018) is included in this exact analytical solution, as demonstrated in the next section. However, for numerical purposes, it is desirable to check the ability of our code to handle very high Lorentz factors above $\gamma =10^{12}$ to look for round-off errors and possible issues related to finite digit precision.

3.5. Motion in the frame where $\boldsymbol {E}$ and $\boldsymbol {B}$ are perpendicular

When electric and magnetic fields are perpendicular, the second invariant vanishes ${\mathcal {I}}_2=0$, there exist always a frame where either the electric or the magnetic field vanishes depending on the sign of the invariant ${\mathcal {I}}_1$. We are then back to the previous cases for a pure electric or magnetic field. Indeed, if ${\mathcal {I}}_1>0$ the constant $\alpha =c^2/E^2$ and therefore the velocity of the frame where $\boldsymbol {B}$ vanishes is

(3.21)\begin{equation} \boldsymbol{V}_B = \frac{\boldsymbol{E} \wedge \boldsymbol{B}}{E^2} c^2 . \end{equation}

If ${\mathcal {I}}_1<0$ the constant $\alpha =1/B^2$ and therefore the velocity of the frame where $\boldsymbol {E}$ vanishes is

(3.22)\begin{equation} \boldsymbol{V}_E = \frac{\boldsymbol{E} \wedge \boldsymbol{B}}{B^2}, \end{equation}

that is the usual electric drift frame. Consequently, the cases ${\mathcal {I}}_1\neq 0$ and ${\mathcal {I}}_2=0$ are included in the previous sections.

3.6. Light-like electromagnetic field

Nevertheless, the case where both electromagnetic invariants vanish ${\mathcal {I}}_1 = {\mathcal {I}}_2 = 0$ must be treated separately because there exists no physical frame moving at a speed strictly less than $c$ where either $\boldsymbol {E}$ or $\boldsymbol {B}$ vanishes. Going back to the equation of motion, let us assume that the electric field is along $\boldsymbol {E} = E \boldsymbol {e}_{y}$ and the magnetic field along $\boldsymbol {B} = B \boldsymbol {e}_{z}$. The electromagnetic tensor then reduces to

(3.23)\begin{equation} F^{ik} =\begin{pmatrix} 0 & 0 & -E/c & 0 \\ 0 & 0 & -B & 0 \\ E/c & B & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} . \end{equation}

The equation of motion therefore becomes

(3.24a)\begin{gather} \frac{\textrm{d} u^0}{\textrm{d}\tau} = \omega_E u^2, \end{gather}

(3.24b)\begin{gather}\frac{\textrm{d} u^1}{\textrm{d}\tau} = \omega_B u^2, \end{gather}

(3.24c)\begin{gather}\frac{\textrm{d} u^2}{\textrm{d}\tau} = \omega_E u^0 - \omega_B u^1, \end{gather}

(3.24d)\begin{gather}\frac{\textrm{d} u^3}{\textrm{d}\tau} = 0 . \end{gather}

The last equation integrates into $u^3=u^3_0 = constant$ and the trajectory is constrained to follow $z-z_0 = u^3_0(\tau -\tau _0)$. Eliminating $u^0$ and $u^1$ from the equation evolving $u^2$ we find

(3.25)\begin{equation} \frac{\textrm{d}^2u^2}{\textrm{d}\tau^2} + (\omega_B^2 - \omega_E^2) u^2 = 0 . \end{equation}

Three cases must be distinguished

1. (i) a dominant electric field for $\omega _E > \omega _B$.

2. (ii) a dominant magnetic field for $\omega _B > \omega _E$.

3. (iii) a light like field if $\omega _E = \omega _B$.

Note that $\omega _E^2 - \omega _B^2$ is related to the relativistic invariant ${\mathcal {I}}_1$ through $({q^2}/{m^2 c^2}) {\mathcal {I}}_1 = \omega _E^2 - \omega _B^2$.

In a light-like electromagnetic field, there is no frame where $\boldsymbol {E}$ or $\boldsymbol {B}$ vanishes. A plane electromagnetic wave propagating in a vacuum is a typical example. In this field, the two Lorentz invariants vanish ${\mathcal {I}}_1 = {\mathcal {I}}_2=0$. The equation for the 4-velocity reduces to

(3.26)\begin{equation} \frac{\textrm{d}^2 u^2}{\textrm{d}\tau^2} = 0, \end{equation}

whose solution is $u^2 = C \tau + \gamma _0v_0^y$, where $C$ is a constant of integration. The other two components are

(3.27a)\begin{gather} u^0 = \omega_E \left( C \frac{\tau^2}{2} + \gamma_0 v_0^y\tau\right) + \gamma_0 c , \end{gather}

(3.27b)\begin{gather}u^1 = \omega_B \left( C \frac{\tau^2}{2} + \gamma_0 v_0^y\tau\right) + \gamma_0 v_0^x. \end{gather}

Substituting in $\textrm {d}u^2/\textrm {d}\tau$ we find $C = \gamma _0 c (\omega _{E} - \omega _{B} \beta _0^x)$ because $\omega _E = \pm \omega _B$ and cannot be simplified further. The solution is therefore

(3.28a)\begin{gather} u^0 = \gamma_0 c \left[ 1 + \omega_{E} (\omega_{E} - \omega_{B} \beta_0^x) \frac{\tau^2}{2} + \beta_0^y \omega_{E} \tau \right], \end{gather}

(3.28b)\begin{gather}u^1 = \gamma_0 c \left[ \beta_0^x + \omega_{B} (\omega_{E} - \omega_{B} \beta_0^x) \frac{\tau^2}{2} + \beta_0^y \omega_{B} \tau \right], \end{gather}

(3.28c)\begin{gather}u^2 = \gamma_0 c [ \beta_0^y + (\omega_{E} - \omega_{B} \beta_0^x) \tau ], \end{gather}

(3.28d)\begin{gather}u^3 = \gamma_0 c \beta_0^z . \end{gather}

This solution allows us to advance the particle velocity in time without taking care of the respective signs of $\omega _{E}$ and $\omega _{B}$. Performing another integration gives the explicit form of the trajectory as

(3.29a)\begin{gather} c (t-t_0) = \gamma_0 c \left[ \tau + \omega_{E} (\omega_{E} - \omega_{B} \beta_0^x) \frac{\tau^3}{6} + \beta_0^y \omega_{E} \frac{\tau^2}{2} \right], \end{gather}

(3.29b)\begin{gather}x - x_0 = \gamma_0 c \left[ \beta_0^x \tau + \omega_{E} (\omega_{E} - \omega_{B} \beta_0^x) \frac{\tau^3}{6} + \beta_0^y \omega_B \frac{\tau^2}{2} \right], \end{gather}

(3.29c)\begin{gather}y - y_0 = \gamma_0 c \left[ \beta_0^y \tau + (\omega_{E} - \omega_{B} \beta_0^x)\frac{\tau^2}{2} \right], \end{gather}

(3.29d)\begin{gather}z - z_0 = \gamma_0 c \beta_0^z \tau . \end{gather}

These solutions must be included in the algorithm whenever both invariants vanish. Tests will be performed in a linearly and a circularly polarized plane electromagnetic wave in § 5.

3.7. Force-free field

In the special case of a force-free field, the electric and magnetic fields are perpendicular according to the expression

(3.30)\begin{equation} \boldsymbol{E} + \boldsymbol{v} \wedge \boldsymbol{B} = \boldsymbol{0}, \end{equation}

where $\boldsymbol {v}$ is the particle velocity in the frame $K$. In that case, by construction $\boldsymbol {E} \boldsymbol {\cdot } \boldsymbol {B} = 0$ and $E < cB$. The privileged frame $K'$ according to the general transformation rule is therefore simply the electric drift frame moving at a velocity

(3.31)\begin{equation} \boldsymbol{V}_{E} = \frac{\boldsymbol{E} \wedge \boldsymbol{B}}{B^2} = \boldsymbol{v} = \boldsymbol{V}_\parallel . \end{equation}

Consequently, the particle velocity $\boldsymbol {v}$ is equal to the drift velocity $\boldsymbol {V}_{E}$ and equal to $\boldsymbol {V}_\parallel$, meaning that in the frame where $\boldsymbol {E}$ and $\boldsymbol {B}$ are parallel (actually $\boldsymbol {E}=0$) the particle is at rest. Thus no Lorentz force acts at all in $K'$ because initially $\boldsymbol {v}'=0$ and $\boldsymbol {E}'=0$, the particle remains at rest forever in $K'$ as long as the force-free condition (3.30) is satisfied, whatever the value of $\boldsymbol {B}'$. This result is analytically exact, there is no numerical error, even to round-off accuracy, contrasting severely with other algorithms showing a slowly breaking of the force-free motion (Ripperda et al. Reference Ripperda, Bacchini, Teunissen, Xia, Porth, Sironi, Lapenta and Keppens2018).

3.8. Low electric and magnetic fields

The expression for particle trajectories given in (3.16) cannot be applied straightforwardly whenever $\omega _{E}$ or $\omega _{B}$ tends to zero because they appear in the denominator of some expressions. However, to include these cases in the general algorithm we designed, we performed a first-order expansion whenever required to smoothly join the vanishing electric or magnetic field case presented previously.

In regions where the electric field is weak according to the condition $\omega _E \tau \ll 1$ the time and $z$ components are integrated according to

(3.32a)\begin{gather} (t-t_0) = \gamma_0 \tau , \end{gather}

(3.32b)\begin{gather}z - z_0 = v_0^z ( t-t_0) . \end{gather}

In regions where the magnetic field is weak according to the condition $\omega _B \tau \ll 1$ the $x$ and $y$ components are integrated according to

(3.32c)\begin{gather} x - x_0 = \gamma_0 c \beta_0^x \tau, \end{gather}

(3.32d)\begin{gather}y - y_0 = \gamma_0 c \beta_0^y \tau . \end{gather}

Taking into account these two limiting cases of vanishing electric or magnetic field, our algorithm can handle any field strength and any geometry of the electromagnetic field as long as these fields are constant in time and uniform in space. Thus its robustness must be tested against spatially and temporally varying fields. This is checked in the next section. Nevertheless, our focus in this paper is about ultra-relativistic regime of particle motion, dealing with a Lorentz factor as high as $\gamma = 10^{12}$. Although the analytical expressions are obviously valid for any $\gamma$, the finite precision of numerical implementation of any algorithm limits the possible range of Lorentz factors achievable. Therefore, we will also check our code on simple test cases such as a purely electric or purely magnetic field in order to assess the stringent limitations of what can be done numerically with ultra-relativistic plasmas.

But before running simulations, we have to properly choose the optimal time step to evolve the particle motion. This crucial problem is addressed in the following section.

4.1. Time step setting in different frames

In a first step, the relation between particle proper time $\Delta \tau$ and the time $\Delta t'$ in frame $K'$, (3.16) can be inverted analytically to give

(4.1)\begin{equation} \omega_E \Delta \tau = \log \left( \frac{\xi + \sqrt{1-(\beta_0^z)^2 + \xi^2 }}{1+\beta_0^z} \right) = - \log \left( \frac{-\xi + \sqrt{1-(\beta_0^z)^2 + \xi^2 }}{1-\beta_0^z} \right), \end{equation}

where we introduced $\xi =\omega _E \Delta t'/\gamma _0 + \beta _0^z$. This expression is optimal for numerical computation, working also when $\beta _0^z$ is nearly +1 for the first expression and nearly $-$1 for the second expression. In a second step, the time $\Delta t'$ must be related to the observer time $\Delta t$ in frame $K$, the one used in the simulation to advance in time one time step; $\Delta t$ and $\Delta t'$ are related by a Lorentz boost similar to (2.1) written most efficiently as

(4.2)\begin{equation} \Delta t = \varGamma \left( \Delta t' + \frac{\boldsymbol{V}\boldsymbol{\cdot} \Delta \boldsymbol{r}'}{c^2} \right), \end{equation}

where $\varGamma = (1-V^2/c^2)^{-1/2}$. The particle advances in position $\Delta \boldsymbol {r}'$ in frame $K'$ known analytically from (3.16) if $\Delta \tau$ is known. Symbolically, we write $\Delta \boldsymbol {r}'(\Delta \tau )$. Moreover, $\Delta \tau$ is found from $\Delta t'$ thanks to (4.1) or symbolically $\Delta \tau (\Delta t')$ thus $\Delta t$ is found from solving the nonlinear scalar equation (4.2), remembering that finally $\Delta \boldsymbol {r}'(\Delta t')$.

This procedure is best understood by a simple example. Let us consider a particle moving along an arbitrary direction in frame $K$ with constant speed $\boldsymbol {v}$. In the frame $K'$, it moves at a constant speed $\boldsymbol {v}'$, $K'$ moving with respect to $K$ at a speed $\boldsymbol {V}$. How to relate then the speeds $\boldsymbol {v}$ and $\boldsymbol {v}'$ knowing the trajectory in $K'$?

In this simple case, the vector position is given in frame $K'$ by $\Delta \boldsymbol {r}' = \boldsymbol {v}' \Delta t'$. Therefore $\Delta t'$ can be solved with respect to $\Delta t'$ as

(4.3)\begin{equation} \Delta t' = \frac{\Delta t}{\varGamma (1 + \boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{v}' / c^2)}, \end{equation}

$\Delta \boldsymbol {r}$ is then related to $\Delta \boldsymbol {r}'$ via the spatial part of the Lorentz transformation. In such a way, the velocity

(4.4)\begin{equation} \boldsymbol{v} = \frac{\Delta \boldsymbol{r}}{\Delta t} = \frac{\boldsymbol{v}'}{\varGamma (1 + \boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{v}' / c^2)} + \left( 1 + \frac{\varGamma}{\varGamma + 1} ( \boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{v}' / c^2 ) \right) \frac{\boldsymbol{V}}{1 + \boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{v}' / c^2} \end{equation}

is computed from the velocity in the frame $K'$. It can be checked that the results agree with the relativistic composition of velocities. The crucial point in this derivation is the analytical inversion of the relation $\Delta t(\Delta t')$ into $\Delta t'(\Delta t)$. This is no longer the case for a particle in an electromagnetic field. We therefore have to resort to numerical inversion by some root finding methods like the Newton–Raphson scheme or other techniques, as described for instance in Press (Reference Press2007).

4.2. Adaptive time steps

The optimal time step strongly depends on the local value of the electromagnetic field and on the particle velocity. Due to the Lorentz transformation of the electromagnetic tensor and because of time dilation, for ultra-relativistic particles, we expect a significant gain in time computation when this time step is adaptively adjusted to the local fields and particle velocities. In order to improve and better control the accuracy of our numerical solution from time step to time step, we implemented also a semi-implicit iterative scheme where the actual constant electromagnetic field employed to advance the particle is the one located midway between the time and position of the particle at the instant $t^n$ and $t^{n+1}$. It resembles the algorithm we already used in Pétri (Reference Pétri2017). If the solution does not converge to a prescribed precision after a limited number of iterations $N_{\textrm {max}}$ (we set it to $N_{\textrm {max}}=10$), the time step is decreased by a factor 2 and the process starts again with at most $N_{\textrm {max}}$ iterations. If necessary, the time step is again diminished by a factor 2 until convergence is reached. Figure 1 summarizes the pseudo-code used to advance the particle position one time step by using the adaptive scheme.

Figure 1.

Pseudo-code summarizing the Picard iteration scheme to advance one time step; $\boldsymbol {F}$ means that it must apply to both fields $\boldsymbol {E}$ and $\boldsymbol {B}$.

Having an algorithm to impose a fixed or prescribed observer time step, we can couple the particle orbit integrator to the field advance like any other existing PIC code. There is no special care required to first advance particles and second to compute the new fields. Let us now look at some particular tests of our algorithm.

5.1. Normalization and adimensionalization

Before showing some numerical results, we normalize the relevant quantities to characteristic values. As we deal with relativistic motion, the speed is conveniently normalized to the speed of light $c$. Next we introduce a characteristic frequency $\omega$ related the electromagnetic frequencies $\omega _{B}$, $\omega _{E}$ or the plasma frequency $\omega _{p}$ depending on the problem studied. This naturally leads to a characteristic length given by $c/\omega$. Moreover, a normalized time is then introduced by $\tilde {t} = \omega t$. In the simulation results shown below, we will always refer to these normalized quantities and plot graphs according to this convention. Concretely, for simulation purposes, we specify the charge $q$ and mass $m$ of the particle in the different tests such that $q=1$ and $m=1$ if not otherwise specified.

5.2. Purely electric field

Let us start with a homogeneous and uniform electric field. Take as an initial condition $t_0=0$ and the position of the particle to be $x_0 = y_0 = z_0 = 0$ with no initial velocity such that $\boldsymbol {\beta } = \boldsymbol {0}$ and thus $\gamma _0=1$. The particle world line is therefore a straight line parametrized with respect to its proper time $\tau$ according to

(5.1a)\begin{gather} \omega_E t = \textrm{sh} (\omega_E \tau), \end{gather}

(5.1b)\begin{gather}x = 0 , \end{gather}

(5.1c)\begin{gather}y = 0 , \end{gather}

(5.1d)\begin{gather}\omega_E z = c \left[ \textrm{ch} (\omega_E \tau) - 1 \right] . \end{gather}

Expressed in terms of the observer time $t$, the trajectory and velocity become

(5.2a)\begin{gather} \omega_E z = c \left[ \sqrt{ 1 + (\omega_E t)^2} - 1 \right], \end{gather}

(5.2b)\begin{gather}\frac{v_z}{c} = \frac{\omega_E t}{\sqrt{ 1 + (\omega_E t)^2}} . \end{gather}

The Lorentz factor grows with time $t$ according to

(5.3)\begin{equation} \gamma = \textrm{ch} (\omega_E \tau) = \sqrt{ 1 + (\omega_E t)^2} . \end{equation}

Note that the only relevant time scale in this problem is the normalized quantity $\tilde {t} = \omega _E t$ or $\tilde {\tau } = \omega _E \tau$.

An example of an accelerating electric field is shown in figure 2 with the position in blue and for the Lorentz factor in red. What matters is not the particle charge and mass, but the quantity $q\,E/m c = \omega _{E}$ which gives the typical time scale for acceleration. The time step in the observer frame $\Delta t$ is imposed by the user. We took an initial value of $\omega _E \Delta t = 10^{-6}$ although there is no restriction on $\Delta t$. Indeed, after each iteration we multiplied it by 2 in order to show the flexibility of adapting the observer time step. Then, the proper time step $\Delta \tau$ is computed according to (4.1). The particle first accelerates in the Newtonian regime with a quadratic increase in position $z$ up to the point where it reaches almost the speed of light. After a time $\omega _E t\gtrsim 1$ it goes at almost constant speed $v_{\textrm {z}}\approx c$. The Lorentz factor then increases almost linearly with time $\gamma \approx \omega _E t$. We let the particle gain energy up to $\gamma = 10^{20}$ to check possible issues related to numerical round off and truncation. No special problems were met for these ultra-relativistic speeds.

Figure 2.

Position $z$ in blue and Lorentz factor in red of an electron in the accelerating electric field $\boldsymbol {E} = E \boldsymbol {e}_{\textrm {z}}$ and with increasing observer time step $\Delta t$.

5.3. Purely magnetic field

Let us go on with a homogeneous and uniform magnetic field. Take as an initial condition $t_0=0$ and the position of the particle in the $xOy$ plane such that $x_0=z_0=0$ and $y_0=\gamma _0 v_0^x/\omega _B = r_\textrm {L}$ with initial velocity such that $\boldsymbol {\beta } = (\beta _0^x,0,\beta _0^z)$. $r_\textrm {L} = \gamma v/\omega _{B}$ is the Larmor radius of the trajectory. The Lorentz factor is constant and given by $\gamma =\gamma _0$. The particle world line is therefore

(5.4a)\begin{gather} t = \gamma_0 \tau, \end{gather}

(5.4b)\begin{gather}x = r_\textrm{L} \sin (\omega_B \tau), \end{gather}

(5.4c)\begin{gather}y = r_\textrm{L} \cos (\omega_B \tau), \end{gather}

(5.4d)\begin{gather}z = v_0^z t . \end{gather}

The particle gyro-frequency is $\omega _B$ in proper time but reduced to $\omega _B/\gamma _0$ in the observer frame as is well known from special relativity. In the observer frame, the relevant normalized time scale is therefore $\omega _{B} t$. For an ultra-relativistic particle with $v_0^x \approx c$, its Larmor radius is $r_\textrm {L} \approx \gamma _0 c/\omega _B$, the expression used in the introduction. An example is shown in figure 3 for $B=1$ and $\gamma = 10^{10}$. The time step is set to $\omega _B \Delta \tau = 10^{-2}$. The radius of the orbit stays at $r_\textrm {L}$ to very high accuracy, more than 15 digits of precision. The Lorentz factor remains constant and equal to $\gamma _0$ as expected.

Figure 3.

Evolution of the error in the orbital radius of an electron in the uniform magnetic field with $\gamma = 10^{10}$.

5.4. Crossed electric and magnetic fields

As another proof of the efficiency of the algorithm, we compute the trajectories in a crossed electromagnetic field ($\boldsymbol {E} \boldsymbol {\cdot } \boldsymbol {B} = 0$) where the average motion is an electric drift in the $\boldsymbol {E} \wedge \boldsymbol {B}$ direction at the electric drift speed $\boldsymbol {v}_{E} = \boldsymbol {E} \wedge \boldsymbol {B}/B^2$. In the frame moving at $\boldsymbol {v}_{E}$ the electric field vanishes and the particle simply follows a helicoidal motion along the constant magnetic field, see the case treated in the previous paragraph. This motion is only allowed for weak electric fields satisfying $E < cB$ thus enforcing $v_{E} < c$. For concreteness, let us assume an electric field directed along $\boldsymbol {e}_{y}$ and a magnetic field directed along $\boldsymbol {e}_{z}$. The electric drift speed becomes $\boldsymbol {v}_{E} = (E/B) \boldsymbol {e}_{x}$. A Lorentz transformation of the electromagnetic field with Lorentz factor $\varGamma _{E} = 1/\sqrt {1-v_{E}^2/c^2}$ along $\boldsymbol {v}_{E}$ shows that in the comoving frame

(5.5a)\begin{gather} \boldsymbol{E}' = 0 , \end{gather}

(5.5b)\begin{gather}\boldsymbol{B}' = \boldsymbol{B}/\varGamma_{E} . \end{gather}

The electric field vanishes as expected and the magnetic field is decreased potentially by a large ratio equal to the Lorentz factor of the comoving frame. As initial conditions for the particle position and velocity, we choose a helicoidal motion as explained in the previous paragraph and corresponding to an evolution in the solely magnetic field $\boldsymbol {B}'$ as seen in the electric drift frame. These quantities are then transformed to the observer inertial frame according to Lorentz transformations for the velocity $\boldsymbol {v}$ of the particle. The algorithm is checked by computing the particle trajectory in the drift frame and the corresponding Lorentz factor of the particle that should remain constant in that frame. Using the Lorentz transformation the coordinates in the drift frame are

(5.6a)\begin{gather} x' = \varGamma_{E} ( x - v_{E} t ), \end{gather}

(5.6b)\begin{gather}y' = y , \end{gather}

(5.6c)\begin{gather}z' = z . \end{gather}

For numerical purposes, the intensity of the electric field is set such that $\varGamma _{E}=10^{3}$, the particle Lorentz factor in this drift frame is $\gamma =10^{10}$ and that of the magnetic field is $B=1$. Typical results are depicted in figure 4 for the trajectory in the electric drift frame which is usually a helicoidal motion and here exactly a circle in the comoving plane $x'O'y'$. The trajectory projected onto the $x'O'y'$ plane remains a circle to very good accuracy with no change in radius within 8 digits, see figure 5. The accuracy seems less good than in the previous case for a purely magnetic field. We lost several digits in the computation of the Larmor radius. This loss of precision is imputed to the procedure used to evaluate the trajectory in the drift frame. Indeed, a Lorentz transform is required for the particle position as explained in (5.6). This coordinate transform induced an additional error in the Larmor radius estimate. We lose several digits during the subtraction for $x'$. We checked that the error does not decrease with decreasing time step.

Figure 4.

Gyromotion of an electron in the electric drift frame with $\varGamma _{E} = 10^{3}$ and $\gamma =10^{10}$. The Larmor radius is $r_\textrm {L}=10^{10}$.

Figure 5.

Relative error in the Larmor radius computed from the transformation in (5.6). The proper time $\tau$ is shown on the $x$-axis normalized to the cyclotron frequency in the electric drift frame $\omega _{\textrm{B}'}$.

The most general and realistic fields are spatially and temporally varying. In these cases, the electromagnetic field has to be determined at some position and time during the particle motion. As severe tests of our algorithm, we study four problems of which three have exact analytical relativistic solutions for the particle trajectory. The relativistic electrostatic Kepler two body problem of an electron orbiting around a fixed positive ion is an interesting test for a spatially varying electric field. It is a central force case. The other two solutions correspond to a particle moving in a plane electromagnetic wave linearly or circularly polarized, a light-like electromagnetic field. A last less trivial example is depicted by a particle drifting in the equatorial plane of a static magnetic dipole. We discuss in depth these regimes in the following sections. The analytical solutions are described in Uzan & Deruelle (Reference Uzan and Deruelle2014) and Gourgoulhon (Reference Gourgoulhon2010). For completeness we recall them in the following paragraphs.

5.5. Central electric force

The two body problem in gravitational physics can be transposed to an equivalent electrostatic problem including a relativistic velocity. In this case, a particle with charge $q$ and mass $m$ orbits around a fixed central particle with charge $Q$. For bounded orbits we require $qQ<0$. Solutions are given by conservation of energy ${\mathcal {E}}$ and angular momentum $L$. The electric force applied to the orbiting particle is

(5.7)\begin{equation} \boldsymbol{f} = \frac{qQ}{4{\rm \pi}\varepsilon_0 r^3} \boldsymbol{r} . \end{equation}

The orbital motion stays in a plane that we choose as the $xOy$ plane. For $|qQ|<4{\rm \pi} \varepsilon _0 L c$, the general solution is

(5.8a)\begin{gather} r(t) = \frac{p}{1 + e \cos(\varOmega_{p} (\varphi(t) - \omega))}, \end{gather}

(5.8b)\begin{gather}\varOmega_{p} = \sqrt{1 - \left( \frac{qQ}{4{\rm \pi}\varepsilon_0L c}\right)^2 } , \end{gather}

(5.8c)\begin{gather}p = \frac{\varOmega_{p}^2}{- \dfrac{qQ}{4{\rm \pi}\varepsilon_0 c^2} \dfrac{{\mathcal{E}}}{L^2}}, \end{gather}

(5.8d)\begin{gather}e^2 = \frac{1}{{\mathcal{E}}^2} \left[ m^2 c^4 + \frac{{\mathcal{E}}^2 - m^2 c^4}{\left( \dfrac{qQ}{4{\rm \pi}\varepsilon_0L c}\right)^2} \right] . \end{gather}

The particle trajectory is plane and therefore described in a cylindrical coordinate system $(r,\varphi )$ by the parametric function $r(\varphi )$. It depends implicitly on time $t$ because of $\varphi$ being a function of time $t$ (the solution for $\varphi$ is not shown here). Note also that $\varOmega _{p}<1$, meaning that the path is not a closed curved but a prograde precessing ellipse; $p$ is called the orbital parameter and is related to the semi-major axis $a$ and eccentricity of the ellipse via $p=a\,(1-e^2)$. The periastron is located at a phase $\omega$ with respect to the $x$-axis. These notations are common when studying stellar orbits in binary systems. This initial phase $\omega$ is deduced from the initial condition $r=r_0$ at $\varphi =\varphi _0$. Explicitly, we find

(5.9)\begin{equation} \varOmega_{p} (\varphi_0 - \omega) = \arccos \left( \frac{p/r_0 - 1}{e} \right) . \end{equation}

An example of relativistic particle trajectory showing the precession of the orbit is given in figure 6. A piece of the exact analytical solution is also shown and matches perfectly the output of the numerical simulations. The total energy ${\mathcal {E}}$ is split into relativistic kinetic energy $\gamma \,m c^2$ and electrostatic potential energy $U$. Inspection of figure 7 demonstrates that the total energy is accurately conserved during time evolution. Angular momentum is also conserved to high accuracy, at least 6 digits, see the relative error evolving in time in figure 8. Finally, the relative error in the total energy $\Delta {\mathcal {E}}/{\mathcal {E}}$ depending on time step is shown in figure 9. We conclude that the scheme is second order in time, decreasing the error according to $\Delta {\mathcal {E}}/{\mathcal {E}} \propto \Delta t^{2}$.

Figure 6.

Motion of an electron in the electric field of a fixed proton, black points. The relativistic precession of the orbit is clearly visible. The minimum and maximum radii of the orbit as predicted by the analytical formulae are shown by two circles tangent to the trajectory. A piece of the exact analytical solution is also shown in red.

Figure 7.

Total energy ${\mathcal {E}}$, relativistic kinetic energy $\gamma m c^2$ and electrostatic potential energy $U$ of an electron in the electric field of a fixed proton; ${\mathcal {E}}$ is accurately conserved, being a constant of motion.

Figure 8.

Relative error in the angular momentum $L$; $L$ is conserved within 11 digits.

Figure 9.

Relative error in the total energy depending on time step $\Delta t$.

5.6. Magnetic drift in a dipole

The magnetic field created by a dipole of magnetic dipole moment $\boldsymbol {\mu }$ is given by the expression

(5.10)\begin{equation} \boldsymbol{B} = \frac{\mu_0}{4{\rm \pi} r^3} \left( \frac{3 (\boldsymbol{\mu}\boldsymbol{\cdot}\boldsymbol{r}) \boldsymbol{r}}{r^2} - \boldsymbol{\mu} \right), \end{equation}

where $\mu _0$ is the magnetic permeability. For a magnetic moment aligned with the $z$-axis, the magnetic field in the equatorial plane $xOy$ is purely vertical and given by

(5.11)\begin{equation} B_{\textrm{z}} = B \frac{R^3}{r^3}, \end{equation}

where $B$ is the field strength at a distance $R$ from the origin.

For non-relativistic particles, the magnetic gradient drift velocity depending on the perpendicular velocity $v_\perp$ is given by (Baumjohann & Treumann Reference Baumjohann and Treumann1996)

(5.12)\begin{equation} \boldsymbol{v}_{\nabla B} = \frac{m v_\perp^2}{2qB} \frac{\boldsymbol{B} \wedge \nabla B}{B^2} . \end{equation}

For the above magnetic dipole geometry, in the equatorial plane we find a magnetic gradient velocity of

(5.13)\begin{equation} \boldsymbol{v}_{\nabla B} = \mp \frac{3}{2} \frac{m v_\perp^2}{qB r} \boldsymbol{e}_{\varphi} = \mp \frac{3}{2} \frac{v_\perp^2}{\omega_{B}r} \boldsymbol{e}_{\varphi} = \mp \frac{3}{2} v_\perp \frac{r_\textrm{L}}{r} \boldsymbol{e}_{\varphi}, \end{equation}

solely directed into the azimuthal direction and decreasing with distance as $1/r$.

Magnetic gradient drift motion is important in fields significantly varying on a length scale comparable to the Larmor radius. We study such motion in the equatorial plane of a magnetic dipole. We know that the particle must gyrate around the origin where the dipole is located. This motion is induced by the magnetic gradient drift in the azimuthal direction. The particle stays within two circles. An example of this drift motion is shown in figure 10 for $B=10^{3}$ and $\gamma \approx 71$. The time step is varied and taken such that $\log _{10} (\omega _{B} \Delta \tau ) = \{0,-1,-2,-3\}$. The Lorentz factor is conserved as is easily checked. The observed drift speed is compared to the expected drift speed in figure 11. For a gyromotion well resolved in the proper time ($\omega _{B} \Delta \tau \ll 1$), the expected trajectory and drift speed are well reproduced.

Figure 10.

Orbit of an electron in the equatorial plane of a magnetic dipole for $B=10^{3}$ and $\gamma \approx 71$.

Figure 11.

Observed and expected drift speed of a positron in the equatorial plane of a magnetic dipole for $B=10^{3}$ and $\gamma \approx 71$.

A second more stringent example is given by $B=10^{8}$ and $\gamma \approx 10^{6}$. The corresponding trajectory is shown in figure 12 for its guiding centre. The particle gyrates at very high frequency, turning millions of rotations before significantly translating its guiding centre a distance of the order of its Larmor radius. To get the true trajectory without stroboscopic effects, we reduce the time step to $\omega_B\,\Delta\tau =10^{-9}$. A zoom into the plot shows clearly the gyration, figure 13. However, in such a case, it is absolutely useless to resolve the gyration frequency. The time step can be increased by several orders of magnitude without losing accuracy. This is the main advantage of our new particle solver, a key feature in solving kinetic problems around strongly magnetized neutron stars. The drifting velocity also corresponds to the expected value, figure 14.

Figure 12.

Orbit of the guiding centre of an electron in the equatorial plane of a magnetic dipole for $B=10^{8}$ and $\gamma \approx 10^{6}$.

Figure 13.

Zoom into the orbit of an electron in the equatorial plane of a magnetic dipole for $B=10^{8}$ and $\gamma \approx 10^{6}$.

Figure 14.

Observed and expected drift speed of a positron in the equatorial plane of a magnetic dipole for $B=10^{8}$ and $\gamma \approx 10^{6}$.

We finish our extensive test of the algorithm by considering two time varying electromagnetic field configurations represented by linearly and circularly polarized plane waves.

5.7. Linearly polarized plane wave

Consider a linearly polarized plane wave propagating along the $\boldsymbol {e}_{x}$ direction such that the 4-vector potential is $A^\alpha =(0,0,({E}/{\omega })\cos \xi ,0)$. The wave vector is therefore $K^\alpha =({\omega }/{c},k,0,0)$ from which we deduce the phase $\xi =\omega t-kx$. The electromagnetic field is then given by

(5.14a)\begin{gather} \boldsymbol{E} = E \sin \xi \boldsymbol{e}_{y}, \end{gather}

(5.14b)\begin{gather}\boldsymbol{B} = \frac{E}{c} \sin \xi \boldsymbol{e}_{z} . \end{gather}

Initially the particle is at rest with a 4-velocity $u^\alpha _0=(c,\boldsymbol 0)$. Introducing the strength parameter of the wave by the ratio

(5.15)\begin{equation} a = \frac{qE}{m c\omega} \end{equation}

the 4-velocity has components

(5.16a)\begin{gather} u^x = \frac{a^2}{2} c (\cos \xi - 1 )^2, \end{gather}

(5.16b)\begin{gather}u^y = -a c (\cos \xi - 1 ), \end{gather}

(5.16c)\begin{gather}u^0 = c + u^x . \end{gather}

More generally, exact analytical solutions for plane electromagnetic waves in a vacuum (linearly or circularly polarized) have been derived. The methodology using 4-vectors and tensors can be found for instance in Uzan & Deruelle (Reference Uzan and Deruelle2014). See also Michel & Li (Reference Michel and Li1999) for typical applications to pulsars. For brevity, we do not reproduce these computations in this work. The mean spatial velocity becomes

(5.17a)\begin{gather} \langle v^x\rangle = \frac{3}{4} \frac{a^2 c}{1+3a^2/4}, \end{gather}

(5.17b)\begin{gather}\langle v^y\rangle = \frac{a c}{1+3a^2/4} . \end{gather}

After integration, assuming the particle starts at rest at the origin at $\xi =0$, we find

(5.18a)\begin{gather} \omega x = \frac{a^2 c}{8} (6\xi - 8 \sin\xi + \sin 2\xi), \end{gather}

(5.18b)\begin{gather}\omega y = a c (\xi - \sin\xi), \end{gather}

(5.18c)\begin{gather}\omega c t = c\xi + \omega x . \end{gather}

The particle motion is completely described by the strength parameter $a$ (disregarding the initial conditions that are not part of the physical parameters). Examples of motion along the $x$ axis are shown in figure 15 for a series of mildly and ultra-relativistic strength parameters $a=10^i$ with $i\in \{0,3,6,9,12,15\}$. The associated Lorentz factor time evolution is shown in figure 16. The numerical integration is compared to the analytical solution depicted by coloured symbols. Both are in perfect agreement. We are able to simulate acceleration to Lorentz factors well above $\gamma = 10^{15}$. This is compulsory to faithfully study lepton acceleration in neutron star magnetospheres. The true motion in this plane wave is perfectly periodic with a period given by

(5.19)\begin{equation} \omega \, T_{\rm lin} = 2 {\rm \pi}\left( 1 + \frac{3}{4} a^2 \right) \end{equation}

the particle returning to rest at each period $T_{\textrm {lin}}$. The maximum Lorentz factor, reached at $t=T_{\textrm {lin}}/2$ is

(5.20)\begin{equation} \gamma_{\textrm{max}} = 1 + 2 a^2. \end{equation}

Arefiev et al. (Reference Arefiev, Cochran, Schumacher, Robinson and Chen2015) performed similar simulations in a strong plane electromagnetic wave and found significant discrepancies between analytical and numerical results already with strength parameters $a\approx 25$. Moreover the time step criterion for the particle pusher they found would be too restrictive in the case of realistic pulsars or magnetars. In order to compare our results with their findings, we computed the invariant quantity (which is another expression of (5.16c))

(5.21)\begin{equation} \gamma m_e c - p_\textrm{x} = m_e c . \end{equation}

We found that this invariant is well conserved for the full time of integration for any strength parameter giving always the value $1$ in normalized units. We show it, however, in a different manner compared to Arefiev et al. (Reference Arefiev, Cochran, Schumacher, Robinson and Chen2015), instead of the above invariant, we plotted the Lorentz factors in figure 16. An estimate of the error is shown in figure 19. Here also our scheme is second order in time.

Figure 15.

Motion of an electron in a linearly polarized plane wave for different strength parameters $a=10^i$ with $i\in \{0,3,6,9,12,15\}$. Symbols correspond to the analytical solution given in (5.18).

Figure 16.

Lorentz factor of an electron in a linearly polarized plane wave for different strength parameters $a=10^i$ with $i\in \{0,3,6,9,12,15\}$. Symbols correspond to the analytical solution (5.16).

5.8. Circularly polarized plane wave

Consider a circularly polarized plane wave propagating in the $\boldsymbol {e}_{x}$ direction such that the vector potential has components $A^\alpha =(0,0,({E}/{\omega })\cos \xi ,({E}/{\omega })\sin \xi )$ and the wave vector $K^\alpha =(\frac {\omega }{c},k,0,0)$ thus the phase $\xi =\omega \,t-k\,x$. The electromagnetic field is then given by

(5.22a)\begin{gather} \boldsymbol{E} = E ( \sin \xi \boldsymbol{e}_{y} - \cos \xi \boldsymbol{e}_{z} ), \end{gather}

(5.22b)\begin{gather}\boldsymbol{B} = \frac{E}{c} ( \sin \xi \boldsymbol{e}_{z} + \cos \xi \boldsymbol{e}_{y} ) . \end{gather}

Initially, the particle is at rest with 4-velocity $u^\alpha _0=(c,\boldsymbol 0)$. The time evolution of the components of this 4-velocity will be

(5.23a)\begin{gather} u^x = a^2 c (1 - \cos \xi ) = a u^y, \end{gather}

(5.23b)\begin{gather}u^y = a c (1 - \cos \xi ), \end{gather}

(5.23c)\begin{gather}u^z = - a c \sin \xi, \end{gather}

(5.23d)\begin{gather}u^0 = c + u^x . \end{gather}

The mean spatial velocity becomes

(5.24a)\begin{gather} \langle v^x\rangle = \frac{a^2 c}{1+a^2}, \end{gather}

(5.24b)\begin{gather}\langle v^y\rangle = \frac{a c}{1+a^2}, \end{gather}

(5.24c)\begin{gather}\langle v^z\rangle = 0 . \end{gather}

After integration, assuming the particle starts at rest at the origin at phase $\xi =0$, we find

(5.25a)\begin{gather} \omega x = a^2 c (\xi - \sin\xi), \end{gather}

(5.25b)\begin{gather}\omega y = a c (\xi - \sin\xi), \end{gather}

(5.25c)\begin{gather}\omega z = a c (\cos\xi-1), \end{gather}

(5.25d)\begin{gather}\omega c t = c\xi + \omega x . \end{gather}

Examples of motion along the $x$ axis are shown in figure 17 for mildly and ultra-relativistic strength parameters $a=10^i$ with $i\in \{0,3,6,9,12,15\}$. The corresponding evolution of the Lorentz factor is given in figure 18. The numerical integration is compared to the analytical solution depicted by coloured symbols. Here again, both are in perfect agreement. We are even able to push the Lorentz factor limit well above $\gamma = 10^{15}$. Circular polarization is more efficient in accelerating particles to ultra-relativistic speeds.

Figure 17.

Motion of a positron in a circularly polarized plane wave for different strength parameters $a=10^i$ with $i\in \{0,3,6,9,12,15\}$. Symbols correspond to the analytical solution (5.25).

Figure 18.

Lorentz factor of a positron in a circularly polarized plane wave for different strength parameters $a=10^i$ with $i\in \{0,3,6,9,12,15\}$. Symbols correspond to the analytical solution (5.23).

The true motion in this plane wave is also perfectly periodic with a period given by

(5.26)\begin{equation} \omega \, T_{\text{circ}} = 2{\rm \pi}\left( 1 + a^2 \right) \end{equation}

the particle returning to rest at each period $T_{\textrm {circ}}$. The maximum Lorentz factor, reached at $t=T_{\textrm {circ}}/2$ is the same as (5.20). Figure 19 again shows the error in the Lorentz factor for a circularly polarized wave, second order in time still holds. It overlaps with the linearly polarized wave and is indistinguishable. The corresponding relative error in the $x$ position of the particle at $t=T_{\textrm {circ}}/2$ is shown in figure 20 again for a linearly and a circularly polarized wave, second order in time still holds, we even observed almost third order for the circularly polarized wave. Note that for the highest strength parameters $a$, in order to increase the precision of the computation, we tried a precision exceeding that of standard built-in types such as ‘double precision’ codes written in C++. We used the boost multiprecision library with cpp_dec_float_50 format with only modest additional computational cost less than a factor two and only small changes in the C++ code.

Figure 19.

Relative error $\Delta \gamma /\gamma$ on the Lorentz factor at time $t=T/2$ corresponding to the maximum Lorentz factor reachable, depending on observer time step $\Delta t$ on log–log scale for the linearly, solid line, and circularly, dashed line, polarized waves and strength parameter $a$. The black solid line triangles correspond to a $\Delta t^{-2}$ slope.

Figure 20.

Relative error $\Delta x/x$ on the particle position at time $t=T/2$ corresponding to the maximum Lorentz factor reachable, depending on observer time step $\Delta t$ on log–log scale for the linearly, solid line, and circularly, dashed line, polarized waves and strength parameter $a$. The black solid line triangles correspond to a $\Delta t^{-2}$ slope and a $\Delta t^{-3}$ slope.

Eventually, we show an implementation of our new pusher in a one-dimensional (1-D) electromagnetic PIC code to demonstrate the feasibility of our approach, switching back to standard double precision built-in types. We emphasize, however, that, for realistic neutron star applications, accuracy for large strength parameters is compulsory and therefore going beyond ‘double precision’ must not be discarded. But in depth analysis of this requirement is left to another work as described by a companion paper (Tomczak & Pétri Reference Tomczak and Pétri2020) specifically dealing with particle acceleration in ultra-strong electromagnetic fields of neutron stars.

6.1. Plasma oscillations

When a charge density $\rho _{s}=n_{s} q_{s}$ grows in a plasma, it tries to relax to its equilibrium state by building a restoring force induced by the electric field. Combined with its inertia, the plasma enters into an oscillatory motion described by plasma oscillations at the relativistic plasma frequency

(6.4)\begin{equation} \omega_{\textrm{p}_{s}}^2 = \frac{n_{s} q_{s}^2}{\langle \varGamma_{s}\rangle \,m_{s}\,\varepsilon_0} , \end{equation}

where $\langle \varGamma _{s}\rangle$ is the average bulk Lorentz factor of the moving plasma slab. For a cold plasma, no wave can propagate but for a warm plasma, propagation is permitted in the form of Langmuir waves.

Plasma oscillations are purely electrostatic problems reducing Maxwell equations to the Maxwell–Gauss part

(6.5)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{E} = \frac{\rho_{s}}{\varepsilon_0} . \end{equation}

In our electrostatic PIC version of the code, we integrate this equation by fast Fourier transform techniques assuming periodic boundary conditions in the $z$ direction. In order to initiate plasma oscillations, we periodically perturb the initial position of the electrons evolving in an otherwise immobile net of ions. The motion becomes relativistic if the electric field energy density $\varepsilon _0 E^2/2$ induced by the perturbation becomes equal to the rest mass energy density $n_{s}m_{s} c^2$.

In figure 21, an example of non-relativistic cold plasma oscillations shows the good agreement between the expectation and the simulations. A second example of a relativistic cold plasma oscillations is shown in figure 22. Here also, the simulations agree with (6.4). The temporal evolution of the 4-velocity $u_{\textrm {z}}$ is fitted with a sinusoidal law

(6.6)\begin{equation} u_{\textrm{z}} = A u_{z\textrm{max}} \sin(\omega t) . \end{equation}

The corresponding average bulk Lorentz factor is $\langle \varGamma _{s}\rangle \approx 31$ whereas its maximum is $\langle \varGamma _{s}\rangle _{\textrm {max}} \approx 50$.

Figure 21.

Non-relativistic plasma oscillations showing the longitudinal velocity $\beta _{\textrm {z}}$ normalized to its maximum value $\beta _{\textrm zmax} \ll 1$, simulations in blue points and theory in red solid line.

Figure 22.

Relativistic plasma oscillations showing the longitudinal 4-velocity $u_{\textrm {z}}$ with a mean factor $\langle \varGamma _{s}\rangle \approx 31$, simulations in blue points and fit in red solid line.

6.2. Two-stream instability

When two beams move through each other, an instability known as the two-stream instability grows, tapping some energy from the particle kinetic energy into the electric field. This instability is ubiquitous when beams of different bulk velocity penetrate each other. The dispersion relation for two counter-propagating streams of bulk velocity $v_{b}$ and $-v_{b}$ is

(6.7)\begin{equation} \frac{\omega_{p_-}^2}{\gamma_{b}^3(\omega-k v_{b})^2} + \frac{\omega_{p_+}^2}{\gamma_{b}^3(\omega+k v_{b})^2} = 1, \end{equation}

with $\omega$ the wave frequency and $k$ its wavenumber. The solution is given for an electron–positron plasma with $\omega _{p_-}=\omega _{p_+}=\omega _{p}$ and $X=\gamma _{b}^3 ({k^2 v_{b}^2}/{\omega _{p}^2})>0$ by

(6.8)\begin{equation} \gamma_{b}^3 \frac{\omega^2}{\omega_{p}^2} = 1 + X \pm \sqrt{1+4X} . \end{equation}

The counter-streaming plasma is unstable whenever $\omega ^2<0$ corresponding to $X<2$. The minimum is reached at $X=3/4$ and it equals $-1/4$. The maximum growth rate $\omega _\textrm {i} = \textrm {Im}(\omega )$ and its corresponding wavenumber are

(6.9a,b)\begin{equation} \frac{\omega_\textrm{i}}{\omega_{p}} = \frac{1}{2 \gamma_{b}^{3/2}}, \quad \frac{k c}{\omega_{p}} = \frac{\sqrt{3}}{2 \gamma_{b}^{3/2} v_{b}/c} . \end{equation}

An example of relativistic two-stream instability is shown in figure 23 for relativistic counter-streaming beams with a Lorentz factor each of $\varGamma _{b}=10$. Simulations results show the total electrostatic energy $\varepsilon _0 E^2/2$ in blue, compared with the expected maximum growth rate in red and a fit of the growth rate in green. The agreement between theory and simulations is good but some discrepancies are observed because our periodic simulation box is not exactly a multiple of the fastest growing wavelength. Therefore the simulation does not pick up rigorously the mode in (6.9a,b).

Figure 23.

Relativistic two-stream instability with $\varGamma _{b}=10$, simulations in blue points, theory in red solid line and fit in green.

6.3. Shock

Finally, for a case evolving the full electromagnetic field equations, we study a relativistic shock in one dimensional by colliding an electron/positron pair plasma onto a wall at rest in the simulation box.

Our simple 1-D electromagnetic PIC code follows the ideas exposed in Birdsall & Langdon (Reference Birdsall and Langdon2005) to solve the Maxwell equations by choosing the grid spacing $\Delta x$ equal to the time step $\Delta t$ times the speed of light $\Delta x = c\Delta t$. This technique has proven very useful to study strongly magnetized and relativistic shocks in pulsar winds, as already done by Lyubarsky (Reference Lyubarsky2005) and Pétri & Lyubarsky (Reference Pétri and Lyubarsky2007).

The initial flow is set up with an electric field $\boldsymbol {E} = E \boldsymbol {e}_{y}$ and a magnetic field $\boldsymbol {B} = B \boldsymbol {e}_{z}$ with $cB>E>0$ inducing an electric drift in the positive $\boldsymbol {e}_{x}$ direction. Particles are reflected on a solid wall located at $x=L$ and there is no incoming flux of particles from the left $x=0$. The simulation frame therefore corresponds to the downstream plasma frame. In this frame, the upstream flow has a velocity $\boldsymbol {V}=E/B \boldsymbol {e}_{x}$ with Lorentz factor $\varGamma =(1-V/c)^{-1/2}$. To compute the time step in the drifting frame, we assume that particles have constant velocity during the integration. For a good guess we can use (4.3) for the relation between observer frame $\Delta t$ and drifting frame $\Delta t'$ time steps. More sophisticated methods could release this assumption but we will show that this approximation already gives accurate results.

The magnetization parameter in the upstream flow defined (within a factor unity) by the ratio between the magnetic energy density over the plasma energy density is given by

(6.10)\begin{equation} \sigma = \frac{B^2}{\mu_0\varGamma nm c^2} . \end{equation}

A strongly magnetized flow implies $\sigma \gg 1$. In this limit, the downstream plasma density, temperature and Lorentz factor are summarized in table A.1 of Pétri & Lyubarsky (Reference Pétri and Lyubarsky2007). As shown in this paper, for strong magnetizations and relativistic shock speeds, the compression ratio between the downstream plasma density $n_2'$ and the upstream plasma density $n_{1/2}$ as measured in the downstream plasma is

(6.11)\begin{equation} \frac{n_2'}{n_{1/2}} = \frac{1+\beta_2}{\beta_2} \approx 2. \end{equation}

Therefore, the compression ratio is almost 2 if the downstream flow remains ultra relativistic with $\beta _2\approx 1$.

A example of such a relativistic shock is shown in figure 24 for the upstream and downstream electromagnetic field $(E_y,B_z)$ and the particle density number $n$, all normalized to their respective values. As a guide to the eye, the lines $y=1$ and $y=2$ are also shown in grey. The flow bulk Lorentz factor is $\varGamma \approx 70$ and the magnetization is $\sigma \approx 70$. In this shock, the magnetohydrodynamics (MHD) jump conditions are satisfied with the magnetic flux conservation property stipulating that $B/n$ is conserved. This is indeed observed in the simulations. Inspecting figure 24, the compression ratio is indeed 2 as expected from an ultra-relativistic strongly magnetized MHD shock. The downstream electric field vanishes because the electric drift must disappear in this plasma, forced to be at rest by the solid wall boundary condition. Such flows are met for instance in pulsar wind nebulae where $\varGamma \gg 1$ and $\sigma \gg 1$ before entering the termination shock (Kirk, Lyubarsky & Petri Reference Kirk, Lyubarsky and Petri2009). The simulation has been stopped when the back propagating electromagnetic wave has reach the middle of the simulation box. This is because there is no incoming particle flow from the left at $x=0$, therefore we observe a sharp jump in the density at $x=L/2$ from $n/n_0=1$ to $n/n_0=0$.

Figure 24.

Upstream and downstream electromagnetic field $(E_y,B_z)$ and particle density number $n$ as a function of position $x$ for an ultra-relativistic strongly magnetized MHD shock with $\varGamma \approx 70$ and $\sigma \approx 70$. All quantities are normalized with respect to their upstream value.