# Application of high order symplectic integration methods with forward integration steps in beam dynamics

K. Skoufaris, J. Laskar, Y. Papaphilippou, Ch. Skokos

► To cite this version:

K. Skoufaris, J. Laskar, Y. Papaphilippou, Ch. Skokos. Application of high order symplectic integration methods with forward integration steps in beam dynamics. Physical Review Accelerators and Beams, 2022, 25, 10.1103/PhysRevAccelBeams.25.034001. insu-03718987

HAL Id: insu-03718987

https://insu.hal.science/insu-03718987

Submitted on 11 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Application of high order symplectic integration methods with forward integration steps in beam dynamics

K. Skoufaris*

European Organization for Nuclear Research (CERN), CH-1211 Geneva 23, Switzerland;
Department of Physics, University of Crete, P.O. Box 2208, GR-71003 Heraklion, Greece;
and Institute of Theoretical and Computational Physics (ITCP), GR-71003 Heraklion, Greece

J. Laskar

IMCCE, CNRS-UMR8028, Observatoire de Paris, PSL University, Sorbonne University,
77 Avenue Denfert-Rochereau F-75014 Paris

Y. Papaphilippou

European Organization for Nuclear Research (CERN), CH-1211 Geneva 23, Switzerland

Ch. Skokos

Nonlinear Dynamics and Chaos Group, Department of Mathematics and Applied Mathematics,
University of Cape Town, Rondebosch 7701, South Africa

<figure><img src="image_2.png" /></figure>

(Received 6 September 2021; accepted 24 February 2022; published 31 March 2022)

The Hamiltonian describing particle motion in an accelerator belongs to a large class of systems, the members of which can be integrated with a new set of high order symplectic integrators. One benefit of these integrators is their strong numerical stability, which results from the inclusion of only forward integration steps, independent of the order of accuracy. Using these integrators, the transfer map of any multipolar accelerator magnet is derived and presented here. From these maps, the Hamiltonian flow in different lattices is simulated and benchmarked against other well established integration schemes in the accelerator community. By comparing quantities such as the linear phase advance and action invariant, the chromaticity, as well as the working point and the tune spread with amplitude, the superiority of the novel symplectic integrators is demonstrated with respect to accuracy and integration cost.

DOI: 10.1103/PhysRevAccelBeams.25.034001

# I. INTRODUCTION

In nonlinear dynamical systems, the accurate tracking of dynamical variables for long integration periods is paramount for a thorough understanding of their evolution. A characteristic example of such a system is the particle motion in an accelerator like the Large Hadron Collider (LHC) [1] and its High Luminosity upgrade (HL-LHC) [2,3]. In particular for such colliders or storage rings, the particles must be tracked for millions of turns in order to estimate the phase space limits where chaotic particle diffusion occurs, defining the dynamic aperture, which is correlated to the beam lifetime and the overall accelerator performance [4,5]. In order to prevent an artificial

“numerical” diffusion to occur during this long time tracking, the integration of the Hamilton equations of motion should be symplectic, i.e., the Hamiltonian form, the fundamental Poisson bracket, and density of states in phase space should be preserved [6–9].

Each electromagnetic element of a beam line can be modeled by a Hamiltonian $\mathcal{H}$, which is generally nonlinear. As the system described by $\mathcal{H}$ cannot be integrated analytically, numerical integration techniques preserving the symplecticity should be employed. Any symplectic integration scheme is a canonical transformation that propagates the dynamical variables of a new Hamiltonian, slightly different from the original one but integrable, from an initial state $[q(0), p(0)]$ to $[q(\lambda), p(\lambda)]$ for some later “time” $\lambda$. Since this mapping $[q(0), p(0)] \rightarrow [q(\lambda), p(\lambda)]$ describes the exact evolution of the new Hamiltonian, the total energy error is kept bounded.

The considerable interest in long term integration for different dynamical systems guided the development of various symplectic integrators which served the peculiarities of each specific system [10–18]. Within the accelerator community, high order symplectic integrators are most often

<footnote>*kyriakos.skoufaris@gmail.com</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
used for tracking particles through the lattice [10,11,19], while for complicated multiparticle effects such as beam-beam and space charge, specialized symplectic maps were developed [20,21], requiring large computational power to improve their accuracy. Suzuki [22] has shown that any high order integrator of order $\zeta > 2$ should include negative time steps, which appear as coefficients of the symplectic map. A problem arising from these negative steps is the need for large integrator coefficients which can cause the lack of good numerical stability when large time integration steps are employed [13,23]. For an equivalent integration cost, the less accurate symplectic integrators with positive steps (e.g., LEAPFROG) are more efficient and numerically stable than the higher order ones. This stability problem can be mitigated for a class of Hamiltonian systems, such as the ones describing the motion in particle accelerators, with the use of the CSABA<sub>m</sub> & CSBAB<sub>m</sub> symplectic integrators [13], which are high order integration schemes with positive integration steps. The performance of this family of integrators as compared to the TEAPOT<sub>m</sub> [19] and the fourth order Yoshida-Forest-Ruth (YFR<sub>3</sub>) [10,11] symplectic integrators is studied in this paper.

In Sec. II, the Lie algebra needed for the construction of the studied integration schemes is defined and later on the CSABA<sub>m</sub> & CSBAB<sub>m</sub> symplectic integrators are described. In Sec. III and in the Appendixes A–C, the maps needed for the use of the CSABA<sub>m</sub> & CSBAB<sub>m</sub> in accelerators are derived. Finally, the performances of the CSABA<sub>m</sub> & CSBAB<sub>m</sub> symplectic integrators in comparison with the TEAPOT<sub>m</sub> and the YFR<sub>3</sub> ones are presented in Sec. IV, while the conclusions of our work are presented in Sec. V.

## II. CSABA<sub>m</sub> & CSBAB<sub>m</sub> SYMPLECTIC INTEGRATORS

For introducing the CSABA<sub>m</sub> & CSBAB<sub>m</sub> symplectic integrators, the Lie algebraic structure [24] of the Hamiltonian mechanics [25–27] should be employed. For a system with N degrees of freedom, the multiplication rule of this Lie algebra is defined by the Poisson bracket, which is given by

$$ \omega(q, p) \cdot v(q, p) = \{\omega(q, p), v(q, p)\} = L_{\omega}v(q, p) \\ = \sum_{n=1}^{N} \frac{\partial \omega}{\partial p_n} \frac{\partial v}{\partial q_n} - \frac{\partial \omega}{\partial q_n} \frac{\partial v}{\partial p_n}, \quad (1) $$

with $L_{\omega}$, the Lie operator. The canonical conjugate variables can be represented by the generalized coordinate

$q = (q_1, q_2, \dots, q_N)$ and momentum $p = (p_1, p_2, \dots, p_N)$ vectors, while $w(q, p)$ and $v(q, p)$ are general functions of the canonical variables. In this framework, and for a function of the canonical variables $\Psi = (q_1, p_1, q_2, \dots, p_N)$, the Hamilton equations of motion and their formal solution are written as

$$ \frac{d\Psi}{ds} = \{\mathcal{H}, \Psi\} = L_{\mathcal{H}}\Psi, \quad (2a) $$

$$ \begin{aligned} \Psi(s^f) &= \mathcal{M}\Psi(s^i) = \exp\left(\int_{s^i}^{s^f} L_{\mathcal{H}} ds\right) \Psi(s^i) \\ &= \exp(\lambda L_{\mathcal{H}})\Psi(s^i) = \sum_{n \ge 0}^{\infty} \frac{\lambda^n}{n!} L_{\mathcal{H}}^n \Psi(s^i), \quad (2b) \end{aligned} $$

where $s$ is the independent variable and $\lambda = s^f - s^i$ is the integration step from the initial $s^i$ to the final $s^f$ point. Using Eq. (2a) and the left equality of Eq. (2b), the differential equation of the symplectic map $\mathcal{M}$ can be obtained and is given by the equation $\frac{d\mathcal{M}}{ds} = \mathcal{M}L_{\mathcal{H}}$. The solution of this differential equation is the Lie transformation $\mathcal{M} = \exp(\int_{s^i}^{s^f} L_{\mathcal{H}} ds) = \exp(\lambda L_{\mathcal{H}})$ and can be obtained from the use of the Magnus expansion [28] with the condition $\{\mathcal{H}(\Psi; s_1), \mathcal{H}(\Psi; s_2)\} = 0$ for any $s_1$ and $s_2$. Even for nonautonomous systems, this solution $[\exp(\int_{s^i}^{s^f} L_{\mathcal{H}} ds) = \exp(\lambda L_{\mathcal{H}})]$ is true if the time depended Hamiltonian $\mathcal{H}(q, p, s)$ is transformed to an autonomous one $\tilde{\mathcal{H}}(\tilde{q}, \tilde{p})$ with the use of a Legendre transformation that is defined by the generating function $F_2 = \sum_{n=1}^{N} \tilde{p}_n q_n + \tilde{p}_{N+1} s$ (the phase space is expanded to $N+1$ degrees of freedom). The Lie transformation $\exp(\lambda L_{\mathcal{H}})$ propagates the initial conditions $\Psi(s^i)$ to the final ones $\Psi(s^f)$ and the map $\mathcal{M} = \exp(\lambda L_{\mathcal{H}}) = \sum_{n \ge 0}^{\infty} \frac{\lambda^n}{n!} L_{\mathcal{H}}^n$ is symplectic (a canonical transformation) only if the summation defined by the previous series is not truncated.

For many nonlinear systems, especially for the non-integrable ones, the application of the Lie transformation $\exp(\lambda L_{\mathcal{H}})$ on the dynamical variables $\Psi$ is not trivial. However, for Hamiltonians of the form $\mathcal{H} = A + eV$ where $A$ and $V$ do not commute (i.e., $\{A, V\} \ne 0$) and are independently integrable, an alternative but equivalent integration scheme to $\exp(\lambda L_{\mathcal{H}})$ can be found with the use of the Campbell-Baker-Hausdorf theorem [24]. This alternative symplectic scheme is described by Hamiltonian splitting in several maps $\exp(c_n\lambda L_A)$ and $\exp(d_n\lambda L_{eV})$ and can be expressed as

$$ \exp(\lambda L_{\mathcal{K}}) = \prod_{n \ge 1} \exp(c_n \lambda L_A) \exp(d_n \lambda L_{eV}) = \exp[\lambda(L_A + L_{eV})] = \exp(\lambda L_{\mathcal{H}}), \quad (3a) $$

$$ \begin{aligned} \mathcal{K} = & k_{1,1}A + \epsilon k_{1,2}V + \lambda \epsilon k_{2,1}\{A, V\} + \lambda^2 \epsilon k_{3,1}\{A, \{A, V\}\} + \lambda^2 \epsilon^2 k_{3,2}\{\{A, V\}, V\} \\ & + \lambda^3 \epsilon k_{4,1}\{A, \{A, \{A, V\}\}\} + \lambda^3 \epsilon^2 k_{4,2}\{A, \{\{A, V\}, V\}\} + \lambda^3 \epsilon^3 k_{4,3}\{\{\{A, V\}, V\}, V\} + \dots \end{aligned} \quad (3b) $$
where “...” in the expression of the effective Hamiltonian $K$ indicates terms involving higher order commutators. The use of parameter $\epsilon$ just facilitates the track of the power of $V$ and therefore can be set to $\epsilon = 1$, without loss of generality. The coefficients $k_{1,1} = \sum_{n=1} c_n$, $k_{1,2} = \sum_{n=1} d_n$, and the rest of $k_{i>1,j}$ coefficients are polynomials of $c_n$ and $d_n$ with degree $i$. The Lie transformation $\exp(\lambda L_K)$ is equivalent to $\exp(\lambda L_H)$ when the coefficients $c_n$ and $d_n$ satisfy the equations $k_{1,1} = k_{1,2} = 1$ and $k_{i>1,j} = 0$. Although the new integration scheme involved Lie transformations $[\exp(c_n\lambda L_A)$ and $\exp(d_n\lambda L_{eV})]$ that are integrable, the upper limit of the $n$ can be very large ($n \to \infty$) for many Hamiltonian systems. Nevertheless, since any Lie transformation is symplectic [25], a truncated form of Eq. (3a) can be used. In that case, the $c_n$ and $d_n$ coefficients are not enough for eliminating all the $k_{i>1,j}$ polynomials (the upper limit of $n$ is finite, $n \to m$) and therefore, the resulted integrator $\exp(\lambda L_K)$ is an approximated but symplectic alternative of the $\exp(\lambda L_H)$ and Eq. (3a) takes the following form:

$$ \exp(\lambda L_K) = \text{Symmetric} \left[ \prod_{n=1}^{m} \exp(c_n \lambda L_A) \exp(d_n \lambda L_{eV}) \right] \\ = \exp(\lambda L_H) + \mathcal{O}(\lambda^{2\zeta}). \qquad (4) $$

where $\zeta$ is the order of the integration. The accuracy of the $\exp(\lambda L_K)$ can be improved if a symmetric placement of the maps with respect to the center of the full integration step is used $\{\text{Symmetric}[\prod_{n=1}^m \exp(c_n \lambda L_A) \exp(d_n \lambda L_{eV})]\}$. As described in [10,13], for symmetric integrators, the effective

Hamiltonian $K$ of Eq. (3b) includes only even powers of $\lambda$ and thus, the symplectic integrator $\exp(\lambda L_K)$ is given by

$$ \exp(\lambda L_K) = \text{Symmetric} \left[ \prod_{n=1}^{m} \exp(c_n \lambda L_A) \exp(d_n \lambda L_{eV}) \right] \\ = \exp(\lambda L_H) + \mathcal{O}(\lambda^{2\zeta}). \qquad (5) $$

The motion of the particles through each element of an accelerator ring can be described by a Hamiltonian of the type $\mathcal{H} = A + \epsilon V$, where one of the two terms is much weaker than the other and can be seen as a perturbation. In general, for perturbative Hamiltonian systems where $\epsilon V \ll A$, a set of symmetric symplectic integrators with only positive integration steps ($c_n, d_n > 0$) is proposed in [12] and all orders of these integrators are derived in [13]. These integrators, named $SABA_m$ & $SBAB_m$ [with $A \equiv \exp(c_n\lambda L_A)$ and $B \equiv \exp(d_n\lambda L_{eV})$], eliminate the terms $k_{i>1,1}$ where the perturbation ($\epsilon V$) appears at first order and keeps the less important ones $k_{i>1,j>1}$ (where the perturbation $\epsilon V$ appears at higher orders). The resultant order of accuracy for these symplectic integration schemes is $\mathcal{O}(\lambda^{2m}\epsilon + \lambda^2\epsilon^2)$. A more accurate version of these integrators can be used for Hamiltonians where the term $\mathcal{F} = \{\{A, V\}, V\}$ is integrable. The authors in [13] proposed the addition of a correction term at the beginning and the end of the $SABA_m$ & $SBAB_m$ in order to eliminate the $k_{3,2}$ term in Eq. (3b). This correction is described by the Lie transformation $\exp(-\frac{1}{2}f_m\lambda^3\epsilon^2L_F)$ where $f_m$ is the corrector coefficient. The resulted symmetric symplectic integrators after the addition of the corrector, named $CSABA_m$ & $CSBAB_m$, are of order $\mathcal{O}(\lambda^{2m}\epsilon + \lambda^4\epsilon^2)$ and they are defined as follows:

for even $m$:

$$ CSABA_m = \exp\left(-\frac{1}{2}f_m\lambda^3\epsilon^2L_F\right) \exp(c_1\lambda L_A) \exp(d_1\lambda L_{eV}) \cdots \\ \times \exp(d_{\frac{m}{2}}\lambda L_{eV}) \exp(c_{\frac{m}{2}+1}\lambda L_A) \exp(d_{\frac{m}{2}}\lambda L_{eV}) \cdots \\ \times \exp(d_1\lambda L_{eV}) \exp(c_1\lambda L_A) \exp\left(-\frac{1}{2}f_m\lambda^3\epsilon^2L_F\right), \qquad (6a) $$

$$ CSBAB_m = \exp\left(-\frac{1}{2}f_m\lambda^3\epsilon^2L_F\right) \exp(d_1\lambda L_{eV}) \exp(c_2\lambda L_A) \cdots \\ \times \exp(c_{\frac{m}{2}+1}\lambda L_A) \exp(d_{\frac{m}{2}+1}\lambda L_{eV}) \exp(c_{\frac{m}{2}+1}\lambda L_A) \cdots \\ \times \exp(c_2\lambda L_A) \exp(d_1\lambda L_{eV}) \exp\left(-\frac{1}{2}f_m\lambda^3\epsilon^2L_F\right), \qquad (6b) $$

for odd $m$:

$$ CSABA_m = \exp\left(-\frac{1}{2}f_m\lambda^3\epsilon^2L_F\right) \exp(c_1\lambda L_A) \exp(d_1\lambda L_{eV}) \cdots \\ \times \exp(c_{\frac{m}{2}+1}\lambda L_A) \exp(d_{\frac{m}{2}+1}\lambda L_{eV}) \exp(c_{\frac{m}{2}+1}\lambda L_A) \cdots \\ \times \exp(d_1\lambda L_{eV}) \exp(c_1\lambda L_A) \exp\left(-\frac{1}{2}f_m\lambda^3\epsilon^2L_F\right), \qquad (6c) $$
$$
CSBAB_m = \exp\left(-\frac{1}{2}f_m\lambda^3\epsilon^2 L_F\right) \exp(d_1\lambda L_{eV}) \exp(c_2\lambda L_A) \cdots \\
\times \exp(d_{\frac{m+1}{2}}\lambda L_{eV}) \exp(c_{\frac{m+3}{2}}\lambda L_A) \exp(d_{\frac{m+1}{2}}\lambda L_{eV}) \cdots \\
\times \exp(c_2\lambda L_A) \exp(d_1\lambda L_{eV}) \exp\left(-\frac{1}{2}f_m\lambda^3\epsilon^2 L_F\right). \qquad (6d)
$$

The parameters $c_n, d_n$, and $f_m$ can be found in [13]. Since all the $c_n$ and $d_n$ are positive and satisfy the condition $k_{1,1} = k_{1,2} = 1$, their value is getting smaller with the increase of $m$ (and so with the integrators order). Due to that and also because the correctors' "strength" $f_m$ is always smaller than the $c_n$ and $d_n$ (for a given integrator $m$) [13], the stability of the $CSABA_m$ & $CSBAB_m$ integrators is guaranteed even for large integration steps. It must be noted that the $CSABA_m$ & $CSBAB_m$ integrators can be still used even for Hamiltonians where $\epsilon$ is not very small ($\epsilon \to 1$) [29] and the resulted integration scheme is up to fourth order [$O(\lambda^{2m} + \lambda^4)$]. The $CSABA_m$ & $CSBAB_m$

symplectic integrators are successfully used in different dynamical problems [13,23,29,30] but never compared with the standard integration schemes like the ones used by MAD-X [31] and SixTrack [32].

### III. INTEGRATION WITH $CSABA_m$ & $CSBAB_m$ SYMPLECTIC INTEGRATORS FOR AN ACCELERATOR LATTICE

The Hamiltonian that describes the particle motion through the electromagnetic fields (elements) of a storage ring or an accelerator with piecewise constant curvature ($\rho$) is described in [25] and is given by

$$ \mathcal{H}(x, p_x, y, p_y, l, \delta; s) = \frac{\delta}{\beta_{r_0}} - (1 + hx) \left[ \sqrt{\left(\delta + \frac{1}{\beta_{r_0}} - \frac{q\phi}{cP_0}\right)^2 - (p_x - a_x)^2 - (p_y - a_y)^2 - \frac{(1 - \beta_{r_0}^2)}{\beta_{r_0}^2} + a_s} \right], \quad (7) $$

where the path length ($s$) along the beam line is the independent variable. $h = 1/\rho$, $\beta_{r_0}$ is the ratio of the reference particle velocity over the speed of light $c$, $q$ is the particle charge, $P_0$ is the reference particle momentum, and $\delta$ is the energy spread. In this work, the lattices that are used include only magnetic multipoles (i.e., no radio frequency cavities or other elements with electric field are simulated) thus, $\phi = 0$. The magnetic fields are assumed to be zero outside the magnets and vary only along the transverse directions (x, y), therefore, the scaled vector potentials $a_x = a_y = 0$ and $a_s = a_s(x, y)$ ("hard edge" approximation). In addition, ultrarelativistic particles are used ($\beta_{r_0} \to 1$) and since their longitudinal momentum is quite larger than the transverse one, the paraxial approximation $\sqrt{(1+\delta)^2 - p_x^2 - p_y^2} \approx (1+\delta) - \frac{p_x^2+p_y^2}{2(1+\delta)}$ can be applied. Using these simplifications, the general form of the studied Hamiltonian is given by

$$ \mathcal{H}(x, p_x, y, p_y, l, \delta; s) = A(x, p_x, p_y, \delta) + V(x, y, \delta) \quad (8a) $$

with

$$ A(x, p_x, p_y, \delta) = (1 + hx) \frac{p_x^2 + p_y^2}{2(1 + \delta)} \quad (8b) $$

and

$$ V(x, y, \delta) = -(1 + \delta)hx - (1 + hx)a_s, \quad (8c) $$

where $A$ includes the kinematic terms ($p_x, p_y, \delta$) and $V$ depends only on the Cartesian coordinates $(x, y)$ and the energy spread $\delta$.

For the majority of multipoles found in an accelerator, the "potential" $V$ in Eq. (8) can be seen as a perturbation of $A$ thus, the $SABA_m$ & $SBAB_m$ symplectic integrators can be used for the study of the particles motion. Given that $V \ll A$ and using Eqs. (8b) and (8c), the quantity $\mathcal{F} = \{\{A, V\}, V\}$ is written as

$$ \mathcal{F} = \{\{A, V\}, V\} = \frac{1 + hx}{1 + \delta} \left[ \left(\frac{\partial V}{\partial x}\right)^2 + \left(\frac{\partial V}{\partial y}\right)^2 \right] = \frac{1 + hx}{1 + \delta} \left[ \left( (1 + hx) \frac{\partial a_s}{\partial x} + h(1 + \delta + a_s) \right)^2 + \left( (1 + hx) \frac{\partial a_s}{\partial y} \right)^2 \right]. \quad (9) $$
It can be seen in Eq. (9) that $\mathcal{F}$ deepens only on the positions and the energy spread [$\mathcal{F} = \mathcal{F}(x, y, \delta)$] therefore is integrable. The integrability of $\mathcal{F}$ is a direct effect of the paraxial approximation where the kinematical terms in Eq. (8b) are of the form $A = A(q^n, p^2)$ with $n = 0, 1$ and because $V = V(\delta^n, q^w)$ with $w = 0, 1, 2, \dots$. When $\mathcal{F}$ is integrable (the paraxial approximation holds), the more accurate $\text{CSABA}_m$ & $\text{CSBAB}_m$ symplectic integrators can be used, whereas, in the opposite case, the $\text{SABA}_m$ & $\text{SBAB}_m$ symplectic integrators are applicable. In the case of a straight element (lack of dipolar component), the bending radius is infinite ($\rho \to \infty$), thus, the curvature $h$ in Eqs. (8b), (8c), and (9) can be set to zero. For a generic $a_s(x, y)$, the symplectic transfer maps (equations of motion) that are generated by the Lie transformations $\exp(-\frac{1}{2}f_m\lambda^3\epsilon^2L_\mathcal{F})$, $\exp(d_n\lambda L_{eV})$, and $\exp(c_n\lambda L_A)$ are given in Appendixes A and B. Using these results, the exact form of the transfer map for the most common linear and nonlinear magnets (specific $a_s$) is also derived and can be found in Appendix C.

## IV. COMPARISON STUDIES

For the comparison of the $\text{CSABA}_m$ & $\text{CSBAB}_m$ against the $\text{TEAPOT}_m$ [19] and the fourth order $\text{YFR}_3$ [10,11] symplectic integrators, which are commonly used in accelerator tracking codes, quantities that can be calculated analytically and others that require numerical simulations are studied. For the Hamiltonians used in the following studies, the performance of the $\text{CSABA}_m$ & $\text{CSBAB}_m$ symplectic integrators is almost identical (due to their construction similarities) and thus, only the results from the $\text{CSABA}_m$ integrator are presented.

### A. Analytical studies: Transfer matrix accuracy and phase advance calculation

A first and quite fundamental study is to estimate the ability of the different integrators to reproduce the linear

properties of a lattice. For that reason, the approximation $A(x, p_x, p_y, \delta) \approx A_0(x, p_x, p_y, \delta) = \frac{p_x^2+p_y^2}{2(1+\delta)}$ is used and a linear periodic lattice that consists of repeated FODO cells that include dipoles (a periodic alternation of drift sections, dipoles, and quadrupoles) is studied. Since the lattice is linear, the transfer map for any section of this lattice can be expressed in a matrix form ($M$) and can be calculated without the use of a symplectic integrator (exact solution). Using small integration steps (i.e., small $\lambda$), the exact transfer matrix can be expanded in a Taylor series. Comparing each $M_{i,j}$ element of this expanded form with the corresponding one in the transfer matrix generated by the symplectic integrators, the integrators accuracy $\zeta$ can be obtained, as described by Eqs. (3b) and (4) through the term $O(\lambda^\zeta)$. Performing this expansion to the exact transfer matrix of a quadrupole and comparing the elements that describe the horizontal motion ($M_{1,1}$, $M_{1,2}$, $M_{2,1}$, and $M_{2,2}$), the order of accuracy $\zeta$ for different symplectic integrators like the LEAPFROG [33], the $\text{CSABA}_m$ and the $\text{TEAPOT}_m$ is plotted in Fig. 1(a). The numbers in the labels inside the parentheses ($n_k, n_{\text{tot}}$) refer to the amount of kicks $n_k$ and the total maps needed to construct the integrator $n_{\text{tot}}$, e.g., (4, 7) the integrator consists of seven maps from which the $n_k = 4$ are kicks and the others are drifts. The $\text{CSABA}_1$ and all the studied TEAPOT integrators can describe the elements of the quadrupole transfer matrix with the same accuracy. On the other hand, the $\text{CSABA}_{m>1}$ integrators can describe the $M_{i,j}$ components better than any of the $\text{TEAPOT}_m$ except the $M_{2,1}$ component where the two integrator families are equivalent. It seems that the TEAPOT splittings improve slowly with the increase of the number of kicks while the overall behavior of the $\text{CSABAs}$ is two orders of magnitude better than any of the TEAPOTs. The only exception is the $\text{CSABA}_1$ which is equivalent to TEAPOTs behavior.

It is also important to study the integrators capabilities to reproduce correctly quantities that describe the dynamics of

<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) The order of accuracy $\zeta$ for each element $M_{i,j}$ of the quadrupole transfer matrix is plotted for various integrators where the numbers in the labels inside the parentheses ($n_k, n_{\text{tot}}$) refer to the amount of kicks $n_k$ and the total maps needed to construct the integrator $n_{\text{tot}}$. The order of accuracy of $\mu$, in a FODO cell that include dipoles, for different $K_Q$ and $L_Q$ with the use of (b) the $\text{CSABA}_2$ and (c) the $\text{TEAPOT}_5$ symplectic integrators. The area below the white dashed curve guarantees stable motion through the FODO cell.</figcaption></figure>
the problem. Such a quantity in this linear periodic lattice is the phase advance. Since there is no coupling between the planes of motion (no skew quadrupoles), the phase advance $\mu$ after one period (one FODO cell), is equal to $\mu_j = \arccot(\frac{M_{1+w,1+w}}{M_{1+w,2+w}}\beta_j - \alpha_j)$, where M is the block diagonal transfer matrix of the repeated cell including both transverse planes. In this respect, the index j corresponds to x or y and w = 0 for j = x, while w = 2 for j = y. The parameters $\beta_j$ and $\alpha_j$ correspond to the usual Courant-Snyder functions [34] at the beginning and end of the repeated periodic cell. The studied FODO cell include one focusing and one defocusing quadrupole with the same length $L_Q$ and the same magnitude of normalized strength $|K_Q| = |k_1|$. Due to these symmetries, the phase advance $\mu$ in both planes is defined by the ratio $\frac{M_{1,1}}{M_{1,2}}$. Starting from a mirror symmetric point in the cell ($\alpha_x = \alpha_y = 0$) and using different symplectic integrators, the $\mu$ for a set of $K_Q$ and $L_Q$ combinations (different lattices) is calculated. For the CSABA₂ and TEAPOT₅ integrators, their order of accuracy defined as $O(\mu) = |\log_{10}(\frac{|M_{1,1}^{(int)}|\beta - |M_{1,1}^{(ex)}|\beta}{|M_{1,2}^{(int)}|\beta})|$ is plotted in Figs. 1(b) and 1(c), respectively. In this formula, the matrix $M^{(ex)}$ is the exact solution and the $M^{(ex)}$ is the one resulted from the studied symplectic integrator. For the configurations under the white dashed curve where the optics is stable, the CSABA₂ is at least two orders of magnitude more accurate than the TEAPOT₅, even if the TEAPOT₅ consists of more maps as compared with CSABA₂ (eleven in contrast to seven). Indeed, the CSABAₘ & CSBABₘ are not only more accurate than TEAPOTₘ but they are also more economical with respect to the integration time.

Another important quantity that can be calculated analytically in a linear lattice like a FODO cell (a periodic alternation of drifts and quadrupoles) is the action invariant $J = \frac{1}{2\beta_j}[j^2 + (\alpha_j j + \beta_j p_j)^2]$ with j = x, y. Using the exact transfer matrices for the different lattice elements, the exact value of the action ($J^{(ex)}$) is calculated for a variate of lattices (different $K_Q$ and $L_Q$ combinations). Applying the CSABA₂ and TEAPOT₃ integration schemes, the action ($J^{(int)}$) is also calculated for the same lattices. Based on these results, the integrators order of accuracy defined as $O(J) = |\log_{10}(|J^{(int)} - J^{(ex)}|)|$ is plotted in Figs. 2(a) and 2(b) for the CSABA₂ and TEAPOT₃ integrations, respectively. With purple color are plotted the least accurate areas while with red the more accurate ones. Interested in the $K_Q$ and $L_Q$ configuration under the white dashed curve where the motion is stable, the CSABA₂ is either similar or more accurate than the TEAPOT₃ by one order of magnitude. This shows that with the same computational cost (both integrators consist of seven steps), the CSABA₂ can preserve better fundamental quantities of integrable motion, which indeed the primary purpose of a symplectic integrator.

<figure><img src="image_4.png" /><figcaption>FIG. 2. The order of accuracy of J in a FODO cell, for different $K_Q$ and $L_Q$ with the use of (a) the CSABA₂ and (b) the TEAPOT₃ symplectic integrators. The area below the white dashed curve guarantees stable motion through the FODO cell.</figcaption></figure>

## B. Analytical studies: Chromaticity and tune spread with amplitude calculation

In order to study the linear chromaticity from sextupoles $\xi$ and the tune spread with amplitude $\frac{\partial Q_j}{\partial J_j}$ with j = x, y, a periodic nonlinear lattice is used. The repeated cell of this lattice is the same as the one used in the previous linear lattice but with the addition of sextupoles that are used to cancel out the natural chromaticity from quadrupoles. The natural chromaticity ($\xi_n$) and the linear chromaticity from sextupoles ($\xi$) are given by

$$ \xi_{nj} = \frac{\Delta Q_j}{\delta} = \frac{(-1)^{\tau}}{4\pi} \int_{s^i}^{s^f} k_1(s) \beta_j(s) ds, \quad (10a) $$

$$ \xi_j = \frac{\Delta Q_j}{\delta} = \frac{(-1)^{\tau+1}}{8\pi} \int_{s^i}^{s^f} k_2(s) \beta_j(s) D_x(s) ds, \quad (10b) $$

where for j = x, $\tau = 1$ and for j = y, $\tau = 2$. $D_x$ is the horizontal dispersion, $k_1$ and $k_2$ are the normalized
<figure><img src="image_5.png" /><figcaption>FIG. 3. The difference in order of accuracy of $\xi_x$ between (a) CSABA$_2$ & YFR$_3$, (b) CSABA$_2$ & TEAPOT$_3$, and (c) CSABA$_4$ & TEAPOT$_5$. The difference in order of accuracy of $\xi_y$ between (d) CSABA$_2$ & YFR$_3$, (e) CSABA$_2$ & TEAPOT$_3$, and (f) CSABA$_4$ & TEAPOT$_5$. The combinations of ($K_Q, L_Q$) in the white area guarantee unstable motion and the black dots are configurations where the CSABA$_\nu$ is accurate up to the 16th decimal digit of $\xi_x$ or $\xi_y$.</figcaption></figure>

strengths of the quadrupole and the sextupole, respectively. Using the symplectic integrators only for the nonlinear elements of the lattice (same optical functions for all the integrators), the $\xi$ and its order of accuracy $O(\xi) = |\log_{10}(|\xi^{(\text{int})} - \xi^{(\text{ex})}|)|$ are calculated for different $K_Q$ and $L_Q$ combinations (different lattices) that guarantee stability. $O(\xi)$ is the accuracy of the studied symplectic integrator, $\xi^{(\text{ex})}$ is the target value of the sextupole chromaticity, and $\xi^{(\text{int})}$ is the one resulting from employing the studied integrator. In these calculations, any contribution from the sextupoles at the optical functions is neglected since it is a second order effect. For the couples CSABA$_2$ & YFR$_3$, CSABA$_2$ & TEAPOT$_3$, and CSABA$_4$ & TEAPOT$_5$, the difference in their order of accuracy $O(\xi)$ at the horizontal plane is presented in Figs. 3(a), 3(b), and 3(c) and at the vertical one in Figs. 3(d), 3(e), and 3(f). As can be seen in these plots, the results for the two planes are very similar for all the compared integrator couples. Focusing at the results in Figs. 3(a) and 3(d), the seven map integrators CSABA$_2$ and YFR$_3$ present very similar performance across the different ($K_Q, L_Q$) configurations. The average difference from all the studied lattices mean [$O(\xi_j)][CSABA_2] - O(\xi_j)][YFR_3]$] equals 0.25, which demonstrates that the CSABA$_2$ is slightly more accurate. Between the CSABA$_2$

and the other seven step integrator TEAPOT$_3$ [Figs. 3(b) and 3(e)], the CSABA$_2$ is at least one order of magnitude more accurate for small ($K_Q, L_Q$) and improves significantly, as the values of $K_Q$ and $L_Q$ are increased. This behavior is true even for higher order integrators, as can be observed in Figs. 3(c) and 3(f), where the order of accuracy difference between the 11 step integrators CSABA$_4$ and TEAPOT$_5$ is plotted. The accuracy degradation of the TEAPOT$_m$ for larger values of ($K_Q, L_Q$) seems to be an inherent problem of that integration scheme. Since these integrators are constructed to reproduce very well the $M_{2,1}$ element of the quadrupole transfer matrix, they start losing in accuracy for higher order multipole magnets. In contrast, any Hamiltonian-centric integrator, such as the CSABA$_m$ & YFR$_m$, is more stable for larger values of ($K_Q, L_Q$). The black dots in Fig. 3 are points where the CSABA$_m$ is accurate up to the 16th decimal digit of $\xi$, i.e., it reaches machine precision.

For long “time” integration studies that exhibit diffusion [35], the high order nonlinear effects (e.g., tune spread with amplitude from sextupoles, nonlinear energy detuning, etc.) should be also considered for an accurate description of the particle motion. For this reason, the tune spread with amplitude $\frac{\partial Q_j}{\partial J_j}$ with $j = x, y$ generated by the sextupoles in the aforementioned nonlinear lattices is studied.
<figure><img src="image_6.png" /><figcaption>FIG. 4. The difference in order of accuracy of $\frac{\partial Q_x}{\partial J_x}$ between (a) CSABA₂ & YFR₃, (b) CSABA₂ & TEAPOT₃, and (c) CSABA₄ & TEAPOT₅. The difference in order of accuracy of $\frac{\partial Q_y}{\partial J_y}$ between (d) CSABA₂ & YFR₃, (e) CSABA₂ & TEAPOT₃, and (f) CSABA₄ & TEAPOT₅. The combinations of ($K_Q, L_Q$) in the white area guarantee unstable motion.</figcaption></figure>

The formulas for the leading order terms of the tune spread with amplitude generated by a sextupole can be found in [36,37] and as before, the symplectic integrators are used only for the sextupolar magnets. Once again, the order of accuracy $O(\frac{\partial Q_j}{\partial J_j}) = |\log_{10}(|\frac{\partial Q_j^{(\text{int})}}{\partial J_j} - \frac{\partial Q_j^{(\text{ex})}}{\partial J_j}|)|$ is calculated for both transverse planes ($j = x, y$) where $\frac{\partial Q_j^{(\text{ex})}}{\partial J_j}$ is the target value of the sextupoles tune spread with amplitude and $\frac{\partial Q_j^{(\text{int})}}{\partial J_j}$ is the one resulting from the use of the studied integrators. For the comparison between CSABA₂ & YFR₃, CSABA₂ & TEAPOT₃ and CSABA₄ & TEAPOT₅, the difference in their order of accuracy $O(\frac{\partial Q_j}{\partial J_j})$ in the horizontal plane is displayed in Figs. 4(a), 4(b), and 4(c) and for the vertical one in Figs. 4(d), 4(e), and 4(f). Among the seven step integrators (CSABA₂, YFR₃, and TEAPOT₃), CSABA₂ can reach higher accuracy in the tune spread with amplitude estimation, for the majority of the ($K_Q, L_Q$) combinations, in both planes. More specific, in the horizontal plane and for all the $K_Q$ and $L_Q$ values, the CSABA₂ is at least one order more accurate than the YFR₃ as can be seen in Fig. 4(a). Similarly, in the vertical plane [Fig. 4(d)], the CSABA₂ is at least one order of magnitude more accurate than the YFR₃ except from a narrow band of ($K_Q, L_Q$) values, where the YFR₃ is more accurate (blue shades).

The difference between CSABA₂ and TEAPOT₃ is not very significant but again the CSABA₂ is slightly more accurate for all the ($K_Q, L_Q$) configurations in both planes (horizontal and vertical), as it is shown in Figs. 4(b) and 4(e). With the 11 step integrators, the CSABA₄ and TEAPOT₅ present very similar results in the vertical plane [Fig. 4(f)] for almost all the $K_Q$ and $L_Q$ values. In the horizontal one [Fiq. 4(c)], the CSABA₄ is at least one order more accurate than the TEAPOT₅ for the majority of lattices (red shades). For the rest, the two symplectic integrators (CSABA₄ & TEAPOT₅) are equally accurate except a narrow band of configurations where the TEAPOT₅ is more accurate (blue shades). In general, the integrator performances are affected by the lattice configuration (mostly by the quadrupole length and strength, i.e., the phase advance between the integrator kicks), as can be seen in Figs. 5(a) and 5(b). In these plots, the accuracy with which the CSABA₄ [Fig. 5(a)] and the TEAPOT₅ [Figs. 5(b)] can describe the sextupole tune spread with amplitude in the horizontal plane ($\frac{\partial Q_x}{\partial J_x}$) for different lattices (different $K_Q$ and $L_Q$ configurations) is presented. For both integrators, the order of accuracy $O(\frac{\partial Q_j}{\partial J_j})$ is smoothly degrading, as the ($K_Q, L_Q$) values approach the unstable motion area (white color). However,
<figure><img src="image_7.png" /><figcaption>FIG. 5. The variation of the $O(\partial Q_i / \partial J_i)$ for differed lattice configurations using (a) the CSABA$_4$ and (b) the TEAPOT$_5$ symplectic integrators.</figcaption></figure>

this smooth degradation of the accuracy is truncated by narrow bands of lattice configurations, for which the integrator performances are significantly improved. This improvement is mostly due to phase advance between the integrator kicks. These “beneficial” phase advances ($K_Q$ and $L_Q$ pairs) are not the same for the two integrators.

### C. Tracking studies: Tune calculation

The results from particle tracking studies are also of greatest importance for a detailed understanding of the system dynamics. Therefore, fundamental quantities like the invariants of motion (e.g., the tune of particles executing regular motion, etc.) as well as quantities like the tune diffusion (used for distinguishing the regular from the irregular motions) should be described as accurately as possible by the symplectic integration schemes. The value of such a quantity is affected by the numerical precision of the tracking algorithms that disturb the symplecticity and from the integration scheme used. Based on that, the more accurate integrators will be affected less.

In order to further evaluate the differences between the CSABA$_m$ and TEAPOT$_m$ integration schemes, a more complicated lattice, similar to the LHC one during the run II, was constructed and analyzed. More specifically, the linear chromaticity is set to 15 units for both transverse planes, the octupoles strength is maximized in order to enhance the nonlinearities and the radio frequency (rf) cavities are switched off (5D tracking with constant energy spread $\delta$). First, a frequency map analysis is performed for studying the particle dynamics in this lattice. Using the CSABA$_m$ or TEAPOT$_m$ symplectic integrators for the propagation through the nonlinear elements and the highly accurate algorithm NAFF (numerical analysis of fundamental frequencies) [38–40] for the calculation of the horizontal and the vertical tunes ($Q_x$, $Q_y$), the tune diffusion resulted from the use of the CSABA$_2$ and TEAPOT$_3$ integrators is plotted in Fig. 6. The tune diffusion is defined by the formula

$$ \log_{10} \left( \sqrt{(Q_x^{\text{in}} - Q_x^{\text{fi}})^2 + (Q_y^{\text{in}} - Q_y^{\text{fi}})^2} \right) \text{ where } Q_x^{\text{in}} \text{ is the tune at the first 1000 revolutions and } Q_x^{\text{fi}} \text{ the tune at the next 1000 revolutions. As can be seen from the results taken with the CSABA}_2 \text{ [Fig. 6(a) in coordinate space and Fig. 6(b) in tune space] and the ones with the TEAPOT}_3 \text{ [Fig. 6(c) in coordinate space and Fig. 6(d) in tune space], the tune diffusion is described in a very similar way by the two integrators.}$$

Despite the similarities seen in Fig. 6, the integrators effectiveness to capture the Hamiltonian (dynamical) properties of the system is better described from the tune study of particles that trace out Kolmogorov-Arnold-Moser (KAM) tori [41–43] in phase space (nonchaotic particles that move on circular trajectories in normal form phase space) which is an invariant of the motion. Given that the lattice used is nonlinear and thus no exact transfer map exists for the calculation of the target values of the invariant (tune), an “exact” symplectic integration scheme is constructed. This integration scheme uses the exact solutions for the linear elements and for the nonlinear ones, a ray tracing LEAPFROG [44] is applied. The general form of the ray tracing LEAPFROG transfer map ($RTLF_\nu$) is written as

$$ RTLF_m = \prod_{w=1}^{m} \left\{ \exp \left[ (2m)^{-1} \lambda L_A \right] \times \exp \left( m^{-1} \lambda L_{eV} \right) \exp \left[ (2m)^{-1} \lambda L_A \right] \right\}_w. \quad (11) $$

With this symplectic integrator ($RTLF_m$), the new Hamiltonian $K$ is given by the original one $H$ as follows $K = H + O(\sum_{w \ge 1} \frac{\lambda^{2w}}{m^{2w}})$. For large $m$ values, the contribution from the leading order ($\lambda^2$) and from the higher ones is very small since they are scaled with $m$. In our case, the density of the kicks used at each nonlinear element is identical to the one needed for the description up to double precision (16th decimal) the deflection (of the conjugate variables) from a thick quadrupole. More specifically, at least 1000 kicks ($m \ge 1000$) are used for each nonlinear
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 6. Tune diffusion resulted from the use of the *CSABA*₂ symplectic integrator in (a) the coordinate space and (b) in the tune space. Tune diffusion resulted from the use of the *TEAPOT*₃ symplectic integrator in (c) the coordinate space and (d) in the tune space.</figcaption></figure>

element. Using the results of the exact integrator, the tune accuracy $\mathcal{O}(Q_j) = |\log_{10}(|Q_j^{(\text{int})} - Q_j^{(\text{ex})}|)$ in both planes $(j=x, y)$ for the seven steps integrators *CSABA*₂ and *TEAPOT*₃ is estimated. Their accuracy difference in the horizontal and vertical planes is presented in Figs. 7(a) and 7(b), respectively. In the horizontal plane [Fig. 7(a)], the *CSABA*₂ is more accurate than the *TEAPOT*₃ [$\Delta\mathcal{O}(Q_x) > 0$] for 86% of the initial conditions, with the average improvement close to 1 order of magnitude [$\text{Mean}(\Delta Q_x) = 0.8$]. In the vertical plane [Fig. 7(b)], the *CSABA*₂ is more accurate than the *TEAPOT*₃ for almost all of the initial conditions (99.7%) and again the

average improvement is close to 1 order of magnitude [$\text{Mean}(\Delta Q_x) = 0.8$]. The *CSABA*₂ is still better than many *TEAPOT* integrators which consists of more than seven maps. A characteristic example can be seen in Figs. 7(c) and 7(d) where the tune accuracy differences between the *CSABA*₂ and the *TEAPOT*₅ (an integrator that consists of 11 maps) are plotted. In the horizontal plane [Fig. 7(c)], the *CSABA*₂ is more accurate for 60% of the initial conditions and, in the vertical one, for 99.6% of the initial conditions. It is worth noticing that, in the horizontal plane, the areas where the *CSABA*₂ is less accurate than the *TEAPOT*₃ [Fig. 7(a)] or the *TEAPOT*₅ [Fig. 7(c)], the tune diffusion
<figure><img src="image_12.png" /><figcaption>FIG. 7. The difference in order of accuracy of $Q_x$ between (a) CSABA$_2$ & TEAPOT$_3$ and (c) CSABA$_2$ & TEAPOT$_5$. The difference in order of accuracy of $Q_y$ between (b) CSABA$_2$ & TEAPOT$_3$ and (d) CSABA$_2$ & TEAPOT$_5$.</figcaption></figure>

is large, as it is shown in Fig. 6(a). The large tune diffusion is an indication that these particles may not be on KAM tori and thus their tune is not an invariant of motion.

V. SUMMARY AND CONCLUSIONS

In this work, the benefits from employing the CSABA$_m$ & CSBAB$_m$ symplectic integrators over TEAPOT$_m$ and the fourth order YFR$_3$, which are widely used in the accelerator community, are presented. This new set of symplectic integrators (CSABA$_m$ & CSBAB$_m$) is ideal for Hamiltonian systems that include a perturbative term, like the ones found in accelerators, and consists of only forward

integration steps that guarantee the stability of the integrator even with large integration steps.

In order to evaluate the effectiveness of the proposed integrators to describe particle motion in an accelerator, quantities that calculated analytically or through numerical simulations are estimated and compared. For the analytical studies, a plethora of lattices was produced and the compared quantities were the phase advance, the linear action invariant, the linear chromaticity, and the tune spread with amplitude. It was clearly shown that for the vast majority of these lattices, the CSABA$_m$ & CSBAB$_m$ symplectic integrators are more accurate. It is to be stressed that, in most of the cases, the accuracy difference is larger than 1 order of magnitude.
Similarly, exceptionally good results were also obtained from the tracking studies, where the CSABA<sub>m</sub> & CSBAB<sub>m</sub> integrators can estimate the tunes of nonchaotic particles with better accuracy as compared to TEAPOT<sub>m</sub>, with the equiv- alent or even larger number of integration steps. These findings indicate that the implementation of these new symplectic integrators in particle tracking codes will be very beneficial for the beam physics community. In this respect, the transfer maps for any combined function magnet with normal or skew symmetry are presented in the Appendixes.

APPENDIX A: GENERAL FORM OF THE SYMPLECTIC TRANSFER MAPS NEEDED FOR THE USE OF THE CSABA<sub>m</sub> & CSBAB<sub>m</sub> SYMPLECTIC IN ACCELERATORS

Using Eqs. (8b), (8c), and (9) in the Lie transformations for the maps exp (c<sub>n</sub>λL<sub>A</sub>), exp (d<sub>n</sub>λL<sub>V</sub>) and the corrector exp (−½f<sub>m</sub>λ<sup>3</sup>L<sub>F</sub>), the solution of the equations of motion (symplectic transfer maps) for a generic vector potential a<sub>s</sub>(x, y) is given by

$$
e^{c_n \lambda L_A} (x, p_x, y, p_y, l, \delta) : \left\{
\begin{aligned}
x^f &= \frac{1}{h} \{(1 + hx^i)[\cos(\eta) + \frac{p_x^i}{p_y^i} \sin(\eta)]^2 - 1\}, \\
p_x^f &= p_y^i \frac{p_x^i - p_y^i \tan(\eta)}{p_y^i - p_x^i \tan(\eta)}, \\
y^f &= y^i + \frac{1+hx^i}{h} \left[ \frac{p_x^2+p_y^2}{p_y^2} \eta + \frac{p_y^2+p_x^2}{p_x^2} \sin(2\eta) + 2 \frac{p_x^i}{p_y^i} \sin^2(\eta) \right], \\
p_y^f &= p_y^i, \\
l^f &= l^i - \eta \frac{(p_x^2+p_y^2)(1+hx^i)}{h p_y^i (1+\delta^i)}, \\
\delta^f &= \delta^i,
\end{aligned}
\right.
\quad (A1)
$$

$$
e^{d_n \lambda L_V} (x, p_x, y, p_y, l, \delta) : \left\{
\begin{aligned}
x^f &= x^i, \\
p_x^f &= p_x^i + d_n \lambda [h(1 + \delta^i) + ha_s|_i + (1 + hx^i) \frac{\partial a_s}{\partial x}|_i], \\
y^f &= y^i, \\
p_y^f &= p_y^i + d_n \lambda (1 + hx^i) \frac{\partial a_s}{\partial y}|_i, \\
l^f &= l^i - d_n \lambda hx^i, \\
\delta^f &= \delta^i,
\end{aligned}
\right.
\qquad (A2)
$$

$$
e^{-\frac{1}{2}f_m \lambda^3 L_x} (x, p_x, y, p_y, l, \delta) : \left\{
\begin{aligned}
x^f &= x^i, \\
p_x^f &= p_x^i + \frac{f_m \lambda^3}{2(1+\delta^i)} \left[ h \left( (h(1+\delta^i + a_s)_i + (1+hx^i) \frac{\partial a_s}{\partial x})_i \right)^2 + \left( (1+hx^i) \frac{\partial a_s}{\partial y} \right)_i^2 \right) \\
&\quad + 2(h(1+hx^i)) \left( (h(1+\delta^i + a_s)_i + (1+hx^i) \frac{\partial a_s}{\partial x})_i \right) \times \left( 2h \frac{\partial a_s}{\partial x} \right)_i + (1+hx^i) \frac{\partial^2 a_s}{\partial x^2} \right)_i \\
&\quad + h(1+hx^i) \left( \frac{\partial a_s}{\partial y} \right)_i^2 + (1+hx^i)^2 \frac{\partial a_s}{\partial y} \left|_i \frac{\partial^2 a_s}{\partial x \partial y} \right|_i \right], \\
y^f &= y^i, \\
p_y^f &= p_y^i + \frac{f_m \lambda^3 (1+hx^i)}{2(1+\delta^i)} \left[ 2(1+hx^i)^2 \frac{\partial a_s}{\partial y} \right]_i \frac{\partial^2 a_s}{\partial y^2} \left|_i \right. \\
&\quad \left. + 2(h(1+\delta^i + a_s)_i + (1+hx^i) \frac{\partial a_s}{\partial x})_i \right) \left( h \frac{\partial a_s}{\partial y} \right)_i + (1+hx^i) \frac{\partial^2 a_s}{\partial x \partial y} \right)_i, \\
l^f &= l^i - \frac{f_m \lambda^3 (1+hx^i)}{2(1+\delta^i)^2} \left[ 2h(1+\delta^i) \left( h(1+\delta^i + a_s)_i + (1+hx^i) \frac{\partial a_s}{\partial x} \right)_i \right. \\
&\quad \left. - (h(1+\delta^i + a_s)_i + (1+hx^i) \frac{\partial a_s}{\partial x})_i \right)^2 - \left( (1+hx^i) \frac{\partial a_s}{\partial y} \right)_i^2, \\
\delta^f &= \delta^i,
\end{aligned}
\right.
\quad (A3)
$$

where λ = s^f − s^i, h = 1/ρ with ρ the bending radius, η = p_y^i c_n λh / (2(1+δ^i)) and c_n, d_n, and f_m are coefficients of the Lie transformations. The detailed derivation of Eq. (A1) can be found in Appendix B, while the derivation of Eqs. (A2) and (A3) is not presented here since are trivial.
APPENDIX B: DERIVATION OF THE
SYMPLECTIC TRANSFER MAP THAT
DESCRIBES THE EFFECT OF THE LIE
TRANSFORMATION exp (c_n λ L_A) ON THE
DYNAMICAL VARIABLES

For $A = (1 + hx) \frac{p_x^2+p_y^2}{2(1+\delta)}$ and $h$ being piecewise constant, the action of $\exp(c_n\lambda L_A)$ on the dynamical variables $(x, p_x, y, p_y, l, \delta)$ is equivalent to solve the following Hamilton equations:

$$
\frac{dx}{ds} = \frac{\partial A}{\partial p_x} = (1 + hx) \frac{p_x}{(1 + \delta)} \quad (\text{B1a})
$$

$$
\frac{dp_x}{ds} = -\frac{\partial A}{\partial x} = -h \frac{p_x^2 + p_y^2}{2(1 + \delta)} \quad (\text{B1b})
$$

$$
\frac{dy}{ds} = \frac{\partial A}{\partial p_y} = (1 + hx) \frac{p_y}{(1 + \delta)} \quad (\text{B1c})
$$

$$
\frac{dp_y}{ds} = -\frac{\partial A}{\partial y} = 0 \Rightarrow \boxed{p_y^f = p_y^i} \qquad (\text{B1d})
$$

$$
\frac{dl}{ds} = \frac{\partial A}{\partial \delta} = -(1 + hx) \frac{p_x^2 + p_y^2}{2(1 + \delta)^2} \quad (\text{B1e})
$$

$$
\frac{d\delta}{ds} = -\frac{\partial A}{\partial l} = 0 \Rightarrow \boxed{\delta^f = \delta^i}. \qquad (\text{B1f})
$$

Inserting Eqs. (B1d) and (B1f) in Eq. (B1b), the evolution of $p_x$ is given by

$$
\begin{align}
\frac{dp_x}{ds} &= -h \frac{p_x^2 + p_y^2}{2(1 + \delta^i)} \Rightarrow \int_{p_x^i}^{p_y^f} \frac{dp_x}{p_x^2 + p_y^2} = - \int_{c_n s^i}^{c_n s^f} \frac{h}{2(1 + \delta^i)} ds \nonumber \\
&\Rightarrow \arctan\left(\frac{p_x^f}{p_y^f}\right) = \arctan\left(\frac{p_x^i}{p_y^i}\right) - c_n \lambda \frac{hp_y^i}{2(1 + \delta^i)} \Rightarrow p_x^f = p_y^i \tan(\sigma - \eta), \tag{B2}
\end{align}
$$

where $\lambda = s^f - s^i$, $\sigma = \arctan(p_x^i)$ and $\eta = p_y^i \frac{c_n \lambda h}{2(1+\delta^i)}$. By applying the identity $\tan(\kappa + \omega) = \frac{\tan(\kappa) - \tan(\omega)}{1 + \tan(\kappa) \tan(\omega)}$ with $\kappa = \sigma$ and $\omega = \eta$, Eq. (B2) is written as

$$
p_x^f = p_y^i \frac{p_x^i - p_y^i \tan(\eta)}{p_y^i - p_x^i \tan(\eta)}. \quad (B3)
$$

Implementing Eqs. (B2) and (B1f) in Eq. (B1a) results in

$$
\begin{align}
\frac{dx}{(1+hx)} &= \frac{p_y^i \tan(\sigma - \eta)}{(1+\delta^i)} ds \Rightarrow \int_{x^i}^{x^f} \frac{dx}{1+hx} = \frac{2}{h} \int_0^\eta \tan(\sigma - \tilde{\eta}) d\tilde{\eta} \nonumber \\
&\Rightarrow \operatorname{Log}_e \left[ \frac{1+hx^f}{1+hx^i} \right] = \operatorname{Log}_e \left[ \frac{\cos(\sigma-\eta)}{\cos(\sigma)} \right] \Rightarrow \boxed{x^f = \frac{1}{h} \left[ (1+hx^i) \left( \cos(\eta) + \frac{p_y^i}{p_x^i} \sin(\eta) \right)^2 - 1 \right].} \tag{B4}
\end{align}
$$

For the evolution of $y$, Eqs. (B4) and (B1d) are used in Eq. (B1c) and the result is expressed as

$$
\begin{align}
\frac{dy}{ds} &= (1 + hx^i) \left( \cos(\eta) + \frac{p_x^i}{p_y^i} \sin(\eta) \right)^2 \frac{p_y^i}{(1 + \delta^i)} \Rightarrow \int_{y^i}^{y^f} dy = \int_0^\eta \frac{2(1 + hx^i)}{h} \left( \cos(\tilde{\eta}) + \frac{p_x^i}{p_y^i} \sin(\tilde{\eta}) \right)^2 d\tilde{\eta} \nonumber \\
&\Rightarrow y^f = y^i + \frac{1 + hx^i}{h} \left[ \frac{p_x^i + p_y^i}{p_y^i} \eta + \frac{p_y^i + p_x^i}{p_y^i} \sin(2\eta) + 2 \frac{p_x^i}{p_y^i} \sin^2(\eta) \right]. \tag{B5}
\end{align}
$$

Finally, for the evolution of $l$, Eqs. (B4), (B2), (B1d), and (B1f) are used in Eq. (B1e) and after the integration over $\eta$, the result is given by

$$
\int_{l^i}^{l^f} dl = \int_{0}^{\eta} -(1 + hx^i) \left( \cos(\tilde{\eta}) + \frac{p_x^i}{p_y^i} \sin(\tilde{\eta}) \right)^2 \frac{p_y^i [1 + \tan^2(\sigma - \tilde{\eta})]}{h(1 + \delta^i)} d\tilde{\eta} \Rightarrow \boxed{l^f = l^i - \eta \frac{(p_x^i + p_y^i)(1 + hx^i)}{h p_y^i (1 + \delta^i)}}. \quad (B6)
$$
## APPENDIX C: SYMPLECTIC TRANSFER MAPS FOR THE LIE TRANSFORMATIONS $\exp(d_n\lambda L_V)$ AND $\exp(-\frac{1}{2}f_m\lambda^3 L_F)$ IN THE CASE OF LINEAR AND NONLINEAR COMBINED FUNCTION MAGNETS

For different combined function multipoles like the linear combined function dipole-quadrupole and the non-linear combined function sextupole-octupole, the action of $\exp(d_n\lambda L_V)$ and $\exp(-\frac{1}{2}f_m\lambda^3 L_F)$ on the dynamical variables $(x, p_x, y, p_y, l, \delta)$ is presented here. From these results, the transfer maps for pure dipole, quadrupole, sextupole, octupole, etc., can be also obtained. The general formulas of V and $\mathcal{F}$ for an arbitrary $a_s$ are given in Eqs. (8c) and (9), respectively. For straight combined function magnets (i.e., except from dipoles where $h^{-1} \neq 0$), the rescaled vector potential $a_s$, in curvilinear coordinates, for normal multipoles is given by

$$a_s = \Re \left( -\sum_{w>1}^{G} \frac{q}{P_0} B_0 \frac{b_{w-1}}{w} r^w e^{iw\theta} \right) = \Re \left( -\sum_{w>1}^{G} \frac{k_{w-1}}{w!} r^w e^{iw\theta} \right) \quad (C1)$$

and for skew ones, it is given by

$$a_s = \Re \left( -\sum_{w>1}^{G} \frac{q}{P_0} B_0 \frac{i\alpha_{w-1}}{w} r^w e^{iw\theta} \right) = \Re \left( -\sum_{w>1}^{G} \frac{i\bar{k}_{w-1}}{w!} r^w e^{iw\theta} \right). \quad (C2)$$

In Eqs. (C1) and (C2), the normalized strength of the normal and skew multipoles is given by $k_w = \frac{g}{P_0} w! B_0 b_w$ and $\bar{k}_w = \frac{g}{P_0} w! B_0 a_w$, respectively, $w = 2, 3, ...$ for the 2wth multipoles, $r^w e^{iw\theta} = (x + iy)^w$, $B_0$ is the magnetic field of the main dipole, $P_0$ the reference particle momentum, and the symbol $\Re$ denotes the real part. Using Eqs. (C1) and (C2), the general form of the function $\mathcal{F}$ for straight combined function magnets is written as

$$\mathcal{F}(x, y, \delta) = \frac{1}{1 + \delta} \sum_{w>1}^{G} \left[ g_{w-1}^2 r^{2(w-1)} + \sum_{u=w+1}^{G} 2g_{w-1}g_{u-1} r^{w+u-2} T_{u-w} \left[ \frac{x}{r} \right] \right], \quad (C3)$$

where $g_w = (\frac{k_w}{w!}, \frac{\bar{k}_w}{w!})$, $G \ge w$, and $T_j[z]$ are the Chebyshev polynomials of the first kind. For a single multipole ($G = w > 1$), Eq. (C3) is reduced to $\mathcal{F} = \frac{1}{1+\delta} g_{w-1}^2 r^{2(w-1)}$. Equation (C3) has the same form for normal and skew multipoles of the same order but with different coefficients $g_w$. It must also be noted that the function $\mathcal{F}$ is only a mathematical object that improves the accuracy of the SABA_m & SBAB_m and not a Maxwellian field despite their similarities. Using Eqs. (A2), (A3), (C1), (C2), and (C3) for straight ($h = 0$) combined function multipoles, the symplectic transfer maps that describe the effect of the Lie transformations $\exp(d_n\lambda L_V)$ and $\exp(-\frac{1}{2}f_m\lambda^3 L_F)$ on the dynamical variables are given by

$$e^{d_n \lambda L_V}(x, p_x, y, p_y, l, \delta): \left\{ \begin{array}{l} x^f = x^i, \\ p_x^f = p_x^i + d_n \lambda \Re(-\sum_{w>1}^G \tilde{g}_{w-1} (x + iy)^{w-1}), \\ y^f = y^i, \\ p_y^f = p_y^i + d_n \lambda \Re(-\sum_{w>1}^G i \tilde{g}_{w-1} (x + iy)^{w-1}), \\ l^f = l^i, \\ \delta^f = \delta^i, \end{array} \right. \quad (C4)$$
$$
e^{-\frac{1}{2}f_m\lambda^3 L_F}(x, p_x, y, p_y, l, \delta): \left\{
\begin{aligned}
& x^f = x^i, \\
& p_x^f = p_x^i + \frac{f_m\lambda^3}{2(1+\delta^i)} \left[ \sum_{w>1}^G 2g_{w-1}^2 x^i(w-1)(x^{i2} + y^{i2})^{w-2} \right. \\
& \qquad \left. + \sum_{u=w+1}^G 2g_{u-1}g_{w-1} \left( (u+w-2)x^i(x^{i2} + y^{i2})^{\frac{u+w-4}{2}} T_{u-w} \left[ \frac{x^i}{\sqrt{x^{i2} + y^{i2}}} \right] \right. \right. \\
& \qquad \left. \left. + (u-w)(x^{i2} + y^{i2})^{\frac{u+w-3}{2}} \left(1 - \frac{x^{i2}}{x^{i2}+y^{i2}}\right) U_{u-w-1} \left[ \frac{x^i}{\sqrt{x^{i2} + y^{i2}}} \right] \right) \right], \\
& y^f = y^i, \\
& p_y^f = p_y^i + \frac{f_m\lambda^3}{2(1+\delta^i)} \left[ \sum_{w>1}^G 2g_{w-1}^2 y^i(w-1)(x^{i2} + y^{i2})^{w-2} \right. \\
& \qquad \left. + \sum_{u=w+1}^G 2g_{u-1}g_{w-1} \left( (u+w-2)y^i(x^{i2} + y^{i2})^{\frac{u+w-4}{2}} T_{u-w} \left[ \frac{x^i}{\sqrt{x^{i2} + y^{i2}}} \right] \right. \right. \\
& \qquad \left. \left. - (u-w)x^i y^i(x^{i2} + y^{i2})^{\frac{u+w-5}{2}} U_{u-w-1} \left[ \frac{x^i}{\sqrt{x^{i2} + y^{i2}}} \right] \right) \right], \\
& l^f = l^i + \frac{f_m\lambda^3}{2(1+\delta^i)^2} \left[ \sum_{w>1}^G g_{w-1}^2(x^{i2} + y^{i2})^{w-1} \right. \\
& \qquad \left. + \sum_{u=w+1}^G 2g_{u-1}g_{w-1}(x^{i2} + y^{i2})^{\frac{u+w-2}{2}} T_{u-w} \left[ \frac{x^i}{\sqrt{x^{i2} + y^{i2}}} \right] \right], \\
& \delta^f = \delta^i
\end{aligned}
\right. \quad (C5)
$$

where $\tilde{g}_w = (\frac{k_w}{w!}, \frac{i\tilde{k}_w}{w!})$ and $U_j[z]$ are the Chebyshev polynomials of the second kind.

1. Combined function dipole-quadrupole

For a combined function dipole-quadrupole magnet with constant inverse bending radius $h = \rho^{-1} = \text{const}$, the function $V$, in curvilinear coordinates, is equal to

$$
V = -(1 + \delta)hx - (1 + hx)a_s = -(1 + \delta)hx + (1 + hx) \left[ hx - \frac{h^2 x^2}{2(1 + hx)} + \frac{k_1}{2} \left( x^2 - y^2(1 + hx) \frac{\sinh^{-1}\left(\frac{hy}{1+hx}\right)}{hy} \right) \right. \\
\left. - \frac{hk_1 x^3}{6(1 + hx)} - \frac{k_1 y^2}{6} \left( \frac{2(1 + hx)^2}{h^2 x^2} + \left(1 - \frac{2(1 + hx)^2}{h^2 x^2}\right) \sqrt{1 + \frac{h^2 x^2}{(1 + hx)^2}} \right) \right], \quad (C6)
$$

where the rescaled vector potential $a_s$, in curvilinear coordinates, can be found in [25]. In order to describe correctly the linear dynamics and the transfer maps of a pure dipole ($h \neq 0, k_1 = 0$) and a pure quadrupole ($h = 0, k_1 \neq 0$), the following approximated forms for $V$ and $F$ are used:
$$V = -hx\delta + \frac{h^2 + k_1}{2}x^2 - \frac{k_1}{2}y^2, \quad (C7a)$$

$$\mathcal{F} = \frac{1 + hx}{1 + \delta} \{[x(h^2 + k_1) - h\delta]^2 + (k_1 y)^2\}. \quad (C7b)$$

The higher-order terms in Eq. (C6) are neglected since they are important mostly for the particles far from the reference trajectory. From Eqs. (C6) and (C7a), the corresponding transfer maps are given by

$$e^{d_n \lambda L_V}(x, p_x, y, p_y, l, \delta): \left\{ \begin{array}{l} x^f = x^i, \\ p_x^f = p_x^i + d_n \lambda [h \delta^i - (h^2 + k_1) x^i], \\ y^f = y^i, \\ p_y^f = p_y^i + d_n \lambda k_1 y^i, \\ l^f = l^i - d_n \lambda h x^i, \\ \delta^f = \delta^i, \end{array} \right. \quad (C8a)$$

$$e^{-\frac{1}{2}f_m \lambda^3 L_F}(x, p_x, y, p_y, l, \delta): \left\{ \begin{array}{l} x^f = x^i, \\ p_x^f = p_x^i + \frac{f_m \lambda^3 h}{2(1+\delta^i)} \{[x^i(h^2+k_1)-h\delta^i]^2 + (k_1 y^i)^2\} + \frac{f_m \lambda^3 (1+h x^i)}{1+\delta^i} (h^2+k_1)[x^i(h^2+k_1)-h\delta^i], \\ y^f = y^i, \\ p_y^f = p_y^i + \frac{f_m \lambda^3 (1+h x^i)}{1+\delta^i} k_1^2 y^i, \\ l^f = l^i + \frac{f_m \lambda^3 (1+h x^i)}{2(1+\delta^i)^2} \{x^i(h^2+k_1)[x^i(h^2+k_1)+2h] - h^2 \delta^i - 2h^2 \delta^i + k_1^2 y^i\}, \\ \delta^f = \delta^i. \end{array} \right. \quad (C8b)$$

In the case of a bending magnet ($h \neq 0$), apart from Eqs. (A1), (C8a), and (C8) which describe the motion of the particle inside the magnet, the effect of the nonperpendicular entrance and exit of the particles in the magnetic field should be taken into account. This effect is described in [45] and can be well approximated by a thin kick that is expressed as

The angle $\theta$ is defined by the face of the magnet and the perpendicular to the beam line plane. When the particle enters the magnet $\theta = \theta_1$ and at the exit $\theta = \theta_2$.

## 2. Combined function sextupole-octupole

For a combined function sextupole-octupole magnet, the expressions for $V$ and $\mathcal{F}$ are written as

$$T(\theta): \left\{ \begin{array}{l} x^f = x^i, \\ p_x^f = p_x^i + hx^i \tan(\theta), \\ y^f = y^i, \\ p_y^f = p_y^i + hy^i \tan(\theta), \\ l^f = l^i, \\ \delta^f = \delta^i. \end{array} \right. \quad (C9)$$

$$V = a_s = \frac{k_2}{6} (x^3 - 3xy^2) + \frac{k_3}{24} (x^4 - 6x^2y^2 + y^4) \quad (C10a)$$

$$\mathcal{F} = \frac{(x^2 + y^2)^2 [(3k_2 + k_3x)^2 + k_3^2 y^2]}{36(1 + \delta)} \quad (C10b)$$

and the corresponding transfer maps are given by

$$e^{d_n \lambda L_V}(x, p_x, y, p_y, l, \delta): \left\{ \begin{array}{l} x^f = x^i, \\ p_x^f = p_x^i - d_n \lambda \left[ \frac{k_2}{2} (x^{i2} - y^{i2}) + \frac{k_3}{6} (x^{i3} - 3x^i y^{i2}) \right], \\ y^f = y^i, \\ p_y^f = p_y^i + d_n \lambda \left[ k_2 x^i y^i + \frac{k_3}{6} (3y^i x^{i2} - y^{i3}) \right], \\ l^f = l^i, \\ \delta^f = \delta^i, \end{array} \right. \quad (C11a)$$
$$e^{-\frac{1}{2}f_m \lambda^3 L_x}(x, p_x, y, p_y, l, \delta): \left\{ \begin{array}{l} x^f = x^i, \\ p_x^f = p_x^i + \frac{f_m \lambda^3}{12(1+\delta^i)} [(x^{i2} + y^{i2})(x^i(2k_2 + k_3x^i)(3k_2 + k_3x^i) + k_3y^{i2}(k_2 + k_3x^i))], \\ y^f = y^i, \\ p_y^f = p_y^i + \frac{f_m \lambda^3}{12(1+\delta^i)} [y^i(x^{i2} + y^{i2})(6k_2^2 + k_3^2(x^{i2} + y^{i2}) + 4k_2k_3x^i)], \\ l^f = l^i + \frac{f_m \lambda^3}{72(1+\delta^i)^2} [(x^{i2} + y^{i2})^2((3k_2 + k_3x^i)^2 + k_3^2y^{i2})], \\ \delta^f = \delta^i. \end{array} \right. \qquad (C11b)$$

The transfer maps of a pure sextupole and a pure octupole are resulted from Eqs. (C10) and (C11) by setting $k_2 = 0$ or $k_3 = 0$, respectively.

[1] O. S. Brüning, P. Collier, P. Lebrun, S. Myers, R. Ostojic, J. Poole, and P. Proudlock, LHC design report, Report No. CERN-2004-003-V-1, CERN Yellow Reports: Monographs, CERN, Geneva, 2004.

[2] G. Apollinari, I. Béjar Alonso, O. Brüning, P. Fessia, M. Lamont, L. Rossi, and L. Tavian, High-Luminosity Large Hadron Collider (HL-LHC): Technical Design Report V. 0.1, Report No. CERN-2020-010, CERN Yellow Reports: Monographs, CERN, Geneva, 2017.

[3] L. Rossi and O. Brüning, Progress with the high luminosity LHC project at CERN, in Proceedings of the 10th International Particle Accelerator Conference, Melbourne, Australia, IPAC-2019 (JACoW, Geneva, 2019), MOYPLM3.

[4] D. Pellegrini, F. Antoniou, G. Arduini, S. Fartoukh, G. Iadarola, N. Karastathis, S. Papadopoulou, Y. Papaphilippou, and G. Sterbini, Incoherent beam-beam effects and lifetime optimisation, in Proceedings of the 2017 Workshop on LHC Beam Operations, Evian Les Bains, France (CERN, Evian Les Bains, France, 2019).

[5] S. Kostoglou, F. Antoniou, I. Efthymiopoulos, G. Iadarola, N. Karastathis, S. Papadopoulou, Y. Papaphilippou, K. Paraschou, D. Pellegrini, G. Sterbini, and G. Trad, Luminosity and lifetime modeling and optimization, in Proceedings of the 9th LHC Operations Evian Workshop, Evian Les Bains, France (CERN, Evian Les Bains, France, 2019).

[6] H. Goldstein, C. Poole, and J. Safko, Classical Mechanics, 3rd ed. (Addison-Wesley, San Francisco, CA, 2002).

[7] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd ed. (Springer, Dordrecht, 2006).

[8] K. Feng and M. Qin, Symplectic Geometric Algorithms for Hamiltonian Systems (Springer, Dordrecht, 2010).

[9] S. Blanes and F. Casas, A Concise Introduction to Geometric Numerical Integration (Routledge & CRC Press, Boca Raton, 2016), 10.1201/b21563.

[10] Haruo Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A 150, 262 (1990).

[11] E. Forest and R. D. Ruth, Fourth-order symplectic integration, Physica (Amsterdam) 43D, 105 (1990).

[12] R. I. McLachlan, Composition methods in the presence of small parameters, BIT 35, 258 (1995).

[13] J. Laskar and P. Robutel, High order symplectic integrators for perturbed hamiltonian systems, Celest. Mech. Dyn. Astron. 80, 39 (2001).

[14] A. C. Petit, J. Laskar, G. Boué, and M. Gastineau, High-order regularised symplectic integrator for collisional planetary systems, Astron. Astrophys. 628, A32 (2019).

[15] S. A. Chin and C. R. Chen, Forward symplectic integrators for solving gravitational few-body problems, Celest. Mech. Dyn. Astron. 91, 301 (2005).

[16] C. Danieli, B. Many Manda, T. Mithun, and Ch. Skokos, Computational efficiency of numerical integration methods for the tangent dynamics of many-body Hamiltonian systems in one and two spatial dimensions, Math. Biosci. Eng. 1, 447 (2019).

[17] B. Senyange and Ch. Skokos, Computational efficiency of symplectic integration schemes: Application to multidimensional disordered Klein-Gordon lattices, Eur. Phys. J. Spec. Top. 227, 625 (2018).

[18] M. Titze, J. Bahrdt, and G. Wüstefeld, Symplectic tracking through straight three dimensional fields by a method of generating functions, Phys. Rev. Accel. Beams 19, 014001 (2016).

[19] H. Burkhardt, R. D. Maria, M. Giovannozzi, and T. Risselada, Improved teapot method and tracking with thick quadrupoles for the LHC and its upgrade, in Proceedings of the 4th International Particle Accelerator Conference, IPAC'13 (JACoW Publishing, Shanghai, China, 2013), pp. 945-947.

[20] K. Hirata, H. W. Moshammer, and F. Ruggiero, A symplectic beam-beam interaction with energy change, Part. Accel. 40, 205 (1992).

[21] J. Qiang, Symplectic multiparticle tracking model for self-consistent space-charge simulation, Phys. Rev. Accel. Beams 20, 014203 (2017).

[22] M. Suzuki, General theory of fractal path integrals with applications to many-body theories and statistical physics, J. Math. Phys. (Cambridge, Mass.) 32, 400 (1991).

[23] L. Nadolski and J. Laskar, Application of a new class of symplectic integrators to accelerator tracking, in Proceedings of the 8th European Particle Accelerator Conference, EPAC 2002, Paris, France, 2002, pp. 1276-1278, http://accelconf.web.cern.ch/AccelConf/e02/PAPERS/WEPLE076.pdf.
[24] B. C. Hall, *Lie Groups, Lie Algebras, and Representations: An Elementary Introduction* (Springer, Berlin, 2015).

[25] A. Wolski, *Beam Dynamics in High Energy Particle Accelerators* (Imperial College Press, London, 2014).

[26] E. Forest, *Beam Dynamics: A New Attitude and Framework* (Harwood, Amsterdam, 1998).

[27] K. Skoufaris, Non-linear dynamics modeling in accelerators with the use of symplectic integrators, Ph.D. thesis, University of Crete, 2021.

[28] W. Magnus, On the exponential solution of differential equations for a linear operator, Commun. Pure Appl. Math. **7**, 649 (1954).

[29] Ch. Skokos, D. O. Krimer, S. Komineas, and S. Flach, Delocalization of wave packets in disordered nonlinear chains, Phys. Rev. E **79**, 056211 (2009).

[30] Y. Papaphilippou and Ch. Skokos, Nonlinear dynamics study of the CLIC damping rings using symplectic integrators, Conf. Proc. C0806233, MOPP061 (2008).

[31] D. Laurent, G. Hans, R. Ghislain, and F. Schmidt, *The MAD-X Program (Methodical Accelerator Design) Version 5.05.02 User's Reference Manual* (CERN, Geneva, Switzerland, 2019), https://mad.web.cern.ch/mad/.

[32] R. De Maria et al., SixTrack V and runtime environment, Int. J. Mod. Phys. A **34**, 1942035 (2019).

[33] R. Ruth, A canonical integration technique, *IEEE Trans. Nucl. Sci.* **30**, 2669 (1983).

[34] E. Courant and H. Snyder, Theory of the alternating-gradient synchrotron, Ann. Phys. (Paris) **3**, 1 (1958).

[35] Instability of dynamical systems with several degrees of freedom, in *Collected Works: Representations of Functions, Celestial Mechanics and KAM Theory, 1957–1965*, edited by A. B. Givental, B. A. Khesin, J. E. Marsden, A. N. Varchenko, V. A. Vassiliev, O. Y. Viro, and V. M. Zakalyukin (Springer-Verlag, Berlin, Heidelberg, 2009), pp. 423–427.

[36] B. Autin, Nonlinear betatron oscillations, AIP Conf. Proc. **153**, 288 (1987).

[37] J. Bengtsson, Non-linear transverse dynamics for storage rings with applications to the low-energy antiproton ring (LEAR) at CERN, Ph.D. thesis, Lund University, Geneva, 1998.

[38] J. Laskar, The chaotic motion of the solar system: A numerical estimate of the size of the chaotic zones, Icarus **88**, 266 (1990).

[39] J. Laskar, C. Froeschlé, and A. Celletti, The measure of chaos by the numerical analysis of the fundamental frequencies. Application to the standard mapping, Physica (Amsterdam) **56D**, 253 (1992).

[40] Y. Papaphilippou, Detecting chaos in particle accelerators through the frequency map analysis method, Chaos **24**, 024412 (2014).

[41] A. A. N. Kolmogorov, Preservation of conditionally periodic movements with small change in the Hamilton function, in *Stochastic Behavior in Classical and Quantum Hamiltonian Systems*, edited by G. Casati and J. Ford (Springer-Verlag, Berlin, 1979), pp. 51–56.

[42] J. Möser, On invariant curves of area-preserving mappings of an annulus, Nachr. Akad. Wiss. Göttingen, II, **1**, 1 (1962), https://cds.cern.ch/record/430015.

[43] Proof of a theorem of A. N. Kolmogorov on the invariance of quasi-periodic motions under small perturbations of the hamiltonian, in *Collected Works: Representations of Functions, Celestial Mechanics and KAM Theory, 1957–1965*, edited by A. B. Givental, B. A. Khesin, J. E. Marsden, A. N. Varchenko, V. A. Vassiliev, O. Y. Viro, and V. M. Zakalyukin (Springer, Berlin, Heidelberg, 2009), pp. 267–294.

[44] A. Chao, Lecture notes on topics in accelerator physics, Report No. SLAC-PUB-9574, 2002.

[45] E. Forest, M. F. Reusch, D. L. Bruhwiler, and A. Amiry, The correct local description for tracking in rings, Part. Accel. **45**, 65 (1994).