A high-order Boris integrator

# A high-order Boris integrator

Mathias Winkel Robert Speck Daniel Ruprecht Institute of Computational Science, University of Lugano, Switzerland. Jülich Supercomputing Centre, Forschungszentrum Jülich, Germany.
###### Abstract

This work introduces the high-order Boris-SDC method for integrating the equations of motion for electrically charged particles in an electric and magnetic field. Boris-SDC relies on a combination of the Boris-integrator with spectral deferred corrections (SDC). SDC can be considered as preconditioned Picard iteration to compute the stages of a collocation method. In this interpretation, inverting the preconditioner corresponds to a sweep with a low-order method. In Boris-SDC, the Boris method, a second-order Lorentz force integrator based on velocity-Verlet, is used as a sweeper/preconditioner. The presented method provides a generic way to extend the classical Boris integrator, which is widely used in essentially all particle-based plasma physics simulations involving magnetic fields, to a high-order method. Stability, convergence order and conservation properties of the method are demonstrated for different simulation setups. Boris-SDC reproduces the expected high order of convergence for a single particle and for the center-of-mass of a particle cloud in a Penning trap and shows good long-term energy stability.

###### keywords:
Boris integrator, time integration, magnetic field, high-order, spectral deferred corrections (SDC), collocation method
journal: Journal of Computational Physics\biboptions

sort&compress

## 1 Introduction

Often when modeling phenomena in plasma physics, for example particle dynamics in fusion vessels or particle accelerators, an externally applied magnetic field is vital to confine the particles in the physical device Stacey2005 (); Stacey2010 (). In many cases, such as instabilities Goldston1995a () and high-intensity laser plasma interaction Gibbon2005b (), the magnetic field even governs the microscopic evolution and drives the phenomena to be studied. Movement of electrically charged particles in an electric and magnetic field is described by the following equations of motion

 dvdt (1a) dxdt =v (1b)

with the particle position , its velocity , the magnetic field , electric field and the charge to mass ratio . In (1a), the well-known Lorentz force depends on and so that a discretizing (1) with a standard velocity-Verlet scheme Verlet1967 (); Birdsall1985 () reads

 xn+1 =xn+Δt(vn+Δt2f(xn,vn)) (2a) vn+1 =vn+Δt2(f(xn,vn)+f(xn+1,vn+1)). (2b)

with an implicit velocity update step (2b). The Boris integration method Boris1970 (); Birdsall1985 () provides a clever way to evaluate (2b) without having to actually solve an implicit system. It has thus become a de-facto standard for the numerical solution of (2) and allows to cheaply integrate the particle trajectory in the presence of electric and magnetic fields.

Being based on the velocity-Verlet scheme, the Boris approach is a second-order method Birdsall1985 (). Whether it is also symplectic is controversial: In Webb2014 () it is claimed that it is while QinEtAl2013 () claim that it is not, but nevertheless shows excellent long term energy stability due to being phase-space volume preserving. Furthermore, it only requires a single evaluation of the right-hand side per time step, making it a cheap numerical method in terms of computational cost Patacchini2009 (). For these reasons, the Boris method is widely used in many Particle-In-Cell-codes (see e. g. Verboncoeur2005 ()), grid-free methods (e. g. Gibbon2010a ()) and Monte-Carlo simulations (e. g. Kirschner2000 ()). Several explicit alternatives to the Boris method have been proposed, compare Patacchini2009 () and references therein. All of them are second-order accurate and apparently no higher-order methods based on the original Boris approach exist. Especially for applications such as trajectory integration in particle accelerators Toggweiler2014 (), space-weather studies Lapenta2012 (), high-intensity laser-plasma interaction Gibbon2005b (), and fusion vessel simulations Bowers2009 (); Kirschner2000 (), where high precision has to be maintained over long physical simulation times, these are desirable, though. In addition, the current development of high-performance computing systems towards high floating point operation rates at stagnating memory data transfer speed favors the use of higher-order methods in essentially all fields of computer simulation Dongarra2014 (). Furthermore, the ability of tuning the order of an integration algorithm adds a new dimension to its parameter space that allows for balancing precision versus runtime.

A number of other high-precision or higher-order methods for (1) have been developed. Examples are methods that use a spatial coordinate instead of time as the independent variable which showed better performance than a fourth-order Runge-Kutta method in beam propagation simulations Stoltz2002 (), a quasi-symplectic Trotter-factorization based scheme that builds upon an explicit-implicit mixture of leap-frog, Verlet, exponential differencing and Boris rotations with a sixth-order rotation angle approximation Bowers2009 () or a Taylor series-based explicit approach with an up to sixth-order replacement for the Boris method using a complex differential operator for the Maxwell fields Quandt2010 (); Quandt2007 (). Essentially, none of these methods are easily tunable for arbitrary order but are formulated for a very specific case. Only the latter, Taylor series-based approach offers this feature but requires a complicated set of appropriate differential operators to be constructed for every order.

In this work, we introduce the high-order Boris-SDC integration method for (1), which is a combination of the classical Boris integrator with the spectral deferred corrections method (SDC). The resulting Boris-SDC method retains much of the simplicity of the Boris integrator (in terms of implementation, alas not derivation) while allowing to easily generate a method of essentially arbitrary order. Based on classical defect correction, SDC has originally been introduced for first-order ODEs as an iterative approach for the generic construction of high-order integration schemes using a low-order base propagator (the “sweeper”) such as implicit or explicit Euler for the correction “sweeps” DuttEtAl2000 (). Several modifications and extensions exist e. g. semi-implicit SDC Minion2003 (), GMRES-accelerated SDC HuangEtAl2006 (), inexact SDC SpeckEtAl2014_DDM2013 (), multi-level SDC SpeckEtAl2014_BIT (), SDC based on DIRK methods Weiser2014 () or the parallel full approximation scheme in space and time (PFASST), a parallel-in-time integrator which exploits the iterative structure of SDC EmmettMinion2012 (). Recently, SDC has been formulated for second-order problems with the standard Verlet integrator as sweeper Minion2014_2ndOrderSDC (). Here, we combine this particular approach with the classical Boris integrator and extend it to a velocity-dependent force of the form (1a).

The derivation of Boris-SDC relies on the interpretation of SDC as a preconditioned Picard iteration for the solution of a collocation problem, see e. g. HuangEtAl2006 (); MinionEtAl2015 (). Collocation methods are based on the integral formulation of an ODE, the approximation of the exact trajectory over a time step by a polynomial and evaluation of the integrals by quadrature. They are a special class of implicit Runge-Kutta methods and, depending on the chosen quadrature nodes, have a number of attractive properties, particularly symplecticness, see e. g. hairer_nonstiff (); hairer_stiff (). The disadvantage of collocation methods is that they require the solution of a very large, possibly nonlinear system of equations to compute the stages. Picard iteration can be used to solve this system, but often requires a too small time step for convergence. SDC can be considered as a preconditioned Picard iteration, where inverting the preconditioner corresponds to ”sweeping” through the quadrature nodes with a low-order method. If sufficiently many sweeps are performed, the advantageous properties of the underlying collocation method are recovered. For e. g. a first-order method such as the implicit Euler as sweeper, SDC formally gains one order per sweep ShuEtAl2007 (), so fixing the number of iterations allows to easily generate a scheme of higher order, up to the order provided by the underlying quadrature. Here, we describe how the classical Boris integrator can be used as a preconditioner to derive an iterative solver for a collocation approximation of (1).

This paper is organized as follows: Section 2 describes collocation methods, briefly discusses their properties and introduces the required notation. In Section 3, we start with spectral deferred corrections based on the velocity-Verlet scheme as base integrator in matrix from. The matrix formulation itself is derived in A. Concentrating on this rather formal notation, these parts are sufficiently general to also be utilized for force expressions other than the Lorentz force in (1a). In the second part of Section 3 we then specialize the formalism to the case of the Lorentz force as the ODE’s right-hand side and derive ready-to-implement expressions for the Boris-SDC method, specifically tailored for problems of the form (1). Section 4 illustrates the properties of Boris-SDC by numerical examples and compares Boris-SDC to the classical Boris integrator. Finally, Section 5 gives a summary and an outlook on possible future directions of research.

## 2 Collocation formulation

In this section, we briefly describe collocation methods and introduce the notation required for the spectral deferred correction approach of Section 3. Note that the notation below is loosely based on the discussion of collocation methods and SDC with a velocity-Verlet integrator as base method for second-order problems in Minion2014_2ndOrderSDC ().

Rewriting equations (1) in Picard formulation for an arbitrary interval with starting value , , , we obtain111We use the following notation here: Vectors and matrices in normal font refer to scalar values at a single node in time. Bold-faced variables indicate aggregation over spatial variables (e. g. particles in more than one dimensions and/or multiple particles). Vectors with capitals are used to denote aggregation over all intermediate steps. Matrices are always denoted with capital letters, slightly abusing our own convention. All matrices in this context, however, refer to aggregated quantities anyway.

 v(t) =v0+∫ttnf(x(s),v(s)) ds, (3a) x(t) =x0+∫ttnv(s) ds. (3b)

Collocation methods are based on the introduction of intermediate nodes

 tn≤τ1<…<τM≤tn+1, M≥1 (4)

and approximating the integrals using quadrature. Details can be found e. g. in (hairer_nonstiff, , II.7). To allow for a more convenient notation below, we set . The quadrature weights are collected in a matrix with entries

 ¯qm,j:=∫τmtnℓj(s) ds,m,j=1,…,M (5)

where , are Lagrange polynomials. Then, in order to account for the initial values , , we further define

 Q:=(000¯Q)∈R(M+1)×(M+1). (6)

For quadrature nodes with , e. g. Gauss-Lobatto nodes, the second row of is zero as well, because .

Let be vectors with approximate values for and at the nodes and with the vector containing the corresponding right-hand side values. Then, the matrix provides approximations of the integral in (3a) with , that is

 QF(X,V)=(Q⊗Id)F(X,V)≈(∫τmtnf(x(s),v(s)) ds)m=0,1,…,M (7)

with being the identity matrix and the standard Kronecker product. Similarly, the term approximates the integrals over . The discrete version of  (3) is then given by the collocation formulation

 V =V0+QF(X,V) (8a) X =X0+QV=X0+QV0+QQF(X,V). (8b)

for

 X0:=(x0,x0,…,x0)T,V0:=(v0,v0,…,v0)T∈R(M+1)d. (9)

We note that .

### 2.1 Collocation system

In order to combine both equations into a single, closed expression, we first note that depends on both and , or, more precisely, each component depends on the tuple . The ordering of Equations (8) (first all the , then all the ), however, is not compatible with the sorting in , where the first entry is depends on , the second on etc. Thus, we need to resort the matrix formulation (8) so that the degrees-of-freedom are ordered as

 U=(u0,u1,…,uM)T:=(x0,v0,…,xm,vm,…,xM,vM)T∈R(M+1)2d. (10)

To this end, we introduce permutation operators , and with

 Ix=(10),Iv=(01)andIxv=(0100) (11)

so that for . The permutation operators and redistribute the entries of and to match the sorting of the degrees of freedom in while the operator reflects the action of a velocity component (entry in the second column) on a position component (first row).

For a matrix we now use for abbreviation , , and . Then, the redistributed version of (8) reads

 U=U0+QxvU0+QQxF(U)+QvF(U) (12)

with , and

 Qv :=Qv⊗Id=Q⊗Iv⊗Id (13a) QQx :=QQ⊗Ix⊗Id (13b)

Equation (12) can be compactly written as a possibly non-linear system of equations

 Mcoll(U)=CcollU0. (14)

with

 (15)

Here, and are matrices, while can in general be a non-linear operator. Setting

 Pcoll :=M−1collCcoll, (16)

the formal update for simply reads . We note that this formalism can be easily extended for higher-order ODEs: For an th-order ODE formulated as first-order system (as done here for ), the permutation operators (11) are simply the unit vectors and is replaced by at set of matrices that couple the components accordingly.

Evaluating requires the inversion of . Only for the sake of notational simplicity, we now focus on linear right-hand side functions , so that is a matrix with

 Mcoll(⋅)=Mcoll=I(M+1)2d−QcollF (17)

and inverse . However, we emphasize that the very same ideas and formulas described in the following apply for the case of non-linear functions, too, as e. g. shown in Section 4. Then, operators like have to be interpreted accordingly, c. f. the discussion at the end of A. For clarity, arguments of the (now linear) mapping are still shown with brackets.

In order to obtain a closed update formula which directly maps the initial data to the final value , we define the linear transfer operators and via

 TPu0=U0 and%  TRU=uM. (18)

Since is not necessarily the final step (if as e. g. for Gauss-Legendre nodes), we make use of the collocation formulation again to obtain approximations to the final values from the full vector or , respectively. We define by the extended vector of quadrature weights over the full interval , where

 ¯qm=∫tn+1tnℓm(τ) dτ. (19)

Then we obtain

 vn+1 =v0+qF(U), (20) xn+1 =x0+qV=x0+qV0+qQF(U), (21)

which can be combined into a single equation again using

 un+1=~CcollU0+~QcollF(U) (22)

with

 ~Ccoll:=TR+qxv and ~Qcoll:=qQx+qv. (23)

Note that if , the vector is equal to the last row of the matrix and computing is equivalent to computing . Now, the complete update formula for reads

 un+1=~Pcoll(u0):=~CcollTPu0+~QcollF(PcollTPu0). (24)

The subsequent parts of this section deal with the numerical properties of this formulation and point towards strategies for efficiently inverting , i. e. solving (14) by an iterative method.

### 2.2 Properties of collocation methods

A collocation method with nodes is equivalent to an -stage implicit Runge Kutta method (IRKM) with a Butcher tableau

 (25)

with being the vector of nodes scaled to the unit interval, see e. g. (hairer_nonstiff, , Theorem 7.7). In this interpretation, equation (14) is a system of equations to be solved for the stages of an IRKM while (22) is the actual update step to be performed once the stages are known.

Collocation methods have a number of attractive numerical properties: They are of optimal order, for Legendre and for Lobatto nodes. For both Gauss-Legendre and Gauss-Lobatto nodes, the resulting method is symmetric because the corresponding nodes are symmetric (hairer_geometric, , Theorem 8.9). Also, for Gauss-Legendre nodes, the resulting method is always symplectic (hairer_nonstiff, , Theorem 16.5) as well as B- and A-stable (hairer_stiff, , Theorem 12.9). In Section 4, we show that collocation methods with Lobatto nodes also have excellent stability properties for the Penning trap example considered there. For Lobatto nodes, however, the method is not necessarily symplectic (hairer_nonstiff, , Table 16.2). Nevertheless, for the cases studied here, the Hamiltonian is given as a quadratic form with a symmetric real matrix for which symmetric methods are also symplectic and vice versa (hairer_geometric, , Theorem 4.9). Hence, for the problems in Section 4, Gauss-Lobatto nodes also yield a symplectic collocation method and because they do not require an additional step to compute the final value at the end of the interval, we focus on Lobatto nodes here. For other cases, Legendre nodes might have to be used to obtain a symplectic method. Note that despite the collocation method being symplectic, energy drift can still emerge due to accumulation of round-off errors, see Hairer2008 ().

### 2.3 Picard iteration and preconditioning

Even in the linear case, the dense structure of calls for an an iterative approach to solve (14). The simplest iteration procedure is a Richardson iteration, see e. g. Kelley1995 (), reading

 U0 =U0 (26a) Uk+1 =(I(M+1)2d−Mcoll)Uk+CcollU0=CcollU0+QcollF(Uk). (26b)

for . Here, superscript denotes the iteration steps. As a measure for convergence, the norm of the residual

 Rk:=CcollU0−McollUk (27)

can be monitored, where . We note that (26) is equivalent to a discretized Picard iteration. Convergence depends on the eigenvalues of the iteration matrix . In the non-linear case, i. e. where is an operator, convergence properties are given by some adequate norm of this operator. Picard iteration typically converges only for very small time steps and is thus usually not an efficient approach, so that more advanced methods based on preconditioners are necessary. For a concise overview of iterative methods including the concept of preconditioning we refer to Kelley1995 (). Introduction of a preconditioner leads to the preconditioned Richardson iteration

 MpcUk+1=(Mpc−Mcoll)Uk+CcollU0. (28)

Here, each iteration requires to solve a linear or nonlinear system of equations determined by the preconditioner . The key is to find a good preconditioner: It has to be easy to invert so that computing (28) is cheap but still provides a sufficiently good approximation of to lead to robust and rapid convergence. In the next section, we will construct such a preconditioner out of the well-known Boris integrator for problems of the form (1).

## 3 Spectral deferred corrections based on the Boris-integrator (Boris-SDC)

The idea to interpret spectral deferred corrections as a preconditioned iterative scheme has been used for different purposes e. g. in HuangEtAl2006 (); Weiser2014 (). Here, we employ it to derive a problem-specific formulation of SDC based on the Boris integration method. The formulation follows the derivation of SDC for second-order problems in Minion2014_2ndOrderSDC (), using the standard velocity-Verlet integrator as preconditioner. There, problems are considered with a force field that depends only on the position . This corresponds to a separable Hamiltonian of a specific form. In general, an that also depends on the velocity leads to an implicit update for in the velocity-Verlet integrator, cf. (2). For the specific form of in (1), however, the Boris integrator provides a trick that essentially allows for a very efficient solution of the implicit system.

In order to use the Boris integrator as a preconditioner, we need a formulation of the velocity-Verlet scheme in matrix form, similar to (14). A concise summary of this rather tedious derivation is given in A.

### 3.1 Velocity-Verlet-based spectral deferred corrections

The standard velocity-Verlet integrator (2) for time steps , , can be written as system of equations

 MvvU=CvvU0, (29)

with system matrix

 Mvv:=I(M+1)2d−QvvF, (30)

see A, in particular (78) and (79), for details (we assume linear right-hand side here as well for notational simplicity). Since both and are lower block-diagonal matrices, (29) can easily be solved by forward substitution. Note that for non-linear functions , can no longer be compactly written in matrix form, but the system itself is still straightforward to solve. Hence, satisfies the conditions necessary for a suitable preconditioner: It can be easily inverted and approximates the original action of , as both systems correspond to equations for approximations of the values of and at the quadrature nodes, a high-order approximation in case of (14) and a composite low-order approximation in case of (29).

In order to precondition the iteration (26), we therefore apply the splitting

 Mcoll=Mcoll−Mvv+Mvv. (31)

Using as a preconditioner in (28) results in the iteration

 U0 =U0 (32a) (I(M+1)2d−QvvF)(Uk+1) =(Qcoll−Qvv)F(Uk)+CcollU0 (32b)

for with iteration matrix

 (33)

As the name suggests, this particular choice of preconditioner leads to the method of spectral deferred corrections (SDC) with velocity-Verlet as base integrator. The “direct” update matrix of (16) is then replaced by the sum of the updates given by the preconditioned iteration, i. e. we have

 Pksdc :=Kksdc+k−1∑l=0KlsdcM−1vvCcoll (34)

so that an approximation to with SDC iterations can be computed through

 (35)

Convergence of the SDC iteration can be monitored via the residual (27).

While the formulation with iteration and update matrices is convenient to analyze, a different approach is needed for an actual implementation to avoid the explicit use and storage of the full right-hand side vector . To this end, is again split according to (31), but the components and are then treated separately, so that (32) becomes

 Xk+1−QxF(Uk+1) =(QQ−Qx)F(Uk)+QV0+X0, (36a) Vk+1−QTF(Uk+1) =(Q−QT)F(Uk)+V0, (36b)

using the definitions of and of (72) and (75) in A. For a component-wise notation (in terms of nodes ), we define

 Q =:(qm,l)m,l=0,…,M, QQ =:(qqm,l)m,l=0,…,M, Qx =:(qxm,l)m,l=0,…,M, QT =:(qTm,l)m,l=0,…,M (37)

and obtain

 vk+1m+1 =v0 +M∑l=0qTm+1,l(f(xk+1l,vk+1l)−f(xkl,vkl)) +M∑l=0qm+1,lf(xkl,vkl), (38a) xk+1m+1 =x0+M∑l=0qm+1,lv0 +M∑l=0qxm+1,l(f(xk+1l,vk+1l)−f(xkl,vkl)) +M∑l=0qqm+1,lf(xkl,vkl) (38b)

for and . These variables are still vectors, i. e. we have . Note that the storage required for the -matrices (37) is usually negligible with respect to the size of . The following observations are useful to simplify these formulas:

• The matrix is a lower diagonal matrix, while is even strictly lower diagonal. Thus, the sum over the difference of the function values can be terminated at in (38a) and in (38b). The formula for is therefore fully explicit, the formula is semi-implicit. In the following section, this semi-implicit update will be reformulated in an explicit way for problems of the form (1) using the Boris ”trick”.

• Since , the first row as well as the first column of is zero. The last sum in both formulas can therefore start at , which is also true for the sum over the initial values in (38b).

• Independently of and , the first term in the summation over the difference of the function values is always zeros, because and . Thus, these sums can start at , too.

• The formulations using the initial conditions for each node are called “-to-node” formulations. More conveniently, this can be reformulated in “node-to-node” form using the matrix instead of , where the th row of is defined as the difference between the th and the th row of (starting with a row of zeros). The matrices and are defined analogously to and . Note that the sum over the th row of is equal to , which defines the factor in front of .

Based on the comments above, taking the difference between (38) for and gives

 vk+1m+1 =vk+1m +Δτm+12(f(xk+1m+1,vk+1m+1)−f(xkm+1,vkm+1)) +Δτm+12(f(xk+1m,vk+1m)−f(xkm,vkm)) +M∑l=1sm+1,lf(xkl,vkl), (39a) xk+1m+1 =xk+1m+Δτm+1v0 +m∑l=1sxm+1,l(f(xk+1l,vk+1l)−f(xkl,vkl)) +M∑l=1sqm+1,lf(xkl,vkl), (39b)

for and . This is the “node-to-node” formulation of SDC with velocity-Verlet integrator as base method. Equation (39) provides the formulation of Boris-SDC that would actually be implemented: Once values from iteration are known, the sums involving the quadrature weights can be computed and the step from to is then essentially a velocity-Verlet step with additional known terms on the right-hand side.

For (39), values for are provided by a simple copy of the initial value to all nodes, see (32). A -times application of these formulas provides approximations and to and . Both are then used to form , which in turn serves as input for the collocation formulation to approximate via

 ukn+1=~CcollTPu0+~QcollF(Uk), (40)

see Eqs. (22) and (35). Using Gauss-Lobatto nodes, i. e.  and , and only a single iteration, the formulas (39) yield the standard velocity-Verlet scheme (2): A brief calculation shows that for this case

 S=⎛⎜ ⎜⎝0000000Δτ22Δτ22⎞⎟ ⎟⎠,Sx=⎛⎜ ⎜⎝0000000(Δτ2)220⎞⎟ ⎟⎠,SQ=⎛⎜ ⎜⎝0000000(Δτ2)24(Δτ2)24⎞⎟ ⎟⎠ (41)

and . The collocation formula (40) is obsolete here (since ) but nevertheless valid with

 q=(0,Δτ22,Δτ22),andQ=S, (42)

see the definitions in (23), so that update and iteration formula result in the same expression. Hence, the first iteration of SDC with velocity-Verlet as base integrator on Gauss-Lobatto nodes is equivalent to applying velocity-Verlet on each node, as long as the initial value at each node is a copy of the initial value.

### 3.2 Boris-SDC

For the specific equations of motion under investigation in this work, the right-hand side as stated in (1a) is given by

 f(x,v)=α[E(x,t)+v×B] (43)

with a constant magnetic field for simplicity (see note on non-constant -fields at the end of this section). Thus, the update formulas for the velocity component in the standard velocity-Verlet integrator (2b) as well as in the SDC iteration (39a) seem to require the solution of an implicit system. For velocity-Verlet integration, the state-of-the-art approach is the Boris integration method Boris1970 (); Birdsall1985 (). Here, (2b) is rewritten as

 vn+1−vnΔt=α[En+12+vn+1+vn2×B]. (44)

The half-step subscript corresponds to the average electric field at times and , i. e. . The idea of the Boris integrator is to separate the electric and magnetic forces. To this end, we define

 v−:=vn+αΔt2En+12andv+:=vn+1−αΔt2En+12, (45)

so that

 v+−v−Δt =α2(v++v−)×B, (46)

which can be shown to correspond to a simple rotation, i. e.  and can be solved for explicitly (cf. Birdsall1985 ()) using with and . Using (45), the new velocity can therefore be computed explicitly. A relativistic generalization of this method is straightforward (Birdsall1985, , Chapter 15–4).

We can use and extend this idea to resolve the implicit dependency in the SDC iteration (39a). More precisely, we rewrite the (seemingly implicit) update for the th component of as

 (47)

Note that the second and third summand of the right-hand side of this equation only depend on values at iteration , i. e. these summands have been computed in the previous iteration. We define

 Δτm+1ck:=−Δτm+12(f(xkm+1,vkm+1)+f(xkm,vkm))+M∑l=1sm+1,lf(xkl,vkl), (48)

so that

 vk+1m+1−vk+1mΔτm+1=f(xk+1m+1,vk+1m+1)+f(xk+1m,vk+1m)2+ck. (49)

With the particular right hand side (43), this yields

 vk+1m+1−vk+1mΔτm+1=α[Ek+1m+12+vk+1m+1+vk+1m2×B]+ck, (50)

Except for the -term (which is known from the previous iteration) this has the very same structure as Equation (44). Extending the idea of (45), we define

 v−:=vk+1m+Δτm+12(αEk+1m+12+ck)andv+:=vk+1m+1−Δτm+12(αEk+1m+12+ck) (51)

to obtain

 v+−v−Δτm+1=α2(v++v−)×B, (52)

which is precisely of the type of (46). As noted before, this can be solved explicitly for , so that (51) can be used to determine . This gives us an explicit solver for the seemingly implicit SDC update (39a) and can be implemented directly into an existing SDC algorithm without further modifications.

We note that this approach can be easily extended to non-constant magnetic fields as follows. Instead of Eq. (44) we now have

 vn+1−vnΔt =α[En+12+12vn+1×Bn+1+12vn×Bn] (53) =α[En+12+vn+1+vn2×Bn+1]+α2vn×(Bn−Bn+1). (54)

The first part has again the same structure as Eq. (44), while the last part does not depend on and can thus be treated separately, i. e. as part of the -term in (51).

## 4 Numerical results

To evaluate the numerical properties of the Boris-SDC integrator, we study particles in a standard Penning trap Penning1936 (). Being confined to a limited volume due to an external magnetic and electric field, the particles’ characteristic properties, such as trajectories in real and phase space, energy conservation, and stability of the integration scheme, can conveniently be analyzed. We consider both the case of a single particle, where an analytic reference solution for the particle’s trajectory is available, as well as the case of many particles. All Boris-SDC runs use Gauss-Lobatto quadrature nodes, see the discussion in 2.2.

We follow the analysis in Patacchini2009 () and choose a constant magnetic field along the -axis with the particle’s charge-to-mass ratio so that

 (55)

The electric field experienced by a particle at position is composed of an ideal quadrupole potential distribution leading to

 Eext(xi)=−ϵω2Eα⎛⎜⎝10001000−2⎞⎟⎠xi (56)

and the inter-particle Coulomb interaction (cgs units)

 Eint(xi)=Nparticles∑k=1k≠iQkxi−xk(|xi−xk|2+λ2)\nicefrac32. (57)

To avoid numerical heating due to spurious close encounters, we follow the standard approach and regularize the Coulomb pole here by means of the smoothing parameter .

In the case , particles are confined in -direction due to the attractive nature of the potential. This setup corresponds to a Penning trap configuration. For , the particle will escape along the -axis but – due to the magnetic field – will still be confined to a limited orbit in - and -direction.

### 4.1 Single particle in a Penning trap

Solving the equations of motion (1) with the magnetic field (55) and the electric field (56) for a single particle inside the Penning trap is a standard textbook task, see e. g. Brown1986 (); Werth2009 (). Movement in direction is a harmonic oscillation, decoupled from the other coordinates,

 z(t) =z(0)cos(~ωt)+vz(0)~ωsin(~ωt), ~ω :=√−2ϵ⋅ωE. (58)

Here, the role of stated before becomes obvious: For , the frequency is purely imaginary and the -trajectory diverges, for the frequency and the trajectory corresponds to a harmonic oscillation. With the definition , the particle movement in the - plane is given by

 w(t) Ω± :=12(ωB±√ω2B+4ϵω2E) (59)
 R− :=Ω+x(0)+vy(0)Ω+−Ω−, R+ :=x(0)−R−, I− :=Ω+y(0)−vx(0)Ω+−Ω−, I+ :=y(0)−I−. (60)

Note that for (which can only happen for ), we have for the revolution frequency. In this case, the physical setup is unstable and the particle escapes from the trap due to a too weak magnetic or too strong electric field. Table 1 lists the physical parameters used in this section; Figure 1 shows a visualization of the particle’s analytical trajectory.

#### 4.1.1 Stability

As a first study, we analyze the stability of the classical Boris integrator, the collocation method, and Boris-SDC for the Penning trap and discuss the relation of the stability of Boris-SDC to the convergence of the corresponding SDC iteration. Note that the numerical stability of SDC has been extensively studied for the generic test equation with a single complex eigenvalue in DuttEtAl2000 (). For a modified version of the original Boris scheme, a numerical stability analysis for an ideal Penning trap has been performed in Patacchini2009 ().

We assess stability by the largest absolute value of the eigenvalues, that is the spectral radius, of the method’s update matrix. For Boris-SDC this is the matrix defined in (35) and for the collocation method defined in (24). A method is stable for a specific configuration, if and only if the spectral radius is smaller than or equal to unity.

Figure 2 shows the resulting stability region for the classical Boris integrator and a collocation method with . The hatched, dark gray area on the left indicates the physical instability of the system, where so that the magnetic field is too weak to confine the particle. For the classical Boris integrator, there is in addition a zone of numerical instability (light gray), where the method is not stable although the physical problem already is. Furthermore, a significant area of numerical instability is present for values of on the right side. In contrast, the stability domain of the collocation method is identical to the domain where the physical setup is stable: There is no region of additional numerical instability. For , the stability domain of the collocation method is identical (not shown). The collocation method is hence stable for every physically stable configuration of the single-particle Penning trap.

##### Convergence of Boris-SDC iteration

For Boris-SDC, the convergence properties of the iteration computing the collocation solution have significant impact on the stability: In cases where the iteration converges poorly or not all, we cannot expect to recover the stability properties of the collocation solution. Convergence of the SDC iteration is governed by the spectral radius of the SDC iteration matrix (not to be confused with the Boris-SDC update matrix). If and only if the spectral radius , the iteration will ultimately converge and the norm of the residual (27) will go to zero. For values close to unity, however, convergence can be unfeasibly slow, resulting in a very large number of required iterations. Figure 3 shows the spectral radius of for and collocation nodes. Small values (blue) indicate fast convergence, values close to unity (yellow and light red) slow convergence and values larger than unity (dark red) indicate divergence of the SDC iteration. For , the area where Boris-SDC shows good convergence does roughly coincide with the stability domain of the classical Boris integrator. Interestingly, for small values of , Boris-SDC also converges well in the region of physical instability. For