Efficient Evaluation of Arbitrary Static Fields For Symplectic Particle Tracking

This article describes a method devised for efficient evaluation of arbitrary static magnetic and electric fields in a source free region needed for long time tracking of charged particles. Field values given on the boundary of the region of interest are reproduced inside by an arrangement of hypothetical magnetic or electric monopoles surrounding the boundary surface. The vector and scalar potentials are obtained by summing the contributions of each monopole. The second step of the method improves the evaluation speed of the potentials and their derivatives by orders of magnitude. This comprises covering the region of interest by overlapping spheres, then calculating the spherical harmonic expansion of the potentials on each sphere. During tracking, field values are evaluated by calculating the solid harmonics and their derivatives inside a sphere containing the particle. Software has been developed to test and demonstrate the method on a small particle accelerator. To our knowledge, there is no other method of this kind, allowing long time symplectic integration in general static fields, without simplification.


INTRODUCTION
Long time tracking of charged particles in electromagnetic fields is an important topic in several areas of physics, and it is fundamental for particle accelerators. There are thousands of particle accelerators in the world, and their number is growing each year. They are used for various purposes in particle physics, radiotherapy, industry, as synchrotron light sources, and other applications. Charged particle tracking in accelerators has a long history, and there are lots of software for this purpose. See [1] for a historical overview. Despite the clear interest, before this work, there was not an efficient method available for long time particle tracking in general magnetic field in circular accelerators. Existing methods for general magnetic fields are limited to short or moderateterm tracking due to reasons we will discuss later. The usual methods for long time tracking are not general. They always use some simplification or idealization of the field. These simplifications are often justified, especially for rings with a large circumference. There are however many important cases when these simplifications are not valid. This is often the case for small rings or rings containing special magnetic elements.
A useful method obtains the field values efficiently, in a way that doesn't introduce errors distorting the long time evolution of the simulation. We describe such a scheme for arbitrary static fields, then confirm it on a small particle accelerator containing eight bending magnets. Although we have not yet implemented the program for static electric fields, we expect it to work for that case and should be easy to include into the present implementation. We do not discuss time-varying fields and fields arising from space-charge effects in this paper.
Errors during tracking of charged particles have two main sources. The first comes from the evaluation of the * Lajos.Bojtar@cern.ch fields and the second is due to the method of integration. Symplectic integrators are the preferred choice for long term tracking. They preserve phase space volume and keep the energy error introduced by the integration bounded. Although the truncation error of symplectic integrators is not smaller in general compared to conventional integrators, their long term behavior is much better. The energy error of symplectic integrators oscillates around the exact value while using conventional integrators the error accumulates and tends to grow. This is because the dominant term in truncation error is dissipative for conventional integrators and it is non-dissipative for symplectic ones [2].
Our study addresses the first source of error, due to the evaluation of the fields. Symplecticity is a property of the integrator. Whatever bad approximations of the potentials are given to a symplectic integrator, the integration step will be symplectic, assuming that the derivatives of the potentials are defined. However, if the magnetic or electric fields described by the potentials do not satisfy certain conditions, the result of the tracking will not correspond to the real dynamics of the system, the potentials will represent something else than one intended to simulate.
Any static magnetic or electric field in a region free of sources has zero divergence, zero curl and the components of the fields are harmonic functions satisfying the Laplace equation.
To our knowledge, there is no 3D grid-based interpolation method satisfying all these conditions for general domains. Grid methods, for example [3], don't satisfy the conditions in Eqns. (1)(2)(3) between the grid points. Attempts were made to overcome these short-comings in [4], but then the C 1 continuity of [3] is lost between the cells of the 3D grid.
In the domain of accelerator physics, one approach to deal with the problem of the representation of 3D fields is to take the Fourier transform of the B field in three variables [5][6][7]. Individual magnets can be represented this way but then remains the problem of how to enter and exit from the fields. Our experience showed that cutting the fields even very far from the magnet leads to discontinuities in the field and important energy drift during tracking. We avoid this problem by expanding the entire beam region, which is not a feasible solution with the existing Fourier transform methods. Another problem is that some fields have fast changing localized features. In those cases, the high number of terms required to represent the field make Fourier transform methods impractical. For specific areas like magneto-spheric physics, similar solutions exist [8], however, it cannot always be adopted to particle accelerators for the same reasons.
Surface methods describe the fields on the boundary surrounding the region of interest. Field values on a closed surrounding surface determine the magnetic or electric fields in a source free region. This is true because these fields are harmonic functions satisfying the Laplace equation which has a unique solution for a given boundary condition. A surface method able to express arbitrary static magnetic fields in a way that they satisfy the relevant conditions in Eqns.(1-3) is given in [9]. Since our method is similar to that stage, we give a summary here.
Based on the Helmholtz decomposition theorem, it is possible to describe a general static magnetic field in a source free region as taking the curl of the sum of two vector potentials. One depends only on the tangential component of the B field and the second depends only on the normal component. Both are given on the boundary surface surrounding the region of interest. where is the vector potential of the Dirac monopole with unit charge, r ′ is the location of the monopole, r is the evaluation point and m is a unit vector pointing toward the direction of the Dirac string. To get A n and A t the bounding surface is partitioned to non-overlapping coordinate patches, then Eqns. (5,6) are integrated by high order cubature formulas for each patch and summed together. The vector potential is related to the B field as:

A. Solid harmonics and Taylor expansion
Direct use of this surface method would require integration on the boundary surface for each element in the accelerator at each time step. In theory, one can use this method directly for long time tracking together with a symplectic integrator, in practice this would be extremely slow. We haven't found any paper describing a direct usage of the method, most probably due to performance reasons.
A more localized description of the fields is needed to speed up the evaluation. One possibility is to expand the potentials into Taylor series around some reference trajectory. That is the approach taken in [9] to construct transfer maps for beamline elements. Particle tracking with transfer maps is fast and useful in beamlines. This approach, however, is not suitable for long time tracking as a remark suggests in A.J. Dragt's vast online book [10] on page 27: "In particular, in the context of Accelerator Physics, the long-term goal of map methods is to be able to describe, predict, and control nonlinear properties with the same facility with which we now handle linear properties. Much has been accomplished in this direction, particularly with regard to single-pass systems and short-to-moderate-term behavior in circulating systems." The problem is how to enter and exit the field of a magnet. As we mentioned earlier, any cut of the fields leads to a non-negligible discontinuity, which causes an energy drift during tracking. Our method avoids this problem by not having cuts at all. This is not possible with transfer maps, because they must start and end somewhere. A sizeable part of the book [10] is dedicated to mitigating this problem on a case-by-case basis.
One might imagine an expansion of the potentials in terms of Taylor series in the entire volume of interest in 3 spacial variables to avoid the problem of field cuts. This would lead to a method enormously inferior compared to ours for reasons explained below.
Spherical harmonics form a complete set of orthogonal functions and can approximate functions on a sphere to any precision [11,12]. The components of the vector potential and the scalar potential are harmonic functions. Inside a sphere, they are uniquely determined by their values on the sphere, therefore spherical harmonics approximations of the potentials on a sphere can also approximate them inside. Spherical harmonics scaled appropriately are called solid harmonics. Regular solid harmonics are the canonical representation for harmonic functions inside a sphere [11].
It is important to understand the differences between multipole expansions in terms of Taylor series and terms of solid harmonics, see [13] for a detailed comparison. Multipole expansion of potentials with Taylor series has (ℓ + 1)(ℓ + 2)/2 coefficients for each degree ℓ. The solid harmonics expansion has only (2 ℓ + 1). That makes a huge difference in the number of calculations. The Taylor expansion needs the partial derivatives of the potentials in 3 variables while the solid harmonic expansion is performed by integration on a sphere only in 2 variables. Taylor expansion can represent any function in a volume, but solid harmonics can represent only harmonics functions, exactly what we need.
To approximate the fields in some general volume, we need to cover the volume by overlapping spheres. Any location where a particle can go during the tracking should be in a sphere. The potentials inside each sphere are continuous and satisfy the conditions given in Eqns. (1)(2)(3). However, when the evaluation algorithm has to switch between the spheres at some step, there is an apparent discontinuity of the potentials and their derivatives between the spheres.

B. The effect of discontinuities
As we mentioned in the introduction, symplecticity is a property of the integrator, and it is not affected by the quality of the approximation of the potentials if the potentials and their derivatives are defined everywhere, which is the case. Any point inside the volume of interest belongs to the interior of a sphere covering the volume and for any point inside a sphere the potentials and their derivatives can be calculated analytically with the formulas given later. The effect of discontinuities is an energy drift during the tracking. The drift is caused by the fact that the approximations of the potentials inside the spheres are exact only when the potentials don't contain higher other components than the maximum degree of the solid harmonics expansion. As we limit the solid harmonics expansion at some degree ℓ max there is a small difference between the exact value of the potentials and the approximated ones. In the overlapping region of two neighboring spheres the approximated values are slightly different, depending on which sphere was chosen.
It is possible to reduce this discontinuity as small as needed within the limitations imposed by the numerical precision of the implementation. The higher the degree of the solid harmonics expansion, the smaller are the discontinuities. That is true by assuming that the amplitudes of the coefficients are decreasing with higher degrees. This is usually the case for particle accelerator magnets. FIG.1 shows an exponential decay of the coefficients for a bending dipole magnet as an example. The exponential decrease of the coefficient amplitudes is also a property of other types of magnets, like quadrupoles, octopoles, solenoids, and other kinds of magnets. This example shows that the coefficients with a degree above 40, contribute little. They are around 10 −16 , the numerical resolution of 64-bit floating-point numbers.
One might ask what the difference between the discontinuities in our method and De Haan's approach described in [4], which we discarded for long time tracking due to the discontinuities is. Both methods must keep Each dot corresponds to a particular order (m in c ℓ,m ). We plotted only the positive orders, for this particular example the coefficients belonging to negative orders are zero.
the discontinuities at a very low level. When the potentials are approximated by solid harmonics expansion, the error of the approximation is about the same order of magnitude as the solid harmonics with degree l max + 1 in the expansion. As FIG.1 shows, the coefficients decrease very fast, exponentially. Therefore the discontinuities between the spheres also decrease exponentially as we go higher with the degree of expansion, l max .
In De Haan's 3D grid method the discontinuities can be reduced by decreasing the distance between the grid points. The number of grid points needed to reduce the distance by some factor goes up with the third power of that factor. The error is also decreasing with the third power of the distance between grid points. Additionally, when the grid size is decreased, the number of discontinuities per unit length is increased. Putting all this together we have a slower than linear improvement of the discontinuity errors as the number of grid points is increased. Compare this to the exponential decrease of the spherical harmonics expansion. To keep the effect of the discontinuities at an acceptable level, the number of grid points needed for De Haan's method is too big to be practical.
When approximating harmonic functions, any 3D grid approach is inferior to ours because harmonics functions are uniquely determined by the boundary conditions. In a 3D grid, only the points on the surface of the volume contain useful information, all the grid points inside the volume are redundant.

C. Evaluation of the potentials
The value of the potentials can be evaluated inside a sphere with radius R with the following expression. In fact, this is the well known formula for regular real solid harmonics inside a unit ball, scaled by R.
We denoted the φ dependent part as Φ(φ; m) given by In the expression above,P m ℓ (cos θ) is the orthonormalized associated Legendre polynomial, given as: P m ℓ (cos θ) is calculated without the Condon-Shortley phase. In Eq.(9), r = |r e − r c |/R is the scaled distance between the sphere center r c and the evaluation point r e .
In the actual implementation, the order of the summation is exchanged as [14] recommends it. Also, the multiplication with the coefficients c ℓm for the negative and positive m values are calculated at the same step. This is the usual technique in most optimized codes [15]. It allows a reduction of the number of operations in the m loop by a factor of two.
The real-valued coefficients c ℓm are pre-calculated by the following integral: where ϕ(θ, φ) is the potential to be approximated, and S 2 denotes the surface of the sphere. There are several ways to do this integration. A good comparison of different methods can be found in [16]. Our priority is numerical accuracy, and secondly to minimize the number of evaluations of the potentials ϕ(θ, φ), because the computational cost of calculating the potentials by integration of Eqns. (5,6) or by the method we will describe later, is much higher than the cost of calculating P m ℓ (cos θ). Usually Eq.(12) is evaluated by some numerical quadrature. An elegant way is, to sum values ϕ(θ, φ) on the surface of the sphere at the t-design points with equal weights. A spherical t-design is a set of N points on the sphere, such that a quadrature with equal weights using these nodes is exact for all spherical polynomials of degree at most t. More on t-designs can be found in [17]. An advantage of integrating on a sphere with t-designs is that the sample points has a uniform distribution, contrary to the often used Gaussian quadrature, which has a denser sampling around the poles. As [18] points out, when the sampling points are distributed uniformly on the sphere, the error in the integral is minimized, assuming the errors in the sample values have a normal distribution.
For a set of N points {x i } in a t-design the following expression is exactly true for all spherical polynomials of maximum degree t or below.
Combining Eqns. (12,13), the spherical harmonics coefficients can be calculated with the following sum: where N is the number of quadrature points in the tdesign. The spherical coordinates θ(x i ) and φ(x i ) are expressed in terms of the vectors x i pointing to the quadrature points. The author of [17] has published a set of files containing t-designs up to degree 325 on his website for download. A subset of these files has been used in our implementation.

D. Evaluation of the derivatives
Apart from the values of the potentials, we also need their first derivatives to integrate the equations of motion. The values and the derivatives are calculated in terms of spherical coordinates first, then converted to Cartesian coordinates. Spherical coordinates are the natural choice in a sphere, but we have many of them, so a global coordinate system common to all spheres is needed. The most convenient is to use Cartesian global coordinates. The first derivatives in the spherical coordinate r is directly obtained from Eq.(9) by differentiation as: (16) The derivative of Φ(φ; m) must be treated for the three case listed in Eq. (10).
The derivatives are expressed in terms of orthonormalized associated Legendre polynomial to achieve good numerical precision.
As [14] points out, calculating P m ℓ (cos θ) directly without normalization, produces much more significant numerical errors than the calculation with the orthonormalized version, therefore it is not advised.
We got the θ derivative of the potentials by the following expression: where the θ derivative ofP m ℓ (cos θ) is calculated for the case m = 0 as: with the pre-calculated constants In Eq. (18) there is a singularity due to the division by sin(θ). This is no singularity in the potentials, only a consequence of the spherical coordinate system. It is eliminated by checking θ for zero. In that case, we replace zero with a tiny number. The θ derivative will be huge, but this apparent near singular behavior disappears when the derivatives are converted to Cartesian coordinates. With this treatment of singularity, Eq. (18) gives a precise value only when m = 0. For the case m = 0, we use a different formula given as: Eq. (20) follows from the fact that Please note the minus sign above. It appears, because P 1 ℓ (cos θ) is calculated without the Condon-Shortley phase. The factor l(l + 1) comes from the normalization.
The potentials and their derivatives are evaluated in the same loop. The most demanding part computationally is the calculation ofP m ℓ (cos θ) values, but it can be reused in Eqns. (9,15,16,18,20).
E. An alternative calculation of the potentials To use the evaluation described above, the coefficients c ℓm must be pre-calculated with Eq. (14) for each sphere covering the volume of interest. This requires the calculation of the exact potentials ϕ(x i ). The magnetic potentials can be calculated in certain cases with Eq.(4), which is the method described in [9], but not always.
Consider a magnet and a vacuum chamber touching the poles of the magnet. The method in [9] calculates the vector potential by integrating the magnetic field on a boundary surface, which is the surface of the vacuum chamber now. Eqns.(5,6) describe a layer of dipoles and monopoles on the surface of the vacuum chamber. We want to cover the entire volume of the chamber with overlapping spheres to approximate the potentials efficiently. The surface of the vacuum chamber will be inside the spheres and also the monopoles and dipoles. This is a problem because the solid harmonics approximation works only if there are no sources inside the spheres. We can use the method in [9] only when there is sufficient space between the vacuum chamber and the iron poles to cover the beam region with the spheres without sources inside.
The solution to this problem we devised is to shift the location of the magnetic monopoles further, toward the exterior of the volume of interest. The shift of the monopoles will change the fields, so the Eqns.(4, 5, 6) cannot be used. The strengths of the monopoles have to be changed such, that the field values on the bounding surface stay the same. The monopoles are sufficient to reproduce the field at the bounding surface. No need for dipoles.
The bounding surface has to be discretized to calculate the strengths of the monopoles. The surface is partitioned into quadrilateral tiles without gaps. Each tile has some number of Gaussian quadrature points. Formulas for placing the quadrature points inside the tiles can be found in [19,20]. There are infinite possibilities to put the monopoles around the volume of interest. If the quadrature points of each tile are translated along the normal vector of the bounding surface, then we have a simple method to place the sources some distance away from the boundary. The placement of the monopoles is illustrated in FIG.2 for a bending dipole.
After the discretization of the bounding surface and the placement of the monopoles we can set up a linear system of n equations to calculate the strengths of n monopoles as: where A is a n by n matrix, x and b are column vectors with n elements. Eq. (22) is solved for x. The elements of the vector b are calculated as the dot product of a unit vector normal to the bounding surface at a quadrature point and the magnetic (or electric) field vector at the same point.
The unit vectorn i is normal to the bounding surface at the ith quadrature point and it points outward. Magnetic field values B i , are given as an input by some magnetic measurement or got from an EM simulation software. The elements a i,j of the matrix A are the dot products of the normal unit vectorn i and the magnetic (or electric) field vector produced by the jth monopole, as: where i, j = 1..n, r i is a vector pointing to the location of the ith quadrature point and r j is a vector pointing to the jth monopole. Eq. (22) is solved by LU decomposition, a standard direct method. Once the strengths of the monopoles are obtained, the vector potential can be calculated anywhere inside the volume of interest and to some extent beyond it by summing the contributions from each monopole. The vector potential of the Dirac monopole with unit charge was already given in Eq. (7).
In case of an electric field, the scalar electric potential is V E (r; r ′ ) = q/(4πǫ 0 r − r ′ ).
One must take care of the singularity in Eq. (7), which is along the half line in the direction of the Dirac string. This direction is a free parameter, but it must be set such that it doesn't intersect the beam region. For more about this topic, see [9] and [21].
The distance of the monopoles from the bounding surface is also a free parameter. It takes a little experimentation to find the best value. If the monopoles are too close, the result will be less smooth. If they are too far, the solution of Eq.(22) will be less accurate. There is another consequence of this parameter, which we found rather important. While the magnetic field on the surface stays the same with different distances, the vector potential will be different. The choice of this distance acts as a gauge transformation. This choice has an influence on the amplitude of the energy oscillation during the tracking.
The placement of monopoles described above and depicted in FIG. 2 works well for all multipoles used in accelerators. In these elements, the longitudinal component of the magnetic or electric field has zero integral around a complete turn of the machine. Let's consider now the field of a solenoid magnet. In this case, the integral of the longitudinal component of the magnetic field is not zero along a complete turn. One can describe the B field of a magnetic monopole not only with the vector potential in Eq.(7) but also with a magnetic scalar potential. Therefore the integral of the longitudinal component of the B field produced by magnetic monopoles along a complete turn of the accelerator must be zero, because ∇ϕ(r) · dr = 0 for any potential ϕ(r). We can conclude that an arrangement of magnetic monopoles alone cannot reproduce the magnetic field of a solenoid magnet. The longitudinal component of the field has to be supplied by electric currents. That can be achieved by adding one or more current loops around the bounding surface as FIG. 3 shows. The geometry of the current loops is not important. However, they should not be too close or too far from the bounding surface. They must give the same integral of the longitudinal component of the B field as the solenoid we want to reproduce. The vector potential is obtained as a sum of the vector potentials of the magnetic monopoles and the current loops.
An arrangement of magnetic monopoles and current loops can describe any static magnetic field. This follows from the Helmholtz decomposition theorem. It states that any smooth rapidly decaying vector field can be described in 3 dimensions as a sum of a curl-free and a divergence-free vector fields. The field to be reproduced is the B field, the divergence-free field corresponds to the curl of the magnetic vector potential of the current loops and the curl-free field can be associated with the gradient of the magnetic scalar potential of the monopoles. We do not use the magnetic scalar potentials of the monopoles, only the vector potentials, but the B field obtained from them is the same apart from the singularity at the Dirac string, which is outside of the beam region.
An arrangement of electric monopoles alone is sufficient to describe any static electric field.
In particle accelerators, an arrangement of monopoles and currents loops depicted on FIG. 2. and FIG. 3. cover most cases. To make the method complete, we consider a simply connected domain with a closed boundary surface around it. In this case, the static magnetic or electric field lines can not have closed loops inside the volume of interest which is free of sources. Therefore a scalar potential is enough to describe the fields. The B or E fields can be reproduced by magnetic or electric monopoles placed outside a closed boundary surface. There is no need for current loops.
Compared to the surface method described in [9], our method is always applicable to calculate the potentials for the solid harmonics expansion, which is not always the case for [9] as we explained earlier. When there is enough space to cover the beam region with spheres without having iron inside the spheres, one is free to choose between the two methods. Ours will give a smoother potential for the same number of quadrature points because the monopoles are located further from the bounding surface. With the method in [9], it is easier to increase the number of quadrature points on the bounding surface at the cost of slower evaluation of ϕ(x i ) during the calculation of coefficients c ℓm in Eq. (14), because there is no need to solve Eq. (22) .   FIG. 4 shows an example for a sector dipole magnet, where the B was reconstructed as the curl of A calculated with the procedure above. In this example, there was no need to use current loops, magnetic monopoles around an open cylinder as boundary surface suffices to obtain a good approximation of the A field.
The accuracy of the reconstructed B field depends on several factors, most importantly the number of monopoles. In this example, we used 4096 monopoles. The relative RMS error between the reference B and ∇ × A is 4 × 10 −3 . The errors were a bit smaller for the quadrupole and for the solenoid magnet we tested. In an actual particle accelerator model, one should use a scaled version of A, such that the integrated curl of A along the reference trajectory matches the magnetic measurement of the accelerator component. After this scaling, a deviation would remain, as illustrated in the magnified part of FIG. 4. The undulation of the blue trace is due to the finite number of monopoles. One can reduce it by increasing the number of monopoles. The amplitude of the undulation also depends on the distance from the reference orbit. The farther the evaluation point is from the reference orbit, the bigger the undulation. However, we haven't found any noticeable effect of this on the beam dynamics during tracking. for a sector dipole magnet halfway between the reference orbit and the maximum aperture. The red curve is the reference field calculated by electromagnetic CAD software, and the blue curve is the curl of the vector potential. On the small sub-picture, the Y-axis is zoomed in by a factor 166.

F. Preparation of the global field map
To test the method, we assembled a small accelerator from eight sector bending magnets. The edge focusing is sufficient to achieve stable optics when the distances between the dipoles are correctly chosen. The workflow to prepare this simple circular accelerator model comprised of the following steps. Each item in the list corresponds to the same task or execution of some command of the software written. Several of these steps have previously been discussed.
1. Obtain the reference orbit through the magnet, including most of the fringe fields. This is usually a tracking result in some EM CAD software such as Opera, CST EM Studio, in the form of a text file containing the coordinates of a reference particle trajectory.
2. Construct a boundary surface around the reference orbit by extruding a profile around it, a circle, ellipse or a rectangle in most cases. Partition this surface into quadrilateral tiles each containing a number of Gaussian quadrature points. Elevate the quadrature points from the boundary surface to some distance. Place monopoles, magnetic ones in this example, to the location of the elevated points and write their coordinates into a text file.
3. Get the field values at those coordinates from the EM simulation software of your choice or some magnetic measurements.
4. Calculate the strengths of the monopoles with the field values mentioned in the previous item by solv-ing Eq. (22). Save the positions of monopoles and their strengths into a text file.
5. Read the file as mentioned above, then translate and rotate the monopoles (and current loops if needed) to the correct location and orientation for each magnet of the accelerator.
6. Set the directions of Dirac strings of the magnetic monopoles such that they don't cross the beam region. That is to avoid problems related to the singularity of the vector potential present in Eq.(7).
7. Cover the entire beam region with overlapping spheres. Then, for each sphere calculate the solid harmonics expansion coefficients with Eq. (14) and write these coefficients along with all supplementary information into a binary file.
After completing this procedure, a binary file is available containing solid harmonics representation of the vector potential for any point inside the beam region. This field map can be used to evaluate the vector potential and the derivatives as it is described in sections II C, II D. Our simple test machine contained only eight instances of the same type of magnetic element. In a real-life scenario, there are many different types of magnets. One would make a global field map for each magnet family driven by the same current circuit. Before simulating with some current setting for the circuits, the spherical harmonics coefficients in each field map must be scaled and summed. This is a rather fast operation. The spherical harmonics expansions of the field maps have to be performed on the same set of spheres covering the beam region, otherwise adding the coefficients together would not make sense. Preparing a global field map this way allows to change current settings for different tracking runs quickly without the recalculation of field maps. The machine in our test contained only magnetic elements, but there is no obstacle to mix static magnetic and electric elements to produce field maps.

G. Symplectic tracking results
Once a global field map is prepared as described above, containing spherical harmonics expansions of the vector and maybe the scalar potentials everywhere in the beam region, symplectic tracking can be performed. The motion of a charged particle calculated by the usual relativistic Hamiltonian where A is the vector potential, V E is the electric scalar potential, p is the canonical momentum and q is the charge of the particle. The equations of motion are six first order ODE's. Integrating these equations in the long term with classical methods like the Runge-Kutta introduces errors, like the drift of energy or non-preservation of phase-space volume. Symplectic integration methods don't suffer from these problems. Until recently, the general belief was that higher order symplectic integrators for general non-separable Hamiltonians have to be implicit, although for specific cases explicit integrators were available [1,22]. An explicit second-order method for general Hamiltonians, symplectic in extended phase space, is described in [23]. We implemented this secondorder explicit symplectic integrator, which has the correct long time behavior. As mentioned earlier, the potentials have discontinuities between the spheres. These appear as fictional sources, and they produce an energy drift during the tracking. This energy drift is a distinct phenomenon from the energy drift caused by a nonsymplectic integration. We performed single particle trackings with several ℓ max (ℓ max is the maximum degree of the solid harmonics expansions) to see the effect of discontinuities. Instead of the energy drift, in FIG. 5 we indicated the relative momentum deviation, more usual in accelerator physics. It can be seen that the drift is decreasing exponentially with the increase of ℓ max . Above a certain value of ℓ max , there is no reduction in the drift. Instead, it increases. This is due to the finite precision of the floating point numbers. With the standard double precision, we got the smallest drift about dp/p = 5 × 10 −8 in 10 5 turns with ℓ max = 42. This is good enough for nearly all simulations in accelerator physics. One can make the drift as small as needed, at the cost of more computation by using higher precision floating numbers.
FIG .6 shows the phase space plot in the horizontal and vertical planes at a position outside the dipoles field. As expected, with some reasonable tunes of the machine, the phase space points are on a thin ellipse. There is no sign of change in the phase space volume. The effect of the energy drift is negligible, and it is not noticeable on the phase space plot. The phase space points of the vertical  plane are on a thicker ellipse, this is the effect of the bigger amplitude motion in the horizontal plane and not due to some artificial phase space volume drift. FIG.7 shows the momentum deviation for a single turn around the machine. The long-term drift is not visible in a single turn. The time scale is too short for that. The highest deviation occurs near the edges of the magnets where the vector potential changes fast. One can reduce this by using smaller time steps or using higher order integrator. The second order integrator we implemented needs four evaluations of the potentials and their derivatives at each time step, which in practice is only three because the first evaluation of the current step and last evaluation of the previous step is the same. Yoshida's method [24] allows the construction of higher order integrators, but these require many more evaluations. For the typical relative momentum deviation required for particle accelerators, the second order method might be faster.

H. Implementation and performance
We implemented the algorithms in this paper and many additional codes for testing different ideas as a Java library named SIMPA. It stands for symplectic integration through monopole arrangements. This code is not yet a finished tool, and it is not ready for release. More work is necessary to turn it into a tracking tool. It is available however from the author on request by e-mail for people interested in experimenting with it. We plan to release it later as open source code if there is interest. Although the ideas in this paper were tested in a particle accelerator, we think they are general enough to be useful in other areas of physics as well.
The performance of the method depends on many factors, and there are several computation intensive tasks mentioned in this paper. Some of these are done only once or a few times as part of the workflow listed in section II F. The most important part from the performance point of view is the evaluation of potentials and their derivatives. As one can see in Eqns. (9,15,16,20,18) the number of operations is proportional to the square of ℓ max , the maximum degree of the solid harmonics expansion. In our test on Intel i7-6700 CPU at 3.40GHz, using a single core with hyper-threading enabled, we got about 58000 evaluations per second of the vector potential and its derivatives with ℓ max = 38. This is about 16 turns per second for a single particle with 2 cm step size in our test machine. The speed can be improved many times by implementing a multi-threaded version of the evaluation.
The spherical harmonics expansion of the entire beam region, which comprises 464 overlapping spheres in our test machine, took 880 seconds using all the eight threads. This task is done only once for each magnet group, connected to the same current circuit. Solving the linear system in Eq. (22) with 4960 variables took 28 seconds. This task has to be done only in the preparation phase of the model, maybe a few times for each magnet type.

III. CONCLUDING REMARKS
We developed a method for symplectic long time tracking of charged particles in static magnetic and electric fields. The main idea is to cover the beam region with overlapping spheres and expand the vector and scalar potentials in terms of solid harmonics for each sphere. The discontinuities between neighboring spheres can be reduced as low as needed to perform long time trackings by increasing the degree of the expansion. Such an expansion of the potentials requires the values of the potentials outside the beam region. We devised and described a surface method in section II E which allows getting the potentials in an analytic way anywhere in the region of interest and to some extent beyond that. This surface method is general. It can describe any combination of static magnetic and electric fields in a source free region.
These two ideas combined allowed us to demonstrate long time symplectic tracking in a circular accelerator with general static fields for the first time.
One should not mislead by the simplicity of the ring used for the demonstration. The method has been already used for modeling a real machine called ELENA [25], taking full advantage of the unique abilities of our approach. For the first time, we managed to model the effect of toroid kicks due to the magnetic system of the electron cooler in long time trackings. The results of these simulations will be the subject of a different paper.
Our method doesn't aim to compete with established tracking codes like MAD-X. It is somewhat complementary to those. Most tracking codes use idealized elements, often they are "kick-codes" and therefore much faster and easier to use. The algorithm described here is useful when those codes can not be used, or their model is not accurate enough due to the idealization of accelera-tor components. This can be the case for small machines with significant fringe fields or machines with special elements. One example of such a machine is ELENA, where the magnetic field of the electron cooler has a considerable influence on the beam dynamics. Despite a significant effort, previous attempts failed to perform long time tracking in ELENA including important effects due to the toroid kicks of the electron cooler [26].
We implemented software to demonstrate and test the algorithm. If there is sufficient interest, our software can be turned into a general purpose library for evaluating static magnetic or electric fields or into a more specific tool for symplectic tracking in accelerators.

ACKNOWLEDGMENTS
I would like to express my gratitude towards Andrea Latina from CERN BE-ABP and Pavel Belochitskii from CERN BE-OP for their valuable comments and encouragements.