{'fraction_non_alphanumeric': 0.029069620092759207, 'fraction_numerical': 0.011455508528271282, 'mean_word_length': 4.69626100765122, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High order time integrators for the simulation of charged particle motion in magnetic quadrupoles

Magnetic quadrupoles are essential components of particle accelerators like the Large Hadron Collider. In order to study numerically the stability of the particle beam crossing a quadrupole, a large number of particle revolutions in the accelerator must be simulated, thus leading to the necessity to preserve numerically invariants of motion over a long time interval and to a substantial computational cost, mostly related to the repeated evaluation of the magnetic vector potential. In this paper, in order to reduce this cost, we first consider a specific gauge transformation that allows to reduce significantly the number of vector potential evaluations. We then analyze the sensitivity of the numerical solution to the interpolation procedure required to compute magnetic vector potential data from gridded precomputed values at the locations required by high order time integration methods. Finally, we compare several high order integration techniques, in order to assess their accuracy and efficiency for these long term simulations. Explicit high order Lie methods are considered, along with implicit high order symplectic integrators and conventional explicit Runge Kutta methods. Among symplectic methods, high order Lie integrators yield optimal results in terms of cost/accuracy ratios, but non symplectic Runge Kutta methods perform remarkably well even in very long term simulations. Furthermore, the accuracy of the field reconstruction and interpolation techniques are shown to be limiting factors for the accuracy of the particle tracking procedures.


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 10 5 ). 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 [2]. 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 z 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 z and harmonics of higher order, see e.g. [7]. 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 [12]. 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 [14] has been carried out in [5] 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 [14] are considered along with im-plicit high order symplectic integrators [8] 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.

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 Since the magnetic field B is irrotational and therefore conservative on simply connected domains, there exist a magnetic scalar potential ψ and a magnetic vector potential A such that B = ∇ × A = ∇ψ. 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 A gauge transformation A = A+∇λ 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 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. [3], [10]. However, in the approach most commonly used in accelerator physics, see e.g. [6], [13], the magnetic vector potential A 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 q(t) = (x, y, z) T and its canonical momentum p(t) = (p x , p y , p z ) T . 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 where m denotes the rest mass of the particle, c is the speed of light and Q the particle charge [11]. The mechanical momenta are denoted by p M = mγv, where v is the particle speed and γ = 1/ » 1 − v 2 /c 2 is the Lorentz factor. The canonical momenta p are related to the mechanical momenta through: Introducing the state vector w = Ä q T , p T ä T , the relativistic Hamilton equations can be written aṡ Here we denote J q A (q) = ∂A/∂q and we set where 0, I are the zero and identity matrix, respectively. For the specific problem at hand, it is convenient to assume that the z axis is the symmetry axis of the quadrupole and to use the longitudinal coordinate z as independent variable instead of time. This change of independent variable leads to the Hamiltonian where now p t = −H is the conjugate momentum of t. The dynamical variables are w = (x, y, t, p x , p y , p t ) T and z 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 z. The magnetic field along the center of a quadrupole is null. Therefore, a particle traveling along the z axis with speed v 0 z is not influenced by the quadrupole. In the following we will indicate the quantities related to the reference particle using the superscript 0. The trajectory described by this particle is called the reference orbit and it is identified by It is convenient to describe the motion of a general particle in terms of deviations with respect of the reference orbit w d = w−w 0 . 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 L, 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 L = 1 m and we will use as a reference the momentum of a proton with rest mass m 0 ∼ 9.38 MeV/c 2 ∼ 1.67 · 10 −27 kg and total energy E = 7 TeV ∼ 1.12 · 10 −6 J. Therefore the reference momentum is given by: The scaled, non dimensional variables are denoted by X, Y, τ and P x , P y , P τ , respectively, where τ = ct d /L, P τ = p d t /(p 0 c). Another important quantity that will be used in this work is the mechanical momentum deviation which is related to P τ by the relation Replacing the canonical pair (τ, P τ ) by ( , δ) , one obtains the Hamiltonian where ‹ A (X, Y, Z) = Q p 0 A (L X, L Y, L Z). 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: 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 Z, this Hamiltonian describes a 6 dimensional non-autonomous system. To obtain an autonomous system, it is possible to introduce a new canonical pair (Z, P z ) and a new independent variable σ: where q = (X, Y, , Z) and p = (P x , P y , δ, P z ). In this case, the resulting Hamilton equations are: Moreover, it can be noticed that the Hamiltonian does not depend on , so that the partial derivative of K with respect to is zero. As a consequence, δ is a constant of motion, equal to the initial value, denoted by the subscript 0, δ 0 . If the evolution of the variable is not needed, the canonical pair ( , δ) can be neglected, considering δ 0 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 w = (X, Y, Z, P x , P y , P z ), the Hamilton equations become: A further simplification can be achieved noticing that P z 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 X, Y , P x and P y .
When this property is preserved by the numerical method, i.e., if the one-step map Φ ∆t : y 0 → Φ ∆t (y 0 ) = y 1 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 . . , s) be real numbers and let c i = s j=1 a ij . A s-stage Runge-Kutta method is given by: 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 a ij = 0 for j ≥ i. The following theorem gives a sufficient condition for a Runge-Kutta method to be symplectic [8].
Theorem 1 If the coefficients of a Runge-Kutta method satisfy: then the method is symplectic.
Gauss methods are particular implicit Runge-Kutta methods, some of which satisfy the condition of theorem (1) 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 [8].
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 H is the Hamiltonian ruling the evolution of an autonomous system, then the exact solution of the Hamilton equations can be formally represented as Here, the notation of equation (4) is used, : : denotes the Lie operator and the exponentiation of a Lie operator is called Lie transformation [6]. 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 [15] and further discussed in [14].
Denoting by M 2 (∆t) the approximation (15) and by M 2n (∆t) an approximation of order 2n, an approximation of order 2n + 2 can be built as follows In this work, methods of order 4 and 6 have been considered, with (α 0 , α 1 ) pairs given by

respectively.
Taking into account the discussion in section 2, the map M 2 , applied to w n yields the following algorithm 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 Z are available. On the other hand, all the methods introduced require the magnetic vector potential evaluation at Z values different from the sampled ones. For some methods, like the midpoint and explicit Runge Kutta method, only the evaluation at Z = Z n + ∆Z/2 is required, so that interpolation of the sampled data can be avoided if a ∆Z 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 A 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 [4]. referred to [6] 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 (ρ, ϕ, z). Due to the periodicity of the field in the angular variable ϕ, it is then possible to expand the angular dependence using Fourier series The field harmonics A m and B m 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 m are different from zero, in particular those with m = 2 (2j + 1) for j ≥ 0. The magnetic scalar potential ψ satisfies the Laplace equation (2). This allows to derive from (18) the representation where C [2 ] m, s (z) (and C [2 ] m, c (z)) are known functions called normal (skew) generalized gradients. The radial component of the harmonics A ρ, m and B ρ, m (see equation (18)), measured at a certain radius R an , also known as radius of analysis, are denoted from now on simply by A m and B m . They can be used to compute the generalized gradients using the following formula Here, I m denotes the modified Bessel function of the first kind, while " A m (R an , k) and " B m (R an , k) denote the Fourier transforms of A m (R an , z) and B m (R an , z), 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 (s) or skew (c). Also the magnetic vector potential be defined using normal and skew terms: Using the generalized gradients, it is possible to derive the expression for a first vector potential gauge, called azimuthal-free gauge, for which A ϕ ≡ 0: If we require that a vector potential is divergence-free ∇ · " A = 0, we obtain the so called Coulomb gauge. The symmetric Coulomb gauge " A belongs to this category and can be expressed as The x and the y components can be written as follows where ‹ m = m + 1 and Finally, via a gauge transformationĀ = " A + ∇λ, a new form of the vector potential can be derived, such thatĀ x ≡ 0. The derivation of this so called horizontal-free Coulomb gauge is described in detail [6] and summarized in the following. As it will be shown later in this section, the propertyĀ x ≡ 0 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).
Notice that a gauge transformation is equivalent to a canonical transformation, see e.g. [9]. In particular, if we consider a Hamiltonian K, similar to (9), with vector potential A and dynamical variables q and p then, using a new vector potential A = A + ∇λ we obtain a Hamiltonian K in the same form of K with dynamical variables: Q = q ; where Q denotes again the particle charge. To derive the horizontal free Coulomb gauge transformation, we build a harmonic function λ as where m, s (z)ρ 2 +m , and the coefficients L [2 ] m, s/c (z) are related to the coefficients C [2 ] m, s/c (z) by the following relations and x leading to the desired horizontal-free Coulomb gauge: where and ∂ z λ is obtained from (26): All the previous vector potential descriptions, when expressed in Cartesian coordinates, take the form where the coefficients a i, j (z) depend on the longitudinal coordinate z. In practical cases, the series expansions in the previous formulae are truncated to a finite number of terms, which depends on the number of harmonics used to describe the field, while the maximum number of generalized gradients derivatives N D determines the range of the indices i and j in (31). Therefore, for a given z = z value, the evaluation time of (31) is proportional to the number of coefficients a i, j ( z) retained, which can therefore be used to estimate the computational cost entailed by each representation.
Using the horizontal-free gauge, all the terms in the previous series expansions corresponding to the x component of the vector potential are null, while the number of terms corresponding to the other two components is similar in the Coulomb and horizontal free gauges. The number of coefficients required by the magnetic vector potential at a specific z location, in different gauges, using the harmonics m ∈ {2, 6, 10, 14} and different generalized gradient derivatives are shown in It can be noticed that the horizontal-free Coulomb gauge requires in general between 20% and 25% less coefficients with respect to the azimuthal-free gauge.
The effective reduction in computational cost achievable by this reformulation is however dependent on the specific features of each numerical method. Indeed, Runge-Kutta methods require the evaluation of the right hand side of Hamilton equations J ∇K, while Lie methods evaluate the map M 2 . Therefore, the speed-up of Runge-Kutta and Lie methods, achieved with the different gauges, is related to the number of vector potential coefficients required by J ∇K and M 2 , respectively. In table 2, the number of evaluations of each vector potential component required by J ∇K and by M 2 are shown. Moreover, in the last columns, the ratios between the number of vector potential coefficients using the horizontal-free Coulomb gauge and the azimuthal-free gauge are shown, when considering harmonics m ∈ {2, 6, 10, 14} and N D = 16. It can be noticed that the map M 2 , and therefore Lie methods, benefit more from the change of gauge, because they require more evaluations of the A x component. The overall efficiency of each method also de-  pends on the total number of evaluations of J ∇K or M 2 required at each step. These are reported in table 3. Moreover, it is important to notice that implicit methods require a certain number of fixed point iterations, between 5 and 8, that increase significantly their computational cost per time step. For a comprehensive efficiency comparison, the total computational cost of each method must be compared to the accuracy level it allows to achieve and to the accuracy level that is actually required for accelerator design. This comparison will be attempted in sections 5 and 6. In the following sections, only normal quadrupoles will be considered (i.e. with null skew generalized gradients), so that the notation will be simplified using C m,s/c (z) . If N D is an even number, all the vector potential gauges produce exactly the same magnetic field, which satisfies the Maxwell equation

Number of evaluations Ratio HFC/AF
where j SP is a spurious current that depends on the generalized gradients only through the generalized gradient derivatives of order N D + 1. This implies that, generally, the largest contributions to these currents will be concentrated in the fringe field. For example, using only two derivatives of C [0] 2 (Z), the spurious currents are: The magnitude of this error can be measured using the following quantity: max where (X 0 , Y 0 ) refer to a fixed position in the transversal plane.
In figure 1, it is reported the behaviour of the error (33) for two different transversal positions, one closer to the quadrupole axis and another further away. The results refer to vector potentials in different gauges, computed considering a second order harmonic B 2 coming from a realistic design of a quadrupole (see section 6). It can be noticed that the error in the different gauges coincides if N D is an even number, coherently with the fact that the magnetic fields are identical. It is unclear to which extent this approximation error, which is far from negligible at larger distances from the quadrupole axis, entails a limitation for the accuracy that tracking algorithms can achieve. Mimetic methods like those presented in [10] guarantee that a discrete equivalent of ∇ × ∇ × A = 0 is exactly satisfied, so that they could represent a potentially interesting alternative to the techniques described in this section. As a final remark we notice that, in the use of equation 20, particular attention has to be devoted to the computation of the input harmonics. If harmonics that do not go to zero sufficiently fast at the boundaries are employed, significant errors in the field description may result. For example, even if the harmonics vanish at the boundaries, the computed generalized gradients, and thus the vector potential, may not vanish as well, causing a discrepancy between the canonical and the mechanical momenta in a region where they should coincide. In our experience, these issues can be avoided by appropriate extension of the harmonics data by sufficiently large regions filled with zero values.

Numerical experiments with an analytic magnetic vector potential
In this section, equations (11) will be solved for the case of a simple vector potential, whose expression is given by a polynomial with just a few non zero coefficients. It is important to notice that this configuration is not completely realistic, since this vector potential does not exactly satisfy the Maxwell equation ∇ × ∇ × A = 0. As discussed in the previous section, this leads to spurious residual currents in the right hand side of this equation, which are however not accounted for in equations (11). In spite of this inconsistency, this benchmark is useful to assess the impact of the interpolation procedure mentioned in section 3, since the simple analytic expression can also be computed cheaply online and comparison between results obtained with and without interpolation can be carried out. A set of scaled vector potentials in different gauges is built using equations (22), (23) and (28) and employing only the C where erf denotes the error function, defined as Using (34), the function C [0] 2 (Z) is then defined as Here, the coefficient α is the plateau value of C to 6 · 10 −4 with L 1 = L 2 = 0.9 and Z 2 = 3.1. The initial conditions for the transversal positions and momenta of the particle are set to w 0 = (0.02, −0.04, 0, 0) and δ 0 = 0. A reference solution is computed using the exact vector potential and the MATLAB ODE solver ode45 with a maximum ∆Z ref = ∆Z data /10 and a relative error tolerance of 10 −13 .
In figure 3, the errors obtained with the different ODE methods are presented, in the case in which no interpolation is used and the magnetic vector potential is computed exactly at each required location. The errors on the X, Y coordinates behave very similarly and the same is true for the corresponding moments, but the errors on transversal momenta are larger than those on positions. As long as the errors are above the tolerance chosen for the reference solution, both second and fourth order solvers behave in agreement with theoretical expectations, while sixth order solvers only seem to display the expected error decay for sufficiently small values of the interval ∆Z. Furthermore, very similar errors are obtained for second and fourth order methods for the largest ∆Z employed. We attribute these facts to the poor resolution of the larger gradient areas at the beginning and at the end of the idealized quadrupole. This effect can also be seen in the numerical results of section 6. It gives a clear indication that there is an upper limit to the value ∆Z that can be employed, independently of the accuracy of the solver employed. In order to study the impact of different interpolation techniques in the more realistic case in which the vector potential is only available as sampled data, the analytic vector potential obtained from (35) are sampled on a fine mesh with ∆Z data = 0.002 over the interval [0, 4]. For each given position along the Z axis where the vector potential is not available, we employ a) the value at the last previously available gridded location (previous), b) the value at the nearest gridded location available (nearest), c) the average of the nearest available potential values (interval) d) a cubic spline interpolation (spline). The results of this comparison are reported in figure 4. Only results obtained with the sixth-order Gauss method are displayed, since the highest order methods are those most af-fected by the accuracy of the field interpolation. It can be noticed that all interpolation methods limit the overall accuracy of the time integrator, to a larger or lesser extent. Ideally, an interpolator of the same order of the time integration method should be employed. On the other hand, spline interpolation seems to be sufficient to achieve errors comparable to those of the exact potential evaluation in most cases. In order to compare the efficiencies of the methods employed, we report in figure 5 the behavior of the error as a function of the CPU time required by each method for a given resolution. Since the interpolation stage is done off-line, the CPU time required does not depend on the interpolation method employed. Among the ODE methods, the fourth-order explicit Runge-Kutta method gives the best results, followed by the sixth-order Gauss method and the fourth-order Lie method. Among the symplectic methods, the fastest method is the second-order Lie one, but it has a relatively low accuracy.
As mentioned in the introduction, the beam stability assessment requires long-term simulations. In order to verify the solution quality in this framework, the phase-space orbits of a particle which travels through a sequence of focusing-defocusing quadrupole cou- ples are measured. This test is carried out considering a sequence of 3000 consecutive quadrupole couples and an integration step of ∆Z = 0.08. All the numerical methods produced stable orbits in this test. As an example, in figure 6 the results for the fourth-order Lie method are reported. Remarkably, the non-symplectic fourthorder Runge-Kutta method gives results entirely analogous to those of the other (symplectic) methods considered, see figure 7.
6 Numerical experiments with a realistic vector potential In this section, equations (11) will be solved for the case of a magnetic vector potential that corresponds to the design of a realistic quadrupole. The harmonics m = 2, 6, 10, 14, at a radius of analysis of 0.05 and sampled at ∆Z data = 0.02 (figure 8) are provided by numerical simulations performed by the FEM/BEM software ROXIE [1], used at CERN to design the accelerator magnets. The presence of connectors on one side of the quadrupole causes asymmetries in the harmonics along the Z axis and the presence of skew harmonics. The generalized gradients are computed using up to 16 derivatives. We will consider the case of a particle quite far from the quadrupole axis: indeed, the initial conditions for the transversal positions and momenta are set to w 0 = (0.02, −0.04, 0, 0) and δ 0 = 0. Using the horizontal-free Coulomb gauge, we obtain a gain in efficiency quite in agreement with the estimates given in section 4, see tables Normal harmonics sampled at ∆z data = 0.02 m. 4 and 5. In the following, only the results with the horizontal-free Coulomb gauge and using cubic spline interpolation of gridded data are reported.
The error comparisons are shown in figure 9. In order to measure the methods' accuracy, the absolute error on the positions and the momenta at the quadrupole exit are used, because these values are useful to to understand if the non-linear effects, mentioned in the introduction, can be described correctly. The trend is similar to those seen in the previous sections, but with a general worsening in accuracy. Longer time integration steps prevent the correct description of the field and consequently the achievement of the theoretical convergence orders for the numerical methods. The final error drops in the high-order methods can be due to the fact that the reference solution has been computed using the same starting vector potential data and only a finer mesh of interpolated values. In order to compare the efficiencies of the methods employed, we plot in figure  10 the behaviour of the error as a function of the CPU time required     In this case, the large ∆Z data of the input data limits the possibility of achieving high accuracies, therefore there is no clear advantage in the use of high order methods and the second-order Lie method is the most efficient.
Finally, we check the beam stability looking at the phase-space orbits of a particle which travels through a sequence of focusingdefocusing quadrupole couples. This test is carried out considering a sequence of 24000 consecutive quadrupole couples and an integration step of ∆Z = 0.04. In this case it is possible to notice (see figures 11 and 12) that the trend in the (X, P x ) phase space appears unstable.
The instability is probably due to the high degree of the vector potential polynomial, obtained using various harmonics and many generalized gradient derivatives. In fact, using only the second order harmonic and two generalized gradient derivatives (figure 13), or considering a motion closer to the center (figure 14), it is possible to obtain stable orbits similar to the ones obtained in the analytical case presented in section 5. It is especially noteworthy that, in spite of its lack of symplectic properties, the classical fourth order RK method does not appear to behave differently with respect to the symplectic ones (figures 15 and 16).
Finally, we check the energy conservation for a particle which travels through 8000 focusing-defocusing quadrupole pairs using an integration step ∆Z = 0.08. Results are shown in figure 17. In this realistic case the magnetic vector potential has a more complex shape and the differences in the K X trend among the numerical methods appear sooner. Nevertheless, the K X oscillations appear to be stable and no significant differences between symplectic and non symplectic methods concerning the energy conservation are observed.   Figure 17: Realistic case II, cubic spline interpolation. Trend of the K X (left) energy component along 8000 focusing-defocusing quadrupole couples for different ODE methods and deviation of K X associated to the different methods with respect to the sixthorder Gauss method (right).

Conclusions and future developments
We have discussed and analyzed several issues that arise in the numerical approximation of charged particles trajectories in magnetic quadrupoles. We have shown that 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. The impact of the interpolation technique employed to compute magnetic vector potential values at arbitrary locations from gridded data has also been assessed, highlighting potential limitations in accuracy induced by insufficiently accurate interpolation methods, if high order time integration techniques are to be applied. However, cubic spline interpolation was found to be sufficient for preserving the accuracy of all the methods considered in this work over a wide range of values for the integration step.
We have then compared several high order integration techniques, which allow to maintain high accuracy even with relatively large integration step values, in order to assess their accuracy and efficiency for long-term simulations. Explicit high order Lie methods have been considered, along with implicit high order symplectic integrators and a more conventional, non symplectic explicit Runge-Kutta methods.
In the case of realistic vector potentials, the errors induced by the vector potential representation and interpolation become significant and reduce the highest possible accuracy that can be attained. Furthermore, since in realistic cases the magnetic vector potential evaluation costlier, numerical methods which require less evaluations, such as the second-order Lie method, appear to be more competitive in terms of efficiency. On the other hand, if these errors could be reduced by different approaches for the vector field representation, higher order methods could be advantageous if a more precise approximation is required. The speed gain obtained by the horizontal-free Coulomb gauge would also allow to use more expensive methods. 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.
All the computations carried out in this work employed vector potential values sampled on a uniform grid along the longitudinal axis of the quadrupole. In order to increase the efficiency of the methods employed, using vector potentials sampled over a non-uniform grid appears to be a straightforward and necessary improvement. In such way, it could be possible to use a smaller step in the fringe field regions, where the vector potential has more complex shape, and a larger step where the field is uniform and has simpler structure.
A particularly interesting aspect of the results obtained was the fact that non symplectic methods appeared to be competitive with symplectic ones, even on relatively long integrations. An even more extreme example of this behaviour was seen in the simulation of a particle crossing 1600000 quadrupole couples of realistic design (32 quadrupoles in 100000 revolutions). An integration step ∆Z = 0.16 was employed. As shown in figure 18, even in this extremely long simulation the behaviour of a symplectic Lie method and of a non symplectic Runge Kutta method was entirely analogous in terms of energy conservation, in contrast to the theoretical expectations. It 0 500,000 1,000,000 1,500,000 2,000,000 2,500,000 3,000,000 3,500,000 4,000,000 4,500,000 5,000,000 5,500,000 6,000,000 6,500,000 7,000,000 0 K X trend: cubic spline interpolation lie4 rk4 Figure 18: Energy conservation of a particle through 1600000 focusing-defocusing quadrupole couples using the fourth-order explicit Runge-Kutta and the fourthorder Lie method.
would be interesting to understand whether more realistic test cases could actually highlight negative effects induced by non-symplectic methods and whether more detailed monitoring of the particle motion could be employed for this purpose. For example, in this work only the energy associated to the transverse variable was used to measure the good behaviour in long-term simulations. Other conserved quantities could also be used to compare the performance of the different methods.