Multiple Boris integrators for particle-in-cell simulation

TL;DR

This paper introduces multiple Boris integrators for particle-in-cell simulations, significantly reducing errors and enabling larger timesteps through a novel combination of Boris procedures and Chebyshev polynomial techniques.

Contribution

The paper develops multiple Boris integrators that enhance accuracy and stability in PIC simulations by combining Boris steps and employing Chebyshev polynomials for one-step formulation.

Findings

Errors reduced by a factor of n^2

Allow larger timesteps without loss of stability

Demonstrated improved performance in numerical tests

Abstract

We construct Boris-type schemes for integrating the motion of charged particles in particle-in-cell (PIC) simulation. The new solvers virtually combine the 2-step Boris procedure arbitrary n times in the Lorentz-force part, and therefore we call them the multiple Boris solvers. Using Chebyshev polynomials, a one-step form of the new solvers is provided. The new solvers give n^2 times smaller errors, allow larger timesteps, and have a long-term stability. We present numerical tests of the new solvers, in comparison with other particle integrators.

Tables

Table 1. Table 1: Field settings for test-particle simulations.

#	$𝑩$	$𝑬$
1	(0,0,0)	(1,0,0)	direct acceleration by E
2	(0,0,0.1)	(1,0,0)	E-dominated
3	(0,0,1)	(1,0,0)	$\left|𝑬\right|=c\left|𝑩\right|$
4	(0,0,1)	(0.1,0,0)	E$\times$B drift
5	(0,0,1)	(0,0,0)	gyration about B

Equations

$$\frac{{\mathbf{x}}^{t+\Delta t}-{\mathbf{x}}^t}{\Delta t}$$

$$m\frac{{\mathbf{u}}^{t+\frac{\Delta t}{2}}-{\mathbf{u}}^{t-\frac{\Delta t}{2}}}{\Delta t}$$

$${\mathbf{u}}^{-}$$

$$\mathbf{t}$$

$${\mathbf{u}}^{\prime }$$

$${\mathbf{u}}^{+}$$

$${\mathbf{u}}^{t+\frac{\Delta t}{2}}$$

$$\alpha \equiv \arctan t$$

$$\sin \alpha =\frac{t}{\sqrt{1+t^2}},\cos \alpha =\frac{1}{\sqrt{1+t^2}},\tan \alpha =t$$

$$\sin 2\alpha =\frac{2t}{1+t^2},\cos 2\alpha =\frac{1-t^2}{1+t^2}$$

$${\mathbf{u}}^{+}$$

$$={\mathbf{u}}^{-}\cos 2\alpha +({\mathbf{u}}^{-}\times \hat{\mathbf{b}})\sin 2\alpha +(1-\cos 2\alpha ){\mathbf{u}}_{\parallel }^{-}$$

$$=({\mathbf{u}}^{-}-{\mathbf{u}}_{\parallel }^{-})\cos 2\alpha +({\mathbf{u}}^{-}\times \hat{\mathbf{b}})\sin 2\alpha +{\mathbf{u}}_{\parallel }^{-},$$

$$\theta\approx 2\alpha=2\arctan\frac{\theta}{2}=\theta\bigg{(}1-\frac{1}{3}\Big{(}\frac{\theta}{2}\Big{)}^{2}+\frac{1}{5}\Big{(}\frac{\theta}{2}\Big{)}^{4}-\cdots\bigg{)}.$$

$$\mathbf{t}\equiv \frac{\theta }{4}\hat{\mathbf{b}}=\frac{q\Delta t}{4m{\gamma }^{-}}\mathbf{B},\alpha \equiv \arctan t$$

$${\mathbf{u}}^{+}$$

$$=({\mathbf{u}}^{-}-{\mathbf{u}}_{\parallel }^{-})\cos 4\alpha +({\mathbf{u}}^{-}\times \hat{\mathbf{b}})\sin 4\alpha +{\mathbf{u}}_{\parallel }^{-}$$

$$\theta\approx 4\alpha=4\arctan\frac{\theta}{4}=\theta\bigg{(}1-\frac{1}{3}\Big{(}\frac{\theta}{4}\Big{)}^{2}+\frac{1}{5}\big{(}\frac{\theta}{4}\Big{)}^{4}-\cdots\bigg{)}.$$

$$\mathbf{t}\equiv \frac{\theta }{2n}\hat{\mathbf{b}}=\frac{q\Delta t}{2nm{\gamma }^{-}}\mathbf{B},\alpha \equiv \arctan t.$$

$${\mathbf{u}}_n^{+}=({\mathbf{u}}^{-}-{\mathbf{u}}_{\parallel }^{-})\cos (2n\alpha )+({\mathbf{u}}^{-}\times \hat{\mathbf{b}})\sin (2n\alpha )+{\mathbf{u}}_{\parallel }^{-}$$

$${\mathbf{u}}_n^{+}$$

$$=c_{n1}\mathbf{u}+c_{n2}(\mathbf{u}\times \mathbf{t})+c_{n3}(\mathbf{u}\cdot \mathbf{t})\mathbf{t}$$

$$\cos nx=T_n(\cos x),\sin nx=\sin x{U}_{n-1}(\cos x)$$

$$p\equiv \cos 2\alpha =\frac{1-t^2}{1+t^2},1+p=\frac{\sin 2\alpha }{t}=\frac{2}{1+t^2}.$$

$$c_{n1}$$

$$c_{n2}$$

$$c_{n3}$$

$$T_{mn}(x)=T_m(T_n(x)),U_{mn-1}(x)=U_{m-1}(T_n(x))U_{n-1}(x)$$

$$c_{n1}$$

$$c_{n2}$$

$$c_{n3}$$

$$c_{11}=\frac{1-t^2}{1+t^2},c_{12}=\frac{2}{1+t^2},c_{13}=\frac{2}{1+t^2},$$

$$c_{21}=\frac{1-6t^2+t^4}{(1+t^2{)}^2},c_{22}=\frac{4(1-t^2)}{(1+t^2{)}^2},c_{23}=\frac{8}{(1+t^2{)}^2},$$

$$c_{31}=\frac{1-15t^2+15t^4-t^6}{(1+t^2{)}^3},$$

$$c_{32}=\frac{2(t^2-3)(3t^2-1)}{(1+t^2{)}^3},c_{33}=\frac{2(t^2-3{)}^2}{(1+t^2{)}^3},$$

$$\displaystyle c_{41}=\frac{\big{(}(1-t^{2})^{2}-24t^{2}\big{)}(1-t^{2})^{2}+16t^{4}}{(1+t^{2})^{4}}~{}~{}$$

$$=\frac{1-28t^2+70t^4-28t^6+t^8}{(1+t^2{)}^4},$$

$$c_{42}=\frac{8(1-t^2)(1-6t^2+t^4)}{(1+t^2{)}^4},c_{43}=\frac{32(1-t^2{)}^2}{(1+t^2{)}^4}$$

$$\theta\approx 2n\alpha=2n\arctan\frac{\theta}{2n}=\theta~{}\bigg{(}1-\frac{\theta^{2}}{12n^{2}}+\frac{\theta^{4}}{80n^{4}}-\cdots\bigg{)}.$$

Full text

††journal: Computer Physics Communications

1 Introduction

Particle-in-cell (PIC) simulation [1, 2] has been extensively used to study complex kinetic phenomena in plasma systems. The Boris solver [3] (or the Boris pusher) is a de fact standard method to advance charged particles in PIC simulation. Although its algorithm is simple, it is robust and provides moderately good accuracy. For these reasons, it has been used about a half century.

Since mid 2010s, there is a renewed interest to develop new particle integrators. In recent PIC simulations, one often solves the system evolution much longer than before. In such a case, we want to solve the particle motion as accurately as possible, because small errors could be accumulated over many timesteps. The Boris solver is proven to be stable [4, 5], however, it does produce a second-order error. For this reason, there is a growing demand for higher-accuracy extensions to the Boris solver.

Umeda [6, 7] has proposed multistep Boris solvers. Repeating a Boris-type procedure multiple steps in the Lorentz-force part, the author managed to set an effective timestep $2^{N-2}$ times smaller. Since the Boris solver is a second-order scheme, $N$-step Boris solvers provide $4^{N-2}$ times higher accuracy for gyration. Zenitani & Umeda [8] have proposed to use a rotation formula in the Lorentz-force part. The solver virtually provides no error in the phase angle in gyration. Because of operator-splitting, their scheme offers second-order accuracy.

In addition to higher numerical accuracy, it is also important to improve computation cost. A fast solver helps us to reduce the total computation time of PIC simulation. The cost depends on many other factors, but the timestep $\Delta t$ is certainly one of them. If the solver allows us to use a larger $\Delta t$, we can further reduce the computation time. Thus it is important to evaluate the constraint on $\Delta t$ for particle solvers.

In this article, we propose new Boris-type integrators to solve the charged-particle motion. The new solvers use a Boris-type procedure multiple times in the Lorentz-force part. We call them “multiple Boris solvers.” They provide higher accuracy and allow larger timesteps than the original Boris solver. In Section 2, we briefly review the classical 2-step Boris solver and discuss its graphical meanings. In Section 3, we introduce our multiple Boris solvers. Using Chebyshev polynomials, we derive coefficients for $n$-tuple Boris solvers, where $n$ is an arbitrary integer. Then we discuss errors, limitations to the timesteps, and the stability of the phase-space volume. In Section 4, we present benchmark results. The new solvers and other particle integrators are briefly compared. Section 5 contains discussion and summary.

2 Boris solver

We describe a popular form of the Boris algorithm [1, 2]. We refer to it as the classical Boris solver. The algorithm solves the particle motion in a time-staggered manner,

[TABLE]

The superscripts ($t,t+\Delta t$, …) indicate time, $\bm{u}=\gamma\bm{v}$ is the spatial part of a 4-vector, and $\gamma=[1-(v/c)^{2}]^{-1/2}=[1+(u/c)^{2}]^{1/2}$ is the Lorentz factor. Other symbols have their standard meanings. Hereafter we denote $\bm{B}^{t}$ as $\bm{B}$ for brevity.

The acceleration part (Eq. (2)) is split into the Coulomb-force part for the first half:

[TABLE]

the middle Lorentz-force part:

[TABLE]

and the Coulomb-force part for the second half:

[TABLE]

In these equations, $\bm{\varepsilon}=({q}/{2m})\bm{E}^{t}$, $\omega_{c}=({qB}/{m}\gamma^{-})$ is the gyro frequency, $\theta=2t=\omega_{c}\Delta t$ is the gyration angle during $\Delta t$, $\bm{\hat{b}}={\bm{B}}/|B|$ is a unit vector in the $\bm{B}$ direction, $\bm{u}^{-}$ and $\bm{u}^{+}$ are two intermediate states.

We focus on the Lorentz-force part in the middle (Eqs. (4a)–(4c)). This 2-step procedure is illustrated in Figure 1(a). The coefficient $2/(1+t^{2})$ is provided such that the particle energy remains constant in the Lorentz-force part, $\gamma^{-}=\gamma^{+}=\gamma^{t}$. A base angle $\alpha$ is given by,

[TABLE]

We will use the following relations for $\alpha$ in subsequent discussions.

[TABLE]

From Eqs. (4b), (4c), and (8), we obtain

[TABLE]

where the subscript $\parallel$ indicates the component along $\bm{B}$. These forms are known as Euler–Rodrigues rotation formula. The equations correspond to a rotation about the magnetic field by an approximate angle of $2\alpha$, instead of the true angle $\theta$. As indicated in Figure 1(a), the half acceleration in the tangent direction ($\bm{u}^{-}\times\bm{t}$; Eq. (4b)) corresponds to the half approximate angle $\alpha$ by the Boris solver. The true rotation angle $\theta$ contains a second-order error of $\delta\theta/\theta\approx-\frac{1}{12}\theta^{2}$.

[TABLE]

3 Multiple Boris solver

Umeda [6] has proposed a novel way to reduce the error ($\delta\theta/\theta$). He employed a half timestep $\Delta t\rightarrow\frac{1}{2}\Delta t$ in the Lorentz-force part (Eqs. (4a)–(4c)) and then repeated the Boris procedure twice. Importantly, the timestep for the position (Eq. (1)) and the Coulomb-force parts (Eqs. (3) and (5)) are unchanged. We call this method the double Boris solver, and it is illustrated in Figure 1(b). Instead of Eq. (4a), we use

[TABLE]

The base angle $\alpha$ is redefined, as indicated in Figure 1(b). Then, we repeat the 2-step Boris procedure (Eqs. (4b) and (4c)) twice. We obtain a half state $\bm{u}^{\prime\prime}$ by the first Boris procedure, which is given by Eq. (11) with Eq. (13). Applying $\cdot\bm{\hat{b}}$ to Eq. (11), one can confirm that the parallel component is preserved, $\bm{u}^{\prime\prime}_{\parallel}=(\bm{u}^{\prime\prime}\cdot\bm{\hat{b}})\bm{\hat{b}}=\bm{u}^{-}_{\parallel}$. Substituting $\bm{u}^{\prime\prime}$ to Eq. (11) again, we obtain the next state $\bm{u}^{+}$,

[TABLE]

The second Boris operation advances the phase angle by another $2\alpha$, and the entire procedure gives gyration by $4\alpha$. In Ref. [6], the author proposed a 3-step procedure to repeat the Boris operation twice. Now we obtain a smaller second-order error in gyration, $\delta\theta/\theta\approx-\frac{1}{48}\theta^{2}$.

[TABLE]

Let us extend this idea. In order to better approximate the gyro motion, we propose to repeat the Boris procedure $n$ times, where $n$ is an arbitrary positive integer. We call the resulting schemes “multiple Boris solvers.” Here below, we describe its procedure in the Lorentz-force part.

We define the tangent vector $\bm{t}$ and the base angle $\alpha$ as follows.

[TABLE]

Note that $\bm{t}$ is divided by $n$. From the discussion on the double Boris solver, repeating the Boris operation $n$ times, one can obtain a rotation angle of $2n\alpha$. The next state $\bm{u}^{+}_{n}$ should be given by

[TABLE]

The subscript $n$ indicates quantities for the $n$-tuple solver. For brevity, we hereafter omit the minus sign from $\bm{u}^{-}$. Since we desire a simulation-friendly expression with $\bm{t}$, we modify Eq. (18) to

[TABLE]

where $c_{n1}$, $c_{n2}$, and $c_{n3}$ are coefficients.

Next, we construct generic expressions for these coefficients. We use Chebyshev polynomials of the first kind $T_{n}(x):T_{0}(x)=1,~{}T_{1}(x)=x,~{}T_{2}(x)=2x^{2}-1,\dots$ and the second kind $U_{n}(x):U_{0}(x)=1,~{}U_{1}(x)=2x,~{}U_{2}(x)=4x^{2}-1,\dots$. They satisfy useful relations:

[TABLE]

We also employ

[TABLE]

With help from Eqs. (7), (8), (21), and (22), we express the coefficients in Eq. (20) in the following way,

[TABLE]

Here, $k$ is a positive integer. We derive Eq. (28) for $n=2k+1$ ($n=3,5,\dots$) in A. Practically, we compute $p$ first, and then calculate Eqs. (23)-(28) by using $p$. Some higher-order polynomials can be calculated via nesting relations [9]:

[TABLE]

We can also express the coefficients with $t$. From the polynomials,

[TABLE]

we obtain explicit forms, for example,

[TABLE]

[TABLE]

It is interesting to see that the $c_{n1}$, $c_{n2}$, and $c_{n3}$ share the same denominator, $(1+t^{2})^{n}$. One can confirm this by inspecting Eqs. (30)–(35).

As a result, the multiple Boris solver advances the velocity by using Eqs. (3), (17), (22)–(28), (20), and (5). One can also use the coefficients in Eqs. (36)–(39) instead of Eqs. (22)–(28).

The Lorentz-force part of the multiple Boris solvers provides a second-order error in phase angle. From Eq. (17), we obtain

[TABLE]

Thus, the multiple Boris solvers provide $n^{2}$ times smaller errors than the single Boris solver. Eq. (40) approaches the exact angle for $n\rightarrow\infty$.

In Figure 2(a), we plot relative accuracy in phase angle, $\theta^{-1}\cdot 2n\arctan({\theta}/{2n})$, for several $n$. In the classical Boris solver, which is equivalent to the $n=1$ case, we typically set $\theta=\omega_{c}\Delta t<0.2$–$0.3$ to keep an error within 1%. Considering ${\theta^{2}}/({12n^{2}})\lesssim 0.01$, we obtain

[TABLE]

The multiple Boris solvers allow $n$-times larger timestep than the single Boris solver.

Next, we evaluate an effective Larmor radius. We consider a simple gyration about $\bm{B}$ in the case of $\bm{E}=0$. From Eq. (1), we obtain

[TABLE]

where $v_{R}$ is the velocity. The two points $\bm{x}^{t+\Delta t}$ and $\bm{x}^{t}$ are located on the same circle with a Larmor radius of $r_{L}^{\prime}$. The two points should be two endpoints of an arc with a central angle $2n\alpha$. Then we expect

[TABLE]

From Eqs. (42) and (43), we obtain

[TABLE]

where $r_{L}=|v_{R}|/\omega_{c}$ is a physical Larmor radius. For $n=1$, we obtain a familiar result of $r^{\prime}_{L}=r_{L}[1+(\theta/2)^{2}]^{1/2}$ from Eq. (7). In general cases, Eq. (44) can be expanded to

[TABLE]

Figure 2(b) shows the relative length of Larmor radius, $r^{\prime}_{L}/r_{L}={\theta}/[2\sin(n\alpha)]$. In the square brackets in Eq. (45), the first term is insensitive to $n$, while the second depends on $n$. Thus, higher $n$ does not allow $n$-times larger timestep. To keep an error within $\approx 1\%$, we need to limit the timestep to

[TABLE]

For $n\rightarrow\infty$, we find $\omega_{c}\Delta t<0.5$. If we set the threshold to $\epsilon\%$, both Eqs. (41) and (46) can be relaxed by a factor of $\sqrt{\epsilon}$. For example, Eq. (46) gives $\omega_{c}\Delta t<1$ for $4\%$.

The multiple Boris solvers have two nice properties for stable solvers. The first one is the time symmetry. As evident in Eqs. (3), (18), and (5), the multiple Boris solvers is symmetric in time. In such a case, it is unlikely that the solvers numerically damp or amplify a motion in one time direction. In Eq. (20), procedures in the reverse-time direction are obtained by reversing the sign of the second coefficient, i.e., $c_{-n1}=c_{n1},c_{-n2}=-c_{n2}$, and $c_{-n3}=c_{n3}$. It is obvious that $c_{01}=1,c_{02}=0$, and $c_{03}=0$ for $n=0$. As a result, we have the coefficients not only for positive integers, but also for any integer $n$.

The second is the volume preservation [4, 5], which is related to the long-term stability. We consider the evolution of a volume in the phase-space from $t-\frac{\Delta t}{2}$ to ${t+\frac{\Delta t}{2}}$. Since the Lorentz-force part of the multiple Boris solver is equivalent to $n$-times repetition of the Boris-type element, we express the equivalent numerical procedures into the following steps [5] of $S_{1},S_{2},\cdots S_{5}$, where $S_{k}$ indicates the $k$-th step. The third step, the Lorentz-force part of the solver, is further divided into $n$ substeps. Boris$(\bm{u},\Delta t)$ indicates the Boris-type procedure (Eqs. (4a)–(4c)) on the vector $\bm{u}$ for a timestep $\Delta t$.

[TABLE]

We check the Jacobian $J_{k}$ for the $k$-th step, $S_{k}$.

[TABLE]

For the first two and the last two steps ($S_{1}$, $S_{2}$, $S_{4}$, and $S_{5}$), we obtain

[TABLE]

In addition, Zhang et al. [5] proved that the Lorentz-force part of the relativistic Boris solver (Boris$(\bm{u},{\Delta t})$) preserves a phase-space volume,111The authors claimed that they developed a new volume-preserving algorithm, but the algorithm appears to be the standard relativistic Boris solver. and the proof holds true for an arbitrary timestep $\Delta t$. Then, for the substeps in the third step ($S_{31},\cdots S_{3m},\cdots S_{3n}$), each Boris-type element (Boris$(\bm{u},\frac{\Delta t}{n})$) preserves the phase-space volume, i.e., $J_{31}=\cdots=J_{3m}=\cdots=J_{3n}=1$. As a result, the Jacobian throughout the entire procedure satisfies $J_{\rm all}=J_{1}J_{2}(J_{31}\cdots J_{3m}\cdots J_{3n})J_{4}J_{5}=1$. This indicates that the entire procedure preserves the phase-space volume. The multiple Boris solvers are volume-preserving.

4 Numerical tests

To evaluate the new solvers, we have conducted four numerical tests. For comparison, the following ten particle integrators are benchmarked:

Classical Boris solver (Eqs. (4a)–(4c)), 2. 2.

Single Boris solver ($n=1$; Eqs. (20), (36)) 3. 3.

Double Boris solver ($n=2$; Eqs. (20), (37)) 4. 4.

Quadruple Boris solver ($n=4$; Eqs. (20), (39)). 5. 5.

Octuple Boris solver ($n=8$; Eqs. (20), (22)–(28)). 6. 6.

$16$x Boris solver ($n=16$; Eqs. (20), (22)–(28)). 7. 7.

$32$x Boris solver ($n=32$; Eqs. (20), (22)–(28)). 8. 8.

Umeda solver ($N=3$) [6] 9. 9.

Umeda solver ($N=4$) [7] 10. 10.

Exact-gyration solver [8]

The first one is the classical Boris solver [1, 2]. The next six are the multiple Boris solvers with $n=1$, $n=2$, $n=4$, $n=8$, $n=16$, and $n=32$. The coefficients are given by Eqs. (36)–(39) or Eqs. (22)–(28). We refer to them as the single, double, quadruple, octuple, $16$x and $32$x Boris solvers, respectively. The next two are the Umeda solvers with $N=3$ steps [6] and $N=4$ steps [7]. We follow procedures in Eqs. (8)–(10) in Ref. [7]. As they contain $2^{N-2}$ Boris elements, they are equivalent to the double and quadruple Boris solvers. The last one is Zenitani & Umeda [8]’s solver, which directly calculates Euler’s rotation formula in the Lorentz-force part:

[TABLE]

By definition, this solver does not produce an error during the gyration phase. We refer to it as the exact-gyration solver. Since it uses trigonometric functions, its computation cost is usually higher than those of polynomial-based solvers. Note that all these solvers share the same Coulomb-force procedures (Eqs: (3) and (5)).

In the first test, we have run test-particle simulations in a static magnetic field $\bm{B}=(0,0,1)$ and $\bm{E}=(0,0,0)$. Parameters are set to $m=1$, $q=1$, and $c=1$. The initial condition is $\bm{u}(t=0)=(1,0,0)$ and the timestep is set to $\Delta t=\pi/60$, $\pi/20$, $\pi/6$, or $\pi/2$. We have calculated errors in $\bm{u}$ from the analytic solutions. We have run the simulations until $t=12\pi$, and then we have evaluated the maximum errors in a 4-vector, $\{|\delta\bm{u}|/|\bm{u}|\}_{\rm max}$. Since $\omega_{c}=qB/(m\gamma)=1/\sqrt{2}$, the interval $12\pi$ corresponds to $\approx 4$ rotations. This test is a subset of numerical tests in our recent study [8]. Here, almost all of the solvers are combinations of the exact solver for the Coulomb-force parts (Eqs. (3) and (5)) and the second-order solver for the Lorentz-force part. These solvers only differ in the Lorentz-force part.

Figure 3 shows the results as a function of $\Delta t$. Errors by the classical Boris solver and the Umeda solvers are not presented, because the errors are identical to those by the single, double, and quadruple Boris solvers. One can see that the multiple Boris schemes give second-order errors with respect to $\Delta t$. The errors also scale like $\propto n^{-2}$. This is because the error simply accumulates in the interval of our interest. From Eq. (40), the errors can be estimated to $\delta\theta\approx(\omega_{c}12\pi)\frac{1}{12}(\omega_{c}\Delta t/n)^{2}=(\pi/\sqrt{2})(\omega_{c}\Delta t/n)^{2}$. The results are in excellent agreement with the prediction. Meanwhile, the exact-gyration solver only gives small errors of $\mathcal{O}(10^{-13})\sim\mathcal{O}(10^{-14})$.

In the second test, we have measured elapse times of test-particle simulations. Since algorithms’ performance relies on CPU architectures, we have run the code on Intel and SPARC processors. We have compiled the program on an Intel Core i7 processor with Intel Fortran compiler 2019. For each solver, we have run the program for $1.2\times 10^{7}$ timesteps. We have measured the elapse time six times for every cases. The black bars in Figure 4 show the results, in normalized units. We have also benchmarked the program on a SPARC processor with FUJITSU Fortran compiler on FX100 supercomputer at the Japan Aerospace Exploration Agency (JAXA). The elapsed times are averaged over three runs, because the results on the SPARC processor have less variations than those on the Intel processor. The results are presented in the gray bars in Figure 4. To emphasize the difference, the leftmost value is set to 0.5. Note that we show total elapse time of test-particle simulations. The programs contain not only the Lorentz-force part, but also the Coulomb-force parts (Eqs. (3) and (5)) and the position part (Eq. (1)). Thus the results are not proportional to, but highly affected by the computation time of the Lorentz-force solvers.

As far as we have tested, the new solvers are adequately fast. On the Intel processor, they are slower than the classical Boris solver, but faster than the Umeda solvers and the exact-gyration solver except for the $n=32$ case. The multiple Boris solvers become slower as $n$ increases. On the SPARC processor, the new solvers are more promising than on the Intel processor. Some of them ($n=1$–$4$) are 50% faster than the exact-gyration solver and even faster than the classical Boris solver. It is noteworthy that the double Boris solver ($n=2$) is significantly faster than the Umeda solver ($N=3$), although the two solvers are equivalent. Similarly, the quadruple Boris solver ($n=4$) is faster than the Umeda solver ($N=4$). Surprisingly, even the 32x Boris solver ($n=32$) looks as good as the Umeda solver ($N=3$), even though it gives $256$ times higher accuracy.

In the third test, we have run test-particle simulations with the single ($n=1$) Boris solver, the quadruple ($n=4$) Boris solver, and the exact-gyration solver in various configurations. We consider five configurations, as shown in Table 1. In cases 1 and 2, the electric field E is dominant: $|\bm{E}|\gg c|\bm{B}|$. In case 3, the electric and magnetic fields are equal: $|{E}|=c|{B}|$. In cases 4 and 5, the magnetic field B is dominant: $|\bm{E}|\ll c|\bm{B}|$. In particular, case 5 stands for the gyration about B, as studied in the first test. Other parameters, timesteps, and the initial conditions are similar to those in the first test. We evaluate maximum values of the normalized error $(|\delta\bm{u}|/|{u}|)_{\rm max}$ during the time interval of $0<t<12\pi$. Here, $\delta\bm{u}$ is deviations from analytic solutions (cases 1 and 5) or numerical reference values (cases 2–4). We obtain the numerical reference values by using the exact-gyration solver with $\Delta t=\pi/240$. These configurations and procedures are the same as in Ref. [8].

Figure 5 shows the errors as a function of the timestep. The dotted, solid, and dashdotted lines indicate the results by the single ($n=1$) Boris solver, the quadruple ($n=4$) solver, and the exact-gyration solver, respectively. In case 1 (in red) by the three solvers and case 5 by the exact-gyration solver, one can essentially see no error. All the other curves suggest the second-order accuracy. In cases 2 and 3, the three solvers return similar results. Meanwhile, in cases 4 and 5, their results are different. The quadruple Boris solver provides 10.9–16.5 times smaller errors than the single Boris solver, as emphasized by the solid arrows in green and blue. The improvement is asymptotic to $16$ for $\omega_{c}\Delta t\ll 1$ in case 5. The exact-gyration solver provides even better results, as indicated by the dashed arrows. The quadruple Boris solver provides similar errors in cases 4 and 5. So does the single Boris solver. It appears that the accuracy in the Lorentz-force solver (case 5) limits the overall accuracy for $|\bm{E}|\ll c|\bm{B}|$ (case 4).

In the fourth test, we study the long-term stability of the particle motion in a non-uniform electromagnetic field. As done in Refs. [4, 8], we have run test-particles in the following field,

[TABLE]

The conditions are set to $\bm{u}(t=0)=(0.1,0.0)$, $\bm{x}(t=0)=(0.9,0,0)$, $m=c=1$, and $\Delta t=\pi/10$. We compare the single ($n=1$) Boris solver, the quadruple ($n=4$) Boris solver, the exact-gyration solver, and the fourth order Runge–Kutta solver. The blue and gray lines in Figure 6 show the trajectories by the quadruple Boris solver and by the Runge–Kutta solver, at (a) an initial stage and at (b) a late stage. The time history of the total energy is presented in Figure 6(c). It combines the particle kinetic energy and the electric potential energy at the midpoint positions ($\bm{x}(t+\frac{\Delta t}{2}$)). The values are normalized by the initial values.

One can see typical trajectories in Figure 6(a). It consists of a fast small-scale gyration and a slow large-scale rotation. The large-scale rotation is due to the $\nabla B$ drift and the E$\times{\bf B}$ drift. In the case of the quadruple Boris solver, the particle keeps gyrating even after a long time (300th turn; $t\approx 2\times 10^{5}$) in Fig. 6(b). The particle energy is excellently conserved. We only recognize a relative error of $\mathcal{O}(10^{-4})$. The trajectories (not shown), the time histories, and the error amplitudes (not shown) by the single ($n=1$) Boris solver and the exact-gyration solver are similar to those of the quadruple ($n=4$) Boris solver. In contrast, as evident in the gray line in Fig. 6(b), the particle gyroradius becomes smaller and smaller in the case of the Runge–Kutta solver. The particle numerically lost 30% of its total energy (Fig. 6(c)).

These results are attributed to the volume-preserving property of the numerical schemes. It is known that the Runge–Kutta solver is not volume preserving. The solver is incapable of preserving the energy of the small-scale oscillation. On the other hand, as discussed in Section 3, the multiple Boris solvers are volume-preserving. The other two solvers are also known to be volume-preserving [4, 5, 8]. These three schemes are capable of preserving the energy of the small-scale oscillation, even in a nonuniform field.

5 Discussion and Summary

We have proposed the multiple Boris integrators that virtually use the Boris procedure $n$ times in the Lorentz-force part. The schemes provide second-order but $n^{2}$-times smaller errors (Eq. (40)) than the classical Boris solver. The coefficients are given by either the generic forms of $p$ (Eqs. (23)–(28)), the generic forms of $t$ (Eqs. (30)–(35)), or the explicit forms (Eqs. (36)–(39)). In our experience, the generic forms with $p$ are useful for large $n$. For small $n$ ($n\leq 4$), we prefer the explicit form, because we can further optimize the code. In our implementation for $n=1\dots 4$, we have calculated the numerators and denominators in Eqs. (36)–(39) separately, because the explicit forms share the same denominator, $(1+t^{2})^{n}$. As $\bm{t}$ contains $1/\gamma$ (Eq. (17)), we multiply the numerators and denominators by $\gamma^{2n}$. By doing so, we reduce the number of division operations to save the total computation cost.

The new solvers are developed in the same spirit as in the multistep Boris solvers [6, 7], which considered $2^{N-2}$-times smaller timesteps in the Lorentz-force part. In the multiple Boris solvers, we virtually divide the rotation angle into $n$ Boris units. This $n$ is arbitrary and not limited to $2^{N-2}$. Thus one can choose the best $n$ for one’s application. In addition, we calculate the vector equation (Eq. (20)) only once, because we prepare the coefficients in advance. The coefficients are rational fractions, and so it is easy to calculate them. Consequently, the program runs efficiently, as evident in the benchmark results.

We have also evaluated numerical accuracy of several solvers in general cases of $\bm{B}\neq 0$ and $\bm{E}\neq 0$. Figure 5 demonstrates the second-order accuracy in most cases, because the second-order Lorentz-force solvers are combined with the Coulomb-force procedures (Eqs. (3) and (5)) in an operator-splitting manner. In this study, since all the Lorentz-force solvers are combined with the same Coulomb-force procedures, the numerical accuracy for $|\bm{E}|\ll c|\bm{B}|$ is essentially controlled by the Lorentz-force solvers, even though the resulting schemes provide the same second-order accuracy. This is evident in cases 4 and 5 in Figure 5 — the quadruple ($n=4$) Boris solver provides $\approx 16$ times better accuracy, and the exact-gyration solver, which corresponds to the multiple Boris solver in the $n\rightarrow\infty$ limit, provides even better accuracy. In general, the multiple Boris solvers have higher accuracy than the classic Boris solver, and then they bring better accuracy in the $|\bm{E}|\ll c|\bm{B}|$ cases. Importantly, we often consider $|\bm{E}|\ll c|\bm{B}|$ in practical applications in PIC simulation.

The new solvers may further reduce the computation cost, by allowing larger timesteps. In PIC simulation, the timestep $\Delta t$ is limited by the plasma frequency $\omega_{p}$ and the gyrofrequency $\omega_{c}$. The ratio of $\omega_{p}$ to $\omega_{c}$ depends on applications, and both frequencies are limited to $\omega\Delta t\lesssim\mathcal{O}(0.1)$ due to the numerical accuracy and $\omega\Delta t<2$ due to the numerical stability [1]. If the gyrofrequency limits the timestep, the multiple Boris solvers can relax the limitations for the accuracy. As evaluated in Eq. (41), the $n$-tuple Boris solver allows $n$-times larger timestep, because it repeats the classical Boris solver with a smaller timestep $(\Delta t/n)$ in the Lorentz-force part. Actually, in the other parts of PIC simulation, the timesteps are unchanged. For example, we update the position (Eq. (1)) for $\Delta t$, not for $\Delta t/n$. This fact makes the timestep condition more complicated. The new condition is given by Eq. (46). This is more restrictive than Eq. (41).

In the nonuniform electromagnetic field, the multiple Boris solvers have a favorable property of the long-term stability. Since they inherit a volume-preservation property from the Boris-type element, the new solvers are capable of preserving a small-scale oscillation, as they should be. As far as we have tested, they conserve the total energy very well. These facts give us further confidence to use the new solvers in PIC simulations for practical applications.

Among the new solvers, considering the accuracy, the timestep, and computation costs, we recommend the quadruple ($n=4$) Boris solver with $\omega_{c}\Delta t\lesssim 0.5$ to the readers. One can develop higher-order solvers ($n\gg 32$) to obtain even better accuracy. In such a case, we recommend the exact-gyration solver [8] to the readers. As shown in Figure 4, the 32x Boris solver ($n=32$) is already more expensive than the exact-gyration solver on Intel processors. It is likely that higher-order solvers are as expensive as the exact-gyration solver on SPARC processors for $n\gg 32$. In any case, since the extra cost is limited, the exact-gyration solver will be another choice for $n>10$. Finally, we note that our recommendation is partly based on the benchmarks on our specific systems. Performance may be different on other systems such as GPU-based systems and FGPA systems. In general, since it is impossible to predict future computing platforms, we can only prepare as many particle solvers as possible. Our solvers will be the latest additions.

In summary, we have proposed the multiple Boris solvers that virtually combine multiple Boris units in the Lorentz-force part. We have derived one-step expressions with polynomial-based coefficients for arbitrary $n$ divisions. The new solvers provide second-order, $n^{2}$-times smaller errors than the classical Boris solver. They also allow larger timesteps. Their performance is promising. They are stable over a long time. We hope that the proposed solvers will be useful in studies with PIC simulations.

Acknowledgements

This work was partially carried out on facilities of the JSS2 system at Japan Aerospace Exploration Agency (JAXA). This work was supported by Grant-in-Aid for Scientific Research (C) 17K05673, (S) 17H06140, and (B) 17H02877 from the Japan Society for the Promotion of Science (JSPS).

Appendix A Third coefficients for $n=2k+1$

We recall the following relation of trigonometric functions.

[TABLE]

We substitute $(2k+1)\alpha$ to $A$, where $k$ is a positive integer, $k\geq 1$. Using Chebyshev polynomials, we obtain

[TABLE]

With help from Eqs. (7), (8), and (22), we obtain

[TABLE]