Modeling the fields of magneto-optical devices , including fringe field effects and higher order multipole contributions , with application to charged particle optics

• A submitted manuscript is the author's version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


I. INTRODUCTION
The transport of charged particle beams through particle optical devices, such as beam-transport lines, particle accelerators, and spectrometer equipment, depends strongly on the shape of their electric and magnetic fields.For this reason, much effort has been put into the derivation of analytical expressions for both these fields and the trajectories of charged particles that pass through said fields.A thorough study on axisymmetric electrostatic and magnetostatic lenses has been performed by El-Kareh and El-Kareh [1].Analytical studies on the trajectory equations for magnetic quadrupole lenses and their solutions have been performed by, among others, Smith [2] (with corrections by Lee-Whiting [3]), Lee-Whiting [4], Matsuda and Wollnik [5], and Nakabushi and Matsuo [6].A derivation of these equations using Hamiltonian theory has been given by Hagedoorn et al. [ 7] and de Leeuw et al. [8].
The solution to the trajectory equations for a beam guiding element is generally presented as a function, called the transfer map of the element, that maps the initial location of a charged particle in phase space on its final location.Three different methods for the calculation of such maps can be distinguished.First, the map can be calculated using aberration coefficients (see Grime et al. [9] or Grime and Watt [10] for their definition).This method has been employed in most of the analytical work mentioned above, and has also been used in various computer codes, e.g., TRANSPORT [11] and TURTLE [12].Second, the map can be calculated using differential algebra, a tool developed by Berz [13], who implemented it in his code COSY INFINITY [14].Finally, Lie methods can be used to calculate the map in such a way that it is always symplectic (i.e., preserving the volume occupied in phase space), regardless of the way it is truncated.These methods have been developed by Dragt et al. [15,16], who also implemented them in their code MARYLIE [16,17].
For the calculation of the transfer map of a magnetooptical element, an accurate description of the magnetic field of the element is essential.This is often done by expanding the scalar potential of the field in a Taylor-Fourier series, as given by, among others, Szilagyi [18], for straight elements having their design orbit along the z axis.The accuracy of this description depends critically on the accuracy at which the z-dependent Taylor coefficients in the expansion, the generalized gradients, can be provided, especially for the so-called fringe fields near the ends of the element, where the z dependence of the generalized gradients is the most obvious.Direct calculation of these gradients using the Biot-Savart law or on-axis field gradients is possible but not very accurate; a more accurate analytical method has been derived by Venturini and Dragt [19].
There are several disadvantages to the Taylor-Fourier expansion mentioned above.First, it leads to a culmination of power series terms in the transfer map calculations, which proves difficult to handle.Moreover, the divergence-free and curl-free nature of the magnetic field is not preserved once the power series is truncated [20].Transfer map computation using series expansions in the magnetic field description usually means a tradeoff between accuracy of the result and a reasonable number of terms in the expressions.
In this paper, we formulate a different approach to the description of the fields of particle optical devices.A more detailed treatment of the work presented here can be found in [21].Parts of it were presented in [22] and [23].The aim of this approach is to overcome the above-mentioned difficulties by using a magnetic field description that is not based on power series expansions.In this way, we do not have to deal with an explosively growing number of higher order terms, which is usually the case in higher order perturbative methods.
It should be noted that this paper is intended to give a detailed description of the new method of field description and does not elaborate on its applications, i.e., on the practical implementation to real-life examples of multipole devices.
As a first step, we give a general description of the magnetic field inside a beam guiding element and its harmonic scalar and vector potentials.This description treats each multipole contribution to the field separately, and gives the total field as a superposition of multipole contributions.
With this done, we assume the magnetic field to be known on a cylindrical surface that has the same axis as the element.We apply the above-mentioned field description to this case, in order to express the magnetic field and its potentials within the cylinder entirely in terms of the boundary values of the field on the surface.As we do not use power series approximations here, and are able to take multipole contributions of any order into account, the accuracy of the field description is limited only by the accuracy of the boundary values.This is true for both the central field and the fringe field region of the device, since this method does not discriminate between z-dependent and z-independent fields.
Since the boundary conditions will usually arise from measurements, the next step is to give expressions for the field, given a series of measurements on either a cylindrical surface, or a median plane, taking the discrete nature of these measurements into account.The resulting field description will be shown to be insensitive to statistical noise in the measurements.Moreover, the field will be a superposition of divergence-free and curl-free terms, so it is always divergence free and curl free, regardless of truncation.As an example, the magnetic field of an existing quadrupole will be reconstructed from magnetic field data.It will be shown that only a fraction of the available data is needed to calculate all of the field.
This method of expressing a magnetic field in terms of measurements knows a wide range of possible applications.One example is the reconstruction of the field of beam guiding elements (application to an existing magnetic quadrupole is presented in Sec.III C).Our method can also be used on clusters of elements with overlapping fringe fields, which are, in fact, to be treated as single elements, and can even be extended to obtain the transfer functions for a complete beam line setup or storage ring.There are also numerous applications in low-energy electron optics.
Once the magnetic vector potential has been described in terms of field measurements, we can use it in particle optics.For this reason, we insert this vector potential in the Hamiltonian equations of motion for a charged particle in a magnetic field, and use a finite difference method to solve the system of equations numerically, in such a way that the steps in the finite difference method match the steps between the measurements.As a consequence, we developed a tool to express the transfer function entirely in terms of these measurements.
The methods presented in this paper are equally suitable for the calculation of generalized gradients of any order, directly from field measurements.This provides an elegant way to combine the results of this paper with existing analytical results.
It should be noted that, although all methods have been devised to describe the fields of magneto-optical devices, they are equally applicable to electrostatic optical devices, since, in the absence of free charges or currents in the interior of the device, the electrostatic potential can be expressed in the same mathematical form as the scalar potential of a magneto-optical device.Once the electrostatic scalar potential has been obtained, it can be introduced into the description of the transfer function of the device, much in the same way as the magnetic vector potential.

A. Basic equations
From Maxwell's equations for the static electromagnetic field, a scalar potential u and a vector potential A exist for a magnetic field B in a region without free charges or currents, satisfying The vector potential A will be chosen such that = ?A 0 (Coulomb gauge).Then A is harmonic: D A 0. These equations are applied to the magnetic field of a magnetic multipole device, within a cylindrical surface.Cylindrical coordinates ͑r, w, z͒ are chosen in such a way that the z axis coincides with the central axis of the device.We assume the scalar potential u to be known at the cylindrical surface r R, and introduce dimensionless coordinates r ‫ء‬ r͞R, z ‫ء‬ z͞R.Then the potential problem for the scalar potential u reads (we drop the stars for convenience) 062401-2 062401-2 PRST-AB 4 MODELING THE FIELDS OF MAGNETO-OPTICAL … 062401 (2001) Here, C 1 and C 2 are constants and U͑w, z͒ satisfies both lim z!6`U ͑w, z͒ C 6 and lim z!6`Uz ͑w, z͒ 0 for all w.If the Fourier series expansion with respect to w of the potential is known to contain no solenoidal (windependent) term, we set C 6 0. The above potential problem has a unique solution for u [24]; from this solution we can derive general expressions for both B and A. These expressions can be used to compute B and A directly from magnetic field measurements, without the need to derive u first.

B. Harmonic potentials.
Since we use cylindrical coordinates, we can expand a solution u͑r, w, z͒ of (4) into a Fourier series: u͑r, w, z͒ a 0 ͑r, z͒ 1 Note that the solutions for a m ͑r, z͒ and b m ͑r, z͒ can be expressed in this form because of the boundary conditions for jzj !`in (4).
In the next section, the smooth (i.e., arbitrarily often differentiable) functions A m ͑z͒ and B m ͑z͒ will be determined from the boundary condition u͑1, w, z͒ U͑w, z͒.For convenience, we write J m for J m ͑r d dz ͒ and introduce the dummy coefficient B 0 0. Then the general solution for u͑r, w, z͒ reads As follows from Fourier integral theory, this general solution will automatically obey the boundary conditions for jzj !`.

Next, we determine a harmonic vector potential
A for the magnetic field B using (1).This vector potential will be expressed solely in terms of the coefficients A m and B m that fix the scalar potential u for B. The vector potential is not unique: adding the gradient of any harmonic scalar field to a harmonic vector potential for B yields another harmonic vector potential.As long as = 3 A B there is no physical reason to prefer one choice for A over another.For convenience, we choose (expressed in cylindrical coordinates) From the expressions for the scalar and vector potentials, we find that both u and A are completely determined by the functions A m and B m , which are in turn uniquely determined by the boundary conditions at r 1.
One should note that we used the gauge freedom to derive a harmonic vector potential.In some cases, however, a different (nonharmonic) vector potential is more convenient.For example, the Hamiltonian equations of motion for a charged particle in a magnetic field contain, when expressed in cylindrical coordinates, a large number of terms containing A w .A vector potential A with A w 0 is usually employed here to get rid of these terms.Although this vector potential is not harmonic, its Fourier coefficients can nevertheless be expressed in terms of the coefficients J m A m and J m B m (see Ref. [19]).

C. Introducing boundary conditions
In this section, we will show how to calculate the Fourier coefficients of u, B, and A, for which we derived formal expressions in the previous section, and their various derivatives directly from given boundary values at the cylindrical surface r 1.Such boundary values result, 062401-3 062401-3 e.g., from direct magnetic field measurements, calculations using the Biot-Savart law, or spinning coil measurements.The function U͑w, z͒, from the remaining boundary condition u͑1, w, z͒ U͑w, z͒, can be expanded into a Fourier series: Since lim z!6`u ͑r, w, z͒ C 6 , we find that lim z!6`V0 ͑z͒ C 6 , while lim z!6`Vm ͑z͒ lim z!6`Wm ͑z͒ 0 for m .0. Inserting the general solution (5) into this boundary condition yields the equations

This allows one to express the coefficient
Note that J m ͑ivr͒ i m I m ͑vr͒ with I m the modified Bessel function of the first kind of order m.
The right-hand side of ( 6) is a convolution product where the basic function g m ͑r, z͒ is given by The notation g m ͑r, ?͒ indicates that evaluation of g m for some value of z is postponed until after the convolution product is taken.Results similar to (7) hold for ͑J m B m ͒ ͑1, z͒.
For fixed z , the function g m ͑r, z 2 z͒ is the solution to (6) in the case that V m ͑z͒ d͑z 2 z͒, and, in this sense, it is a fundamental solution.This fundamental property of g m ͑r, z͒ will be used to the full extent in the next section.
The basic result ( 7) is also obtained in Ref. [19], through slightly different methods, and used in a different context.An expansion of (7) into powers of r is derived in this paper in order to obtain power series expansions of the coefficients a m ͑r, z͒ and b m ͑r, z͒.Expressions similar to (7) describing the field of axisymmetric lenses ͑m 0͒ can also be found in Ref. [1].
It can be shown that g m ͑r, z͒ is strictly positive [21,26], so Z Since J m A m depends linearly on V m , this allows one to calculate the effect of errors in V m on the calculation of J m A m .We find, using (8), where d͑J m A m ͒ ͑r, ?͒ and dV m denote the maximum errors in J m A m ͑r, ?͒ (as a function of r) and V m , respectively, on the interval 2`, z , `(see also Ref. [21]).This result will prove useful in the next section, where we show how to obtain V m from measurements.The coefficients J k A m and J k B m with k different from m, that occur in the transverse components of A, and their partial derivatives with respect to r and w can be calculated in the same way as J m A m and J m B m .This will prove useful in Sec.IV B, where the components of A and their derivatives will occur in the Hamiltonian equations of motion for a charged particle in a magnetic field.In general, we write for k $ m 2 1, where It should be noted that the definition of g k m contains no physical parameters, since the input of magnetic field data is performed completely via the boundary values V m and W m .Therefore, the values of g k m can be calculated in advance, which greatly simplifies the calculation of J k A m for a given V m .
As has been done for g m , it can be shown that R 2`j g k m ͑r, z͒j dz , `for 0 # r , 1 and k .m. Analogous to (9), we find that thus relating the accuracy of There are indications that C k ͑r͒ ϳ O͑r k ͒, although this has not yet been proven.For k m 2 1, however, it can be shown that g k m is not integrable on 2`, z , `, so Eq.(11) does not apply to this case.
In the case that the field description given in this section is used in charged particle optics, knowledge of both the partial derivatives (occurring in the trajectory equations) and the integrals (occurring in the solutions to these equations) of J k A m is vital (see also Sec.IV C).Their calculation is rather straightforward.The partial derivatives of J k A m are given by the convolution product of the corresponding derivative of g k m ͑r, ?͒ with V m .Integration of J k A m with respect to z can also be done by integrating either g k m or V m with respect to z. m .0, V m ͑z͒ tends to zero for jzj !`, integration by parts of the right-hand side of (10) yields where Using the basic functions g m m and g m11 m , one can calculate the components of A and their derivatives directly from the boundary values V m and W m , without need to calculate the magnetic field B or the scalar potential u first.In fact, since the functions A m and B m are fully determined by V m and W m , any quantity related to u can be calculated directly from the boundary values, if the appropriate basic function is used.
In practice, boundary values originate either from the potential u (e.g., when produced by a three-dimensional code) or from quantities such as the magnetic field B (e.g., when obtained from measurements).Determining V m and W m in the second case takes a little extra work as compared to the first.As an example, we show how to determine the functions V m and W m from measured values of (a component of) B. We have where the components of B are given in dimensional (unscaled) form.Using the boundary condition , where a r, w, z, we find for m 0, 1, 2, . . .
͒ .These relations, combined with (10), allow one to calculate J k A m and J k B m for 0 # r , 1 and all z if one component of B is known at r 1.Note that, in the case of B z being known, one needs to replace g k m by its primitive G k m , as in (12), while in the case of B r being known, we obtain the correct basic function using (6), where the de-nominator I m ͑v͒ has to be replaced by vI 0 m ͑v͒, resulting in a different basic function.
In general, many field-related quantities can be calculated directly by taking the convolution product of (the Fourier coefficients of) the boundary values at r 1 with an appropriate basic function.This will prove to be a powerful tool.In the next section, we will treat the case of discrete boundary values at r 1 and we will show how to reconstruct these boundary values if the available magnetic field data originate from a surface other than r 1, e.g., a median plane.

III. CALCULATING THE FIELD FROM MEASUREMENTS A. Using measurements at the boundary
In practice, we can measure (the components of) B at r 1 for a discrete set of w and z values and approximate them by interpolating the measurements.Measurements of B w and B z are more convenient than measurements of B r , since the former provide direct approximations of ͑V m , W m ͒ and ͑V 0 m , W 0 m ͒, respectively, while the latter do not.Since piecewise constant or piecewise linear interpolations are almost always employed, the special cases of V m being piecewise constant or linear will be considered.As the measured values for B will be negligible for sufficiently large z, we assume that, for m .0, V m ͑z͒ 0 for z sufficiently large.On the other hand, since lim z!6`V0 ͑z͒ C 6 , a solenoidal contribution to u should be treated by fitting B z ≠u ≠z rather than by fitting u.By using (15), B z can be calculated from V 0 0 , which is considered to be zero for sufficiently large z.
If we assume that V m is piecewise constant, there are pairs ͑l i , z i ͒, with P i l i 0, such that V 0 m ͑z͒ P i l i d͑z 2 z i ͒.By using (12), we find where G k m ͑r, z͒ is defined in the previous section.If, on the other hand, we assume that V m is piecewise linear, then V 0 m is piecewise constant, and V 00 m P i l i d͑z 2 z i ͒, where both P i l i 0 and P i l i z i 0. In this case we have where Gk m ͑r, z͒ Following from the conditions for the l i in the above cases, we find that the terms on the right-hand side of both ( 16) and ( 17) cancel each other for jzj !`, so J k A m tends to zero for jzj !`in both cases.
Since expressions such as ( 16) and ( 17) are derived in a straightforward fashion from Eq. ( 10) for specific instances of V m , the accuracy of J k A m is completely determined by that of V m through Eq. ( 11), and convergence of J k A m with increasing number of z subdivisions is directly related to that of V m .
From the integral expressions in the previous section, we find that related quantities can be derived by replacing the functions G m and Gm by basic functions corresponding to these quantities, while retaining the pairs ͑l i , z i ͒.In fact, these pairs determine the corresponding multipole contribution completely.
At this point, we show how to determine the various multipole contributions to the magnetic field from the values of B z at r 1, as given by (15).Assume B z j r1 has been measured at the points with coordinates ͑w, z͒ ͑w j , w i ͒.Since V 0 m is the Fourier coefficient of cos͑mw͒ in the Fourier series representation of B z , V 0 m ͑w i ͒ is obtained from ("Џ" denotes a numerical approximation) where P j ͑w j11 2 w j ͒ 2p.(Note that the simple midpoint scheme used in this example can always be replaced by a more sophisticated scheme, if so desired.)We choose the piecewise constant approximation for V 0 m ͑z͒, i.e., V 00 m Џ P i l i d͑z 2 z i ͒, with z i Then ͑J k A m ͒ ͑r, z͒ is determined by (17).Since B z ͑z͒ ഠ 0 for jzj sufficiently large, so is V 0 m ͑z͒.Since Gm ͑r, z͒ is a known function, we see that ͑J k A m ͒ ͑r, z͒ can directly be related to the field measurements at r 1.This provides a convenient and flexible way to calculate the magnetic field and related quantities from the field measurements of any desired particle optical device.The nature of the methods outlined in this section also allows the calculation of the field of a cluster of devices, i.e., multipoles, which can in fact be treated as a single device.
It should be noted that one needs measurements performed at 2m different angles w i at least in order to determine the 2m-pole contribution.In other words, the number of different angles at which measurements are taken determines the highest order multipole contributions that can be approximated from a given set of measurements.

B. Using measurements not at the boundary
In the previous sections, the multipole coefficients J k A m were calculated under the assumption that the values of the magnetic field or the magnetic scalar potential at the surface r 1 were available.In many cases, however, the measurements have not been performed at the surface r 1, but at a plane containing the z axis or on the z axis itself instead.We will show that it is possible to obtain the various multipole contributions to the magnetic field in such cases by means of a least squares method, which will be outlined below.
We will treat the case that measurements of B w ͑r, w, z͒ were taken in P ͑P $ 2͒ planes containing the z axis, i.e., at the points ͑r k , w i , w j ͒, k 1, . . ., M, i 1, . . ., P, j 1, . . ., N. In this case, we are able to fit at most M different multipole contributions; in most cases, the 2m-pole contributions corresponding to m 1, . . ., M will be fitted.If less than M multipole coefficients are fitted, the remaining data can be used to improve the statistical properties of the fit.
We assume that 0 , r k , 1, k 1, . . ., M; if this is not the case, then the values r k and z j should be suitably scaled.
As an example, we assume that B w originates from the field of a realistic magnetic quadrupole, where the quadrupole and sextupole contributions are dominant, and other contributions are negligible.In this case, B w is approximated by B w ͑r, w, z͒ 2 2 Rr a 2 ͑r, z͒ sin͑2w͒ 2 3 Rr a 3 ͑r, z͒ sin͑3w͒ .
The coefficients a 2 and a 3 are then approximated by where z l 1 2 ͑w l 1 w l11 ͒.This corresponds to a piecewise constant approximation of the multipole coefficients of B w at r 1.We denote the measured value of B w at ͑r k , w l , w j ͒ by f kij , and define the quantity M͑l, m͒ by In order for B w ͑r, w, z͒ to vanish for jzj !`, l, m must satisfy The optimal values for l, m are then obtained by minimizing M͑l, m͒ under the conditions (18).

062401-6 062401-6
There are a few remarks to be made concerning this example.
(i) In the above example, the assumption was made that the multipole contributions were all normal oriented, i.e., B w did not contain any cos͑mw͒ terms.In this case, it is sufficient to have measurements for one angle w 1 .In general, the orientation of the multipole contributions will be unknown, so B w should contain both sin͑mw͒ and cos͑mw͒ terms, and we need at least two different angles w 1 and w 2 , which are such that mjw 0 2 w 1 j͑͞2p͒ is not an integer for any 2m-pole contribution that is assumed to occur in the magnetic field.
(ii) In the case that the l l and m l are obtained directly from boundary values, as in the previous section, the resulting solution automatically satisfies the boundary conditions for jzj !`in (4).However, this is not true for the least squares method presented in this section; for this reason, the condition (18) has to be imposed explicitly.More generally, if a piecewise nth degree polynomial approximation is employed at the boundary, the conditions have to be imposed.
(iii) In theory, one can determine the various multipole contributions from measurements taken at only two different angles, but in practice one can use measurements at more angles in order to reduce the influence of random errors.The same can be said about the number of different values r k , if there are M different values r 1 , . . ., r M available, and one wishes to fit the coefficients of less than M multipole contributions.As for the different z values, one could use a smaller number of pairs ͑l i , z i ͒ for the fitting process than the number of available values w j , but this reduces the z range where B w ͑1, w, z͒ is assumed to be nonzero.For this reason, one should use all pairs ͑l i , z i ͒ corresponding to the values z i where B w ͑1, w, z͒ is assumed to be nonzero, and use the remaining data (at the z values where statistical noise is assumed to be dominant) to improve on statistics.
(iv) Further improvement on the statistics can be obtained by using a weighted average for the sum of squares in the definition of M͑l, m͒ instead of the unweighted arithmetic mean used in the example.The inverse square of the relative error in the measurement f kij can be used as the weight for the term ͓f kij 2 B w ͑r k , w i , w j ͔͒ 2 .
As shown in the previous section, the pairs ͑l l , z l ͒ fully determine the corresponding Fourier coefficient of any field-related quantity, e.g., u or B w .Knowledge of these pairs allows one therefore to fit the corresponding multipole contribution to any field-related quantity, using the corresponding basic function.The method presented in this section allows one to obtain these pairs from measurements on the z axis or in the plane w w 1 instead of on the surface r 1.This will be very useful in cases where it is not possible to measure on a cylindrical surface.
Comparing the methods developed in this and the previous section, we find that the method of the previous section gives better approximations for the individual multipole coefficients, while the method of this section gives a better overall approximation for the total magnetic field.This effect is caused by the way these methods deal with higher order multipole contributions that are assumed to be zero but, in fact, are not.The method of the previous section yields accurate approximations of the lower order multipole coefficients and completely neglects any higher order multipole contributions that might exist in the field; the method of this section, however, "distributes" the total contribution of higher order terms among the lower order coefficients.Although this decreases the accuracy of the approximation of the lower order coefficients, it improves the accuracy of the overall approximation of the magnetic field.For this reason, the latter method for field approximation might be preferred when using this approximation for calculating a transfer function numerically, as described in Sec.IV B.

C. Experimental test of the presented theory
The theory developed in the previous sections has been verified using actual field measurements for a magnetic quadrupole, performed by Brooijmans [27].Part of the measurements have been used as input for the calculation of the complete quadrupole field, and the outcome has been compared to the remaining measurements.
The quadrupole used for the measurements was normal orientated (i.e., antisymmetric with respect to the planes x 0 and y 0), and the component B y on the plane y 0 was measured, for a number of equidistant x and z values.Since in the center of the quadrupole, B y turned out to depend linearly on x, it has been assumed at first that any higher order multipole contributions could be neglected with respect to the quadrupole field.Therefore, B y has been fitted using basic functions for m 2 only.Since B y B w for x .0 and B y 2B w for x , 0, the basic function Ĝ2 ͑r, z͒ 1 r G 2 ͑r, z͒ 1 pr has been employed.Furthermore, the method of Sec.III A (measurements at the boundary) has been used to fit B y ; the outermost row of measurements provided a piecewise constant approximation of V 2 ͑z͒.This resulted in the following expression for B y at y 0 and 2R , x , R: where B k denotes the measured value for B y at ͑x, z͒ ͑x 1 , z k ͒ and z k11͞2 ͑z k 1 z k11 ͒͞2.Finally, B y has been 062401-7 062401-7 From this figure, we find that there is a good agreement between measured and calculated values of B y .The small differences between calculations and measurements, most obvious around the physical ends of the quadrupole at z 155 and z 455, arise either from parasitic dipole and/or sextupole contributions to the magnetic field or from misalignment of the device's design orbit with respect to the grid followed by the Hall probe.A more thorough field calculation, in which these effects are included from the beginning, will take care of this.
In short, the above example shows that the methods developed here can be used to calculate the complete magnetic field of a multipole device, while only a limited amount of measurements is needed as input.By the same methods, other field-related quantities, such as the vector potential A, can be calculated easily from the same measurements.

IV. APPLICATION TO CHARGED PARTICLE OPTICS
In this section, we derive the Hamiltonian for the motion of a charged particle in a magnetic field, with z as the independent coordinate instead of the time t, and insert the expressions for A obtained in Sec.II C. We outline a numerical method for solving the resulting system of first order nonlinear ordinary differential equations.This method will take both the Hamiltonian nature of the equations and the special nature of the components of A, which are superpositions of shifted basic functions, into account.We will use this method to show that the total transfer function of the beam guiding device under consideration can be expressed completely in terms of the field measurements at r R.
Finally, we show how to incorporate the given field description in existing analytical methods that express the transfer function in terms of aberration coefficients.We outline how these coefficients can be obtained, up to any desired degree, directly from field measurements.

A. Charged particle Hamiltonian
The motion of a particle with mass m 0 and charge Q in the field of a magnetic multipole is governed by the following Hamiltonian: Since the location of the particle in the z direction is well defined, while the location in the t direction is not, it is often desirable to make z the independent variable instead of t.This is possible for all parts of the particle trajectory where dz͞dt fi 0. To obtain a Hamiltonian system for which z is the independent variable, a new canonical momentum p t is introduced, and the Hamiltonian K for this system is chosen such that the action integral remains invariant: This way, the canonical form of the equations of motion is preserved.K͑q k , p k , t, p t , z͒ is then obtained by solving the equation p t 1 H 0 for p z : where the sign of the square root has been chosen such that dt͞dz .0. Note that in the case that ≠H͞≠t 0, the order of the system has decreased by 2 as a result of the exchange.A complete treatise on this method can be found in, among others, Refs.[15,28].
From the Hamiltonian K, the equations of motion for a charged particle in a magnetic field can be derived, with z as the independent variable.The magnetic field description presented in this paper enters the equations through the vector potential A. Methods for solving these equations will be presented below.

B. Calculating transfer functions
In this section, we show how to obtain the transfer function from numerical integration of the equations of motion.The system of equations derived from the Hamiltonian (19) will take the following form: Here, q ͑r, w, t͒, p ͑p r , p w , p z ͒, x ͑r, w, z͒, and f is a given function.The prime (9) denotes total differentiation with respect to z.
Note that if A does not explicitly depend on t, p t is a constant of motion, and the particle trajectories can be calculated without solving the equations for t and p t .
As derived in Sec.II C, the components of A and their partial derivatives all take the form The l i originate from measurements at the points z w i .The z i were defined by z i 1 2 ͑w i 1 w i11 ͒.This provides good approximations for A and its derivatives in the region 0 # r , 1, 2p , w , p, w 0 , z , w n .Now we proceed to solving the system (20) by means of a finite difference method.The discrete version of ( 20) is given by We apply the initial condition ͑q, p͒ ͑w 0 ͒ ͑q 0 , p 0 ͒, and calculate ͑q, We then find ͑q f , p f ͒ ͑q, p͒ ͑w n ͒ by repeated application of (21).Through repeating this procedure for a number of initial locations, and interpolating between the corresponding final locations, the complete transfer function can be obtained.
It should be noted that the description for A used in the differential equations is accurate for all r # 1, and not only for small r, as in the case of a Taylor expansion.Therefore, a transfer function obtained using the above method is accurate for all r # 1, and its use is not limited to the paraxial region.
Since the steps of the finite difference method are located precisely at the points z w i , where the field measurements were performed, we find that there are no interpolation errors in the values for A used in the calculations.By optimizing the interpolation of the boundary 062401-9 062401-9 values at r 1, such that not only the boundary values but also their z derivatives are matched at the points z w i , we can also remove interpolation errors from the values of the partial derivatives of A for better results.A comparison between the method presented in this section and common perturbative methods reveals the following.For small to moderate values of r, a perturbative method produces an accurate result with less effort than is needed when using our method.It is in situations where nonlinear terms dominate the field description where the present method is most useful, e.g., for large values of r, or close to the edge of a multipole device, or in a cluster of devices, where there may be nonlinear coupling of overlapping fields.In such cases, the number of terms needed in a perturbative method, and hence its complexity, increases considerably, while our method can be applied just as easily as to a paraxial case.Also, our method does not have any restrictions on the size of the radius r, where a perturbative method of a given order is usually only applicable at radii smaller than a well-defined maximum.
However, even in cases where a perturbative method is the most logical choice, the field description presented in this paper can be used to obtain the coefficients of the power series terms.This will be demonstrated in the next section.

C. Incorporating the field description into existing results
As mentioned in Sec.I, many methods have been developed to obtain analytical expressions for the transfer function of a given magneto-optical device.They usually start from the Hamiltonian for the motion for a charged particle in the magnetic field of the device; then both the magnetic vector potential and the Hamiltonian are approximated by a truncated power series.The truncated Hamiltonian is then used to derive approximated, but still Hamiltonian, equations of motion, for which solutions are obtained using successive substitution.These solutions are mostly given in the form of power series approximations with respect to the transverse coordinates x and y, and the coefficients of the various terms, the so-called aberration coefficients, are expressed in terms of the generalized gradients of the vector potential.
These generalized gradients are in fact the on-axis radial derivatives of A, which can hardly be obtained from direct measurements, but can, on the other hand, easily be obtained by expanding the basic functions in the description of A into powers of r (see also Ref. [19]): The expansion of g k m into powers of r is then inserted into (10) in order to obtain expansions of J k A m , and finally of A, into powers of r.From ( 16), we find that the z-dependent coefficients in this expansion, i.e., the generalized gradients, are obtained in the same way as the coefficients J k A m .In other words, we can easily derive accurate approximations for the desired generalized gradients, in term of measurements far from the axis, up to any desired order.This provides a very convenient way to incorporate the field descriptions outlined in this paper into existing analytical work.

V. CONCLUSIONS
The magnetic field inside a magneto-optical device, and its harmonic scalar and vector potentials, has been explored in the area 0 # r , R, 2`, z , `.The various multipole contributions to these quantities have been fitted using field measurements at the boundary r R and shifted basic functions.The same set of measurements and shifts can be used to fit various field-related quantities.
All the approximations of the field and its scalar and vector potentials satisfy Maxwell's equations; the approxi-mated potentials are harmonic, which allows one to apply harmonic potential theory whenever necessary.
As the definitions of the basic functions corresponding to the quantities to be fitted do not contain any physical parameters, these functions can be calculated in advance to the desired accuracy, which greatly simplifies the fitting procedure.
The developed procedure is independent of the exact form of the boundary conditions and can be used to fit the field of one device or a cluster of various consecutive devices.
The procedure works for any order multipole contribution, but will be the most useful for lower order multipole contributions, since higher order multipole contributions are more difficult to obtain from measurements, while their effect on particle trajectories will often be small.
Once the vector potential A for the field of a multipole device has been calculated, it can be inserted in the Hamiltonian equations of motion for a charged particle passing through the device.The description of A in terms of z shifts of basic functions can conveniently be combined with the numerical integration of these equations.By expanding A into powers of r, descriptions of the z-dependent generalized gradients in terms of z shifts of basic functions are obtained.These descriptions are much more accurate than descriptions in terms of on-axis field derivative measurements, and are an excellent way of combining the methods outlined in this paper with much of existing analytical work on particle optics.
The methods developed in this paper can also be used for the calculation of the fields and transfer functions of static electro-optical devices, since the electric field of such a device can be written in the same form that was used for the field of a magneto-optical device, and the calculation of the transfer function is similar to the magnetooptical case.
X m1 ͓a m ͑r, z͒ cos͑mw͒ 1 b m ͑r, z͒ sin͑mw͔͒ .The terms corresponding to a certain value of m represent the 2m-pole contribution to u: m 0 corresponds to a solenoid, m 1 to a dipole, m 2 to a quadrupole, etc.For convenience, we adopt the following notation for this Fourier series: u͑r, w, z͒ ͑ ͑ ͑a m ͑r, z͒, b m ͑r, z͒͒ ͒ ͒ m0 , using the dummy coefficient b 0 ͑r, z͒ 0. Having inserted this expansion into the Laplace equation for u, we find second-order partial differential equations for the coefficients a m ͑r, z͒ and b m ͑r, z͒.Their formal solutions are given by a m ͑r, z͒ J m : F 21 ͓J m ͑ivr͒ 3 ͑F A m ͒ ͑v͔͒ ͑z͒ , where J m denotes the Bessel function of the first kind of order m and F denotes the classical Fourier integral transformation with respect to z: ͑F f͒ ͑v͒ : Z 2`f ͑z͒e 2ivz dz .

FIG. 1 .
FIG. 1. (Color) Comparison of the measured and calculated B y as a function of z for various values of x.The black dots represent the measurements, the red curves indicate the calculated field.The top row of measurements served as input for the calculations.