Advanced numerical techniques for time integration of relativistic equations of motion for charged particles

Advanced numerical techniques for solving the relativistic equations of motion for charged particles are provided. A new fourth-order integrator is developed by combining the Taylor series expansion of the numerical angle of relativistic gyration and the fourth-order Runge–Kutta method for integrating the Lorentz factor. The new integrator gives the exact relativistic E-cross-B drift velocity, but has a numerical accuracy much higher than the classic fourth-order Runge–Kutta integrator.


Background
The relativistic equations of motion for charged particles are important for planetary and space science studies, especially for relativistic particle accelerations at collisionless shocks (Umeda and Yamazaki 2006;Nakanotani et al. 2018) and in planetary radiation belts (Katoh 2014;Katoh and Omura 2016;Hiraga and Omura 2020;Omura 2021;Fujiwara et al. 2022).The classic fourth-order Runge-Kutta integrator (RK4) (Runge 1895;Kutta 1901) has been used over many years for solving the relativistic equations of motion for charged particles.Although RK4 has a high numerical accuracy, it does not satisfy any of conservation laws.
As an alternative to RK4, the Boris integrator (Boris 1970) has been used over fifty years in particle-in-cell simulations.The Boris integrator has the second-order accuracy in time.The Boris integrator conserves the kinetic energy during the gyration of charged particles in a magnetic field.This property makes the phase-space trajectories of particles with the Boris integrator more accurate than with the classic RK4 in a long-term time integration (Qin et al. 2013;He et al. 2015).The Boris integrator also provides the exact non-relativistic E × B drift velocity.It has been known, however, that the Boris integrator has a numerical error in the gyration angle per time step.
It is also known that the Boris integrator has a large numerical error in the drift velocity of the relativistic E × B drift (Vay 2008; Zhang et al. 2015; Higuera and  Cary 2017; Ripperda et al. 2018).There have been several attempts for obtaining a more accurate relativistic E × B drift velocity, such as the Vay integrator (Vay 2008) and the Higuera-Cary integrator (Higuera and Cary 2017).These explicit integrators reduce numerical errors in the phase-space trajectories of particles when a particle velocity vector is close to the guiding-center velocity vector only.For a velocity vector far from the guiding-center velocity vector, a numerical error in the phase-space trajectories of particles with the Boris integrator is smaller than that with the Vay or the Higuera-Cary integrators.
Based on the analytic solution to the relativistic E × B drift motion, a new explicit integrator for the relativistic equations of motion for a charged particle has been developed recently (Umeda 2023).The new integrator (hereafter the Umeda integrator) has the second-order accuracy in time and provides the exact relativistic E × B drift velocity.However, the Umeda integrator has a numerical error in the gyration angle per time step as the Boris integrator.A purpose of the present study is to increase the numerical accuracy of the Umeda integrator.

Brief introduction to the Umeda integrator
Let us start with the relativistic equations of motion for a charged particle with a position vector r , a velocity vector v , a mass m, and a charge q, which expresses the accelera- tion of charged particles by electric field E and magnetic field B : where γ is the relativistic Lorentz factor given as: with c being the speed of light.
To simplify Eq.( 1b), a motion of a charged particle is separated into motions in the directions parallel and perpendicular to a magnetic field.The velocity components parallel and perpendicular to a magnetic field are given as follows: Then, Eq.( 1b) is separated into two components parallel and perpendicular to the magnetic field as follows: where u ≡ γ v is a momentum vector and v E is the E × B drift velocity vector with E and E ⊥ being the electric field components paral- lel and perpendicular to the magnetic field, respectively, (1a) Let us take a coordinate system such that B z and E y (i.e., E × B � x ).For the relativistic E × B drift (i.e., c|B| > |E| ) with a constant electromagnetic field, the perpendicular component of the momentum vector in the present coordinate system, i.e., u ⊥ = (u x , u y ) moves along the following elliptical trajectory (Friedman and Semon 2005): where and γ B is a correction factor due to the gradient-B-type drift caused by dγ /dt = 0 , which is also known as a "boosted" Lorentz factor (Ripperda et al. 2018): The boosted Lorentz factor is constant on the elliptical trajectory in Eq.( 5), which is obtained by taking the time derivative of γ B as follows: ( Here, a similar approach to Zenitani and Umeda (2018) for the Boris integrator is used.The elliptical trajectory in the u x − u y space is written by using the elliptical rota- tion matrix as follows: where ω c ≡ q|B|/m is the gyro frequency and τ is the proper time in relativity, This is rewritten by using the half-angle formula of trigonometric functions together with Eq.( 6) as follows: The gyration angle per proper time step �τ is approxi- mated as: Equation ( 8) is then rewritten in the following vector form: where The time-development equation for the momentum vector u is obtained by adding the momentum component parallel to the ambient magnetic field u to Eq.( 9): The proper time step is approximated by: The previous study (Umeda 2023) chose which has been used in the Boris integrator (Boris 1970).Note that the Boris integrator (Boris 1970) is also rewritten with the same form as Eq.( 10): where (10) The Umeda integrator (Umeda 2023) corresponds to the Boris integrator (Boris 1970) if E = 0 .The Umeda inte- grator also approaches to the non-relativistic Boris integrator as c → ∞ (i.e., γ → 1).

Construction of higher-order integrator
For an arbitrary time step t , the momentum vector u obtained with the Umeda integrator (Umeda 2023) is always on the elliptical trajectory, which is the analytic solution to the relativistic E × B drift in a uniform and constant electromagnetic field.That is, the Umeda integrator satisfies dγ B /dt = 0 and gives the exact relativ- istic E × B drift velocity.However, a numerical error in the gyration angle becomes larger as the time step t becomes larger.Below, numerical techniques for reducing numerical errors are given.

Taylor series expansion of tangent function
The Taylor series expansion of the tangent function in Eq.( 8) with Eq.( 11) is given as follows (Kato and Zenitani 2021): The Boris integrator (Boris 1970) and the previous study (Umeda 2023) used the first (i.e., t ) term only.The time-development equation for the momentum vector u is written in the following form by using Eq.( 13): (13) Here, an operator F with arguments of 1/ Ŵ and t is defined for simplicity.The operator F has three forms as described below.
Taylor series expansion: If γ 2 E > 1 , i.e., |v E | < c , the numerical error of the Taylor series expansion from the tangent function becomes larger as the argument q�t 2mγ E Ŵ |B| approaches to π/2 , which appears as a numer- ical error in the gyration angle.The operator F with the Taylor series expansion has an upper-bound in the numerical gyration angle per time step, , γ E disappears from the numera- tor of the operator F with the Taylor series expansion.Hence, the operator F is computationally stable.If γ 2 E < 0 , i.e., |v E | > c ; however, the operator F with the Taylor series expansion becomes unstable numerically for β U < 0 , i.e., q�t Although higher-order terms of the Taylor series relax the numerically unstable condition, a safety factor given in the previous study (Umeda 2023) is necessary to avoid the numerical instability. (14d) Tangent function: If γ 2 E > 1 , i.e., |v E | < c , the tangent function leads to overflow or underflow at q�t 2mγ E Ŵ |B| = π/2 with some compilers.It is recom- mended to use cosine and sine functions instead.If 14) is replaced with −ι tanh(ιθ) , where ι is the imag- inary unit.The operator F with the hyperbolic tangent function is computationally stable, since tanh 2 (ιθ) < 1.
Cosine and sine functions: , |v E | < c , the cosine and sine functions do not lead to overflow or underflow for any argument.14) , |v E | > c , sin θ and cos θ in Eq.( 14) is replaced with −ι sinh(ιθ) and cosh(ιθ) , respectively.The operator F with the hyperbolic cosine and sine func- tions is computationally stable as well.
RK4 integrator is widely used in scientific computing in various fields, which has the fourth-order accuracy in time (Kutta 1901).The momentum vector at time t 1 ≈ t is estimated with the Euler integrator at the first step.The momentum vector at time t 2 ≈ t is re-estimated with the mid-point rule at the second step.The momentum vector at time t 3 ≈ t + �t/2 is estimated with the mid-point rule at the third step.Then, the time integration is performed based on the Simpson integration rule by using 1/γ t− �t 2 , 1/γ t 1 and 1/γ t 2 (at time t), and 1/γ t 3 (at time t + �t/2 ).Finally, the momentum vector at time t + �t/2 is obtained at the fourth step as follows: (15a)

Numerical tests
Numerical errors of the present integrators from the theoretical solution to the relativistic E × B drift in a uni- form and constant electromagnetic field (Friedman and Semon 2005; Umeda 2023) are examined.For evaluating the relativistic gyration angle, the t term only, up to t 3 terms, up to t 5 terms of the Taylor series, and the tan- gent function are used.For the integration of 1/ Ŵ , Euler, (15d)  The order of accuracy for various combinations of the Taylor series terms and the multi-stage integrators is summarized in Table 1.The t term, up to t 3 terms, and up to t 5 terms of the Taylor series have the second-, fourth-, and sixth-order accuracy in time.However, the Euler integrator, i.e., 1/ Ŵ = 1/γ t− �t 2 has the first-order accuracy in time independently of the number of the Taylor series terms.Note that the previous Umeda integrator with the t term only and 1/ Ŵ = 1/γ − has the second-order accuracy in time (Umeda 2023), the choice of which is not so bad.Both of the mid-point rule and the trapezoidal rule have the second-order accuracy in time independently of the number of the Taylor series terms.The accuracy of the trapezoidal rule is slightly higher than that of the mid-point rule (not shown).The Heun3 and RK3 integrators have the third-order accuracy in time with the Taylor series terms higher than t but have the second-order accuracy in time with the t term only.The accuracy of Heun3 is almost the same as that of RK3 (not shown).The RK4 integrator and the Kutta-3/8 rule have the fourth-order accuracy in time with the Taylor series terms higher than t but have the second-order accuracy in time with t term only.The accuracy of the Kutta-3/8 rule is almost the same as that of RK4 (not shown).These results suggest that the order of accuracy is determined by the lowest order of accuracy of either the Taylor series terms or the multi-stage integrators.
Numerical tests are performed with a uniform, constant and time-independent electromagnetic field E = (0, E y , 0) and B = (0, 0, B z ) .The E × B drift veloc- ity v E = E y /B z = 0.8c and the initial velocity vector v 0 = (0.5c, 0, 0) are used.The time step is varied from ω c �t = 2 3 to 2 −9 , where ω c = qB z /m is the non-relativ- istic gyro frequency, The numerical results at ω c t = 24 is compared against the theoretical solution.
Figure 1 shows the relative error η of the momentum vector u = (u x , u y ) and position vector r = (x, y) from the theoretical solution as a function of t with the com- bination of RK4 and different number of Taylor series terms.The circle, triangle, square, and x-mark show the results with t term only, up to t 3 terms, up to t 5 terms, and tangent function, respectively.It is clearly shown that the order of accuracy is second in time with t term only.The order of accuracy is fourth in time with more than t 3 terms.As the number of Taylor series terms increases, the numerical accuracy approaches to that with the tangent function.The numerical error with the Taylor series expansion is enhanced for an extremely large time step, ω c �t γ E Ŵ ≥ π , because of a large numerical error in the gyration angle.In Fig. 1, the maximum argument is ω c �t = 8 with γ E = 5/3 and 1.15 < Ŵ < 2.7 .The computational cost is also measured on a single compute core of the Intel Xeon Gold 6342 processor with the Intel oneAPI compiler ver.2021.5.0.The compile options are "-ipo -ip -O3 -xICELAKE-SERVER." The elapsed time per particle and per time step with t term only, up to  (1) directly), and the present study with the combination of RK4 and the tangent function, respectively t 3 terms, up to t 5 terms, and tangent function is 67.7 ns, 73.3 ns, 76.2 ns, and 94.2 ns, respectively.
Figure 2 shows relative error η of the momentum vec- tor u = (u x , u y ) and position vector r = (x, y) from the theoretical solution as a function of t with the com- bination of the tangent function and different integrators.The circle, triangle, square, and x-mark show the results with Euler, trapezoidal, Heun3, and RK4 integrators, respectively.It is clearly shown that the order of accuracy is first, second, third, and fourth in time with Euler, trapezoidal, Heun3, and RK4 integrators, respectively.
Figure 3 shows the relative error η of the momentum vector u = (u x , u y ) and position vector r = (x, y) from the theoretical solution as a function of t for the previ- ous and present integrators.The circle, triangle, square, and x-mark show the results with Boris (1970), Umeda (2023), fourth-order Runge-Kutta integrators (RK4 to Eq.( 1) directly), and the present study with the combination of RK4 and the tangent function, respectively.It is shown that the numerical error of direct RK4 is much lower than that of both Boris (1970) and Umeda (2023) integrators.The numerical error of the present study is two orders lower than that of direct RK4.The elapsed time per particle and per time step with the Umeda integrator is 58.4 ns.Hence, the numerical cost of the present study is much cheaper than that of the previous study.
A long-term numerical test is also performed with the same electromagnetic field E = (0, E y , 0) and B = (0, 0, B z ) , the E × B drift velocity v E = E y /B z = 0.8c , and the initial velocity vector v 0 = (0.5c, 0, 0) .The time step is fixed to ω c �t = 0.1 , and the numerical test is per- formed up to ω c t = 10 7 .Figure 4 shows the relative error η of the momentum vector u = (u x , u y ) , the position vec- tor r = (x, y) , the constant for the elliptical trajectory C in Eq.( 5), and the boosted Lorentz factor γ B in Eq.( 6) from the theoretical solution as a function of time for the previous and present integrators.
For a short-term numerical test, the numerical errors in the momentum and position vectors with the secondorder (Boris and Umeda), RK4, and present integrators are ∼ 10 −4 , ∼ 10 −8 , and ∼ 10 −10 , respectively, which is consistent with Fig. 3.The numerical error in the constant for the elliptical trajectory C with the Boris integrator varies between 10 −6 and 10 −2 at a certain period ( T ≈ 20/ω c in the present case), which depends on γ E and γ B .This result indicates that the momentum vector moves along an elliptical trajectory with a wrong relativistic E × B drift velocity.The numerical error in the Fig. 4 Relative error η of the momentum vector u = (u x , u y ) , the position vector r = (x, y) , the constant for the elliptical trajectory C in Eq. ( 5), and the boosted Lorentz factor γ B in Eq.( 6) from the theoretical solution as a function of time for the previous and present integrators.(a) 0 ≤ ω c t ≤ 100 and (b) 0 ≤ ω c t ≤ 10 7 .The circle, triangle, square, and x-mark show the results with Boris (1970), Umeda (2023), fourth-order Runge-Kutta integrators (RK4 to Eq.( 1) directly), and the present study with the combination of RK4 and the tangent function, respectively boosted Lorentz factor γ B is related to the numerical error in C, because the center of the elliptical trajectory is given as (γ B γ E v E , 0) .The numerical error in the constant for the elliptical trajectory C with RK4 is ∼ 10 −8 .On the other hand, the numerical errors in the constant for the elliptical trajectory C with the Umeda and present integrators are ∼ 10 −14 .
For a long-term numerical test, the numerical errors in the momentum vectors with the Boris and Umeda integrators vary between 10 −2 and 10 0 because of the accu- mulation of the numerical error in the gyration angle.The numerical error in the momentum vectors with the present integrator is much lower than that with other integrators.Note that oscillations in the numerical errors in Fig. 4b are due to a data sampling rate.The numerical errors in the position vector with the Boris, Umeda, and present integrators are ∼ 10 −4 , ∼ 10 −5 and ∼ 10 −8 , respectively.The numerical error in the constant for the elliptical trajectory C with RK4 increases to 10 −2 , which indicates that the momentum vector deviates from the theoretical trajectory.On the other hand, the numerical error in C with the Boris integrator varies between 10 −6 and 10 −2 , which indicates that the momentum vector is on the elliptical trajectory for a long term.The numerical errors in C and γ B with the Umeda and present integra- tors are ∼ 10 −12 .The result clearly shows that the opera- tor F in Eq.( 14) is well-designed for conserving both C and γ B .
It is noted that the operator F(1/ Ŵ, �t) is Eq.( 14) is replaced with F[1/ Ŵ, E t (r t ), B t (r t ), �t] with the posi- tion vector obtained with the multi-stage numerical integrator, if electromagnetic fields are non-uniform and time-dependent.

Conclusions
Instead of the classic Runge-Kutta (RK4) integrator (Runge 1895;Kutta 1901), the relativistic equations of motion for charged particles have been conventionally solved with the Boris integrator (or the Boris push) (Boris 1970) in particle-in-cell (PIC) plasma simulations.The conventional Boris integrator (Boris 1970), the previous Vay (2008) and Higuera and Cary (2017) integrators conserve the kinetic energy of charged particles during the gyration but have large numerical errors in the relativistic E × B drift velocity.Recently, a new integrator (the Umeda integrator) has been developed, which provides the exact relativistic E × B drift velocity (Umeda 2023).However, the Umeda integrator has the second-order accuracy in time and has a numerical error in the relativistic gyration angle as well.
To reduce the numerical error of the Umeda integrator and to make its order of accuracy higher, two advanced numerical techniques are adopted in the present study.One is the Taylor series expansion of the numerical gyration angle (Kato and Zenitani 2021), and the other is multi-step (multi-stage) Runge-Kutta method (Butcher 1996) for the integration of the Lorentz factor.With the combination of the trigonometric function for the gyration angle and RK4 for the integration of the Lorentz factor, a new fourth-order integrator has been developed.The new integrator conserves both of the kinetic energy of charged particles during the gyration and the relativistic boosted Lorentz factor, but has a numerical accuracy much higher than the classic RK4.The computational cost of the new integrator is much cheaper than the Umeda integrator.Although the order of accuracy is extended up to fourth in the present study, higher-order multi-stage integrators for the integration of the Lorentz factor could make the order of accuracy higher in time.

Euler integrator
The Euler integrator is known as an explicit time-forwarding numerical procedure with the first-order accuracy in time.The momentum vector is updated with 1/ Ŵ = 1/γ t− �t 2 as follows:

Mid-point rule
The numerical integrator based on the mid-point rule has second-order accuracy in time.The momentum vector at time t is estimated with the Euler integrator at the first step.Then, by assuming that 1/ Ŵ = 1/γ t is constant in time from t − �t/2 to t + �t/2 , the momentum vec- tor at time t + �t/2 is obtained at the second step as follows: (A.1a)

Trapezoidal rule
The numerical integrator based on the trapezoidal rule has second-order accuracy in time.The momentum vector at time t 1 ≈ t + �t/2 is estimated with the Euler integrator at the first step.By assuming that 1/γ changes linearly from 1/γ t−�t/2 to 1/γ t 1 , 1/ Ŵ is obtained with the numerical integration based on the trapezoidal rule.Then, the momentum vector at time t + �t/2 is obtained at the second step as follows:

Heun3 integrator
The third-order Heun (Heun3) integrator is known as an explicit time-forwarding numerical procedure with the third-order accuracy in time (Heun 1900).The momentum vector at time t − �t/6 is estimated with the Euler integrator at the first step.Then, the momentum vector at time t + �t/6 is estimated with the mid-point rule at the second step.Finally, the momentum vector at time t + �t/2 is obtained at the third step as follows: (A.2b)

RK3 integrator
The third-order Runge-Kutta (RK3) integrator is known as an explicit numerical procedure with the third-order accuracy in time (Kutta 1901).The momentum vector at time t is estimated with the Euler integrator at the first step.The momentum vector at time t 1 ≈ t + �t/2 is esti- mated at the second step.By assuming that 1/γ changes in a parabolic manner from 1/γ t−�t/2 to 1/γ t 1 via 1/γ t , 1/ Ŵ is obtained with the numerical integration based on the Simpson integration rule.Finally, the momentum vector at time t + �t/2 is obtained at the third step as follows: (A.4a) (A.5a) (A.5b)

Kutta-3/8 rule
There is another version of the fourth-order Runge-Kutta integrator, which is known as the Kutta-3/8 rule (Kutta 1901).The momentum vector at time t − t 6 is estimated with the Euler integrator at the first step.The momentum vector at time t + t 6 is estimated at the second step.The momentum vector at time t 1 ≈ t + �t/2 is estimated at the third step.Then, the time integration is performed based on the numerical integration of a cubic polynomial function passing through four points of 1/γ t− �t 2 , 1/γ t− �t 6 , 1/γ t+ �t 6 , and 1/γ t 1 (at time t + �t/2 ).Finally, the momentum vector at time t + �t/2 is obtained at the fourth step as follows: (A.5c) (A.6a) (A.6b) (A.6c) (A.6d)

Fig. 1
Fig.1Relative error η of the momentum vector u = (u x , u y ) and position vector r = (x, y) from the theoretical solution as a function of t for the present integrators with the combination of RK4 and different number of Taylor series terms.The circle, triangle, square, and x-mark show the results with t term only, up to t 3 terms, up to t 5 terms, and tangent function, respectively

Fig. 2
Fig. 2 Relative error η of the momentum vector u = (u x , u y ) and position vector r = (x, y) from the theoretical solution as a function of t for the present integrators with the combination of the tangent function and different integrators.The circle, triangle, square, and x-mark show the results with Euler, trapezoidal, Heun3, and RK4 integrators, respectively

Fig. 3
Fig.3Relative error η of the momentum vector u = (u x , u y ) and position vector r = (x, y) from the theoretical solution as a function of t for the previous and present integrators.The circle, triangle, square, and x-mark show the results withBoris (1970),Umeda (2023), fourth-order Runge-Kutta integrators (RK4 to Eq.(1) directly), and the present study with the combination of RK4 and the tangent function, respectively

Table 1
Order of accuracy for various combinations of Taylor series terms and multi-stage integrators