\documentclass[11pt]{article}
% \usepackage{multicol}
\usepackage{geometry}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{bbold}
\usepackage{hyperref}
\usepackage[capitalize]{cleveref}%to use \cref
%\usepackage{color,soul}
\geometry{margin=3cm}
%Math stuffs
\DeclareMathOperator{\EE}{E}
\DeclareMathOperator{\EK}{K}
\newcommand{\dd}{{\rm{d}}}
% Bold vectors, only works for single characters
\def\*#1{\mathbf{#1}}
\def\.#1{\phantom{}^#1}
% Bolt mathcal
\DeclareMathAlphabet\mathbfcal{OMS}{cmsy}{b}{n}
\title{Symplectic and Self-Consistent Algorithms for Particle Accelerator Simulation.}
\author{Thomas Planche, Paul M.\ Jung}
\begin{document}
\maketitle
\abstract{
%The focus of this paper is given to particle accelerator simulation, even though the formalism used is often borrowed from the plasma physics community.
This paper is a review of algorithms, applicable to particle accelerator simulation, which share the following two characteristics:
%In a first part we review existing algorithms that share two properties:
(1) they preserve to machine precision the symplectic geometry of the particle dynamics, and
(2) they track the evolution of the self-field consistently with the evolution of the charge distribution.
This review includes, but is not limited to, algorithms using a Particle-in-Cell discretization scheme.
At the end of this review we discuss to possibility to derived algorithms from an electrostatic Hamiltonian.}
\section{Introduction}
The conventional approach to simulating charged particle dynamics is to start from equations of motion, such as Newton's law and Poisson's equation, and solve them approximately using some standard ordinary or partial differential equation solver. The truncation errors often lead to non-physical artifacts, such as the non-conservation of phase space volume, or the violation of conservation laws resulting from symmetries of the system (Noether's theorem).
By contract, symplectic integrators produce exactly stationary solutions to an approximate action. Solutions are exact to machine precision. Approximations are made up front, when choosing the approximate Lagrangian or Hamiltonian, and the corresponding approximate discrete action. Once the physical description of the system is chosen, there is no more arbitrariness in the arcane of the algorithm.
In accelerator physics, symplectic integrators are primarily used to study long-term stability of orbits in storage rings~\cite{Forest-2006}. But their properties, the first of which is the lack of arbitrariness after the choice of physics, make them desirable from the study of all conservative processes in particle accelerators.
Self-consist algorithms are, on the other hand, essential to study betatron resonances, and by extension dynamical aperture, in the presence of space-charge forces~\cite{baartman-1998-shelter}.
\section{Notations}
Throughout this paper a bold character always denotes a vector, a vector field, or a matrix. As is the convention in classical field theory, we use a dot to denote a {\it partial} derivative with respect to time $t$. All formulas are given in SI units, and we use $c$, $\epsilon_0$, and $\mu_0$ to denote respectively: the speed of light, the vacuum permittivity, and permeability.
\section{From a Single-Particle Hamiltonian}
A first class of symplectic and self-consistent algorithms may be derived from a single-particle Hamiltonian. For the sake of demonstration we consider a set of (macro-)particles whose dynamics is governed by the following Hamiltonian:
%The first approach consists in modelling the beam as an ensemble of macro-particles whose dynamics is govern by a single-particle Hamiltonian. For the sake of simplicity we choose the following Hamiltonian:
%Let us start from a single-particle Hamiltonian and model the beam as an ensemble of macro-particles whose dynamics is govern by it.
%For the sake of simplicity let us consider the following Hamiltonian:
\begin{align}
H(\*x, \*P;t)=\frac{\*P^2}{2 m}+q\phi(\*x)+q\psi(\*x),
\end{align}
where $m$ is the mass of the particle, $q$ is its charge, $\*x$ and $\*P$ are its coordinate vector and associated canonically conjugated momentum. The space-charge force derives from the self-potential $\phi$. The external focusing forces derive from the scalar potential $\psi$.
Since this Hamiltonian has no explicit dependence on the independent variable $t$, and is the sum of terms depending on either position or momentum alone, the particle motion can be numerically integrated using a concatenation of jolt maps~\cite{Forest-4rthorder-1990}:
\begin{align}
e^{-\Delta t:H:}\approx\left(I-\frac{\Delta t}{2}:\frac{\*P^2}{2 m}:\right)\left(I- \Delta t\, q:\phi+\psi:\right)\left(I-\frac{\Delta t}{2}:\frac{\*P^2}{2 m}:\right)\,,
\end{align}
This approximate map, accurate to second order in $\Delta t$, is
symplectic if $\phi$ and $\psi$ are functions of class $C^2$. This
is shown by proving that the Jacobian matrix of the map is symplectic~\cite{Qiang-2017}. Higher order integrators may be derived from this second order one using Yoshida's method~\cite{Yoshida-1990}.
With this approach the numerical method for solving the equation of motion for the self-potential -- namely Poisson's equation -- is not obtained from a least action principle. This leads to a certain level of arbitrariness in the way the self-potential is to be computed.
An algorithm of this type, which uses a spectral Poisson solver, has demonstrated excellent long-term stability for 2 and 3-dimensional simulations~\cite{Qiang-2017}.
%Another approach consists in deriving the algorithm entirely from a least-action principle. This can be done in a number of ways which will shall attempt to review.
\section{From a Discretized Lagrangian}
Let us now consider methods based on variational integrators derived from a Lagrangian. We will see that with these methods all the dymanics -- the evolution of the particle distribution as well as the evolution of the self-field -- are obtained from Hamilton's principle of stationary action.
\subsection{Low's Lagrangian}
To illustrate this approach we choose to start from the Lagrangian for non-relativistic collisionless plasma proposed by Low~\cite{Low-1958}. In the electrostatic limit, where the self-field derives solely from a scalar potential, it writes:
\begin{align}\label{eq:Low}
L(\*x,\dot{\*x},\phi;t) = \int f(\*x_0, \dot{\*x}_0) L_P(\*x(\*x_0, \dot{\*x}_0,t),\dot{\*x}(\*x_0, \dot{\*x}_0,t);t)\,\dd \*x_0 \dd \dot{\*x}_0
+\frac{\epsilon_0}{2}\int |\nabla\phi(\bar{\*x},t)|^2\dd \bar{\*x}.
\end{align}
where $L_P$ has the form of a a single-particle Lagrangian:
\begin{align}
L_P(\*x,\dot{\*x};t)=\frac{m}{2}|\dot{\*x}|^2-q\phi(\*x,t)\,,
\end{align}
$\*x$ and $\dot{\*x}$ are vector fields that map the initial coordinates $(\*x_0, \dot{\*x}_0)$ to the corresponding coordinates at time $t$. $f$ is the initial plasma density function. $\bar{\*x}$ is a dummy integration variable.
Note that the single particle Lagrangian can be made more general by adding external fields term, and by replacing the non-relativistic kinetic energy term by $-m\sqrt{1-|\dot{\*x}|^2/c^2}$. For the sake of simplicity, and without much loss of generality, we choose to put aside these refinements.
\subsection{Discretized Lagrangian}
Let's now discretize our system, i.e.\ approximate the continuous system by one with a finite number of degrees of freedom. The choice of the discretization scheme, for both the phase-space distribution $f$ and real-space potential $\phi$, is arbitrary. For the sake of illustration we choose the following particular Particle-in-Cell (PIC) scheme:
\begin{align}\label{eq:PICdiscr}
f(\*x_0, \dot{\*x}_0)&=\sum_i w^i \,\delta(\*x_0-\*x^i_0)\delta(\dot{\*x}_0-\*v^i_0)\,,\\
\phi(\*x,t)&=\sum_j\phi^j(t) K(\*x-\*x^j)\,,
\end{align}
where $w^i$ is the weight of the $i^{\rm th}$ macro-particle, $\*x^i_0$ and $\*v^i_0$ are its initial coordinates.
$\*x^j$ is the position of $j^{\rm th}$ node of the PIC grid, and $\phi^j(t)$ is the potential assigned to this node. $\delta$ is the Dirac function. $K$ is an interpolation kernel function of class $C^2$ which satisfies the requirement of norming~\cite{grigoryev2012numerical}:
\begin{align}
\sum_j K(\*x - \*x^j) = 1\,,
\end{align}
for all $\*x$.
The choice of the kernel function is arbitrary. It is usually chosen among positively defined even functions. A noticeable example of kernel function is the Gaussian wavelet used in {\tt COSY INFINITY} (see section on `General Particle Optical Elements' in~\cite{COSYman}). Suitable kernel functions may also be constructed out of piecewise polynomials~\cite{Qin-PIC-2016}.
Note that this discretization scheme is similar, although not identical, to the one used in Ref.~\cite{shadwick2014variational}.
%The purpose of the numerical simulation is to compute the time evolution of each individual degree of freedom, namely $\*x_i(t)$, $\dot{\*x}_i(t)$, and $\phi_i$.
Combining~\cref{eq:Low,eq:PICdiscr} leads to the discretized Lagrangian:
\begin{equation}
\begin{aligned}
L_D(\*x,\dot{\*x},\phi;t)&=\frac{m}{2}\sum_i w^i |\dot{\* x}^i(t)|^2 -q\sum_i \sum_j \phi^j w^i K(\*x^i(t)-\*x^j)\\
&+\frac{\epsilon_0}{2}\int \bigg(\sum_j\phi^j\nabla K(\bar{\*x}-\*x^j)\bigg)^2\,\dd \bar{\*x}\,,
\end{aligned}
\end{equation}
where $\*x^i(t)=\*x(\*x^i_0,\*v^i_0,t)$.
%Several other discretization schemes have been proposed and tested: for a similar approach based on a spectral descretization of real-space see~\cite{webb-spectral-2016}; for spectral discretization of phase-space to solve steady-state problems see~\cite{gold-2018}; for a discretization of the phase-space on a regular grid, see the extensive literature on Eulerian Vlasov codes. %The choice of the most appropriate discretization scheme depends on the specificity of the problem one is trying to solve.
\subsection{Discretized Action and Equations of Motion}
%We are now seeking to obtain equations of motion from the discretized Lagrangian using Hamilton's principle of least action.
The action $S=\int L_D \dd t$ can be approximated to first order using a Riemann sum:
\begin{align}\label{eq:Riemann}
S\approx S_{D}= \Delta t\sum_n \, L_D(\*x_{n},\frac{\*x_{n+1}-\*x_{n}}{\Delta t},\phi_{n};t)\,,
\end{align}
where the subscript $n$ denotes an evaluation at time $t=n\,\Delta t$. Minimization of the action follows from the discrete Euler-Lagrange equation (see section 1.1.2 of~\cite{marsden-2001}) which leads to the following equations of motion:
\begin{align}\label{eq:9}
m\frac{\*x^i_{n+1}-2\,\*x^i_{n}+\*x^i_{n-1}}{\Delta t}&=-q\sum_j \phi_{n}^j \nabla K(\*x^i_{n}-\*x^j)\,,\\\label{eq:10}
\sum_k \phi_{n}^k \mathcal{M}^{jk} &=-\frac{q}{\epsilon_0}\rho_{n}^j\,,
\end{align}
where:
\begin{align}
\mathcal{M}^{jk}=\int K(\bar{\*x}-\*x^j)\nabla^2 K(\bar{\*x}-\*x^k) \,\dd \bar{\*x}\,,
\end{align}
and
\begin{align}
\rho_{n}^j=\sum_i w^i K(\*x_n^i-\*x^j)\,.
\end{align}
Equation~\ref{eq:9} is a discrete form of Newton's equation with the Lorentz force.
To solve it numerically one may split this second-order equation into two first order equations~\cite{webb-spectral-2016}.
Equation~\ref{eq:10} is obtained after integrating by parts the outcome from the Euler-Lagrange equation, and dropping the boundary term. It is a discrete form of Poisson's equation. It defines a linear relation between all $\phi_{n}^j$ and $\rho_{n}^j$ and can be solved by inverting the square matrix $\mathbfcal{M}$.
%exactly using conventional linear algebra techniques.
Equations~\ref{eq:9} and~\ref{eq:10} constitute a complete numerical integration scheme. Numerical integration leads to an exact (to machine precision) solution of an approximate action, which makes it a variational integrator. Variational integrators are symplectic integrators~\cite{marsden-2001}. As a matter of fact this particular one is a symplectic Euler.
Higher-order variational integrators can be obtained using higher-order approximations of $S$. A second order variational inegrator using a spectral discretization of $\phi$ has been tested in one, two, and three-dimensional, and has demonstrated excellent long-term stability~\cite{webb-spectral-2016}.
\section{From a Discretized Hamiltonian}
Symplectic integrators are more commonly obtained from a Hamiltonian~\cite{Forest-2006}.
% Unfortunately the electrostatic Low's Lagrangian in~\cref{eq:Low} cannot be Legendre transformed into a Hamiltonian, since it has no explicit dependence on $\dot{\phi}$. In this section we will discuss two ways to overcome this degeneracy to obtain a Hamiltonian from Low's Lagrangian.
Unfortunately the electrostatic Low Lagrangian in~\cref{eq:Low} is degenerate: it contains no explicit dependence on $\dot{\phi}$. This makes the application of a Legendre transformation to this Lagrangian, if not impossible, at least beyond the abilities of the authors.
In this section we will discuss two ways to overcome this issue and obtain a Hamiltonian from different versions of Low's Lagrangian.
\subsection{Electromagnetic Hamiltonian}
Let us choose to use the Weyl gauge (also referred to as temporal gauge) $\phi=0$.
%Under this condition the electric field derives entirely from a vector potential:
% \begin{align}
% \mathbfcal{E}=-\dot{\*A}\,.
% \end{align}
Low's Lagrangian becomes:
\begin{align}\label{eq:WeylLow}
L(\*x,\dot{\*x},\*A,\dot{\*A};t)&=\int f(\*x_0,\dot{\*x}_0)\left(\frac{m}{2} |\dot{\*x}(\*x_0,\dot{\*x}_0,t)|^2+q\,\dot{\*x}\cdot \*A(\*x,t)\right) \,\dd \*x_0 \dd \dot{\*x}_0\\
&+\frac{\epsilon_0}{2}\int |\dot{\*A}(\bar{\*x},t)|^2 - |c\nabla\times \*A(\bar{\*x},t)|^2 \,\dd \bar{\*x}\,,
\end{align}
\subsection{Discretization and Legendre Transformation}
Let's use a PIC discretization scheme identical to~\cref{eq:PICdiscr}, only replacing $\phi$ by $\*A$. For compactness we write the discretized Lagrangian in matrix form:
\begin{align}
L_D&=\frac{m}{2} \dot{\* x}^\intercal\, \*w\, \dot{\* x} +q \dot{\* x}^\intercal\, \*w \*K\,\*A+\frac{\epsilon_0}{2}\dot{\*A}^\intercal\,\mathbfcal{K}\,\dot{\*A}-\frac{1}{2 \mu_0}\*A^\intercal\,\mathbfcal{K}_{_\times}\*A
\end{align}
% \begin{align}
% L_D&=\frac{m}{2}\sum_i w^i |\dot{\* x}^i(t)|^2 +q\sum_i \sum_j w^i W(\*x^i(t)-\*x^j) \, \* A^j \cdot \dot{\* x}^i(t) \\
% &+\frac{\epsilon_0}{2}\int \bigg|\sum_j \dot{\*A}^j(t) W(\bar{\*x}-\*x^j) \bigg|^2\,\dd \bar{\*x}
% % -\frac{1}{2 \mu_0} \int \bigg|\sum_j \mathbb{W}(\bar{\*x}-\*x^j) \*A^j(t) \bigg|^2\,\dd \bar{\*x}
% -\frac{1}{2 \mu_0} \int \bigg|\sum_j \nabla\times\left(\*A^j(t) W(\bar{\*x}-\*x^j)\right)\bigg|^2\,\dd \bar{\*x}\,.
% \end{align}
where $\*x$ and $\*A$ are now vectors with components $\*x^i$ and $\*A^j$ respectively;
$\*x$ contains $N_p$ elements (the number of macro-particles), and $\*A$ contains $N_g$ elements (the number of grid nodes). $\*w$ is a diagonal matrix with components $w^i$. The components of the other matrices are:
\begin{align}
K^{ij}&=K(\*x^i(t)-\*x^j)\\
\mathcal{K}^{jk}&=\int K(\bar{\*x}-\*x^j)K(\bar{\*x}-\*x^k) \,\dd \bar{\*x}\\
\mathcal{K}_{_\times}^{jk}&=\int [\nabla K(\bar{\*x}-\*x^j)]_{_\times}^\intercal\, [\nabla K(\bar{\*x}-\*x^k)]_{_\times} \dd \bar{\*x}\,,
\end{align}
where $i$ and $l$ go from 1 to $N_p$, while $j$ and $k$ go from 1 to $N_g$.
Superscript $^\intercal\,$ refers to the transpose operation. $[\quad]_{_\times}$ denotes a skew matrix used to express the cross product as a matrix multiplication ($[\*a]_{_\times} \*b=\*a\times \*b$).
The Legendre transformation writes:
\begin{align}\label{eq:LegendreEM}
H_D= \dot{\* x}^\intercal\,\*P+ \dot{\*A}^\intercal\, \*Y-L_D\,,
\end{align}
The components of the canonical momentum vectors are $\*P^i=\frac{\partial L_D}{\partial \dot{\* x}^i}$ and $\*Y^j=\frac{\partial L_D}{\partial \dot{\*A}^j}$, which are explicitly as:
\begin{align}
\*P&=m \*w \dot{\* x}+q\*w\*K \*A\\
\*Y&=\epsilon_0 \mathbfcal{K} \dot{\*A}\\
\end{align}
The discretized Hamiltonian becomes:
\begin{align}\label{eq:HDEM}
H_D=\left( \*P-q \*w\*K \*A \right)^2+\frac{1}{2\epsilon_0} \*Y^\intercal\,\mathbfcal{K}^{\intercal^{-1}}\,\*Y+\frac{1}{2 \mu_0}\*A^\intercal\,\mathbfcal{K}_{_\times}\*A\,,
\end{align}
and the associated canonical Poisson bracket writes:
\begin{align}
\{F,G\}=\sum_i \frac{\partial F}{\partial\*x^i}\frac{\partial G}{\partial\*P^i}-\frac{\partial F}{\partial\*P^i}\frac{\partial G}{\partial\*x^i} +\sum_j \frac{\partial F}{\partial\*A^j}\frac{\partial G}{\partial\*Y^j}-\frac{\partial F}{\partial\*Y^j}\frac{\partial G}{\partial\*A^j}\,.
\end{align}
Since $H_D$ is the sum of exactly solvable terms (for an explicit solution of the $\left( \*P-q \*w\*K \*A \right)^2$ term, see~\cite{Wu-2003}), one can build a second order symplectic integrator by concatenating maps~\cite{Forest-4rthorder-1990}.
A similar symplectic integrator derived from the Morrison-Marsden-Weinstein electromagnetic Hamiltonian~\cite{morrison-1980,weinstein1981comments,marsden1982hamiltonian} has been tested~\cite{Qin-PIC-2016}. Note that a corresponding variational integrator had previously been tested~\cite{Qin-variational-2008}.
\subsection{Electrostatic Hamiltonian}
In most accelerator physics problems particles do not move at relativistic speeds with respect to each other. In such a case a scalar potential is sufficient to describe the self-field. Keeping track of the three components of a vector potential is wasteful.
We have already discussed the fact that an electrostatic Hamiltonian cannot be obtained from~\cref{eq:Low}.
In this section we show that it is however possible to obtain an electrostatic Hamiltonian after changing the independent variable.
The action associated with~\cref{eq:Low} writes:
\begin{align}\label{eq:Actiont}
S = \iint f\left(\frac{m}{2}|\dot{\*x}|^2-q\phi(\*x,t) \right)\,\dd \*x_0 \dd \dot{\*x}_0\dd t
+\frac{\epsilon_0}{2}\iint |\nabla\phi|^2\dd \bar{\*x}\dd t.
\end{align}
We proceed in the first integral to a change of variable using the substitution function:
\begin{align}
g(x_0,y_0,t_0,x_0',y_0',t_0',z)=(x_0,y_0,z_0,\dot{x}_0,\dot{y}_0,\dot{z}_0,t)\,,
\end{align}
where primes $'$ denote partial derivative w.r.t.\ $z$. The determinant of the Jacobian matrix $\det\left( D_g \right)=t'/(t_0'^5)$,
% \begin{align}
% \det\left( D_g \right)=\frac{t'}{t_0'^5}\,
% \end{align}
where $t_0'=\frac{\partial t}{\partial z}\big|_{z=0}$.
Similarly we proceed in the second integral to the change of variable given by:
\begin{align}
h(\bar{x},\bar{y},\bar{t},z)=(\bar{x},\bar{y},\bar{z},t)\,,
\end{align}
The determinant of the Jacobian of $h$ is 1.
The action now writes:
\begin{align}\label{eq:Actionz}
S &= \int \left( \int \hat{f}\left(m\frac{x'^2+y'^2+1}{2 t'}-q\,t'\hat{\phi}\right) \,\dd x_0 \dd y_0 \dd t_0 \dd x_0' \dd y_0' \dd t_0'
+\frac{\epsilon_0}{2}\int |\nabla\hat{\phi}|^2\,\dd \bar{x} \dd \bar{y} \dd \bar{t} \right)\dd z\\
&= \int L_z(x,y,t,\hat{\phi},x',y',t',\hat{\phi}';z) \, \dd z\,,
\end{align}
%where $L_z(x,y,t,\phi,x',y',t',\phi';z)$ is the Low's Lagrangian with the space coordinate $z$ as independent variable.
where:
\begin{align}
&\hat{f}(x_0,y_0,t_0,x_0',y_0',t_0',z)=f(x_0,y_0,z_0,\dot{x}_0,\dot{y}_0,\dot{z}_0,t)/t_0'^5, \quad {\rm and}\\
&\hat{\phi}(x_0,y_0,t_0,x_0',y_0',t_0',z)=\phi(x_0,y_0,z_0,\dot{x}_0,\dot{y}_0,\dot{z}_0,t)\,.
\end{align}
To simplify notations we drop the hats.
One notices that $L_z$ depends explicitly on $\phi'$, $x'$, $y'$, and $t'$, enabling us to define the following canonical momentum densities:
\begin{align}
\Pi&=\epsilon_0 \phi'\\
P_x&=mf\frac{x'}{t'}\\
P_y&=mf\frac{y'}{t'}\\
-E&=-mf\frac{x'^2+y'^2+1}{2t'^2}-qf\phi\\
\end{align}
which in turn enables us to perform a Legendre transformation and obtain the following electrostatic Hamiltonian:
\begin{align}\label{eq:Hamiltonianz}
H_z&=\int \sqrt{2mf(E-qf\phi)-P_x^2-P_y^2} \,\, \dd x_0 \dd y_0 \dd t_0 \dd x_0' \dd y_0' \dd t_0'
+\frac{1}{2}\int \frac{\Pi^2}{\epsilon_0} -\epsilon_0 \nabla_{\perp}\phi\,\, \dd \bar{x} \dd \bar{y} \dd \bar{t}\,.
\end{align}
%As far as the authors know, no algorithm derived an electrostatic Hamiltonian has been tested yet.
An attempt to implement an algorithm based on a relativistic version of this Hamiltonian is presented Ref.~\cite{paulICAP18}.
%$m^3m^3s^{-3}s=m^6s^{-2}$\\
%$m^2s^2$\\
%$m^4s^{-4}$
%The single particle Lagrangian can be made more general by including external scalar and vector potential terms, and by replacing the kinetic energy term $\frac{m}{2}\dot{\*x}^2$ by $-mc^2\sqrt{1-\frac{\dot{\*x}^2}{c^2}}$.
\section{Conclusion}
A variety of symplectic and self-consistent multi-particle algorithms have been developed by both the accelerator physics and the plasma physics community. They are superior to most algorithm derived from equations of motion as they guaranty that the symplectic nature of the particle dynamics is conserved to machine precision. Some of them guaranty the conservation of the symplectic nature of the self-field dynamics as well.
The algorithms discussed in this paper were all derived from a collision-less picture, and as such are unable to model non-Hamiltonian processes such as intra-beam scattering.
\section{Acknowledgement}
The authors would like to thank Rick Baartman and Kyle Gao for their corrections and comments on the original manuscript.
\bibliographystyle{elsarticle-num}
\bibliography{./bib}
\end{document}