1 Introduction
Relativistic magnetized flows arise in many astrophysical contexts, especially around compact objects such as neutron stars and black holes. Plasma regimes are typically collisionless and comprise a substantial fraction of electron–positron pairs. Numerical techniques to investigate the behaviour of such plasmas require a kinetic description. Unfortunately solving the full Vlasov equations in six dimensions (3 space and 3 velocity directions) on a computer remains out of reach with current technology. Therefore particle in cell (PIC) methods are preferred and well suited for numerically solving collisionless plasma physics problems where kinetic effects prevail. Nevertheless, particle methods are inherently prone to numerical noise. Moreover for ultra-strong magnetic fields such as those anchored into pulsars and magnetars, the time scale of gyration is much too small compared to the long-term evolution of the plasma rendering computations on dynamical time scales impossible.
The magnetospheres of compact objects, black holes and neutron stars, are filled with relativistic electron–positron pairs travelling in a ultra-strong magnetic field close to or even above the critical magnetic field of
$4.4\times 10^{9}$
T where quantum electrodynamics processes become dominant. Moreover the plasma evolves in an almost perfectly collisionless regime where Coulombian collisions can be neglected. These plasmas are prone to several instabilities and reconnection phenomena as well as radiative processes that need to be studied within a kinetic description of the fluid.
The Vlasov approach represents the best mean to accurately investigate the long-term photon–plasma bath evolution. However it requires integration in a six-dimensional phase space, computationally very expensive and unfortunately impossible to treat on contemporary computers. Another less demanding alternative is given by PIC simulations that sample the phase space into millions or billions of individual particles evolving in time via macroscopic electromagnetic interactions (no particle–particle interactions taken into account in the simplest approach). Both techniques require a numerical integration of the particle equation of motion in an external electromagnetic field. Performance, advantages and drawbacks of PIC methods are extensively discussed in the standard textbooks on particle simulations of plasmas such as Hockney & Eastwood (Reference Hockney and Eastwood1988) and Birdsall & Langdon (Reference Birdsall and Langdon2005). In any case, it is usually difficult to solve this problem in the ultra-relativistic regime where particles can reach Lorentz factors larger than several millions or billions
$\unicode[STIX]{x1D6FE}>10^{9}$
. At these energies, radiation reaction is expected to play a prominent role on the long-term dynamics. In the non-relativistic regime the Boris method (Boris Reference Boris1970) and the leapfrog method are classical explicit schemes often quoted. Boris algorithm has been extended to relativistic motion by Vay (Reference Vay2008). A volume-preserving relativistic integrator was recently developed by Higuera & Cary (Reference Higuera and Cary2017). Splitting techniques for high-order symmetric volume-preserving methods can be found in He et al. (Reference He, Sun, Zhang, Wang, Liu and Qin2016). A relativistic version of the leapfrog method was also implemented in the code of Melzani et al. (Reference Melzani, Winisdoerffer, Walder, Folini, Favre, Krastanov and Messmer2013) but with possible issues for ultra-relativistic speeds. Realistic simulations have to take care of exact particle energy conservation which is a concern for all available methods (Lapenta & Markidis Reference Lapenta and Markidis2011). Flaws already encountered in non-relativistic plasma simulations are enhanced in relativistic regimes. Vay & Godfrey (Reference Vay and Godfrey2014) recently reviewed about these relativistic PIC methods.
Simulations of plasmas around neutron stars are particularly demanding because of the extreme conditions in the vicinity of such stars. Indeed a typical ratio between the cyclotron frequency
$\unicode[STIX]{x1D714}_{B}$
and the stellar rotation frequency
$\unicode[STIX]{x1D6FA}$
is
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D714}_{B}}{\unicode[STIX]{x1D6FA}}=\frac{eB}{m\unicode[STIX]{x1D6FA}}\approx 10^{16-19}.\end{eqnarray}$$
Moreover the ratio between neutron star radius
$R$
and Larmor radius
$r_{B}$
is also large where
$e$
is the electron charge,
$m$
its mass and
$B$
the magnetic field strength.
$$\begin{eqnarray}\frac{R}{r_{B}}=\frac{R\unicode[STIX]{x1D714}_{B}}{c}=\frac{\unicode[STIX]{x1D714}_{B}}{\unicode[STIX]{x1D6FA}}\frac{R}{r_{L}}\approx \frac{10^{15}}{\unicode[STIX]{x1D6FE}}\frac{R}{r_{L}},\end{eqnarray}$$
where
$r_{L}=c/\unicode[STIX]{x1D6FA}\approx 10-10.000R$
is the light-cylinder radius, i.e. another important length scale indicating the transition between the static quasi-stationary and the wave zone induced by rotating the magnet anchored into the neutron star.
$\unicode[STIX]{x1D6FE}$
is the particle Lorentz factor. Leptons extracted from the neutron star surface are expected to start with modest kinetic energy,
$\unicode[STIX]{x1D6FE}\gtrsim 1$
, thus the most severe constrain on this ratio is obtained for
$\unicode[STIX]{x1D6FE}=1$
. Such huge gaps need to be taken into account for realistic modelling of pulsars. In this communication, we show how to efficiently and accurately solve for this motion for any Lorentz factor up to at least
$\unicode[STIX]{x1D6FE}=10^{9}$
in some cases. We first remind how to implicitly solve the equation of motion in the non-relativistic and relativistic case and how update implicitly the particle position. Then we test our algorithm in some special configurations with constant external electromagnetic fields and with spatially and temporally varying fields.
2 Equation of motion
The equation of motion for a relativistic charged particle of mass
$m$
and charge
$q$
, evolving in an external electromagnetic field with
$\boldsymbol{E}$
the electric part and
$\boldsymbol{B}$
the magnetic part, is given by the Lorentz force according to
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{r}}{\text{d}t}=\boldsymbol{v} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{p}}{\text{d}t}=q(\boldsymbol{E}+\boldsymbol{v}\wedge \boldsymbol{B}). & \displaystyle\end{eqnarray}$$
$\boldsymbol{p}$
is deduced from the spatial three-velocity
$\boldsymbol{v}$
or from the four-velocity
$\boldsymbol{u}$
by 
$$\begin{eqnarray}\boldsymbol{p}=\unicode[STIX]{x1D6FE}m\boldsymbol{v}=m\boldsymbol{u}.\end{eqnarray}$$
The particle position is denoted by the vector
$\boldsymbol{r}$
and
$t$
is the time as measured by an observer attached to a frame at rest. In order to show the complication arising in the relativistic regime, we first remind the non-relativistic implicit algorithm to solve for the particle velocity. Then we extend to the relativistic case and show how to deal with the Lorentz factor.
2.1 Non-relativistic motion
In the non-relativistic limit, the Lorentz factor is
$\unicode[STIX]{x1D6FE}=1$
and the momentum simplifies to
$\boldsymbol{p}=m\boldsymbol{v}$
. The equation of motion therefore becomes
$$\begin{eqnarray}\frac{\text{d}\boldsymbol{v}}{\text{d}t}=\frac{q}{m}(\boldsymbol{E}+\boldsymbol{v}\wedge \boldsymbol{B}).\end{eqnarray}$$
We look for an implicit method to solve this equation. The discretized version of the Lorentz force (2.3) is
$$\begin{eqnarray}\frac{\boldsymbol{v}^{n+1/2}-\boldsymbol{v}^{n-1/2}}{\unicode[STIX]{x0394}t}=\frac{q}{m}\left(\boldsymbol{E}^{n}+\frac{\boldsymbol{v}^{n+1/2}+\boldsymbol{v}^{n-1/2}}{2}\wedge \boldsymbol{B}^{n}\right).\end{eqnarray}$$
We rearrange terms in order to collect expressions containing only the unknown velocity and the next time step
$\boldsymbol{v}^{n+1/2}$
by writing
$$\begin{eqnarray}\boldsymbol{v}^{n+1/2}-\boldsymbol{v}^{n+1/2}\wedge \frac{q\unicode[STIX]{x0394}t\boldsymbol{B}^{n}}{2m}=\boldsymbol{v}^{n-1/2}+\boldsymbol{v}^{n-1/2}\wedge \frac{q\unicode[STIX]{x0394}t\boldsymbol{B}^{n}}{2m}+\frac{q\unicode[STIX]{x0394}t\boldsymbol{E}^{n}}{m}.\end{eqnarray}$$
This system can be solved explicitly in terms of
$\boldsymbol{v}^{n+1/2}$
by introducing two constant parameters such that
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}=\frac{q\unicode[STIX]{x0394}t}{m};\quad p=\frac{\unicode[STIX]{x1D6FC}}{2}. & & \displaystyle\end{eqnarray}$$
We also define the matrix for the magnetic rotation
$$\begin{eqnarray}\unicode[STIX]{x1D64D}^{n}=\left[\begin{array}{@{}ccc@{}}0 & pB_{z}^{n} & -pB_{y}^{n}\\ -pB_{z}^{n} & 0 & pB_{x}^{n}\\ pB_{y}^{n} & -pB_{x}^{n} & 0\end{array}\right]\end{eqnarray}$$
and the matrix for the electric acceleration
$$\begin{eqnarray}\unicode[STIX]{x1D64E}^{n}=\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D6FC}E_{x}^{n}\\ \unicode[STIX]{x1D6FC}E_{y}^{n}\\ \unicode[STIX]{x1D6FC}E_{z}^{n}\end{array}\right].\end{eqnarray}$$
In matrix notation, the discretized equation of motion becomes
$$\begin{eqnarray}(I_{3}-\unicode[STIX]{x1D64D}^{n})V^{n+1/2}=(I_{3}+\unicode[STIX]{x1D64D}^{n})V^{n-1/2}+\unicode[STIX]{x1D64E}^{n}.\end{eqnarray}$$
This system is solved for the unknown vector
$V^{n+1/2}$
to deduce the velocity at the next time
$t^{n+1/2}$
by
$$\begin{eqnarray}V^{n+1/2}=(I_{3}-\unicode[STIX]{x1D64D}^{n})^{-1}[(I_{3}+\unicode[STIX]{x1D64D}^{n})V^{n-1/2}+\unicode[STIX]{x1D64E}^{n}].\end{eqnarray}$$
The inverse matrix of
$(I_{3}-\unicode[STIX]{x1D64D}^{n})$
is easily found to be
$$\begin{eqnarray}(I_{3}-\unicode[STIX]{x1D64D}^{n})^{-1}=\frac{1}{1+p^{2}B^{2}}\left[\begin{array}{@{}ccc@{}}1+B_{x}^{2}p^{2} & B_{x}B_{y}p^{2}+B_{z}p & B_{x}B_{z}p^{2}-B_{y}p\\ B_{x}B_{y}p^{2}-B_{z}p & 1+B_{y}^{2}p^{2} & B_{y}B_{z}p^{2}+B_{x}p\\ B_{x}B_{z}p^{2}+B_{y}p & B_{y}B_{z}p^{2}-B_{x}p & 1+B_{z}^{2}p^{2}\end{array}\right].\end{eqnarray}$$
$V^{n+1/2}$
can be computed for any time step at the price of a costly matrix multiplication. The inversion of the matrix in front of
$V^{n+1/2}$
is the crucial step to an implicit integration of the equation of motion. In the relativistic case, the procedure is far from trivial because the Lorentz factor is a priori not known at the next time step. However, as we show in the next section, it is possible to find the Lorentz factor by solving a biquadratic equation.
2.2 Relativistic motion
The relativistic equation of motion is treated in a similar manner as before. However, a complication stems from the indeterminacy of the particle Lorentz factor at the next time
$t^{n+1/2}$
. Indeed, the discretized version of the Lorentz force (2.1b
) is
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FE}^{n+1/2}\boldsymbol{v}^{n+1/2}-\unicode[STIX]{x1D6FE}^{n-1/2}\boldsymbol{v}^{n-1/2}}{\unicode[STIX]{x0394}t}=\frac{q}{m}\left(\boldsymbol{E}^{n}+\frac{\boldsymbol{v}^{n+1/2}+\boldsymbol{v}^{n-1/2}}{2}\wedge \boldsymbol{B}^{n}\right).\end{eqnarray}$$
Written in terms of the four-velocity we get
$$\begin{eqnarray}\frac{\boldsymbol{u}^{n+1/2}-\boldsymbol{u}^{n-1/2}}{\unicode[STIX]{x0394}t}=\frac{q}{m}\left[\boldsymbol{E}^{n}+\frac{1}{2}\left(\frac{\boldsymbol{u}^{n+1/2}}{\unicode[STIX]{x1D6FE}^{n+1/2}}+\frac{\boldsymbol{u}^{n-1/2}}{\unicode[STIX]{x1D6FE}^{n-1/2}}\right)\wedge \boldsymbol{B}^{n}\right].\end{eqnarray}$$
This version has the advantage of being fully symmetric in time, therefore conserves energy and is moreover Lorentz invariant. This is the procedure chosen by Vay (Reference Vay2008). The unknown velocity at time
$t^{n+1/2}$
is separated from the known velocity at time
$t^{n-1/2}$
such that
$\boldsymbol{u}^{n+1/2}$
has to be found according to the relation
$$\begin{eqnarray}\boldsymbol{u}^{n+1/2}-\frac{\boldsymbol{u}^{n+1/2}}{\unicode[STIX]{x1D6FE}^{n+1/2}}\wedge \frac{q\unicode[STIX]{x0394}t\boldsymbol{B}^{n}}{2m}=\boldsymbol{u}^{n-1/2}+\frac{\boldsymbol{u}^{n-1/2}}{\unicode[STIX]{x1D6FE}^{n-1/2}}\wedge \frac{q\unicode[STIX]{x0394}t\boldsymbol{B}^{n}}{2m}+\frac{q\unicode[STIX]{x0394}t\boldsymbol{E}^{n}}{m}.\end{eqnarray}$$
Contrary to the non-relativistic case, we cannot solve explicitly for
$\boldsymbol{u}^{n+1/2}$
because of the Lorentz factor
$\unicode[STIX]{x1D6FE}^{n+1/2}$
that remains unknown. However it is linked to the sought velocity
$\boldsymbol{u}^{n+1/2}$
by
$$\begin{eqnarray}(\unicode[STIX]{x1D6FE}^{n+1/2})^{2}=1+(u^{n+1/2}/c)^{2}.\end{eqnarray}$$
Fortunately, this factor can be determined explicitly analytically following the method outlined in the next lines. Let us introduce the same variables as before defined by (2.6) and two new parameters taking into account the actual Lorentz factor
$$\begin{eqnarray}p^{\pm }=\frac{p}{\unicode[STIX]{x1D6FE}^{n\pm 1/2}}\end{eqnarray}$$
as well as the matrices
$\unicode[STIX]{x1D64D}^{n}$
and
$\unicode[STIX]{x1D64E}^{n}$
. Compared to the non-relativistic situation, we need to replace the matrix
$\unicode[STIX]{x1D64D}^{n}$
by introducing the Lorentz factor to get the useful matrices as
$\unicode[STIX]{x1D64D}^{n}/\unicode[STIX]{x1D6FE}^{n+1/2}$
for
$U^{n+1/2}$
and as
$\unicode[STIX]{x1D64D}^{n}/\unicode[STIX]{x1D6FE}^{n-1/2}$
for
$U^{n-1/2}$
. The four-velocity at time
$t^{n+1/2}$
is solution of the matrix system
$$\begin{eqnarray}\left(I_{3}-\frac{\unicode[STIX]{x1D64D}^{n}}{\unicode[STIX]{x1D6FE}^{n+1/2}}\right)U^{n+1/2}=\left(I_{3}+\frac{\unicode[STIX]{x1D64D}^{n}}{\unicode[STIX]{x1D6FE}^{n-1/2}}\right)U^{n-1/2}+\unicode[STIX]{x1D64E}^{n}.\end{eqnarray}$$
We invert the matrix
$(I_{3}-\unicode[STIX]{x1D64D}^{n}/\unicode[STIX]{x1D6FE}^{n+1/2})$
to obtain
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(I_{3}-\frac{\unicode[STIX]{x1D64D}^{n}}{\unicode[STIX]{x1D6FE}^{n+1/2}}\right)^{-1}=\frac{1}{1+(p^{+})^{2}B^{2}}\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\left[\begin{array}{@{}ccc@{}}B_{x}^{2}(p^{+})^{2}+1 & B_{x}B_{y}(p^{+})^{2}+B_{z}p^{+} & B_{x}B_{z}(p^{+})^{2}-B_{y}p^{+}\\ B_{x}B_{y}(p^{+})^{2}-B_{z}p^{+} & B_{y}^{2}(p^{+})^{2}+1 & B_{y}B_{z}(p^{+})^{2}+B_{x}p^{+}\\ B_{x}B_{z}(p^{+})^{2}+B_{y}p^{+} & B_{y}B_{z}(p^{+})^{2}-B_{x}p^{+} & B_{z}^{2}(p^{+})^{2}+1\end{array}\right]\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \det \left(I_{3}-\frac{\unicode[STIX]{x1D64D}^{n}}{\unicode[STIX]{x1D6FE}^{n+1/2}}\right)=1+(p^{+})^{2}B^{2}. & & \displaystyle\end{eqnarray}$$
The matrix cannot be singular because its determinant is always larger than unity. The four-velocity is derived at the next time step
$t^{n+1/2}$
according to
$$\begin{eqnarray}U^{n+1/2}=\left(I_{3}-\frac{\unicode[STIX]{x1D64D}^{n}}{\unicode[STIX]{x1D6FE}^{n+1/2}}\right)^{-1}\left[\left(I_{3}+\frac{\unicode[STIX]{x1D64D}^{n}}{\unicode[STIX]{x1D6FE}^{n-1/2}}\right)U^{n-1/2}+\unicode[STIX]{x1D64E}^{n}\right].\end{eqnarray}$$
The Lorentz factor
$\unicode[STIX]{x1D6FE}^{n+1/2}$
remains unknown but it can be deduced from the solution of a biquadratic equation. Indeed, from the relation between four-velocity and Lorentz factor as given in (2.15), equation (2.20) allows one to write the Lorentz factor as a root of the biquadratic polynomial defined by
$$\begin{eqnarray}(\unicode[STIX]{x1D6FE}^{n+1/2})^{4}+a_{1}(\unicode[STIX]{x1D6FE}^{n+1/2})^{2}+a_{2}=0.\end{eqnarray}$$
The coefficients are given by
$$\begin{eqnarray}\displaystyle a_{2}=-p^{2}[(B^{n})^{2}+\{(\unicode[STIX]{x1D64E}^{n}+\boldsymbol{u}^{n-1/2})\boldsymbol{\cdot }\boldsymbol{B}^{n}\}^{2}] & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle a_{1} & = & \displaystyle p^{2}(B^{n})^{2}-1-(\unicode[STIX]{x1D64E}^{n}+\boldsymbol{u}^{n-1/2})^{2}-p^{2}(\boldsymbol{v}^{n-1/2}\wedge \boldsymbol{B}^{n})^{2}+2p(\boldsymbol{v}^{n-1/2}\wedge \unicode[STIX]{x1D64E}^{n})\boldsymbol{\cdot }\boldsymbol{B}^{n}\nonumber\\ \displaystyle & = & \displaystyle p^{2}(B^{n})^{2}-1-(\boldsymbol{u}^{n-1/2})^{2}-2\unicode[STIX]{x1D64E}^{n}\boldsymbol{\cdot }\boldsymbol{u}^{n-1/2}-(p\boldsymbol{v}^{n-1/2}\wedge \boldsymbol{B}^{n}+\unicode[STIX]{x1D64E}^{n})^{2}.\end{eqnarray}$$
There are formally four solutions for the Lorentz factor, two of them being negative and rejected, and only one satisfying the physical condition
$\unicode[STIX]{x1D6FE}^{n+1/2}\geqslant 1$
. It is explicitly given by
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}^{n+1/2}=\sqrt{\frac{\displaystyle \sqrt{a_{1}^{2}-4a_{2}}-a_{1}}{2}}.\end{eqnarray}$$
This Lorentz factor is inserted into the matrix
$(I_{3}-(\unicode[STIX]{x1D64D}^{n}/\unicode[STIX]{x1D6FE}^{n+1/2}))^{-1}$
to deduce the four-velocity at the next time step
$t^{n+1/2}$
.
This analytical procedure to solve for the implicit discretization of the momentum equation has already been given by Vay (Reference Vay2008) although in a different form. In order to retrieve his results, we conform to his notation and introduce the following quantities
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D749}=p\boldsymbol{B}^{n} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{t}=\frac{\unicode[STIX]{x1D749}}{\unicode[STIX]{x1D6FE}^{n+1/2}}=p^{+}\boldsymbol{B}^{n} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\prime }=\boldsymbol{u}^{n-1/2}+\frac{\boldsymbol{u}^{n-1/2}}{\unicode[STIX]{x1D6FE}^{n-1/2}}\wedge \unicode[STIX]{x1D749}+\unicode[STIX]{x1D64E}^{n}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\boldsymbol{u}^{n+1/2}=\boldsymbol{u}^{\prime }+\boldsymbol{u}^{n+1/2}\wedge \boldsymbol{t}.\end{eqnarray}$$
Taking the dot and cross product with
$\boldsymbol{t}$
we get
$$\begin{eqnarray}\boldsymbol{u}^{n+1/2}=\frac{1}{1+t^{2}}(\boldsymbol{u}^{\prime }+\boldsymbol{u}^{\prime }\wedge \boldsymbol{t}+(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{t})\boldsymbol{t})\end{eqnarray}$$
as found by Vay (Reference Vay2008). In matrix form this reduces to
$$\begin{eqnarray}\boldsymbol{u}^{n+1/2}=\frac{1}{1+t^{2}}\left[\begin{array}{@{}ccc@{}}1+t_{x}^{2} & t_{x}t_{y}+t_{z} & t_{x}t_{z}-t_{y}\\ t_{x}t_{y}-t_{z} & 1+t_{y}^{2} & t_{y}t_{z}+t_{x}\\ t_{x}t_{z}+t_{y} & t_{y}t_{z}-t_{x} & 1+t_{z}^{2}\end{array}\right],\end{eqnarray}$$
which is the matrix form (2.20). Adapting our notation to expressions given by Vay (Reference Vay2008), our biquadratic equation for the new Lorentz factor reduces exactly to the same expression. Consequently, the analytic expressions we found, although formulated in a matrix form, are in the end exactly the same as in Vay (Reference Vay2008) who uses a vector form instead.
2.3 Position update and fully implicit scheme
The Vay (Reference Vay2008) algorithm treats the velocity part implicitly but evolves the position explicitly in a leapfrog fashion. A more judicious scheme would treat both position and velocity updates implicitly, with a similar discretization. Such modification gets rid of the staggered grids in time where positions are evaluated at full time steps whereas velocities are computed at half-time steps or vice versa. Fully implicit schemes are usually computationally expensive but in the present work we are mostly concerned with precision, accuracy and stability of the integration technique. This is particularly relevant for particles moving in ultra-strong electromagnetic fields like those met in pulsar magnetospheres.
The new scheme we propose is able to handle spatially and temporally varying fields to high accuracy. The position and velocity updates are given similarly by
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\boldsymbol{r}^{n+1}-\boldsymbol{r}^{n}}{\unicode[STIX]{x0394}t}=\boldsymbol{v}^{\ast } & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\unicode[STIX]{x0394}t}=\frac{q}{m}[\boldsymbol{E}(\boldsymbol{r}^{\ast },t^{\ast })+\boldsymbol{v}^{\ast }\wedge \boldsymbol{B}(\boldsymbol{r}^{\ast },t^{\ast })] & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{r}^{\ast }=\frac{\boldsymbol{r}^{n+1}+\boldsymbol{r}^{n}}{2} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}^{\ast }=\frac{\boldsymbol{v}^{n+1}+\boldsymbol{v}^{n}}{2}=\frac{1}{2}\left(\frac{\boldsymbol{u}^{n+1}}{\unicode[STIX]{x1D6FE}^{n+1}}+\frac{\boldsymbol{u}^{n}}{\unicode[STIX]{x1D6FE}^{n}}\right) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle t^{\ast }=\frac{t^{n+1}+t^{n}}{2}. & \displaystyle\end{eqnarray}$$
$t^{\ast }$
between
$t^{n}$
and
$t^{n+1}$
. The external electromagnetic field is evaluated at the mean position
$\boldsymbol{r}^{\ast }$
between
$\boldsymbol{r}^{n}$
and
$\boldsymbol{r}^{n+1}$
.
A direct solution of this nonlinear system is usually too costly. Another technique to solve it uses the Picard iteration method. Starting from a first guess of the solution given by
$\boldsymbol{r}^{n+1}=\boldsymbol{r}^{n}$
and
$\boldsymbol{u}^{n+1}=\boldsymbol{u}^{n}$
, the position and velocity are updated at time
$t^{n+1}$
according to (2.28) until convergence to a prescribed accuracy. When the difference between two successive updates becomes less than a user prescribed tolerance, convergence is achieved. The number of iterations required is usually less than ten.
In the next section we provide some numerical examples of particle trajectories in constant and spatially/temporally varying electromagnetic fields to check the efficiency and accuracy of our algorithm.
3 Special cases
The solution of the system simplifies drastically in special cases where only a magnetic or an electric field is present. For the purpose of numerical tests, we use normalized units and follow the motion of an electron such that its mass is
$m=1$
and its charge
$q=-1$
in these units. The strength of the electric or magnetic field is also normalized to values taken to be in the set
$\{10^{-3},1,10^{3}\}$
. The matrices in the algorithm contain products between different Cartesian components of the electromagnetic field. In order to test extensively the code, we chose to integrate numerically along each coordinate axis
$\{x,y,z\}$
separately but also along a combination of two or even the three coordinate axes in order to get the most general matrices with all elements computed to give non-trivial values. Comparing the simulation outputs with exact analytical expressions is easy in a constant and homogeneous electromagnetic field. Interestingly, some exact expressions are also known for a monopolar electric field and for polarized plane waves. Comparison between numerical and analytical solutions are given below for special electromagnetic field configurations.
3.1 Constant electric field
Let us assume a constant electric field along the
$x$
axis
$\boldsymbol{E}=E\boldsymbol{e}_{x}$
and a particle initially at position
$x=0$
at time
$t=0$
with zero initial velocity. The particle is subject to a constant acceleration (in its rest frame) in the
$x$
direction with intensity
$g=qE/m$
. Moreover its trajectory can be exactly integrated giving
$$\begin{eqnarray}\displaystyle & \displaystyle x=\frac{c^{2}}{g}\left[\sqrt{1+\left(\frac{gt}{c}\right)^{2}}-1\right] & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v=\frac{gt}{\sqrt{1+\displaystyle \left(\frac{gt}{c}\right)^{2}}}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}(t)=\sqrt{1+\left(\frac{gt}{c}\right)^{2}}.\end{eqnarray}$$
At very late times
$t\gg c/g$
the trajectory can be approximated by
$$\begin{eqnarray}\displaystyle & \displaystyle x\approx \text{sign}(g)ct & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v\approx \text{sign}(g)c & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}\approx \frac{|g|t}{c}. & \displaystyle\end{eqnarray}$$
$y$
and
$z$
directions are easily found with similar expressions. For an arbitrary direction of the electric field, if its strength remains constant, the final Lorentz factor remains the same for all simulations if the final time measured in units of
$c/g$
remains identical.
As a diagnostic of the performances of the algorithm, we compute the relative error in the final value of the distance left from the origin and the final Lorentz factor for different time steps
$\unicode[STIX]{x0394}t$
and different electric field strengths. For concreteness, the final time is fixed to
$10^{9}\,c/g$
such that the final Lorentz factor is close to
$\unicode[STIX]{x1D6FE}_{f}=10^{9}$
. The time step is also measured in units of
$c/g$
. The final distance from the origin is estimated from the trajectory and equal to
$\{10^{12},10^{9},10^{6}\}$
for respectively intensities of
$E=\{10^{-3},1,10^{3}\}$
. A summary of final errors for an electric field equally directed along all three axes is given in table 1. The relative error is insensitive to the normalized electric field strength. It is computed for any quantity
$X$
according to
$$\begin{eqnarray}\unicode[STIX]{x1D716}_{X}=\frac{|X_{sim}-X_{th}|}{|X_{th}|},\end{eqnarray}$$
where
$X_{sim}$
is the value extracted from the simulations and
$X_{th}$
the theoretical value. In the table 1 we only show the results for an electric field strength
$E=1$
. For any value of the time step even when
$\unicode[STIX]{x0394}t\gg c/g$
the final Lorentz factor is computed to high accuracy with 11 digits for the largest time step. Decreasing
$\unicode[STIX]{x0394}t$
will increase the number of operations to be done and therefore lower the precision. However the final position is less sensitive to
$\unicode[STIX]{x0394}t$
given with 7–8 digits of precision except for
$\unicode[STIX]{x0394}t=10^{-3}c/g$
where only 5 digits of precision are computed. Obviously, reducing the time step increases the accuracy of the position but the number of iterations also increases and the relative error in the Lorentz factor increases too. For
$\unicode[STIX]{x0394}t=10^{-3}$
the distance is given with 5 digits and the Lorentz factor is also thus 3 digits less for distance and 6 digits less for velocity. Our algorithm is able to catch particles with very high Lorentz factors even with large time steps. A compromise between precision in the position and velocity space and computational cost requires an appropriate choice of the time step.
Table 1. Relative precision for the final position of the particle and for the Lorentz factor in an uniformly accelerating normalized electric field directed along the axis
$\boldsymbol{e}_{x}+\boldsymbol{e}_{y}+\boldsymbol{e}_{z}$
. The particle was initially at rest and its final Lorentz factor is
$\unicode[STIX]{x1D6FE}_{f}=10^{9}$
.
To conclude on pure electric acceleration, an example is provided in figure 1 where the Lorentz factor is shown to increase linearly with time. There are no bounds for the particle energy and the code will at some stage fail to reproduce a Lorentz factor above
$10^{15}$
corresponding to double precision representation of floating numbers.
Figure 1. Acceleration of an electron in a constant electric field along the
$x$
axis with
$E=10^{3}$
. The Lorentz factor is plotted versus time in a log–log scale.
Figure 2. Gyromotion of an electron in a constant magnetic field along the
$z$
axis with
$v_{\bot }=0.999c$
and
$\unicode[STIX]{x1D6FE}=10^{6}$
.
3.2 Constant magnetic field
The second instructive example assumes a constant magnetic field along the
$z$
axis. The relativistic trajectory of a charged particle is well known and given by an helicoidal motion corresponding to a constant speed
$v_{z}$
in the
$z$
direction and a circular path in the plane perpendicular to the magnetic field lines. Therefore the motion is described by
$$\begin{eqnarray}\displaystyle & \displaystyle x=r_{B}\cos \unicode[STIX]{x1D714}_{B}t & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle y=r_{B}\sin \unicode[STIX]{x1D714}_{B}t & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z=v_{z}t, & \displaystyle\end{eqnarray}$$
$r_{B}=v_{\bot }/\unicode[STIX]{x1D714}_{B}$
, the synchrotron frequency is
$\unicode[STIX]{x1D714}_{B}=qB/\unicode[STIX]{x1D6FE}m$
and appropriate initial conditions for the position and velocity should be used. The Lorentz factor is obviously constant because the magnetic field does not work and is equal to 
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}=\left(1-\frac{v_{\bot }^{2}+v_{z}^{2}}{c^{2}}\right)^{-1/2}.\end{eqnarray}$$
The solutions for a magnetic field directed along the
$x$
and
$y$
axes are obtained by a permutation of the components of the position vector. The most general configuration requires magnetic field components along all three axis.
For the numerical tests, the time step
$\unicode[STIX]{x0394}t$
is given in units of
$2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{B}$
and the particle makes 100 turns. For the numerical parameters, we used the same field intensities and time steps as for the constant electric field case. The perpendicular speed is
$v_{\bot }=0.999c$
and the constant Lorentz factor is set to
$\unicode[STIX]{x1D6FE}=10^{6}$
from which we deduce
$v_{z}$
according to (3.6). An example of gyromotion is shown in figure 2. The trajectory remains perfectly circular, there is no loss of energy and no shrinking of the orbit after 100 turns. The Lorentz factor is shown in figure 3 and is perfectly constant during the full time of the run and equal to its initial value of
$10^{6}$
to very good accuracy. As a more precise check of the algorithm, we compute the relative errors in the distance travelled along the magnetic axis and the final Lorentz factor as shown in table 2. In this run, we assume a constant magnetic field along the vector
$\boldsymbol{e}_{x}+\boldsymbol{e}_{y}+\boldsymbol{e}_{z}$
. Whatever the time step, the position of the guiding centre is given with almost full accuracy of double precision, 13 digits. The Lorentz factor also remains constant within 13 digits of precision. This is achieved for any time step. Longer integration time does not degrade this accuracy.
Figure 3. Lorentz factor of an electron in a constant magnetic field along the
$z$
axis with
$v_{\bot }=0.999c$
and
$\unicode[STIX]{x1D6FE}=10^{6}$
. The Lorentz factor is plotted versus time in a log–log scale and remains constant.
Table 2. Relative precision for the final position of the guiding centre and for the Lorentz factor in a constant magnetic field directed along the axis
$\boldsymbol{e}_{x}+\boldsymbol{e}_{y}+\boldsymbol{e}_{z}$
. The Lorentz factor is initially
$\unicode[STIX]{x1D6FE}=10^{6}$
.
3.3 Cross 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 an 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 certainty, 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_{y}/B_{z})\boldsymbol{e}_{x}$
. A Lorentz transformation of the electromagnetic field with Lorentz factor
$\unicode[STIX]{x1D6E4}_{E}=1/\sqrt{1-v_{E}^{2}/c^{2}}$
along
$\boldsymbol{v}_{E}$
shows that in the comoving frame
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{E}^{\prime }=0 & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}^{\prime }=\boldsymbol{B}/\unicode[STIX]{x1D6E4}_{E}. & \displaystyle\end{eqnarray}$$
$\boldsymbol{B}^{\prime }$
as seen in the electric drift frame. These quantities are then transformed to the observer inertial frame according to Lorentz transformations for the three-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 
$$\begin{eqnarray}\displaystyle & \displaystyle x^{\prime }=\unicode[STIX]{x1D6E4}_{E}(x-v_{E}t) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle y^{\prime }=y & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z^{\prime }=z. & \displaystyle\end{eqnarray}$$
$E=0.999999$
and that of the magnetic field is
$B=1$
leading to a Lorentz factor of the electric drift motion of
$\unicode[STIX]{x1D6E4}_{E}\approx 707.107$
. The Lorentz factor in the electric drift frame is fixed to
$\unicode[STIX]{x1D6FE}=100$
. Typical results are depicted in figure 4 for the trajectory in the electric drift frame which is usually an helicoidal motion and here exactly a circle in the comoving plane
$x^{\prime }0^{\prime }y^{\prime }$
. The trajectory projected onto the
$x^{\prime }O^{\prime }y^{\prime }$
plane remains a circle to very good accuracy with a relative change in radius less than 1 %. Figure 5 supports this fact. The Lorentz factor is shown in figure 6. In the observer frame it reaches values up to
$\unicode[STIX]{x1D6E4}\approx 10^{5}$
whereas in the electric drift frame it is computed according to
$\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6E4}_{E}\unicode[STIX]{x1D6E4}(1-\unicode[STIX]{x1D737}_{E}\boldsymbol{\cdot }\unicode[STIX]{x1D737})$
. As a check of our algorithm, we plot these Lorentz factors on the same figure 6 and indeed retrieve
$\unicode[STIX]{x1D6FE}=100$
within 8 digits of precision.
Figure 4. Gyromotion of an electron in the electric drift frame with
$\unicode[STIX]{x1D6E4}_{E}\approx 707.107$
and
$\unicode[STIX]{x1D6FE}=100$
in the electric drift frame. Ten rotations are plotted.
Figure 5. Radius of the orbit of an electron in the electric drift frame with
$\unicode[STIX]{x1D6E4}_{E}\approx 707.107$
.
Figure 6. Lorentz factor of an electron in the observer frame and in the drift frame with
$\unicode[STIX]{x1D6E4}_{E}\approx 707.107$
.
The radius of the circular motion suffers from a slight oscillation amplitude increase with time. The electric drift regime in the relativistic limit represents a severe test for any relativistic particle pusher. Indeed, for lower electric drift speed
$v_{E}$
, the circular orbit is better conserved like for instance taking
$\unicode[STIX]{x1D6E4}_{E}\approx 70.71$
(
$E=0.9999$
), see figure 7.
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 Keplerian two body problem of a 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 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 paragraphs. The analytical solutions are described in Uzan & Deruelle (Reference Uzan and Deruelle2014) and for completeness we recall them in the following paragraphs.
Figure 7. Radius of the orbit of an electron in the electric drift frame with
$\unicode[STIX]{x1D6E4}_{E}\approx 70.71$
.
3.4 Central electric force
The two body problem in gravitational physics can be transposed in an equivalent electrostatic problem including relativistic corrections to the 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
$E$
and angular momentum
$L$
. The electric force applied to the orbiting particle is
$$\begin{eqnarray}f=\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}r^{3}}\boldsymbol{r}.\end{eqnarray}$$
The orbital motion stays in a plane that we choose as the
$xOy$
plane. In polar coordinates
$(r,\unicode[STIX]{x1D719})$
the conservation laws for angular momentum and energy are translated into
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}\dot{\unicode[STIX]{x1D719}}=\frac{L}{mr^{2}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}mc^{2}=E-\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}r}. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D6FE}$
from these expressions, the equations of motion in polar coordinates become 
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D719}}=\frac{Lc^{2}}{r^{2}\displaystyle \left(E-\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}r}\right)} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{r}}^{2}=c^{2}\left[1-\frac{L^{2}c^{2}/r^{2}+m^{2}c^{4}}{\displaystyle \left(E-\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}r}\right)^{2}}\right]. & \displaystyle\end{eqnarray}$$
${\dot{r}}=0$
and
$\dot{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D6FA}$
. Equality between centrifugal force and electrostatic attraction induces an angular rotation rate satisfying 
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}^{4}+\left(\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}mr}\right)^{2}\frac{\unicode[STIX]{x1D6FA}^{2}}{r^{2}c^{2}}-\left(\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}mr}\right)^{2}\frac{1}{r^{4}}=0.\end{eqnarray}$$
The solution is
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}=-\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}mr^{2}c\sqrt{2}}\left[-1+\sqrt{1+\frac{4m^{2}c^{4}}{\displaystyle \left(\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}r}\right)^{2}}}\right]^{1/2}.\end{eqnarray}$$
In the general case of an elliptic bound orbit, dividing
${\dot{r}}^{2}$
by
$\dot{\unicode[STIX]{x1D719}}^{2}$
and setting
$u=1/r$
we find a ordinary differential equation
$$\begin{eqnarray}\frac{\text{d}^{2}u}{\text{d}\unicode[STIX]{x1D711}^{2}}+\left[1-\left(\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}Lc}\right)^{2}\right]u=-\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}c^{2}}\frac{E}{L^{2}}.\end{eqnarray}$$
For
$|qQ/4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}Lc|<1$
, the general solution is
$$\begin{eqnarray}\displaystyle & \displaystyle r=\frac{p}{1+e\cos (\unicode[STIX]{x1D6FA}_{p}(\unicode[STIX]{x1D711}-\unicode[STIX]{x1D714}))} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{p}=\sqrt{1-\left(\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}Lc}\right)^{2}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p=\frac{\unicode[STIX]{x1D6FA}_{p}^{2}}{-\displaystyle \frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}c^{2}}\frac{E}{L^{2}}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle e^{2}=\frac{1}{E^{2}}\left[m^{2}c^{4}+\frac{E^{2}-m^{2}c^{4}}{\displaystyle \left(\frac{qQ}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}Lc}\right)^{2}}\right]. & \displaystyle\end{eqnarray}$$
$E$
is split into relativistic kinetic energy
$\unicode[STIX]{x1D6FE}mc^{2}$
and electrostatic potential energy
$U$
. Inspection of figure 9 demonstrates that the total energy is exactly conserved during time evolution.
3.5 Magnetic dipole
In the equatorial plane of a magnetic dipole, the particle gyrates around the origin where the dipole is located. This motion is induced by the magnetic gradient drift in the azimuthal direction. The particle stay within two circles. An example of magnetic drift motion is shown in figure 10 for
$\unicode[STIX]{x1D6FE}=10^{6}$
. The Lorentz factor is conserved as checked in figure 11.
We finish our extensive test of the algorithm by considering two time varying electromagnetic fields represented by linearly and circularly polarized plane waves.
Figure 8. Motion of an electron in the electric field of a fixed proton. The relativistic precession of the orbit is clearly visible. The minimum and maximum radius 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 yellow.
Figure 9. Total energy
$E$
, relativistic kinetic energy
$\unicode[STIX]{x1D6FE}mc^{2}$
and electrostatic potential energy
$U$
of an electron in the electric field of a fixed proton.
$E$
is perfectly conserved, being a constant of motion.
Figure 10. Relativistic drift motion of an electron in a magnetic dipolar field.
Figure 11. Lorentz factor of the drift motion of an electron in a magnetic dipolar field.
3.6 Linearly polarized plane wave
Consider a linearly polarized plane wave propagating along the
$\boldsymbol{e}_{x}$
direction such that the vector potential is
$A^{\unicode[STIX]{x1D6FC}}=(0,0,(E/\unicode[STIX]{x1D714})\cos \unicode[STIX]{x1D709},0)$
. The wave vector is therefore
$(\unicode[STIX]{x1D714}/c,k,0,0)$
from which we deduce the phase
$\unicode[STIX]{x1D709}=\unicode[STIX]{x1D714}t-kx$
. The electromagnetic field is then given by
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{E}=E\sin \unicode[STIX]{x1D709}\boldsymbol{e}_{y} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}=\frac{E}{c}\sin \unicode[STIX]{x1D709}\boldsymbol{e}_{z}. & \displaystyle\end{eqnarray}$$
$u_{0}^{\unicode[STIX]{x1D6FC}}=(c,\mathbf{0})$
. Introducing the strength parameter of the wave by the ratio 
$$\begin{eqnarray}a=\frac{qE}{mc\unicode[STIX]{x1D714}}\end{eqnarray}$$
the four-velocity will have components
$$\begin{eqnarray}\displaystyle & \displaystyle u^{x}=\frac{a^{2}}{2}c(\cos \unicode[STIX]{x1D709}-1)^{2} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u^{y}=-ac(\cos \unicode[STIX]{x1D709}-1) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u^{0}=c+u^{x}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle v^{x}\rangle =\frac{3}{4}\frac{a^{2}c}{1+3a^{2}/4} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle v^{y}\rangle =\frac{ac}{1+3a^{2}/4}. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D709}=0$
, we find 
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}x=\frac{a^{2}c}{8}(6\unicode[STIX]{x1D709}-8\sin \unicode[STIX]{x1D709}+\sin 2\unicode[STIX]{x1D709}) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}y=ac(\unicode[STIX]{x1D709}-\sin \unicode[STIX]{x1D709}) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}ct=c\unicode[STIX]{x1D709}+\unicode[STIX]{x1D714}x. & \displaystyle\end{eqnarray}$$
$x$
axis are shown in figure 12 for a mildly
$a=1$
and an ultra-relativistic
$a=10^{3}$
strength parameter. The mean motion with average velocity
$\langle v_{x}\rangle$
is also shown. The associated Lorentz factor evolution is shown in figure 13. The numerical integration is compared to the analytical solution depicted by coloured symbols.
Figure 12. Motion of an electron in a linearly polarized plane wave for different strength parameters
$a=1,10^{3}$
. Symbols correspond to the analytical solution (3.20).
Figure 13. Lorentz factor of an electron in a linearly polarized plane wave for different strength parameters
$a=1,10^{3}$
. Symbols correspond to the analytical solution (3.18).
3.7 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^{\unicode[STIX]{x1D6FC}}=(0,0,(E/\unicode[STIX]{x1D714})\cos \unicode[STIX]{x1D709},(E/\unicode[STIX]{x1D714})\sin \unicode[STIX]{x1D709})$
and the wave vector
$(\unicode[STIX]{x1D714}/c,k,0,0)$
thus the phase
$\unicode[STIX]{x1D709}=\unicode[STIX]{x1D714}t-kx$
. The electromagnetic field is then given by
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{E}=E(\sin \unicode[STIX]{x1D709}\boldsymbol{e}_{y}-\cos \unicode[STIX]{x1D709}\boldsymbol{e}_{z}) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}=\frac{E}{c}(\sin \unicode[STIX]{x1D709}\boldsymbol{e}_{z}+\cos \unicode[STIX]{x1D709}\boldsymbol{e}_{y}). & \displaystyle\end{eqnarray}$$
$u_{0}^{\unicode[STIX]{x1D6FC}}=(c,\mathbf{0})$
. The time evolution of the components of this four-velocity will be 
$$\begin{eqnarray}\displaystyle & \displaystyle u^{x}=a^{2}c(1-\cos \unicode[STIX]{x1D709})=au^{y} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u^{y}=ac(1-\cos \unicode[STIX]{x1D709}) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u^{z}=-ac\sin \unicode[STIX]{x1D709} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u^{0}=c+u^{x}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle v^{x}\rangle =\frac{a^{2}c}{1+a^{2}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle v^{y}\rangle =\frac{ac}{1+a^{2}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle v^{z}\rangle =0. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D709}=0$
, we find 
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}x=a^{2}c(\unicode[STIX]{x1D709}-\sin \unicode[STIX]{x1D709}) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}y=ac(\unicode[STIX]{x1D709}-\sin \unicode[STIX]{x1D709}) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}z=ac(\cos \unicode[STIX]{x1D709}-1) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}ct=c\unicode[STIX]{x1D709}+\unicode[STIX]{x1D714}x. & \displaystyle\end{eqnarray}$$
$x$
axis are shown in figure 14 for a mildly
$a=1$
and an ultra-relativistic
$a=10^{3}$
strength parameter. The mean motion with average velocity
$\langle v_{x}\rangle$
is also shown. The corresponding evolution of the Lorentz factor is given in figure 15. The numerical integration is compared to the analytical solution depicted by coloured symbols.
Figure 14. Motion of an electron in a circularly polarized plane wave for different strength parameters
$a=1,10^{3}$
. Symbols correspond to the analytical solution (3.24).
Figure 15. Lorentz factor of an electron in a circularly polarized plane wave for different strength parameters
$a=1,10^{3}$
. Symbols correspond to the analytical solution (3.22).
4 Conclusion
We developed an efficient and fully implicit algorithm to solve the Lorentz force equation for ultra-relativistic particles reaching Lorentz factors higher than
$10^{9}$
. Although it requires matrix multiplication, this additional computational cost is completely compensated by the implicit nature of the scheme and its very good behaviour and accuracy for ultra-relativistic motion on large time scales for spatially and temporally varying fields. These features are compulsory for any study of plasmas in high-energy astrophysics where such extreme conditions are easily met such as for instance in neutron stars and especially pulsars for which the ratio between the Larmor frequency and the neutron star spin frequency is huge, spanning more than 10 decades. Several more important tests have been envisaged such as motion in a Coulombian electric field, a kind of relativistic Kepler two-body problem thus in an inhomogeneous electric field. Motion in a linearly or circularly polarized plane wave represents also a challenging test of the algorithm in extreme intense electromagnetic fields. Applications are not restricted to high-energy astrophysics, interaction of high-intensity laser fields with plasmas is another fruitful domain to investigate.
Our implicit algorithm updates the velocity from an evaluation of the electromagnetic field at some point in space–time, chosen as
$t^{\ast },\boldsymbol{r}_{\ast }$
, equation (2.28). Other schemes could be tested such as an interpolation of the Lorentz force in such a way that
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\boldsymbol{r}^{n+1}-\boldsymbol{r}^{n}}{\unicode[STIX]{x0394}t}=\frac{\boldsymbol{v}^{n+1}+\boldsymbol{v}^{n}}{2}=\frac{1}{2}\left(\frac{\boldsymbol{u}^{n+1}}{\unicode[STIX]{x1D6FE}^{n+1}}+\frac{\boldsymbol{u}^{n}}{\unicode[STIX]{x1D6FE}^{n}}\right) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\unicode[STIX]{x0394}t}=\frac{q}{2m}[\boldsymbol{E}(\boldsymbol{r}^{n+1},t^{n+1})+\boldsymbol{v}^{n+1}\wedge \boldsymbol{B}(\boldsymbol{r}^{n+1},t^{n+1})+\boldsymbol{E}(\boldsymbol{r}^{n},t^{n})+\boldsymbol{v}^{n}\wedge \boldsymbol{B}(\boldsymbol{r}^{n},t^{n})]. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Following the discussion presented in the work by Qiang (Reference Qiang2017), it is possible to use higher-order integration schemes to reduce the computational costs to reach a prescribed accuracy. Qiang (Reference Qiang2017) technique can be applied to any time reversible second-order relativistic integrator algorithm.
Ultra-relativistic particles in ultra-strong magnetic fields are prone to damping via radiation reaction. This viscous force needs to be implemented in the equation of motion to catch the full dynamics of relativistic plasmas subject to intense radiation fields. This additional frictional force can be implemented straightforwardly with the fully implicit technique exposed in this paper but is left for future work.
Acknowledgements
I am grateful to the referees for their constructive comments and suggestions. This work has been supported by the French National Research Agency (ANR) through the grant no. ANR-13-JS05-0003-01 (project EMPERE).
