Splitting methods for time integration of trajectories in combined electric and magnetic fields
The equations of motion of a single particle subject to an arbitrary electric and a static magnetic field form a Poisson system. We present a second-order time integration method which preserves well the Poisson structure and compare it to commonly used algorithms, such as the Boris scheme. All the methods are represented in a general framework of splitting methods. We use the so-called functions, which give efficient ways for both analyzing and implementing the algorithms. Numerical experiments show an excellent long term stability for the new method considered.
pacs:02.70.-c, 45.10.-b, 52.65.Cc, 52.65.Yy
The time integration of a particle motion due to electric and magnetic fields is a common task in several fields of science, from molecular dynamics to plasma physics and astrophysics. Various integration methods have been established within different scientific areas. A prominent example is the leap frog method, which is also called Boris pusher in the context of particle-in-cell (PIC) codes.
We consider here time-independent magnetic fields which, however, may be nonuniform in space. The electric field on the other hand is not subject to any constraints. A classical Coulomb system subject to a magnetic field serves as an example. If the particles within this system are moving sufficiently slow, the pairwise interaction can be described purely electrostatically and only the time invariant external magnetic field has to be taken into account.
In this article we review some second-order methods from the literature under a unifying framework: the methods are interpreted as splitting methods and the solutions of the subflows are given in terms of the so-called functions. This offers an economic way of stating, analyzing and implementing the methods. We also present a new second-order method based on the symmetric splitting of the Hamiltonian of the system. Further we discuss the structure preserving character of the methods regarding time symmetry, volume preservation and preservation of the system’s underlying Poisson structure.
The article is organized as follows. In Section II we formulate the equations of motion of a particle subject to an electric and a static magnetic field as a Poisson system, and we state some important properties of the exact flow of the system. In Section III we shortly discuss the corresponding properties we desire from the discrete flow given by the numerical integrator. In Section IV we give formulas for the maps corresponding to the subflows of different terms of the Poisson system and shortly discuss their structural properties. In Section V we present four time integrators, which are constructed using the formulas of Section IV and certain functions. In Section VI we give numerical results in which the integrators are compared on four different cases. The needed formulas for the matrix functions are given in the appendix (Section VIII).
Ii Equations of motion and their description as a Poisson system
The equations of motion for a particle with mass and charge subject to an electromagnetic field are given by the Lorentz force
where the electric field and the magnetic field are determined by a scalar potential and a vector potential , respectively, as
In the Lagrangian formalism the Lorentz force (1) can be derived from the Lagrangian
which leads to the usual Hamiltonian formulation
However, we use a different formulation of the equations of motion, which exhibits both the Hamiltonian structure of the system and the splitting of the right-hand side into efficiently computable subflows. We use for the momentum and state the equations of motion as
where is the electric force and
corresponds to the cross product of the magnetic force as
We use the matrix formulation of the cross product as it will simplify the analysis and derivation of the time integration methods.
is called a Poisson system, if is a skew-symmetric matrix which satisfies the identity
In this case, is also called a Poisson tensor. The exact flow of any Poisson system is a Poisson map (also called noncanonically symplectic map), i.e., the Jacobian matrix of the flow satisfies
Moreover, it is time symmetric, i.e. .
and the skew-symmetric structure matrix
by setting . We will denote from now on the conjugate variables as or , whichever is more convenient.
For the system matrix (9) we have the following. Let be a differentiable function, and let the matrix (4) be defined by . Then a direct calculation using the relation (6) shows that as in (9) defines a Poisson tensor if and only if . This is satisfied as Maxwell’s laws are assumed to hold for the magnetic field (e.g., we do not consider magnetic monopoles). Thus the ODE (3) can be written as a Poisson system (5).
The flow of a Hamiltonian system is volume preserving, i.e., it holds that
as . Since is a continuous function of , and since , (10) holds also for the exact flow of the Poisson system.
As we consider only the motion of a single particle, the electrostatic force will be time-independent. In case of many particles and pairwise forces this is not the case. However, the Poisson system (3) generalizes from dimension to in a trivial way: is replaced by a block diagonal matrix consisting of blocks.
Iii Structure preserving integrators
It is desirable that the one-step map given by the numerical flow satisfies the same structural properties as the exact flow of the ODE. The preservation of time symmetry or the Poisson structure implies good long time behavior for the numerical solution (see the backward error analysis of HairerLubichWanner ()). The properties we want from the integrators are:
Time symmetry: .
The integrator gives a Poisson map, i.e., it holds that
The preservation of volume:
Whether (11) is satisfied depends both on the integrator and the Poisson system considered. Moreover, the integrator has to respect the Casimirs (the first integrals of the Poisson system), which in general is difficult to verify for a nonlinear system matrix .
However, in some special cases one can show that these properties hold. That is the case for example if the integrator is a composition method resulting from a Hamiltonian splitting. This will be discussed in Section IV.
We will study analytically in Section V whether the considered integrators satisfy these properties.
Iv Splitting methods
To have a common description for the various time integrators we state them in the framework of splitting methods and use the functions for representing the subflows containing the term . To this end, we write the Poisson system (3) as
We denote the exact flows of the subsystems corresponding to , , , etc. by , , , etc. This means that provides the exact solution of the system at time , and so on. We further denote by the exact flow of the ODE with right-hand side , where . This means that the magnetic field is fixed at , i.e. . One easily verifies the following solutions for the exact flows:
We will also need the solutions of the subsystems corresponding to the fields and . The flow corresponding to is given by the exact solution of the ODE
which can be expressed by the matrix exponential as
Using (42) we see that
On the other hand, the flow corresponding to the field is given by the exact solution of the ODE
This solution is given by the variation-of-constants formula (44) as
We also state the solution for the flow corresponding to in the particular case when both and are constant. This is presented also in Chin () using another formulation. The ODE corresponding to the flow is now
and its solution is given again by the variation-of-constants formula and by formula (42) as
We note that all of these substeps are time symmetric, i.e., , and so on.
where for and for . Its determinant is given by
As the Cayley transform of a skew-symmetric matrix is a unitary matrix, and as is a continuous function of and , we find for the maps and the relation
Therefore both are volume preserving. One easily verifies from (12) that and are also volume preserving.
Splitting the Hamiltonian (8) gives subsystems corresponding to the right-hand sides
Both of these constitute Poisson systems of the form (5) with the system matrix (9) and with Hamiltonians and , respectively. Therefore, the subflows and are Poisson maps, as well as their compositions , etc.
We will next state four time integration methods. They are formulated using the functions and the solutions of the subflows given in Section IV. Three of them can be found in the literature: the Boris–Buneman scheme BirdsallLangdon (), the method of Chin Chin () and the method of Spreiter and Walter SpreiterWalter (). We also state a new method which is nearly symmetric and nearly preserves the Poisson structure (method C below).
v.1 Boris–Buneman scheme
The so-called Boris pusher BirdsallLangdon () is undoubtedly the most common particle trajectory integrator used for PIC codes. It is a second-order accurate centered difference leap frog scheme, given by
PIC codes usually define the position with an offset of half a time step relative to the momentum. In (17) is updated alternatingly with .
To formulate this scheme as a one-step method, we denote , , and so on. Using this notation the linear drift operation corresponding to can be split as
The momentum change in between these drift operations is defined by the second equation of (17). Solving it for gives
where is the Cayley transform
From this one can easily read off that the mid-step can be split as
where , and is given by
Since the Cayley transform of a skew-symmetric matrix is a unitary matrix, the subflow is also a volume preserving map (see the analysis of Section IV). And as also the maps and are volume preserving, the method (20) is volume preserving as a composition of volume preserving maps. The volume preservation property of the Boris scheme was also shown in Qinetal ().
Replacing the Cayley transform by the exponential function in the mid-step, i.e., considering the splitting
gives also a symmetric and volume preserving method. Both the Cayley transform and the matrix exponential rotate the momentum around the vector in the mid-step. This can be seen as follows. Since , we may deduce from the power series representations that
On the other hand, since is skew-symmetric both and are unitary matrices, so that
v.2 Chin’s schemes
Since the substeps are symmetric maps, both (24) and (25) give symmetric methods. Moreover, as the methods (24) and (25) are compositions of volume preserving maps (see Section IV), both of them are volume preserving as well.
We note that the methods of Chin become exactly energy preserving when Chin ().
v.3 Symmetric splitting of the Hamiltonian
As a new method, we propose a symmetric splitting of the Hamiltonian for the Poisson system (5). This means that the right-hand side of the (5) is split into the flows (15) and (16) according to the Hamiltonian (8). Then, Strang splitting is applied to this decomposition giving the symmetric one-step method
The mid-step of (26), i.e., the solution
can be computed efficiently using the formula (13) if the magnetic field is uniform. If the magnetic field is position dependent, some approximations have to be made.
We infer from (26) that the last substep changes only . As the magnetic field depends only on , we see from (13) that by using a symmetric approximation for (27), we get a symmetric approximation for (26).
The approaches we consider for the mid-step are the following ones:
The symmetric method
combined with fixed-point iteration. This approach can be found in EinkemmerOstermann ().
The symmetric method
combined with fixed-point iteration.
In the methods above the numerical strategy is to perform a fixed-point iteration for the implicit mid-step. This approach is also taken in EinkemmerOstermann (). For the approach (28) this means that after the initial value
is computed, the iteration goes on as
Here can be computed explicitly using the formula (13).
The implementation for (29) goes analogously. After the initial value
is computed, the iteration goes on as
As shown in EinkemmerOstermann () the method (26) with the approach (28) will be symmetric up to an order that equals the number of fixed-point iterations. It is easy to verify that the same holds true for (29). By the results of EinkemmerOstermann () the same holds true for the third approach.
It is easily verified that the high-order composition scheme will also be a Poisson integrator up to an order of the scheme used for the mid-step (i.e., the relation (11) will be satisfied up to this order). Although the number of magnetic field evaluations increases considerably, the term is still evaluated only once per time step, which may be beneficial in situations where the calculation of the electric field is more expensive than that of the magnetic field.
Since the flows of the subsystems resulting from a Hamiltonian splitting are Poisson maps, so is the composition (26); see also (HairerLubichWanner, , Ch. 7). As the method (26) is symmetric and a Poisson integrator, it has a modified Hamiltonian (see the backward error analysis of HairerLubichWanner ()) and preserves the first integrals for exponentially long times.
The strategies 1-3 above give symmetric (or Poisson) maps up to a high order and we expect a good preservation of first integrals numerically.
v.4 The method of Spreiter and Walter
In SpreiterWalter () a second-order accurate scheme is proposed which incorporates a homogeneous magnetic field in -direction, and which overcomes the stability condition , where is the cyclotron frequency. We repeat here the derivation of this method using the functions which shortens the derivation considerably. As this derivation shows, the method can be seen as an exponential version of the velocity Verlet scheme. Integrating the Lorentz force law
with respect to gives an ODE for ,
Under the assumption that , we can approximate the ODE (30) up to second order as
By the variation-of-constants formula (44) the solution for this ODE at is given by
This gives the time stepping formula for . Integrating the Lorentz force law and using the variation-of-constants formula (44) once more yields
Using the fact that
we find that
Since , and as , we see that in the limiting case the velocity Verlet scheme is obtained. We note that the method of Spreiter and Walter can be interpreted also as an exponential Taylor method KoskelaOstermann ().
Despite of its similarity with the velocity Verlet scheme, we have not found any structure preserving properties for the method of Spreiter and Walter.
Vi Numerical comparison of the integrators
In the first experiment, we compare the integrators for a particle motion in an inhomogeneous magnetic field in the absence of an electric field. Then, results for experiments including combined electric and magnetic fields are provided. The motion of a charged particle in a Penning trap Kretzschmar () provides a suitable numerical test in the case of combined fields. In this setting the particle trajectory is a periodic orbit, formed by a superposition of three harmonic oscillators. In addition to the ideal Penning trap with a uniform magnetic field we consider two variations of this setup with static but nonuniform magnetic fields.
vi.1 2d particle motion in an inhomogeneous magnetic field
In the first experiment we consider an example provided in Chin (): a particle of charge and mass at initial position with initial velocity moves in a magnetic field , . The particle motion is composed of a periodic motion in the -plane and a drift in -direction with drift velocity . The drift reduced trajectory forms an ellipse centered at with . is given by
One full orbit has the period .
The scalar potential is set and the vector potential with . The Lagrangian takes the form . As it is independent of , the -component of the generalized momentum,
gives another integral of motion.
Since now , the flow vanishes and the integrators used in the previous example simplify as follows.
Boris–Buneman scheme with Cayley transformation:
Chin-b (equals Boris–Buneman scheme with matrix exponential):
The method of Scovel gave poor results as it is not symmetric when applied to a nonuniform field. Therefore it was discarded in this experiment. The symmetric methods (28) (denoted as Impl. Strang) and (29) (denoted as Impl. mp.) were used with five (for Impl. Strang) and six (for Impl. mp.) fixed-point iteration steps.
We compare these methods by performing a numerical integration over 20000 cycles of period with step sizes varying from to . The initial velocity is set to in all test cases. Figures (a)a and (b)b show the time step length plotted against the relative error of the invariant (as defined in (33)) and the absolute position error using a logarithmic scale. Both iterative methods provide slightly lower errors of the invariant, while only the implicit Strang splitting method has a lower position error than the method of Chin. The method of Boris preserves well first-order drifts. When the step size increases almost about one order of magnitude there is only a little impact on the position error for the method of Boris. The position error stays around , which is the diameter of the ellipse in –direction. For small time steps, however, the implicit Strang splitting method is the most accurate.
The error behavior becomes more visible when considering separately the error of the drift velocity and the error of the period length as shown in Figures (a)a and (b)b. The Boris pusher is the best method regarding the conservation of the drift velocity but the worst regarding the period length.
vi.2 Ideal Penning trap with uniform magnetic field
The ideal Penning trap consists of a static quadrupole electric field and a uniform magnetic field given by
where depends on the geometry and the voltage of the electrodes. The analytic solution for the motion of a single particle in an ideal Penning trap can be found in Kretzschmar (). Choosing the trap parameters suitably gives a stable periodic orbit consisting of a fast gyromotion perpendicular to with the frequency , a slower axial motion along with the frequency solely caused by the electrostatic fields and an even slower circular drift motion around the -axis with the magnetron frequency . We set , , , and the initial coordinates of the particle and . The particle trajectory of this example is shown in Figure (a)a.
As the magnetic field is uniform in space the central step of the symmetric Hamiltonian splitting (26) can be computed exactly as is done in Scovel’s method. We compare this method to the method of Chin (25), the method of Spreiter and Walter (denoted by SpW), the Boris pusher (20) and its matrix exponential variation (21) (denoted by Boris exp).
We compare the methods by their ability to deliver a stable orbit over a large number of magnetron cycles (full circles around the -axis), by the relative energy error (Figure (a)a) and by the absolute error of the position (Figure (b)b). As both types of errors are heavily oscillating, the maximum value of the error along a trajectory with a total simulation time of (a complete magnetron cycle) is plotted against the time step length, using logarithmic scales for both quantities. The time step size ranges from 0.002 to 2.4 cyclotron cycles , violating the stability criteria for some methods. A further reference line with a slope corresponding to a second-order method has been added to both plots.
With the exception of Spreiter and Walter all methods feature a bounded energy error and a position error increasing linearly in time on a closed orbit, as long as the time step size is considerably smaller than . The method of Spreiter and Walter shows a linearly increasing energy error until its orbit becomes unstable, so with longer simulation time its errors become worse relative to the other methods.
For larger time steps the method of Chin, the Hamiltonian splitting (Scovel) and the Boris pusher in its matrix exponential variation become unstable for , integer, while the Boris pusher using the Cayley transform preserves a closed orbit and a bounded energy error.
The Boris pusher provides in both variations a smaller energy error than the method of Chin and the Hamiltonian splitting, but shows a larger position error for small time steps. Figure (b)b demonstrates an interesting behavior of the Boris pusher: for a large range of time step sizes the maximum position error remains constant at , the diameter of the gyromotion. This can be explained by the fact that the Boris pusher can reproduce the correct first-order gyrocenter drift motion ParkerBirdsall ().
vi.3 Penning trap with magnetic bottle
The third experiment features a nonuniform magnetic field: a Penning trap with a magnetic bottle determined the magnetic field
Here is as given in (35). A physical implementation of a Penning trap featuring such a magnetic field can be found in Rodegheri (). We set and all the other parameters are as in the last example. The strength of the magnetic field and the cyclotron frequency vary only a little along the trajectory. The particle trajectory of this example is shown in Figure (b)b.
When using the symmetric splitting method (26), approximations have to be made for the mid-step since now the magnetic field is nonuniform in space. Figure (a)a shows the effect of the number of employed fixed-point iterations on the energy preservation of this method in a long time simulation of a trapped particle using a time step size of , where is the cyclotron frequency at the initial position. As the position error is bounded by the closed orbit we choose the relative error of (defined as the angle of motion around the -axis) as a measure for the trajectory’s error in Figure (b)b. Both errors are heavily oscillating, so the maximum value over each consecutive time steps is drawn.
The mid-step is either approximated by the implicit method (28) (denoted as Impl. Strang) or by (29) (denoted as Impl. mp.). For 5 fixed-point iterations (in the case of Impl. Strang) and 6 iterations (in the case of Impl. mp.) both alternatives seem to become approximately energy preserving and provide a stable orbit. This is expected from the theory, since the order of symmetry of these methods equals the number of fixed-point iterations employed EinkemmerOstermann () and since a good long time behavior can be expected for symmetric methods HairerLubichWanner ().
We compare the methods of Chin, the Boris pusher and the symmetric splitting method (26) using the implicit schemes as described above. Furthermore, we apply a high-order composition of (28) or (29) (denoted as Impl. mp. comp.). The composition is performed using an 8th-order symmetric scheme of Suzuki and Umeno SuzukiUmeno () consisting of 15 substeps. Each substep uses 16 fixed-point iterations. Consequently, such a composition scheme requires a large number of magnetic field evaluations per time step: 240 for the method (29), or 255 when using the method (28) due to the calculation of the initial value within each substep. Nevertheless, a single electric field evaluation suffices for each time step. As the results of the implicit schemes do not differ considerably when the composition is applied, only the ones of (29) are drawn in the plots.
As before we provide two log-log plots showing the dependence of the energy and the position errors on the time step size in the range of 0.002 to 0.4 cyclotron cycles . Lacking an analytic solution to determine the position errors we use the trajectory of a high-order method with a small time step (0.00025 cyclotron cycles) as a reference: the fourth-order method of McLachlan () applied on and (denoted as M in Chin ()) in combination with a composition scheme Suzuki_Fractals () increasing the order to 6.
All variations of the symmetric Hamiltonian splitting method give an energy error similar to Chin’s scheme (Figure (a)a). However, the smallest position error is obtained by using the composition scheme for the mid-step of the Hamiltonian splitting.
vi.4 Penning trap with asymmetric magnetic field
As the Penning trap with magnetic bottle is still rotationally invariant, we also consider an example with an asymmetric magnetic field given by
where . Notice that . is again as given in (35). The particle trajectory of this example is shown in Figure (c)c. Along the trajectory the magnetic field strength increases by a factor of 1.8. Therefore time step sizes between 0.001 and 0.1 are investigated in a simulation lasting for approximately 24 orbit cycles. Figures (a)a and (b)b show the energy and position errors for the methods giving similar results as in the previous example. The reference trajectory for the position errors has been computed with the method described in the previous example. Using the largest time step size shows unbounded energy errors for all the methods with the exception of the Boris pusher with the Cayley transform and the fixed point iteration schemes of the symmetric Hamiltonian splitting in combination with the high-order composition method of the mid-step. Therefore another long time integration over time steps of length 0.1 has been performed using both methods and the fourth-order method of McLachlan (); Chin () (denoted by Chin M). For this example we measure also the magnetic moment defined by
It is a so-called adiabatic invariant of the system HairerLubichWanner () and should be conserved approximately. The magnetic moment as well as the energy error are oscillating within some range, so the maximum value on consecutive intervals of time steps is drawn in Figures (a)a and (b)b.
While the energy error of the Boris pusher remains smaller than that of the symmetric Hamiltonian splitting for over time steps, the evolution of the magnetic moment indicates the instability long before.
An initial decrease of is followed by a rapid growth after time steps for the Boris scheme. The Boris scheme with 2 times smaller time steps (denoted Boris h/2) shows a similar behavior with a delay. The fourth-order method denoted by Chin M exhibits a quickly growing error in energy but a smoother change of . Only the symmetric Hamiltonian splitting in combination with a high-order approximation of the mid-step is able to provide a stable orbit.
We have considered here time integrators for particle trajectories in combined electric and magnetic fields. We have reformulated and analyzed existing methods found from literature under a common framework and additionally proposed a new method which deals the effect of the magnetic field implicitly. This is done by splitting the right-hand side of the equations of motion into three terms and by using the matrix function approach. The matrix functions are used for the magnetic field term, and the so-called Rodrigues formula simplifies their evaluation. We represent the equations of motion as a Poisson system and analyze the structure preserving properties of the integrators.
The existing methods we found in literature are: Boris–Buneman scheme, Chin’s schemes, the method of Spreiter and Walter, and the Scovel’s method. Additionally we have proposed a new implicit scheme. This scheme is the same as the Scovel’s method when the magnetic field is homogeneous.
The Boris–Buneman scheme and the Chin’s schemes are shown to be both symmetric and volume preserving. For the Spreiter and Walter’s scheme we did not find any structure preserving properties. The new method we propose is fully structure preserving up to an order equaling the number of fixed-point iterations used for the implicit mid-step. In numerical experiments we illustrate the effect of the number of fixed-point iterations on the preservation of the structure. In our experiments 6 iterations were enough to reduce the error negligible.
As a first numerical comparison, we tested all the schemes except Spreiter and Walter’s method in a case where there is no electric field. Due to the time symmetry, all the methods seemed to preserve invariants. For the position error the newly proposed implicit splitting method seemed to be the most accurate for small time steps, whereas for larger time steps the method of Boris was the most accurate. When the motion was separated into a periodic and a drift term, it was found that the method of Boris preserved the drift motion the best and the periodic motion the worst, whereas this was the opposite for the implicit splitting method.
To compare schemes in combined electric and magnetic fields, we integrated particle trajectories within a Penning trap. As the first example of combined fields we considered a constant homogeneous magnetic field, so Scovel’s method could be used producing energy and solution errors very similar to the method of Chin: a lower solution error than the Boris pusher for step-sizes considerably smaller than the inverse cyclotron frequency at the cost of a slightly larger but likewise bounded energy error. Spreiter and Walter’s method on the other hand did not feature a bounded energy error, so the orbit of the particle was unstable.
As a second Penning trap example we have an inhomogeneous but radially symmetric magnetic field. Replacing Scovel’s scheme by the new method still gave the same energy error like Chin’s scheme. The solution error differed depending on the method used for the implicit mid-step iteration. Applying a high order composition scheme on the mid-step provided the lowest solution error of all methods tested.
In the last numerical example of a Penning trap, we performed long time integration of a trajectory in an asymmetric magnetic field and compared three methods: one of the Chin’s schemes, the Boris–Buneman scheme and the newly proposed implicit method. We followed the preservation of the energy and of the magnetic moment. As expected from the analysis of Section V, the implicit method turned out to be the most stable. The Boris–Buneman scheme, being the second most stable, was more unstable than the implicit scheme even for 4 times smaller time steps.
Viii Appendix: Matrix functions
The following identities are easily verified for the skew-symmetric matrix (see (4)):
where . Using the second identity and the Taylor series of and , we get the Rodrigues formula for the exponential of :
We list here the necessary definitions and properties of the functions. These functions are defined by
They satisfy and the recurrence relation