On the Boris solver in particle-in-cell simulation

A simple form of the Boris solver in particle-in-cell (PIC) simulation is proposed. It employs an exact solution of the Lorentz-force part, and it is equivalent to the Boris solver with a gyrophase correction. As a favorable property for stable schemes, this form preserves a volume in the phase space. Numerical tests of the Boris solvers are conducted by test-particle simulations and by PIC simulations. The proposed form provides better accuracy than the popular form, while it only requires few additional computation time.


I. INTRODUCTION
Particle-in-cell (PIC) method [1][2][3][4] is one of the most important techniques in modern plasma simulations.Solving Lagrangian motions of many charged particles and the temporal evolution of Eulerian electromagnetic fields, it can simulate various kinetic phenomena, whose length scale is larger than the grid size.
The particle integrator, an algorithm to advance charged particles, is a fundamental element of PIC simulation.Since the particle integrator is used for all the particles, its accuracy, stability, and computational cost has a big impact on those of the entire simulation run.One of the most common integrators is the Boris solver 5 , also known as the Buneman-Boris solver.It solves the particle motion in a leap-frog manner, and the acceleration part is split into several parts, as will be shown later.Owing to its simplicity and reliability, the Boris solver has been used for nearly 50 years.
In addition to the Boris solver, various solvers have been developed.For earlier ones, we refer the readers to classic textbooks 1,2 and references therein.][8][9][10][11][12][13] Vay 6 has developed a solver to carefully deal with the force balance.By splitting the integrator into the explicit first half and the implicit second half, his solver better deals with the E×B drift at a relativistic speed, and therefore it has drawn growing attention in laser physics and in astrophysics.Pétri 7 has proposed a relativistic implicit solver, which iteratively uses a matrix-form of the Vay solver for non-staggerd timesteps.Qiang 8 has presented a fast Runge-Kutta relativistic integrator, ready for the force balance problem.Umeda 9 has proposed a multi-step extension of the Boris integrator.His solver effectively deals with the gyration at a half timestep.
From the theoretical viewpoint, Qin et al. 10 have recently pointed out that the nonrelativistic Boris solver is stable, because it preserves a volume in the phase space every timestep.This property can be examined by a simple Jacobian matrix.Presently, the volume preservation is regarded as a key property for stable solvers.Zhang et al. 11 split the scheme into several substeps to discuss the volume preservation of their solver, which appears to be another expression of the Boris solver.Higuera & Cary 12 have proposed a relativistic volume-preserving solver, which employs Vay's characteristic velocity in the Lorentz-force part of the Boris solver.Ripperda et al. 13 extensively compared selected set of particle solvers.
In this contribution, we propose a simple form of the Boris solver, based on an exact solution for the Lorentzforce part.We further examine the volume preservation of the proposed form of the Boris solver.Then we will present numerical tests of the Boris solvers, followed by discussion and summary.

II. BORIS SOLVER
First, we outline the Boris algorithm. 5It handles the particle motion in the following discrete forms, where the superscripts (n, n + 1, ...) indicate timesteps, u = γv is the spatial part of a 4-vector, γ = [1 − (v/c) 2 ] −1/2 = [1 + (u/c) 2 ] 1/2 is the Lorentz factor, and v is an effective velocity.Other symbols have their standard meanings.The action part is split into the Coulomb force for the first half timestep, the Lorentz force for the entire timestep, and the Coulomb force for the second half: where ε = (q/2m)E n+ 1 2 , u − and u + are two intermediate states.Hereafter we denote the field quantities E n+ 1 2 , B n+ 1 2 at n = n + 1 2 as E, B for brevity.Since Eq. ( 4) is an energy-conserving rotation in the momentum space, the Lorentz factor is set to be constant during the operation, γ − = γ + .The phase angle in the rotation part is The rotation is solved in the following way where b = B/|B| is a unit vector.There are two choices in Eq. (7).One can advance u n to u n+1 by using either of the two equation sets.Boris 5 presented a procedure of Eqs. ( 3), ( 6), (7a), ( 8), (9), and (5) in his original article.The subsequent textbooks 1,2 described a simplified procedure [Eqs.(3), (7b), ( 8), (9), and ( 5)], by replacing Eq. (7a) with Eq. (7b).We call the two procedures the Boris-A solver and the Boris-B solver, respectively.Although we do not describe the detail, the Boris-A solver accurately handles the rotation.As a result of the replacement (Eq.(7a) → Eq. (7b)), the Boris-B solver approximates the rotation, as illustrated in gray in Figure 1.There is always a delay in the gyrophase angle, Thus, the Boris-B solver is an approximate form of the original Boris solver (the Boris-A solver).Despite this, owing to its simplicity and computational cost, the Boris-B solver is widely used.Scientists often indicate the Boris-B solver simply by "the Boris solver."Eqs.(7a) and (7b) differ by a factor of tan( θ 2 )/( θ 2 ).Therefore, the Boris-A solver is sometimes referred to as the Boris solver with a gyrophase correction factor. 1,2< l a t e x i t s h a 1 _ b a s e 6 4 = " N W X i f H 6 D o U O O s 0 R n p A g + 1 3 Y z r j I = " > A A A D X H i c h V L L a t t Q E D 2 2 n E e d 5 l k I h W 5 M T a A r c x U C 2 Y Z k 0 2 V C a z s Q G y P J N 4 5 q v Z C u H Y L Q D 2 T V b f s H I T + S H 8 i i m + 6 7 z q K b L n p 0 L S 9 c 4 / o K S f M 4 Z 2 b u z N i R 5 y Z K i B + l s l F Z W V 1 b f 1 X d e L 2 5 t b 2 z u 9 d K w l H s y K Y T e m F 8 a V u J 9 N x A N p W r P H k Z x d L y b U + 2 7 e F Z 7 m + P Z Z y 4 Y f B Z 3 U W y 6 1 u D w L 1 2 H U v R 1 O q o G 6 m s 3 k 5 d N I Q + t X n B L I Q 6 i n M e 7 p Z q 6 K C P E A 5 G 8 C E R Q F H 2 Y C H h c w U T A h F t X a S 0 x Z R c 7 Z f I U C V 3 R J Q k w q J 1 y O + A 2 l V h D a j n M R P N d p j F 4 x u T W c O B e B Y P 4 k U 8 i U f x S / x Z G C v V M f J a 7 v i 3 J 1 w Z 9 b b v 3 3 7 6 T d Z i n k 1 8 h g P i b Z 2 5 T 8 m n V + F m a b Y p 7 n / 8 W N 9 Q 4 l b f z N e + g K i U v j H t D t l 5 / 7 r a M o 0 x v U e u p a j T n x W d X B R r w E w W 5 U x r e U 9 t P Z 9 s a R V 9 z m o 8 x 8 s n 7 u i Z q y W Z 4 w I z y + / o H Z j M l X z u n P n v h s 0 L r c O G K R r m x V H 9 5 L T Y v n W 8 w 3 t 8 Y A e O c Y K P O E e T N X z B V 3 z D 9 / J P o 2 J s G J s T a L l U c N 5 g 5 h j 7 f w E 9 w r G p < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " N W X i f H 6 D o U O O s 0 R n p A g + 1 3 Y z r j I = " > A A A D X H i c h V L L a t t Q E D 2 2 n E e d 5 l k I h W 5 M T a A r c x U C 2 Y Z k 0 2 V C a z s Q G y P J N 4 5 q v Z C u H Y L Q D 2 T V b f s H I T + S H 8 i i m + 6 7 z q K b L n p 0 L S 9 c 4 / o K S f M 4 Z 2 b u z N i R 5 y Z K i B + l s l F Z W V 1 b f 1 X d e L 2 5 t b 2 z u 9 d K w l H s y K Y T e m F 8 a V u J 9 N x A N p W r P H k Z x d L y b U + 2 7 e F Z 7 m + P Z Z y 4 Y f B Z 3 U W y 6 1 u D w L 1 2 H U v R 1 O q o G 6 m s 3 k 5 d N I Q + t X n B L I Q 6 i n M e 7 p Z q 6 K C P E A 5 G 8 C E R Q F H 2 Y C H h c w U T A h F t X a S 0 x Z R c 7 Z f I U C V 3 R J Q k w q J 1 y O + A 2 l V h D a j n M R P N d p j F 4 x u T W c O B e B Y P 4 k U 8 i U f x S / x Z G C v V M f J a 7 v i 3 J 1 w Z 9 b b v 3 3 7 6 T d Z i n k 1 8 h g P i b Z 2 5 T 8 m n V + F m a b Y p 7 n / 8 W N 9 Q 4 l b f z N e + g K i U v j H t D t l 5 / 7 r a M o 0 x v U e u p a j T n x W d X B R r w E w W 5 U x r e U 9 t P Z 9 s a R V 9 z m o 8 x 8 s n 7 u i Z q y W Z 4 w I z y + / o H Z j M l X z u n P n v h s 0 L r c O G K R r m x V H 9 5 L T Y v n W 8 w 3 t 8 Y A e O c Y K P O E e T N X z B V 3 z D 9 / J P o 2 J s G J s T a L l U c N 5 g 5 h j 7 f w E 9 w r G p < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " N W X i f H 6 D o U O O s 0 R n p A g + 1 3 Y z r j I = " > A A A D X H i c h V L L a t t Q E D 2 2 n E e d 5 l k I h W 5 M T a A r c x U C 2 Y Z k 0 2 V C a z s Q G y P J N 4 5 q v Z C u H Y L Q D 2 T V b f s H I T + S H 8 i i m + 6 7 z q K b L n p 0 L S 9 c 4 / o K S f M 4 Z 2 b u z N i R 5 y Z K i B + l s l F Z W V 1 b f 1 X d e L 2 5 t b 2 z u 9 d K w l H s y K Y T e m F 8 a V u J 9 N x A N p W r P H k Z x d L y b U + 2 7 e F Z 7 m + P Z Z y 4 Y f B Z 3 U W y 6 1 u D w L 1 2 H U v R 1 O q o G 6 m s 3 k 5 d N I Q + t X n B L I Q 6 i n M e 7 p Z q 6 K C P E A 5 G 8 C E R Q F H 2 Y C H h c w U T A h F t X a S 0 x Z R c 7 Z f I U C V 3 R J Q k w q J 1 y O + A 2 l V h D a j n M R P N d p j F 4 x u T W c O B e B Y P 4 k U 8 i U f x S / x Z G C v V M f J a 7 v i 3 J 1 w Z 9 b b v 3 3 7 6 T d Z i n k 1 8 h g P i b Z 2 5 T 8 m n V + F m a b Y p 7 n / 8 W N 9 Q 4 l b f z N e + g K i U v j H t D t l 5 / 7 r a M o 0 x v U e u p a j T n x W d X B R r w E w W 5 U x r e U 9 t P Z 9 s a R V 9 z m o 8 x 8 s n 7 u i Z q y W Z 4 w I z y + / o H Z j M l X z u n P n v h s 0 L r c O G K R r m x V H 9 5 L T Y v n W 8 w 3 t 8 Y A e O c Y K P O E e T N X z B V 3 z D 9 / J P o 2 J s G J s T a L l U c N 5 g 5 h j 7 f w E 9 w r G p < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " N W v c G s 5 r X 9 L g 7 g j < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " v c G s 5 r X 9 L g 7 g j < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " v c G s 5 r X 9 L g 7 g j < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Schematic diagram for the Lorentz-force part of the Boris solvers.Procedures by the Boris-C solver (in black) and by the Boris-B solver (in gray) are illustrated.There is an error in the phase angle (δθ) in the Boris-B solver.

III. EXACT GYRATION SOLVER
To avoid the phase error in the Boris-B solver, one can use the Boris-A solver.Here, we propose another expression.We propose the following rotation procedure based on an analytic solution: This is illustrated in black in Figure 1.Then, one can advance the particle by Eqs. ( 3), ( 6), ( 11), (12), and ( 5) in a time-reversible manner.We call this the Boris-C solver.
In practice, we use a threshold b in |B| 2 = max(B 2 , b ), to avoid dividing by zero.This does not cause a significant error in the B → 0 limit.The Boris-C solver provides second-order accuracy.The splitting into the Coulomb-force part (Eqs.( 3) and ( 5)) and the Lorentz-force part (Eqs.( 11)-( 12)) is equivalent to an operator splitting, also known as the Strang splitting. 14,15The Strang splitting gives second-order accuracy, if each operator maintains more than secondorder accuracy.In this case, since both of the two parts give exact solutions, the resulting scheme should have second-order accuracy.
One can also prove the second-order accuracy in a straightforward manner.In the constant fields, the second-order expansion of the new state u n+1 is where Following the Boris procedure, we find In Eq. ( 17), the Lorentz-force part is expanded for u − , similarly as Eqs.( 13)-( 15), but with E = 0.Although the Boris-C solver accurately solves the third and higherorder terms in Eq. ( 17), we focuses on the terms up to second-order.Importantly, Eqs. ( 16) and ( 18) do not contain second or higher-order terms, O(∆t 2 ), because the solver already gives an exact solution for the Coulomb-force part.From Eq. ( 16), we obtain first-order expansions of the following variables, Substituting Eqs. ( 19)-( 20) into Eqs.( 16)-( 18), one can obtain an expanded form of the Boris-C solver, ).The first and second-order coefficients are identical to those of the Talyor expansion (Eq.( 13)-( 15)) for u = u n .This proves that the Boris-C solver has second-order accuracy.
As a theoretical property, we examine the volume preservation 10 of the Boris-C solver.This discussion holds true for the Boris-A solver as well, because both methods accurately solve the rotation.We consider the temporal evolution of a phase-space volume, by splitting the Boris-C solver to the following substeps. 11We consider half-adjusted positions x n and x n+1 .R stands for a 3-D rotation matrix.
For comparison, we examine the volume preservation for the nonrelativistic fourth-order Runge-Kutta method in Appendix B. The phase-space volume is not preserved in this case, J RK4 = 1.

IV. NUMERICAL TESTS
In order to see the accuracy, stability, and performance of the Boris solvers, we have carried out four numerical tests.First, to check the accuracy, we have carried out test-particle simulations in a uniform electromagnetic field.Physics parameters are set to m = 1, q = 1, c = 1, and u(t = 0) = (1, 0, 0).The time interval, the timestep, and the numerical threshold are set to 0 < t < 12/π, ∆t = π/6, and b = 10 −20 .We consider five configurations, as shown in Table I.The first case corresponds to direct acceleration by the electric field.A weak magnetic field is imposed in the second case, but the electric field is still dominant.The third is a special case of E ⊥ B and |E| = |B|.The fourth case corresponds to the E×B drift.The drift is slightly modulated by the relativistic effect.The last case corresponds to gyration about the magnetic field.In each cases, we evaluate an error in 4vector to a reference value, δu = u − u ref .For cases 1 and 5, we employ analytic solutions as reference values.For cases 2-4, we refer numerical results with a small timestep ∆t = π/240, because we do not know analytic solutions as a simple function of t.The numerical reference values are checked by analytic solutions in other forms. 16,17igure 2(a) shows temporal evolution of errors |δu|/|u| in our test-particle simulations.The dotted lines represent results by the Boris-B solver, whereas the solid lines represent results by the Boris-C solver.The results by the Boris-A solver are not shown, because they are essentially the same as those by the Boris-C solver.In case 1 (in red), the Boris-B and Boris-C solvers use the same parts (Eqs.( 3) and ( 5)) and their results are identical.The two solvers give accurate results, because Eqs. ( 3) and ( 5) give an exact solution for the linear acceleration by E. In cases 2 and 3 (in magenta and in orange), the two solvers give similar results.The curves drop, because |δu| does not grow and because |u| increases in time.In case 4 of the E×B drift (in green), the Boris-C solver drastically improves the results.It reduces the error by two orders-of-magnitude, as indicated by the green arrow.The error by the Boris-B solver grows in time, be-   cause the phase error |δu| ≈ |u|δθ = δθ accumulates in time.On the other hand, the error by the Boris-C solver remains small.The repeated drop corresponds to the gyroperiod.Since the Boris-C solver exactly solves the phase, the solid curve repeatedly show the same pattern, and it does not grow further.In case 5 of gyration (in blue), the error by the Boris-B solver linearly grows in time (∝ t) because of the accumulation of the phase error.In contrast, the Boris-C solver gives accurate results, as we have expected.
Figure 2(b) presents maximum errors (|δu|/|u|) max during 0 < t < 12π as a function of the timestep, ∆t = π/60, π/20, π/6, and π/2.The results demonstrate the second-order accuracy of the solvers.As evident in cases 1-3, the Boris-C solver is as good as the Boris-B solver when |E| |B|.As the electric field dominates, the errors decrease, and then they give accurate results in case 1 (red line).When |B| |E|, the two solvers give different results.In cases 4 and 5, the Boris-B solver remains moderately good (the green and blue dotted lines).From Eq. ( 10) and using γ = √ 2, one can estimate the curve in case 5, |δu|/|u| ≈ (12π/∆t)δθ = (π/ √ 8)∆t 2 , in agreement with the blue dotted line.In contrast, the Boris-C solver produces drastically smaller errors than the Boris-B solver, even though it maintains the secondorder accuracy.As the magnetic field dominates, the errors decrease, and then the solver gives accurate results in case 5 (blue solid line).In the bottom part in Figure 2(b), the errors in the two exact cases probably come from a machine error.For example, for case 5, if an error of O(10 −15 ) in double precision floating numbers is accumulated every timestep, the total error would be |δu|/|u| 10 −15 × (12π/∆t) ≈ 10 −13.5 ∆t −1 , in consistent with the blue line.The error in case 1 (red line) is even smaller, because of larger |u|.
As a second numerical test, we evaluate the long-term stability of the particle solvers.We have run a code in the following field, as was done in Ref. 10, The initial conditions are u(t = 0) = (0.1, 0.0), x(t = 0) = (0.9, 0, 0), m = c = 1, and ∆t = π/10.Figure 3 shows the trajectories at (a) an initial stage and at (b) a late stage by the Boris-C solver (in blue) and by the fourth order Runge-Kutta solver (in gray).From the beginning, the particle undergoes a fast small-scale gyration and a slow large-scale rotation due to the ∇B drift and the E×B drift.The large-scale drifts keep the particle in the same domain, and then we inspect the evolution of the trajectory.After a long time (300th turn; t ≈ 2 × 10 5 ), as evident in Fig. 3(b), the Runge-Kutta solver dissipates the small-scale motion and the relevant kinetic energy.This is because the Runge-Kutta solver is not volume-preserving, at least in the uniform fields in the nonrelativistic regime (Appendix B), and quite likely so in generic relativistic cases.In contrast, the Boris-C solver is free from the numerical damping, similarly as the Boris-B solver (not shown).This numerical experiment demonstrates that the Boris-C solver does preserve the phase-space volume for the small-scale gyration over a long time.
(a) 1st turn (b) 300th turn x < l a t e x i t s h a 1 _ b a s e 6 4 = " A B s 7 V I a u n x b w y 4 0 m 8 U 7 t Z P P s k a z z O I j 7 F F v K E y N y k 5 9 E a 4 m Z h t i P u P H 6 g b S t y q m z n K 5 x L V p 6 9 L u 0 l 2 0 r + 6 s g V I a u n x b w y 4 0 m 8 U 7 t Z P P s k a z z O I j 7 F F v K E y N y k 5 9 E a 4 m Z h t i P u P H 6 g b S t y q m z n K 5 x L V p 6 9 L u 0 l 2 0 r + 6 s g V I a u n x b w y 4 0 m 8 U 7 t Z P P s k a z z O I j 7 F F v K E y N y k 5 9 E a 4 m Z h t i P u P H 6 g b S t y q m z n K 5 x L V p 6 9 L u 0 l 2 0 r + 6 s g V I a u n x b w y 4 0 m 8 U 7 t Z P P s k a z z O I j 7 F F v K E y N y k 5 9 E a 4 m Z h t i P u P H 6 g b S t y q m z n K 5 x L V p 6 9 L u 0 l 2 0 r + 6 s g As a third numerical test, to compare computational costs, we have run the case 5 in Table I until t = 2 × 10 7 π (1.2 × 10 8 steps), by only advancing the particle velocity.Average elapse times are proportional to 1.46 : 1 : 1.24 for the Boris-A, B, and C solvers on our PC (Intel processor) and 1.46 : 1 : 1.26 on FX100 computer (SPARC processor) at Japan Aerospace Exploration Agency (JAXA).While performance depends on implementation, compilers, and CPU architectures, the Boris-C solver runs faster than the Boris-A solver, and it runs only 25% slower than the Boris-B solver.
As a final test, we have carried out PIC simulations of a physics problem on magnetic reconnection.The settings are documented in Ref. 18, and therefore we do not repeat them here.Comparing average elapse times, we have found that the Boris-C solver is slower than the Boris-B solver by 1.2% on XC50 supercomputer (Intel processor) at National Astronomical Observatory of Japan and by 3.3% on FX100 computer (SPARC processor) at JAXA.

V. DISCUSSION AND SUMMARY
As already discussed, the Boris-C solver is a combination of the two exact solvers for the Coulomb-force part and for the Lorentz-force part, while the Boris-B solver is based on the exact Coulomb-force solver and a second-order solver for the Lorentz-force part.From the viewpoint of operator-splitting, in both cases, the combined solver maintains second-order accuracy, because each part has second-order accuracy.This is evident in Figure 2(b).Since the Boris-C solver better deals with the Lorentz-force part, the amplitude of the second-order error is much smaller than in the Boris-B solver, in particular in the |B| |E| cases.The Boris-B solver gives a second-order error in phase (Eq.( 10)).Other secondorder solvers, such as the Higuera-Cary solver and the Umeda solver, could suppress the total error, by better solving the phase. 9,12,13In this line, the Boris-C and A solvers provide the best possible results, because they give no phase error in phase.
We have formally proved that the Boris-C solver preserves a volume in the phase space during the temporal evolution, and then we have confirmed that it preserves the small-scale gyration after a long-time computation.The volume preservation of the popular Boris-B solver has been extensively studied, [10][11][12] however, it has never been evaluated in the Boris-A or C solvers.Since the Boris-C solver is simple, the proof is given straightforwardly.This gives further confidence to the reliablity of the Boris-C and A solvers.
Many scientists prefer the Boris-B solver to the Boris-A solver because of the low computation cost.In fact, according to the numerical test, the Boris-A solver is 46% more expensive than the Boris-B solver.Since the Boris-C solver employs simple expressions, it is more favorable than the Boris-A solver.The Boris-C solver gives 25% in test-particle simulations.Obviously the conventional solvers (the Boris-A and Boris-B solvers) tried to avoid trigonometric functions which were expensive at that time.Presently, these functions are not so expensive as they used to be, and therefore the Boris-C solver runs adequately fast.In PIC simulations, the Boris-C solver runs slower only by 1-3% in PIC simulations.Typically, the particle integrator and the electric current calculator account for the most of the computation time, and they consume equal amount of time.Then largely due to memory-access performance, the particle integrator does not run at full speed.If we assume that it runs at 50% efficiency, the particle integrator is responsible for 25% of the computation time.Then, the 25% slow-down in the particle integrator should result in 6% slow-down in the entire run.We have observed 1-3% in our PIC simulations, probably because our code runs less efficiently, but the results are reasonable.The computational cost can be easily compensated by employing a larger ∆t.
From the viewpoint of the stability, we usually keep the timestep small, ω max ∆t < 2, where ω max is the maximum frequency to resolve. 2 Nevertheless, in PIC simulation, the magnitude of the magnetic field may instantly approach or exceed the qB γmc ∆t = 2 criteria, and so it is useful to check the stability for a large ∆t limit.The Boris-B solver delays the gyrophase (Eq.( 10)), and then the angle never exceed π, i.e., θ < π.The particle eventually moves back-and-forth. 19The Boris-A and Boris-C solvers have no such limitation.They simply allow particles to gyrate more than π, although it leads to an opposite gyration.Meanwhile, the Boris-A solver needs some care for θ ≈ π, where Eq. (7a) is undefined or approaches ±∞.Therefore the Boris-B and Boris-C solvers are safer choices among the three.
1][22][23] In fact, the symplectic algorithms [24][25][26] are favored in many applications, owing to their long-term accuracy.However, the symplectic schemes are often implicit and computationally expensive, while the number of explicit symplectic solvers for relativistic charged particles in arbitrary electromagnetic fields is limited. 22,23In addition, it was reported that volume-preserving solvers are sometimes more accurate than symplectic solvers. 27Considering these issues, it will take some time before symplectic solvers will be popular in relativistic PIC simulations.Meantime, it is important to improve the second-order Boris solvers, which have an advantage in computational efficiency and which are proven to be reliable.
In summary, we have proposed a simple form of the Boris solver, based on the exact solution for the Lorentzforce part.It is equivalent to the Boris solver with a gyrophase correction.It has a favorable property of preserving the phase-space volume, and therefore it appears to be stable.The proposed form gives much more accurate results than the popular form (the Boris-B solver), while it only requires few additional computation time.We hope that the proposed form will be a good alternative to the conventional Boris solvers.
5 w b e i h t y 5 W 1 y I p s e e 4 H U n c M W 5 4 Z 1 / u p / 6 w r g 9 D y 3 5 w b e i h t y 5 W 1 y I p s e e 4 H U n c M W 5 4 Z 1 / u p / 6 w r g 9 D y 3 5 w b e i h t y 5 W 1 y I p s e e 4 H U n c M W 5 4 Z 1 / u p / 6 w r g 9 D y 3 5 w b e i h t y 5 W 1 y I p s e e 4 H U n c M W 5 4 Z 1 / u p / 6 w r g 9 D y 3

9 e 8 M
W j c O X P U / 0 v D f b 3 d 2 9 5 v a t 4 x m e 4 w U 7 s I N d 7 O M A f e 7 h E t f 4 g q 9 O 2 x H O j v N q R l 1 r N Z o n u L G c v b 8 h z c P h < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 0 0 y X n P b n 0 O y B 6 k T B J 6 U 5 y y s O m S V m K + m 0 s W W + C 4 + i 1 / i m 7 g S P 8 W f p b l K m 6 P e y z n / w U y r 8 u HG p 6 d v f 1 O 1 X B e Q X 2 G L / M B W H t F K G D V 4 v 7 L a n P c / v b Y n V P h g T 5 b Y W E p W y d i U e E h 1 3 b + B R e Y 5 5 u e o v R J d x q u m k 8 t y j V l J 0 q 6 s V / c 0 s P O p V u 5 i x F l N F 3 T 1 x E M 7 c 7 O i s m 4 4 N / W + v Q O z u V L P O + f 9 e 8 M W j c O X P U / 0 v D f b 3 d 2 9 5 v a t 4 x m e 4 w U 7 s I N d 7 O M A f e 7 h E t f 4 g q 9 O 2 x H O j v N q R l 1 r N Z o n u L G c v b 8 h z c P h < /l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 0 0 y X n P b n 0 O y B 6 k T B J 6 U FIG. 2. (a) Time evolution of numerical errors |δu|/|u| in test-particle simulations.The dotted lines indicate the results by the Boris-B solver, whereas the solid lines indicate the results by the Boris-C solver.The timestep is fixed to ∆t = π/6.(b) Maximum errors, (|δu|/|u|)max during 0 < t < 12π as a function of ∆t.
9 1 m B e J G l V S n 5 A t / o R 4 o / 4 A y 7 6 C a 5 d u B H 0 5 D Z d a K m 9 I c k 8 z p mZ O z O G b 1 t h J E Q v k 8 1 N T E 5 N 5 2 d m 5 + Y X F p e W C y u n o d c O T F k 1 P d s L a o Y e S t t y Z T W y I l v W / E D q j m H L M + P 6 I P G f d W Q Q W p 5 7 E t 3 5 s u H o L d e 6 t E w 9 o q l y e 7 F c E m W h T n F Y 0 F K h h P Q c e Y V M E e d o w o O J N h x I u I g o 2 9 A R 8 q l D g 4 B P W w N d 2 g J K l v J L x J g l t 0 2 U J E K n 9 Z r f F r V 6 a n W p J z F D x T a Z x e Y b k F n E p n g T z + J D v I o X 8 S 6 + R s b q q h h J L X f 8 G 3 2 u 9 C + W 7 t e O P 8 k a z T O I j 7 F J v K E y N y k 5 9 E a 4 G p t t g P u P H 6 g b S t y o m z n K 5 x L V p a 9 D u 0 l 2 0 r + G s g x i D O 6 R a F 2 U 6 I / T T o 6 K 1 W I m n X K s t K S n h p p P P L a K J m f V G e I l E z f V z K M x m Y M U 8 5 t / r n a g P 1 f y u X P a 3 w 0 b F k 6 3 y 5 o o a 5 W d 0 t 5 + u n 1 5 r G M D W + z A L v Z w i C N U W Y P E A x 7 x l O 1 l v 3 N T u X w f m s 2 k n F X 8 O r n C D 8 gg s H w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " A B s 7 8 p J f Q 0 9 1 m B e J G l V S n 5 A t / o R 4 o / 4 A y 7 6 C a 5 d u B H 0 5 D Z d a K m 9 I c k 8 z p mZ O z O G b 1 t h J E Q v k 8 1 N T E 5 N 5 2 d m 5 + Y X F p e W C y u n o d c O T F k 1 P d s L a o Y e S t t y Z T W y I l v W / E D q j m H L M + P 6 I P G f d W Q Q W p 5 7 E t 3 5 s u H o L d e 6 t E w 9 o q l y e 7 F c E m W h T n F Y 0 F K h h P Q c e Y V M E e d o w o O J N h x I u I g o 2 9 A R 8 q l D g 4 B P W w N d 2 g J K l v J L x J g l t 0 2 U J E K n 9 Z r f F r V 6 a n W p J z F D x T a Z x e Y b k F n E p n g T z + J D v I o X 8 S 6 + R s b q q h h J L X f 8 G 3 2 u 9 C + W 7 t e O P 8 k a z T O I j 7 F J v K E y N y k 5 9 E a 4 G p t t g P u P H 6 g b S t y o m z n K 5 x L V p a 9 D u 0 l 2 0 r + G s g x i D O 6 R a F 2 U 6 I / T T o 6 K 1 W I m n X K s t K S n h p p P P L a K J m f V G e I l E z f V z K M x m Y M U 8 5 t / r n a g P 1 f y u X P a 3 w 0 b F k 6 3 y 5 o o a 5 W d 0 t 5 + u n 1 5 r G M D W + z A L v Z w i C N U W Y P E A x 7 x l O 1 l v 3 N T u X w f m s 2 k n F X 8 O r n C D 8 gg s H w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " A B s 7 8 p J f Q 0

2 u 9 CFIG. 3 .
FIG. 3. Particle trajectories by the Boris-C solver during (a) the initial stage and (b) the late stage.A trajectory by the RK4 solver is overplotted in gray in (b) the late stage.