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

TL;DR

This paper introduces a new second-order, structure-preserving integration method for relativistic charged particle trajectories in electromagnetic fields, improving energy conservation and velocity accuracy over existing methods.

Contribution

A novel second-order integration method that preserves volume and the correct  imes \u00B5 B velocity, differing from the relativistic Boris Push by its Lorentz factor calculation.

Findings

The new method conserves energy in absence of electric fields.

It preserves volume and  imes B5 B velocity, unlike some existing methods.

Numerical tests show improved energy stability near resonant orbits.

Abstract

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 [1] only in the method for calculating the Lorentz factor. This is shown to have the same second-order accuracy in time as that (Boris Push) [1] found by splitting the electric acceleration and magnetic rotation and that [2] found by averaging the velocity in the magnetic rotation term. All three methods are shown to conserve energy when there is no electric field. The Boris method and the current method are shown to be volume-preserving, while the method of [2] and the current method preserve the $\vec{E} \times \vec{B}$ velocity. Thus, of these second-order relativistic momentum integrations, only the integrator introduced here both preserves volume and gives the correct $\vec{E} \times \vec{B}$ 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 [2] develops energy errors near resonant orbits of a test problem that neither volume-preserving integrator does. [1] J. Boris, Relativistic plasma simulation-optimization of a hybrid code, in: Proc. Fourth Conf. Num. Sim. Plasmas, Naval Res. Lab, Wash. DC, 1970, pp. 3-67. [2] J.-L. Vay, Simulation of beams or plasmas crossing at relativistic velocity, Physics of Plasmas (1994-present) 15 (5) (2008) 056701.

Equations

$$\frac{d\vec{x}}{dt}=\vec{v}(\vec{u})=\vec{u}/\gamma ,$$

$$\frac{d\vec{u}}{dt}=(q/m)(\vec{E}+\vec{v}(\vec{u})\times \vec{B}).$$

$${\vec{x}}_f$$

$$\Delta \vec{x}$$

$${\vec{u}}_f$$

$$\Delta \vec{u}$$

$${\bar{v}}_{new}\equiv \vec{v}\left(\frac{{\vec{u}}_i+{\vec{u}}_f}{2}\right).$$

$${\bar{v}}_v\equiv \frac{\vec{v}({\vec{u}}_i)+\vec{v}({\vec{u}}_f)}{2}$$

$${\bar{v}}_b\equiv \frac{\vec{v}\left(u_i+\overrightarrow{ϵ}\right)+\vec{v}\left(u_f-\overrightarrow{ϵ}\right)}{2}.$$

$$\overrightarrow{ϵ}\equiv \frac{q\vec{E}}{2m}\Delta t.$$

$${\bar{v}}_v$$

$${\vec{u}}_f$$

$${\vec{u}}_{new}$$

$${\vec{u}}_i$$

$$\Delta {\vec{u}}_{new}=2\overrightarrow{ϵ}+\frac{{\bar{u}}_{new}}{{\gamma }_{new}}\times 2\vec{\beta },$$

$$\vec{\beta }\equiv \frac{q\vec{B}}{2m}\Delta t,$$

$${\bar{u}}_{new}\left(1+\frac{{\beta }^2}{{\gamma }_{new}^2}\right)={\vec{u}}_{-}-\frac{\vec{\beta }\times {\vec{u}}_{-}}{{\gamma }_{new}}+\frac{\vec{\beta }\vec{\beta }\cdot {\vec{u}}_{-}}{{\gamma }_{new}^2},$$

$${\vec{u}}_{-}\equiv {\vec{u}}_i+\overrightarrow{ϵ}.$$

$$({\gamma }_{new}^2-1)\left({\gamma }_{new}^2+{\beta }^2\right)={\gamma }_{new}^2|{\vec{u}}_{-}{|}^2+|\vec{\beta }\cdot {\vec{u}}_{-}{|}^2$$

$${\gamma }_{new}^2=\frac{1}{2}\left({\gamma }_{-}^2-{\beta }^2+\sqrt{({\gamma }_{-}^2-{\beta }^2{)}^2+4({\beta }^2+|\vec{\beta }\cdot {\vec{u}}_{-}^2|)}\right),$$

$${\gamma }_{-}\equiv \gamma ({\vec{u}}_{-}).$$

$${\vec{u}}_i+\overrightarrow{ϵ}$$

$${\vec{u}}_f-\overrightarrow{ϵ}$$

$${\vec{u}}_{+}-{\vec{u}}_{-}=({\vec{u}}_{+}+{\vec{u}}_{-})\times \frac{\vec{\beta }}{{\gamma }_{new}}$$

$${\gamma }_{-}\ge {\gamma }_{new}$$

$$\theta ={\Omega }_c{\gamma }^{-1}\Delta t\left(1-\frac{({\Omega }_c{\gamma }^{-1}\Delta t{)}^2}{12}+\mathcal{O}(({\Omega }_c{\gamma }^{-1}\Delta t{)}^4)\right)$$

$$\theta ={\Omega }_c{\gamma }^{-1}\Delta t\left(1+\left[\frac{1}{8}\left(1-\frac{1}{{\gamma }^2}\right)-\frac{1}{12}\right]({\Omega }_c{\gamma }^{-1}\Delta t{)}^2+\mathcal{O}(({\Omega }_c{\gamma }^{-1}\Delta t{)}^4)\right).$$

$$\bar{v}=\vec{E}\times \vec{B}/|\vec{B}{|}^2.$$

$${\bar{v}}_f={\bar{v}}_i=\vec{E}\times \vec{B}/|\vec{B}{|}^2.$$

$$\frac{\partial ({\vec{u}}_f)}{\partial ({\bar{u}}_{new})}=\mathbf{I}-\mathbf{\Omega }+\frac{{\vec{\beta }}_{new}\times {\bar{u}}_{new}\otimes {\bar{u}}_{new}}{{\gamma }_{new}^3},$$

$$\mathbf{\Omega }\cdot \vec{V}\equiv \vec{\beta }\times \vec{V}/{\gamma }_{new}.$$

$$J_{f,new}\equiv \det \left(\frac{\partial ({\vec{u}}_f)}{\partial ({\bar{u}}_{new})}\right)=\det (\mathbf{I}-\mathbf{\Omega })\det \left(1+\frac{{\bar{u}}_{new}\cdot (\mathbf{I}-\mathbf{\Omega }{)}^{-1}\cdot (\vec{\beta }\times {\bar{u}}_{new})}{{\gamma }_{new}^3}\right).$$

$$\det (\mathbf{I}-\mathbf{\Omega })=1+{\beta }^2/{\gamma }_{new}^2.$$

Full text

I Introduction

Computing trajectories of classical charged particles in an electromagnetic field is critical to many areas of physics. In beam physics (cf cary2016select and references therein) particle tracking through accelerator lattices is needed to determine whether particle 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 birdsall2004plasma ; hockney1988computer ; nieter2004vorpal . 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 goldstein2014classical . Symplectic integrators are those that also preserve the Poincaré invariants. Symplectic integrators allow one to invoke the 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 ruth1983canonical ; feng1986difference ; forest1990fourth ; candy1991symplectic ; yoshida1990construction relevant to particle tracking are now in common use, and there are also explicit symplectic integrators cary1993explicit 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 preserve 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 boris1970relativistic (also known as the Boris push) is volume-preserving (more detailed discussion in light of qin2013boris ; zhang2015volume below) and is found to have good properties when used to compute trajectories. In particular, the spatial Boris push stoltz2002efficiency 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 ${\vec{E}\times\vec{B}}$ drift velocity, and this was found vay2008simulation to lead to problems in some special cases. Vay proposed using a different integrator that does preserve the ${\vec{E}\times\vec{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 ${\vec{E}\times\vec{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 the following section 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 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 vay2008simulation properly capture the ${\vec{E}\times\vec{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 vay2008simulation 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,

[TABLE]

where $\gamma=\sqrt{1+|\vec{u}/c|^{2}}$, and the Lorentz force,

[TABLE]

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

[TABLE]

and then one is left with solving the Lorentz force equation Eq. (2) to second order for a time step, $\Delta 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 $|\vec{E}|<c|\vec{B}|$, should be constant velocity in the third direction of magnitude $|\vec{E}|/|\vec{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,

[TABLE]

where $\bar{v}$ is an average of the initial and final values of $\vec{v}$. There are multiple choices for how to do the average. Here we introduce a new choice,

[TABLE]

while vay2008simulation made the choice,

[TABLE]

and the Boris choice is,

[TABLE]

where

[TABLE]

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. E.g.,

[TABLE]

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 vay2008simulation . We first write the new integrator as the composition of two integrators.

[TABLE]

where

[TABLE]

with

[TABLE]

and $\gamma_{new}\equiv\gamma(\vec{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 vay2008simulation , the explicit solution comes from taking the dot and cross products of Eq. (14) and clearing terms to obtain

[TABLE]

in which

[TABLE]

Eq. (17) gives $\bar{u}_{new}$ explicitly provided $\gamma_{new}$ can be computed. Squaring Eq. (17) gives the biquadratic polynomial,

[TABLE]

which can be solved to obtain

[TABLE]

where

[TABLE]

Equation (17) is needed to derive Eq. (20) for $\gamma_{new}$, but, once $\gamma_{new}$ is known, $\bar{u}_{new}$ is not necessary to obtain $\vec{u}_{f}$. Using the definition of $\Delta\vec{u}_{new}$ and Eqs. (12-14)

[TABLE]

Defining then and $\vec{u}_{+}=\vec{u}_{f}-\vec{\epsilon}$, using Eq. (18), and subtracting Eq. (22) from Eq. (23) yields the familiar Boris rotation equation:

[TABLE]

The new integrator’s implementation differs from the Boris integrator’s in the use of Eq. (20) instead of Eq. (21) to calculate $\gamma$ for the magnetic rotation. The two prescriptions differ by terms second-order in $\Omega_{c}\Delta t$. One can furthermore show that the relativistic factor used by the Boris integrator always exceeds or equals that of the present integrator.

[TABLE]

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

[TABLE]

It can be shown that, for the same problem, the new integrator rotates the particle’s momentum through an angle

[TABLE]

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 $\gamma=\sqrt{3}$). The new integrator’s third-order error term is therefore always smaller than the Boris integrator’s for electron gyro-motion with $\Omega_{c}\Delta t\ll 1$.

For perpendicular electric and magnetic fields, a stationary solution, $\Delta\vec{v}=0$, must have, from Eq. (6),

[TABLE]

Stationarity also implies $\vec{u}_{f}=\vec{u}_{i}$, which used in either of Eqs. (7,8) implies that

[TABLE]

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 vay2008simulation , 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 $\vec{u}_{-}$ exactly enough that the second half-acceleration exactly cancels the first.

V Volume-preservation

It has been surmised qin2013boris that the excellent conservation properties of the nonrelativistic Boris push, as, e.g., observed in stoltz2002efficiency , 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, $\gamma$, which is constant for that part of the transformation) shows that the relativistic Boris push is also volume-preserving. (The reference zhang2015volume 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,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,

[TABLE]

where $\bm{\Omega}$ is the matrix defined by

[TABLE]

By the determinant lemma, this matrix has determinant,

[TABLE]

Explicit computation gives

[TABLE]

Explicit computation can also be used to find

[TABLE]

by noting that

[TABLE]

can be solved to give

[TABLE]

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

[TABLE]

By an identical process, one can show that

[TABLE]

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 $u_{i}/\gamma_{i}$ followed by an implicit step using $u_{f}/\gamma_{f}$. As a result the full Jacobian for the Vay integrator is

[TABLE]

where

[TABLE]

Consequently, after N steps, the differential volume element in the Vay integrator is

[TABLE]

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.

VI Numerical results

We choose a test problem of the form,

[TABLE]

for which

[TABLE]

The Hamiltonian for this system is

[TABLE]

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,

[TABLE]

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,p_{y})$ 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.

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 $A_{z}(y)$. We then choose the particular fields,

[TABLE]

which corresponds to a sheared magnetic field with the reversal at $y=0$. For this system,

[TABLE]

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, $\Delta t=1/40$, which corresponds to $80\pi\approx 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 give nested surfaces in this plane, in essence showing that they are all preserving the invariants of the problem.

Next we use a time step, $\Delta t=1/10$ that is more typically used in simulations. This corresponds to $20\pi\approx 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 nested surfaces in this plane, showing that they are accurately representing the topology of the trajectories. However, for the Vay integrator with the initial condition $p_{y}\approx 1.7$ a resonance island is seen. (This is a two-fold degenerate resonance; trajectories started in the top island do not visit the bottom island.) Thus, a finite volume of trajectories is trapped in this resonance. 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. 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 ${\vec{E}\times\vec{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 vay2008simulation .

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., yoshida1990construction 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.

Acknowledgements

This work is supported by DOE/NSF Grant No. DE-SC0012584.