Explicit symplectic integrator for particle tracking in s-dependent static electric and magnetic fields with curved reference trajectory

We describe a method for symplectic tracking of charged particles through static electric and magnetic fields. The method can be applied to cases where the fields have a dependence on longitudinal as well as transverse position, and where the reference trajectory may have non-zero curvature. Application of the method requires analytical expressions for the scalar and vector potentials: we show how suitable expressions, in the form of series analogous to multipole expansions, can be constructed from numerical field data, allowing the method to be used in cases where only numerical field data are available.


I. INTRODUCTION
The magnetic and (in some cases) electric fields used to guide particles in an accelerator are often arranged so that particles ideally follow a curved trajectory. In simple cases, for example a magnetic dipole field, standard expressions can be used to calculate the path of a particle through both the main field and the fringe field regions of the relevant element. However, in more complex cases, calculating particle trajectories can be challenging: such cases include, for example, situations where quadrupole or higher-order multipole fields are included by design within a dipole field, or where account needs to be taken of multipole components occuring from systematic or random errors within the element. In general, the problem of particle tracking can be broken down into two parts. First, an accurate description of the field is needed; and second, the equations of motion through the field must be integrated to find the path followed by a given particle. It is often possible to use a numerical field map to describe the field; then, standard integration algorithms (for example, Runge-Kutta algorithms) can be used to integrate the equations of motion. However, an approach such as this can be computationally expensive, both in terms of the memory needed to store the field data, and in terms of the processing involved in integrating the equations of motion. Furthermore, if there are specific constraints or requirements for the trajectories, then additional challenges can occur. For example, if the tracking must obey the symplectic condition, then an explicit Runge-Kutta integration algorithm cannot be used. Symplectic Runge-Kutta algorithms do exist, but are implicit in the sense that each step requires the solution of a set of algebraic equations that can add significantly to the computation time.
Regarding the description of the field, an alternative approach to a numerical field map is to represent the field as a superposition of a number of "modes". Given a set of coefficients, the field can be calculated at any position by summing the functions describing the different modes. This is the approach generally taken for multipole fields, for example, where the horizontal and vertical magnetic field components B x and B y (respectively) are given by: The upper limit of the sum, m max is chosen to provide the accuracy required for the field. The advantages of this approach over a numerical field map are first, that the data describing the field are contained in a relatively small set of coefficients, and second, that the calculation of the field at an arbitrary point does not need interpolation between grid points, which can be an issue in some circumstances for a numerical field map. The field represented by the multipole expansion (1) is independent of the distance along the reference trajectory, and so is appropriate for the main field region within an accelerator element. Depending on the situation being considered, fringe fields may be neglected altogether (as in the "hard edge" approximation), or may be represented using appropriate expressions based, for example, on generalised gradients [1] or formulae representing solutions to Maxwell's equations with appropriate limiting behaviour [2]. A semi-analytical field description such as (1) has a further advantage over a purely numerical description in the context of particle tracking. In some cases, it is possible to construct explicit transfer maps parameterised, for example, in terms of the mode coefficients and element length: the transfer maps then offer the possibility of greater computational efficiency over numerical integration techniques, such as Runge-Kutta algorithms. Furthermore, if the transfer maps are constructed in an appropriate way, then the tracking can satisfy requirements such as symplecticity. An explicit symplectic integrator for general s-dependent static magnetic fields, in systems with a straight reference trajectory, has been presented by Wu, Forest and Robin [3]. Application of the integrator requires the derivatives of the vector potential; it is therefore convenient to have a semi-analytical field description, which allows the derivatives to be expressed in terms of appropriate modes in the same way as the potential itself, thus avoiding the need for taking derivatives numerically.
In elements designed to bend the beam trajectory, it is usually convenient to use a reference trajectory that follows the intended curvature of the path followed by the beam. In such cases, the standard multipole expansion (1) must be modified to give a field that satisfies Maxwell's equations. For completeness, we would like to have a set of modes that can be used to describe threedimensional electric and magnetic fields in a co-ordinate system based on a curved reference trajectory, and an efficient method for integrating the equations of motion for particles moving through these fields. In this paper, we present a suitable set of modes for static electric and magnetic fields, and an explicit symplectic integrator for tracking particles through a given field (i.e. a field represented by a certain set of coefficients). The mode decomposition that we use is based on solutions to Laplace's equation in toroidal co-ordinates; the explicit symplectic integrator is developed following the method of Wu, Forest and Robin [3].

II. DEFINITIONS
We consider a particle of charge q moving (at a relativistic velocity v) through a static electromagnetic field described by a scalar potential Φ and a vector potential A = (A x , A y , A s ). The Hamiltonian for the motion of the particle is [4]: where a particle with the chosen reference momentum P 0 has velocity β 0 c and relativistic factor γ 0 = (1 − β 0 ) − 1 2 , and the scaled vector potential a = (a x , a y , a s ) = qA/P 0 . The independent variable for the system is s, corresponding to distance along a reference trajectory. The reference trajectory follows the arc of a circle (in the plane perpendicular to y) with radius ρ = 1/h. At any point along the reference trajectory, the co-ordinates x and y describe (respectively) the horizontal and vertical position of the particle in a plane perpendicular to the reference trajectory. The longitudinal co-ordinate is defined: where the particle arrives at position s along the reference trajectory at time t (and we can assume that for the reference particle, s = 0 at time t = 0). The momenta conjugate to the co-ordinates x and y are: where γ is the relativistic factor of the particle, m is the mass, and v x and v y are the components of the velocity parallel to the x and y axes. The longitudinal conjugate momentum is: where E = γmc 2 + qΦ is the total energy of the particle. To simplify some of the formulae, we introduce the "scaled" scalar potential φ, defined by:

III. DERIVATION OF THE SYMPLECTIC INTEGRATOR
Our method follows the technique of Wu, Forest and Robin [3]. We first extend phase space by introducing a new independent variable σ, so that s is now a dynamical variable with conjugate momentum p s . The Hamiltonian describing the motion of a particle through an electrostatic field with scaled potential φ = φ(x, y, s) and magnetic field described by a scaled potential a = (a x , a y , a s ) is now: We shall consider the special case where the magnetic field has a uniform vertical field component, which can be represented by a component of the vector potential: where k 0 = qB 0 /P 0 for a magnetic field of strength B 0 . If the field is correctly matched to the curvature of the reference trajectory (so that the reference trajectory is a possible physical trajectory of a particle with momentum P 0 ), then h = k 0 . Other components of the magnetic field can be included in the components a x and a y of the vector potential. We assume that the dynamical variables take small values, so that we can approximate the Hamiltonian by expanding the square root to some order in the dynamical variables. In the conventional paraxial approximation, the expansion is made to second order. Here, we expand to third order, and obtain: where: Viewed as a Hamiltonian in its own right, the term H 1s is integrable, but this is not the case for the other terms, H 1y , H 1x , H 2 or H 3 . However, by making appropriate canonical transformations to new variables, we can express H 1y , H 1x and H 2 in integrable form. H 3 is of order 3 (or higher) in the dynamical variables; we assume we can drop this term (with some loss of accuracy in the solution to the equations of motion). We can then construct an explicit symplectic integrator as follows: where: Continuing the process: e − ∆s 8 :H1s with the transformations of all other variables (not shown explicitly) corresponding to the identity. Now consider H 1y . To find an explicit form for the transformation generated by H 1y , we first consider a transformation to new variables, defined by a mixedvariable generating function: where X i = (X, Y, Z, S) are the new co-ordinates, p i = (p x , p y , δ, p s ) are the original momenta, and I Y is defined by: In Goldstein's nomenclature [5] F y (X i , p i ; σ) is a mixedvariable generating function of the third kind. The new co-ordinates (X, Y, Z, S) are identical to the original coordinates (x, y, z, s), since: and similarly for y, z and s. The new momenta are: and: In terms of the new variables, H 1y can be written: Viewed as a Hamiltonian, H 1y is integrable. The transformations (generated by H 1y ) of the dynamical variables are: Again, the transformations of all other variables (i.e. for those variables not shown explicitly, above) are given by the identity transformation. To apply the transformation e − ∆s 4 :H1y: , we first transform from the original variables to a set of new variables using (24)-(26); we then apply the transformations (29)-(31), and finally transform back to the original variables using the inverse of the transformations (24)-(26). Note that although the new momenta do not change under the transformation generated by H 1y , the change in the Y co-ordinate leads to a change in p x , p y and p s because the inverse of transformations (24)-(26) have to be calculated at a different point from the original transformations. Thus: e − ∆s 4 :H1y: where y 0 and y 1 correspond to the initial and final values of the co-ordinate y under the transformation e − ∆s 4 :H1y: . There is also a change in p s ; but this has no effect on the dynamics. In summary, to apply the transformation e − ∆s 4 :H1y: we need to evaluate a y (at the initial value of the co-ordinate y = y 0 , and at the final value of the coordinate y = y 1 ), and the integral (with respect to y) of the derivative of a y (with respect to x).
In Section IV B we give analytical expressions for the components of the vector potential, based on a threedimensional "multipole" decomposition of a magnetic field in a region with a curved reference trajectory. It is also possible to write down expressions for the derivatives of the vector potential; however, the integral in (32) needs to be performed numerically. Although this will make a significant contribution to the computational cost for each step in the tracking calculation, in most cases the integral should converge reasonably quickly given that the derivative of the potential (which is related to the field strength) should vary slowly over the range of the integral (corresponding to the change in the y co-ordinate over the tracking step).
The transformation with generator H 1x may be handled in a similar way to that generated by H 1y , by first transforming to new variables. For the case of H 1x , we use the mixed-variable generating function: where: Note that the new variables in this case (co-ordinates X, Y , Z and S, and momenta P X , P Y , P Z and P S ) are formally different from the variables in the previous case; but to avoid introducing further notation, we use the same symbols. The transformations (generated by H 1x ) of the dynamical variables are: where: and x 0 and x 1 are the values of x before and after the transformation, respectively. The variables y and δ are unchanged by the transformation. Finally, we find explicit expressions for the transformation with generator H 2 by again first transforming to new variables. In this case, we use a mixed-variable generating function: where X i = (X , Y , Z , S ) are the new co-ordinates, and p i = (p x , p y , δ, p s ) are the original momenta. The new coordinates are identical to the original co-ordinates, since: and similarly for y, z and s. The new momenta are: and: In terms of the new variables, H 2 can be written: which is an integrable Hamiltonian, leading to the trans-formations: e −∆s :H2: Again, transformations of the variables not given explicitly above, are equal to the identity.

IV. s-DEPENDENT FIELDS IN TOROIDAL CO-ORDINATES
Applying the symplectic integrator described in Section III involves derivatives of the scalar potential, and derivatives and integrals of the vector potential. It is therefore helpful to have analytic representations of the scalar and vector potentials, from which expressions for the derivatives and integrals may be found. In practice, however, only a purely numerical representation of the potentials may be available (giving, for example, the values of the potentials on a grid of discrete points over some region of space). With a straight reference trajectory (h = 0), it is possible to fit the coefficients of series representations of the potentials, for example using generalised gradients [1]; the series representation gives the functional dependence of the potential on the coordinates, and this therefore provides a suitable representation for applying the integrator.
A similar approach is possible in the case that the reference trajectory has some non-zero curvature. Expressions for "curvilinear multipoles" (multipole fields around curved reference trajectories) have been given by McMillan and others [7][8][9][10], and have been implemented in the tracking code Bmad [11]. However, the available expressions are not ideal for use where the potential is given in purely numerical form. In much of the previous work, the multipoles are expressed in terms of the transverse Cartesian co-ordinates, x and y: obtaining the multipole coefficients then involves fitting polynomials to the numerical data along either the x or y axis [12]. The nature of the potential (which satisfies Laplace's equation) is such that residuals to the fit will grow exponentially with distance from the line along which the fit is performed. A more robust approach is based on fitting to a surface bounding some region of space enclosing the reference trajectory: within the surface, the residuals decrease exponentially with distance from the surface. Although the residuals will still grow exponentially outside the region enclosed by the surface, if the surface is chosen appropriately then the enclosed region will cover the volume of interest for particle tracking.
To obtain a multipole decomposition based on fitting numerical data on a surface, it is convenient in the case of a curved reference trajectory to work in toroidal coordinates [13,14]. The co-ordinates in the transverse plane are illustrated in Fig. 1. The toroidal co-ordinates u and v are related to the accelerator co-ordinates x and y (Cartesian co-ordinates in a plane perpendicular to the reference trajectory) by: where ρ = 1/h is the radius of curvature of the reference trajectory. The longitudinal co-ordinate s (the distance along the reference trajectory) is related to the toroidal co-ordinate θ by: A surface enclosing the reference trajectory can be defined by specifying a fixed value u ref for the co-ordinate u: a surface defined by u = u ref for 0 ≤ v < 2π and 0 ≤ θ < 2π resembles a torus. If numerical field data are available for the scalar and vector potentials on such a surface, then it is possible to fit the coefficients of series expansions for the scalar and vector potentials (up to some desired order) to the data. This produces expressions that are suitable for use in the explicit symplectic integrator described in Section III. We first discuss the case of the scalar potential, and then extend the results to the vector potential.
A. Scalar potential in toroidal co-ordinates In terms of the toroidal co-ordinates, an harmonic potential (such that ∇ 2 φ = 0) may be written [13,15]: where the f mn are coefficients representing the strength of a multipole component φ mn . The multipole components are given by: where P µ ν (ξ) is an associated Legendre polynomial of the first kind, and: An algorithm for computation of the associated Legendre polynomials with positive µ has been presented by Segura and Gil [16]; values for negative µ are readily obtained using [17]: where Q µ ν (ξ) is an associated Legendre polynomial of the second kind. Note that for integer µ (which is the case of interest here), the term in Q µ ν (ξ) in (58) vanishes. We shall show in Section IV C that each component φ mn has properties that may be expected of a multipole of order m, with m = 1 corresponding to a dipole, m = 2 a quadrupole, and so on. Note that a normal dipole deflects a particle horizontally, whereas a skew dipole deflects a particle vertically.
Given numerical data for a potential φ(u, v, θ), the coefficients f mn may be obtained from: where u ref is a fixed value of u that defines the surface (enclosing the reference trajctory, x = y = 0) on which the fit to the numerical data is performed, and N mn is a normalising factor: As an alternative to calculating the coefficients f mn from the scalar potential, they may be calculated from the electric field components. The electric field is derived from the potential by: The E v component of the field (tangential to a line defined by fixed values of u and θ) is given by: The coefficients f mn can then be found from the values of E v on a surface u = u ref : where: To apply the symplectic integrator described in Section III, we need the derivatives of the potential with respect to the Cartesian co-ordinates. The derivates can be obtained from: For a given multipole component (56) the derivatives with respect to the toroidal co-ordinates u and v are: and: Finally, we need the derivatives of the toroidal coordinates (u, v) with respect to the Cartesian coordinates (x, y). The toroidal co-ordinates can be expressed in terms of the Cartesian co-ordinates as follows: We then find: and: The derivatives of the potential with respect to the Cartesian co-ordinates can be found by using equations (68), (69), (71) and (72) in equations (65) and (66). Tracking a particle through a field described by a scalar potential can then be achieved by using the potential and its derivatives (with respect to x and y) in the symplectic integrator described in Section III.

B. Vector potential in toroidal co-ordinates
To apply the explicit symplectic integrator to a particle moving through a magnetic field, we need expressions for the components of the vector potential. Since we address the case of a curved reference trajectory, we assume that the magnetic field has a (normal) dipole component derived from the longitudinal component a s of the vector potential (8). Other components of the magnetic field (corresponding to quadrupole, or higher-order multipole components) may be derived from the transverse components of the vector potential. In toroidal co-ordinates, these components may be expressed as follows: where the functions φ mn are given by (56). In the case that a θ = 0 (i.e. the longitudinal component of the vector potential is zero, so that k 0 = 0 in (8)), and α mn = f mn for all m, n, it is found that: with φ given by (55). Hence, the magnetic field derived from the vector potential a = (a u , a v , 0) with components (in toroidal co-ordinates) given by (73) and (74) has the same form as the electric field derived from the scalar potential φ given by (55).
To apply the symplectic integrator described in Section III, we require the components of the vector potential in Cartesian co-ordinates, and their derivatives. Given the components (a u , a v ) in toroidal co-ordinates, the components (a x , a y ) in Cartesian co-ordinates are obtained from: where the normalising factor N is: .
The derivatives of a x and a y with respect to the Cartesian co-ordinates x and y can be expressed in terms of the derivatives with respect to the toroidal co-ordinates u and v: Given (73) and (74), the derivatives of a x and a y with respect to the toroidal co-ordinates may be found from the second derivatives of the scalar potential: where: (2u)) cos(v) cosh(u) + 4n(6(n − 1) + (n cosh(2u) + n − 2) cos(2v)) +(5 + 4n(7n − 3)) sinh(u) 2 + (1 + 2n) 2 sinh(u) sinh(3u), (84) C. Examples of multipole potentials in toroidal co-ordinates To illustrate the scalar potential given by (55), we consider the case that the potential is independent of the longitudinal co-ordinate, θ: as a consequence, we need to include only a single longitudinal mode, n = 0 in the summation in (55). With a straight reference trajectory (h = 0), we expect a multipole potential to take the form: where the real and imaginary parts of the coefficient C m determine the strengths of the normal and skew components of the field. Hence, in a normal multipole field of order m the potential varies along the x axis as: and along the y axis as: With a curved reference trajectory, we expect to see similar behaviour in the dependence of the potential for a given order of multipole on the x and y co-ordinates, but with some difference from the dependence given in (89) arising from the curvature. One way to show a similarity between multipoles with straight and curved reference trajectories would be to expand the potential in the case of a multipole with curved reference trajectory as a series in x and y; unfortunately, the fact that the limit x → 0, y → 0 corresponds to u → ∞ makes it problematic to obtain the appropriate series. However, we can plot the potential for a given order of (normal or skew) multipole as a function of x and y: plots for dipoles, quadrupoles and sextupoles are shown in Fig. 2 (normal multipoles) and Fig. 3 (skew multipoles).
From Fig. 2 (top), for example, we see that for a normal dipole the potential has an approximately linear dependence on x. With a straight reference trajectory, we would expect the potential to be independent of y; however, the curvature of the reference trajectory introduces a second-order dependence of the potential on y. In the case of a normal quadrupole (Fig. 2, middle), the potential has a (roughly) quadratic dependence on both x FIG. 2. Scalar potential in normal multipoles with a curved reference trajectory. Each row shows (top to bottom) the potential in a multipole of order n = 1 (dipole), order n = 2 (quadrupole) and order n = 3 (sextupole). The left-hand and middle plots in each row show respectively the potential (black line) as a function of horizontal position x, with y = 0, and as a function of vertical position y, with x = 0. The red lines in the left-hand plots show curves φ ∝ x n . The red lines in the middle plots show curves φ ∝ y n+1 for odd n, and φ ∝ y n for even n. The right-hand plot in each row shows contours of constant potential in the plane perpendicular to the reference trajectory. and y: this again corresponds to the behaviour that we would expect in the case of a straight reference trajectory. Because the curvature of the reference trajectory breaks the symmetry between positive and negative values of x, the effect of the curvature is more evident in the dependence of the potential on x, than in the dependence of the potential on y. For a skew quadrupole (Fig. 3, middle), the potential with a straight reference trajectory is exactly zero along the x and y axes. With a curved reference trajectory, the potential is zero along the x axis (as required by symmetry); but there is a relatively weak fourth-order dependence of the potential on y (with x = 0). Other cases demonstrate the general behaviour we would expect for a multipole potential in a straight co-ordinate system, but with some differences arising from the curvature of the reference trajectory.

V. TEST CASES
To illustrate application of the explicit symplectic integrator presented in Section III, we consider three test cases: a curvilinear magnetic skew sextupole, a curvilinear electrostatic quadrupole, and the fringe field region of an electrostatic quadrupole in the g-2 storage ring [18][19][20][21]. The first two cases are "artificial" in the sense that they are based on fields described by a small number of components; the third case is more realistic, and uses field component coefficients fitted to numerical data obtained from a modelling code. In each case, we track a particle with some chosen initial conditions through the field using the explicit symplectic integrator. For comparison, we also integrate numerically the (Hamiltonian) FIG. 3. Scalar potential in skew multipoles with a curved reference trajectory. Each row shows (top to bottom) the potential in a multipole of order n = 1 (dipole), order n = 2 (quadrupole) and order n = 3 (sextupole). The left-hand and middle plots in the top row (dipole) show respectively the potential as a function of x, with y = 0, and as a function of y, with x = 0 (black line). In the middle and bottom rows (quadrupole and sextupole), the left-hand and middle plots show respectively the potential as a function of x, with y = x tan(π/2n), and as a function of x, with y = −x tan(π/2n) (black lines). The red lines in the left-hand and middle plots show curves φ ∝ x n (or φ ∝ y n in the top row, middle plot). The right-hand plot in each row shows contours of constant potential in the plane perpendicular to the reference trajectory.
equations of motion derived from the exact Hamiltonian (2). All calculations are performed in Mathematica 5.0 [6]; for numerical integration of the equations of motion derived from the Hamiltonian (2), we use the NDSolve function with default settings; although this provides a non-symplectic integration, it should achieve good accuracy.

A. Curvilinear magnetic skew sextupole
As a first illustration of the explicit symplectic integrator presented in Section III we consider the motion of a particle in an electric field with (scaled) magnetic scalar potential given by: The field derived from this potential has the characteristics of a skew sextupole field, as shown in Fig. 4. We choose the field strength such that φ 0 = 5×10 4 , and use a radius of curvature for the reference trajectory ρ = 5 m. A dipole magnetic field is included, represented by the longitudinal component of the vector potential (8), but with k 0 = 1.05/ρ so that there is a slight mismatch between the field and the curvature of the reference trajec- tory.
For the reference particle, we choose β 0 = 0.8, and the initial conditions for the particle to be tracked are: (x, p x , y, p y , z, δ) = (1 mm, 4 × 10 −3 , 1 mm, −0.1 × 10 −3 , 0, 0.02). (93) We track the particle using the explicit symplectic integrator presented in Section III, from s = 0 to s = s max = π 6 ρ, with a step size of ∆σ = s max /10. The integration required in (32) is approximated by Simpson's rule: where the derivative is evaluated in each case at the appropriate (fixed) values of x and s, and at the indicated value of y. A similar approximation is made for the integration in (38). Although these approximations will lead to some symplectic error, this should be small for small step size. In cases where symplecticity is important, more accurate integration routines can be used, though at greater computational cost.
The tracking results are shown in Fig. 5. There is good agreement between the two integration methods.

B. Curvilinear electrostatic quadrupole
As a second illustration of the explicit symplectic integrator presented in Section III we consider the motion of a particle in an electric field with (scaled) scalar potential given by: This represents the potential for a "curvilinear" electrostatic quadrupole, with a strength that varies with longitudinal position along the reference trajectory. The transverse and longitudinal variation of the field are described by m = 2 and n = 12 (respectively) in Eq. (55). The potential is illustrated in Fig. 6. We choose the field strength φ 0 = 200, and use a radius of curvature for the reference trajectory ρ = 5 m. We include a magnetic field, represented by the vector potential (8), but we introduce a small mistmatch between the field and the curvature of the reference trajectory by setting k 0 = 1.05/ρ. For the reference particle, we choose β 0 = 0.8, and the initial conditions for the particle to be tracked are: (x, p x , y, p y , z, δ) = (2 mm, 0, 1 mm, −1.1 × 10 −3 , 0, 0.02). (96) We track the particle using the explicit symplectic integrator presented in Section III, from s = 0 to s = s max = FIG. 5. Results of tracking a particle through a magnetic curvilinear skew sextupole, described by the magnetic scalar potential given in Eq. (92). The black points show the results from the explicit symplectic integrator presented in Section III. The red lines show the results of numerical integration of the equations of motion derived from the Hamiltonian (2).
FIG. 6. Variation of the electrostatic potential (95) in a curvilinear quadrupole, as a function of the co-ordinates s (left-hand plot, for x = 10 mm and y = 0), x (middle plot, for y = 0 and s = π 12 ρ) and y (right-hand plot, for x = 0 and s = π 12 ρ). π 6 ρ, with a step size of ∆σ = s max /40. For comparison, we also integrate numerically the (Hamiltonian) equations of motion derived from the exact Hamiltonian (2). The tracking results are shown in Fig. 7, and again we see good agreement between the two integration methods.
C. g-2 storage ring electrostatic quadrupole As a final example of application of the symplectic integrator, we consider the fringe field regions of the electrostatic quadrupoles in the g-2 storage ring [18][19][20][21]. Values for the potential were calculated (using an FEA code) at points on a uniform Cartesian grid; the values of the potential on a surface defined (in toroidal co-ordinates) by u = u ref = 5.76 were then obtained by (spline) interpolation. On the surface u = u ref , we used 120 grid points in v, with 0 ≤ v < 2π, and 80 grid points in θ, with 0 < θ ≤ 2 • (such that the ends of the quadrupole electrodes are at approximately θ = 1 • ). The reference radius for the co-ordinate system is taken to be the radius of curvature of the reference trajectory in the g-2 storage ring, ρ = 7.112 m: this is the radius of the closed orbit for muons with momentum 3.094 GeV/c. The value of u = 5.76 then corresponds, for v = 0, to a point with x = 0.045 m and y = 0, in the conventional accelerator co-ordinate system, with the origin for the x and y co-ordinates on the reference trajectory.
Based on equation (55), coefficients f mn were calculated so that the potential on any grid point can be found from: (coth(u)) cos(mv) sin(n θ), where n = n 0 (2n + 1), with n 0 = 45 (so that n = 0 corresponds to a sine function with quarter period equal to 2 • , i.e. the range of θ over which values for the potential are given). The values of f mn are obtained essentially by a discrete Fourier transform of the potential on the given grid points. Mode numbers 0 ≤ m ≤ 10 and 0 ≤ n ≤ 79 are used. The truncation in the azimuthal mode number m (compared to the number of data points available) means that the data are not fitted perfectly; however, the contribution of modes (multipoles) of order m > 10 is found to be small. Note that the dominant multipole is the quadrupole component, m = 2. The potential as a function of θ (at u = u ref and v = 0) is shown in Fig. 8, and as a function of v (at θ = 2 • and at θ = 0.25 • , with u = u ref in both cases) in Fig. 9. In the lower plot in Fig. 9, we see that the variation of the potential with the "azimuthal" co-ordinate v in the fringe-field region (about 30 mm from the ends of the electrodes) is significantly distorted from a simple sine wave, indicating the presence of higher-order multipoles.
Using the coefficients f mn we can calculate the potential at any point within the surface on which the fit is performed. As an example, Fig. 10 shows the potential as a function of θ (for v = 0) and as a function of v (for θ = 2 • ). In each plot, the black line shows the potential at u = u ref = 5.76 and the red line shows the potential at u = 6.11: the larger value of u corresponds to a value of x that is a factor of √ 2 smaller than the value of x at u = u ref , so that the potential (for a pure quadrupole) is expected to be smaller by a factor of two. The expected behaviour of the potential (as a function of u) is indeed what we observe.
Tracking a particle through the fringe field of an electrostatic quadrupole using the symplectic integrator described in Section III requires the derivatives of the potential with respect to the accelerator co-ordinates, x, y and s. The derivatives can be calculated (at any point within the surface used to fit the coefficients f mn for the given potential) using equation (97), together with (71) and (72). Some example results from tracking a muon through the fringe field are shown in Fig. 11. The black points in Fig. 11 show the muon trajectory calculated using the symplectic integrator for the detailed fringe-field model, i.e. the model based on the numerical data for the scalar potential. The red line shows the results of an integration using a (non-symplectic) adaptive Runge-Kutta integration of the equations of motion in the same field. The blue line shows the results of a Runge-Kutta integration of the equations of motion through a region with the same magnetic field, but with a "hard-edge" model for the electric field. The hard-edge model is constructed so that the scalar potential is zero up to a point s = s 1 , and is given simply by φ = 1 2 k 1 (x 2 − y 2 ) for s > s 1 . The value of k 1 is chosen to correspond to the focusing potential in the body of the quadrupole found from the numerical data for the scalar potential. The point s 1 is chosen so that the integrated gradient, smax 0 k 1 ds in the hard-edge model is equal to the integrated gradient in the fringe-field model.
There is good agreement between the symplectic integrator and the Runge-Kutta integrator for the detailed fringe-field model. There is little difference between the detailed fringe-field model and the hard-edge model for the horizontal motion, which is dominated by the magnetic field (that is the same in both cases). There is some small but observable difference between the detailed fringe-field model and the hard-edge model for the vertical motion. The change in the vertical momentum after integrating through the full region is approximately the same in both cases: this is expected, since the length of the quadrupole field in the hard-edge model was chosen to give the same integrated focusing strength as the detailed fringe-field model. However, the fact that the change in the vertical momentum occurs at a discrete point in the hard-edge model leads to a slightly larger difference between the models in the vertical co-ordinate at the end of the integration. It is unclear what impact this may have on the beam dynamics in the storage ring, FIG. 7. Results of tracking a particle through the field of a curvilinear electrostatic quadrupole. The potential is given by Eq. (95). The black points show the results from the explicit symplectic integrator described in Section III. The red lines show the results from numerical integration of the equations of motion derived from the Hamiltonian (2).
but it is possible that it may lead to an observable effect over a sufficiently large number of turns. FIG. 10. Scalar potential in an electrostatic quadrupole in the g-2 storage ring. The potential is plotted as a function of toroidal co-ordinate θ at v = 0 (left) and as a function of v at θ = 2 • (right). In each plot, the black line shows the potential at u = u ref = 5.76, and the red line shows the potential at u = 6.11. At the larger value of u, the value of the co-ordinate x is reduced by a factor of √ 2 compared to the value of x at u = u ref ; the potential is a factor of two smaller at the larger value of u, as expected for a quadrupole field.
FIG. 11. Trajectory of a muon through the fringe field region of an electrostatic quadrupole in the g-2 storage ring. The electrostatic potential is shown in Figs. 8 and 9. The reference momentum is 3.094 GeV/c, and the reference trajectory is the arc of a circle with radius 7.112 m, determined by the magnetic field strength, B ≈ 1.45 T. The initial co-ordinates (x, px, ypy, z, δ) of the muon are (10 mm, 5 × 10 −4 , 10 mm, −2 × 10 −6 , 0, −0.02). The black points show the results from the symplectic integrator, with step size 12.4 mm, i.e. a total of 20 steps. The red line shows the results of an integration using a (non-symplectic) adaptive Runge-Kutta integration of the equations of motion in the same field. The blue line shows the results of a Runge-Kutta integration of the equations of motion through a region with the same magnetic field, but with a "hard-edge" model for the electric field.