21pt 15pt

# Structure-preserving techniques in accelerator physics

Dan T. Abella and Alex J. Dragtb CONTACT A. J. Dragt Email:
<EMAIL_ADDRESS>aRadiaSoft LLC, Boulder, Colorado, USA  bUniversity of
Maryland Physics Department (Professor Emeritus)

###### Abstract

To a very good approximation, particularly for hadron machines, charged-
particle trajectories in accelerators obey Hamiltonian mechanics. During
routine storage times of eight hours or more, such particles execute some
${10}^{8}$ revolutions about the machine, ${10}^{10}$ oscillations about the
design orbit, and ${10}^{13}$ passages through various bending and focusing
elements. Prior to building, or modifying, such a machine, we seek to identify
accurately the long-term behavior and stability of particle orbits over such
large numbers of interactions. This demanding computational effort does not
yield easily to traditional methods of symplectic numerical integration,
including both explicit Yoshida-type and implicit Runge-Kutta or Gaussian
methods. As an alternative, one may compute an approximate one-turn map and
then iterate that map. We describe some of the essential considerations and
techniques for constructing such maps to high order and for realistic magnetic
field models. Particular attention is given to preserving the symplectic
condition characteristic of Hamiltonian mechanics.

###### keywords:

symplectic map, symplectic jet, geometric integration, Lie algebra, Poincaré
generating function, Cremona map, Cremona symplectification

††articletype: Book Chapter

*

###### Contents

  1. 1 Introduction
  2. 2 Lagrangians and Hamiltonians
  3. 3 Transfer Maps, the Symplectic Condition, and Symplectic Integrators
  4. 4 Lie Algebraic Concepts and Tools
  5. 5 Symplectic Completion of Symplectic Jets
  6. 6 Symplectic Completion Using Generating Functions
  7. 7 Symplectic Completion Using Cremona Maps
  8. 8 Concluding Discussion
  9. Acknowledgements

##  1 Introduction

Particle accelerators consist of arrays of magnets and radio-frequency (rf)
cavities. Their purpose is to produce intense beams of high-energy charged
particles including electrons, positrons (the antimatter counterpart of
electrons), protons, antiprotons, and various ions. The magnets provide
magnetic fields that bend, focus, and exert various nonlinear effects on beam
particle orbits; and the rf cavities provide electric fields that accelerate
and longitudinally bunch particles.

Successful accelerator design requires the accurate calculation/simulation of
particle orbits over long periods of time. The most computationally
challenging are orbits in _storage rings_ , essentially circular machines in
which particles continually circulate for long periods of time. Electron (or
positron) storage rings are used to produce intense high-energy X-rays.
Electron (or positron) storage rings, as well as proton (or antiproton or ion)
storage rings, are used in pairs to produce colliders. In a collider, the beam
from one ring collides head-on with the counter-circulating beam in a second
ring.111But, even with the highest achievable beam densities, beam-beam
collisions are sufficiently rare that the stored beams are only partially
depleted over the storage time. However, the colliding beams can have
significant dynamical (both linear and nonlinear) effects on each other, a
complication that must be understood/managed but lies beyond the scope of this
chapter. For example, the Large Hadron Collider (LHC) in CERN (near Geneva,
Switzerland), which collides protons on protons, has a circumference of
$27\text{\,}\mathrm{k}\mathrm{m}$, and each ring has some $19\,200$ elements
(magnets, rf cavities, and intervening drift spaces). Protons moving at
essentially the speed of light are stored for about
$8\text{\,}\mathrm{h}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{s}$, and during
this time make approximately $4\text{\times}{10}^{8}$ turns around the ring,
approximately $8\text{\times}{10}^{12}$ element passages, and approximately
$20\text{\times}{10}^{9}$ (betatron) oscillations about the design orbit.
Following this number of oscillations is comparable to following the earth’s
orbit about the sun from the time of the Big Bang. While the number of
betatron oscillations in electron (or positron) storage rings is comparable,
they need not be followed for as long because those oscillations are damped by
the energy loss associated with X-ray emission. Correspondingly,
electron/positron storage-ring orbits are less computationally challenging
than proton/antiproton/ion storage-ring orbits. This chapter is devoted to the
most challenging problem of calculating/simulating particle orbits over long
periods of time in proton (or antiproton or ion) storage rings.

##  2 Lagrangians and Hamiltonians

In Cartesian coordinates, the relativistic Lagrangian $L$ for the motion of a
particle of mass $m$ and charge $q$ in an electromagnetic field is given by
the expression

$L(\boldsymbol{r},\boldsymbol{v},t)=-mc^{2}(1-v^{2}/c^{2})^{1/2}-q\psi(\boldsymbol{r},t)+q\boldsymbol{v}\cdot\boldsymbol{A}(\boldsymbol{r},t).$
(1)

Here $\boldsymbol{r}$ is the particle position at time $t$,
$\boldsymbol{v}=\mathop{}\\!\mathrm{d}{\boldsymbol{r}}/\mathop{}\\!\mathrm{d}{t}$
is the particle velocity, and $c$ is the speed of light. The quantities $\psi$
and $\boldsymbol{A}$ are the scalar and vector potentials defined in such a
way that the electromagnetic fields $\boldsymbol{E}$ and $\boldsymbol{B}$ are
given by the standard relations

$\displaystyle\boldsymbol{B}$
$\displaystyle=\operatorname{\nabla\times}\boldsymbol{A},$ (2a)
$\displaystyle\boldsymbol{E}$
$\displaystyle=-\operatorname{\nabla}\psi-\partial{\boldsymbol{A}}/\partial{t}.$
(2b)

This formulation ignores spin, radiation reaction (X-ray emission, also
referred to as synchrotron radiation), and quantum effects [LM:init]Much of
the background material for this chapter is most easily found on the Web in a
draft book: A. J. Dragt, _Lie Methods for Nonlinear Dynamics with applications
to Accelerator Physics_ , http://www.physics.umd.edu/dsat/dsatliemethods.html.
It, in turn, provides numerous additional references. In subsequent citations
it will be referred to as _LM_. For a discussion of Lagrangians and
Hamiltonians for charged-particle motion in electromagnetic fields, see _LM_ ,
Sections 1.5–1.7.. These effects may be important over long times for lighter
particles such as electrons and positrons, but they are significantly less
important for heavier particles such as protons and antiprotons and ions.

For the Lagrangian (1) the _canonical_ momentum in Cartesian coordinates is
given by the equation

$\boldsymbol{p}^{\text{can}}=\partial{L}/\partial{\boldsymbol{v}}=m\boldsymbol{v}/(1-v^{2}/c^{2})^{1/2}+q\boldsymbol{A}.$
(3)

Here we use the superscript _can_ to emphasize that (3) defines the
_canonical_ momentum. Note that the first term in (3) is just the relativistic
_mechanical_ momentum,

$\boldsymbol{p}^{\text{mech}}=m\boldsymbol{v}/(1-v^{2}/c^{2})^{1/2}=\gamma
m\boldsymbol{v},$ (4)

where $\gamma$ is the standard relativistic factor

$\gamma=1/(1-v^{2}/c^{2})^{1/2}.$ (5)

Consequently, the relation (3) may also be written in the forms

$\boldsymbol{p}^{\text{can}}=\boldsymbol{p}^{\text{mech}}+q\boldsymbol{A}\quad\text{and}\quad\boldsymbol{p}^{\text{mech}}=\boldsymbol{p}^{\text{can}}-q\boldsymbol{A}.$
(6)

Upon implementing the standard procedure that relates Lagrangians and
Hamiltonians, one finds that the Hamiltonian $H$ associated with the
Lagrangian $L$ specified by (1) is given by the expression

$H=[m^{2}c^{4}+c^{2}(\boldsymbol{p}^{\text{can}}-q\boldsymbol{A})\cdot(\boldsymbol{p}^{\text{can}}-q\boldsymbol{A})]^{1/2}+q\psi=[m^{2}c^{4}+c^{2}(\boldsymbol{p}^{\text{can}}-q\boldsymbol{A})^{2}]^{1/2}+q\psi.$
(7)

In the usual Hamiltonian formulation (as in the usual Lagrangian formulation)
the time $t$ plays the distinguished role of an _independent_ variable, and
all the coordinates $q$ and momenta $p$ are _dependent_ variables.222We are
embarrassed by the custom of also using the symbol $q$ to denote the _charge_
of the particle in question. That is, the canonical variables are viewed as
functions $q(t)$, $p(t)$ of the independent variable $t$. In some cases, it is
more convenient to take some _coordinate_ to be the independent variable
rather than the time, in which case the time becomes a dependent variable. So
doing may facilitate the use of _transfer maps_ , as described in the next
section. For example, consider the passage of a collection of particles
through a rectangular-shaped beam-line element such as a magnet or an rf
cavity. In such a situation, particles with different initial conditions will
require different times to pass through the beam-line element. If the
quantities of interest are primarily the locations and momenta of the
particles as they leave the exit face of the beam-line element, then it would
clearly be more convenient to use for an independent variable a coordinate
that measures the progress of a particle through the beam-line element. With
such a choice, the relation between entering coordinates and momenta and
exiting coordinates and momenta could be treated as a transfer map.
Remarkably, this goal can be achieved within a Hamiltonian framework
[LM:S.1.6]_LM_ , Section 1.6..

######  Theorem 2.1.

Suppose $H(q,p,t)$ is a Hamiltonian for a system having $n$ degrees of
freedom. Suppose further that $\dot{q}_{1}=\partial{H}/\partial{p_{1}}\neq 0$
for some interval of time $T$ in some region $R$ of the phase space described
by the $2n$ variables $(q_{1},\dotsc,q_{n})$ and $(p_{1},\dotsc,p_{n})$. Then,
in this region and time interval, $q_{1}$ can be introduced as an independent
variable in place of the time $t$. Moreover, the equations of motion with
$q_{1}$ as an independent variable can be obtained from a Hamiltonian that
will be called $K$. To construct $K$, define a quantity $p_{t}$ by the rule

$p_{t}=-H(q,p,t).$ (8)

Suppose that this relation is solved for $p_{1}$ to give a relation of the
form

$p_{1}=-K(t,q_{2},\dotsc,q_{n};p_{t},p_{2},\dotsc,p_{n};q_{1}).$ (9)

Such an inversion is possible according to the inverse function theorem
because $\partial{H}/\partial{p_{1}}\neq 0$ by assumption. Then, as the
notation is intended to suggest, $K$ is the desired new Hamiltonian. In this
formulation, $t$ is treated as a coordinate like the remaining $q_{2}$ …
$q_{n}$, and $p_{t}$ is its conjugate momentum.

As an example, let us use this construction to find the Hamiltonian $K$
corresponding to the Hamiltonian $H$ given by (7) when the $z$ coordinate is
taken to be the independent variable. Assume that $\dot{z}>0$ for the
trajectories in question. Then one finds the result

$K=-[(p_{t}+q\psi)^{2}/c^{2}-m^{2}c^{2}-(p_{x}-qA_{x})^{2}-(p_{y}-qA_{y})^{2}]^{1/2}-qA_{z}.$
(10)

Here the quantities $p_{x}$ and $p_{y}$ denote _canonical_ momenta. Note that,
according to (8), $p_{t}$ is usually negative. For the example at hand, one
finds that

$p_{t}=-[m^{2}c^{4}+c^{2}(\boldsymbol{p}^{\text{mech}}\cdot\boldsymbol{p}^{\text{mech}})]^{1/2}-q\psi=-\gamma
mc^{2}-q\psi.$ (11)

There is yet another Hamiltonian formulation that is of interest. In the
spirit of relativity, and following the insight of Hermann Minkowski
(1864–1909), it is reasonable to try to treat space and time on a similar
footing [LM:S.1.6e]_LM_ , Section 1.6 and Exercise 1.6.7. Let us review some
of the mathematical machinery of Special Relativity. Suppose the world-line of
a particle through space-time is parameterized in terms of some parameter
$\tau$ by specifying four functions $x^{\mu}(\tau)$ that, taken together, form
a 4-vector with four contravariant components $x^{\mu}$. We adopt the
convention that the first three components of $x^{\mu}$ are the spatial
coordinates of the particle, and the fourth (with a factor of $c$) is its
temporal coordinate. Specifically (for $\mu=1,2,3,4$ and with $x^{4}=ct$) we
write

$x^{\mu}=(x,y,z,ct)=(\boldsymbol{r},ct).$ (12)

In addition, let $(x^{\prime})^{\mu}$ denote the four derivatives defined by
the equations

$(x^{\prime})^{\mu}=\mathop{}\\!\mathrm{d}{x^{\mu}}/\mathop{}\\!\mathrm{d}{\tau}.$
(13)

Under the assumption that the parameterization is unchanged by a Lorentz
transformation, $(x^{\prime})^{\mu}$ is evidently also a 4-vector, which will
be called the 4-velocity. The 3-velocity, $\boldsymbol{v}$, of a particle is
given by the ratio
$\boldsymbol{v}=(\mathop{}\\!\mathrm{d}{\boldsymbol{r}}/\mathop{}\\!\mathrm{d}{\tau})/(\mathop{}\\!\mathrm{d}{t}/\mathop{}\\!\mathrm{d}{\tau})$.
Since the speed of a massive particle must be less than $c$,
$||\boldsymbol{v}||<c$, it follows that (for physical particles) the
4-velocity must satisfy the condition

$x^{\prime}\cdot x^{\prime}=(x^{\prime})^{\mu}(x^{\prime})^{\nu}g_{\mu\nu}>0.$
(14)

Here $g_{\mu\nu}$ denotes the metric tensor, and we have employed the usual
Einstein convention that repeated indices are to be summed over. In Cartesian
coordinates and for flat space-time, only the diagonal entries of $g$ are
nonzero, and we take them to have the values

$g_{11}=g_{22}=g_{33}=-1,\;g_{44}=1.$ (15)

That is, the space-time interval $\mathop{}\\!\mathrm{d}{s}$ is taken to be
given by the relation

$\mathop{}\\!\mathrm{d}{s}^{2}=g_{\mu\nu}\mathop{}\\!\mathrm{d}{x}^{\mu}\mathop{}\\!\mathrm{d}{x}^{\nu}=c^{2}\mathop{}\\!\mathrm{d}{t}^{2}-(\mathop{}\\!\mathrm{d}\boldsymbol{r})^{2}.$
(16)

We remark that the notation $\mathop{}\\!\mathrm{d}{s}^{2}$ appearing in (16),
although universally employed, can be misleading since, depending on
circumstances, $\mathop{}\\!\mathrm{d}{s}^{2}$ can be negative, zero, or
positive, and is therefore not necessarily the square of anything. But note
that $\mathop{}\\!\mathrm{d}{s}^{2}>0$ for time-like displacements. Space-time
endowed with the metric (15) is sometimes called _Minkowski_ space.

The metric tensor can be used to raise and lower indices. For example, there
are the relations

$x_{\mu}=g_{\mu\nu}x^{\nu}.$ (17)

In particular, $x_{\mu}$ has the entries

$x_{\mu}=(-x,-y,-z,ct)=(-\boldsymbol{r},ct).$ (18)

Finally, we define a 4-potential $A^{\mu}$ with entries

$A^{\mu}=(A_{x},A_{y},A_{z},\psi/c)=(\boldsymbol{A},\psi/c).$ (19)

We are now ready to employ some of this mathematical machinery. Consider the
_relativistic_ Lagrangian $L_{R}$ defined by the relation

$L_{R}=\tfrac{1}{2}mc\,(x^{\prime})^{\mu}(x^{\prime})^{\nu}g_{\mu\nu}+q\,(x^{\prime})^{\mu}A^{\nu}g_{\mu\nu}.$
(20)

It has the pleasing property that it is algebraically simple and treats space
and time on a similar footing. In particular, $L_{R}$ is evidently a Lorentz
scalar. That is, it is invariant under Lorentz transformations.333The quantity
$(x^{\prime})^{\mu}A^{\nu}g_{\mu\nu}$ is a scalar under Lorentz
transformations provided the 4-potential $A^{\nu}$ actually transforms as a
4-vector. See [LM:S.1.6fn]_LM_ , footnotes in Exercise 1.6.7. for a discussion
of the contrary case.

The _canonical_ momentum $p_{\mu}$ is given by the relation

$p_{\mu}=\partial{L_{R}}/\partial{(x^{\prime})^{\mu}}=mc\,(x^{\prime})_{\mu}+qA_{\mu},$
(21a) which can also be written in the form
$p_{\mu}=p^{\text{mech}}_{\mu}+qA_{\mu},$ (21b) where the _mechanical_
momentum is given by $p^{\text{mech}}_{\mu}=mc\,(x^{\prime})_{\mu}.$ (21c)

According to (21c), the mechanical momentum transforms like a 4-vector under
Lorentz transformations because $(x^{\prime})_{\mu}$ transforms like a
4-vector. From (21b) we see that the canonical momentum also transforms like a
4-vector to the extent that the 4-potential does so.444Again see [LM:S.1.6fn]
for a discussion of the contrary case.

Again implementing the standard procedure that relates Lagrangians and
Hamiltonians, one finds that the relativistic Hamiltonian $H_{R}$ associated
with the Lagrangian $L_{R}$ specified by (20) is given by the expression

$\displaystyle H_{R}$ $\displaystyle=$
$\displaystyle\tfrac{1}{2}mc\,(x^{\prime})^{\mu}(x^{\prime})^{\nu}g_{\mu\nu}=\tfrac{1}{2mc}(p^{\mu}-qA^{\mu})(p^{\nu}-qA^{\nu})g_{\mu\nu}$
(22) $\displaystyle=$
$\displaystyle\tfrac{1}{2mc}(p^{\mu}-qA^{\mu})(p_{\mu}-qA_{\mu}).$

Observe that $H_{R}$, like $L_{R}$, is Lorentz invariant. Note also that the
phase space associated with world-lines is eight-dimensional with canonical
coordinates $x^{\mu}$ and $p_{\nu}$.

Let us see what can be said about the phase-space trajectories generated by
$H_{R}$. Evidently $H_{R}$, as given by (22), does not depend explicitly on
$\tau$,

$\partial{H_{R}}/\partial{\tau}=0.$ (23)

It follows that $H_{R}$ is a constant (and integral) of motion. Moreover, from
(22), we see that the quantity
$\mathop{}\\!\mathrm{d}{s}^{2}/(\mathop{}\\!\mathrm{d}{\tau})^{2}$ defined by

$\mathop{}\\!\mathrm{d}{s}^{2}/(\mathop{}\\!\mathrm{d}{\tau})^{2}=g_{\mu\nu}(x^{\prime})^{\mu}(x^{\prime})^{\nu}=(x^{\prime})^{\mu}(x^{\prime})_{\mu}=x^{\prime}\cdot
x^{\prime}$ (24a) is a constant (and integral) of motion.
$\mathop{}\\!\mathrm{d}{s}^{2}/(\mathop{}\\!\mathrm{d}\tau)^{2}=\text{constant}.$
(24b)

Suppose we restrict our attention to those solutions that satisfy the relation

$x^{\prime}\cdot x^{\prime}=1.$ (25)

From (21c) and (25) we see that for these solutions $(p^{\text{mech}})^{\mu}$
satisfies the _mass-shell condition_

$p^{\text{mech}}_{\mu}(p^{\text{mech}})^{\mu}=(p^{\text{mech}})\cdot(p^{\text{mech}})=m^{2}c^{2}.$
(26)

From (22) and (25), we find that for these solutions $H_{R}$ has the value

$H_{R}=(mc/2).$ (27)

Moreover, for those solutions that satisfy (25), we have the result
$\mathop{}\\!\mathrm{d}{s}^{2}>0$ and may therefore select, in accord with
(16), (24a), and (25), a parameterization such that

$\mathop{}\\!\mathrm{d}{s}/\mathop{}\\!\mathrm{d}{\tau}=1.$ (28)

We have introduced three Hamiltonians, namely $H$, $K$, and $H_{R}$. It can be
shown that they all describe the same physics, and in this sense are
equivalent.555When using $H_{R}$, we employ only solutions that obey (25).
Which is to be employed depends on context. In what follows, we will use the
Hamiltonian $K$ with some reference to the Hamiltonian $H_{R}$.

##  3 Transfer Maps, the Symplectic Condition, and Symplectic Integrators

To proceed, it is convenient to introduce some terminology and definitions.
Suppose we are working with a $2n$-dimensional phase space. Let the symbol $z$
denote the collection of canonical phase-space variables arranged in the form

$z=(z_{1},z_{2},\dotsc,z_{2n})=(q_{1},\dotsc,q_{n};p_{1},\dotsc,p_{n}).$ (29)

Let $f(z,t)$ and $g(z,t)$ be any two (possibly “time-dependent”) functions of
$z$.666Here, by the “time”, we mean whatever has been selected to be the
independent variable. Note also that the symbol $z$ now no longer refers to
the third component of $\boldsymbol{r}$, but rather to the collection of
phase-space variables. Define their _Poisson bracket_ , $[f,g]$, by the rule

$[f,g]=\sum_{j}\frac{\partial{f}}{\partial{q_{j}}}\frac{\partial{g}}{\partial{p_{j}}}-\frac{\partial{f}}{\partial{p_{j}}}\frac{\partial{g}}{\partial{q_{j}}}.$
(30)

From this definition one may compute the _fundamental_ Poisson brackets

$[z_{a},z_{b}]=J_{ab},$ (31)

where $J$ is the $2n\times 2n$ matrix

$J=\begin{pmatrix}0&I\\\ -I&0\end{pmatrix}.$ (32)

Here $0$ and $I$ denote $n\times n$ zero and identity blocks, respectively.
The matrix $J$ is sometimes called the _Poisson matrix_.

There is an _existence and uniqueness theorem_ to the effect that sets of
first-order ordinary differential equations have solutions, and each solution
is uniquely specified by its initial conditions [LM:S.1.3]_LM_ , Section 1.3..
Hamilton’s equations of motion are first order. Now suppose a charged particle
enters a beam-line element (or collection of beam-line elements) with
_initial_ conditions $\boldsymbol{z}^{i}$ and subsequently exits with _final_
conditions $\boldsymbol{z}^{f}$. Then, by the existence and uniqueness
theorem, $\boldsymbol{z}^{f}$ is uniquely specified by $\boldsymbol{z}^{i}$.
Thus, there is a map $\mathcal{M}$, called a _transfer_ map, that sends
$\boldsymbol{z}^{i}$ to $\boldsymbol{z}^{f}$, and we write

$\boldsymbol{z}^{f}=\mathcal{M}\boldsymbol{z}^{i}.$ (33)

This relation between $\boldsymbol{z}^{i}$ and $\boldsymbol{z}^{f}$ is
illustrated by the picture shown in figure 1.

Figure 1: The transfer map $\mathcal{M}$ sends initial conditions
$\boldsymbol{z}^{i}$ to final conditions $\boldsymbol{z}^{f}$.

Next suppose small changes $\mathop{}\\!\mathrm{d}\boldsymbol{z}^{i}$ are made
in the initial conditions. The result will be associated small changes
$\mathop{}\\!\mathrm{d}\boldsymbol{z}^{f}$ in the final conditions. These
small changes will be connected by the relations

$\mathop{}\\!\mathrm{d}\boldsymbol{z}^{f}=M\mathop{}\\!\mathrm{d}\boldsymbol{z}^{i},$
(34)

where $M$ is the _Jacobian_ matrix with entries

$M_{ab}=\partial{z^{f}_{a}}/\partial{z^{i}_{b}}.$ (35)

It can be shown that if $\mathcal{M}$ is the result of integrating Hamilton’s
equations of motion, then its associated Jacobian matrix $M$ will satisfy the
condition

$M^{T}JM=J,$ (36)

where $M^{T}$ denotes the transpose of $M$. A matrix that satisfies (36) is
said to be _symplectic_ ; correspondingly $\mathcal{M}$ is called a
_symplectic map_.777In what follows, the letters $sp$ and $Sp$ are used as
abbreviations for _symplectic_. Note that in general $M$ depends on
$\boldsymbol{z}^{i}$. But $J$ does not. Therefore (36), since it must hold for
all $\boldsymbol{z}^{i}$, places strong (linear and nonlinear) restrictions on
$\mathcal{M}$ [LM:S.6.4.1]_LM_ , Subsection 6.4.1..

Suppose the final conditions $\boldsymbol{z}^{f}$ are to be determined by
integrating Hamilton’s equations of motion _numerically_ , and we also wish to
satisfy (36). A numerical integrator with this property is called a
_symplectic integrator_. In the case of a storage ring, such as the LHC, we
must integrate through thousands of beam-line elements to integrate through
even a single turn, and we wish to integrate through a large number of turns.
Therefore, even if we wish to integrate for only a small number of turns, we
would like to be able to use an _explicit_ symplectic integrator because the
numerous iterations required for an _implicit_ symplectic integrator would
make the computation extremely slow. Here we envision that one must iterate
the implicit solve to the point where convergence has been achieved to machine
precision in order to achieve symplecticity to machine precision. Our concern
is that the effort involved in iterating implicit Runge-Kutta or implicit
Gauss to machine precision will exceed the effort required for one step of
some explicit symplectic method for $K$ or $H_{R}$ if such an explicit
symplectic method can be found.

Explicit symplectic integrators are available if the Hamiltonian has the form
$T(p)+V(q)$. The Hamiltonian $K$ given by (10) is of this form if we make the
approximations

$\psi=0,\ A_{x}=0,\text{ and }A_{y}=0,$ (37)

for it then takes the form

$K=-\bigl{(}p_{t}^{2}/c^{2}-m^{2}c^{2}-p_{x}^{2}-p_{y}^{2}\bigr{)}^{1/2}-qA_{z}.$
(38)

However, it can be shown that this approximation excludes magnetic fringe-
field effects, which, by the Maxwell equations, _must_ occur at entry and exit
of every magnetic beamline element. It also excludes transverse electric
fringe fields, which must occur at entry and exit of all rf cavities. To make
accurate calculations that include fringe-field effects (which can be
important when nonlinear and even some linear effects are considered), one
would like to have an _explicit_ symplectic integrator that does not make the
approximations in (37).

Remarkably, there is an explicit symplectic integrator for the Hamiltonian
$H_{R}$ [LM:S.12.9]_LM_ , Section 12.9.. But there is a caveat: One of the
advertised features of symplectic integrators is that they can be used with a
rather large step size (thereby reducing computation time) since they at least
exhibit the qualitative nature of solutions exactly. However, there is a
theorem to the effect that (for any finite step size) symplectic integrators
do not preserve the Hamiltonian, even if it has no explicit dependence on the
independent variable ($\tau$ in the case of $H_{R}$) [8]. This may not be a
serious problem in some applications of symplectic integrators. In the case of
$H_{R}$, however, (27) tells us that failure to preserve $H_{R}$ means the
particle mass is _not_ preserved. Therefore, to preserve the particle mass to
good accuracy, which would seem highly desirable, it is necessary to employ a
sufficiently small step size, thereby making symplectic integration in this
situation relatively slow.

The rest of this chapter is devoted to exploring other possible approaches to
satisfying the symplectic condition (36) while at the same time achieving
improved computational speed and taking into account, through some desired
order, all linear and nonlinear effects associated with realistic
electromagnetic fields, including fringe fields and high-order multipole
fields.

##  4 Lie Algebraic Concepts and Tools

According to (30), Poisson brackets obey the _antisymmetry_ property

$[g,f]=-[f,g].$ (39)

It can be verified that Poisson brackets also satisfy the _Jacobi identity_.
Let $f$, $g$, and $h$ denote any three functions on phase-space. Then there is
the identity

$[f,[g,h]]+[g,[h,f]]+[h,[f,g]]=0.$ (40)

As a consequence, the Poisson bracket satisfies all the requirements for a
_Lie product_. The set of all phase-space functions therefore constitutes a
Lie algebra with the Poisson bracket as the Lie product [LM:S.5.1]_LM_ ,
Section 5.1..

Given any function $f(z,t)$, define an associated _differential_ operator,
denoted by $\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$ and called a _Lie operator_ , by the rule

$\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits=\sum_{j}\frac{\partial{f}}{\partial{q_{j}}}\frac{\partial{\
}}{\partial{p_{j}}}-\frac{\partial{f}}{\partial{p_{j}}}\frac{\partial{\
}}{\partial{q_{j}}}.$ (41)

Then, if $g(z,t)$ is any other phase-space function, the action of
$\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$ on $g$ is defined by writing

$\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits
g=\sum_{j}\frac{\partial{f}}{\partial{q_{j}}}\frac{\partial{g}}{\partial{p_{j}}}-\frac{\partial{f}}{\partial{p_{j}}}\frac{\partial{g}}{\partial{q_{j}}}=[f,g].$
(42)

Thus, a Lie operator may be viewed as a Poisson bracket waiting to happen.

In general, Lie operators do not commute. However, the commutator
$\\{\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits,\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits\\}$ of any two Lie operators
$\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$ and $\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$ is again a Lie operator. Indeed, as a
consequence of the Jacobi identity (40), the commutator
$\\{\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits,\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits\\}$ may be written in terms of the Poisson
bracket of the two underlying functions $f$ and $g$ according to the relation

$\\{\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits,\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits\\}=\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits-\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits=\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{[f,g]}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits.$ (43)

The relation (42) defines the action of $\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$. Powers of $\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$ can be defined by the rules

$\displaystyle\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits^{0}={\cal{I}}\Leftrightarrow\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits^{0}g=g,$ (44a)
$\displaystyle\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits^{1}g=[f,g],$ (44b)
$\displaystyle\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits^{2}g=[f,[f,g]],$ (44c)

and so on. Now that powers of the Lie operator $\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$ have been defined, one may also define power
series in $\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$. Of particular interest is the power series
associated with the exponential function by the rule

$\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)=\mathrm{e}^{\mspace{1.5mu}\mspace{-1.5mu}\mathop{\text{\raise
0.2411pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.2411pt\hbox{{:}}}}\nolimits}=\sum_{m=0}^{\infty}\frac{1}{m!}\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits^{m}.$ (45)

This operator $\mathrm{e}^{\mspace{1.5mu}\mspace{-1.5mu}\mathop{\text{\raise
0.2411pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.2411pt\hbox{{:}}}}\nolimits}$, called a _Lie transformation_ , therefore
acts on $g$ according to the relation

$\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)g=g+[f,g]+\frac{1}{2!}[f,[f,g]]+\dotsb.$ (46)

In this context $\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$ (and sometimes $f$) is called a _Lie
generator_. Its importance for us lies in the fact that _any Lie
transformation generates a symplectic map_.

Poisson brackets are invariant under symplectic maps [LM:S.6.1.2]_LM_ ,
Section 6.1.2.. An important consequence of this fact is the extremely useful
similarity relation [LM:S.8.2]_LM_ , Section 8.2.

$\mathcal{L}\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\mathcal{L}^{-1}=\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{\mathcal{L}f}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits),$ (47)

where $\mathcal{L}$ denotes any symplectic map.

Suppose the transfer map $\mathcal{M}$ has the property that it maps the
origin into itself. In other words, we assume that $\mathcal{M}z$ has a Taylor
expansion of the form

$z_{a}^{f}=\sum_{b}R_{ab}z_{b}^{i}+\sum_{bc}T_{abc}z_{b}^{i}z_{c}^{i}+\sum_{bcd}U_{abcd}z_{b}^{i}z_{c}^{i}z_{d}^{i}+\dotsb,$
(48)

which has no constant term. This can always be accomplished by the use of
_deviation_ variables. If $\mathcal{M}$ is symplectic, then $R$ must be a
symplectic matrix. Moreover, the Taylor coefficients $T$, $U$, $\cdots$ cannot
be arbitrary, but are constrained by complicated nonlinear relations that
follow from the symplectic condition (36). A _truncated_ Taylor expansion of a
symplectic map is called a symplectic _jet_. Finally, the series (48) _cannot_
in general be truncated without violating the symplectic condition. Therefore
a symplectic jet is generally not a symplectic map.

However, there is a _factorization theorem_ [LM:S.7.6]_LM_ , Section 7.6. to
the effect that $\mathcal{M}$ can also be written in the Lie product form

$\mathcal{M}=\mathcal{R}\,[\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{3}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{4}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\dotsb]\,\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{1}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits).$ (49)

Here $\mathcal{R}$ is the linear symplectic map associated with $R$, and the
$f_{m}$ are _homogeneous_ polynomials of degree $m$. Translations/deviations
from the origin, which correspond to constant terms added if desired to the
Taylor series (48), are described by $f_{1}$; and the $f_{3}$, $f_{4}$
$\cdots$ describe the nonlinear terms in (48). Unlike the Taylor coefficients,
there are no restrictions imposed on the $f_{m}$ by the symplectic condition.

Any analytic symplectic map is uniquely specified by a first-order polynomial
$f_{1}$, a symplectic matrix $R$, together with a collection of homogeneous
polynomials $f_{3}$, $f_{4}$, … that describe the nonlinear part of the map.
And the converse also holds. In addition, the (in principle infinite) product
appearing in square brackets in (49) can be truncated at any stage without
violating the symplectic condition. It can be shown that each factor in (49)
is a symplectic map, and the product of any number of symplectic maps is also
a symplectic map. In what follows it is also convenient to write (49) in the
form

$\mathcal{M}=\mathcal{R}\mathcal{N}\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{1}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)$ (50)

where $\mathcal{N}$, the _nonlinear_ part of $\mathcal{M}$, is given by

$\mathcal{N}=\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{3}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{4}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\cdots.$ (51)

##  5 Symplectic Completion of Symplectic Jets

We know that the Lie algebra of all Lie operators, which we will call
$\operatorname{ispm}(2n,\mathbb{R})$, is infinite dimensional. Correspondingly
$\operatorname{ISpM}(2n,\mathbb{R})$, the group of symplectic maps, is
infinite dimensional.888The letters $m$ and $M$ that appear in this and the
previous sentence are abbreviations for $map$. The letters $i$ and $I$ are
abbreviations for $inhomogeneous$. By $inhomogeneous$ it is meant that the
possibility of constant terms appearing in (48) is included. Indeed, the
factorization (50) gives a representation of the general analytic symplectic
map. We see that the specification of a symplectic map generally requires an
infinite number of parameters. This fact produces an awkward situation for
human beings and computers, which can work only with a finite number of
quantities (and often only with finite precision).

An optimistic perspective on the experimental and theoretical situation, for
example in the field of accelerator physics, might be stated as follows: We
know that a beam transport system, accelerator, storage ring, or any portion
thereof may be described by a symplectic transfer map. However, because we
cannot measure or control electromagnetic fields exactly, we are unsure of and
unable to control exactly what this map is. Also, since it is impossible to
perform computations with an infinite number of parameters/variables and to
infinite precision, it is necessary to develop various approximation schemes.
Thus, we are able to study computationally (and probably theoretically) the
detailed properties of only a subset of all symplectic maps. The hope is that
if two symplectic maps are in some sense nearly the same, then their behavior
[including, in some cases, long-term (repeated iteration) behavior] will be in
some important ways nearly the same.999Note that a similar optimism is shared
by practitioners of symplectic integration. Were that not true from an
experimental standpoint, it would be impossible to build satisfactory storage
rings and the like. Were that not true from a theoretical standpoint, it would
be impossible to design storage rings and the like with any assurance of
satisfactory performance.

Suppose, as an approximation, the product appearing in the square brackets of
(49), namely the map $\mathcal{N}$, is _truncated_ at $m=\text{maxm}$ to
produce the map $\mathcal{N}^{\text{trunc}}$ given by

$\mathcal{N}^{\text{trunc}}=\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{3}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{4}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\dotsb\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{\text{maxm}}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits).$ (52)

In analogy with (50) we also make the definition

$\mathcal{M}^{\rm{trunc}}=\mathcal{R}\mathcal{N}^{\text{trunc}}\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{f_{1}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits).$ (53)

The map $\mathcal{M}^{\rm{trunc}}$, while exactly symplectic, requires only a
finite number of parameters for its specification. We only need store a first-
order polynomial $f_{1}$, a symplectic matrix $R$, and a collection of
homogeneous polynomials $f_{3}$, $f_{4}$, …$f_{\rm{maxm}}$. For example, in
the case of a 6-dimensional phase space, $3002$ parameters are required when
$\text{maxm}=8$.

What can be said about the _accuracy_ of $\mathcal{M}^{\text{trunc}}$? Suppose
the Taylor expansion (48) is terminated by retaining only terms through degree
$(\text{maxm}-1)$, thereby forming a symplectic jet that we will call
$\mathcal{J}$. An examination of the proof of the factorization theorem (49)
shows that a knowledge of the coefficients in $\mathcal{J}$ supplies just
enough information to determine the ingredients of
$\mathcal{M}^{\text{trunc}}$, and vice versa. In other words, a knowledge of
$\mathcal{J}$ supplies just enough information to determine $\mathcal{R}$, the
polynomial $f_{1}$, and the polynomials $f_{3}$ through $f_{\text{maxm}}$, and
vice versa. (In particular, $f_{3}$ contributes only to $T$ terms and to terms
beyond second order, $f_{4}$ contributes only to $U$ terms and to terms beyond
third order, _etc._) Thus, while exactly symplectic, the map
$\mathcal{M}^{\rm{trunc}}$ is guaranteed accurate only through terms of degree
$(\text{maxm}-1)$. With regard to memory requirements, the storage of a jet
requires more locations because it does not exploit the symplectic condition.
For example, $10\,296$ locations are required in the case of 6-dimensional
phase space when $\text{maxm}=8$.

At present there are Lie-algebraic results and Truncated Power Series (TPSA)
methods that make it possible to compute in principle the ingredients in
$\mathcal{M}^{\rm{trunc}}$ as given by (53), with $\text{maxm}=8$, for any
beam-line element or collection of beam-line elements (including a full ring)
based on field data provided numerically on a grid [LM:S.17ff]_LM_ , Chapters
17–25..

What has been accomplished here? Given a symplectic jet $\mathcal{J}$, we have
found a map $\mathcal{M}^{\text{trunc}}$ that is guaranteed symplectic and
whose Taylor expansion agrees with $\mathcal{J}$ through terms of degree
$(\text{maxm}-1)$. This Taylor expansion will in general contain terms of
degrees beyond $(\text{maxm}-1)$. It will be called the _Lie symplectic
completion_ of $\mathcal{J}$.

But there is a problem: Suppose we wish to evaluate
$\mathcal{M}^{\text{trunc}}\boldsymbol{z}^{i}$ for some initial condition
$\boldsymbol{z}^{i}$. Because of the infinite series that appears in the
definition (45) of a Lie transformation, the symplectic completion of
$\mathcal{J}$ provided by $\mathcal{M}^{\text{trunc}}$ will in general contain
an _infinite_ number of terms. As a consequence, the evaluation of
$\mathcal{M}^{\text{trunc}}\boldsymbol{z}^{i}$ will generally involve the
summation of infinite series, a task that generally lies beyond numerical
methods, or is at best numerically intensive unless the series converges
rapidly.

We have seen that in principle the symplectic completion of a symplectic jet
is possible. What we would like are other symplectic jet completions whose
actions on $\boldsymbol{z}^{i}$ can be computed rapidly and to machine
precision. Two such methods will be described in the next two sections of this
chapter.

##  6 Symplectic Completion Using Generating Functions

We have described how in general the computation of the action of a Lie
transformation on phase space involves the summation of an _infinite_ series
if the symplectic condition is to be honored. Before exploring a particular
method to deal with this problem, we begin this section by studying a simple
example of what happens if only a finite number of terms in the series
expansion are employed.

As such an example consider, for a two-dimensional phase space, the symplectic
map $\mathcal{M}$ given by the relation

$\mathcal{M}=\mathcal{R}\mathcal{N}$ (54)

with linear part

$\mathcal{R}=\exp(-\tfrac{\theta}{2}\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{p^{2}+q^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)$ (55)

and nonlinear part

$\mathcal{N}=\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{qp^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits).$ (56)

This map may be viewed as a toy model for the one-turn map of a storage ring.

In this case the infinite series for $\mathcal{R}$ can be summed exactly to
give the results

$Q=\mathcal{R}q=q\cos\theta+p\sin\theta,\quad
P=\mathcal{R}p=-q\sin\theta+p\cos\theta.$ (57)

And the map $\mathcal{N}$ can also be evaluated exactly, since the exponential
of any monomial Lie operator can be evaluated exactly [LM:S.monom]_LM_ ,
Subsections 1.2.3, 1.4.1, 1.4.2, and Exercise 1.4.3.. For the case at hand
there is the result

$Q=\mathcal{N}q=q(1-p)^{2},\quad P=\mathcal{N}p=p/(1-p).$ (58)

Therefore $\mathcal{M}$ can also be evaluated exactly. Note that the result
(58) for $P$, and therefore the map $\mathcal{N}$, has a pole on the phase-
space surface $p=1$.

Figure 2 shows the result of applying $\mathcal{M}$ repeatedly to seven
initial conditions for the case $\theta/2\pi=0.22$. In other words, seven
initial conditions have been selected, and their orbits have been found under
the repeated action of $\mathcal{M}$. One initial condition lies very near the
origin, and its orbit appears to lie on a closed curve that is nearly
elliptical. (It would be nearly circular had the horizontal and vertical
scales been equal.) This is to be expected because the effect of the nonlinear
part $\mathcal{N}$ is so small on such orbits that they are essentially those
of the rotation map $\mathcal{R}$. By contrast, the other initial conditions
lie successively farther from the origin, and the effect of $\mathcal{N}$
becomes ever more significant. Their orbits appear to lie on closed curves
that, the farther they lie from the origin, the more noticeably do
nonlinearities distort them from being circular. The origin itself is an
elliptic fixed point corresponding to a one-turn closed orbit.

Figure 2: Phase-space portrait, in the case $\theta/2\pi=0.22$, resulting from
applying the map $\mathcal{M}$ repeatedly (2000 times) to the seven initial
conditions $(q,p)=(0.01,0)$, $(0.1,0)$, $(0.15,0)$, $(0.2,0)$, $(0.25,0)$,
$(0.3,0)$, and $(0.35,0)$ to find their orbits.

Now suppose the nonlinear map $\mathcal{N}$ is _truncated_ to form the map
$\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$ by retaining only the first _two_
terms in its Taylor expansion. In Lie form we have the result

$\mathcal{N}^{\mspace{1.0mu}\text{tr2}}=\mathcal{I}+\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{qp^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits.$ (59)

This truncated map $\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$ has the effect

$\displaystyle Q$
$\displaystyle=\mathcal{N}^{\mspace{1.0mu}\text{tr2}}q=(\mathcal{I}+\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{qp^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)q=q+[qp^{2},q]=q-2qp\Leftrightarrow Q-q=-2qp,$
(60a) $\displaystyle P$
$\displaystyle=\mathcal{N}^{\mspace{1.0mu}\text{tr2}}p=(\mathcal{I}+\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{qp^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)p=p+[qp^{2},p]=p+p^{2}\Leftrightarrow
P-p=p^{2}.$ (60b)

Evidently, $\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$ is a symplectic jet map
that retains only terms through degree 2. Indeed, one finds using (60) the
result

$[Q,P]=1-4p^{2}\neq 1,$ (61)

and therefore $\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$, while a symplectic jet
map, is as expected _not_ a symplectic map.101010It follows from (36) that a
symplectic map must preserve Poisson brackets, and vice versa.

Next define a corresponding map $\mathcal{M}^{\text{tr2}}$ by writing

$\mathcal{M}^{\text{tr2}}=\mathcal{R}\mathcal{N}^{\mspace{1.0mu}\text{tr2}}.$
(62)

The left-hand graphic in Figure 3 shows the orbits of
$\mathcal{M}^{\text{tr2}}$ for two initial conditions: one near the origin,
and one quite far away. Inspection of the figure shows that orbits are no
longer distorted circles, but instead appear to spiral into the origin. This
motion toward the origin occurs because
$\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$, and consequently
$\mathcal{M}^{\text{tr2}}$, is _not_ symplectic.

We could also retain the next term in the Taylor series for
$\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{qp^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)$ to form the symplectic jet map
$\mathcal{N}^{\mspace{1.0mu}\text{tr3}}$,

$\mathcal{N}^{\mspace{1.0mu}\text{tr3}}=\mathcal{I}+\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{qp^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits+\tfrac{1}{2}\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{qp^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits^{2}.$ (63)

It retains terms through degree 3. It is still nonsymplectic, but more nearly
symplectic than $\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$. Again define a
corresponding map $\mathcal{M}^{\text{tr3}}$ by writing

$\mathcal{M}^{\text{tr3}}=\mathcal{R}\mathcal{N}^{\mspace{1.0mu}\text{tr3}}.$
(64)

The phase-space portrait for $\mathcal{M}^{\text{tr3}}$ is found to be
somewhat more like that for $\mathcal{M}$ than that provided by
$\mathcal{M}^{\text{tr2}}$, because $\mathcal{M}^{\text{tr3}}$ is more nearly
symplectic. However, there is still substantial/disastrous nonsymplectic
spiraling, in this case out of the origin. We see that violation of the
symplectic condition can convert the origin, initially an elliptic fixed
point, into a nonlinear attracting or a nonlinear repelling fixed point.

Suppose we know that the behavior of some dynamical system is describable by a
symplectic map. This system may be a beam-line element, some collection of
beam-line elements, or even the one-turn map for a storage ring. And suppose
that a truncated Taylor expansion (symplectic jet $\mathcal{J}$) is known for
this map through terms of some order. As stated earlier, such knowledge is in
fact computable using Lie algebraic and TPSA algorithms [LM:S.10.5]_LM_ ,
Section 10.5 and Chapter 39. Some authors refer to TPSA as _Differential
Algebra_ (DA). For an exposition of DA, see M. Berz, _Modern Map Methods in
Particle Beam Physics_ , vol. 108 of _Advances in Imaging and Electron
Physics_ , Academic Press, San Diego, 1999.. According to the factorization
theorem, as far as nonlinear effects are concerned, such knowledge is
equivalent to the knowledge of a set of homogeneous polynomials $f_{3}$, …,
$f_{n}$. What we would like to find is a map that is symplectic, has the jet
$\mathcal{J}$ through terms of degree $(n-1)$, and is relatively easy to
compute. Because it is symplectic, its Taylor expansion must in general also
have terms beyond degree $(n-1)$. In some way that is not yet clearly
described, we would like these additional terms to be as small as possible
while remaining consistent with the symplectic condition. For example, their
effect over the phase-space region of interest should not be appreciably
larger than the extent to which $\mathcal{J}$ violates the symplectic
condition.

How can we find such symplectic maps? We need a supply of relatively easily
computed symplectic maps. It is known that such maps can be produced with the
aid of generating functions.111111Symplectic completion of symplectic jets
using a generating function was first implemented—in the context of
Accelerator Physics—in the Lie-algebra based accelerator design code MaryLie.
Four types of generating functions, commonly called $F_{1}(q,Q)$,
$F_{2}(q,P)$, $F_{3}(p,Q)$, and $F_{4}(p,P)$, are usually presented in
graduate Classical Mechanics text books. What is less familiar is that, for a
$2n$ dimensional phase space, there is in fact a $2n(4n+1)$ dimensional
_family_ of types of generating functions: There is a type for every $4n\times
4n$ symplectic matrix [LM:S.ntype]_LM_ , Section 6.7.. Among these types, we
have found the so called _Poincaré_ generating function, which we denote as
$F_{+}$, to be particularly attractive [LM:S.Poin]_LM_ , Chapter 34. See also
the article B. Erdélyi and M. Berz, “Optimal symplectic approximation of
Hamiltonian flows,” _Phys. Rev. Lett._ , vol. 87, 114302, Aug. 2001. Currently
we do not find persuasive their invocation of the Hofer metric, but do agree
(for other reasons) with their conclusion that use of the Poincaré generating
function has several desirable features..

To describe the use of $F_{+}$, it is useful to make some additional
definitions. In addition to the definition of $z$ given by (29), which
describes the collection of what we might call _old_ variables, we introduce
the symbol $Z$ to denote a collection of what we may call _new_ variables:

$Z=(Z_{1},Z_{2},\dotsc,Z_{2n})=(Q_{1},\dotsc,Q_{n};P_{1},\dotsc,P_{n}).$ (65)

We shall also need the sums and differences defined by

$\Sigma=Z+z\text{ and }\Delta=Z-z.\Leftrightarrow
z=\tfrac{1}{2}(\Sigma-\Delta)\text{ and }Z=\tfrac{1}{2}(\Sigma+\Delta).$ (66)

Specifically, for future use and in the case of a two-dimensional phase space,
the relations (66) take the equivalent forms

$\displaystyle\Sigma_{1}$ $\displaystyle=Q+q\text{ and
}\Delta_{1}=Q-q\Leftrightarrow q=\tfrac{1}{2}(\Sigma_{1}-\Delta_{1})\text{ and
}Q=\tfrac{1}{2}(\Sigma_{1}+\Delta_{1}),$ (67a) $\displaystyle\Sigma_{2}$
$\displaystyle=P+p\text{ and }\Delta_{2}=P-p\Leftrightarrow
p=\tfrac{1}{2}(\Sigma_{2}-\Delta_{2})\text{ and
}P=\tfrac{1}{2}(\Sigma_{2}+\Delta_{2}).$ (67b)

Finally, we will need a collection of $2n$ _auxiliary_ variables, which we
will call $u$:

$u=(u_{1},u_{2},\dotsc,u_{2n}).$ (68)

Now let $F_{+}(u,t)$ be any function of $u$ and perhaps the independent/time
variable $t$. We define its action on phase space by the rule

$\Delta=J\partial_{u}F_{+}|_{u=\Sigma}.$ (69)

Note that this definition specifies a relation between $\Delta$ and $\Sigma$,
which in turn, when the relations on the right-hand side of (66) are taken
into account, specifies a relation between $Z$ and $z$. By the general theory
of generating function machinery, this relation between $Z$ and $z$ is
guaranteed to be a symplectic map for _any_ choice of $F_{+}(u,t)$.

To see how this mathematical machinery works in some detail, let us apply it
to a simple example in two-dimensional phase-space. Suppose $F_{+}$ is the
cubic monomial

$F_{+}(u)=F^{3}_{+}(u)=-\tfrac{1}{4}u_{1}{u_{2}}^{2}.$ (70)

We then compute

$J\partial_{u}F^{3}_{+}=\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}\begin{pmatrix}-\tfrac{1}{4}{u_{2}}^{2}\\\\[3.01385pt]
-\tfrac{1}{2}u_{1}u_{2}\end{pmatrix}=\begin{pmatrix}-\tfrac{1}{2}u_{1}u_{2}\\\\[3.01385pt]
\,\tfrac{1}{4}{u_{2}}^{2}\end{pmatrix},$ (71)

and hence

$J\partial_{u}F^{3}_{+}|_{u=\Sigma}{\;}=\begin{pmatrix}-\tfrac{1}{2}\Sigma_{1}\Sigma_{2}\\\\[3.01385pt]
\tfrac{1}{4}{\Sigma_{2}}^{2}\end{pmatrix}.$ (72)

From (69) and (72) it follows that

$\Delta_{1}=-\tfrac{1}{2}\Sigma_{1}\Sigma_{2},\text{ and
}\Delta_{2}=\tfrac{1}{4}{\Sigma_{2}}^{2}.$ (73)

Finally, employ the relations (67) in (73) to obtain the relations

$\displaystyle Q-q=-\tfrac{1}{2}(Q+q)(P+p)$ $\displaystyle\Leftrightarrow
Q=q-\tfrac{1}{2}(Q+q)(P+p),$ (74a) $\displaystyle P-p=\tfrac{1}{4}(P+p)^{2}$
$\displaystyle\Leftrightarrow P=p+\tfrac{1}{4}(P+p)^{2}.$ (74b)

The two equations on the right-hand sides of (74) specify a relation between
$Z$ and $z$.

As is the case with the use of any generating function, the above relation
between $Z$ and $z$ is _implicit_. We can _begin_ to make it _explicit_ by
seeking a Taylor expansion using iteration: For the zeroth iteration, make the
Ansatz

$\displaystyle Q=q\text{ and }P=p.$ (75a) Now substitute this Ansatz into the
right-hand sides of (74) to yield for the first iteration the results
$\displaystyle Q=q-2qp\text{ and }P=p+p^{2}.$ (75b)

Observe that the results (75b) agree with the jet results (60). We have found
a symplectic map whose jet through terms of second order agrees with
$\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$.

Figure 3: Phase-space portraits, in the case $\theta/2\pi=0.22$. Left: The
resulting of applying the map $\mathcal{M}^{\text{tr2}}$ repeatedly ($1000$
times) to the two initial conditions $(q,p)=(0.01,0)$ and $(0.4,0)$. These
orbits spiral into the origin. Right: The result of applying the map
$\mathcal{M}^{\text{psc2}}$ repeatedly ($2000$ times) to the seven initial
conditions $(q,p)=(0.01,0)$, $(0.1,0)$, $(0.15,0)$, $(0.2,0)$, $(0.25,0)$,
$(0.3,0)$, and $(0.35,0)$. Light gray curves in the background indicate the
exact result.

It happens that for this example the implicit equations on the right-hand
sides of (74) can be made explicit by algebraic manipulation. The equation on
the right-hand side of (74b) is quadratic in $P$, and on choosing the negative
square root (the solution for which $p$ vanishing implies that $P$ also
vanishes), we obtain the relation

$P=-(p-2)-2\sqrt{1-2p}.$ (76a) Once $P$ is known, the equation on the right-
hand side of (74a) is linear in $Q$ and has the immediate solution
$Q=\frac{q\sqrt{1-2p}}{2-\sqrt{1-2p}}.$ (76b)

We see that, in this case, the use of $F_{+}$ produces a map that has a
square-root branch point on the phase-space surface $p=1/2$ and a pole at
$p=-3/2$. Note that the branch-point singularity lies closer to the origin
than does the pole of the exact map, which we have noted lies on the surface
$p=1$.

Let $\mathcal{N}^{\mspace{1.0mu}\text{psc2}}$ be the map given by (76). One
can verify by direct computation that it has the property

$[Q,P]=1,$ (77)

and hence the map $\mathcal{N}^{\mspace{1.0mu}\text{psc2}}$ is exactly
symplectic, as desired and expected. The reader can also verify that the terms
through degree $2$ in the Taylor expansions of (76) agree with the terms in
(60). This is just the result (75b) that we have already found by
iteration.121212Surprisingly, the terms through degree $3$ in the Taylor
expansions of (76) _agree_ with the terms generated by applying (63) to
$(q,p)$. This happens due to the second relation in (80) and our tacit
assumption that $f_{4}=0$. We may therefore say that
$\mathcal{N}^{\mspace{1.0mu}\text{psc2}}$ is the _Poincaré symplectic
completion_ of the degree-two symplectic jet map
$\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$ given by (60). Correspondingly,
suppose we define the associated map $\mathcal{M}^{\text{psc2}}$ by the
relation

$\mathcal{M}^{\text{psc2}}=\mathcal{R}\mathcal{N}^{\mspace{1.0mu}\text{psc2}}.$
(78)

We expect it to be symplectic because it is the product of two symplectic
maps.

The right-hand graphic in Figure 3 shows the result of applying
$\mathcal{M}^{\text{psc2}}$ repeatedly to seven initial conditions for the
case $\theta/2\pi=0.22$. Note that orbits generated by
$\mathcal{M}^{\text{psc2}}$ exhibit no spurious spiraling towards or away from
the origin. Moreover, comparison with the background light gray curves showing
the exact result reveals that the orbits closely agree in shape, but there is
some difference in phase advance. It seems remarkable that the relatively
meager information about $\mathcal{N}$ present in
$\mathcal{N}^{\mspace{1.0mu}\text{tr2}}$ suffices, after Poincaré symplectic
completion has been performed, to give such good agreement.

The reader may wonder how we knew to make the inspired choice (70) for
$F_{+}^{3}(u)$. Suppose we make for $F_{+}(u)$ the expansion

$F_{+}(u)=\sum_{m=3}^{\text{maxm}}F_{+}^{m}(u),$ (79)

where the $F_{+}^{m}(u)$ are homogeneous polynomials of degree $m$. Then there
are formulas that determine the $F_{+}^{m}(u)$ in terms of the $f_{m}(u)$. For
example, there are the relations

$F_{+}^{3}(u)=-\tfrac{1}{4}f_{3}(u)\text{ and
}F_{+}^{4}(u)=-\tfrac{1}{8}f_{4}(u).$ (80)

Observe that if we assume $F_{+}(u)$ has _only_ a $F_{+}^{3}(u)$ component,
then the associated map will have an $f_{3}$ given by (80) and a _vanishing_
$f_{4}$. (This is one of the virtues of the Poincaré generating function.) In
general the $f_{>4}$ will not vanish. Results for the $F_{+}^{m}(u)$ in terms
of the $f_{m}(u)$ are known through order $m=8$, but become increasingly
complicated as $m$ increases beyond $m=4$ [LM:S.34.4.]_LM_ , Section 34.4..
However, it is also possible to proceed _without_ these formulas, thereby
bypassing their complications. We will next illustrate how to do so for our
simple example.

Suppose only the functions $f_{3}$, $f_{4}$, …, $f_{\text{maxm}}$ are known.
For our example we know from (56) that $f_{3}=qp^{2}$ and $(\text{maxm}-1)=2$.
Next compute the set of Taylor series through terms of degree
$(\text{maxm}-1)$ for the jet
$\mathcal{N}^{\mspace{1.0mu}\text{tr(\text{maxm}-1)}}$. For our example this
set is given by (60). In this set, replace $Q,q$ and $P,p$ by their
representations in terms of $\Sigma$ and $\Delta$ using, for this example, the
relations given on the right-hand sides of (67). So doing, for the relations
appearing on the far-right sides of (60), yields the results

$\displaystyle\Delta_{1}$
$\displaystyle=-2\,\bigl{[}\tfrac{1}{2}(\Sigma_{1}-\Delta_{1})\,\tfrac{1}{2}(\Sigma_{2}-\Delta_{2})\bigr{]}=-\tfrac{1}{2}(\Sigma_{1}-\Delta_{1})(\Sigma_{2}-\Delta_{2}),$
(81a) $\displaystyle\Delta_{2}$
$\displaystyle=\bigl{[}\tfrac{1}{2}(\Sigma_{2}-\Delta_{2})\bigr{]}^{2}=\tfrac{1}{4}(\Sigma_{2}-\Delta_{2})^{2}.$
(81b)

Approximately solve these equations for the quantities $\Delta$ in terms of
the quantities $\Sigma$. Do so in the form of a Taylor series in $\Sigma$
truncated beyond terms of degree $(\text{maxm}-1)$, which can be done by
iteration: For the zeroth iteration make the Ansatz

$\displaystyle\Delta_{1}=0$ $\displaystyle\text{ and }\Delta_{2}=0.$ (82a) Now
substitute this Ansatz into the right-hand sides of (81) to obtain, through
terms of degree $2$, the Taylor expansion
$\displaystyle\Delta_{1}=-\tfrac{1}{2}\Sigma_{1}\Sigma_{2}$
$\displaystyle\text{ and }\Delta_{2}=\tfrac{1}{4}{\Sigma_{2}}^{2}.$ (82b)

In this case the iteration process is finished because
$(\text{maxm}-1)=3-1=2.$ Observe that the relations (82b) agree with the
relations (73)! The transition from (73) to the right-hand sides of (74) now
proceeds as before. We have found, _directly_ from the jet
$\mathcal{N}^{\mspace{1.0mu}\text{tr(maxm-1)}}$, the relations that would have
flowed from the use of the related $F_{+}$. And we have therefore identified
(in implicit form) the desired symplectic map. Note also that the operations
we have just performed involve only well-defined polynomial manipulations, and
therefore can be performed on a computer using TPSA routines.

There remains the problem of converting the implicit results given by (69) to
explicit results for $Z$ in terms of $z$. For our simple example we were able
to do so by solving a quadratic equation. For most applications, however, the
equations to be solved are much more complicated, and must be handled
numerically. That is, given $z$ as a collection of numbers, we would like to
find the associated collection of numbers $Z$. This can be done by simple
iteration or by use of Newton’s method. In either case the process can be
started using jet results.

As an example of the use of simple iteration, the equations on the right-hand
sides of (74) may be converted into the iteration rule

$\displaystyle Q^{[n+1]}$
$\displaystyle=q-\tfrac{1}{2}(Q^{[n]}+q)(P^{[n]}+p),$ (83a) $\displaystyle
P^{[n+1]}$ $\displaystyle=p+\tfrac{1}{4}(P^{[n]}+p)^{2}.$ (83b)

Using the jet results (60), one may begin the iteration optimally with the
values

$Q^{[0]}=q-2qp\text{ and }P^{[0]}=p+p^{2}.$ (84)

Observe that in general the quantities to be evaluated numerically at each
step are polynomials, and therefore this evaluation is quite fast.

For example, consider the case

$q=-0.3\text{ and }p=-0.2,$ (85)

which is a point near the boundary of figure 2. Then use of (76) shows that in
this case we hope to find the results

$Q^{\infty}=$-0.434\,588\,297\,681\,520\,63$\cdots\text{ and
}P^{\infty}=$-0.166\,431\,913\,239\,846\,35$\cdots.$ (86)

Table 1: Convergence of $Q^{[n]},P^{[n]}$ and $\Delta Q^{[n]},\Delta P^{[n]}$ as a function of $n$ using simple iteration. $n$ | $Q^{[n]}$ | $P^{[n]}$ | $\Delta Q^{[n]}$ | $\Delta P^{[n]}$
---|---|---|---|---
$0$ | $-0.420\,000\,000\,000\,000\,0$ | $-0.160\,000\,000\,000\,000\,0$ | $-1.459\text{\times}{10}^{-2}$ | $-6.432\text{\times}{10}^{-3}$
$1$ | $-0.429\,599\,999\,999\,999\,9$ | $-0.167\,600\,000\,000\,000\,0$ | $-4.988\text{\times}{10}^{-3}$ | $1.168\text{\times}{10}^{-3}$
$2$ | $-0.434\,100\,480\,000\,000\,0$ | $-0.166\,217\,560\,000\,000\,0$ | $-4.878\text{\times}{10}^{-4}$ | $-2.144\text{\times}{10}^{-4}$
$3$ | $-0.434\,420\,243\,290\,214\,4$ | $-0.166\,471\,174\,686\,911\,6$ | $-1.681\text{\times}{10}^{-4}$ | $3.926\text{\times}{10}^{-5}$
⋮ | ⋮ | ⋮ | ⋮ | ⋮
$17$ | $-0.434\,588\,297\,681\,512\,6$ | $-0.166\,431\,913\,239\,848\,3$ | $-7.994\text{\times}{10}^{-15}$ | $1.943\text{\times}{10}^{-15}$
$18$ | $-0.434\,588\,297\,681\,520\,0$ | $-0.166\,431\,913\,239\,846\,1$ | $-6.661\text{\times}{10}^{-16}$ | $-2.776\text{\times}{10}^{-16}$
$19$ | $-0.434\,588\,297\,681\,520\,4$ | $-0.166\,431\,913\,239\,846\,5$ | $-1.665\text{\times}{10}^{-16}$ | $1.388\text{\times}{10}^{-16}$
$20$ | $-0.434\,588\,297\,681\,520\,7$ | $-0.166\,431\,913\,239\,846\,4$ | $5.551\text{\times}{10}^{-17}$ | $5.551\text{\times}{10}^{-17}$
Table 2: Convergence of $Q^{[n]},P^{[n]}$ and $\Delta Q^{[n]},\Delta P^{[n]}$ as a function of $n$ using Newton’s method. $n$ | $Q^{[n]}$ | $P^{[n]}$ | $\Delta Q^{[n]}$ | $\Delta P^{[n]}$
---|---|---|---|---
$0$ | $-0.420\,000\,000\,000\,000\,0$ | $-0.160\,000\,000\,000\,000\,0$ | $-1.459\text{\times}{10}^{-2}$ | $-6.432\text{\times}{10}^{-3}$
$1$ | $-0.434\,534\,931\,789\,995\,8$ | $-0.166\,440\,677\,966\,101\,7$ | $-5.337\text{\times}{10}^{-5}$ | $8.765\text{\times}{10}^{-6}$
$2$ | $-0.434\,588\,297\,975\,149\,4$ | $-0.166\,431\,913\,256\,077\,6$ | $2.936\text{\times}{10}^{-10}$ | $1.623\text{\times}{10}^{-11}$
$3$ | $-0.434\,588\,297\,681\,520\,8$ | $-0.166\,431\,913\,239\,846\,4$ | $1.110\text{\times}{10}^{-16}$ | $8.327\text{\times}{10}^{-17}$
$4$ | $-0.434\,588\,297\,681\,520\,8$ | $-0.166\,431\,913\,239\,846\,4$ | $1.110\text{\times}{10}^{-16}$ | $5.551\text{\times}{10}^{-17}$

For the case (85), table 1 displays simple iteration results near the
beginning and end of the iteration process. We see that the point
$Q^{\infty},P^{\infty}$ is, as desired, an _attractor_ for the simple
iteration process. Convergence to machine precision has been achieved after
about 20 iterations. However, the simple iteration process converges rather
slowly.131313The iteration process does converge more rapidly for $q,p$ values
that lie closer to the origin, but the improvement is not impressive. For
example, after reducing the distance from the origin by a factor of $5$, at
least $12$ iterations are still required to achieve convergence to machine
precision. Let us examine the convergence rate. Define errors $\Delta
Q^{[n]},\Delta P^{[n]}$ by the rules $\Delta Q^{[n]}=Q^{[n]}-Q^{\infty}$ and
$\Delta P^{[n]}=P^{[n]}-P^{\infty}$. Examination of the error columns in table
1 shows that there are the results $|\Delta Q^{[n+1]}/\Delta Q^{[n]}|\sim 0.2$
and $|\Delta P^{[n+1]}/\Delta P^{[n]}|\sim 0.2$, and therefore the errors
decrease _geometrically_ with each iteration by factors that are not
significantly less than $1$. The number of correct digits increases roughly
_linearly_ with the number of iterations

We would like an iteration process that is more rapidly convergent. Solution
of the equations on the right-hand sides of (74) can be converted into a
fixed-point problem, and this problem can be solved by Newton’s method. The
results of so doing for the problem at hand are shown in table 2. Evidently
with Newton’s method convergence to machine precision has been achieved with
$3$ iterations. And for the errors we find the results $|\Delta
Q^{[n+1]}/(\Delta Q^{[n]})^{2}|\sim 0.2$ and $|\Delta P^{[n+1]}/(\Delta
P^{[n]})^{2}|\sim 0.1$. Therefore the convergence is _quadratic_ as expected
for Newton’s method. The number of correct digits roughly _doubles_ from one
iteration to the next.

##  7 Symplectic Completion Using Cremona Maps

The previous section illustrated how one may begin with a symplectic jet $\cal
J$ truncated at terms beyond degree $(\text{maxm}-1)$ and then, by use of a
suitable Poincaré generating function, add to $\cal J$ terms of degree maxm
and higher to produce an exactly symplectic map. In effect, that approach
usually adds an _infinite_ number of terms beyond degree $(\text{maxm}-1)$,
because the resulting map generally contains singularities. However, as seen
for phase-space regions of physical interest, these additional terms seem to
have little (no deleterious) effect beyond achieving symplectification.

In this section we describe how one may achieve symplectification by adding
only a _finite_ number of terms beyond degree $(\text{maxm}-1)$. The result is
a map that is both _polynomial_ and exactly symplectic. We call such maps
_Cremona_ maps. For simplicity, as was done in the previous section, we
restrict our discussion to the case of one degree of freedom
[TMCA:init]Detailed background material for this section is most easily found
on the Web: D. T. Abell, _Analytic Properties and Cremona Approximation of
Transfer Maps for Hamiltonian Systems_ , PhD dissertation, University of
Maryland, College Park, 1995. %https://www.radiasoft.net/wp-
content/uploads/2021/01/thesis_2e.pdf. In subsequent citations it will be
referred to as _TMCA_. See also the chapter “Symplectic maps and computation
of orbits in particle accelerators” by A. J. Dragt and D. T. Abell, in the
book _Integration Algorithms and Classical Mechanics_ (J. E. Marsden, G. W.
Patrick, and W. F. Shadwick, eds.), vol. 10 of _Fields Inst. Comm._ ,
(Providence, Rhode Island), pp. 59–85, American Mathematical Society, 1996.
See further S. Blanes, “Symplectic maps for approximating polynomial
Hamiltonian systems,” _Phys. Rev. E_ , vol. 65, 056703, May 2002. For a
history of the term _Cremona maps_ , and for a fuller discussion of the
concepts of _kicks_ and _jolts_ than we present here (including the extension
to two and three degrees of freedom), see _TMCA_ , Chapters 10 and 11..

As Cremona maps are symplectic, one can use them to approximate the behavior
of Hamiltonian systems; as they are polynomial, one can compute them rapidly
and exactly. But how might one construct such maps? To answer this question,
we consider first the Lie transformation of a polynomial function of $q$
alone; in other words, a map of the form $\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g(q)}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)$ with $g$ an arbitrary polynomial in $q$.
Using (46) and (30), one may compute the action of such a map to be

$\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g(q)}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\begin{pmatrix}q\\\
p\end{pmatrix}=\begin{pmatrix}q\\\ p+\partial{g}/\partial{q}\end{pmatrix},$
(87)

which is necessarily a polynomial symplectic map. Because this map changes
only the momentum, we call it a _kick map_ and refer to the corresponding Lie
generator $g$ as a _kick_. Now suppose we choose any linear symplectic map
$\mathcal{L}$ and use it to form a more general map $\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{\mathcal{L}g(q)}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)$. With the aid of (47) one may write

$\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{\mathcal{L}g(q)}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)=\mathcal{L}\exp(\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{g(q)}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits)\mathcal{L}^{-1}.$ (88)

Then the linearity of $\mathcal{L}$ assures us that maps of the form (88) are
also necessarily polynomial symplectic maps. We call such a generalized kick
map a _jolt map_ and its Lie generator, $\mathcal{L}g$, a _jolt_. What we have
learned here is that jolt maps—defined by an $\mathcal{L}$ and a $g(q)$—can
supply us with an endless stream of Cremona maps.

Later in this paper, we shall make use of jolt maps to construct a Cremona
symplectification for our example map (54). To do so, we first need to develop
some concepts and tools.

In the vector space of dynamical polynomials (_i. e._ all polynomials on phase
space), we define a set of general basis monomials of degree $l$ by the rule

$G^{(l)}_{r}(z)=\frac{q^{l-r}p^{r}}{\sqrt{(l-r)!r!}}.$ (89)

And for the basis monomial in $q$ alone, we write $Q^{(l)}=G^{(l)}_{0}$; thus,

$Q^{(l)}(q)=\frac{q^{l}}{\sqrt{l!}}.$ (90)

In addition, we introduce a (very special) inner product $\braket{\,,\,}$
defined by the rule

$\braket{G^{(l)}_{r},G^{(l^{\prime})}_{r^{\prime}}}=\delta_{ll^{\prime}}\delta_{rr^{\prime}}.$
(91a) We see that with respect to this inner product, the $G^{(l)}_{r}$
constitute an orthonormal basis for the space of dynamical polynomials. Now
suppose we have dynamical polynomials $f=\sum_{lr}f_{lr}G^{(l)}_{r}$ and
$g=\sum_{lr}g_{lr}G^{(l)}_{r}$. We extend the inner product (91a), to the
entire vector space of dynamical polynomials by defining
$\braket{f,g}=\sum_{lr}f_{lr}\,g_{lr}.$ (91b)

For our purposes, the essential feature of what we call the _invariant scalar
product_ , (91), is that any transformation belonging to the $\mathrm{U}(1)$
subgroup of $\mathrm{Sp}(2,\mathbb{R}^{\mspace{2.0mu}\mspace{-2.0mu}})$ leaves
this inner product unchanged [LM:S.7.3]_LM_ , section 7.3.. That subgroup is
the group of plane rotations, and hence

$\braket{\mathcal{L}f,\mathcal{L}g}=\braket{f,g}$

for any $\mathcal{L}$ having the form, (55), of a plane rotation. Because the
jolt $\mathcal{L}g(q)$ is generally a polynomial in both $q$ and $p$ (recall
(57)), we may hope that a set of jolts $\\{\mathcal{L}_{j}Q^{(l)}(q)\\}$ can
span the space of relevant dynamical polynomials.

We have two other concepts to introduce. The first is that of _sensitivity
vectors_ $\sigma^{r}$, which have components141414The $\sigma^{r}_{j}$ do
depend on $l$, but we have suppressed this index to avoid notational clutter.

$\sigma^{r}_{j}=\braket{G^{(l)}_{r},\mathcal{L}_{j}Q^{(l)}}.$ (92)

These components measure the content of each jolt $\mathcal{L}_{j}Q^{(l)}$
within the given monomial $G^{(l)}_{r}$. The second concept is that of the
_Gram matrix_ $\Gamma(l)$, which has components

$\Gamma(l)_{rs}=\frac{1}{N}\sum_{j=1}^{N}\sigma^{r}_{j}\,\sigma^{s}_{j}\equiv\\{\sigma^{r},\sigma^{s}\\},$
(93)

where we have introduced a weighted scalar product denoted
$\set{\,,\,}$.151515In the case of two or three degrees of freedom, this
becomes $\sum_{j}w_{j}\sigma^{r}_{j}\sigma^{s}_{j}$, with weights $w_{j}$
differing from $1/N$. This symmetric matrix, which measures the uniqueness, or
linear independence, of the different sensitivity vectors, depends only on
one’s choice of $\mathcal{L}_{j}$.

With the above concepts and tools in hand, let us return to the example
nonlinear Lie generator $qp^{2}$ of (56). We ask ourselves, “How can one
express this generator as a linear combination of jolts?” Or the slightly more
general question: How do we determine $N$ jolts $\mathcal{L}_{j}Q^{(3)}$
together with associated _jolt strengths_ $a^{(3)}_{j}$ so as to obtain a
_jolt decomposition_ ,

$f_{3}(q,p)=\frac{1}{N}\sum_{j=1}^{N}a^{(3)}_{j}\mathcal{L}_{j}Q^{(3)},$ (94)

for any homogeneous dynamical polynomial $f_{3}$ of degree three? Because the
monomials $G^{(3)}_{r}$ form a basis for such polynomials, we may, with the
use of (91a), write

$f_{3}=\sum_{r=0}^{3}c^{(3)}_{r}G^{(3)}_{r},\text{ where
}c^{(3)}_{r}=\braket{G^{(3)}_{r},f_{3}}.$ (95)

Now insert the jolt decomposition (94) into the latter equality. We find that

$c^{(3)}_{r}=\frac{1}{N}\sum_{j=1}^{N}a^{(3)}_{j}\braket{G^{(3)}_{r},\mathcal{L}_{j}Q^{(3)}}=\frac{1}{N}\sum_{j=1}^{N}a^{(3)}_{j}\sigma^{r}_{j},=\\{a^{(3)},\sigma^{r}\\}.$
(96)

Since no component of $a^{(3)}$ that lies orthogonal to $\sigma^{r}$ can
contribute to $c^{(3)}_{r}$, we make the _Ansatz_ that the vector of jolt
strengths, $a^{(3)}$, must be a linear combination of the sensitivity vectors;
thus $a^{(3)}=\sum_{s}\alpha^{(3)}_{s}\sigma^{s}$. Inserting this expansion
into (96), we obtain

$c^{(3)}_{r}=\sum_{s=0}^{3}\alpha^{(3)}_{s}\\{\sigma^{s},\sigma^{r}\\}=\sum_{s=0}^{3}\Gamma(3)_{rs}\,\alpha^{(3)}_{s};\text{
hence }c^{(3)}=\Gamma(3)\,\alpha^{(3)}.$ (97)

This result tells us that in order to compute $\alpha^{(3)}$—and hence the
jolt strengths $a^{(3)}$—we require a non-singular Gram matrix $\Gamma(3)$. We
may then compute $\alpha^{(3)}=\Gamma(3)^{-1}c^{(3)}$.

One may parameterize the $\mathcal{L}_{j}$ by rotation angles $\theta_{j}$. It
is then possible to compute the sensitivity vectors analytically, with result

$\sigma^{r}_{j}=\binom{3}{r}^{1/2}c_{j}^{3-r}s_{j}^{r},$ (98)

where $c_{j}=\cos\theta_{j}$ and $s_{j}=\sin\theta_{j}$. One then obtains the
Gram matrix elements in the explicit form [TMCA:S.16.1.1]_TMCA_ , Section
16.1.1

$\Gamma(3)_{rs}={\biggl{[}\binom{3}{r}\binom{3}{s}\biggr{]}}^{1/2}\frac{1}{N}\sum_{j}c_{j}^{6-(r+s)}s_{j}^{r+s}.$
(99)

The questions that remain are (i) what is the best choice of angles
$\theta_{j}$, and (ii) how many do we need? It seems reasonable to be
democratic about our choice of angles, and analysis has indeed shown evenly-
spaced angles to be optimal [TMCA:S.16.1.3]_TMCA_ , Section 16.1.3. A simple
dimension-counting argument tells us that for the case $l=3$, the number of
jolts cannot be less than $4$. However, an analysis of (99)
[TMCA:S.16.1.4ff]_TMCA_ , Sections 16.1.4 and 16.1.6, or direct numerical
computation, shows that four evenly-spaced angles yield a singular
$\Gamma(3)$; but use of five evenly-spaced angles does not.

Let us review our progress so far: We have the third-order Lie generator
$qp^{2}$, which we want to decompose into a linear combination of jolts, as in
(94). Using the evenly-spaced angles $\theta_{j}=2\pi{}j/5$, we construct,
_cf._ (55), the five rotation maps $\mathcal{L}_{j}=\mathcal{R}(\theta_{j}),$
with $j\in\set{0,\dotsc,4}$. We also construct the Gram matrix (99),
sensitivity vectors (98), and the jolt strengths
$a^{(3)}=\sum_{r}\alpha^{(3)}\sigma^{r}$, where
$\alpha^{(3)}=\Gamma(3)^{-1}c^{(3)}$. This allows us, finally, to construct
the jolt decomposition

$qp^{2}\approx-0.6532\,\mathcal{L}_{0}Q^{(3)}+0.8936\,\mathcal{L}_{1}Q^{(3)}-0.5670\,\mathcal{L}_{2}Q^{(3)}-0.5670\,\mathcal{L}_{3}Q^{(3)}+0.8936\,\mathcal{L}_{4}Q^{(3)}.$

The last step is to split this decomposition into five individual jolt maps.
We thereby achieve an approximation to $\mathcal{N}$ in the form of a Cremona
map,

$\mathcal{N}^{\mspace{1.0mu}\text{cr}}=\mathcal{L}_{0}\mathrm{e}^{\mspace{1.5mu}\mspace{-1.5mu}\mathop{\text{\raise
0.2411pt\hbox{{:}}}\mspace{1.0mu}{-0.6532Q^{(3)}}\mspace{1.0mu}\text{\raise
0.2411pt\hbox{{:}}}}\nolimits}\mathcal{L}_{0}^{-1}\,\mathcal{L}_{1}\mathrm{e}^{\mspace{1.5mu}\mspace{-1.5mu}\mathop{\text{\raise
0.2411pt\hbox{{:}}}\mspace{1.0mu}{0.8936Q^{(3)}}\mspace{1.0mu}\text{\raise
0.2411pt\hbox{{:}}}}\nolimits}\mathcal{L}_{1}^{-1}\,\mathcal{L}_{2}\mathrm{e}^{\mspace{1.5mu}\mspace{-1.5mu}\mathop{\text{\raise
0.2411pt\hbox{{:}}}\mspace{1.0mu}{-0.5670Q^{(3)}}\mspace{1.0mu}\text{\raise
0.2411pt\hbox{{:}}}}\nolimits}\mathcal{L}_{2}^{-1}\\\
\times\mathcal{L}_{3}\mathrm{e}^{\mspace{1.5mu}\mspace{-1.5mu}\mathop{\text{\raise
0.2411pt\hbox{{:}}}\mspace{1.0mu}{-0.5670Q^{(3)}}\mspace{1.0mu}\text{\raise
0.2411pt\hbox{{:}}}}\nolimits}\mathcal{L}_{3}^{-1}\,\mathcal{L}_{4}\mathrm{e}^{\mspace{1.5mu}\mspace{-1.5mu}\mathop{\text{\raise
0.2411pt\hbox{{:}}}\mspace{1.0mu}{0.8936Q^{(3)}}\mspace{1.0mu}\text{\raise
0.2411pt\hbox{{:}}}}\nolimits}\mathcal{L}_{4}^{-1}.$ (100)

We make two observation about this form. First, the error made by this
approximation is, according to the Baker-Campbell-Hausdorff theorem
[LM:S.3.7.3]_LM_ , Section 3.7.3, determined by commutators and multiple
commutators of the different jolts, which in this case have degree $4$ and
higher, _i. e._ , starting at one degree higher than the jolts themselves. The
hope is that those higher-degree terms effectively added to $\mathcal{N}$ by
our Cremona factorization do not damage the dynamics. Second, the order of the
factors in (100) is not prescribed. We may reorder them without changing the
degree of approximation.

Figure 4 shows the result of applying
$\mathcal{M}^{\text{cr}}=\mathcal{R}\mathcal{N}^{\mspace{1.0mu}\text{cr}}$
repeatedly to seven initial conditions for the case $\theta/2\pi=0.22$. The
two graphics in that figure correspond to different orderings of the factors
in (100). To facilitate comparison, light gray curves in the background
indicate the corresponding exact results. Evidently, the order of factors can
have a dramatic impact on the result’s absolute accuracy. In the left-hand
graphic, the original outer curve is replaced by five islands, and the next
curve inwards has angular corners. While this effect remains a topic of
research, we point out how one may—absent knowledge of the correct
result—address the associated uncertainty as to which is the more accurate
result.

Figure 4: Phase-space portrait, in the case $\theta/2\pi=0.22$, resulting from
applying the map
$\mathcal{M}^{\text{cr}}=\mathcal{R}\mathcal{N}^{\mspace{1.0mu}\text{cr}}$
repeatedly (2000 times) to the seven initial conditions $(q,p)=(0.01,0)$,
$(0.1,0)$, $(0.15,0)$, $(0.2,0)$, $(0.25,0)$, $(0.3,0)$, and $(0.35,0)$ to
find their orbits. Light gray curves in the background indicate the exact
result. In the left-hand graphic, the jolt maps are applied in numerical
order, _i. e._ $\set{0,1,2,3,4}$. In the right-hand graphic, the jolt maps are
instead applied in the order $\set{3,4,0,1,2}$.

Suppose we rewrite the map $\mathcal{RN}$ in the form

$\mathcal{RN}=\mathcal{R}\mathcal{N}^{1/2}\mathcal{N}^{1/2}.$ (101)

The reduced nonlinear content of the square root, $\mathcal{N}^{1/2}$, will
make a Cremona symplectification of this map more accurate. In words, this
process means we (i) scale the Lie generator $\mathop{\text{\raise
0.34444pt\hbox{{:}}}\mspace{1.0mu}{qp^{2}}\mspace{1.0mu}\text{\raise
0.34444pt\hbox{{:}}}}\nolimits$ by $1/2$; (ii) split $\mathcal{N}^{1/2}$, as
we did $\mathcal{N}$ in (100), to obtain
$\mathcal{N}^{\mspace{1.0mu}\text{rtcr}}=(\mathcal{N}^{\mspace{1.0mu}\text{cr}})^{1/2}$;
and (iii) square this map to obtain an approximation to $\mathcal{M}$, which
we call $\mathcal{M}^{\text{rtcr}}$, in the form

$\mathcal{M}^{\text{rtcr}}=\mathcal{R}\mathcal{N}^{\mspace{1.0mu}\text{rtcr}}\mathcal{N}^{\mspace{1.0mu}\text{rtcr}}.$
(102)

(Note that in our example we now must apply ten successive jolt maps.) This
process, which we call the _root trick_ , can also improve results obtained by
use of generating functions [LM:S.34.3.5]_LM_ , Section 34.3.5. The general
technique of _scaling, splitting, and squaring_ has application to a broad
range of problems in computational physics [LM:S.4.1]_LM_ , Sections 4.1 and
10.8.

The left-hand graphic in Figure 5 shows the result of applying the root trick
to the ordering $\set{3,4,0,1,2}$. Note the improved accuracy as compared to
the right-hand graphic in Figure 4. Not shown is the result of applying the
root trick to the ordering $\set{0,1,2,3,4}$, which exhibits greatly improved
accuracy. In particular, gone are the offensive islands and angular corners
present in the left-hand graphic of Figure 4.

Figure 5: Cremona symplectification using five-plus-five jolts with the root
trick. The left-hand graphic shows the result of applying the root trick to
the map used for the right-hand graphic in Figure 4. And the right-hand
graphic here shows the result of applying the _symmetrized_ root trick to that
same map.

Moreover, we can improve accuracy still further by _symmetrizing_ the factors:
By this we mean that one obtains an improved result by approximating
$\mathcal{N}$ as

$\mathcal{N}^{\mspace{1.0mu}\text{rtscr}}=\mathcal{N}^{\mspace{1.0mu}\text{rtcr}}\mathcal{N}^{\mspace{1.0mu}\text{rtcr}}_{\text{rev}}$
(103)

where the subscript ‘rev’ means that the given map’s factors should be applied
in the _reverse_ order. It can be shown that this symmetrization of the
nonlinear factors leads to an automatic reduction of many spurious higher-
order terms including the _canceling_ of _all_ $f_{4}$ terms. (This is the
same desirable feature found for the use of a Poincaré generating function
described in the previous section.) The full map $\cal{M}$ now has the
approximation

${\cal{M}}^{\rm{rtscr}}=\mathcal{R}{\cal{N}}^{\rm{rtscr}}.$ (104)

The right-hand graphic in Figure 5 shows the result of using
$\mathcal{M}^{\text{rtscr}}$, with the ordering $\set{3,4,0,1,2}$ for
$\mathcal{N}^{\mspace{1.0mu}\text{rtcr}}$, to track particle trajectories.
Evidently, there is now near perfect agreement with exact results, agreement
that is comparable to that found with the use of a Poincaré generating
function.

We have illustrated how to find, over a substantial region of phase space, a
Cremona approximation for a nonlinear map (in one degree of freedom) having
the form ${\cal{N}}=\exp(:f_{3}:)$. In an analogous (but substantially more
complicated) manner, with the use of more jolt maps, one may find, in two and
three degrees of freedom, suitable Cremona approximations to maps of the form
(52). Key to the construction of jolt maps for an $\mathcal{N}$ having
$f_{>2}$ in its factored product form, and acting on phase spaces
corresponding to two or three degrees of freedom, is _an optimal choice of the
linear maps $\mathcal{L}_{j}$_ used in constructing the jolt maps. It is known
that an optimal choice of linear maps is related to the construction of
suitable _cubature_ formulas for various manifolds.161616Cubature formulas are
higher-dimensional analogs of quadrature formulas. How to do so in the case of
one degree of freedom is well understood for all $f_{>2}$; and the case of two
degrees of freedom is reasonably well understood for all $f_{<15}$. Much work
remains for the case of three degrees of freedom. In particular, one would
like to have, for that case, cubature formulas for the manifold $SU(3)/SO(3)$.
Using this approach, a suitable set of 108 $\mathcal{L}_{j}$ has been found
for all $f_{<7}$ [TMCA:PtII]See _TMCA_ , Part II. Also see D. T. Abell, E.
McIntosh, and F. Schmidt, “Fast symplectic map tracking for the CERN Large
Hadron Collider,” _Phys. Rev. ST Accel. Beams_ , vol. 6, 064001, June 2003..

##  8 Concluding Discussion

By design a storage ring has a (one-time-around) _closed_ orbit. Even if the
design is not perfectly executed, there is a fixed-point theorem to the effect
that there is still a closed orbit that is near the design closed orbit.
Consider the passage of particles near the closed orbit through individual
beam-line elements or through collections of successive elements (called
_lumps_) or once around the entire ring. Each such passage is described by a
symplectic map whose jet is computable/known to some order $(maxm-1)$. Each of
these jets can be symplectified to produce symplectic maps that can be used to
propagate particles around the ring by letting them act in succession thereby
producing in effect a _net_ one-turn symplectic map. This operation is called
_tracking_. The slowest, but presumably most accurate, procedure would be to
track element-by-element. A faster procedure would be to track lump-by-lump.
Its accuracy could be checked by comparing its result to element-by-element
tracking results. Even faster and more daring would be to perform full turn-
by-turn tracking using the symplectic map produced by symplectifying the jet
for the full one-turn map. Whatever method is selected, it can be applied
repeatedly a large number of times to simulate the effect of a large number of
turns while being exactly symplectic (to machine precision) and having
accuracy through order $(maxm-1)$. Since evaluation of the action of these
symplectified maps (obtained either by generating function or Cremona
symplectification) on phase space is fairly fast, it is possible to track in
these ways for a relatively large number of turns with relatively modest use
of computer time.

For example, let us consider the case of the LHC, for which a particle passes
through approximately 19,000 elements per turn. One approach that has been
used in the past is to track particles element by element using the
approximation given by (37) and (38). That is, all fringe-field effects are
neglected. Moreover, the maps for drift spaces, bending magnets (dipoles), and
focusing/defocusing magnets (quadrupoles) are approximated by linear (matrix)
maps. Finally, the nonlinear effects of higher multipole magnets (sextupoles,
octupoles, etc.), as well as multipole errors in dipoles and quadrupoles, are
simply treated as kicks. This approach is often referred to as _direct_ /“
_exact_ ” tracking, although what it actually does is equivalent to
implementing a relatively crude, but exactly symplectic, one-turn map.

Since generating function or Cremona tracking begins with symplectic jets, and
these jets can in principle be computed for _realistic_ electromagnetic
fields, generating function or Cremona tracking can in principle be expected
to give more accurate results for realistic machines. However, in order to
assure the Accelerator Physics community of their reliability, generating
function or Cremona tracking should also be able to reproduce the results of
direct tracking. That is, based on the assumptions made for direct tracking,
relatively crude but symplectic jets can be computed for each beam line
element, and these jets can be concatenated to form symplectic jets for lumps
or one-turn maps. These jets can then be symplectified using generating
function or Cremona methods, and their tracking results compared with those
obtained by direct tracking. Preliminary studies/comparisons of this kind for
the LHC, prior to its completed construction, were carried out for various
nonlinear imperfection models. They show that, even for relatively large
betatron amplitudes where nonlinear effects are expected to be important,
there is good agreement in dynamic aperture (phase-space region of long-term
storage) between direct and one-turn Cremona map tracking results using jets
containing generators $f_{<8}$. Moreover, in these studies Cremona tracking is
approximately 20 times faster than direct tracking. And, if one wishes to
simulate orbits with smaller betatron amplitudes, amplitudes associated with
normal LHC operation, then use of generators with $f_{<6}$ appears to be
adequate, in which case Cremona tracking is about 60 times faster than direct
tracking. Finally, the same Cremona map tracking speeds can be achieved to
produce accurate results for realistic machines. Generating function and
Cremona symplectification methods are therefore worth further study,
development, and implementation.

## Acknowledgements

We are grateful to the U. S. Department of Energy Office of Science for
research support over the years on the use of Map and Lie-Algebraic methods in
Accelerator Physics. In addition, we thank RadiaSoft LLC for partial support
provided to one of us (DTA) during the preparation of this paper.

## References

  * [1] Much of the background material for this chapter is most easily found on the Web in a draft book: A. J. Dragt, _Lie Methods for Nonlinear Dynamics with applications to Accelerator Physics_ , http://www.physics.umd.edu/dsat/dsatliemethods.html. It, in turn, provides numerous additional references. In subsequent citations it will be referred to as _LM_. For a discussion of Lagrangians and Hamiltonians for charged-particle motion in electromagnetic fields, see _LM_ , Sections 1.5–1.7.
  * [2] _LM_ , Section 1.6.
  * [3] _LM_ , Section 1.6 and Exercise 1.6.7.
  * [4] _LM_ , footnotes in Exercise 1.6.7.
  * [5] _LM_ , Section 1.3.
  * [6] _LM_ , Subsection 6.4.1.
  * [7] _LM_ , Section 12.9.
  * [8] G. Zhong and J. E. Marsden, “Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators,” Phys. Lett. A, vol. 133, pp. 134–139, Nov. 1988.
  * [9] _LM_ , Section 5.1.
  * [10] _LM_ , Section 6.1.2.
  * [11] _LM_ , Section 8.2.
  * [12] _LM_ , Section 7.6.
  * [13] _LM_ , Chapters 17–25.
  * [14] _LM_ , Subsections 1.2.3, 1.4.1, 1.4.2, and Exercise 1.4.3.
  * [15] _LM_ , Section 10.5 and Chapter 39. Some authors refer to TPSA as _Differential Algebra_ (DA). For an exposition of DA, see M. Berz, _Modern Map Methods in Particle Beam Physics_ , vol. 108 of _Advances in Imaging and Electron Physics_ , Academic Press, San Diego, 1999.
  * [16] _LM_ , Section 6.7.
  * [17] _LM_ , Chapter 34. See also the article B. Erdélyi and M. Berz, “Optimal symplectic approximation of Hamiltonian flows,” _Phys. Rev. Lett._ , vol. 87, 114302, Aug. 2001. Currently we do not find persuasive their invocation of the Hofer metric, but do agree (for other reasons) with their conclusion that use of the Poincaré generating function has several desirable features.
  * [18] _LM_ , Section 34.4.
  * [19] Detailed background material for this section is most easily found on the Web: D. T. Abell, _Analytic Properties and Cremona Approximation of Transfer Maps for Hamiltonian Systems_ , PhD dissertation, University of Maryland, College Park, 1995. https://www.radiasoft.net/wp-content/uploads/2021/01/thesis_2e.pdf. In subsequent citations it will be referred to as _TMCA_. See also the chapter “Symplectic maps and computation of orbits in particle accelerators” by A. J. Dragt and D. T. Abell, in the book _Integration Algorithms and Classical Mechanics_ (J. E. Marsden, G. W. Patrick, and W. F. Shadwick, eds.), vol. 10 of _Fields Inst. Comm._ , (Providence, Rhode Island), pp. 59–85, American Mathematical Society, 1996. See further S. Blanes, “Symplectic maps for approximating polynomial Hamiltonian systems,” _Phys. Rev. E_ , vol. 65, 056703, May 2002. For a history of the term _Cremona maps_ , and for a fuller discussion of the concepts of _kicks_ and _jolts_ than we present here (including the extension to two and three degrees of freedom), see _TMCA_ , Chapters 10 and 11.
  * [20] _LM_ , section 7.3.
  * [21] _TMCA_ , Section 16.1.1.
  * [22] _TMCA_ , Section 16.1.3.
  * [23] _TMCA_ , Sections 16.1.4 and 16.1.6.
  * [24] _LM_ , Section 3.7.3.
  * [25] _LM_ , Section 34.3.5.
  * [26] _LM_ , Sections 4.1 and 10.8.
  * [27] See _TMCA_ , Part II. Also see D. T. Abell, E. McIntosh, and F. Schmidt, “Fast symplectic map tracking for the CERN Large Hadron Collider,” _Phys. Rev. ST Accel. Beams_ , vol. 6, 064001, June 2003.