Numerical computation of high-order transfer maps for rf cavities

Modern map-based accelerator beam-dynamics codes model magnetic elements so as to include nonlinear effects and realistic fringe ﬁelds, but they persist in modeling rf cavities as either energy kicks or linear maps. This work presents a method for including the nonlinear effects of rf cavities in a map-based code.


I. INTRODUCTION
Techniques now exist for constructing accurate transfer maps for very general magnetic elements, including fringes and overlapping fields [1][2][3]. For radio-frequency (rf) cavities, however, current modeling practices usually compute transfer maps as either energy kicks that vary with cavity phase, e.g. [4], or linear maps [5]. Both approaches omit significant physics: Using just an energy kick ignores the transverse particle dynamics important in, for example, high-gradient rf linear accelerators [6,7], rf cavities that contain axial asymmetries [8,9], or rf cavities designed for purposes other than acceleration [10 -13]. Using a linear map in a z-based code implies the absence of second-order energy variation -important when accelerating near peak voltage. The work described here aims to place the treatment of rf cavities on the same footing as the present treatment of magnetic elements in modern nonlinear beam-dynamics codes.
Rosenzweig and Serafini [6] presented a generalized matrix approach -generalized in the sense that phase information is included in the matrix entries -for ultrarelativistic particles in an axisymmetric and periodic rf cavity. Van Zeijts [14] brought to bear the power of Hamilton's machinery [15]. The advantages of this approach include the fact that computing transfer maps to high order can be automated; other fields-from nearby magnets or other cavity modes-can easily be superposed; and the inclusion of axial asymmetries is straightforward.
To compute transfer maps for rf cavities, one must know the vector potential; in particular, one must have a transverse expansion of the vector potential at many longitudinal locations. One may construct the coefficients of such an expansion from the on-axis field and its derivatives. For a realistic cavity, however, the on-axis field is known only from experimental measurement or electromagnetic simulation, and numerically computed derivatives become increasingly suspect as the order increases. This paper describes a method for computing robustly and to high order the coefficients in the transverse expansion of the vector potential for a given rf cavity. Moreover, the technique uses Fourier transforms (as opposed to Fourier series), so it makes no implicit assumption about the periodicity of the rf structure and therefore correctly models the fringe fields.
Section II outlines very briefly the computation of transfer maps. Section III describes a general expansion for the fields of an rf cavity and how to obtain expansions from surface field data. Section IV presents a model for a realistic rf cavity, including fringe fields, for which fields may be computed by a straightforward quadrature. It then uses the model to test the method of this paper. Section V shows results based on data from an electromagnetic simulation. This paper concludes with a summary. A sequel to this paper will present the results of computing transfer maps and using them to describe particle motion in rf cavities.

II. COMPUTATION OF TRANSFER MAPS
For only a handful of dynamical systems can one compute a transfer map in closed form. In all other cases one must resort to approximations. For Hamiltonian systems there exists a well-developed calculus, using methods of Lie algebra, for computing transfer maps to high order [15,16]. It rests on the fact that maps for such systems must be symplectic. Here we give a very brief overview.
First, define for any function f of dynamical variables, coordinates q and momenta p, an associated Lie operator :f: by the rule that it acts on any other such function by taking a Poisson bracket: @f @q i @g @p i ÿ @f @p i @g @q i : Then powers of this Lie operator, naturally enough, we define by iteration: :f: 0 g g; :f: 1 g f; g; :f: 2 g f; f; g; :f: 3 g f; f; f; g : Finally, the corresponding Lie transformation is the operator e :f: These operators have the very special property that for any dynamical function f, the Lie transformation e :f: constitutes a symplectic map. Moreover, the degree of the map generator f is associated with the degree of nonlinearity in the map: A homogeneous polynomial of degree two, commonly denoted f 2 , gives rise to a linear map; while a homogenous polynomial of degree n 3, an f n , gives rise to a map that in Taylor series form possesses nonlinear terms of degree n ÿ 1 and higher. The Dragt-Finn factorization theorem [17] shows that, for Hamiltonian systems, a transfer map M possessing a convergent Taylor series expansion may also be written as a (usually infinite) product of Lie transformations, where each succeeding generator is of successively higher degree: e :f 2 : e :f 3 : e :f 4 : . Much effort has gone into the building of computer codes that can construct, concatenate, analyze, and apply transfer maps in this form [4].
which we solve with the initial condition M0 equals the identity. This equation amounts to a set of coupled ordinary differential equations to be solved for the Lie generators f n [18]. The computer codes MARYLIE [4] and MARYLIE/ IMPACT [19] each contain a set of routines -the GENMAP routines-which solve the equation (1) for a given Hamiltonian to produce a transfer map in Dragt-Finn form.

III. HAMILTONIAN FOR A RELATIVISTIC CHARGED PARTICLE IN AN RF CAVITY
One may write the Hamiltonian for a relativistic charged particle in electromagnetic fields, described by scalar and vector potentials and A, respectively, in the standard form H m 2 c 4 p ÿ eA 2 c 2 q e : To use the longitudinal coordinate z, rather than the time t, as the independent variable, one notes that the ''momentum'' conjugate to t is p t ÿH and then solves for p z to determine the new Hamiltonian [15]: H ÿ p t q 2 =c 2 ÿ mc 2 ÿ p x ÿ qA x 2 ÿ p y ÿ qA y 2 q ÿ qA z : For the purposes of constructing a transfer map, we shall require an expansion, to arbitrary order in the transverse variables, of this Hamiltonian, and hence of the potentials and A. The electric and magnetic fields derive from the potentials and A according to E ÿr ÿ @ t A; B r A: One may, of course, use an arbitrary gauge function to transform these potentials to 0 ÿ @ t ; A 0 A r without altering the fields E and B [20]. In particular, one may choose a gauge such that the scalar potential vanishes. In this gauge the electric and magnetic fields derive solely from the vector potential: Moreover, assuming a simple-harmonic time dependence for the fields, one converts the electric field E to the vector potential A by a simple division plus a 90 phase shift. We therefore seek an arbitrary-order expansion for the electric field E in terms of the transverse variables. To this end, we first develop a completely general Bessel function expansion for the spatial component of a given cavity mode. We then show how one may use surface field data -from, say, three-dimensional electromagnetic simulations -to determine accurately the coefficients in the Bessel function series. After expanding the Bessel functions in the transverse variables, we will obtain the desired power series expansion of E, hence also of A.

A. General representation of fields in an rf cavity
Assuming the vacuum form of Maxwell's equations applies to our rf cavity, we know the electric field satisfies the homogeneous wave equation Here theẽ m k andf m k denote arbitrary functions of k which describe the particular mode in the particular rf cavity. The radial function R m is either a regular or a modified Bessel function, depending on the relative values of k l and the separation constant k. Define s l k sgnk 2 ÿ k 2 l ; 2 l jk 2 ÿ k 2 l j; (9) so that k 2 ÿ k 2 l s l k 2 l . Then R m k; J m l ; s l k < 0; I m l ; otherwise: For future reference, we note that R m k; has the derivative obeys the recursion relation s l kR m1 k; ÿ 2m l R m k; R mÿ1 k; ; (11b) and satisfies the differential equation In addition, R m has the series expansion Because we are working in cylindrical coordinates, the remaining two components of the vector Helmholtz equation, (7a) and (7b), constitute a pair of coupled partial differential equations for E and E . We shall write E and E in the general forms where F cm , F sm , F cm , and F sm denote radial functions to be determined. Inserting (12) into (7a) and (7b) yields the following two pairs of coupled differential equations: Consider first the case m 0, for which we need only F c0 and F c0 . Here, noting (9) and (11c), we immediately identify both radial functions as proportional to R 1 k; . In addition, the m 0 component of the divergence equation (6) tells us that It follows from (11a) that Since the equation for F c0 has the same form as that for F c0 , we write F c0 , by analogy, in the form Then an azimuthally symmetric mode has the general representation For accelerating modes, of course,f 0 0 and E vanishes.
For the remainder of this section, we assume m 1. To decouple the equations (13), we insert (8) and (12) into the divergence equation (6), obtaining the constraints Using these in (13a) and (13c), we obtain separate nonhomogeneous differential equations for the radial functions F cm and F sm that define E : A straightforward computation shows that ÿ ik lẽ m kR m1 k; and ÿ ik lf m kR m1 k; constitute, respectively, particular solutions to (17a) and (17b), and that R m k; l constitutes a solution to the homogeneous versions of (17) which is regular at the origin. We therefore obtain for F cm k; and F sm k; the general solutions where m k and m k denote arbitrary functions to be determined for each particular cavity mode. To determine the remaining two radial functions, F cm and F sm , we use (16) and (18) to obtain Putting together (8), (12), (15), (18), and (19), we obtain the following completely general representation for the spatial component of the electric field for a particular standing-wave mode in an rf cavity: Note that the functionsẽ m ,f m , m , and m completely determine the electric field and are independent of radius.
(Moreover, because E is simply related to the vector potential A, they also determine the associated magnetic field.) We turn now to the task of determining those functions.

B. Determination ofẽ m ,f m , m , and m
Since the functionsẽ m ,f m , m , and m do not vary with radius, one may determine them using any convenient radius. Given E on the surface of a cylinder of radius R about the z-axis, an azimuthal Fourier decomposition of the z component yields For E and E we write exactly similar expressionsexcept that in the case of E we shall write E cm cosÿE sm sin. Comparing (21) with (20c), equating coefficients of cosm, and inverting the Fourier integral yieldsẽ for m 2 f0; 1; 2; . . .g. Similarly, the coefficients of sinm yieldf for m 2 f1; 2; . . .g. To determinef 0 k, we require the ''constant'' term in the azimuthal Fourier decomposition of E : It remains for us to determine m k and m k for m 2 f1; 2; . . .g. This one may do using either E or E together with the already determinedẽ m k andf m k. For example, from the cosm component of E we determinẽ From the sinm component we determinẽ An alternative expression for m k derives from the cosm component of E : To obtain the corresponding alternate expression for m k, replace in (25b)Ẽ cm andf m byẼ sm andẽ m . Although and do not appear in the m 0 fields (15), it is convenient -and consistent with the above expressions [22]-to define 0 k 0 k 0: (26)

C. Expansion of E in the transverse variables
To expand the electric field E in the transverse variables, we first expand in terms of the radial variable by inserting (11d) into (20). After some algebra, we obtain the following expansions, valid for all m 0: [Note that C c00 z and C c00 z both vanish, so that in the case m 0, the j 0 term makes no contribution to the expressions for E c0 and E c0 .] The corresponding expansions for E sm , E sm , and E zsm and C smj , C smj , and C zsmj (m 1) are obtained by interchangingẽ m and m withf m and m , respectively. For the important special case of an azimuthally symmetric accelerating field-m 0,f 0 0-it is well known that the electric field is determined entirely by its on-axis behavior. We can confirm this from the above expansions: In this case, according to (26), m and m make no contribution. Then (20c), (27), and (28), with (9), yield In particular, for the longitudinal field on axis, 0 and only the j 0 term contributes: in agreement with (20c). It follows that one may rewrite (29) in the form demonstrating that, for the case m 0, the on-axis field does indeed determine the entire electric field. In addition, a comparison of the above expressions with (29) leads us to call the functions C mj z generalized gradients.
According to the results in (31), one may obtain the coefficients in the transverse expansion of a given electric field by computing various derivatives of the on-axis component E z 0; z. While straightforward, this technique suffers from the fact that repeated numerical differentiation is an unstable operation, exquisitely sensitive to the presence of even small amounts of noise. On the other hand, (29) shows that one may compute the very same coefficients by performing various integrals of the cavity's characteristic functionẽ 0 k, essentially, see (22), the Fourier transform of the longitudinal field. This latter approach, as a consequence, yields a much more robust technique, capable of obtaining the expansion coefficients, or generalized gradients, to very high order. (See Sec. IV for an explicit example.) In addition, even if one performs the integrations of (28) rather crudely, the fact that integration is a linear operation means the computed fields can still satisfy Maxwell's equations through the order of the computation. In particular, this will be the case whenever taking the longitudinal derivative commutes with the operation of converting from integral to quadrature formula.
To obtain expansions in terms of the transverse variables x and y, one makes the replacements To complete the process, one then also converts the transverse fields to Cartesian components according to the rules E x E cos ÿ E sin; E y E sin E cos:

D. Expansion of A in the transverse variables
To construct the Hamiltonian, we require not the electric field E, but the vector potential A. According to (2) and (4), this is just (for a given mode) As an explicit example, we write down the spatial component of the vector potential, in Cartesian coordinates, for the azimuthally symmetric accelerating cavity: The ÿi out front serves simply to convert a cos! l t l modulation of the electric field [cf. (4)] to a ÿ sin! l t l modulation of the vector potential.

IV. EXAMPLE RF CAVITY
To see how well the approach outlined above may work in practice, we examine how errors (or noise) in surface field data affect the computation of interior fields. First we develop a simple, though nontrivial, model rf cavity for which one may compute the fields analytically. Then we generate appropriate discrete surface field ''data'' and use it to investigate the sensitivity of the interior field computations. According to (31), the on-axis derivatives of the longitudinal field compose the generalized gradients, the coefficients in the transverse expansion (29) of the electric field. We therefore examine first the computation of those derivatives, and then the actual fields.

A. Analytic cavity
Consider an azimuthally symmetric -mode cavity with five equally spaced identical cells and the property that at the bore radius a the longitudinal field E z has fixed magnitude E g in the gaps (alternating in sign from gap to gap) and 0 elsewhere (see Fig. 1). In this case one may compute L g the characteristic functionẽ 0 k analytically: Using (22), and assuming the central cavity has bore field E g , one obtains where g, L, and z c denote, respectively, the gap size, cell length, and longitudinal position of the cavity midpoint. (For a number of cells different from five, only the form factor, in square brackets, differs-though it remains a function of kL.) Then, according to (30) and (34), the on-axis derivatives are kg=2 1 ÿ 2 coskL 2 cos2kLk n ( ÿ1 n=2 coskz ÿ z c ; neven; ÿ1 n1=2 sinkz ÿ z c ; n odd: [For derivatives off axis at radius , one must include an additional factor of R 0 k; in the integrand above.] Figure 2 shows the result of using (35) to compute the on-axis longitudinal field E z 0; z as a function of z for the case L 0:2125 m, g L=2, and a 0:4L 0:085 m, which has a corresponding cavity frequency c=2L 705 MHz. (These parameter values are very close to those of the five-cell cavities being designed for the electron LINAC required by the electron-cooling system currently under investigation for the RHIC luminosity upgrade.) Similarly, Fig. 3 shows the first six derivatives of the field in Fig. 2, required for computing particle motion correct through fifth order in the transverse variables. For these two figures, the integrals were computed using the built-in numerical integration routines of MATHEMATICA [23].
The indentations seen in the peaks of the first derivative (top-left plot in Fig. 3) derive from the very slight flattening of the electric field between one gap and the next. This flattening, though not visible in Fig. 2 (see instead Fig. 4 below), makes its presence felt increasingly with each successive derivative.

B. Using surface field data
To simulate the process of computing fields from numerical data, we first generated data using (34) in (20c), with m 0, to compute field values for discrete points at . To obtain accurate values for fields at this radius, one must set the upper integration limit to a k max of nearly 1400. The field values are shown in Fig. 4. We then used this data -with varying amounts of added error-to compute the characteristic functionẽ 0 k, (22) with m 0, and then the on-axis field values and longitudinal derivatives, based on (30): Knowing the asymptotic behavior ofẽ 0 k-from (34), / k ÿ1=2 e ÿka -one may for any desired accuracy identify a maximum k at which to truncate this (and other) integrals. For the code MARYLIE/IMPACT we shall need C zcoj for j 2 f0; 1; 2; 3g. This implies that we need E 6 z 0; z, hencẽ e 0 k for k as high as 800. The numeric integrals in (22) and (36) were computed using Filon integration [24]. Figures 5-7 compare the results obtained with 10% noise-meaning surface field values (the data) were multiplied by uniformly distributed random numbers in the range 0:9; 1:1-with the analytic results obtained earlier.
In particular, Fig. 7 shows as functions of z, and for n 2 f0; . . . ; 6g, the relative errors between the values,Ê n z 0; z, obtained from noisy data and the analytic values, E n z 0; z, obtained earlier of the nth longitudinal derivative of the on-axis field. We see that,  even with a relative noise level of 10%, the on-axis field values remain good to 1%; and not until the fourth or fifth derivative does the relative error approach 10%. In general, the relative error in the on-axis field is about a factor of 10 smaller than that of the surface field data; and only for the fourth or fifth derivative does the relative error approach that of the surface field data. This behavior results largely from the computation ofẽ 0 k in (22): in particular, the integration over z and the presence of R 0 k; R in the denominator both serve to damp errors.
We turn now to the actual computation of the field or vector potential -both on and off the axis: Figs. 8 and 9 show, as functions of z, and for several different radii, the error made in using the first three terms of (29) to compute the longitudinal and radial electric fields, respectively. In Fig. 8 two of the curves-those representing the errors on axis and at one-tenth the bore radius-have a character very different from the others: An examination of the truncation error shows that this dominates the curves for radii 0:2a. (See Table I for the amplitudes of the generalized gradients.) The curve for 0:1a exhibits some truncation error in the central 30 cm of the cavity, but it is otherwise dominated by the accuracy of the computed (reconstructed) on-axis field, as shown by the curve for 0 which necessarily has zero truncation error. For the radial fields, see Fig. 9, all the curves shown are dominated by truncation error.   Fig. 2, and an approximate computation, the first three terms of (29b). Results are shown for various fractions of the bore radius a: 0 (black), 0:1a (red), 0:2a (orange), 0:3a (green), 0:4a (light blue), 0:5a (dark blue), and 0:6a (violet). If errors are present in the surface field data, then these will propagate (somewhat damped) into the generalized gradients. For the longitudinal field, the jth generalized gradient goes-cf. (29b) and (31b)-as the 2jth derivative of the on-axis field. To the extent that the first term in (29b) dominates the computation, we expect the error to be dominated by the error in the on-axis field. This claim is illustrated in Fig. 10. The surface field data used to construct that figure includes 0.1% relative noise; and a comparison with Fig. 8 shows that the errors at the smaller radii have increased to roughly the 0.01% level (recall the factor of 10 damping), while the errors at the larger radii -already larger than 0.01%-remain dominated by truncation error and hence see no significant increase. Similar considerations apply to the radial field, and the radial version of Fig. 10 we show in Fig. 11.

V. SIMULATED RF CAVITY
In general, the data one would use to implement the technique described in this paper would come from an electromagnetic simulation of the rf cavity in question.
Here we have used the results of a SUPERFISH simulation of a five-cell rf cavity having cell length 21.3 cm, bore radius about 9.5 cm, and frequency 703.75 MHz. One can extract surface field data from that simulation, use it to compute interior fields, and compare with the SUPERFISH data at those interior points. The results, based on data at radius 9 cm, are shown in Figs. 12 and 13, which have the same format as Figs. 8 and 9.
Based on the same sorts of observations as in the last paragraph of the previous section-only in reverse-Figs. 12 and 13 suggest that the SUPERFISH data is accurate to about four digits. In fact, the SUPERFISH   radii up to about 4.5 cm, according to Table I, keeping terms in (29) through j 3 results in a truncation error comparable to the accuracy of the SUPERFISH results. To achieve this, however, we must know the sixth derivative of E z 0; z-the dominant contribution to C zco3 z-to at least two digits, a demand that is difficult to meet using numerical differentiation, but straightforward using (22) and (36).

VI. SUMMARY
To construct transfer maps for rf cavities, one must know the vector potential, which in the gauge employed here differs from the electric field simply by a factor i! l . The coefficients in the transverse expansion of the field/vector potential are, see (27) or (29), the generalized gradients (28). This paper describes a robust technique for computing those generalized gradients to the high order demanded by modern map-based tracking codes. This technique is not restricted to azimuthally symmetric cavities, or even to accelerating cavities. Moreover, one may superpose other fields by adding their vector potentials and thereby explore the effects of higher-order modes or overlapping magnetic fields.