Accurate Transfer Maps for Realistic Beamline Elements: Part I, Straight Elements

The behavior of orbits in charged-particle beam transport systems, including both linear and circular accelerators as well as final focus sections and spectrometers, can depend sensitively on nonlinear fringe-field and high-order-multipole effects in the various beam-line elements. The inclusion of these effects requires a detailed and realistic model of the interior and fringe fields, including their high spatial derivatives. A collection of surface fitting methods has been developed for extracting this information accurately from 3-dimensional field data on a grid, as provided by various 3-dimensional finite-element field codes. Based on these realistic field models, Lie or other methods may be used to compute accurate design orbits and accurate transfer maps about these orbits. Part I of this work presents a treatment of straight-axis magnetic elements, while Part II will treat bending dipoles with large sagitta. An exactly-soluble but numerically challenging model field is used to provide a rigorous collection of performance benchmarks.


I. INTRODUCTION
For the design of high-performance linear accelerators, synchrotrons, and storage and damping rings it is essential to have realistic electric and magnetic field information for the various beam-line elements in order to compute accurate design orbits and accurate high-order transfer maps about the design orbits. There are similar requirements for other charged-particle beam transport systems such as final focus sections, high-resolution spectrometers, and high-resolution electron microscopes.
Realistic field data can be provided on a grid with the aid of various 3-dimensional finite element codes, sometimes spot checked against measured data. But the computation of high-order transfer maps based on this data appears to pose an insurmountable problem: the calculation of high-order transfer maps requires a knowledge of high derivatives of the field data. The direct calculation of high derivatives based only on grid data is intolerably sensitive to noise (due to truncation or round-off) in the grid data [1]. We will see that this problem can be solved by the use of surface methods. The effect of numerical noise can be overcome by fitting field data on a bounding surface far from the beam axis and continuing inward using the Maxwell equations. (We recall that interior fields are uniquely specified by their values on bounding surfaces.) While the process of differentiation serves to amplify the effect of numerical noise, the process of continuing inward using the Maxwell equations is smoothing. This smoothing is related to the fact that harmonic functions take their extrema on boundaries. When using surface methods, all fits are made on such boundaries. Therefore, if these fits are accurate, interior data * Work supported by U.S. Department of Energy Grant DE-FG02-96ER40949.
† Electronic address: cemitch@umd.edu ‡ Electronic address: dragt@umd.edu based on these fits will be even more accurate. In this paper we will devote our attention to magnetic beam-line elements. Static electric beam-line elements, and static electric and magnetic beam-line elements such as velocity selectors (Wien filters), could be treated in a similar way. For a treatment of RF cavities, see [2].
There are two magnetic cases that it is convenient to handle separately: straight and curved. For straight magnetic elements such as solenoids, quadrupoles, sextupoles, octupoles, etc., and wigglers, it is convenient to employ cylindrical surfaces. For the case of curved magnetic elements, such as dipoles with large design-orbit sagitta, it is convenient to employ the surface of a bent box with straight ends. In all cases the bounding surface will surround the design orbit within the beam-line element and will extend into the fringe-field regions outside the beam-line element, thus taking into account all fringe-field effects as well as all effects within the body of the beam-line element.
For the case of straight beam-line elements it is convenient to describe the magnetic field in terms of a magnetic scalar potential ψ. Then, if one wishes to compute transfer maps in terms of canonical coordinates, one can proceed with the aid of an associated vector potential A computed from ψ. Alternatively, if one wishes to integrate noncanonical equations employing the magnetic field B, it can be obtained from the relation B = ∇ψ.
For the case of curved beam-line elements it is convenient to work directly with the vector potential. Its use in the case of canonical coordinates is then immediate. If instead one wishes to integrate noncanonical equations employing the magnetic field B, it can be obtained from the relation B = ∇ × A.
In this paper we will treat the case of straight beamline elements. For this purpose we will employ the surface of a (virtual) cylinder of uniform cross-section (circular, elliptical, or rectangular) surrounding the beam, and which lies within all iron or other magnetic sources. In these cases, a Green function can be determined for the geometry of the fitted domain as a series composed of known, orthogonal, and complete functions. In the case of bent-bore magnetic elements, such as dipoles with large sagitta, it is not possible to obtain a suitable Green function. New techniques have therefore been developed for treating cases of such general geometries, and these techniques will be presented in a subsequent paper.
Our emphasis will be on the accurate representation of fields in terms of scalar and vector potentials. Once these representations are available, there are a variety of methods for computing the associated design orbits and transfer maps about these orbits. An Appendix summarizes how this can be done when canonical coordinates are employed and the associated Lie-algebraic structure is exploited.
A brief outline of this paper is as follows: Section II reviews circular cylinder harmonic expansions and introduces a collection of functions, known as on-axis gradients, which uniquely characterize the magnetic scalar potential. It also describes how the vector potential can be obtained from the scalar potential. The on-axis gradients themselves are generally unspecified functions of z. In some simple cases they can be found analytically. However, in general they must be determined numerically. Section III describes how this can be done using known magnetic field values determined at points on some regular 3-dimensional grid. In Section IV, we treat an analytically soluble model problem, that of a magnetic monopole doublet, which is used to benchmark these methods. In Section V, we discuss the advantages of these surface-fitting methods that result from numerical smoothing. In Section VI, we use these methods to study a proposed ILC damping ring wiggler. The paper concludes with a summary and an Appendix. Further detail may be found in [3,4].

A. Scalar Potential
In a current-free region the magnetic field B is curl free, and can therefore be described most simply in terms of a magnetic scalar potential ψ. Because B is also divergence free, ψ must obey the Laplace equation Since we wish to describe straight beam-line elements, it is convenient to work initially in circular cylindrical coordinates ρ, φ, and z with We note for future use that (2) can be written in the form from which it follows that, for non-negative integers l and m, We see that even powers of ρ and the combinations ρ m cos mφ and ρ m sin mφ are analytic (in fact, polynomial) functions of x and y.
A general solution ψ satisfying the Laplace equation and analytic near the axis ρ = 0 takes the form Note that the term containing G 0,s (k) does not contribute to the above sum, and we may set G 0,s = 0 without loss of generality. By utilizing the Taylor series of the modified Bessel function I m , we may write ψ in the form of a circular cylinder harmonic (multipole) expansion: The functions C [n] m,α , known as on-axis gradients [3,4], are defined by We will require an expansion of (6) in the transverse variables x and y. Define the polynomials for integer l ≥ 0 and m ≥ 0. In Cartesian coordinates we then have: Note that each polynomial F l,m α is homogeneous of degree 2l + m.

B. Vector Potential
To determine an associated vector potential A, we must find a solution to the coupled system of equations ∇ × A = ∇ψ, where ψ is given by the series (11). We also must select some particular gauge. And even after a particular type of gauge has been chosen, say a Coulomb gauge, there is still some remaining freedom [4]. One convenient Coulomb gauge choice is given by the rules m,s (z) and C

[n]
m,c (z) describe normal and skew components, respectively. For example, in the body of a long normal dipole (m = 1), we expect C 1,s (z) will be nearly constant (independent of z) and therefore C [n] 1,s (z) 0 for n > 0. Correspondingly, in the body of a long normal dipole, use of (12) gives the results A x A y 0 and Similarly, in the body of a long normal quadrupole (m = 2), use of (12) gives the results A x A y 0 and For any set of analytic functions C [n] m,α (z) employed in (11), the vector potential defined by (12) satisfies the conditions: 1) ∇ × A = ∇ψ = B, 2) ∇ × ∇ × A = 0, and 3) ∇·A = 0. Note that both Maxwell's equations and the Coulomb gauge condition are satisfied by construction. In the following section, we will see how the coefficient functions C m,α (z) can be numerically determined.

III. SURFACE METHODS
There are cases in which Taylor expansions of the form (11) and (12) can be found analytically. In general, however, we have available only measured or numerical threedimensional magnetic field data on a discrete mesh of points distributed throughout the region of interest. The required high derivatives of ψ or A cannot be reliably computed directly from this data by numerical differentiation due to numerical noise, whose effect becomes worse with the order of derivative desired. The effect of numerical noise, and its amplification by numerical differentiation, can be overcome by fitting on a bounding surface far from the axis and interpolating inward using the Maxwell equations. Surface fitting methods have several advantages, including: 1. Only functions with known (othonormal) completeness properties and known (optimal) convergence properties are employed.
2. The Maxwell equations are exactly satisfied.
3. The results are manifestly analytic in all variables.
4. The error is globally controlled. Fields that are harmonic (fields that satisfy the Laplace equation) take their extrema on boundaries. Both the exact and computed fields are harmonic. Therefore their difference, the error field, is also harmonic, and must take its extrema on the boundary. But this is precisely where a controlled fit is made. Thus, the error on the boundary is controlled, and the interior error must be even smaller.

5.
Because harmonic fields take their extrema on boundaries, interior values inferred from surface data are relatively insensitive to errors/noise in the surface data. Put another way, the inverse Laplacian (Laplace Green function), which relates interior data to surface data, is smoothing. It is this smoothing that we seek to exploit. In general, the sensitivity to noise in the data decreases rapidly (as some high inverse power of distance) with increasing distance from the surface, and this property improves the accuracy of the high-order interior derivatives needed to compute high-order transfer maps.
Let us briefly compare this approach to that of on-axis or midplane fitting. In the case of on-axis fitting, it is common to use various analytic model profiles, such as Enge functions, and then differentiate them repeatedly to achieve objective 2 by continuing outward using the Maxwell equations. However, these functions do not have completeness properties, item 1. And there is no smoothing, item 5, to overcome the amplification of noise due to numerical differentiation.
In the case of midplane fitting, one approach would be to attempt to employ expressions that relate the onaxis gradients to midplane data and its derivatives. For example, in the case of midplane symmetry, there are the relations By repeatedly differentiating these relations with respect to z, one can obtain the C

[n]
m,s (z) for n > 0. In general, the determination of C m,s (z) requires the computation of m + n − 1 derivatives. Although this approach achieves objective 2, since all relevant quantities are subsequently computed in terms of on-axis gradients, in the case where data is available only at grid points it presupposes the ability to compute very high-order derivatives by high-order numerical differentiation. This is generally impossible, due to the high noise sensitivity associated with high-order numerical differentiation, because there is no intrinsic smoothing, item 5.
Another approach is to use an analytic functional form, with free parameters, that is known to satisfy the the 3dimensional Laplace equation for all parameter values. These parameters can be adjusted so that the field derived from this representation well approximates the field at various grid points. (These grid points could be in the midplane, but could be out of the midplane as well.) This representation can then be repeatedly differentiated to provide the required field derivatives. However, commonly this fitting procedure has no known completeness/convergence property, item 1. In some cases Fourier series representations with known completeness properties are used. But with Fourier series representations, an artificial periodicity is imposed in the transverse horizontal direction. As a result, the Fourier coefficients for the field expansions, call them a n , can fall off at best as (1/n 2 ) [4]. Correspondingly, the Fourier coeficients in the associated expansion for ψ fall off at best as (1/n 3 ). As a result, repeated differentiation produces nonconvergent series, and there is no analyticity, item 3. Whatever representation is used, there is again no intrinsic smoothing to overcome the amplification of noise due to numerical differentiation.

A. Use of Field Data on Surface of Circular Cylinder
All three-dimensional electromagnetic codes calculate all three components of the field on some threedimensional grid. Also, such data is in principle available from actual field measurements. In this subsection we will describe how to compute the on-axis gradients from such field data using the surface of a circular cylinder [5]. Once these gradients are known, we may use (11) and (12) to compute the associated scalar and vector potentials.
Consider a circular cylinder of radius R, centered on the z-axis, fitting within the bore of the beam-line element in question, and extending beyond the fringefield regions at the ends of the beam-line element. The beam-line element could be any straight element such as a solenoid, quadrupole, sextupole, octupole, etc., or it could be wiggler with little or no net bending. See Fig.  1, which illustrates the case of a wiggler. Suppose the magnetic field B(x, y, z) = B(ρ, φ, z) is given on a grid, and this data is interpolated onto the surface of the cylinder using values at the grid points near the surface. Next, from the values on the surface, compute B ρ (R, φ, z), the component of B(ρ, φ, z) normal to the surface. The major remaining task is to compute the on-axis gradients from a knowledge of B ρ (R, φ, z). See Fig. 2. At this point we note that the functions exp(ikz) sin(mφ) and exp(ikz) cos(mφ) form a complete set over the surface of the circular cylinder.
LetB ρ (R, φ, k) be the Fourier transform of B ρ (R, φ, z) given by the integral Because B decays rapidly in the fringe field region, Calculation of realistic design orbit z d and its associated realistic transfer map M based on data provided on a 3-dimensional grid for a real beam-line element. Only a few points on the 3-dimensional grid are shown. In this illustration, data from the 3-dimensional grid is interpolated onto the surface of a cylinder with circular cross section, and this surface data is then processed to compute the design orbit (trajectory) and the associated transfer map about the design orbit. The use of other surfaces is also possible, and may offer various advantages. therefore its Fourier transform is well defined. Next define the functionsb m,s andb m,c bỹ for m ≥ 1 and from which it follows, using the representation (5), that Now substitute (19) into the right sides of (17) and perform the indicated integrations to get the results from which it follows that This relation for G m,α (k) can be employed in (8) to give the result We have found an expression for the on-axis gradients in terms of field data (normal component) on the surface of the cylinder. Equation (22) may be viewed as the convolution of Fourier surface datab m,α (R, k) with the inverse Laplacian kernel k n+m−1 /I m (kR). Moreover, this kernel has a very desirable property. The Bessel functions I m (kR) have the asymptotic behavior Since I m (kR) appears in the denominator of (22), we see that the integrand is exponentially damped for large |k|. Now suppose there is uncorrelated point-to-point noise in the surface data. Such noise will result in anomalously large |k| contributions to theb m,α (R, k). But, because of the exponential damping arising from I m (kR) in the denominator, the effect of this noise is effectively filtered out. Moreover, this filtering action is improved by making R as large as possible. This filtering, or smoothing, feature will be discussed in more detail in Section V. We close this subsection with the remark that if one wishes to extract the C [n] 0,c (z) (monopole) on-axis gradients from field data, as is required for example in the case of a solenoid, it may be preferable to use the longitudinal component B z (R, φ, z) on the surface of the cylinder rather than the normal component B ρ (R, φ, z) [4].
B. Use of Field Data on Surface of Elliptical Cylinder

Background
In the previous subsection we employed a cylinder with circular cross section, and observed mathematically that it is desirable for error insensitivity to use a cylinder with a large radius R. Physically, this is because we want the field data points employed to be as far from the axis as possible since the effect of inhomogeneities (noise) in the data decays with distance from the inhomogeneity. Evidently the use of a large circular cylinder is optimal for beam-line elements with a circular bore. However, for dipoles or wigglers with small gaps and wide pole faces, use of a cylinder with elliptical cross section should give improved error insensitivity. See Fig. 3. In this subsection we will set up the machinery required for the use of elliptical cylinders, and apply it to the calculation of on-axis gradients based on field data [3,4].

Elliptic Coordinates
Elliptic coordinates in the x, y plane are described by the relations [6] x = f cosh(u) cos(v), Contours of constant u, with u ∈ [0, ∞], are nested ellipses with common foci located at (x; y) = (±f ; 0). Contours of constant v, with v ∈ [0, 2π], are hyperbolae. Together these contours form an orthogonal coordinate system. See Fig. 4. Data is to be interpolated onto the ellipse whose cross section is that of the elliptical cylinder of Fig. 3. See Fig. 5.
EllipticCylindricalCoord_1000.gif (GIF Image, 475x412 pixels) http://mathworld.wolfram.com/images/eps-gif/EllipticCylindricalCoord_1000.gif For our work we will need the unit vectorê u , the unit vector (outwardly) normal to the surface of the elliptical cylinder. Write Then, by definition, we have the result It is also convenient to employ the complex variables and In these variables, the relations (24) can be written in the more compact form Form differentials of both sides of (28). Doing so gives the result and the complex conjugate result Now form the product of (29) and (30) to get the transverse line-element relation From this relation we infer the results

Mathieu Equations
Let us seek to construct harmonic functions of the form where the functions P and Q are yet to be determined. Employing the Ansatz (33) in Laplace's equation and use of (32b) yields the requirement We also observe that there is the trigonometric identity so that the requirement (34) can be rewritten in the form Upon dividing both sides by P Q, (36) becomes from which it follows that Therefore, there is a common separation constant a such that and Correspondingly, P and Q must satisfy the ordinary linear differential equations Equation (40b) for Q is called the Mathieu equation, and Equation (40a) for P is called the modified Mathieu equation. For our purposes, we will need solutions Q(v) of (40b) that are periodic with period 2π. Such solutions exist only for certain characteristic values of the separation constant a. These values are denoted a n (q) for n = 0, 1, 2, 3, · · · and b n (q) for n = 1, 2, 3, · · · . The solutions associated with the separation constants a = a n (q) are denoted ce n (v, q). They are even functions of v and, in the small q limit, are proportional to the functions cos(nv). The solutions associated with the separation constants a = b n (q) are denoted se n (v, q). They are odd functions of v and, in the small q limit, are proportional to the functions sin(nv). The functions ce n (v, q) and se n (v, q) form a complete orthogonal set over the interval v ∈ [0, 2π] and are normalized so that With regard to the solutions of the modified Mathieu equation, note that (40b) is transformed into (40a) under v → iu. As a result, corresponding (real-valued) solutions to (40a) are defined by Ce n (u, q) = ce n (iu, q) and Se n (u, q) = −i se m (iu, q). We refer the reader to [3,4] and [8] for a detailed treatment of the Mathieu functions and their properties.

Elliptic Cylinder Harmonic Expansion and On-Axis Gradients
The stage is now set to describe the expansion of any harmonic function ψ in terms of Mathieu functions. The general harmonic function that is analytic in x and y near the origin can be written in the coordinates (24) in the form where the functions c n (k) and s n (k) are arbitrary. We will call (42) an elliptic cylinder harmonic expansion.
To exploit this expansion, suppose the magnetic field B(x, y, z) is interpolated onto the surface u = U of an elliptic cylinder using values at the grid points near the surface. See Fig. 5. Let us employ the notation B(x, y, z) = B(u, v, z) so that the magnetic field on the surface can be written as B(U, v, z). Next, from the values on the surface, compute B u (U, v, z), the component of B(x, y, z) normal to the surface. Our aim will be to determine the on-axis gradients from a knowledge of B u (U, v, z). At this point we note that the functions exp(ikz)se n (v, q) and exp(ikz)ce n (v, q) form a complete set over the surface of the elliptical cylinder.
Let us begin by solving (32a) for (∂ψ/∂u). We find, using (26), the result, We see that the right side of (43) is a well-behaved function F (u, v, z) whose values are known for u = U , Moreover, using the representation (42) in (43) and (44), we may also write Next multiply both sides of (45) by exp(−ik z) and integrate over z. So doing gives the result 1 2π Now, employ the orthogonality properties of the Mathieu functions to obtain the relations In view of (47), define the functionF (v, k) by the rulẽ and define functionsF c r (k) andF s r (k) by the rules Note the similarity to (17), where cos(rφ) and sin(rφ) are replaced by ce r (v, q) and se r (v, q). We will call the functionsF α r (k) Mathieu coefficient functions in analogy to the Fourier coefficients that arise in Fourier analysis.
With these definitions, the relations (47) can be rewritten in the form Finally, employ (50) in (42). So doing gives the result We have obtained an elliptical cylinder harmonic expansion for ψ in terms of surface field data.
Of course, what we really want are the on-axis gradients. Again, once these gradients are known, we may use (11) and (12) to compute the associated scalar and vector potentials. The gradients can be found by employing two remarkable connections (identities) between elliptic and circular cylinder functions of the form [3] Ce r (u, q) ce r (v, q) = ∞ m=0 α r m (k)I m (kρ) cos(mφ), (52a) For further reference, we will call the quantities α r m (k) and β r m (k) Mathieu-Bessel connection coefficients [3,4]. Using these results, (51) can be rewritten in the form Upon comparing (53) with (5), we conclude that there are the relations and Finally, in view of (8) and (54), we have the desired results We have found expressions for the on-axis gradients in terms of field data (normal component) on the surface of an elliptic cylinder.

C. Use of Field Data on Surface of Rectangular Cylinder
A similar procedure has been developed for computing the on-axis gradients in terms of field data provided on the surface of a rectangular cylinder. In this case, each on-axis gradient C [n] m,α may be written as the sum of four contributions, Each contribution is determined by the integration of the normal component of the field against an appropriate kernel over one of the four faces (Top, Bottom, Left, Right) of the rectangular cylinder. The resulting expressions are quite lengthy, and we therefore refer the reader to [3,4] for further details. The remainder of this paper will focus on the circular and elliptical cylinder cases.

A. Monopole Doublet
In this section, we develop an exactly-soluble but numerically challenging model field to be used to numerically benchmark the procedures developed in Section III. Suppose two magnetic monopoles having strengths ±g are placed at the (x, y, z) locations r + = (0, a, 0), See Fig. 6, which also shows a circular cylinder with radius R (the surface ρ = R). These monopoles generate a scalar potential ψ(x, y, z) described by the relation Correspondingly, they produce a magnetic field B = ∇ψ having the components off-axis field components B x (ρ = 1/2, φ = π/4, z) and B z (ρ = 1/2, φ = π/4, z).
Due to the symmetries of the field, the functions C m,s = 0 for m even. It can be shown that the nonvanishing on-axis gradients are given for m odd by [3]: For fixed z, the domain of convergence for the poly-   nomial series (11)(12) representing the field in terms of the functions (60) is given by the condition x 2 + y 2 < √ z 2 + a 2 . In particular, the domain of convergence is a region of circular cross-section whose radius increases as we move longitudinally away from the location of the monopoles at z = 0.
Suppose we wish to calculate a transfer map through 7 th order. Then, as shown in the Appendix, to do so requires knowledge of the C m,α (z) with (m + n) ≤ 8 when m is even.
Graphs of a selected few of these functions, for the monopole doublet in the case that a = 2.5 cm and g = 1 Tesla-(cm) 2 , are shown in Figs. 11 through 13. In these plots z has units of centimeters. Evidently the C m,s become ever more highly peaked with increasing m. Fortunately, when working through some fixed degree, we need fewer derivatives with increasing m. Note that we expect that the function C m,s (z) should have n zeroes. This is indeed the case.

B. Circular Cylinder Results
The procedure discussed in Section III A has been benchmarked using the field of a monopole doublet in the case that a = 2.5 cm and g = 1 Tesla-(cm) 2 . We set up a regular grid in x, y, z space, where we let each variable range over the intervals The known values of the three components of the field are computed using (59) at each grid point. Consider a cylinder of radius R = 2 cm and length 600 cm. We use bicubic interpolation to interpolate B at these grid points onto 49 selected angular points on the cylinder, for each of the 4801 selected values of z. The angular integration in (17) is performed using a Riemann sum with N = 49. (This is necessary to ensure sufficient convergence of the angular integrals to within 10 −4 .) We evaluate the Fourier transform at 401 values of k in the range [−K c , K c ] with K c = 20, using a spline-based Fourier transform algorithm [3,4]. We use these same points in k space to evaluate the inverse Fourier transform, providing a set of numerically determined functions C [n] m,α (z). A comparison between the exact on-axis gradients (60) and those obtained from grid data is provided in Figs. 11-13 for the functions C 7,s . Evidently the agreement is excellent. Fig. 14 illustrates the difference between the on-axis gradient C 1,s as obtained from grid data and the exact on-axis gradient obtained from (60) with m = 1. The maximum error attained relative to peak is 1.7 × 10 −4 . Further detailed study shows that numerical results agree with exact results to within relative errors less than a few parts in 10 4 for all the relevant C [n] m,α (z). By using exact data on the cylinder rather than interpolating off the grid onto the cylinder, we have also verified that the error due to interpolation onto the cylinder is comparable to that produced by numerical integration [4]. Finally, all these small errors can be further reduced with the aid of a finer grid [3,4].  1,s (z). Exact results are shown as a solid line, and numerical results are shown as dots.

C. Elliptical Cylinder Results
The procedure discussed in Section III B has been benchmarked using grid values identitical to those described in the previous section. Consider an elliptical cylinder of semiminor axis of y max = 2 cm and semimajor axis x max = 4 cm. In this case, we evaluate the angular integrals (49) using a Riemann sum with N = 120. (This is necessary to ensure sufficient convergence of the angular integrals to within 10 −4 .) Doing so requires interpolation off the grid onto the elliptical cylinder at 120 m,α (z).

Results for the functions C
[n] m,α (z) are similar to those found in the circular cylinder case [3,4], and have comparable accuracy. In the following section, however, we illustrate that functions obtained using an elliptical cylinder are significantly more robust against numerical noise in the original grid values.

V. SMOOTHING
In this section, we investigate the smoothing of numerical noise that results from the use of the surface fitting techniques described in Section III. Consider computing the on-axis gradients C , 300] cm used in Section IV for fitting the field of the monopole doublet, with mesh points indexed by j = 1, · · · N . Let B y (0, 0, z) denote the value of the monopole-doublet on-axis vertical field at longitudinal location z, as determined from (59). At each point (x j , y j , z j ) we set Here the δ x (j) and δ y (j) are uniformly distributed random variables taking values in the interval [−1, 1], and = 0.01. After interpolating these values onto the surface of a circular cylinder with R = 2 cm and z ∈ [−300, 300] cm, we use the procedure described in Section IIIA to compute the on-axis gradients (22).
In Fig. 15, we have displayed the computed quantityb m,s (R, k) appearing in (17) for the case m = 1. It is a function of the spatial frequency k, having random variations of approximately uniform variance over the interval [−20, 20] cm −1 . Fig. 16 displays the kernel k m−1 /I m (kR) multiplyingb m,s (R, k) in (22) for the case m = 1 and R = 2 cm. Note the rapid decay of this function for large |k|. Finally, Fig. 17 displays the product of these functions, illustrating the dramatic suppression of high-k contributions to the Fourier integral appearing in (22).
A similar phenomenon occurs when fitting is performed using an elliptical cylinder. In this case, a sequence of kernels contributes to (55) for each fixed m. In Fig. 18 we have displayed the kernels contributing to the case m = 1, with x max = 4 cm and y max = 2 cm. Kernels take their maxima at k = 0, and these maxima decrease monotonically with increasing index r. All kernels decrease rapidly with increasing |k| [3,4].
To study the effect of noise on the on-axis gradients, twelve distinct random fields were generated on a mesh according to (61). Figs. 19-20 illustrate the on-axis gradients C [6] 1,c and C [0] 7,c as computed using these field values. The solid line in Fig. 19 illustrates the rms value of the on-axis gradient C [6] 1,c , as computed using a circular cylinder of radius R = 2 cm according to (22). The dashed line in Fig. 19 illustrates the rms value of the on-axis gradient C [6] 1,c , as computed using an elliptical cylinder of semiminor axis 2 cm and semimajor axis 4 cm according to (55). In Fig. 20, similar results are shown for the on-axis gradient C m,s are produced by fields that are predominantly in the vertical (y) direction, and for such fields there is relatively less difference between circular and elliptical surface fitting because the semi-minor axis of the elliptical cylinder is the same as the radius of the circular cylinder. Note also that in both cases only the component of the field normal to the surface is used. Horizontal fields drive primarily the C m,c . However, in the horizontal direction, the semi-major axis of the elliptical cylinder is substantially larger than the radius of the circular cylinder. We therefore expect the advantage of using elliptical cylinders compared to circular cylinders will be most apparent in the C  Figs. 19 and 20 show. The effects of errors in the surface data are suppressed more when the bounding surface is farther away from the field observation point. As a result of this suppression, it is advantageous to use a fitting surface that is as far away as possible.
Suppose we compare the C m,s shown in Figs. 12 and 13. Remarkably, we find that, as a result of smoothing, the effect of 1% noise in the field data produces, on average, only on the order of 0.01% changes in the on-axis gradients.

A. ILC Damping Ring Wiggler Fields
A less stringent test of the accuracy of surface methods is that the magnetic field, as computed from surface data using field values on a 3-dimensional mesh, should reproduce the magnetic field at the interior mesh points. (This is also a test of the quality of the magnetic data on the mesh.) We computed such an interior fit, and the associated transfer map, for the modified CESR-c design of the Cornell wiggler, which has been adopted as the design prototype for use in International Linear Collider studies [9,10]. Cornell provided data obtained from the   The interior field was computed using the on-axis gradients through terms of degree 6 in x, y over the domain of the original data. This solution for the interior field was then compared to the original data at each grid point. Fig. 21 displays the vertical field B y off-axis at The field data (points) are shown along with computed values (solid line). Note that the fitted field captures the fringe-field behavior. The relative error was found to satisfy the bound δ|B data − B f it |/|B| peak ≤ 3.5×10 −4 . We observe that this error is comparable to that found for the monopole-doublet benchmark. Presumably it is due to errors in numerical integration, errors in interpolating onto the elliptic cylinder, errors arising from neglecting terms beyond degree 6, etc., as well as possible failure of the OPERA-3d data to be Maxwellian. Fig. 22 illustrates the horizontal roll-off of the vertical field at y = 0.1 cm, z = 104.2 cm for 0 ≤ x ≤ 1 cm. Note the discrete jumps in the original data, reflecting the number of digits retained in the output of the numerical computation. Despite the small variation of B y in x, the fit goes through the interior data. Finally, Fig. 23 illustrates the longitudinal field B z , again at (x, y) = (0.4, 0.2) cm along the wiggler. Note that no information about B z was used to generate this field, since only the component of B normal to the elliptic cylinder surface was used to generate the interior solution.
The error for B y on-axis lies in the range 0.1-0.2 G along the length of the wiggler, increasing slightly near the end poles. A plot of residuals in the plane y = 0 is displayed in Fig. 24. Note that the error is within 0.6 G over this region of the x-z plane. The error begins to increase rapidly at about x = 2 cm; this may be due to retaining only terms through degree 6 in the on-axis gradients, or perhaps a finite domain of convergence of the associated power series for B y (x, y, z). For all |x| ≤ 2 cm, the peak error is 0.3 G. We remark that this peak error amounts to a relative error of less than 2 parts in 10 5 , which is remarkably small compared to the error for the data of Fig. 21. We expect the error to behave  (11) or (12). Dots represent numerical data provided by OPERA-3d. like a harmonic function, and therefore it should grow as one approaches the boundary. Conversely, it should be the smallest on the center line. The observed error follows this pattern. Finally, this phenomenon may also be a factor in the observed increase in the error at and beyond x = 2 cm. The design orbit is displayed in Fig. 25, and Table I lists its initial and final conditions. Table II displays the matrix R 2 that describes the linear part of the transfer map in (62). Phase-space coordinates are arranged in the order (x, p x , y, p y , τ, p τ ). The first few Lie generators f m of the nonlinear part of the transfer map are listed in Table III.

VII. CONCLUDING SUMMARY
Surface methods provide a reliable and numerically robust technique for extracting transfer maps from numerical field data. By benchmarking them against a nu-   m,α (z) required to compute transfer maps through 7 th order, and that the results obtained were remarkably insensitive to noise. Moreover, these small errors can be further reduced if desired with the aid of a finer grid.
Subsequently we applied surface methods to compute the interior field for the proposed ILC Damping Ring wigglers. Consistent with the accuracy displayed by the monopole-doublet benchmark results, excellent fits were demonstrated for interior fields. We also illustrated the computation of the design orbit and its associated transfer map based on surface methods.
In summary, the use of surface methods makes it possible, for the first time, to compute for straight beam-line elements realistic transfer maps for real magnets includ-ing all fringe and high-order multipole field effects.
In many cases, however, we are interested in magnetic elements with significant sagitta, such as dipoles with large bending angles. In these cases, it is not possible in general to surround the design orbit with a cylindrical surface that lies interior to all iron or other magnetic sources. Part II of this paper will describe an alternative, but more computationally intensive, method suitable for general geometries including that of a bent box with straight ends. See Fig. 26. In this case we obtain simple, geometry-independent kernels for computing the interior vector potential and its derivatives. All the advantages demonstrated in this paper for surface methods will be retained.

APPENDIX
MaryLie/IMPACT is a hybrid code that utilizes Lie-algebraic methods for computing and manipulating charged-particle transfer maps through 5 th order, while space-charge effects are treated using Particle-in-Cell methods [11]. We made use of its Lie algebraic capabilities. In the Lie algebraic approach, maps are computed and manipulated in Lie algebraic form. Each map describes the transformation of the full six-dimensional phase-space coordinates of a particle as it passes through a given beam-line element. Because of the symplectic nature of Hamiltonian motion, through aberrations of order (n − 1) such a map has the Lie representation where is the charged-particle Hamiltonian expressed in terms of deviation variables about the design orbit and expanded in a homogeneous polynomial series. The deviation variable Hamiltonian H is determined in turn by the Hamiltonian K for which path length is the independent variable. In Cartesian coordinates and with z taken as the independent variable, K is given by the relation Here Φ and A are the electric scalar and magnetic vector potentials, respectively. We conclude that (in the case of no electric fields, Φ = 0) what is needed are Taylor expansions for the components of A in the deviation variables x and y. In the straight-element case these Taylor expansions can be found using the circular cylinder harmonic based expansion (12). Suppose, for example, we wish to retain in the expansion of the Hamiltonian H appearing in (64) homogeneous polynomials through degree 8. That is what is required to compute transfer maps through 7 th order. If the design orbit lies on the z axis, as will be the case for any straight multipole such as a solenoid, quadrupole, sextupole, octupole, etc., this expansion is straightforward because in this case the Cartesian coordinates x, y are already deviation variables. We see from (65) that we must retain homogeneous polynomials in the variables x, y through degree 7 in the expansions of A x and A y , and homogeneous polynomials in the variables x, y through degree 8 in the expansion of A z . Inspection of (12) shows that for the cases m = 0 or m odd we then need the C m,α (z) with (m + n) ≤ 8.
In the case of a wiggler, the design orbit oscillates around the z axis. See Fig. 25. Now the requirement that the Hamiltonian H appearing in (64) be an expansion in deviation variables about the design orbit is more involved: The components of A must be expanded about the design orbit. As they stand in (12), they are expanded about the z axis. If we are to retain terms in H through degree 8 when a re-expansion is made about the design orbit, then coefficients C m,α beyond those listed in the previous paragraph must in principle be included in the calculation. We say that higher-degree terms produce lower-degree terms due to feed down. However, the effect of feed down is small if the oscillations are of small amplitude, as they are for the proposed ILC wiggler. In this example we have found that the approximation of retaining only the terms in A x , A y , and A z through degree 8 produces changes in the design orbit and the transfer map about the design orbit that are comparable to or smaller than the errors found in the benchmark studies of Section IV. We also remark that this feed-down problem does not arise when the general geometry methods to be described in Part II are employed.