Symplectic tracking through straight three dimensional fields by a method of generating functions

For simulating single-particle trajectories, the derivation of final coordinates from known initial coordinates through arbitrary electromagnetic fields is of key interest in accelerator physics. We address this task in the case of straight stationary magnetic fields, using generating functions via a perturbative ansatz for the corresponding Hamilton-Jacobi equation. Such an approach is always symplectic, independent of the expansion order. We set up the Hamiltonian by static fields, represented by Fourier series, and outline this approach for the correct and complete set of 3D-multipole fields. Different types of multipoles can be treated with the same formalism, combining them with a specific table of Fourier coefficients characterizing their fields. The resulting particle-tracking routine maps the multipole in a single step. Results are compared with analytical estimations and high-resolution integration methods.


I. INTRODUCTION
Particle tracking is a well-established tool to estimate the dynamic aperture in particle accelerators.Many thousands of turns around a circular machine need to be simulated, and the results depend very sensitively on the correct modeling of nonlinear devices such as insertion devices or fringe fields in end poles of multipoles.Despite the high computing capacity nowadays, efficient tracking algorithms are of great importance in multiparameter scans.Among other methods, the particle tracking based on generating functions (GFs) is extremely efficient in the case of undulators [1].
Having an analytic field description of the device under consideration, analytic GFs can be established [2].These descriptions can be derived for typical magnetic fields by fitting a Fourier series to the magnetic fields, which has been demonstrated for permanent magnet undulator structures [1].The method is extremely useful for the complicated magnetic structure of APPLE II type undulators [3] which produce arbitrary polarization and, thus, have a huge parameter space.The usual kick-map approach can be rather tedious for these devices, because the kick maps have to be established for each gap and each state of polarization individually.Contrarily, the analytic GF approach permits a complete parametrization with respect to an arbitrary gap and state of polarization which simplifies the tracking studies significantly.The energy dependence of the particle path is simply included by scaling the field strength.Other approaches use numerical fitted GFs; see the discussion in Refs.[4][5][6][7].
An analytic GF for particle tracking passing APPLE undulators was implemented into the code ELEGANT which is freely available [8].It was demonstrated that the new method outperforms a symplectic integration tracking scheme by at least one order of magnitude in terms of CPU time [9].
Based on this success, the idea arose to adopt the method of analytic GF-based particle tracking to 3D multipoles.In third-generation storage rings, the multipoles (higher than dipoles) are often treated as 2D devices.With the advent of diffraction-limited light sources with multibend achromats, compact magnet designs with enhanced fringe field effects are getting more and more interesting [10].In future designs, fringe field and cross talk effects will play an important role, and efficient codes for a realistic assessment of the beam performance via particle tracking are needed.
In this paper, we construct and present the use of GFs to set up a mapping routine through arbitrary stationary magnetic fields with a straight reference axis, which includes cases like fringe fields of 3D multipoles, solenoids, combined function magnets, or taper.Hereby we provide the complete set of differential equations to arbitrary orders for the expansion functions, which is not covered in Refs.[1,11], in Sec.II.Magnets with dipole components leading to curved reference orbits are not considered here; they require a more refined calculation.
In our application to multipoles, we apply the method to the cylindrical field descriptions given in Refs.[12,13].We prefer a cylindric geometry for a simplification of the equations; in this case, the coordinates can be separated, and the longitudinal dependency is described as a Fourier series.This is not necessarily the case in Cartesian coordinates: For example, in the analytic Cartesian formulas for fringe fields in Ref. [14], the longitudinal "falloff" functions still involves the transversal coordinates.The drawback of cylinder coordinates is, however, the appearance of a singularity at r ¼ 0, which must be handled carefully.Regardless of this approach, we will see that most of the formulas can also be represented in a Cartesian system.
The 3D-field descriptions are derived by fitting Fourier coefficients with respect to an axis parallel to the longitudinal axis with a slight offset.The Maxwell equations will then determine the entire field, as discussed in Refs.[12,13].This approach can be used to approximate an "Enge-type field falloff" as discussed in Ref. [15]; see also Fig. 2.
For our first demonstration code, it turns out that the step size can be in the decimeter range for reasonable field strengths, as long as the slope of the particle path remains small along the device.In the presented example, a 3D quadrupole is mapped by a single step of around 40 cm.We will also see that the principle of our method is not limited in the order of the slopes by which the Hamiltonian was developed: Higher-order results will depend on higher orders in the transversal derivatives of the fields.All equations for the implementation into a numeric code are given.
In the following, we call our demonstration code STGFM, which means symplectic tracking based on generating functions applied to multipoles.STGFM is not optimized for efficiency; rather, it serves for a thorough comparison in different expansion orders with results from direct integration.In Sec.IV, numeric results of STGFM are presented, and some concluding remarks are given in Sec.V.

A. Hamilton-Jacobi equations
Our starting point will be the well-known Hamiltonian H for time-independent magnetic fields of a linear section.We express H in Cartesian coordinates ðx; y; zÞ, in which the longitudinal coordinate z will play the role of the "time" parameter, x and y are the transversal coordinates, and p x and p y are their corresponding conjugate momenta: Let us denote by x 0 the transverse slope or angle or kick of a moving particle with respect to the x coordinate.It is given by and similarly y 0 ¼ p kin;y =p kin;z .If the total (constant) kinetic momentum p kin ≔ jp kin j is large in comparison to p kin;x and p kin;y , which is the case if the slopes are small, then the square root can be expanded.
In the following, we will assume that the slopes are small enough along the device so that an expansion of H up toand including-third order is sufficient for the physics we are dealing with.Later on, it will become clear, however, that there is no restriction on this order; the method can be developed completely analogous for higher orders.
With the normalized quantities H ≔ H=p kin , px ≔ p x =p kin , and py ≔ p y =p kin , it holds that and similarly for the y components.So H is a Hamiltonian with respect to these quantities.Let us denote the truncated Hamiltonian of H by H, and let us also drop the tilde symbols from now on: Here A x , A y , and A z are corresponding dimensionless time-independent vector potentials; i.e., they emerge from the usual vector potential by multiplication with e=p kin , and the functions A r ¼ cosðφÞA x þ sinðφÞA y and A φ ¼ − sinðφÞA x þ cosðφÞA y are the corresponding normalized potentials in cylinder coordinates.We will switch back and forth between Cartesian and cylinder coordinates whenever it suits best; i.e., we use the transformation rules for cylinder coordinates, which propagate to their derivative counterparts, if necessary: For 3D multipoles, cylinder coordinates are the best choice, because the variables are well separated and partial derivatives are simple.This is important if the code will be generated automatically.On the other hand, some formulas like Eq. ( 6) look more appealing in Cartesian coordinates.
We are looking for a generating function F, which will be a function of r, φ, z and the cyclic momenta v r , v φ , satisfying and The reason is that if such a generating function can be constructed, it will provide a canonical transformation between the original generalized coordinates ðx; y; p x ; p y Þ and a set of comoving (cyclic) coordinates ðu x ; u y ; v x ; v y Þ [16,17].For every particle, these comoving coordinates will, in general, be different; however, they stay constant along its motion.If one can establish the relation between the original coordinates and this new system of cyclic coordinates, one obtains the sought symplectic map from known (original) coordinates to unknown (original) coordinates, by expressing the unknown ones in dependency of the known ones, using said relation.We will show this explicitly in Sec.II B. Inserting (3a) into (2b) yields which is the Hamilton-Jacobi equation we will be dealing with.We introduced a parameter ϵ ∈ ½0; 1 by redefinition of the potential functions.We can tune ϵ globally to express solutions in terms of the order in the field strength in analogy to Refs.[1,11].

B. The recursive differential system
In order to solve (4), we make the following perturbative ansatz: , where F free is the solution (8) in the drift case and f ijk are functions of the parameters r, φ, and z.Let i þ j þ k be the total order of f ijk .Then we see from the Hamilton-Jacobi equation (4) that the products forcing an increase of this order; thus, f ijk;z will be determined recursively by all lower-order f ijk 's and the potential functions.The same argument holds if H emerged from an expansion in higher orders in the slopes.To make this precise in our case, we insert the above ansatz into (4) and perform a coefficient comparison, yielding the following system of differential equations (values with respect to negative indices are treated as zero) [17]: where Z qpk ≔ f qpk;z − δ q0 δ p0 δ k1 A z , X qpk ≔ f qpk;x − δ q0 δ p0 δ k1 A x , and Y qpk ≔ f qpk;y − δ q0 δ p0 δ k1 A y , using Kronecker symbols.The first couple of equations of (6), belonging to the lowest total orders, look similar to Eqs. (29a)-(29k), if one replaces A x with A x − f 001;x and A y with A y − f 001;y together with the additional equation f 001;z ¼ A z .From Eqs. (3a), we obtain with X ij ≔ P ∞ k¼1 f ijk;x ϵ k and Y ij ≔ P ∞ k¼1 f ijk;y ϵ k .Take an interval I ¼ ½I 1 ; I 2 of length L along the z axis and z 0 and z f two positions of interest in the interior of I. Let us understand coordinates labeled without an index (or a zero) as the initial or known coordinates and the ones with an "f" as the final coordinates.Note that in this paper we do not require that z 0 < z f , so the final coordinates may also be some unknown starting coordinates, if the particle travels in the positive z direction.For example, in the case of no fields, the so-called "drift case", one can easily verify that Eq. ( 4) is solved by F free ðr; φ; v r ; v φ ; zÞ where α 0 is a constant, leading to the usual drift solution Since the system of equations ( 6) determines f qpk only up to a function h qpk of r and φ, we can set these functions so that f qpk j z f ≡ 0 holds.Then, we take (7a) and (7b) at the final position to get v r ¼ p xf and v φ ¼ p yf , i.e., the cyclic momenta equal the final Cartesian canonical momenta.
To obtain these momenta p xf and p yf , we then have to consider (7a) and (7b) at the initial position; i.e., we have to compute X ij and Y ij at the initial coordinates, and then invert the equations, for example, by a Newton routine in which the initial momenta p x0 and p y0 are given by the initial slopes x 0 0 and y 0 0 and the corresponding values of the potentials there.We choose a Newton routine, because it will be sufficient for a first demonstration code and is straightforward to be SYMPLECTIC TRACKING THROUGH STRAIGHT THREE … Phys.Rev. Accel.Beams 19, 014001 (2016) implemented.For a first-order expansion in the canonical momenta, Eqs.(7a) and (7b) can be solved explicitly and yield a similar form as Eq. ( 27) in Ref. [1].
Note that since we have by the conditions above, Eq. ( 5) turns out to be in fact an expansion with respect to these final slopes x 0 f and y 0 f and the field strength parameter ϵ.Once the canonical momenta p xf and p yf are determined, we compute the offsets x f and y f using the second set of equations (3b) for the generating function: After the offsets are computed, we obtain the slopes x 0 f and y 0 f at the final position from the previous equations, by inserting the offsets into the potential functions.
In practice, we have to stop the summations at a given index, which means to force higher-order f ijk 's to be zero in contradiction to (6).This error does not affect symplecticity, because the transformation obtained by such a truncated generating function belongs to a correspondingly modified Hamilton-Jacobi equation.The convergence towards the physical values depend on the normalized field strength parameter ϵ and the step size, whereas only the latter one can be freely chosen for real devices.

III. IMPLEMENTATION
A. An iterative solution using Fourier decompositions

Let us assume a Fourier decomposition
of the (transversal) field components within I along the z axis, with k n ≔ 2πn=L and in which b n and bn are complex functions of r and φ.We further assume a gauge in which A z ≡ 0. Then by integration of the normalized vector potentials A r and A φ have the form By Eq. ( 6), higher-order f ijk 's are then computed by successively differentiating lower-order terms with respect to x and y and then integrating them with respect to z.
During that process, the terms with b 0 and b0 introduce more and more mixed expressions involving the quantities R z k expðik n zÞ due to the products in Eq. ( 6).Although this does not constitute an essential obstacle, we have chosen to first expand z into a series in order to keep the resulting expressions simpler.This will introduce certain errors, as we will see, and we also do not recommend it for efficient application.The basic principles, however, can be demonstrated.
Without loss of generality, assume that I is symmetric around z ¼ 0 and consider a Fourier decomposition z ¼ lim K→∞ T K ðzÞ with respect to I, where Then we obtain in which d 001 n and d001 n are complex functions of r and φ, related to the Fourier coefficients b n and bn of the magnetic fields by (n ≠ 0): and d001 0 ¼ 0 ¼ d 001 0 , with e=p kin as the normalization factor.By introducing corresponding Fourier decompositions of f ijk , X ijk , and Y ijk , where c ≔ cosðφÞ and s ≔ sinðφÞ, we obtain in cylinder coordinates [17] where for ðq; p; kÞ ≠ ð0; 0; 1Þ dqpk n e ik n z f else: and with w j ≔ e ik j z − e ik j z f , using the boundary conditions f qpk ðr; φ; z f Þ ≡ 0. In order to maintain a high numeric accuracy, we choose z 0 and z f accordingly so that z f0;fg ¼ T K ðz f0;fg Þ is satisfied; see Fig. 1.Equations ( 16), (17a), and (17b) show how higher total orders of the f ijk 's, and therefore higher-order results in the computations of the coordinates, depend on higher orders in the derivatives of the field coefficients.This goes completely analogous if we would have included higher orders in the slopes.The integration of the Fourier-approximated functions will, however, introduce specific errors: For example, the next integral will have extrema at the zeros of its integrand.Nonetheless, we implemented this method in STGFM in order to verify our ideas, and it turns out that the errors are within reasonable ranges; see Sec.IV.To find precise values (which must be computed only once) for z 0 and z f , we used a Newton iteration beforehand.Equations (7a) and (7b) can be written in vectorial form ðp x ; p y Þ tr ≕ gðp xf ; p yf Þ, which has to be evaluated at the initial coordinates and inverted for ϵ ¼ 1.This was done using another Newton iteration piþ1

B. Multipole fields
Let us take a closer look on the 3D-multipole fields.According to Ref. [12], the analytic formulas for these multipoles can be written as where 2 ≤ m ∈ N is the order of the multipole (m ¼ 2 quadrupole, m ¼ 3 sextupole, etc.) and for p ≥ 0 The function G m;0 ≕ P ∞ n¼−∞ c n e ik n z is determined by the field distribution along the z axis; see below.From (20a) and (20b), it follows for p ≥ 1 that In practice, we have to truncate the infinite summations.
With A z ¼ 0 this yields for n ≠ 0 For a given truncation order P, the c n 's can be determined, for example, by the form of one field component along the z axis, say, B φ ðr 0 ; 0; zÞ (see Fig. 2).Note that if we introduce the Fourier coefficients bn ∈ C by and b0 ¼ 0. This B z component emerged out of a truncated potential using Eq. ( 10) and, thus, does not coincide with the corresponding truncated formula of (19c).In summary, the truncated formulas to (19a)-(19c) do not satisfy ∇ B ¼ 0, whereas the B field emerging from Eqs. (22a), (23b) and (25) do not satisfy ∇ × B ¼ 0. The error eðzÞ due to truncation is on the order of This error cannot be avoided; however, it shrinks with the expansion order P drastically.Even more important, the non-Maxwellian character of the fields does not affect the symplecticity of the whole ansatz.
Note that for P → ∞ we can express the partial derivatives of b n and bn again in terms of b n and bn .One can easily verify that for all n ∈ Z it holds [17] that These formulas can be used to derive general formulas for arbitrary high derivatives using only the original functions [17].We have seen that, for higher orders in the slopes, we require higher orders in the derivatives of the multipole field at the starting position.For small values of P it is, however, recommended to compute the derivatives explicitly.

C. Low-order expansion
In STGFM we implemented a generating function F of order up to 2 in the field strength parameter ϵ and of total order 2 in the momenta p xf and p yf .Recall from Sec. III A that we can write Eqs.(7a) and (7b) in vectorial form.With the given order above, the Newton routine to obtain p xf and p yf thus involves a matrix of the form ½g 0 ða; bÞ Here D is the determinant of g 0 .Since we make the computations up to second order in the field strength parameter ϵ, we have to compute f 101 , f 102 , f 011 , f 012 , f 111 , f 112 , f 201 , f 021 , f 202 , f 022 , and f 002 .The first equations of ( 6), written down explicitly, are Note that in Ref. [1], a subset of these equations are presented.From Eq. ( 16), we obtain The partial derivatives of the η qpk j 's with respect to r and φ are computed from these equations, as well as all higherorder d qpk n 's and dqpk n 's using Eqs.(17a) and (17b).By this method we obtain the entries X ij and Y ij in the matrix above and thus are able to compute the final canonical momenta p xf and p yf .
The computed values have been cross-checked numerically by the model of an extended 2D quadrupole of constant field gradient K.With B ¼ ðKy;Kx;0Þ tr we obtain A x ¼ −KxΔz and A y ¼ KyΔz, where Δz ≔ z f − z.It follows that The values η qpk j have been computed from these equations by inverting (15a): Similarly, the values d qpk j and dqpk j have been checked using To compute the offsets x f and y f (respectively, r f and φ f ), we insert the known final canonical momenta p xf and p yf into Eqs.(9a) and (9b) at ϵ ¼ 1: Note that all higher derivatives of d001

IV. TRACKING RESULTS
The tracking scheme is based on a new method which is not implemented yet in any other tracking code available.Various numeric tests have been performed with respect to other references: The first set of tests compares a full 3D tracking based on the new code with the results derived from numeric integration utilizing the magnetostatic code RADIA [18,19].In a second test run, the new code is compared to analytic expressions which are exact in second-order expansion in the magnetic field for specific initial conditions.In a third step, particle tracking with up to 3 × 10 4 turns has been performed.It is demonstrated that specific resonances are driven only if the fringe fields are switched on.We want to discuss these tests in more detail.
A FORTRAN code has been written for the numeric tests.The various expansions with respect to the canonic momenta and the field strength are programed explicitly which results in a rather bulky code which cannot easily be expanded to higher orders.Furthermore, it is not optimized for high-speed number crunching.On the other hand, the explicit programing is well structured, and, thus, it is suited for all types of numeric checks.The code treats 3D multipoles parametrically; thus, all kind of multipoles discussed in Refs.[12,13] are implemented.The code is freely available [20].As explained in Sec.II, STGFM involves a Fourier decomposition of the function T K within the expansion interval of the transverse fields up to an order K of our choice.
The numeric tests were performed for a quadrupole.In a first step, the Fourier coefficients c n for the field (or vector potential) needs to be determined.For this purpose, a pure permanent magnet Halbach-type quadrupole [21] consisting of 256 radial slices per 360 degrees is evaluated with RADIA.The hard edge length is 200 mm, and the inner and outer radii are 20 and 40 mm, respectively.A single transverse field component is evaluated at three radii of r 0 ¼ 3, 9, and 18 mm over a length of 400 mm (100 mm beyond both ends of the magnet).
Each field distribution is decomposed into an individual set of 25 harmonics a n;r 0 where B φ ðr 0 ; 0; zÞ ¼ P ∞ n¼0 a n;r 0 cosðk n zÞ, with k n ¼ 2πn=L and L ¼ 0.4 m.Equations ( 22b) and (23b), taken at r ¼ r 0 and φ ¼ 0, were then inverted to obtain the c n 's.Because of symmetry considerations, the sinlike terms are set to zero to exclude a possible source of numeric noise in this case.We emphasize that sin terms are implemented into the code as well for the description of nonsymmetric field profiles, as was shown in Sec.II.The Fourier decomposition is well suited for the field description, and even close to the magnet pole tips (radius ¼ 18 mm) the longitudinal field distribution is reproduced with high accuracy (Fig. 2).
In the following, we will use the coefficients as derived from the fields at the reference radius of r 0 ¼ 9 mm.Because of the finite number of segments (256) in the RADIA simulations, the radial field dependence is nonlinear further off axis (Fig. 3).Similarly, in a real magnet, saturation effects of an iron yoke or permeability of the permanent magnets different from one would lead to nonlinearities, which are often approximated by 2D fields, independent of the longitudinal position.Fringing fields are strongly dependent on the longitudinal position and require a 3D-field description.Both types of fields are correctly treated by the applied multipole expansion.Nonlinearities are automatically introduced in the model via the higherorder derivatives of the longitudinal field distribution: P ¼ 1, 2, 3 introduces third-, fifth-, and seventh-order dependencies in r.The third-order terms are also called pseudooctupoles.P can freely be chosen in the FORTRAN implementation.

A. Comparison STGFM versus RADIA
In the first set of tests, we tracked 63 particles with different initial coordinates in a single step of 400 mm through the complete quadrupole including the fringe fields.All 25 Fourier coefficients are included.Different choices for the initial displacement and angle distributions were chosen: (i) Initial displacements and angles are equally distributed either on a horizontal or on a vertical phase space ellipse.The two phase spaces are still uncoupled behind the quadrupole (Fig. 4).(ii) The horizontal and the vertical initial displacements are equally distributed on an ellipse whereas the initial angles are set to zero (Fig. 6).(iii) The horizontal and the vertical initial angles are equally distributed on an ellipse whereas the initial displacements are set to zero (Fig. 7).
In the last two cases the coupling terms (i.e., terms which contain products of the coordinates of both planes) are tested.The results are compared to results from RADIA which uses a Runge-Kutta method.For a better visibility, the drift related to the initial angle and the linear focusing term have been subtracted from the results.The linear focusing term is evaluated as Figures 4-7 give the results for two orders of expansion: The red solid curves belong to orders taken as in Sec.III C, whereas the red dotted curves belong to total order one in the canonical momenta p xf and p yf and order two in the field strength parameter ϵ.The blue curves give the corresponding differences to RADIA (black, solid line).For a better visibility of these differences, a scaling factor of 50 has been used for the solid lines and a factor of 5 for the dotted lines.For the red solid curves, the accuracy is on the percent level in all cases.Higher accuracy is expected for higher expansion orders or for a coordinate transformation in several steps through the magnet.

B. Comparison STGFM versus analytic formulas
In a second test, we checked the code against analytic expressions where the level of accuracy is known.The analytic formulas assume an infinitely long quadrupole, characterized by only the lowest-order Fourier coefficient c 0 in STGFM.The analytic formulas come from the wellknown linear transfer matrix, expanded to second order in field strength.Note that, although both the analytic formulas and the code constitute comparable expansions in field strength, they do not necessarily agree: The Newton routine, for example, introduces higher products of field strength because of the division by the determinant.In our test, c 0 was derived from an 800 mm long quadrupole modeled with RADIA.For an initial displacement or angle in FIG. 4. Final coordinates after tracking through a quadrupole magnet as evaluated with RADIA (black, solid line) and with STGFM for different expansion orders (see the text).All particles are initially equally distributed on the rim of a horizontal phase space ellipse with displacement.Top: ðx 0 ;x 0 0 ;y 0 ;y 0 0 Þ ¼ ð9 mm;10 mrad;0;0Þ.Bottom: ðx 0 ; x 0 0 ;y 0 ;y 0 0 Þ ¼ ð0;0;9 mm;10 mradÞ.
SYMPLECTIC TRACKING THROUGH STRAIGHT THREE … Phys.Rev. Accel.Beams 19, 014001 ( FIG. 5.The same initial conditions, color code, and line type as Fig. 4, this time computed for two steps.We see that, in comparison to Fig. 4, the errors in case of the loworder expansion (dashed blue curves) are significantly reduced.The errors with more expansion orders (solid blue curves) improve only slightly.one plane (the other three initial coordinates are zero), the final displacement and angle in the same plane can be derived analytically.We summarize the results up to second-order expansion: For the starting coordinates x 0 ¼ 0, x 0 0 ≠ 0, y 0 ¼ 0, and y 0 0 ¼ 0, the final coordinates are evaluated from For starting coordinates x 0 ≠ 0, x 0 0 ¼ 0, y 0 ¼ 0, and y 0 0 ¼ 0, we get where K ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c 0 =ðBρr 0 Þ p .Similar equations hold for the vertical phase space.
The analytic expressions can conversely be obtained from the theory section when inserting (31a)-(31i) and the corresponding expressions for the canonical momenta into (34a) and (34b) together with the expansion of the inverse of the determinant (28) for the lower orders.The simulations have been performed for different electron energies to check the convergence of the series (Fig. 8).The following terms y are plotted: where val STGFM and val analytic model are the corresponding final coordinates evaluated with the two methods.At low energies, y scales with an order equal to or higher than three.This is explained by the fact that the analytic expansions are truncated at the second order.At first glance, a third-order scaling is expected to remain for higher electron energies.Interestingly, the scaling is lower than third order.Since the inverse of the determinant is involved in our expressions, we have especially first-and second-order terms in the energy in any expansion order.The local minima in several plots (e.g., lower left in Fig. 8) are due to a sign change of the ratio in the expression of y.Thus, the convergence behavior with energy is well understood.

C. Tracking with FODO lattice
In a third test, we examined the tracking suitability of the new code and compared it against a conventional integration routine WAVE [22,23].For multiturn tracking, simple FODO-like lattices have been set up besides the extended 3D quadrupole.
In the first case, one defocusing kick quadrupole at either side of the 3D-focusing quadrupole is introduced.Second, we added one more defocusing and focusing kick quadrupole on both sides (Fig. 9).The integrated strength of the individual kick quadrupoles equals half the strength of the 3D quadrupole, and the sign is reversed.The simple lattice is closed with a drift of 200 mm between the kick quadrupoles.The observation plane is placed in the center between the kick quadrupoles.We tracked 10 3 turns and analyzed the displacements at the observation plane.The strength of the 3D quadrupole and the two kick quadrupoles can be scaled by a common factor afact.
We started from the model 3D quadrupole depicted in Fig. 2 including 12 Fourier coefficients.This quadrupole has an integrated strength of 2.12=m (normalized to the particle energy).In order to move the tune close to a quarter resonance, all Fourier coefficients and also the strength of the two kick quads are multiplied by the same factor afact ¼ 2.27.The orbit displacement versus turn number in Relative errors between the results of STGFM and those from analytic equations which contain only a linear and a quadratic term in the field and energy (for details, see the text).Blue line, displacement; red line, angle.The data are given for different starting coordinates.Top left: ðx 0 ; x 0 0 ; y 0 ; y 0 0 Þ ¼ ð0; 10 mrad; 0; 0Þ; top right: (0, 0, 0, 10 mrad); bottom left: (9 mm, 0, 0, 0); bottom right: (0, 0, 9 mm, 0).FIG. 9. Top: Simple lattice with one extended focusing quadrupole (EFQ) treated with STGFM and two defocusing kick quads (DFQ) which together have the same absolute strength as the extended quad.Bottom: Extended quadrupole with additional focusing and defocusing quads.The same scaling factor is applied to all elements for adjusting the tune close to 0.25.the symmetry plane is Fourier transformed to determine the tune and the dynamic behavior of the particles (Fig. 10).The positions of the additional lines for P ¼ 1 are of no further relevance in this example; they depend on the initial particle coordinates and distance to the resonance.Switching off the nonlinearity by setting P ¼ 0, a single pronounced tune is observed, describing the oscillation in a parabolic potential.For P ¼ 1 an additional third-order radial field dependency in the transversal components is introduced.
In the next step, the tune is shifted close to the quarter resonance of 0.25 in order to drive the nonlinearity, using a scaling factor of afact ¼ 2.370.Now, the orbit plotted versus the number of turns follows a long period beating (not shown), the Fourier transform of the position variables shows additional lines, and the appearance of a fourfold island structure is shown in Figs.11 and 12 in red.For comparison, results of the WAVE integration routine are shown in black.For the tracking, the subroutine related to the GF was replaced by the WAVE integration routine, applying the same B-field routines and lattice setup.The results show the island structure symmetrically placed in the x − x 0 plane, as expected.The discrepancies to STGFM indicated the limited accuracy of the approximation used.The scaling factor for the quadrupoles was chosen to be afact ¼ 2.370, which corresponds to an integrated quadrupole strength of 5.02=m.When the tune is close to 0.25, the fringe fields as described in STGFM drive a fourth-order resonance as demonstrated in the plot.Switching off the fringe fields (i.e., taking only c 0 ), the resonances disappears (not plotted).The results of STGFM are given in red, and results of WAVE are shown in black.To achieve good agreement, the focal strength in the case of WAVE was reduced by 7%.

SYMPLECTIC TRACKING THROUGH STRAIGHT THREE …
Phys.Rev. Accel.Beams 19, 014001 (2016) 014001-13 The tilt of the pattern in Fig. 11 is due to the finite expansion order, which comes into effect at high fields: The integrated 3D-quadrupole strength is 5.02=m, while the strongest BESSY II quadrupole has 1.3=m.The rotation nearly vanished when the 3D-quadrupole strength is reduced by roughly a factor of 2 as depicted in Fig. 12 (strength 2.58=m), indicating a much better agreement of STGFM.Here, the extended lattice of Fig. 9 (bottom) has been applied.All multipoles are scaled with the same factor such that the tune is close to 0.25, as before.It must be emphasized that all tracking results are evaluated with a single step of about 0.4 m along the 3D quadrupole.
To compare the integration routine with results of STGFM, the field strength factor afact was reduced by about 1% for the integration routine to achieve a better agreement of the island size and position of both routines.

V. CONCLUSION
We presented a GF-based tracking scheme through a complete set of 3D multipoles.For test purposes we implemented the algorithm into a FORTRAN code.Various tests with this code have been performed to confirm the validity of the algorithm.In this paper, a quadrupole magnet has been studied as a test case.The method is applicable to any kind of stationary magnetic fields where the reference orbit is described by a straight line.
Today, the implementation of the code is limited in the expansion (total) order which is "hard wired," and it is not optimized for fast computing.In the future, we plan to examine generalizations to comoving coordinates in storage rings and time-dependent fields.Furthermore, CPU tests with an improved code will be performed for a comparison with other codes and step-size tests.This may also involve a review of the effects using different gauges.
It is expected that the new code will be extremely useful for sophisticated magnet designs such as combined function magnets or compact magnet designs with significant cross talk as well as long periodic structures with fringing fields.These sophisticated magnet designs are already now part of the upcoming diffraction-limited light sources [24].

ACKNOWLEDGMENTS
The first author thanks the sponsoring by the Wolfgang Gentner Program of the Federal Ministry of Education and Research (BMBF) and Atoosa Meseck of Helmholtz-Zentrum Berlin (HZB) for fruitful discussions.We also thank Michael Scheer for providing us with a WAVE subroutine which does the tracking part of his program.12. Tracking results in the horizontal plane for 10 3 turns with the lattice as depicted in Fig. 9 (bottom).Starting coordinates are x 0 ¼ 5, 10.0, 15.0, 20.0, 25.0, 30, and 35 mm, while x 0 0 ¼ 0 mrad in all cases.The scaling factor for the quadrupoles was chosen to be afact ¼ 1.218, which corresponds to an integrated quadrupole strength of 2.58=m.The rotation of the pattern is greatly reduced as compared to Fig. 11, since the quadrupole field strength is about 2 times smaller compared to Fig. 11.The results of STGFM are given in red, and results of WAVE are shown in black.As before, the focal strength was reduced by 1% in the case of WAVE to achieve good agreement of both methods.

5 FIG. 1 .
FIG.1.Enlarged part of the Fourier decomposition of the function T K ðzÞ − z with K ¼ 25 (red curve), 15 (blue curve), and 5 (magenta curve) Fourier coefficients.The function is expanded over the same interval where the Fourier decomposition of the reference field from RADIA is done (Sec.IV).The integration is conducted from one zero crossing to another one, where the zero crossings are determined with the Newton method.The integration interval starts at z 0;25 , z 0;15 , and z 0;5 depending upon the number of Fourier coefficients utilized (25, 15, and 5, respectively). b

FIG. 2 .
FIG.2.Thick lines: Longitudinal B φ -field distribution in the midplane of a quadrupole at horizontal off-axis distances of r 0 ¼ 3, 9, and 18 mm; thin lines: enhanced residuals (same color coding) between the RADIA field and a Fourier representation with 25 coefficients.

n and d 001 n
can be expressed in terms of b n and bn by (27a)-(27d).

FIG. 3 .
FIG.3.Distribution of the field component B φ in the midplane of a permanent magnet quadrupole as described in the text.The residuals represent the terms with nonlinear radial dependency in our test model.

FIG. 6 .
FIG.6.The same color code and line type as Fig.4, but with different initial coordinates (see the text): starting with particles equally distributed on a circle with displacements x 0 ¼ 9 mm and y 0 ¼ 9 mm (bottom: x 0 0 ¼ 0 mrad and y 0 0 ¼ 0 mrad).

1 FIG. 10 .FIG. 11 .
FIG.10.Logarithmic plot of the Fourier transform of the orbit excursion for 10 3 turns.The maximum of the peaks corresponds to the tunes.P ¼ 0 corresponds to purely linear fields; P ¼ 1 introduces a third-order field contribution.Parameters: 1.7 GeV, afact ¼ 2.27, start amplitudes x 0 ¼ 40 mm and y 0 ¼ 0 mm.The tracking through the thick element was done in a single step with 12 Fourier We use the same ordering as described at the beginning of Sec.III C.
FIG.12.Tracking results in the horizontal plane for 10 3 turns with the lattice as depicted in Fig.9(bottom).Starting coordinates are x 0 ¼ 5, 10.0, 15.0, 20.0, 25.0, 30, and 35 mm, while x 0 0 ¼ 0 mrad in all cases.The scaling factor for the quadrupoles was chosen to be afact ¼ 1.218, which corresponds to an integrated quadrupole strength of 2.58=m.The rotation of the pattern is greatly reduced as compared to Fig.11, since the quadrupole field strength is about 2 times smaller compared to Fig.11.The results of STGFM are given in red, and results of WAVE are shown in black.As before, the focal strength was reduced by 1% in the case of WAVE to achieve good agreement of both methods.