Magnetic ﬁeld distribution of injection chicane dipoles in Spallation Neutron Source accumulator ring

We have performed 3D computing simulations to study the magnetic ﬁeld distribution of the injection chicane dipoles in the SNS accumulator ring. The simulations yield the performance characteristics of the magnets and generate the magnetic ﬁeld data in three dimensional grids for further beam tracking study. Based on the simulation data, a 3D multipole expansion of the chicane dipole ﬁeld, consisting of the generalized gradients and their derivatives, has been made. The harmonic and pseudoharmonic components in the expansion give much insight into the magnet physics and can ﬁt directly into theoretical frame work of beam optics. The expansion is quasianalytical by ﬁtting numeric data into interpolation functions. A 5th-order representation of the magnetic ﬁeld is generated, and the effects of even higher-order terms on the ﬁeld representation are discussed.


I. INTRODUCTION
In the Spallation Neutron Source (SNS) under construction at Oak Ridge National Laboratory (ORNL), a 1 GeV H ÿ beam from its linac is injected into an accumulator ring to produce a proton beam by stripping off electrons in low magnetic field. In the injection region there are four DC magnetic dipoles (D1 to D4) in a chicane structure to produce a horizontal orbit bump, as shown in Fig. 1. The striping foil is located at the left edge of the second dipole (D2), whose field is very critical to determine the stripping efficiency, beam loss, and the trajectories of stripped electrons and circulating beam. The dipoles D2 and D3 are the C-shaped magnets designed with sophisticated pole tips to achieve high field uniformity and good field quality in order to minimize uncontrolled beam loss and to dump stripped electrons to their collectors. The third dipole (D3) is close to D2 and their fringe field overlaps. In addition, D2 and D3 have complementary pole tip structures, resulting in significant cancellation of the integrated harmonics through these two dipoles. In order to assess the effects of the dipoles on the dynamics of the circulating beam and on the striped electron trajectories, it is necessary to know the detailed magnetic field distribution of the chicane dipoles, especially D2 and D3.
The injection chicane dipoles were designed and developed at Brookhaven National Lab (BNL) [1]. Tremendous efforts including 3D computing simulations were devoted to these critical devices. In order to establish magnetic field data files for beam dynamics analysis at ORNL, we initiated an independent simulation study of the D2 and D3 properties. Our simulations also yield performance characteristics of the dipoles and generate magnetic field distribution in 3D grids for beam tracking study. However, the data file occupies excessively large memory, and so man-aging the data in subsequent work becomes difficult. More importantly, discrete 3D data points do not provide any insight into the magnet physics. They do not fit directly into the existing frame work of beam optics [2], where it is usually required that the magnetic field is given by the components of multipoles. Therefore, we go further to expand the chicane field into 3D multipoles based on the simulation data. The 3D multipole expansion describes the magnetic field in terms of harmonic and pseudoharmonic components. The expansion is analytical in nature if we fit a few generalized gradients of numerical data into interpolation functions. The field at any point within a cylindrical volume on the z axis of magnets can be computed by approximate formulas such as a 5th-order representation presented in this paper. Although the 3D multipole expansion is illustrated with the chicane dipoles in this work, the theory and techniques presented here can be equally applied to other accelerator magnets, such as quads, sextupoles, etc.
In Sec. II, we describe briefly 3D simulations of the chicane dipoles D2 and D3. The field distribution and other properties of the magnets are reported and compared with that from the BNL simulation. Section III is devoted to the 3D multipole expansion of magnetic field from the simulation data. The theory is briefly reviewed and the expansion techniques in our work are described. The validity of the expansion is demonstrated by the on-axis generalized gradients and a 5th-order representation in comparison with the simulation data. The effects on the magnetic field representation by even higher-order terms up to the harmonic number m 10 in the expansion are discussed. A short summary of this work is presented in Sec. IV.

II. 3D COMPUTING SIMULATIONS
The simulation environment used in this work is OPERA3D/TOSCA [3]. We first simulate D2 and D3 sepa-rately and then combine the two dipoles together in the final model, as shown in Fig. 2. The simulation models are built by the OPERA package ''MODELLER,'' rather than by ''Pre-Processor'' as was done at BNL. The MODELLER makes it easier to simulate two or more magnets together. The mechanical dimensions of the dipoles are taken from the BNL design drawings and parameters [4]. Both the dipoles D2 and D3 have a nominal gap height of 24.7904 cm, pole width of 50.038 cm, and 14 turns per coil. The nominal pole length is 70.104 cm for D2 and 47.9552 cm for D3. The distance between the D2 gap center and the D3 gap center is 181.4068 cm. The origin of the coordinates in the simulations of the two dipoles together is located at the D2 gap center. In the model, the dipoles D2 and D3 are energized at 2140 and 1690 A, respectively. These currents should be close to the operation values for a 1 GeV proton beam.
The magnetic field on the z axis from the simulation is plotted in Fig. 3. The dominant component is the dipole field B y , which has a minimum at z 93:87 cm between D2 and D3. The axial field B z is also rather large, as plotted by the dashed red curve. This is an indication of significant fringe fields produced by the dipoles. The integrated field and the integrated harmonic contents are analyzed by a rotating Cartesian patch, similar to an unbucked winding in a Halbach search coil in accelerator magnet measurements [5]. Figure 4 shows the integrated harmonic distribution evaluated with a reference radius of 8 cm. The integrated  harmonic amplitudes are expressed in units with one unit as 1E-4 of the integrated dipole field. Note that the vertical axis is in logarithmic scale. The blue and green bars are calculated for the D2 region and D3 region, respectively, with a separation point at z 93:87 cm. The red bars are computed for the entire length of both D2 and D3. It is clear that the overall integrated strength of the quadrupole, sextupole, and octupole terms through the two dipoles is much smaller than that in each individual dipole. This is due to opposite phases of these harmonics in the two dipoles as produced by the complementary pole tips. The stripping foil is located at x 0, y 2 cm and z 30:7 cm (note that the origin of the coordinates is at the center of the D2 gap), which is around the left edge of D2 in Fig. 2. It is important to know the field in the stripping foil area since the motion of stripped electrons is deter-  mined by the magnetic field lines in that region. In Table I we list these parameters and compare the simulation results between ORNL and BNL [6]. The agreement looks good. From the D2 and D3 simulation we have generated the magnetic field distribution in 3D grids. The data file stores the magnetic field within a rectangular bar, which has a cross section of 20 20 cm 2 and a length from z ÿ200 cm to z 400 cm. The step size is 0.5 cm in all dimensions. Thus, the file contains more than 2 10 6 space points, each of which is specified by 3 Cartesian coordinates (x; y; z) and 3 field components (B x , B y , B z ). The file size is about 250 MB. If a step size of 0.25 cm is desired, it would produce a data file of more than 2 GB. This would be very difficult to handle in subsequent applications. More importantly, the existing frame work of beam optics usually requires the expansion of all the physical parameters including the magnetic field into various orders. The discrete 3D magnetic data on grids do not satisfy this requirement. These issues have motivated us to develop the 3D multipole expansion as described below.

A. Review of theory
The fundamental theory of the 3D multipole expansion of the magnetic field in cylindrical coordinates has long been established [7]. Some early work on this topic in accelerator magnet analyses and measurements can be found in Refs. [8][9][10][11]. A more recent and comprehensive treatment of the subject, based on magnetic field data obtained through measurement or computation, was completed by Venturini [12,13]. In this brief review, we follow his approach and limit our discussion to the 3D multipole expansion of the magnetic scalar potential, rather than the vector potential.
It is well known that in a current-free region the magnetic scalar potential m satisfies the Laplace equation In cylindrical coordinates (r; ; z) the solution of Eq. (1) can be expanded in terms of the eigenfunctions of the operator @ 2 =@ 2 as Equation (2a) is usually referred to as the 3D multipole expansion for the magnetic scalar potential. The ''sinelike'' and ''cosinelike'' terms are the normal and skew components. The integer m is the order of the multipoles. For example, m 0 corresponds to a pure solenoid, m 1 to a dipole, and m 2 to a quadrupole, etc. The solenoidal field is described by the cosinelike term only, i.e., 0;c r; z. The physical meaning of the functions b m k and a m k in Eqs. (2b) and (2c) will become clear later. I m is the modified Bessel function of the first kind of order m, which can be expressed by a Taylor expansion With the help of Eq. (3), we can rewrite Eqs. (2b) and (2c) as where s or c, and the functions C 2' m; zare defined as It is easy to verify that C 2'2 m; z d 2 dz 2 C 2' m; z. Therefore, if C 0 m; z C m; z is known, all the expansion coefficients in Eq. (4) can be found by successive differentiation. The functions C m; z are called the gen- eralized gradients. Their derivatives constitute the socalled pseudoharmonics. The magnetic field can be found by the gradients of the magnetic scalar potential:B r m . The convention here is to use a sign. The field components in cylindrical coordinates are Note that the solenoidal field (m 0) does not contribute to the azimuthal field B , although we keep the formulas in the same format. Recall that in a 2D multipole expansion the radial dependence of the transverse field harmonics of order m is r mÿ1 . However, in a 3D multipole expansion this applies only to the terms associated with the generalized gradients C m; z (' 0). There are still many more terms of the pseudoharmonics that vary radially as r 2'mÿ1 (' is not equal to 0). For example, in a quadrupole (m 2) the component C 2 2; r 3 (' 1) is called a pseudooctupole since its azimuthal variation follows a quadrupole field while its radial dependence resembles an octupole. The pseudoharmonics are produced by the fringe field. If the field is uniform in the z direction, then the derivatives of the generalized gradients vanish and the pseudoharmonics disappear, and therefore we recover the 2D multipoles.
According to the uniqueness theorem, the magnetic field in a volume of current-free region is uniquely determined by the boundary conditions. If the magnetic field is known on the surface of a very long cylinder coaxial with the magnet z axis, we should be able to find the field at any point within the cylinder by the method of the 3D multipole expansion. Suppose the radial field B r on a cylindrical surface of radius R is found from simulation or experiment, we first analyze the field via Fourier expansion to obtain where B m R; z and A m R; z are the amplitudes of the normal and skew components of B r R; ; z harmonics of order m. By comparison of B r obtained from Eqs. (2a) and (7), we find (9a) HereB m R; k andÃ m R; k are the Fourier transforms of B m R; z and A m R; z In summary, a 3D multipole expansion of the magnetic field starts with the field data on a cylindrical surface, e.g. B r R; z from simulation in our discussion. We then decompose B r R; z into its Fourier components B m R; z and A m R; z, as indicated in Eq. (7). The next step is to obtaiñ B m R; k andÃ m R; k by the Fourier transforms (10a) and (10b). The inverse Fourier transforms ofB m R; k and A m R; k, weighted by a function k mÿ1 =I 0 m kR, yield the generalized gradients C m; z according to Eq. (9). Finally, the magnetic field can be obtained by Eq. (6).

B. Expansion techniques
We first need to calculate the surface field using OPERA. There are two basic field calculation methods in OPERA3D post processor: nodal interpolation and integration. The former, which is more commonly used, offers fast computation but less accuracy. The latter is slow but should give an order of magnitude improvement in accuracy, since it evaluates the iron magnetizations from the solution and then computes the field by integration of r M r1=r [14].
In the 3D multipole expansion of the D2 and D3 field, we calculate the radial field B r on the surface of a cylinder of radius R 10 cm, starting from z ÿ240 cm and ending at z 420 cm in a step size of 0.25 cm. There are a total of 2641 data points in the z direction. At each z position, B r is computed along a circle of radius R 10 cm at a step of 0:5 , which yields 720 data points for the Fourier decomposition. The B r harmonic components up to m 10 (20th pole) are computed. Figure 5 shows the dominant term B 1 versus z, which is the normal term for m 1. To some extent, B 1 resembles the dipole field component B y on the z axis as shown in Fig. 3. However, it contains highly nonlinear terms, i.e., the pseudoharmonics, since it is evaluated at R 10 cm. The other two important components of the B r harmonics are the skew term A 1 for m 1 and the solenoidal term A 0 as plotted in Figs. 6 and 7. The former resembles B x on the z axis while the latter has connection with B z on the z axis, which will be made clearer in the next subsection. The next step is to perform Fourier transforms of B m R; z and A m R; z to arrive atB m; R; k and A m; R; k as described in Eq. (10). Because of their complex nature, the Fourier transforms produce two terms for each harmonic component of B r R; ; z: one associated with the sine term and the other with the cosine term. The wave number k is chosen from 0 to 2 (1=cm) in 4000 steps. The transforms yield spectrumlike distributions of the fields. The dominant termsB 1;s R; k andB 1;c R; k are plotted in Fig. 8 and the other two componentsÃ 1; R; k andÃ 0; R; k are shown in Figs. 9 and 10. Having done this, the generalized gradients can be computed using the inverse Fourier transforms ofB m; R; k andÃ m; R; k, weighted by a function k mÿ1 =I m 0 kR as expressed in Eq. (9). Note that Eq. (9) has a singularity at k 0 for m 0, which is the case of a solenoidal field. There are three options to avoid this problem. First, we can make k sufficiently close but not equal to zero. The magnitude of the resulting C 0 z C 0;c z depends on the minimum k chosen, but its first derivative C 0 0 z remains approximately the same. As will be explained later, it is C 0 0 z rather than C 0 z that has physical meaning. Second, we can calculate C 0 00 z in the inverse Fourier transform and then to get C 0 0 z via numerical integration. This will be described in the next paragraph. Some errors will still occur in numerical integration. The third option is an easy method in practice and will be described in Sec. III C.
With the generalized gradients calculated above, the remaining work is to perform successive differentiation to obtain the desired pseudoharmonics. This has been done by a MATHEMATICA program [15], which first fits the numerical data of C m; z into interpolation functions of a certain order. Although numerical differentiation in MATHEMATICA is straightforward, numerical errors are always introduced and numerical oscillations from higherorder differentiation are often large enough to invalidate the results. A way to remedy this problem is to calculate the even derivatives C m; 2' z according to Similarly, we can obtain the odd derivatives C m; 2'1 z by This method is very effective in eliminating numerical oscillations and has been proven to work well for all the trials up to 2' 10 in this practice. The second option to calculate C 0 0 z in the previous paragraph follows this technique.

C. On-axis gradients
The generalized gradients C 0 0 z, C 1;s z, and C 1;c z have special physical meanings. They express the on-axis magnetic field components B z z, B y z, and B x z at x y 0, respectively. Since these gradients are obtained via a number of mathematical manipulations, it would be interesting to know if they agree with the on-axis field data directly computed from the OPERA3D simulation. This is a way to verify the theory and check all the procedures of the mathematical calculations. In Fig. 11(a) we compare the generalized gradient C 1;s z with the field data B y z at x y 0 obtained directly from OPERA3D. The agreement looks remarkably good. In order to see the discrepancies of the two sets of data more clearly, we plot their differences in Fig. 11(b). The maximum difference is about 0.017 G out of a total field range of 3009.36 G. This corresponds to a relative error in 10 ÿ6 level, which is attributed to numerical errors. Indeed, the generalized gradient C 1;s z is an accurate representation of the dipole field B y z on the z axis. Figure 12 compares the generalized gradient C 0 0 z and the OPERA3D data B z z at x y 0, and Fig. 13 compares the generalized gradient C 1;c z and the OPERA3D data B x z at x y 0. Again, the agreements are excellent. In principle, the 3D multipole expansion is insensitive to the raw data errors produced in simulations or experiments. This is due to a basic property of the solutions of the Laplace equation: the value of the magnetic scalar potential m at some interior point is an appropriately weighted average of its values over any surrounding boundary.
The above results suggest that we can in fact obtain the generalized gradients C 0 0 z, C 1;s z, and C 1;c z by the on-axis field components B z z, B y z, and B x z at x y 0 directly from the 3D simulation data rather than using complicated expansion procedures. This method is easier, faster, and more accurate. This is especially true for C 0 0 z since the other two options involve more numerical errors.

D. A 5th-order representation
To calculate the field at any point off the z axis within the cylinder, Eq. (6) has to be used. In practice, a cut to a certain order of the azimuthal harmonic number m and the radial power number l has to be made. A 5th-order representation of the magnetic field covers the multipoles up to the regular dodecapoles and pseudododecapoles. In this case, a total of 13 generalized gradients are required and they can be stored in the computer as interpolation functions for subsequent calculations. In this way much less memory is required and the field calculations become quasianalytical. The explicit expressions for this representation are given below: We verify these equations by comparing their results with the field data obtained directly in the OPERA3D simulation. Figures 14 and 15 plot the field components B y versus y at x z 0, and B x versus x at y z 0, respectively. The agreement is very good. In Figs. 16 and 17 we plot B y versus x at y z 0, and B x versus y at x z 0. Again, we see good agreement in general between the 5th-order representation and the simulation data, especially in the near-axis region. The deviation between the two becomes large far away from the axis. This is due to the cut off of higher-order terms. A 9th-order (up to 20th pole) approximation can yield much better agreement, as shown by the green curves in Figs. 16 and 17. This higher-order effect will be explained in more detail in the next subsection.

E. Higher-order effects
Figrues 16 and 17 reveal that the results calculated from the 5th-order representation deviate more from the simulation data when space points are far away from the magnet axis. This is caused by cutting off higher-order terms (m > 6 and beyond r 5 ) and can be corrected. In order to see this more clearly, we write down explicitly a 9th-order formula for B y at y z 0 versus x when x is on the positive axis: The green curve for x > 0 in Fig. 16 is generated by Eq. (13), which agrees with the simulation data very well. It would be interesting to know the contributions to B y from each harmonic term up to m 10. This is listed in Table II for point A, which is indicated by a cross at x 10 cm in Fig. 16. The summations on the rightmost column give the contributions from different orders in r, while the summations on the bottom rows give the contributions from different harmonic number m for three different orders in r. A full 9th-order representation yields B y 3023:23 G at x 10 cm and y z 0. As compared to the simulation data of 3023.13 G at the same point, the difference is only 30 ppm. Therefore, we have demonstrated that we can approximate the simulation data to very high accuracy by higher and higher-order representations in the expansion.
In summary, the results from the 3D multipole expansion are very good and achieve the level of accuracy desired. In general, a 5th-order approximation can make a fairly accurate analytical representation of the real field distribution. The accuracy can be improved even more for the field points far away from the axis by adding higherorder terms of m > 6.

IV. SUMMARY
3D computing simulations with OPERA3D/TOSCA have been performed to study the magnetic field of the injection chicane dipoles D2 and D3 in the SNS ring. The simulations have produced the field distribution in 3D grids that can be used for further beam tracking study. The simulations have also provided other useful information on the performance characteristics of the magnets.
Based on the simulated radial field on the surface of a long cylinder coaxial with the magnet z axis, the magnetic field inside the cylinder has been expanded into the 3D multipoles, consisting of the generalized gradients and their derivatives. The 3D multipole expansion has given much insight into the magnet physics in terms of harmonics and pseudoharmonics. The 3D multipole expansion is analytical in nature when one fits the generalized gradients of numerical data into interpolation functions. This feature will prove to be very useful in further beam dynamics study. A 5th-order representation of the chicane dipole field from the 3D multipole expansion has been presented. The effects of even higher-order terms on the field are discussed. Overall, the agreements between the expansion data and the simulation results are very good. The theory and techniques presented in this paper can be equally applied to experimental data on a cylindrical surface, as well as to other accelerator magnets, such as quads, sextupoles, octupoles, etc.

SNS is managed by UT-Battelle, LLC, under Contract
No. DE-AC05-00OR22725 for the U.S. Department of Energy. SNS is a partnership of six national laboratories: Argonne, Brookhaven, Jefferson, Lawrence Berkeley, Los Alamos, and Oak Ridge. The author would like to thank W. Z. Meng of BNL, who performed the 3D simulation of D2 and D3 a few years ago and provided us with his OPERA input file as a reference. The help from J. Simkin of Vector Field, England, in many details of the OPERA code and simulation techniques is highly appreciated. We also would like to thank M. Venturini of LBNL, the author of Refs. [12,13], for his help and comments on this work.