## 1 Introduction

Magnetic quadrupoles are key components of particle accelerators that are used to focus particle beams. In high energy circular accelerators, the quality of their magnetic fields can influence the overall beam stability. The beam stability is measured in terms of the dynamic aperture, defined as the region in the phase space outside which a particle is considered as lost from the beam. The dynamic aperture is estimated by solving numerically the Hamilton equations describing the particle trajectories for a large number of accelerator revolutions (typically, more than

). Therefore, it is important to use very efficient time integration methods which also guarantee good long-term preservation of dynamical invariants of motion. The action of a quadrupole can be approximately described using a linear combination of the position and momenta of the particles at the inlet that yields the position and momenta at the outlet Carey (1988). This approach yields a good approximation if the particles travel near the quadrupole center (small apertures) and it is exact in regions where the magnetic field is constant along the longitudinal axis and only the main quadrupole field is present. In realistic cases, however, the magnetic field has a more complex structure which involves non-uniformity along and harmonics of higher order, see e.g. Forest (1998). These inhomogeneities of the field at the quadrupole ends, known as fringe field, can lead to a non-linear dependency of the position and momenta of the particles at the outlet from the position and momenta at the inlet. In this paper, we focus on the numerical problems encountered when modelling these non linear dependencies in an accurate and efficient way, as necessary for the design of the large aperture quadrupoles foreseen for the HL-LHC project Rossi (2011). The accurate numerical solution of the complete Hamilton equations is mandatory in this case. A preliminary study on the applicability of second-order methods based on the Lie algebra integrators proposed in Wu et al. (2003) has been carried out in Dalena et al. (2014) for the case of a realistic quadrupole.In this work, we first consider a specific gauge transformation that allows to reduce by approximately 50% the computational cost of each vector potential evaluation, thus significantly enhancing the efficiency of any numerical approximation method employed for the particle trajectory simulation. We then compare several high order integration techniques, which allow to maintain sufficiently high accuracy even with relatively large integration step values, in order to assess their accuracy and efficiency for these long-term simulations. Explicit high order Lie methods Wu et al. (2003) are considered along with implicit high order symplectic integrators Hairer et al. (2006) and more conventional, non symplectic explicit Runge-Kutta methods.

In the case of realistic vector potentials, the errors induced by the vector potential reconstruction and interpolation become significant and reduce the highest possible accuracy that can be attained. Furthermore, since in realistic cases the magnetic vector potential evaluation is more costly, numerical methods which require less evaluations, such as the second-order Lie method, appear to be more competitive. On the other hand, experiments with idealized fields show that, if these errors could be reduced, higher order methods could be advantageous and the speed gain obtained with the horizontal-free Coulomb gauge would enhance their efficiency. In particular, the explicit fourth-order Runge-Kutta appears to be the most efficient method and the fourth-order Lie the most efficient among symplectic methods. A particularly interesting aspect of the results obtained is the fact that non symplectic methods appear to be competitive with symplectic ones even on relatively long integrations, when stability of the computed trajectories and energy conservation are considered. Indeed, the spurious energy losses appear to be more closely related to the errors in the representation of the magnetic vector potential than to those introduced by the time integrators. This unexpected result warrants more detailed investigation.

A detailed outline of the paper is as follows. In section 2, we introduce the Hamiltonian that describes the motion of a charged particle inside a magnetic quadrupole and the corresponding Hamilton equations. In section 3, the numerical methods used in this work are briefly reviewed, along with the specific issues that arise when the vector potential values are only available at given sampling intervals, so that appropriate interpolation procedures must be employed if high order methods are to be applied. In section 4, we review the approach employed to represent the magnetic vector potential. We also show how the non-uniqueness of the vector potential can be exploited to identify a gauge that allows to reduce the number of vector potential evaluations. Numerical results for a simple vector potential that can be expressed analytically are presented in section 5. An assessment of time integration methods on a realistic case is presented in section 6. Finally, some conclusions are drawn in section 7, where we also discuss possible developments of this work.

## 2 Charged particle motion in magnetic quadrupoles

Magnetic quadrupoles are devices in which a stationary magnetic field with cylindrical symmetry is generated. Near the quadrupole center, where the particles travel, the field is a solution of the Maxwell equations in the vacuum and in absence of charges and currents

(1) | ||||||

Since the magnetic field is irrotational and therefore conservative on simply connected domains, there exist a magnetic scalar potential and a magnetic vector potential such that The magnetic scalar potential is defined up to constants and the magnetic vector potential is defined up to gradients of scalar functions, so that given a scalar function and defining , one has A gauge transformation does not change the magnetic field but can yield a more convenient representation of the vector potential for specific purposes. The second Maxwell equation can be rewritten as

(2) |

which means that the scalar potential satisfies the Laplace equation. Many high accuracy numerical methods are available to solve this equation starting from appropriate boundary conditions and taking into account also the real geometry of the accelerator, see e.g. Corno et al. (2016), Lipnikov et al. (2011). However, in the approach most commonly used in accelerator physics, see e.g. Dragt (1997), Venturini and Dragt (1999), the magnetic vector potential is represented by a power series in the transversal coordinates, whose coefficients are functions of the longitudinal coordinate. A detailed presentation of this approach is given in section 4.

The motion of a charged particle in a magnetic quadrupole is described in terms of its position and its canonical momentum . In the case of a high energy particle accelerator, the particles speed is very close to the speed of light. The relativistic Hamiltonian is given by

(3) |

where denotes the rest mass of the particle, is the speed of light and the particle charge Panofsky and Phillips (1962). The mechanical momenta are denoted by where is the particle speed and is the Lorentz factor. The canonical momenta are related to the mechanical momenta through:

Introducing the state vector , the relativistic Hamilton equations can be written as

(4) | ||||

Here we denote and we set

(5) |

where

are the zero and identity matrix, respectively. For the specific problem at hand, it is convenient to assume that the

axis is the symmetry axis of the quadrupole and to use the longitudinal coordinate as independent variable instead of time. This change of independent variable leads to the Hamiltonianwhere now is the conjugate momentum of . The dynamical variables are and is now the independent variable. It can be noticed that this new Hamiltonian does not describe an autonomous system, due to the fact that the magnetic vector potential depends on the independent variable .

The magnetic field along the center of a quadrupole is null. Therefore, a particle traveling along the axis with speed is not influenced by the quadrupole. In the following we will indicate the quantities related to the reference particle using the superscript . The trajectory described by this particle is called the reference orbit and it is identified by , and . It is convenient to describe the motion of a general particle in terms of deviations with respect of the reference orbit . This change of variables is a canonical transformation which leads to the new Hamiltonian

Moreover, it is convenient to introduce a specific scaling of the deviation variables. In particular, all the position variables will be scaled by a fixed length usually denoted as the bunch length, while the momentum variables are scaled by the module of the reference mechanical momentum

In this work, we will use a reference length and we will use as a reference the momentum of a proton with rest mass and total energy . Therefore the reference momentum is given by:

(6) |

The scaled, non dimensional variables are denoted by and respectively, where Another important quantity that will be used in this work is the mechanical momentum deviation

which is related to by the relation

Replacing the canonical pair by one obtains the Hamiltonian

(7) |

where . For relativistic particles, the momenta in the transversal plane are much smaller than the total momentum module, i.e.

Therefore, the so called paraxial approximation can be introduced, that amounts to substituting the square root in equation (7) with its Taylor expansion truncated at the first order:

(8) |

where the constant terms can be neglected because they do not influence the Hamilton equations. In general, due to the fact that the vector potential depends on the independent variable , this Hamiltonian describes a dimensional non-autonomous system. To obtain an autonomous system, it is possible to introduce a new canonical pair and a new independent variable :

(9) |

where and . In this case, the resulting Hamilton equations are:

(10) |

Moreover, it can be noticed that the Hamiltonian does not depend on so that the partial derivative of with respect to is zero. As a consequence, is a constant of motion, equal to the initial value, denoted by the subscript , If the evolution of the variable is not needed, the canonical pair can be neglected, considering as a parameter and reducing again the size of the phase space. In this case, the Hamiltonian is still given by (9) but, since the dynamical variables are now , the Hamilton equations become:

(11) |

A further simplification can be achieved noticing that is decoupled from the other dynamical variables, so that its computation can be neglected if we are only interested in the dynamics of the transversal variables, reducing the number of equations (11) to the four ones associated to , , and .

## 3 Review of high order numerical methods for ODE problems

In this section, some high order numerical methods for the solution of a first order system will be reviewed, in view of their application to the solution of the Hamilton equations. A more detailed presentation of the relevant numerical methods can be found for example in Hairer et al. (2006).

An important feature of Hamiltonian flows is their symplectic property, which can be defined more precisely as follows. A differentiable map , where is an open set, is called symplectic if:

where is the Jacobian matrix of and is the matrix (5). When this property is preserved by the numerical method, i.e., if the one-step map is symplectic, quadratic invariants of motion are preserved, thus ensuring in principle a good behaviour for long-term simulations. The first symplectic techniques that will be considered are Runge-Kutta methods, which can be written in general as Let be real numbers and let . A s-stage Runge-Kutta method is given by:

(12) |

Runge-Kutta methods are often summarized via the so called Butcher tableau, in which all the coefficients are arranged as:

A Runge-Kutta method is explicit if for The following theorem gives a sufficient condition for a Runge-Kutta method to be symplectic Hairer et al. (2006). If the coefficients of a Runge-Kutta method satisfy:

(13) |

then the method is symplectic. Gauss methods are particular implicit Runge-Kutta methods, some of which satisfy the condition of theorem (3) and are thus symplectic. The midpoint method can be interpreted as the second-order Gauss method, characterized by the Butcher tableau

The fourth-order Gauss method, considered in this paper, is characterized by the Butcher tableau

The sixth-order Gauss method is instead characterized by the Butcher tableau

Implicit methods require the solution of a nonlinear system of equations at each time step. This can be done using either the Newton or the fixed-point method. Even though Newton’s method is usually superior, numerical results show that the latter is faster. This behaviour can be justified by the fact that, for the problems at hand, both methods require a small number of iterations to achieve convergence, but the Newton method implies higher initial costs related to the evaluation of Jacobian matrix, see e.g. the discussion in Hairer et al. (2006). Also the best known fourth-order explicit Runge-Kutta method has been considered in this work. This method is not symplectic and it is characterized by the Butcher tableau

If is the Hamiltonian ruling the evolution of an autonomous system, then the exact solution of the Hamilton equations can be formally represented as

(14) |

Here, the notation of equation (4) is used, denotes the Lie operator and the exponentiation of a Lie operator is called Lie transformation Dragt (1997). The methods based on Lie algebra techniques most widely applied in accelerator physics employ a second-order approximation of the Lie transformation. Higher order Lie methods are then built using the procedure introduced by Yoshida in Yoshida (1990) and further discussed in Wu et al. (2003). The first step is to split the Hamiltonian in solvable parts

such that can be computed exactly for . This is true if is nilpotent of order two (i.e. for ), because in this case the exponential series reduces to a finite sum. A second order approximation is then given by

(15) | ||||

Denoting by the approximation (15) and by an approximation of order , an approximation of order can be built as follows

(16) |

where and . In this work, methods of order and have been considered, with pairs given by

respectively.

Taking into account the discussion in section 2, the map , applied to yields the following algorithm

(17) |

In the case of particle motion inside a magnetic quadrupole, the ODE system is given by (11). The magnetic vector potential is written in the form (31) and, in many practical applications, only its sampled values at equally spaced locations in are available. On the other hand, all the methods introduced require the magnetic vector potential evaluation at values different from the sampled ones. For some methods, like the midpoint and explicit Runge Kutta method, only the evaluation at is required, so that interpolation of the sampled data can be avoided if a is employed for computation that is twice that of the data. However, in general an interpolation is needed in order to provide the magnetic vector potential evaluated at the points needed by each specific ODE solver. This introduces a further source of error whose quantification is not a straightforward task. Some proposals to compute intermediate values will be compared in section 5, extending the preliminary results in Dalena and Pugnat (2015).

## 4 Representation of the magnetic vector potential

In this section, the approach used in this work for the reconstruction of the magnetic vector potential will be introduced. The reader is referred to Dragt (1997) for a complete presentation of this technique. Due to the geometry of the quadrupole, it is natural to describe its magnetic field using cylindrical coordinates . Due to the periodicity of the field in the angular variable , it is then possible to expand the angular dependence using Fourier series

(18) |

The field harmonics and are the basis of the vector potential approximation. Exploiting the quadrupole symmetries, it is possible to show that only the harmonics associated to certain values of are different from zero, in particular those with for The magnetic scalar potential satisfies the Laplace equation (2). This allows to derive from (18) the representation

(19) |

where (and

) are known functions called normal (skew) generalized gradients. The radial component of the harmonics

and (see equation (18)), measured at a certain radius also known as radius of analysis, are denoted from now on simply by and They can be used to compute the generalized gradients using the following formula(20) | ||||

Here, denotes the modified Bessel function of the first kind, while and

denote the Fourier transforms of

and , respectively. It is possible to express different quantities, such as the harmonics or the magnetic potentials, using either the normal or the skew generalized gradients. Depending on which generalized gradient is used, these quantities are labelled as normal () or skew (). Also the magnetic vector potential be defined using normal and skew terms:(21) |

Using the generalized gradients, it is possible to derive the expression for a first vector potential gauge, called azimuthal-free gauge, for which :

(22) |

If we require that a vector potential is divergence-free we obtain the so called Coulomb gauge. The symmetric Coulomb gauge belongs to this category and can be expressed as

(23) |

The and the components can be written as follows

(24) | ||||

where and

Finally, via a gauge transformation a new form of the vector potential can be derived, such that The derivation of this so called horizontal-free Coulomb gauge is described in detail Dragt (1997) and summarized in the following. As it will be shown later in this section, the property implies that using this representation for the vector potential leads to a significant reduction in the computational cost of each right hand side evaluation in the numerical solution of system (11).