High Order Numerical Integrators for Relativistic Charged Particle Tracking

TL;DR

This paper develops higher order, time-reversible numerical integrators for relativistic charged particle tracking, improving accuracy and computational efficiency in electromagnetic field simulations.

Contribution

The paper introduces a fourth order algorithm for relativistic particle tracking, extending existing integrators to higher order accuracy with explicit testing and validation.

Findings

Fourth order integrator achieves higher accuracy.

High order integrators reduce computational cost.

Explicit numerical examples demonstrate effectiveness.

Abstract

In this paper, we extend several time reversible numerical integrators to solve the Lorentz force equations from second order accuracy to higher order accuracy for relativistic charged particle tracking in electromagnetic fields. A fourth order algorithm is given explicitly and tested with numerical examples. Such high order numerical integrators can significantly save the computational cost by using a larger step size in comparison to the second order integrators.

Equations

$$\frac{d\mathbf{r}}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

$$\frac{dt}{ds}$$

$$\frac{d\mathbf{r}}{ds}$$

$$\frac{d\mathbf{p}}{ds}$$

$$\frac{d\zeta }{ds}$$

$$A$$

$$\zeta (\tau )$$

$$B$$

$$C$$

$$\zeta (\tau )$$

$$t(\tau /2)$$

$$\mathbf{r}(\tau /2)$$

$${\mathbf{p}}_{-}$$

$${\gamma }_{-}$$

$${\mathbf{p}}_{+}-{\mathbf{p}}_{-}$$

$$\mathbf{p}(\tau )$$

$${\gamma }_0$$

$${\mathbf{p}}_{-}$$

$${\mathbf{p}}_{+}$$

$${\gamma }_1$$

$$\mathbf{t}$$

$$\lambda$$

$$\sigma$$

$${\gamma }_2$$

$${\mathbf{t}}^{\ast }$$

$$s$$

$$\mathbf{p}(\tau )$$

$${\mathbf{p}}_{-}$$

$${\gamma }_{-}$$

$$\mathbf{t}$$

$${\gamma }_{new}$$

$${\mathcal{M}}_{2nd}(\tau )$$

$${\mathcal{M}}_{4th}(\tau )$$

$${\mathcal{M}}_{2n+2}(\tau )$$

$$E_x$$

$$E_y$$

$$E_z$$

$$B_x$$

$$B_y$$

$$B_z$$

$$E_x$$

Full text

I Introduction

Numerical tracking charged particle in electric and magnetic fields has many applications in beam physics and plasma physics. It normally involves solving the Lorentz force equations numerically with external electromagnetic fields. In previous studies, a second order, time reversible numerical algorithm, known as Boris integrator boris has been widely used in plasma and beam physics simulations. However, for numerical simulation of relativistic charged particle motion in electromagnetic fields, this integrator can produce large error vay ; qiangx ; qiangy . A new time reversible second order integrator was proposed in reference vay that avoids this problem. Recently, another time reversible second order integrator was proposed by Higuera and Cary hc . Besides working well for particle tracking in electromagnetic fields with large relativistic factor, this algorithm also preserves phase space volume.

The above three time reversible numerical integrators are second order accuracy of integration step size. In some numerical simulations, a higher order numerical integrator can be more effective in attaining the desired numerical accuracy. So far, the extension of these second order integrators to higher order accuracies have not been reported in literature. Meanwhile, in the area of symplectic numerical integrator study of Hamiltonian systems, high order numerical integrators have been reported by using a split-operator method forest1 ; yoshida . In this paper, after reformating the original Lorentz force equations, we observed that higher order numerical integrators can be obtained from these symmetric, time reversible second order integrators for relativistic charge particle tracking.

The organization of this paper is as follows: after the introduction, we present the high order numerical integrator in section II; We present numerical tests of the fourth order integrator in Section III and draw conclusions in Section IV.

II High Order Numerical Integrators

The Lorentz equations of motion for a charged particle subject to electric and magnetic fields can be written as:

[TABLE]

where ${\bf r}=(x,y,z)$ denotes the particle spatial coordinates, ${\bf p}=(p_{x}/mc,p_{y}/mc,p_{z}/mc)$ the particle normalized mechanic momentum, $m$ the particle rest mass, $q$ the particle charge, $c$ the speed of light in vacuum, $\gamma$ the relativistic factor defined by $\sqrt{1+{\bf p}\cdot{\bf p}}$, $t$ the time, ${\bf E}(x,y,z,t)$ the electric field, and ${\bf B}(x,y,z,t)$ the magnetic field. Instead of using the time $t$ as an explicit independent variable, we write the above equations using $s$ as independent variable:

[TABLE]

Letting $\zeta(t,{\bf r},{\bf p}:s)$ denote a vector of coordinates, the above equations of motion can be rewritten as:

[TABLE]

where the matrix $A$ is a given as:

[TABLE]

A formal solution for above equation after a single step $\tau$ can be written as:

[TABLE]

The matrix $A$ can be written as a sum of two terms $A=B+C$, where

[TABLE]

and

[TABLE]

Using the Baker-Campbell-Hausdorff theorem baker ; campbell ; haus , a second order approximation for above single step solution can be obtained as:

[TABLE]

Letting $\exp(\frac{1}{2}\tau B)$ define a transfer map ${\mathcal{M}}_{1}$ and $\exp(\tau C)$ a transfer map ${\mathcal{M}}_{2}$, for a single step, the above splitting results in a second order numerical integrator for the original equation as:

[TABLE]

From definitions of the matrices $B$ and $C$, it is seen that the transfer map ${\mathcal{M}}_{1}$ corresponds to the solutions of Eqs. 3 and 4 for half step, and transfer map ${\mathcal{M}}_{2}$ corresponds to the solution of Eq. 5 for one step. The solutions of transfer map ${\mathcal{M}}_{1}(\tau/2)$ is straightforward and can be written as:

[TABLE]

The solution for ${\mathcal{M}}_{2}(\tau)$ can have different forms depending on different ways of approximation. In the Boris algorithm, ${\mathcal{M}}_{2}(\tau)$ is given as:

[TABLE]

where ${\bf p}_{+}$ can be solved analytically from the linear equation Eq. 17. The Boris algorithm is time reversible and has been widely used in numerical plasma and beam physics simulations. However, it was found that the Boris algorithm could have large numerical error for charged particle tracking when the particle relativistic factor is largevay,qiangx,qiangy. The source of this numerical error might result from the momentum update in separate steps from the electric field and from the magnetic field. This becomes especially a serious problem to simulate a relativistic charged particle beam including space-charge effects, where the electric field and the magnetic field cancel each other significantly in the laboratory frame and results in $1/\gamma^{2}$ decrease of the transverse space-charge effects. The new time reversible solution that avoids this problem in ${\mathcal{M}}_{2}(\tau)$ was proposed in reference vay as:

[TABLE]

This algorithm does not have the problem of the Boris algorithm and works well for relativistic particle tracking. Recently, another time reversible and structure-preserving algorithm for ${\mathcal{M}}_{2}(\tau)$ was proposed in reference hc . This algorithm is similar to the Boris algorithm by replacing the $\gamma_{-}$ in Eq. 17 with the following $\gamma_{new}$:

[TABLE]

This algorithm also works well for particle tracking with large relativistic factor and preserves phase space volume.

So far, all these three time reversible algorithms have a second order accuracy of integration step size. The second order one step map can be rewritten as:

[TABLE]

Since this one-step map is symmetric and time reversible, one can follow the exact steps of the reference yoshida and construct the same fourth order accuracy numerical integrator as:

[TABLE]

where $\alpha=1-2^{1/3}$, and $s=\tau/(1+\alpha)$. An arbitrary even order accuracy integrator can also be obtained following that reference. Assume that ${\mathcal{M}}_{2n}$ denotes a transfer map with an accuracy of order $2n$, the transfer map ${\mathcal{M}}_{2n+2}$ with $(2n+2)$th order of accuracy can be obtained from the recursion equation yoshida :

[TABLE]

where $z_{0}=1/(2-2^{1/(2n+1)})$ and $z_{1}=-2^{1/(2n+1)}/(2-2^{1/(2n+1)})$.

III Numerical Tests

We tested the above $4^{th}$ order extension of the Boris algorithm, the Vay algorithm, and the Higuera-Cary algorithm using two numerical examples. In the first example, we considered an electron moving inside static electric and magnetic fields. These fields are given as:

[TABLE]

where $\gamma_{0}$ is the relativistic factor of the moving beam, $\beta_{0}=\sqrt{1-(1/\gamma_{0})^{2}}$, and the constant $E_{0}=9\times 10^{6}$. The above external fields correspond to the space-space fields generated by a moving infinitely long transverse uniform cylindrical positron beam. First, we assume that both the initial electron kinetic energy and the moving positron beam kinetic energy are $2$ MeV. Figure 1 shows the electron trajectory evolution as a function of time from the $4^{th}$ order extension of the Boris integrator (magenta), the Vay integrator (green), and the Higuera-Cary integrator (blue) with a step size of $0.25$ ns (around $0.01$ oscillation period).

It is seen that three numerical integrators agree with each other very well in this case.

Figure 2 shows the relative numerical errors at the end of integration as a function of step size from the $4^{th}$ order extension of the Boris integrator (magenta), the Vay integrator (green), and the Higuera-Cary integrator (blue) together with a power $4$ polynomial. It is seen that all three numerical integrators have nearly the same relative errors and converge as $4^{th}$ power with respect to the step size. Next, we assumed that both the initial electron and the moving positron beam have a kinetic energy of $50$ MeV. Figure 3 shows the electron trajectory evolution as a function of time from the $4^{th}$ order extension of the Boris integrator (magenta), the Vay integrator (green), and the Higuera-Cary integrator (blue) with a step size of $4$ ns (around $0.008$ oscillation period).

It is seen that at beginning all three integrators agree with each other well. After $2$ ns, the Boris integrator starts to deviate from the other two numerical integrators due to the cancellation errors from the electric field and the magnetic field. Figure 4 shows the relative numerical errors at the end of integration as a function of step size from the $4^{th}$ order extension of the Boris integrator (magenta), the Vay integrator (green), and the Higuera-Cary integrator (blue) together with a power $4$ polynomial. It is seen that all three numerical integrators converge as $4^{th}$ power of the step size. However, the $4^{th}$ order extension of the Boris integrator still shows much larger relative errors than the other two $4^{th}$ order integrators.

In the second example, we assume that a $10$ MeV electron transports through a standing wave radio-frequency (RF) cavity with time-dependent electromagnetic fields. The electromagnetic fields are given as:

[TABLE]

with $r^{2}=x^{2}+y^{2}$ and

[TABLE]

where $\omega$ is the RF angular frequency of the cavity, $\theta$ is the initial driven phase of the cavity, and $e_{0}(z)$ is the on-axis longitudinal electric field. In this case, the RF frequency is $1.3$ GHz, and the initial phase is $224$ degree, and the on-axis electric field $e_{0}(z)$ is shown in Fig. 5.

The electron kinetic energy evolution through the RF cavity is shown in Fig. 6.

The electron is accelerated from the initial $10$ MeV to the final about $37$ MeV at the exit of the cavity. Figure 7 shows the electron trajectory evolution through the cavity.

The electron is focused from the transverse electromagnetic forces through the cavity.

Figure 8 shows the relative numerical errors at the exit of the cavity as a function of step size from the $4^{th}$ order extension of the Boris integrator (magenta), the Vay integrator (green), and the Higuera-Cary integrator (blue) together with a plot of the $4^{th}$ power polynomial. It is seen that in this example, all three $4^{th}$ order numerical integrators have nearly the same relative errors and converge as $4^{th}$ power of the step size.

IV Conclusion and Discussion

In the paper, three second order, time reversible numerical integrators were extended to $4^{th}$ and arbitrary even order accuracy integrators following the split-operator method. These high order numerical integrators have the potential to significantly save computational cost with a given numerical error tolerance. The two high order relativistic integrators can also be used to track the charged particle with large relativistic factor in electromagnetic fields. These integrators implemented in some modern beam dynamics simulation code impact-t will be a useful tool for high brightness electron beam dynamics study.

The extension to higher order accuracy presented in this paper is not limited to the above three integrators. The same extension can be applied to the other time reversible second order relativistic integrators qin ; petri , which were noticed by the author after this work had been done.

ACKNOWLEDGEMENTS

Work supported by the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research used computer resources at the National Energy Research Scientific Computing Center.