Time-centered, hence second-order, methods for integrating the relativistic momentum of charged particles in an electromagnetic field are derived. A new method is found by averaging the momentum before use in the magnetic rotation term, and an implementation is presented that differs from the relativistic Boris Push only in the method for calculating the Lorentz factor. This is shown to have the same second-order accuracy in time as that found by splitting the electric acceleration and magnetic rotation (Boris Push) and that found by averaging the velocity in the magnetic rotation term (Vay's method) [J.-L. Vay, Phys. Plasmas 15, 056701 (2008)]. All three methods are shown to conserve energy when there is no electric field. The Boris Push and the current method are shown to be volume-preserving, while Vay's method and the current method preserve the velocity. Thus, of these second-order relativistic momentum integrations, only the integrator introduced here both preserves volume and gives the correct velocity. While all methods have error that is second-order in time, they deviate from each other by terms that increase as the motion becomes relativistic. Numerical results show that Vay's method develops energy errors near resonant orbits of a test problem that neither volume-preserving integrator does.
Computing trajectories of classical charged particles in an electromagnetic field is critical to many areas of physics. In beam physics (cf. Ref. 1 and references therein), particle tracking through accelerator lattices is needed to determine whether particles will remain in, e.g., storage rings for times long enough for significant reactions to occur. In self-consistent plasma physics computations, the computation of particle trajectories is one part of kinetic modeling of plasmas using the Particle-In-Cell (PIC) method.2–4 Particle tracking is also important for understanding the limits of precision in spectrometers.
The underlying system is Hamiltonian, and so has a complete set of Poincaré geometric invariants.5 Symplectic integrators are those that also preserve the Poincaré invariants. Symplectic integrators allow one to invoke the Kolmogorov-Arnold-Moser (KAM) theorem, so that the integrated motion shadows the real motion (except for the very slow Arnold diffusion in systems with 2.5 or more degrees of freedom). Explicit symplectic integrators6–10 relevant to particle tracking are now in common use, and there are also explicit symplectic integrators11 for certain plasma systems. Unfortunately, there is at present no known explicit symplectic integrator for charged particle motion in arbitrary electromagnetic fields.
A less restrictive condition is that the integrator be volume-preserving, i.e., that it preserves the last Poincaré invariant. This prevents attractors and repellers in the integrated system, important as they do not exist in the underlying system. The relativistic Boris integrator12 (also known as the Boris push) is volume-preserving (more detailed discussion in light of Refs. 13 and 14 below) and is found to have good properties when used to compute trajectories. In particular, the spatial Boris push15 has been found excellent for tracking studies in accelerators, as it eliminates numerical cooling or heating, which can mask physical effects due to small dissipative terms.
On the other hand, the relativistic Boris push does not correctly compute the drift velocity, and this was found16 to lead to problems in some special cases. Vay proposed using a different integrator that does preserve the drift velocity at the cost of a slightly more complicated computation of the new relativistic factor. Unfortunately, Vay's integrator is not volume-preserving, as we show here, and so it lacks an important feature of the Boris integrator. The lack of volume-preservation is shown to lead to unphysical behavior in the computed trajectories.
Here, we construct an integrator that preserves both the drift velocity and phase-space volume. Like both the Boris and Vay integrators, it also conserves energy in the absence of an electric field. This new integrator has computational cost comparable to that of the Vay integrator.
The organization of this paper is as follows. In Sec. II, we introduce the leap-frog separation of the spatial and momentum integrations that reduce this to second-order integration of the momentum change equation. In the same section, we discuss also the various centerings that can be used to obtain a second-order-accurate momentum integration, deriving a new integrator and showing that the previous Boris and Vay integrators follow from this principle and have the same order of accuracy. Section III discusses implementation. In Sec. IV, we show that all three integrators conserve energy when there is no electric field, but that only the current integrator and the integrator of Ref. 16 properly capture the drift. In Sec. V, we show that only the current and Boris integrators are volume-preserving. In Sec. VI, we show that the new and Boris integrators eliminate un-physical behavior displayed by the integrator of Ref. 16 in a test problem. Finally, in Sec. VII, we summarize our results and conclude.
II. SECOND-ORDER CHARGED PARTICLE INTEGRATORS
Charged particles in the electromagnetic field follow the dynamics of relativistic translation:
where , and the Lorentz force
The first step to an efficient, second-order integration of these equations is use of the Leap-Frog method, in which the finite-time-step relativistic translation becomes
and then one is left with solving the Lorentz force equation Eq. (2) to second order for a time step, , with the assumption of constant electric and magnetic fields.
Ideally one would like to solve this equation in a way that preserves as many properties of the underlying differential equations as possible. In this paper, we consider the following properties: (1) Energy should be conserved in the absence of an electric field. (2) The static solution for crossed electric and magnetic fields, with , should be constant velocity in the third direction of magnitude . (3) The differential volume, which is preserved by any solution of the differential equation, should be preserved by the finite-time-step solution.
The standard way to obtain a second-order solution is by time centering. That is, one uses the solution
where is an average of the initial and final values of . There are multiple choices for how to do the average. Here, we introduce a new choice
while16 made the choice
and the Boris choice is
The Boris push is not usually written in this way, but one can see that Eq. (9) is equivalent to a magnetic rotation by the velocity found after an initial half electric acceleration (the first term in the numerator) and before a final half acceleration (the second term in the numerator).
All of these integrators have second-order error in the time step, as one can see from Taylor expansion. For example,
is equivalent to Eq. (7) to second order. That the Boris push has the same order of error follows from a similar calculation.
III. EXPLICIT EVALUATION
The new integrator (5)–(7) looks implicit, as the final momentum on the left side of (5) is involved in its definition through (7). However, it can be explicitly computed by methods similar to that of Ref. 16. We first write the new integrator as the composition of two integrators
which can be solved to obtain
The new integrator's implementation differs from the Boris integrator's in the use of Eq. (20) instead of Eq. (21) to calculate γ for the magnetic rotation. The two prescriptions differ by terms second-order in . One can furthermore show that the relativistic factor used by the Boris integrator always exceeds or equals that of the present integrator
IV. PRESERVATION OF LIMITING SOLUTIONS
Here, we analyze the above integration methods with regard to whether they preserve properties of limiting solutions, (1) no electric field and (2) correct value of perpendicular velocity in crossed electric and magnetic fields. In the first case, the differential equations have no change of energy. This is well-known to be true for the Boris push. It follows that it is true for the other integrators, as, without an electric field, the relativistic factor is preserved even for the finite difference solution, and so both the new integrator (7) and the Vay integrator (8) preserve energy when there is no electric field.
When there is no electric field and the magnetic field is static and uniform, the Boris integrator rotates the particle's momentum through an angle
It can be shown that, for the same problem, the new integrator rotates the particle's momentum through an angle
In the non-relativistic case, Eq. (27) reduces to Eq. (26), with the coefficient of the third-order error term increasing monotonically to a limit of in the ultra-relativistic case (vanishing at ). The new integrator's third-order error term is therefore always smaller than the Boris integrator's for electron gyro-motion with .
For perpendicular electric and magnetic fields, a stationary solution, , must have, from Eq. (6)
Hence for either the new integrator or the Vay integrator, the electric drift velocity has the same value as for the differential equations. For the Boris integrator, as noted in Ref. 16, this is not the case. Inequality (25) explains why: the Boris integrator always rotates the momentum less than the new integrator, which is shown here to rotate exactly enough that the second half-acceleration exactly cancels the first.
It has been surmised13 that the excellent conservation properties of the nonrelativistic Boris push, as, e.g., observed in Ref. 15, are due to the fact that it preserves phase space volume, like the underlying differential equations. This follows straightforwardly from the observation that the nonrelativistic Boris push is a sequence of sheared (spatial change) translations, unsheared (electric acceleration) translations, and rotations (magnetic acceleration). An additional observation, that the relativistic Boris push simply changes the rotation to being sheared (due to the dependence of the rotation on the relativistic factor, γ, which is constant for that part of the transformation) shows that the relativistic Boris push is also volume-preserving. (Reference 14 shows that an integrator that is not actually the Boris integrator does not preserve volume but then ultimately gives a more detailed proof of volume preservation by what is in fact the Boris integrator.)
While the two transformations (12) and (13) are not inverses of each other (otherwise the integrator would be the identity), we will show that their two Jacobian determinants are reciprocals of one another and so the transformation for the new integrator is volume-preserving.
Differentiating (12) gives the Jacobi matrix
where is the matrix defined by
By the determinant lemma, this matrix has determinant
Explicit computation gives
Explicit computation can also be used to find
by noting that
can be solved to give
By an identical process, one can show that
Thus, the Jacobian of the first half step equals the inverse of the Jacobian for the second half step, and so their product is unity, and the new integrator is volume-preserving, just like the Boris integrator.
The Vay integrator can be analyzed in this same way. It can likewise be composed into two half steps, but with the first being an explicit step using followed by an implicit step using . As a result, the full Jacobian for the Vay integrator is:
Consequently, after N steps, the differential volume element in the Vay integrator is
Because the two Jacobians in any of the fractions of Eq. (42) depend on different variables (the initial and final momenta), while having the same functional form, the Vay integrator is not generally volume-preserving. At subsequent steps, the differential volume grows or shrinks by a similar factor, but now evaluated at the new position. The function, J(x, u) varies over space but is bounded for bounded regions. Hence, generally the factors in Eq. (42) are variously greater and less than unity. While one could imagine a trajectory that conspires to have all factors either greater or less than unity, we have not so far been able to construct such a case. In the case where the magnetic field is constant in space and time, the series telescopes and then, the boundedness of prevents the existence of attractors or repellers.
VI. NUMERICAL RESULTS
We choose a test problem of the form
The Hamiltonian for this system is
in units of This Hamiltonian is independent of z, so pz is an invariant of the motion. It is also time independent, so H is an invariant. Further
is an invariant of the motion and in involution with the other invariants. Hence, this system of three degrees of freedom has three invariants in involution and so is integrable.
To determine the degree to which the various integrators preserve these invariants, we use the Poincaré surface of section technique. We follow the trajectories for initial conditions of the same energy and same pz and plot the points in the plane x = 0 when that plane is crossed with positive px. (This does require interpolation to that plane, when that plane is crossed in some time step.) We can then check how well the invariant (47) is preserved by the integration: if all of the other invariants associated with the problem are fixed, then the surface of section should consist of non-intersecting curves.
For numerical integration, we have to be more specific. We choose pz = 0 without loss of generality, as a different value of pz is equivalent to choosing a different function . We then choose the particular fields
which corresponds to a sheared magnetic field with the reversal at y = 0. For this system
and the surface of section should show the family of curves given by
We choose the energy to be moderately relativistic to expose the effects of the varying volume element (40), and we choose a = 1 and b = 2, i.e., of order unity.
We first use all integrators with a very small time step, , which corresponds to time steps per period of the oscillation in x for the given potential. The results are shown in Fig. 1. At this temporal resolution, all integrators are seen to produce non-intersecting curves in this plane, in essence showing that they are all preserving the invariants of the problem.
Next we use a time step, that is more typically used in simulations. This corresponds to time steps per period of the oscillation in x for the given potential. The results are shown in Fig. 2. At this temporal resolution, the Boris and new integrators continue to show non-intersecting curves in this plane, showing that they are accurately representing the topology of the trajectories. However, for initial conditions near , where the periods of the x and y oscillations are equal and the overall motion periodic, the Vay integrator exhibits resonance islands: closed curves confined to the upper and lower half-planes, clearly not belonging to the family of curves given by Eq. (52). Moreover, the sections for different trajectories are seen to cross each other. This is an indication that the other invariants (either energy or pz) are less well-preserved, as crossing can occur only if the trajectories are for different values of the other invariants. This is demonstrated in Figure 3, which plots the error in the energy for initial conditions near the border of the region where the Vay integrator's surface-of-section curves are confined to the upper and lower half-planes. The likely cause of this phenomenon is the near-periodicity of these orbits—the Vay integrator's errors accumulate constructively over several similar orbits.
These are the orbits where the Vay integrator exhibits the largest error; it exhibits smaller errors for initial conditions far from the points . These errors are nevertheless bounded and periodic. Figure 4 shows in more detail how the errors in energy and Iy vary for orbits with different initial values of in the neighborhood of the problematic value . The spread of the energy errors (Fig. 4(a)) in the solutions of the Boris and new integrators is smaller than that of the Vay integrator's solution by a factor where the Vay integrator's errors are largest. The errors in Iy have more comparable spreads, but the Boris and new integrators' mean errors in Iy (Fig. 4(b)) are smaller by factors reaching 101 and 102, respectively.
Thus, in this example, with a reasonable choice for the time step, use of the Vay integration leads to unphysical consequences.
VII. SUMMARY AND FUTURE DIRECTIONS
We have derived a new integrator for charged particle motion in arbitrary electromagnetic fields that is both volume-preserving (like the Boris push) and also correctly computes the drift velocity (like the Vay push). This new integrator has been tested numerically and compared with the other integrators. It is found not to introduce sizable resonances at reasonable values for the time step in contrast with the integrator of Ref. 16.
A number of new directions deserve attention in this area. What are the consequences of composing integrators for which the volume alternately grows and shrinks, as it does for the Vay integrator? Can one have attractors and repellers in this case? How does one extend these methods to higher order, as was done for symplectic integrators by, e.g., Ref. 10 and others. More importantly, what is gained? For beam tracking, one needs a spatial version of the current integrator. Combining higher-order with spatial integration could be very powerful for beam tracking. Finally, the extension of these concepts to self-consistent simulations, with higher-order integrators, could be explored.
This work was supported by DOE/NSF Grant No. DE-SC0012584.