Medical therapy and imaging fixed-field alternating-gradient accelerator with realistic magnets

NORMA is a design for a normal-conducting race track fixed-field alternating-gradient accelerator (FFAG) for protons from 50 to 350 MeV. In this article we show the development from an idealised lattice to a design implemented with field maps from rigorous two-dimensional (2D) and three-dimensional (3D) FEM magnet modelling. We show that whilst the fields from a 2D model may reproduce the idealised field to a close approximation, adjustments must be made to the lattice to account for differences brought about by the 3D model and fringe fields and full 3D models. Implementing these lattice corrections we recover the required properties of small tune shift with energy and a sufficiently-large dynamic aperture. The main result is an iterative design method to produce the first realistic design for a proton therapy accelerator that can rapidly deliver protons for both treatment and for imaging at up to 350 MeV. The first iteration is performed explicitly and described in detail in the text.


I. INTRODUCTION
A. Proton therapy and the need for imaging Modern radiotherapy is a mainstay of cancer treatment today, and in the developed world around half of cancer patients will receive radiotherapy as part of their treatment. The bulk of radiotherapy treatments are given by using x-rays from a low-energy linac (c. 10-20 MeV kinetic energy), where the method of intensity-modulated radiotherapy (IMRT) utilizes numerous treatment field directions and a moveable set of multi-leaf collimators that may obtain a highly-conformal dose to a prescribed treatment volume [1][2][3][4][5][6][7]. Whilst accurate pre-treatment patient imaging (such as computed tomography) are of course crucial to creating and delivering the planned treatment, the inherent near-exponential reduction of xray intensity with depth makes x-ray treatment comparatively less sensitive to errors in the patient density that is inferred from the patient imaging process.
Proton therapy is an alternative method of delivering a radiotherapeutic treatment to a patient, already known since 1947 [8] to potentially offer an inherently more precise delivered dose since the proton stopping transfers those particles' kinetic energy to deposited dose according to the Bethe-Bloch equation [9]. A well-known initial proton energy coupled with a well-known tissue density allows one to place the maximum dose -at the Bragg peak at the end of the particle range -at a desired position within the patient. The overlap of numerous Bragg peaks from protons with differing initial energy allows a conformal dose in the treatment volume whilst potentially better sparing the surrounding tissue from un- * sam.tygier@manchester.ac.uk † kiril.marinov@stfc.ac.uk ‡ robert.appleby@manchester.ac.uk wanted dose, particularly important for nearby organs at risk. That said, proton therapy is not thought to be advantageous over x-ray therapy for all radiotherapy treatments, and is often prescribed for particular complex or pediatric treatments where the unwanted ancillary dose may cause a later induction of secondary cancers as a treatment side-effect. The treatment advantages of proton therapy have led to the creation of the more than fifty operating centers around the world today [10,11]. Supplanting early facilities at research laboratories, today's hospital-based centers predominantly utilize cyclotrons although synchrotrons are also used. Modern cyclotrons -particularly superconducting ones -offer a number of advantages in terms of simplicity, capital cost and possible dose rate at the patient; 1 Gy may be accurately delivered to a patient by such a source in less than a minute using an average current of <1 nA [9]. However, higher-intensity cyclotrons are typically limited to around 230-250 MeV kinetic energy due to relativistic effects, and their fixedenergy extraction requires the use of a mechanical degrader (typically graphite) to lower the proton energy for shallower proton dose delivery. 230 MeV protons are sufficient for treating adult patients (this energy corresponds to about a 33 cm proton range in water), but some treatments would benefit from more rapid (and therefore finer) variation of the energy than is readily achieved using cyclotrons. The latter benefit has led to the proposal both of linacs and of FFAGs (fixed-field alternatinggradient accelerators) for particle therapy; each offers the possibility of obtaining a rapid change of the delivered proton energy at rates up to 1 kHz. Various FFAG designs have therefore been put forward, a notable example being the French RACCAM study that showed how such an accelerator could operate with multiple treatment rooms [12].
Irrespective of the accelerator technology used to deliver the particles, the benefits of proton therapy are today limited by inadequacies in the patient imaging used to estimate the patient tissue density [9]. This determination is crucial to deriving a suitable treatment plan, and is further complicated by the complex scattering effects around inhomogeneities that demands Monte-Carlo dose estimation [13]. X-ray computed tomography is seen as somewhat inadequate as the patient density and composition derived from the Hounsfield measurement can lead to several millimeters of error in the resulting proton range [14][15][16]. Proton tomography is a better, more direct, measure of the desired proton stopping power and therefore several groups are developing clinical proton tomography instruments that typically track millions of individual protons to assemble a three-dimensional image with a resolution approaching a millimeter [17,18]. But, proton tomography requires protons of sufficient energy to pass through the part of the patient to be imaged, implying significantly-higher incident energies than those that would be used to deliver a (stopped beam) treatment in that same volume. A proton source of 250 MeV could be used for imaging through smaller thicknesses, but patients requiring treatment with 230 MeV protons of course require imaging with much higher proton energies -perhaps as high as 330 MeV or more depending on the image resolution required. No commercial cyclotron today offers this higher energy, and whilst linacs and synchrotrons both in principle could offer such an energy only ProTom has offered a 330 MeV system commercially [19]. This lack has previously motivated both the PAMELA design study [20] -which examined a combined accelerator system offering both protons and carbon ions -and our NORMA study [21,22]. Both designs aimed at providing proton energies suitable for tomography, but in the latter proton-only NORMA study the additional goal has been to offer a simple, robust design.

B. FFAGs for proton therapy and imaging
In our previous work on NORMA we concluded that an FFAG design utilizing an FDF (focusing-defocusingfocusing) arrangement of normal-conducting gradient dipoles was the best way to achieve our desired 350 MeV. The maximum field of around 1.6 T necessitates a somewhat larger circumference than could be obtained with superconducting magnets, but allows for an easier control of the tune during proton bunch acceleration using the higher-order field components in the magnets. Injection of a single proton bunch would be obtained from a cyclotron at an energy of at least 50 MeV, and both injection and extraction would be via a conventional pulsed kicker/septum combination that may benefit from lengthening two of the straight sections between the 10 FDF cells to give a racetrack layout; the acceleration cycle of around 1 ms per extracted single bunch would enable rapid bunch-by-bunch variation of the delivered energy at the patient. A more detailed discussion of the accelerator magnet lattice design is given in Garland et al. [21], and is summarized in section II. In this previ-ous article both the round and racetrack configurations of NORMA are considered; however, here we focus on just the round variant.
The optics design of NORMA presented in [21] was performed with idealized magnet models, to enable the overall beam dynamics to be studied and optimized. This is a common approach that allows optimization using a reduced number of variables, and gives a tractable simulation time. Our idealized magnet models use analytic expressions for the magnetic field within the body of the magnet, and a simple analytic expression for the fringe field fall-off; this is a compromise between simplicity and realism. However, once a design has been obtained with idealized magnets it is of course important to check that the approximations used do not have a significant effect on the dynamics of the accelerator. For example, the presence, size and shape of the realistic fringe fields will affect the focusing of a magnet and the fields from the physical magnets will not completely match analytic models. Their shape can introduce higher-order effects that are not expressed in their ideal analytic form and manufacturing tolerances will cause deviations from the ideal field. This is an issue that must in particular be addressed in FFAG design -where the magnets are large in aperture and inherently nonlinear in nature -before one can be confident that a realistic and therefore buildable design has been obtained.
In this article we build on our previous work [21], and discuss the detailed 2D and 3D magnet modelling carried out to improve the realism of the NORMA design. Realistic magnet fields introduce perturbations on the idealized dynamics such as tune shift and a reduction in the dynamic aperture. We show how this can be mitigated -by re-matching and re-optimizing -in order to recover acceptable dynamical properties; we demonstrate the recovery of a sufficiently flat tune profile and sufficient dynamic aperture for injection and acceleration. We introduce the magnet models in three stages, so as to methodically understand the importance of different effects.
The layout of this paper is as follows. In section II we describe the original NORMA lattice and the methods used for modelling complex elements and performing long-term stability studies. This is followed in section III by the magnet design process. In section IV NORMA is modelled using the radial profile from the 2D magnet simulations, with realistic transverse fields but with a simplified fringe model; the effects on the dynamical properties are shown to be small. Then in section V we show the effects of changing the parameters used for the fringe fields, and how these can be mitigated by rematching the overall lattice; the results show that all realistic likely fringe field perturbations to the dynamical properties can be absorbed by retuning the strength and field index of the magnets. Finally, in section VI we model NORMA with a full set of 3D magnet models. Here the effects on the dynamics are more significant, and so we show the necessary corrections to the field profile to re-store the dynamical properties and crucially of the dynamic aperture. The overall result is a method for the iterative design of realistic magnets for a medical FFAG capable of delivering 350 MeV protons for imaging, and a demonstration of the first step of the iteration.

II. THE NORMA ACCELERATOR
In this section we introduce the nominal NORMA design, as described in detail in [21], as well as the code Zgoubi [23] which is used for tracking particles through the lattice.
A. The NORMA lattice The round NORMA lattice design variant is a FFAG consisting of 10 identical FDF triplets. It accelerates proton bunches to 350 MeV, with injection from a cyclotron with at least 50 MeV kinetic energy. NORMA utilizes normal-conducting sector magnets with a scaling FFAG field achieved by pole-face shaping. Within the 36 • cell, each magnet has a sector angle of 6 • ; within the FDF triplet the magnets are spaced by 1.8 • . The scaling FFAG magnets result in a flat tune (tune shift below 10 −3 ) over the full energy range from 50 to 350 MeV. The magnet field fall-off is modelled using Enge-like fringe fields with a 6 cm extent [24,25].
A NORMA triplet cell is shown in Fig. 1. The full ring composed of 10 cells is shown in Fig. 2. Between each triplet is 2.4 m of magnet-free drift space. The parameters for the round NORMA lattice are given in table I. In this article we refer to this lattice as the nominal NORMA design and use it as the baseline for further study. Note that -during the early stages of this workwith idealized and 2D magnets we considered an energy range down to 30 MeV. Some of the dynamic aperture simulations in sections IV and V were carried out at 30 MeV, however this is always a tougher requirement than at higher energy due to adiabatic damping.

B. Tracking simulations
Zgoubi, which was used for tracking simulations, is a charged-particle tracking code widely used for designing and studying FFAGs. It uses a stepwise ray-tracing method; the particle is propagated in small steps and at each of these the magnetic field and its derivatives are evaluated. This gives accurate results for particles over a wide range of momenta and trajectories when moving through large-aperture magnets. Zgoubi features a range of magnet descriptions that together are capable of simulating the complex magnets which are typical of FFAGs. We use the PyZgoubi [22,26] framework around Zgoubi to expand its capability, allowing advanced scripting and optimization.
Zgoubi features both analytic and field-map-based magnet descriptions. In each, the magnet description is used to calculate the field and its derivatives at each in-tegration step along the charged-particle's trajectory. In our studies we have used several magnet models: FFAG, an analytic model of a sector scaling FFAG magnet; DIPOLES, a sector dipole with optional higher multipole fields; POLARMESH, a 2D polar mid-plane field map with an out-of-plane expansion.

Analytic Scaling FFAG
Zgoubi offers an idealized scaling FFAG sector magnet, which we use as the reference model for the nominal error-free lattice. The field is composed of the product of the radial scaling law where r 0 is the reference radius and k is the field index, and a longitudinal fringe function where s is the distance from the effective field boundary, λ is the fringe field extent and C i are the well-known Enge coefficients [24]. In the nominal design we use C 1 = 2.24 and λ = 4 cm with other coefficients set to zero. Zgoubi allows any field overlap between magnets to be modelled by assuming linear superposition. We use the FFAG element for initial optimization of the lattice; however it is limited in its flexibility for error studies as only the position and strength can be adjusted.

Multipole expansion
Zgoubi's DIPOLES element can be used to model a sector dipole, but also allows additional multipole components in the field. This can be used to describe combinedfunction magnets, but also to approximate more complex fields such as a scaling FFAG. It allows the radial field profile to be expressed as where B i are the multipole components. Again this is multiplied by the longitudinal fringe field given in Eq. 2 to give the mid-plane field. There are two methods to find the multipole coefficients to use with the DIPOLES element. The ideal FFAG scaling field can be Taylor expanded about a given radius, usually halfway between the minimum and maximum orbits. Alternatively a multipole expansion can be fitted over the required good field region. A Taylor expansion of B(r) around R 0 gives which can be equated with Eq. 3 to obtain the coefficients. Whilst the Taylor expansion gives the correct field and derivatives about the expansion point, using a fit gives a lower maximum deviation in field for any given expansion order, as illustrated in Fig. 3. Use of DIPOLES with a fit up to 9 th order gives a very good agreement to the dynamics of the FFAG element: the mean closed orbit deviation is below 10 −9 cm and the mean fractional deviation of the tune is below 10 −10 horizontally and below 10 −5 vertically.

Midplane field maps
Zgoubi's POLARMESH element allows a magnet to be defined in terms of a 2D polar field map in the mid-plane, with field maps generated in external magnet simulation codes such as OPERA [27]. This allows simulation of the deviations from an ideal field that are likely to occur in a real magnet. It is also possible to generate field maps representing the ideal field by evaluating the analytic equations for a scaling FFAG magnet at the grid points. This can be used to distinguish simulation effects due to an interpolation step from actual effects due to the imperfections of a realistic magnet. It also allows investigation into the resolution requirements for the field map. We find that for the nominal NORMA lattice a 1000×1000 grid of mesh points per half cell (from the start to the center of the triplet) is sufficient to give a mean closed orbit deviation of 10 −6 cm and mean fractional tune deviations of 10 −5 horizontally and 10 −4 vertically. Increasing the mesh density to 2000×2000 does not significantly improve agreement with the analytic models.

C. Dynamic aperture
Dynamic aperture (DA) is a measure of the stable area of phase-space for particles circulating in an accelerator. A particle within the stable area will survive a large number of turns though the accelerator lattice, while a particle outside this region will be lost after a small number of turns. In practice a finite number of turns must be simulated in order to predict the stability. For NORMA we consider a coordinate stable if a particle starting there will survive 1000 turns, as this is representative of the length of the acceleration cycle of around 1 ms. We use a strict definition for dynamic aperture where -for any given amplitude -a set of particles with a range of phase-space angles are tested and must all be stable, as described in [22].
A large dynamic aperture is required to transport the injected bunches through the accelerator with a low loss. A large dynamic aperture increases overall transmission efficiency (from injection to extraction) and therefore reduces radiation and activation. The injected bunches from the cyclotron will have a typical normalized emittance of less than 10 mm mrad [28]. We therefore require that the normalized dynamic aperture is kept around 50 mm mrad or greater, a specification considered sufficient for this application [20].
It is useful to see the effect on the dynamic aperture of the nominal lattice design from using the differing tracking methods, before studying modifications to the lattice. This allows us to distinguish DA changes that are due to using a different field from those due to the choice of tracking method. Figure 4 shows the dynamic aperture over a range of real space angles in x and y, i.e. where 0 • is the horizontal and 90 • is the vertical dynamic aperture, using FFAG, DIPOLES and POLARMESH elements. As before, the DIPOLES shows good agreement to the FFAG element. However, with POLARMESH we see a significant reduction in DA. Increasing the number of mesh points used in the POLARMESH does not improve agreement with the analytic models. We believe that this is due to small errors from the interpolation of the field map building up when tracking for a large number of turns. Therefore, we consider only FFAG and DIPOLES elements to be suitable for DA calculations.

D. Field Errors
Deviations of a given field from an ideal field can be measured in a number of ways, and it is common to specify maximum deviations from the ideal field or gradient. In general the deviation of each multipole component from the ideal field can be measured. In a synchrotron a multipole expansion about the magnet center can be used, but for a magnet that accepts a wide range of orbits one needs to be careful about where the multipole components are measured. This is important because introducing a multipole error of a given order at one lo- cation will change the lower-order multipole values everywhere else; for example a sextupole error at a given orbit radius causes a quadrupole and dipole shift across the magnet. For a field map defined by field strengths on a regular grid, the multipole components can be found either by fitting a polynomial to the whole map or a subsection of it, or by repeated numerical differentiation (e.g. the quadrupole component is proportional to the 1 st derivative of the field with respect to the radius).

III. MAGNET DESIGN AND OPTIMIZATION
This section describes the finite element models of the lattice magnets, the strategy chosen for their optimization and the main steps of its implementation.

A. 2D magnet models
The required nominal radial field profile B N (r) is given by Eq. 1. The values of the parameters r 0 , B 0 , k and the extent of the good-field region required are specified for each magnet. The first step of the magnet design and optimization process is to fit Eq. 1 to a polynomial within the good field region. By following the standard procedure (see e.g. [29]) the coefficients of the fitting polynomial can be used to obtain the two-dimensional scalar potential Ψ = Ψ (x, y) such that B = ∇Ψ in the air gap of the two magnets. Figure 5 shows the result. Note, that the 2D magnet modelling and optimization is performed in a coordinate system with an origin located at the center of the good-field region of each magnet and the x-axis points along the machine radius.
If the permeability of the magnet yoke material is infinitely high the lines of constant scalar potential in Fig. 5 coincide with the faces of the magnet poles that generate FIG. 5: Lines of constant scalar potential for the D-magnet (blue). The good-field region required is represented by the rectangle (red) and the origin of the coordinate system is the center of the good-field region.
The vectors (green) are proportional to the local gradient of the scalar potential.
the required field. The good-field region (plus any reasonable contingencies) must fit within the air gap of the magnet and this determines uniquely the actual magnet pole shapes and the nominal value of the scalar potential that the magnet will be operated at. Indeed, as Fig. 5 shows, each pole configuration is uniquely determined by the absolute value of its scalar potential |Ψ 0 |. Finally, the nominal magnet current I N = 2 |Ψ 0 |/µ 0 (Ampere-turns) can be obtained from Ampère's law as usual, where µ 0 is the vacuum permeability.
However, in reality the permeability of ferromagnetic steel normally employed in accelerator magnet construction is finite. Therefore, strictly speaking the magnet pole shape obtained by means of the multipole expansion technique described in the preceding paragraph is an approximation only. In addition, the onset of magnetic steel saturation lowers the permeability even further and this leads to field errors that cannot be neglected, particularly in the high-field regions of the magnets. If the field is the actual field generated by the magnet, is too high then the pole shape y = y(x) must be replaced with a new pole shapeỹ = y(x) + ∆y(x) such that the magnitude of the field error B N − B M it generates is acceptable. The approximate relationship between ∆B and ∆y(x) can be obtained as follows.
where ∆y (x) is the correction to the pole profile. In Eq. 5 terms of the order of ∆y 3 and higher have been neglected. An approximate solution to Eq. 5 is easy to obtain and the result is In Eq. 6 terms of the order of (B N − B M ) 3 and higher have been neglected. Equivalently Eq. 6 can be derived by formally inverting the function B M = B M (y) and Taylor-expanding the result y = y (B).
It is clear that dB M dx must not be zero as it is present in the denominators of both the expressions for α 1 and α 2 . However, in deriving Eq. 6 it has been assumed that the term of the order of (B N − B M ) 2 is a small perturbation compared to the term of the order of (B N − B M ), (or equivalently, the term proportional to α 2 in Eq. 6 is much smaller than one) and so the series can be truncated without any significant loss of accuracy. This assumption may not be valid if the value of dB M dx is sufficiently small. This means that in the vicinity of local extrema of B M Eq. 6 becomes inaccurate. In addition dy dx must be nonzero in order to ensure the existence of the inverse pole profile function x = x(y). If however, dy dx → ∞ (i.e. dy dx has a pole (see Fig. 5) then in the vicinity of that pole Eq. 6 may no longer be small compared to one. Hence, the applicability of Eq. 6 relies upon the following two conditions: (i) dy dx is non-zero and not too large and (ii) is not too small. A magnet pole optimization procedure based on Eq. 6 must be applied iteratively as Eq. 6 itself is an approximation. The zeroth-order approximation to the pole profile y 0 (x) is obtained by means of the multipole expansion technique and the field distribution B  M (x) and y 0 (x) are then substituted into the right-hand side of Eq. 6 to obtain the corrected pole profile y 1 = y 0 +∆y 0 which is in turn used to obtain the field distribution B M (x) and so on. The implementation of this scheme on a computer uses MATLAB [30] to calculate the pole profile correction according to Eq. 6 from field distribution data generated by OPERA 2D and the result is then fed to OPERA 2D, which in turn calculates the updated field distribution. This process is repeated until the desired maximum field error is reached. Figure 6 shows the result from a test run. The initial FIG. 6: A test of the pole optimization procedure performed on the 2D F-magnet model. Distribution of the relative field error obtained from three different iterations. The inset shows the 2D flux density distribution as obtained with OPERA 2D illustrating the low-field and high-field regions of the magnets. Red color corresponds to high flux density.
(zeroth-order) pole profile and the current strength were obtained with the multipole expansion technique. In order to test the optimization algorithm the coil current was intentionally increased so that the field in the main part of the magnet (which operates far from saturation) is 2% higher than the nominal field. Within 7 iterations the peak relative field error was reduced by a factor of over 200: from 2% to less than 10 −4 . Figure 7 shows the good match for the dipole, quadrupole and sextupole components between the 2D field profiles and the ideal nominal magnet field. The first four iterations were performed in the region −0.3 m < x < 0.3 m and the remaining three in the region −0.25 m < x < 0.25 m while the specified good field region for this test was −0.23 m < x < 0.23 m. It was noted that the convergence in the central region of the magnet was faster than in the low-field and high-field regions. Indeed, as Fig. 6 shows, the extent to which the error is reduced by the first iteration is greater in the central part of the magnet than towards the magnet ends. A possible explanation for this is that in the lowfield region dy dx gradually increases towards the pole rolloff area whilst in the high-field region dB M dx approaches zero near the peak of the magnetic field. As pointed out earlier these are precisely the conditions when Eq. 6 is no longer accurate. Figure 8 shows the correction to the pole profile obtained after 7 iterations. As can be seen, ∆y (x) ≥ 0 in the region −0.23 m < x < 0.23 m, which corresponds to increasing the magnet gap and decreasing the field strength in full accordance with Fig. 6. The notch formed near x = 0.3 m is most a likely a result of Eq. 6 becoming less accurate close to the region where dB M dx → 0.

B. 3D magnet models
The 3D magnet models were created by extruding the 2D magnet shapes, Fig. 9 shows the OPERA 3D model of the focusing (F) magnet. Two clamping plates were added on each side of the magnets to restrict the extent of their fringe fields and eliminate possible cross-talk between adjacent magnets. In order to save time the overall height of each magnet, the size of its back leg, the thickness of the clamping plates and the geometry of its roll-off in the high-field region were not optimized. The optimization in 3D was performed on the integrated field strength only and no attempt was made to adjust the main field (generated in the air gap of the magnet) and the fringe field (generated outside the air gap) separately. The integrated field strength I(r) is defined as where the origin of the cylindrical coordinate system is located in the machine center, and r, ϕ and y are the radial, azimuthal and axial coordinates, respectively, and B y (r, ϕ, y) is the vertical (axial) magnetic field component. The integration limit Γ 0 > 0 was chosen such that B y (r, |ϕ| > Γ 0 , y = 0) → 0. Both magnets were specified as six-degree sector magnets. It was found that the integrated field strength distribution generated by the 3D six-degree sector magnet models with straight edges deviates considerably from the ideal, "hard-edge", six-degree sector magnetic field distribution given by Eq. 1 (see Fig. 10). This means that pole edges of the magnets need to be adjusted accordingly. To achieve this the new pole edge shape ϕ = ϕ(r) where ϕ 0 = 6 • is the initial approximation to the magnet edge profile and I (r) is the specified value of the integrated field strength obtained from Eqs. 1 and 7. The ratio ∆I ∆ϕ was obtained by creating another magnet model with ϕ 1 = 5 • sector angle and subtracting the integrated field strengths obtained from both models. To simplify the task a discrete version of the goodfield region of each magnet was considered. The radial coordinates r 1 , r 2 , ..., r 5 split the good-field region into four intervals of equal length and the corresponding azimuthal coordinates ϕ(r 1 ), ϕ(r 2 ), ..., ϕ(r 5 ) were obtained from Eq. 8. The corrected pole edge shape was recovered by fitting a polynomial to the calculated nodes {r i , ϕ i }, 1 ≤ i ≤ 5. As Fig. 10 shows this procedure reduces the peak relative integrated field error by an order of magnitude. Further improvement can be obtained by "tweaking" the nodes individually until the desired accuracy goal has been reached. Figure 11 shows the pole edge profile of the F-magnet obtained from Eq. 8 compared to a 6 • -edge and a 5 • -edge.
Whilst the radial profile in the final 3D magnets deviates from that of the nominal design, the integrated field along each arc of constant radius is well matched by this adjustment of the edge geometry. Figure 12  The 2D and 3D magnet models can now be used to obtain field maps for particle tracking. In following sections we analyze the beam dynamics with these realistic magnets models.

IV. BEAM DYNAMICS WITH 2D MAGNETS
In this section we model NORMA using the field profiles from the 2D magnet models, to compare the dynamics with those due to the nominal fields.

A. Implementing the 2D field maps in Zgoubi
The 2D magnet design is output from OPERA as a table of the vertical component of the field, B y , along a radial line through the magnet. The B x and B z are zero by symmetry. These radial profiles can be used in Zgoubi with either the POLARMESH or DIPOLES elements.
To track the 2D model with POLARMESH a midplane field map must be generated. For each magnet the radial profile is combined with a longitudinal Enge fringe function (Eq. 2) to give the field at each point on a 2D mesh. We use the same Enge parameters as in the nominal model. These are summed to give the whole cell midplane field map. The process is shown in Fig. 13. The multipole fit for DIPOLES can be a global fit across the whole magnet giving a single model that can be used at every energy. Alternatively, a region around each orbit can be selected and fitted. The latter gives a higher accuracy, avoiding any smoothing from the fit, but means that the magnet description must be refitted at every energy. For simulations, we used 9 th -order local fits about each energy.

B. Dynamics and DA with 2D magnets
The 2D model gives good agreement with the nominal design. Figure 15 shows the horizontal and vertical tunes for the 2D models against the nominal design. The mean tune is within 10 −4 in both planes and the tune excursion grows from 6.4×10 −5 and 9.1×10 −4 to 7.9×10 −4 and 1.5×10 −3 between the nominal and the POLARMESH model. Agreement between the POLARMESH and DIPOLES simulations methods are very good, demonstrating that the multipole expansion is an effective method to use for tracking. An increased tune excursion may cause accelerated particles to cross resonances that reduce their stability. In order to confirm that this is not the case for the 2D magnets we determine the 1000-turn DA. Using DIPOLES with fits to the 2D radial field profile at each orbit, we find a 45 • DA (i.e. where particle amplitude is increased equally in both transverse planes) that is very close to the nominal lattice over the full range of energies; this is shown in Fig. 16. The figure shows that deviations from the ideal field in the 2D model do not cause a significant drop in DA. As with the nominal design, the DA is kept above 50 mm mrad over the full range of energies in the accelerator.
We can also add random multipole errors to the magnets. By using errors of a similar magnitude to those seen in the 2D field we can obtain some understanding of how the worst-case errors might affect the DA. However, applied multipole errors tend to cause big field changes away from the point that they are applied, which would not be seen in a well-designed real magnet. For a given error size, the ideal field is expanded around the mean position of a given orbit; Gaussian-distributed errors are then applied to a given multipole component. Figure 17 shows how the DA is reduced as quadrupole and sextupole errors applied around the 30 MeV orbit increase, using 20 seeds per error size; the dotted vertical lines Overall, the deviations of the 2D profile from the nominal design do not affect the dynamics enough to cause problems. We find that the increased tune shifts as a function of energy are tolerable. Random quadrupole and sextupole errors of similar magnitudes can also be added without causing the DA to drop significantly below our 50 mm mrad requirements.
Whilst the 2D simulations are a good indicator that it is feasible to produce the magnetic field profile speci-fied in the nominal lattice design, the following sections show that the longitudinal and full 3D fields are needed for a sufficient understanding of the lattice dynamics. 2D modelling is however a useful step in the design and optimization stage, as it can be performed more rapidly than with full 3D models.

V. BEAM DYNAMICS WITH FRINGE FIELDS
In this section the effect of the fringe extent and radial dependence is investigated. For the initial lattice design an Enge fringe field was used, with a 4 cm extent constant with respect to radius. In the real magnet, the fringe extent depends on the pole shapes. It should be expected that adjustments to the magnet strengths and/or field profile will be needed to account for the realistic fringe fields.
The fringe-field extent has a significant effect on the edge focusing of the magnets and therefore contributes to the tune of the cell. Changes in the fringe-field extent as a function of radius will therefore have an effect on the tune as a function of energy.
To show the effect on the tune, the fringe extent was varied while all other lattice and magnet parameters were held fixed and no rematching was performed. Figure 18 shows the change in tunes as a function of fringe extent; the original value of 4 cm in the nominal lattice is highlighted. As the tune shifts, the working point moves and approaches resonances. At an extent of around 10 cm the vertical tune approaches 0.25; the black points in Fig. 19 show how this causes a large drop in DA as the working point approaches a fourth-order resonance. The fringe-field extent can be varied by extending the fringe-field description from Eq. 2 with the variable κ and making λ a function of r, so that λ = λ 0 (r 0 /r) κ . A positive κ gives a fringe field with a larger extent at smaller radii, as would be expected due to the larger pole gap. Figure 20 shows the tune as a function of energy for a range of κ values. These tune shifts with energy can cause problems as the beam will cross resonances during acceleration, so it is important that they can be corrected. The tune shift due to a change in the fringe field can be compensated by changing the magnet body fields. For a given fringe extent the working point of the lattice can be rematched to recover the original tunes by adjusting the magnet strengths B 0 of the F and D magnets, and the field index k shared by both F and D magnets. The match is constrained to keep the outer closed-orbit position constant. DA is most critical at in-jection energy before the emittance is reduced somewhat by adiabatic damping, so the rematching was performed at 30 MeV. PyZgoubi's optimization feature, making use of the Nelder-Mead downhill simplex method, was used for the rematching. Figure 19 shows how the rematch recovers some of the lost DA for lattices where the change in fringe field caused the tune to approach a resonance.
We find that the original working point is not necessarily optimal for maximising DA, as the shape of the fringe field can affect the relative strength of the resonances. It is therefore important to re-optimize the working point to find the largest DA region. Re-optimization is performed by scanning the lattice through a range of working points and calculating the dynamic aperture at each point. Figure 21 shows how the working point is rematched for a fringe extent of 6 cm. Rematching the tune to the original values recovers some DA, but the scan is required to recover a DA above 50 mm mrad. The extent and shape of the fringe fields play a large role in the dynamics of an accelerator. The lattice must be re-optimized to account for their effects. This shows that 3D magnet simulation, which will give the fringe field, is a necessary step for realistic lattice design.

VI. BEAM DYNAMICS WITH 3D MAGNETS
In this section we model NORMA using field maps from 3D magnet simulations. These contain not just the radial field profile but also how the field falls off outside the magnet body.
A. Implementing the 3D field maps in Zgoubi 3D magnets are output from the OPERA simulation as midplane field maps, i.e. the B y component at grid points on a horizontal midplane; the B x and B z are zero by symmetry. The midplane is sufficient to fully define the field over the vacuum region of the magnet including the fringe field, while significantly reducing the computational resources required to generate and store the field map. Figure 22 shows the midplane fields for each magnet. The dynamics in the full magnet design will be strongly influenced by both the body field and the fringe fields. For example, the reduced peak field in the F magnet is compensated by an extended field length such that the integrated field is close to the nominal design. The 3D field maps must therefore be used as a whole, as the radial and fringe parts cannot be independently combined with the nominal field.
The midplane maps from the 3D magnet simulation can be used in Zgoubi with either the POLARMESH or DIPOLES elements. The F and D magnets are designed independently, as this allows simple boundary conditions and reduction in computation time due to symmetry. The midplane maps for the magnets are combined assuming a linear superposition to create the map for the full cell, which is then read into the POLARMESH element as shown in Fig. 23. At the field crossover point in the overlap between the F and the D magnets, the residual field from the magnets is about 5-10% of the body field, greatest at the low-energy orbits. Figure 24 shows the overlap at 50 MeV. A further step would be to combine the magnets within the FEM simulations so that the overlap region is fully calculated. If this is found to cause a significant effect, then either clamp plates can be used to reduce the fringe extent or rematching as described below can be used to account for it.
Alternatively, the 3D magnet can be modelled with a DIPOLES element, as shown in Fig. 25. This makes rematching simpler by reducing the description to a smaller set of variables, gives more reliable long-term tracking and improves performance compared to tracking in the map directly. The DIPOLES element has several parameters that allow it to model a range of magnets. The radial profile is specified as a 9th-order polynomial, the entrance and exit boundaries can be moved and rotated independently, the fringe can be specified with up to 6 Enge coefficients and the fringe extent can be a function of radius.
We fit the DIPOLES using only the magnetic field and then use tracking to verify that the solution is good. For the initial fit we minimised the difference between the map and the DIPOLES at each grid point with in the good field region, using a Nelder-Mead downhill simplex, allowing the above parameters to vary. This however did not give a good match for the dynamics. Much better dynamic agreement was found by including the differences of the integrated field and integrated gradient along the particle trajectories in the objective function.  The difference between the nominal fields and the 3D maps cause significant changes to the tune of the lattice. Figures 26 and 27 show the tune for POLARMESH in green and fitted DIPOLES in red, compared to the nominal FFAG in blue. It can be seen that the fit does a good job of reproducing the dynamics of the field map. At low energies the vertical tune with the 3D magnets crosses the quarter-integer resonance so we expect a large drop in DA below 100 MeV.  We simulate the 1000-turn DA in the 3D magnet using both the POLARMESH and fitted DIPOLES magnet elements. Figure 28 shows a significant reduction compared to the nominal FFAG design, especially for energies below 100 MeV where the vertical tune crosses below the quarter-integer resonance. It is clear that the uncorrected map does not result in a sufficiently large stable region.
In order to improve the DA the magnets must be re- This rematching is carried out by varying the multipole components of the fitted analytic field, and optimizing to obtain the nominal tune at all energies. The shape of the fringe falloff is held fixed as it is expected that a small change to the body field will not have a dramatic effect on the fringe, so that after another iteration of the magnet design the tune flatness will be retained. We found that a rematch varying multipole coefficients up to decapole in the bodies of the magnets was sufficient to flatten the tune such that it no longer crossed harmful resonances. The light blue lines (3D rematched) in Fig. 26 and 27 show the improvement in tune flatness when the magnets are rematched up to decapole order. Figure 28 shows the DA for the rematched magnets in light blue. The rematched tune is sufficient to increase the DA to around 50 mm mrad in the critical region below 100 MeV and for it to remain above 50 mm mrad at higher energies.
The main result from this work is the corrected (rematched) field profile for each magnet that recovers the original beam dynamics obtained with the idealistic (formula) fields. These field profiles together with the original (Map) field profiles obtained from the 3D Opera models are plotted in Figs. 29 and 30 for the F and Dmagnets, respectively. The plots labelled "Fitted" are practically identical to the "Map" fields and represent an intermediate result needed to obtained the corrected (rematched) fields. In the F magnet the difference between the fitted field and the rematched field is less than 5%. In the D at small radius there is a 20% reduction in the field, falling to 1% at large radius. These modifications to the body field are sufficient to rematch the magnets that would otherwise give significantly altered dynamics to the original design.
We now take the difference between the fitted and rematched field profiles, and add it to the original field  profile given by Eq. 1 in order to obtain the updated pole shapes of the two magnets that generate the rematched field profiles shown in Figs. 29 and 30. The first step in this process is to obtain the updated lines of constant scalar potential. These are shown in Fig. 31 for the D-magnet (cf. Fig. 5) and represent the zerothorder approximation to the pole shape of the rematched D-magnet. At this stage the procedure described in Section III can be implemented to yield the revised magnet designs that generate the rematched field profiles followed by an analysis of the beam dynamics produced by these updated magnet designs. The process is repeated until the desired level of performance is reached. Table II shows the parameters of the F and D magnets achieved by our design process.

VII. CONCLUSION
The NORMA design for a medical proton accelerator has been previously demonstrated with idealized magnet modelling. More detailed modelling using realistic magnet designs is important to show that the design is robust, and to understand the retuning needed to maintain the required dynamics and stability of the original design. In this paper we have shown that the NORMA lattice can be modelled with magnet designs from 2D and 3D FEM simulations, taking the design from an idealized lattice to a detailed study with realistic magnets.
The F and D magnets were designed in 2D and 3D using OPERA. The field maps from the these models were imported into the PyZgoubi tracking code so that dynamics of the proton bunches could be compared to the idealized lattice. For the 2D models -which give a radial profile -no significant effect on the dynamics was found. Random multipole errors, with similar-sized deviations as between the 2D profile and the analytic field, were found to have only a small effect on the stability of the lattice. Performing the magnet modelling and optimization work entirely in 3D is not optimal because of the complexity of the problem. As shown in the text, highly accurate 2D pole shapes for the two magnets are easy to obtain. These solutions allow a quick and realistic assessment of the available vertical and horizontal magnet apertures, the attainable good field region and, equivalently, the energy range of the machine to be made. In addition, the 2D pole solutions are a convenient starting point for the actual design work in 3D. The 3D models in turn provide a realistic representation of the actual magnetic fields and allow a detailed study and mitigation of field errors, fringe field and magnet cross-talk effects.
To see the importance of the fringe fields, simulations with altered fringe extents were carried out. It was found that changing the fringe length had a significant effect on the focusing of the lattice, and that for any given fringe extent the lattice must be retuned. By rematching the field strength and index in the magnets the working point could be re-optimized and sufficient DA recovered. 3D magnet models predict greater shifts from the nominal fields. However, we show that it still possible to make small adjustments to the body fields that recover a sufficient DA over the energy range from injection to extraction; further iterations of the 3D design may be used to refine it. More detailed modelling of the overlap region between the magnets is also needed, though it is expected that any changes to the dynamics can be accounted for by similar re-optimization.
This work demonstrates the importance of relying on rigorous 3D magnet simulations to model an FFAG accelerator, in order to obtain a realistic assessment of the overall machine performance. It also demonstrates methods of retuning an accelerator design to account for these effects, and to recover the original dynamic properties and beam stability.