%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    INSTITUTE OF PHYSICS PUBLISHING                                   %
%                                                                      %
%   `Preparing an article for publication in an Institute of Physics   %
%    Publishing journal using LaTeX'                                   %
%                                                                      %
%    LaTeX source code `ioplau2e.tex' used to generate `author         %
%    guidelines', the documentation explaining and demonstrating use   %
%    of the Institute of Physics Publishing LaTeX preprint files       %
%    `iopart.cls, iopart12.clo and iopart10.clo'.                      %
%                                                                      %
%    `ioplau2e.tex' itself uses LaTeX with `iopart.cls'                %
%                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% First we have a character check
%
% ! exclamation mark    " double quote  
% # hash                ` opening quote (grave)
% & ampersand           ' closing quote (acute)
% $ dollar              % percent       
% ( open parenthesis    ) close paren.  
% - hyphen              = equals sign
% | vertical bar        ~ tilde         
% @ at sign             _ underscore
% { open curly brace    } close curly   
% [ open square         ] close square bracket
% + plus sign           ; semi-colon    
% * asterisk            : colon
% < open angle bracket  > close angle   
% , comma               . full stop
% ? question mark       / forward slash 
% \ backslash           ^ circumflex
%
% ABCDEFGHIJKLMNOPQRSTUVWXYZ 
% abcdefghijklmnopqrstuvwxyz 
% 1234567890
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\documentclass[]{iopart}
% Uncomment next line if AMS fonts required
%\usepackage{iopams}  
\usepackage[]{epsfig}  
\begin{document}

\title{Boris push with spatial stepping}


\author{G. Penn\footnote[3]{gpenn@socrates.berkeley.edu}, P. H. Stoltz, 
J. Wurtele\dag, and J. R. Cary\ddag}

\address{\dag\ UC Berkeley Dept. of Physics, Berkeley, CA 94720}

\address{\ddag\ Tech-X Corporation, Boulder, CO 80301}


\begin{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.
\end{abstract}

%Uncomment for PACS numbers title message
%\pacs{00.00, 20.00, 42.10}

% Uncomment for Submitted to journal title message
%\submitto{\JPA}

% Comment out if separate title page not required
%\maketitle

\section{The spatial Boris push algorithm}

Plasma physics simulations often use the Boris push algorithm\cite{boris}
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\cite{stoltz} 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\cite{icool}.

For the spatial Boris scheme, the equations for particle motion
are written in the form:
\begin{eqnarray}
{dp_x \over dz}&=&q\left({E_x \over v_z} + {v_y B_z \over v_z} - B_y\right),
\ \ 
{dp_y \over dz}=q\left({E_y \over v_z} - {v_x B_z \over v_z} + B_x\right),
\ \ \cr
{dU/c \over dz}&=&q\left({v_x E_x \over v_z c} + {v_y E_y \over v_z c} 
+ {E_z \over c}\right).
\end{eqnarray}
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
\begin{equation}
{dw\over dz} = M w + b, 
\label{eq:evolve}
\end{equation}
where $w$ is the vector consisting of $p_x$, $p_y$, and $U/c$;
the vector and matrix terms are
\begin{equation}
b = q \left(\begin{array}{c} -B_y\\ B_x\\ E_z/c\end{array}\right),\ \ 
M = {q\over p_z} 
\left(
\begin{array}{c} 0\\ -B_z\\ E_x/c\end{array}
\begin{array}{c} B_z\\ 0\\ E_y/c\end{array}
\begin{array}{c} E_x/c\\ E_y/c\\ 0\end{array}
\right).
\end{equation}

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.~(\ref{eq:evolve}) 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.~(\ref{eq:evolve}) 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 step-centered scheme will be second-order accurate.
To evolve the momenta from $w^n$ to $w^{n+1}$,
this implies solving the equations:
\begin{equation}
w^{-} = w^n + {\Delta \over 2} b, \\ \\
w^{+} - w^{-} = \Delta M {w^{+} + w^{-} \over 2}, \\ \\
w^{n+1} = w^{+} + {\Delta \over 2} b,
\end{equation}
where $\Delta$ is the step size.
Solving for the intermediate value $w^{+}$ yields:
\begin{equation}
w^{+} = \left(I - M \Delta /2\right)^{-1} \left(I + M \Delta /2\right) 
w^{-} \equiv (I+R) w^{-},
\end{equation}
where
\begin{equation}
R = \alpha
\left(
\begin{array}{c}
\delta (E_x^2-c^2B_z^2) \\
-c B_z+\delta E_xE_y   \\
E_x-\delta c B_zE_y
\end{array}
\begin{array}{c}
c B_z+ \delta E_xE_y     \\
\delta (E_y^2-c^2B_z^2) \\
E_y+\delta c B_zE_x     
\end{array}
\begin{array}{c}
E_x+\delta cB_zE_y \\
E_y-\delta cB_zE_x  \\
\delta (E_x^2+E_y^2)
\end{array}
\right),
\end{equation}
$\delta \equiv q \Delta /2 p_z c$,
and 
$\alpha = 2 \delta
[1+\delta^2(c^2B_z^2-E_x^2-E_y^2)]^{-1}$.


\section{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\cite{dullweber}.
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 $\Delta$,
the errors, $e$ , scale as
\begin{equation}
e_{RK} \sim R_0 \left({\Delta_{RK} \over L}\right)^5 {z \over L}, \ \ 
e_{B} \sim B_0 \left({\Delta_{B} \over L}\right)^2.
\end{equation}
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_\perp$ is conserved and in fact is exactly
preserved by the Boris push scheme.  However, $p_\perp$ decays exponentially
for the Runge-Kutta scheme.
With a 10 T field, $\pi$ 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 above the cutoff for this channel.
After 200 meters, the crossover point is 
$e \approx 3 \times 10^{-5}$,
well below fluctuations for an ensemble of $\sim 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.

\begin{figure}
\begin{center}
\epsfysize=4.0 true cm
\epsfbox{penn_fig1.eps}
\end{center}
\caption{\label{fig:fofoboris}Relative errors for 
spatial Boris push in FOFO lattice,
for step sizes of 5 and 10 cm.}
\end{figure}

\begin{figure}
\begin{center}
\epsfysize=4.0 true cm
\epsfbox{penn_fig2.eps}
\end{center}
\caption{\label{fig:fofork}Relative errors for Runge-Kutta algorithm
for different step sizes, in a uniform solenoid and FOFO lattice after 200 m.}
\end{figure}

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.

\section*{References}

\begin{thebibliography}{1}

\smallskip
\bibitem{boris}
J. Boris, in
{\it Proc.~of the 4th Conf.~on Numerical Simulation of Plasmas} (NRL, 1970),
pp.~3-67.

\smallskip

\smallskip
\bibitem{stoltz}
P.H.~Stoltz, G.~Penn, J.~Wurtele, and J.R.~Cary, 
Phys.~Rev.~S.T.-A.B.~5, 094001 (2002).
\smallskip

\smallskip
\bibitem{icool}
R.C.~Fernow, in {\it Proc. of the 1999 Particle Accelerator
Conference, New York} (IEEE, Piscataway, NJ, 1999),
p.~3020.
\smallskip

\smallskip
\bibitem{dullweber}
A.~Dullweber, B.~Leimkuhler, and R.~McLachlan, 
J.~Chem.~Phys.~15 (1997) 5840.
\smallskip

\end{thebibliography}


\end{document}


