Boris push with spatial stepping

G. Penn§, P. H. Stoltz, J. Wurtele†, and J. R. Cary‡

† UC Berkeley Dept. of Physics, Berkeley, CA 94720

‡ Tech-X Corporation, Boulder, CO 80301

Abstract. The Boris push is commonly used in plasma physics simulations because of its speed and stability. It is second order accurate, requires only one field evaluation per time step, and has good conservation properties. However, for accelerator simulations it is convenient to propagate particles in z down a changing beamline. A “spatial Boris push” algorithm has been developed which is similar to the Boris push but uses a spatial coordinate as the independent variable, instead of time. This scheme is compared to the 4th order Runge-Kutta algorithm, for two simplified muon beam lattices: a uniform solenoid field, and a “FOFO” lattice where the solenoid field varies sinusoidally along the axis. Examination of the canonical angular momentum, which should be conserved in axisymmetric systems, shows that the spatial Boris push improves accuracy over long distances.

1. The spatial Boris push algorithm

Plasma physics simulations often use the Boris push algorithm[1] because of its efficiency and stability. It is a leapfrog scheme where the particle is moved, then half of the energy changed is applied. After this, the momentum is rotated by the magnetic field, and the rest of the energy change applied. The rotation is automatically energy conserving, and the algorithm is symmetric to time reversal, which improves the performance. This algorithm is second order accurate and requires only one field evaluation per time step. However, the Boris push is rarely used for accelerator simulations, where it is desirable to propagate particles in space instead of time. The spatial Boris push[2] adapts this algorithm to use a spatial independent variable by interchanging the roles of energy U and forward momentum p_{z}, as well as t and z. An implementation of the new scheme is available in the ICOOL simulation code[3].

For the spatial Boris scheme, the equations for particle motion are written in the form:

\frac{dp_{x}}{dz} = q (\frac{E_{x}}{v_{z}} + \frac{v_{y}B_{z}}{v_{z}} − B_{y}) , \frac{dp_{y}}{dz} = q (\frac{E_{y}}{v_{z}} − \frac{v_{x}B_{z}}{v_{z}} + B_{x}) ,
\frac{dU/c}{dz} = q (\frac{v_{x}E_{x}}{v_{z}c} + \frac{v_{y}E_{y}}{v_{z}c} + \frac{E_{z}}{c} ) .  (1)

In the temporal Boris scheme, the electric fields modify U, while magnetic fields conserve energy. For the spatial Boris push, the terms that directly change p_{z} are the field components E_{z}, B_{x}, and B_{y}. Noting that v_{z} = p_{z}c^{2}/U, we can write

\frac{dw}{dz} = Mw + b, (2)

§ gpenn@socrates.berkeley.edu

where w is the vector consisting of p_{x}, p_{y}, and U/c; the vector and matrix terms are

b = q ( \begin{matrix}−B_{y} \\ B_{x} \\ E_{z}/c \end{matrix} ) , M = \frac{q}{p_{z}} ( \begin{matrix}0 \\ −B_{z} \\ E_{x}/c \end{matrix} \begin{matrix}B_{z} \\ 0 \\ E_{y}/c \end{matrix} \begin{matrix}E_{x}/c \\ E_{y}/c \\ 0 \end{matrix} ) . (3)

The integration cycle consists of a leap-frog scheme in which the particle positions are shifted a half-step, where the fields are evaluated. The vector term of Eq. (2) is propagated a half-step and the forward momentum is recalculated. Then, the matrix term is calculated for a full step, and the vector term for the second half-step. Finally, the particle positions are moved another half-step using the new particle momenta.

In this scheme, the vector and matrix terms of Eq. (2) are integrated separately. To preserve the symmetry for propagating the particles backwards, it is necessary to solve the matrix portion of the evolution implicitly; because the matrix M by itself conserves p_{z}, within this step M can be treated as a constant matrix. Thus, a stepcentered scheme will be second-order accurate. To evolve the momenta from w^{n} to w^{n+1}, this implies solving the equations:

w^{−} = w^{n} + \frac{Δ}{2}b, w^{+} − w^{−} = ΔM \frac{w^{+} + w^{−}}{2} , w^{n+1} = w^{+} + \frac{Δ}{2}b, (4)

where Δ is the step size. Solving for the intermediate value w^{+} yields:

w^{+} = (I − MΔ/2)^{−1} (I + MΔ/2) w^{−} ≡ (I + R)w^{−}, (5)

where

R = α ( \begin{matrix}δ(E_{x}^{2} − c^{2}B_{z}^{2}) \\ −cB_{z} + δE_{x}E_{y} \\ E_{x} − δcB_{z}E_{y} \end{matrix} \begin{matrix}cB_{z} + δE_{x}E_{y} \\ δ(E_{y}^{2} − c^{2}B_{z}^{2}) \\ E_{y} + δcB_{z}E_{x} \end{matrix} \begin{matrix}E_{x} + δcB_{z}E_{y} \\ E_{y} − δcB_{z}E_{x} \\ δ(E_{x}^{2} + E_{y}^{2}) \end{matrix} ) , (6)

δ ≡ qΔ/2p_{z}c, and α = 2δ[1 + δ^{2}(c^{2}B_{z}^{2} − E_{x}^{2} − E_{y}^{2})]^{−1}.

2. Results

The spatial Boris push is compared with the Runge-Kutta scheme, with a focus on canonical angular momentum, which should be conserved in an axisymmetric system. Runge-Kutta is approximately fifth order accurate but is known to produce articial damping or growth of conserved quantities[4]. The Boris push is only second order, but is simpler to calculate, requires fewer field calculations, and preserves invariants. For a characteristic length scale L and step size Δ, the errors, e , scale as

e_{RK} ∼ R_{0} (\frac{Δ_{RK}}{L} )^{5} \frac{z}{L}, e_{B} ∼ B_{0} (\frac{Δ_{B}}{L} )^{2} . (7)

Note that for the Runge-Kutta the exponent = 5, not 4, because of the axisymmetry. For the Boris push, the maximum error is reached in a single betatron oscillation. With 4 times more field calculations per step for the Runge-Kutta algorithm, the spatial Boris push is more efficient when it is acceptable to have e > B_{0}^{5/3}R_{0}^{−2/3}(L/z)^{2/3}.

Two examples are considered, both with 200 MeV/c muons: for a uniform solenoid, p_{⊥} is conserved and in fact is exactly preserved by the Boris push scheme. However, p_{⊥} decays exponentially for the Runge-Kutta scheme. With a 10 T field, π phase advance corresponds to 40 cm, and particles were tracked for 200 m.

A FOFO lattice was also considered, with a sinusoidal field on axis having a 1 m half-period and peak field B = 2 T. The muon momentum, p_{z} = 200 MeV/c, is well

Figure 1. Relative errors for spatial Boris push in FOFO lattice, for step sizes of 5 and 10 cm.

Figure 2. Relative errors for Runge-Kutta algorithm for different step sizes, in a uniform solenoid and FOFO lattice after 200 m.

above the cutoff for this channel. After 200 meters, the crossover point is e ≈ 3×10^{−5}, well below fluctuations for an ensemble of ∼ 10^{6} particles. If the minimum step size is fixed by the geometry or other concerns, the Boris step can be up to 4 times faster and may still be more accurate.

In summary, the spatial Boris push speeds up particle tracking; fields are only evaluated once, where the momentum kick is applied, compared to 4 field evaluations for the Runge-Kutta. The local nature of the kick also simplifies the calculation. By being space-symmetric, the spatial Boris push exhibits the same conservation properties as the temporal Boris push, and the errors, though second order in step size, tend to average out. In the Runge-Kutta scheme, by contrast, errors will slowly accumulate. For muon cooling channels, ionization energy loss and scatter are still applied at the end of the step, requiring the assumption that the energy loss per step is small.

References

[1] J. Boris, in Proc. of the 4th Conf. on Numerical Simulation of Plasmas (NRL, 1970), pp. 3-67.

[2] P.H. Stoltz, G. Penn, J. Wurtele, and J.R. Cary, Phys. Rev. S.T.-A.B. 5, 094001 (2002).

[3] R.C. Fernow, in Proc. of the 1999 Particle Accelerator Conference, New York (IEEE, Piscataway, NJ, 1999), p. 3020.

[4] A. Dullweber, B. Leimkuhler, and R. McLachlan, J. Chem. Phys. 15 (1997) 5840.