New method for generating linear transfer matrices through combined rf and solenoid fields

We present a new method for computing the transverse transfer matrix for superimposed axisymmetric rf and solenoid field maps. The algorithm constructs the transfer matrix directly from one-dimensional rf and solenoid field maps without computing numerical derivatives or eigenfunction expansions of the field map data. In addition, this method accurately describes the dynamics of low energy particles starting from a solenoid-immersed cathode, allowing the method to simulate transport through both rf and electrostatic guns. Comparison of particle tracking with the transfer matrix, and direct integration of the equations of motion through several field setups, shows excellent agreement between the two methods.


I. INTRODUCTION
Linear transfer matrices continue to serve as important simulation tools in the design, commissioning, and operation of modern accelerators.Common examples of their use include the computation of linear centroid motion, beam-based alignment of optical elements, and orbit feedback.In addition, they feature prominently in the theory of round-to-flat beam transforms in rf and DC guns [1,2], and are used for emittance measurements for beams without space charge [3].Because of this utility, analytic expressions for the transfer matrix through many beam line elements, such as magnets with constant fields, are well known and are widely used.In contrast to these simple elements, the fields of many beam line elements in modern accelerators have no analytic form and may overlap each other.For example, in high brightness electron sources, solenoid fields used for emittance compensation may overlap the accelerating fields at or near the cathode.To properly describe the dynamics in these machines, the transfer matrix through superimposed rf and solenoid fields must be constructed.In general, to model these elements, one must use numerically computed electromagnetic field maps.Unfortunately, no closed form solutions for the transfer matrix through such elements exist.
Nonetheless, a significant amount of work has gone into developing both semianalytic and numerical techniques to compute these matrices.In general, these techniques require some form of manipulation of the field map data and often have a limited range of validity.For example, the widely used rf transfer matrix given by Rozensweig and Serafini [4] requires a Floquet expansion of the on-axis rf field map, and is only valid for ultrarelativistic particles.Other methods expand the field map data in terms of the general solution to the homogeneous Maxwell equations in cylindrical coordinates [5,6].From this expansion, the vector potential can be computed, allowing the Hamiltonian for the overlapping solenoid and rf fields to be constructed.When combined with differential and Lie algebra techniques, this method is quite powerful and can be used to generate arbitrary order maps.Despite this, it requires significant overhead in setting up, and it may not be suitable for online modeling purposes.A simpler solution has been put forth in [7].In this approach, the particle energy and rf fields are assumed constant over a small time step.The second order linear transverse equations of motion can then be solved ''exactly,'' and the transfer matrix for the time step constructed.While this method works well for relativistic particles, it does have several drawbacks.First, the assumption that the energy is constant requires extremely small time steps for very low energy particles like those emitted from a photocathode.Second, the determinant of the resulting transfer matrix is only correct to first order in the step size.Another simple approach is given in [8,9], and makes use of the ''equivalent field'' concept.While this method guarantees the correct determinant of the transfer matrix, it assumes a constant velocity over each integration step and is therefore not valid for very low particle energies.Alternatively, one can also compute the transfer matrix elements from the equations of motion directly using numerical integration [10].In general, this requires computing numerical derivatives of the field maps (or differentiating orbits).Additionally, in order to properly capture the dynamics of the problem, some care must be taken to ensure the integrator used truly maintains symplecticity [11].
To our knowledge, there is currently no simple, inherently symplectic method of computing the transfer matrix directly from superimposed rf and solenoid field maps for low particle energies.Because no analytic solution for the transfer matrix through the field of an entire rf cavity or solenoid exists, any attempt at constructing a simple algorithm for low energy dynamics will require multiple steps.Given the availability of numerical integrators, any new algorithm must satisfy the following three requirements in order to be useful: (i) using reasonable step sizes, the method should be able to describe the low energy dynamics found in rf or DC guns; (ii) the method should generate matrices with the correct determinant, regardless of the step size; and (iii) the method for constructing the transfer matrix should be easy to implement.Based on a new analytic solution to the equations of motion, we derive a transfer matrix algorithm with all three of these qualities.
The layout of this work is as follows.First, the longitudinal equations of motion are solved for a single small step in the independent variable.Then, the transfer matrix over the same step is derived for electrostatic and solenoid fields.Building on this result, the matrix for combined rf and solenoid fields is computed.This matrix is then tested with tracking through a DC gun, a superconducting rf cavity, and a rf gun with a solenoid-immersed cathode.Excellent agreement between the transfer matrix and direct integration of the equations of motion is demonstrated in all three cases.

II. DERIVATION OF THE TRANSFER MATRIX
In order to treat both standing wave and traveling wave rf structures, the rf fields are written in complex phasor form: E ¼ Ẽðx; y; zÞe i!t and B ¼ Bðx; y; zÞe i!t .In these and all subsequent expressions, tildes are used to denote phasor quantities.For notational simplicity, the phase factor e i!t and the real symbol Re½. . .are suppressed in all but the final results for the transfer matrix.The beam line axis is taken to point along the positive z axis.Derivatives with respect to z are denoted with a prime: f 0 df=dz.The equations of motion for the longitudinal phase space variables t and are Here the constant E e ¼ mc 2 =e gives the (signed) rest energy of the electron in (eV).Assuming cylindrically symmetric fields, the linearized transverse equation of motion in the Larmor frame takes the form [12,13] In this expression L stands for both the transverse Larmor coordinates x L and y L , p ¼ is the normalized reference particle momentum, and Á L is the Larmor angle, defined by Note the use of a negative sign in front of the integral.With this definition, a positive solenoid field creates a positive change in the Larmor angle for electrons.Because all of the following derivations are carried out in the Larmor frame, the subscript L on the transverse variables is suppressed in the remainder of this work.As stated before, general solutions to the differential equations in Eqs. ( 1)-(3) do not exist for arbitrary cavity and solenoid field maps.As a result, any method for computing the transfer matrix analytically requires some form of approximation to these equations.The approach taken in this work is to find an exact solution to the equations of motion for a step in z and t that is small enough so that the field profiles Ẽz ðzÞ and Bz ðzÞ, as well as the rf phase, do not change appreciably.The solution is exact in the sense that the particle energy changes correctly over the course of the step.The change in the field map profile as well as the rf phase are then included with the use of edge-focusing matrices.Slicing the field maps and consecutively multiplying the matrices for each step gives the total transfer matrix: The first step in constructing the transfer matrix for one step is to solve the longitudinal equations of motion in Eq. ( 1).To do so the electric and magnetic fields are assumed constant and equal to their average value over the step from z i to z f : For a constant electric field, 0 ffi hE z i=E e is constant, and the normalized energy, momentum, and velocity are given by Using these expressions, the derivatives of tðzÞ and Á L ðzÞ in Eqs. ( 1) and ( 3) can be directly integrated: In the last line, the constant b is the normalized solenoid field defined as b ¼ p Á Á 0 L ¼ ÀehB z i=2mc [14].Note that the expression for time given here is essentially the same as Eq. ( 7) in [15].For simplicity and speed, we take the average values here to be equal to the fields evaluated at the midpoint of the step.In addition to defining the transformation between the lab and Larmor coordinates, the function Á L plays an important role in the derivation of the transfer matrix for both the electrostatic and rf field cases.

A. Overlapping electrostatic and solenoid fields
The derivation of the transfer matrix for superimposed rf and solenoid fields is based in part on the method used to derive the transfer matrix for static fields [3].In this section, a detailed derivation of the static field result is given.The same techniques are then modified and used to derive the rf matrix in the following section.For static fields, the equation of motion is found by taking !!0 in Eq. ( 2): Note that p 0 / hE z i and Á 0 L / hB z i in this expression, implying depends on both the accelerating field and its gradient, as well as the solenoid field.The transfer matrix from z i to z f ¼ z i þ Áz is derived in a three-step process.Over the interval ½z i ; z f , the electric and magnetic fields are approximated as rectangular step functions.Formally the fields are written as where ðzÞ is the Heaviside step function, and ðzÞ is the Dirac delta function.The transfer matrix is then found by solving the transverse equation of motion piecewise from z i to z f .First, the equations of motion are integrated across the rising edge of the electric field at z i .Because the rising edge is approximated as an instantaneous step, the particle's position does not change: Integrating the equation of motion gives the kick delivered to the particle's trajectory [4,14]: The corresponding transfer matrix for the rising edge, R E , takes the form Next, the equation of motion is solved across the interval ðz i ; z f Þ, where both the electric and magnetic field are approximately constant.In this region, the equation of motion reduces to where Á 0 L ðzÞ / 1=pðzÞ.With the electric field held constant, , p, , and Á L are given by Eqs. ( 6) and (8).With these functions, the differential equation can be solved by assuming ¼ ðÁ L Þ. Plugging this into Eq.( 12) gives Using Eq. ( 8), it is possible to show Á 00 ¼ À 0 Á 0 = 2 , canceling the second term in the above equation.Assuming Á 0 Þ 0, the first term in this expression must also vanish.
Completing the initial value problem for this solution determines the transfer matrix for the step from z i to z f : In this and following expressions, C cos½Á L ðzÞ, and S sin½Á L ðzÞ.The last step in constructing the full matrix for the interval ½z i ; z f is to evaluate the transfer matrix for the falling edge of the accelerating field.The result is essentially the same as before, except now the derivative of the electric field has the opposite sign.This allows the transfer matrix for the falling edge to be written as Combining the three matrices for each region gives the full transfer matrix for the step Áz: One important thing to note about the matrix in Eq. ( 16) is that it has the correct determinant for the phase space variables chosen: detðÁM dc x;x 0 Þ ¼ p i =p f .The transfer matrix for the canonical phase space variables x and p x can be found by applying the transformation: It follows from this expression that the matrix ÁM dc x;p x satisfies the symplectic condition detðÁM dc x;p x Þ ¼ 1.In addition to this, the transfer matrix also has the convenient feature that the derivative of the accelerating field never has to be calculated, bypassing the need to compute derivatives numerically.It is important to note that the risingedge matrix should not be included in Eq. ( 16) when the NEW METHOD FOR GENERATING LINEAR TRANSFER . . .Phys.Rev. ST Accel.Beams 15, 024002 (2012) 024002-3 electric field is nonzero at the initial position of the reference particle.Including the edge matrix in this case is equivalent to a particle seeing the field rise from zero to the actual value at the initial reference position.For particles starting from a cathode, this is not physical.Similarly, the falling edge matrix should not be included when tracking of the reference particle stops in a region of nonzero electric field.

B. Overlapping rf and solenoid fields
With the results for the electrostatic field determined, it is now possible to construct the transfer matrix for rf fields.The approximation used here is similar to that used in the electrostatic case: the field profiles Ẽz and B z , as well as the rf phase, are assumed constant over the step Áz.This implies ðz; tÞ is also constant, allowing the solutions to the longitudinal equations of motion in Eq. ( 6) to be used.Over the interval ðz i ; z f Þ, the general equation of motion in Eq. ( 2) reduces to Because the rf magnetic focusing term proportional to ! scales as p À1 and not p À2 / ðÁ 0 L Þ 2 , the solution to the equation of motion used in the electrostatic case is no longer valid.In the constant field and phase approximation, this term introduces a square root of a quadratic in z into the equation of motion: p À1 ¼ ½ð i þ 0 ÁzÞ 2 À 1 À1=2 .The presence of this factor makes this equation difficult to solve.To our knowledge no analytic solution exists.
In order to proceed, the rf magnetic focusing term must be removed from Eq. (18).Changing variables to reduced coordinates, defined by ¼ ffiffiffiffiffiffiffi p , provides a clue as to how to remove this term.To see this requires transforming the general equation of motion in Eq. ( 2) first, and then making the constant field and phase approximation.The general equation of motion for the reduced variables is [12,14] Unfortunately, the rf magnetic focusing term, now proportional to p À3 , still contains a square root in the denominator, and the reduced equation of motion does not have an analytic solution.It is important to note, however, that switching to the reduced variables effectively removes the rf electric contribution (the term proportional to Ẽ0 z ) to the focusing function [14].This implies that a similar variable transformation can be used to eliminate the rf magnetic focusing term.It turns out that such a transformation is possible, and requires switching the independent variable from longitudinal position to time.Doing so yields the equivalent equation of motion: Next, the coordinates are transformed to the new reduced variables defined by " ¼ ffiffiffi ffi p .In matrix form, this transformation is written as As before, the equation of motion for the reduced variables takes the form of Hill's equation: Comparing Eqs. ( 19) and ( 22) shows that switching the independent variable and using the new reduced coordinates effectively exchanges the roles of the functions and p.Additionally, switching between the two sets of reduced coordinates allows one to choose which of the two rf focusing terms to hide in the variable transformation.
With the new definition of the reduced coordinates, all that remains now is to solve Eq. ( 22) in the region where the solenoid and accelerating fields are constant.This requires knowing the functions pðtÞ and ðtÞ.The momentum is easily found by rearranging Eq. ( 7): pðtÞ ¼ p i þ c 0 ðt À t i Þ.The normalized energy is then given by ðtÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p 2 ðtÞ þ 1 p . Inserting these expressions into the equations of motion in the constant field region yields Note that by approximating the fields as constant after changing variables, the rf electric focusing term proportional to Ẽ0 z À3 in Eq. ( 22) is set to zero in the transformed equation of motion.This eliminates the presence of any square roots in the resulting focusing function, and effectively breaks the normal equivalence between using the two sets of independent variables and reduced coordinates.
Equation ( 23) can be solved with the function " The transfer matrix is found by completing the initial value problem for this solution.The resulting matrix can be written in the compact form In this expression _ f ¼ c f 0 .The matrix in Eq. ( 24) correctly describes the evolution of the reduced variables.In order to get the transfer matrix for the usual phase space variables, it must be transformed back: COLWYN GULLIFORD AND IVAN BAZAROV Phys.Rev. ST Accel.Beams 15, 024002 (2012) 024002-4 It is clear from this expression that the effects of the rf magnetic focusing must be contained in the matrices left multiplying M dc i!f .Combining these matrices together gives In the limit that Át ¼ t À t i is small, this matrix reduces to From Eq. ( 2), it is clear that this is nothing but a thin lens approximation to the focusing delivered by the rf magnetic field.The complete matrix for the step Áz is found by including the effects of the field edges.The effect of adding the time dependence to the edge matrices is minimal; the rising-edge matrix remains the same.For the falling edge matrix, the change in the rf phase is included in the electric field: 0 f ¼ 0 e i!Át .Adding the edge matrices to Eq. (25) gives This is the main result of this work.As a reminder, 0 ¼ e Ẽz ðz i Þ=mc 2 , where Ẽz is the complex electric field map.Additionally, the expressions for , p, and Át can be found in Eqs. ( 6)-( 8), respectively.For clarity we leave the result in the above factorized matrix form.This allows several limiting cases to be easily evaluated.First, in the limit that both the rf and solenoid fields vanish, b, 0 !0, and Eq. ( 28) reduces to a drift matrix.For vanishing rf fields, 0 !0, and Eq. ( 28) reduces to a hard edge solenoid matrix.In the limit that !!0, R rf reduces to I 2Â2 , and the total rf matrix reduces to the previous electrostatic result in Eq. ( 16).In addition to having the correct limiting behavior, the rf transfer matrix also has the correct determinant: det½ÁM rf x;x 0 ¼ det½ÁM dc x;x 0 ¼ p i =p f .This follows directly from the fact that det½R rf ¼ 1.The factorized matrix form in Eq. (28) also clearly demonstrates the dynamics of a step through overlapped rf and solenoid fields: focusing from the change in the accelerating field, rotation from the solenoid, compression from acceleration, as well as focusing from the rf magnetic field.Additionally, because the analytic form depends only on the reference trajectory defined by tðzÞ, ðzÞ, Á L ðzÞ, and the field maps at the same position, these quantities do not have to be computed using Eqs.( 6)-( 8).This allows one to construct a transfer matrix directly from the reference and field data output from any simulation code if desired.

III. TESTING THE TRANSFER MATRIX
To test the validity of our approach, the energy gain and transfer matrix are calculated through three different field setups and compared to direct integration of the equations of motion using a fourth order Runge-Kutta algorithm.The three field setups used are a DC gun with overlapping solenoid, a superconducting radio frequency (SRF) cavity, and a rf gun with a solenoid-immersed cathode.To check that the transfer matrix correctly describes the transverse dynamics in each case, all four transfer matrix elements are compared with Runge-Kutta integration.To do so, the two principle trajectories through each field setup are computed.These trajectories are defined by the initial phase FIG. 2. Comparison of Runge-Kutta integration (blue) and the tracking using the transfer matrix (green) through the DC gun and solenoid fields.FIG. 3. The field map for the SRF injector cavity and the corresponding energy gain.The cavity voltage is 3 MV, the initial kinetic energy is 1 MeV, and the phase is on-crest.space coordinates u 1 ¼ ð1 mm; 0Þ T , and u 2 ¼ ð0; 1 mradÞ T .Table I shows the various settings used in each simulation.
For the first simulation, we use the field maps for the high voltage DC gun and first emittance compensation solenoid of the Cornell ERL injector prototype.Figure 1(a) shows the field maps corresponding to gun voltage and solenoid field strength given in Table I. Figure 2 shows how the energy gain computed from the constant field solution in Eq. ( 6) compares to the energy gain computed using Runge-Kutta integration.The step size used for the constant field solution is 1 mm.The agreement between the two methods demonstrates that the constant field solution works very well for the longitudinal variables, even for low initial kinetic energies.Figure 2 shows the results of tracking the principle trajectories using both the rf transfer matrix with !¼ 0, and Runge-Kutta integration.As with the longitudinal variables, the agreement between both methods of tracking is excellent.In addition to these results, the expression for the electrostatic transfer matrix has been experimentally verified in [3].Next, the two principle trajectories are computed through the field map of the 1.3 GHz Cornell ERL injector SRF cavity.Figure 3(a) shows the on-axis electric field map of the SRF cavity model with a 3 MV cavity voltage.The corresponding energy gain through the cavity, computed using the constant FIG. 4. Comparison of direct integration (blue) and tracking using the transfer matrix (green) of the two principle trajectories.FIG. 5.The fields energy gain for the rf gun setup.The cavity field is scaled so that the cavity voltage is 1 MV, and the phase is set to the on-crest value for a 1 MeV electron.4. As in the electrostatic case, the agreement is very good.Finally, the principle trajectories are computed through the rf gun setup.To simulate an rf gun, the last 1.5 cells of the injector cavity field map are used.The solenoid field is positioned so that the maximum value of the solenoid field occurs at the cathode.In order to make sure that the rf phase is constant over each step, a simple adaptive step size algorithm is included.This algorithm adjusts the step size so that the change in rf phase over the step is less than a user-defined tolerance.Figure 5(a) shows field maps for the rf gun setup.The corresponding energy gain through the gun is shown in Fig. 5(b).The accelerating field is scaled so that the rf cavity voltage is 1 MV and the phase is set for maximum acceleration.The tracking results for the two principle trajectories are shown in Fig. 6.From the figure, it is clear that the transfer matrix works well in the low energy case.The average step size for the simulation was roughly Áz ¼ 0:1 mm.

IV. CONCLUSION
We have derived and tested a new method for calculating the 4 Â 4 transfer matrix through superimposed rf and solenoid fields.The algorithm computes the transfer matrix directly from the field data without computing eigenfunction expansions or numerical derivatives.Comparison with numerical integration demonstrates that this new method works for low energy beams starting from a solenoidimmersed cathode.Additionally, because the algorithm relies on analytic solutions to the equations of motion, it is simple to implement and guarantees the correct value for the determinant of the transfer matrix.One limitation to this approach is the assumption (inherent in the derivation) that the fields display cylindrical symmetry.For many applications this is a reasonable assumption; however, previous work shows that asymmetric focusing from input power couplers may be noticeable when heavy beam loading is present [11].In addition, when tracking ultrarelativistic particles, the algorithm takes steps typically on the order of a few millimeters, and therefore may not be the best choice for computational speed.In this case, one may still choose to use the Rosenzweig-Serafini matrix.Nonetheless, the matrix algorithm given here strikes an appropriate balance between accuracy, speed, and simplicity not previously achieved.
NEW METHOD FOR GENERATING LINEAR TRANSFER . . .Phys.Rev. ST Accel.Beams 15, 024002 (2012) 024002-7 field solution and Runge-Kutta integration, are shown in Fig. 3(b).The results of tracking the two principle trajectories through the cavity, with a fixed step size of 2 mm, are shown in Fig.

TABLE I .
Simulation parameters.The fields and energy for a 500 kV gun voltage, 0.04 T maximum solenoid field setting, and 1 eV initial kinetic energy.NEW METHOD FOR GENERATING LINEAR TRANSFER . ..Phys.Rev.ST Accel.Beams 15, 024002 (2012)