A highorder Boris integrator
Abstract
This work introduces the highorder BorisSDC method for integrating the equations of motion for electrically charged particles in an electric and magnetic field. BorisSDC relies on a combination of the Borisintegrator 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 loworder method. In BorisSDC, the Boris method, a secondorder Lorentz force integrator based on velocityVerlet, 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 particlebased plasma physics simulations involving magnetic fields, to a highorder method. Stability, convergence order and conservation properties of the method are demonstrated for different simulation setups. BorisSDC reproduces the expected high order of convergence for a single particle and for the centerofmass of a particle cloud in a Penning trap and shows good longterm energy stability.
keywords:
Boris integrator, time integration, magnetic field, highorder, spectral deferred corrections (SDC), collocation methodsort&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 highintensity 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
(1a)  
(1b) 
with the particle position , its velocity , the magnetic field , electric field and the charge to mass ratio . In (1a), the wellknown Lorentz force depends on and so that a discretizing (1) with a standard velocityVerlet scheme Verlet1967 (); Birdsall1985 () reads
(2a)  
(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 defacto 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 velocityVerlet scheme, the Boris approach is a secondorder 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 phasespace volume preserving. Furthermore, it only requires a single evaluation of the righthand 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 ParticleInCellcodes (see e. g. Verboncoeur2005 ()), gridfree methods (e. g. Gibbon2010a ()) and MonteCarlo simulations (e. g. Kirschner2000 ()). Several explicit alternatives to the Boris method have been proposed, compare Patacchini2009 () and references therein. All of them are secondorder accurate and apparently no higherorder methods based on the original Boris approach exist. Especially for applications such as trajectory integration in particle accelerators Toggweiler2014 (), spaceweather studies Lapenta2012 (), highintensity laserplasma 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 highperformance computing systems towards high floating point operation rates at stagnating memory data transfer speed favors the use of higherorder 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 highprecision or higherorder 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 fourthorder RungeKutta method in beam propagation simulations Stoltz2002 (), a quasisymplectic Trotterfactorization based scheme that builds upon an explicitimplicit mixture of leapfrog, Verlet, exponential differencing and Boris rotations with a sixthorder rotation angle approximation Bowers2009 () or a Taylor seriesbased explicit approach with an up to sixthorder 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 seriesbased 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 highorder BorisSDC integration method for (1), which is a combination of the classical Boris integrator with the spectral deferred corrections method (SDC). The resulting BorisSDC 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 firstorder ODEs as an iterative approach for the generic construction of highorder integration schemes using a loworder base propagator (the “sweeper”) such as implicit or explicit Euler for the correction “sweeps” DuttEtAl2000 (). Several modifications and extensions exist e. g. semiimplicit SDC Minion2003 (), GMRESaccelerated SDC HuangEtAl2006 (), inexact SDC SpeckEtAl2014_DDM2013 (), multilevel SDC SpeckEtAl2014_BIT (), SDC based on DIRK methods Weiser2014 () or the parallel full approximation scheme in space and time (PFASST), a parallelintime integrator which exploits the iterative structure of SDC EmmettMinion2012 (). Recently, SDC has been formulated for secondorder 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 velocitydependent force of the form (1a).
The derivation of BorisSDC 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 RungeKutta 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 loworder method. If sufficiently many sweeps are performed, the advantageous properties of the underlying collocation method are recovered. For e. g. a firstorder 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 velocityVerlet 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 righthand side and derive readytoimplement expressions for the BorisSDC method, specifically tailored for problems of the form (1). Section 4 illustrates the properties of BorisSDC by numerical examples and compares BorisSDC 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 velocityVerlet integrator as base method for secondorder problems in Minion2014_2ndOrderSDC ().
Rewriting equations (1) in Picard formulation for an arbitrary interval with starting value , , , we obtain^{1}^{1}1We use the following notation here: Vectors and matrices in normal font refer to scalar values at a single node in time. Boldfaced 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.
(3a)  
(3b) 
Collocation methods are based on the introduction of intermediate nodes
(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
(5) 
where , are Lagrange polynomials. Then, in order to account for the initial values , , we further define
(6) 
For quadrature nodes with , e. g. GaussLobatto 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 righthand side values. Then, the matrix provides approximations of the integral in (3a) with , that is
(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
(8a)  
(8b) 
for
(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 degreesoffreedom are ordered as
(10) 
To this end, we introduce permutation operators , and with
(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
(12) 
with , and
(13a)  
(13b) 
Equation (12) can be compactly written as a possibly nonlinear system of equations
(14) 
with
(15) 
Here, and are matrices, while can in general be a nonlinear operator. Setting
(16) 
the formal update for simply reads . We note that this formalism can be easily extended for higherorder ODEs: For an thorder ODE formulated as firstorder 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 righthand side functions , so that is a matrix with
(17) 
and inverse . However, we emphasize that the very same ideas and formulas described in the following apply for the case of nonlinear 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
(18) 
Since is not necessarily the final step (if as e. g. for GaussLegendre 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
(19) 
Then we obtain
(20)  
(21) 
which can be combined into a single equation again using
(22) 
with
(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
(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 GaussLegendre and GaussLobatto nodes, the resulting method is symmetric because the corresponding nodes are symmetric (hairer_geometric, , Theorem 8.9). Also, for GaussLegendre nodes, the resulting method is always symplectic (hairer_nonstiff, , Theorem 16.5) as well as B and Astable (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, GaussLobatto 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 roundoff 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
(26a)  
(26b) 
for . Here, superscript denotes the iteration steps. As a measure for convergence, the norm of the residual
(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 nonlinear 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
(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 wellknown Boris integrator for problems of the form (1).
3 Spectral deferred corrections based on the Borisintegrator (BorisSDC)
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 problemspecific formulation of SDC based on the Boris integration method. The formulation follows the derivation of SDC for secondorder problems in Minion2014_2ndOrderSDC (), using the standard velocityVerlet 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 velocityVerlet 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 velocityVerlet scheme in matrix form, similar to (14). A concise summary of this rather tedious derivation is given in A.
3.1 VelocityVerletbased spectral deferred corrections
The standard velocityVerlet integrator (2) for time steps , , can be written as system of equations
(29) 
with system matrix
(30) 
see A, in particular (78) and (79), for details (we assume linear righthand side here as well for notational simplicity). Since both and are lower blockdiagonal matrices, (29) can easily be solved by forward substitution. Note that for nonlinear 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 highorder approximation in case of (14) and a composite loworder approximation in case of (29).
In order to precondition the iteration (26), we therefore apply the splitting
(31) 
Using as a preconditioner in (28) results in the iteration
(32a)  
(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 velocityVerlet 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
(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 righthand side vector . To this end, is again split according to (31), but the components and are then treated separately, so that (32) becomes
(36a)  
(36b) 
using the definitions of and of (72) and (75) in A. For a componentwise notation (in terms of nodes ), we define
(37) 
and obtain
(38a)  
(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 semiimplicit. In the following section, this semiimplicit 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 “tonode” formulations. More conveniently, this can be reformulated in “nodetonode” 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
(39a)  
(39b) 
for and . This is the “nodetonode” formulation of SDC with velocityVerlet integrator as base method. Equation (39) provides the formulation of BorisSDC 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 velocityVerlet step with additional known terms on the righthand 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
(40) 
see Eqs. (22) and (35). Using GaussLobatto nodes, i. e. and , and only a single iteration, the formulas (39) yield the standard velocityVerlet scheme (2): A brief calculation shows that for this case
(41) 
and . The collocation formula (40) is obsolete here (since ) but nevertheless valid with
(42) 
see the definitions in (23), so that update and iteration formula result in the same expression. Hence, the first iteration of SDC with velocityVerlet as base integrator on GaussLobatto nodes is equivalent to applying velocityVerlet on each node, as long as the initial value at each node is a copy of the initial value.
3.2 BorisSDC
For the specific equations of motion under investigation in this work, the righthand side as stated in (1a) is given by
(43) 
with a constant magnetic field for simplicity (see note on nonconstant fields at the end of this section). Thus, the update formulas for the velocity component in the standard velocityVerlet integrator (2b) as well as in the SDC iteration (39a) seem to require the solution of an implicit system. For velocityVerlet integration, the stateoftheart approach is the Boris integration method Boris1970 (); Birdsall1985 (). Here, (2b) is rewritten as
(44) 
The halfstep 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
(45) 
so that
(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 righthand side of this equation only depend on values at iteration , i. e. these summands have been computed in the previous iteration. We define
(48) 
so that
(49) 
With the particular right hand side (43), this yields
(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
(51) 
to obtain
(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 nonconstant magnetic fields as follows. Instead of Eq. (44) we now have
(53)  
(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 BorisSDC 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 BorisSDC runs use GaussLobatto 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 chargetomass ratio so that
(55) 
The electric field experienced by a particle at position is composed of an ideal quadrupole potential distribution leading to
(56) 
and the interparticle Coulomb interaction (cgs units)
(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,
(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
(59) 
(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 BorisSDC for the Penning trap and discuss the relation of the stability of BorisSDC 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 BorisSDC 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 singleparticle Penning trap.
Convergence of BorisSDC iteration
For BorisSDC, 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 BorisSDC 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 BorisSDC shows good convergence does roughly coincide with the stability domain of the classical Boris integrator. Interestingly, for small values of , BorisSDC also converges well in the region of physical instability. For