Structure-preserving second-order integration of relativistic charged particle trajectories in electromagnetic fields

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.⁵ 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 integrators^6–10 relevant to particle tracking are now in common use, and there are also explicit symplectic integrators¹¹ 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 integrator¹² (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 push¹⁵ 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 $E→×B→$ drift velocity, and this was found¹⁶ to lead to problems in some special cases. Vay proposed using a different integrator that does preserve the $E→×B→$ 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 $E→×B→$ 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 $E→×B→$ 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.

Charged particles in the electromagnetic field follow the dynamics of relativistic translation:

where $γ=1+|u→/c|2$⁠, 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, $Δt$⁠, 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 $|E→|<c|B→|$⁠, should be constant velocity in the third direction of magnitude $|E→|/|B→|$⁠. (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 $v¯$ is an average of the initial and final values of $v→$⁠. There are multiple choices for how to do the average. Here, we introduce a new choice

while¹⁶ made the choice

and the Boris choice is

where

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.

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

where

with

and $γnew≡γ(u→new)$⁠. The first equation (12) corresponding to the last half update is explicit. So we need only solve Eqs. (14)–(16) to obtain an explicit solution.

In a process similar to that of Ref. 16, the explicit solution comes from taking the dot and cross products of Eq. (14) and clearing terms to obtain

in which

Eq. (17) gives $u¯new$ explicitly provided γ_new can be computed. Squaring Eq. (17) gives the biquadratic polynomial

which can be solved to obtain

where

Equation (17) is needed to derive Eq. (20) for γ_new, but, once γ_new is known, $u¯new$ is not necessary to obtain $u→f$⁠. Using the definition of $Δu→new$ and Eqs. (12)–(14)

Defining then and $u→+=u→f−ϵ→$⁠, using Eq. (18), and subtracting Eq. (22) from Eq. (23) yields the familiar Boris rotation equation

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 $ΩcΔt$⁠. One can furthermore show that the relativistic factor used by the Boris integrator always exceeds or equals that of the present integrator

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 $+1/24$ in the ultra-relativistic case (vanishing at $γ=3$⁠). The new integrator's third-order error term is therefore always smaller than the Boris integrator's for electron gyro-motion with $ΩcΔt≪1$⁠.

For perpendicular electric and magnetic fields, a stationary solution, $Δv→=0$⁠, must have, from Eq. (6)

Stationarity also implies $u→f=u→i$⁠, which is used in either of Eqs. (7) and (8) implies that

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 $u→−$ exactly enough that the second half-acceleration exactly cancels the first.

It has been surmised¹³ 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

which can be used in Eq. (34) and ultimately with Eqs. (32) and (33) to obtain

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 $ui/γi$ followed by an implicit step using $uf/γf$⁠. As a result, the full Jacobian for the Vay integrator is:

where

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 $J(x,u)$ prevents the existence of attractors or repellers.

We choose a test problem of the form

for which

The Hamiltonian for this system is

in units of $q=m=c=1$ This Hamiltonian is independent of z, so p_z 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 p_z and plot the points $(y,py)$ in the plane x = 0 when that plane is crossed with positive p_x. (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 p_z = 0 without loss of generality, as a different value of p_z is equivalent to choosing a different function $Az(y)$⁠. 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 $H=4.$ 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, $Δt=1/40$⁠, which corresponds to $80π≈250$ 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, $Δt=1/10$ that is more typically used in simulations. This corresponds to $20π≈60$ 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 $py≈1.7$⁠, 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 p_z) 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 $δH(t)=γ+12x(t)2−4$ 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 $(y,py)=(0,±1.7)$⁠. These errors are nevertheless bounded and periodic. Figure 4 shows in more detail how the errors in energy and I_y vary for orbits with different initial values of $py∝Iy1/2$ in the neighborhood of the problematic value $py,0≈1.7$⁠. 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 $≈102$ where the Vay integrator's errors are largest. The errors in I_y have more comparable spreads, but the Boris and new integrators' mean errors in I_y (Fig. 4(b)) are smaller by factors reaching 10¹ and 10², respectively.

Thus, in this example, with a reasonable choice for the time step, use of the Vay integration leads to unphysical consequences.

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 $E→×B→$ 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.