Analytical expressions for fringe fields in multipole magnets

Fringe fields in multipole magnets can have a variety of effects on the linear and nonlinear dynamics of particles moving along an accelerator beamline. An accurate model of an accelerator must include realistic models of the magnet fringe fields. Fringe fields for dipoles are well understood and can be modelled at an early stage of accelerator design in such codes as MAD8, MADX or ELEGANT. However, usually it is not until the final stages of a design project that it is possible to model fringe fields for quadrupoles or higher order multipoles. Even then, existing techniques rely on the use of a numerical field map, which will usually not be available until the magnet design is well developed. Substitutes for the full field map exist but these are typically based on expansions about the origin and rely heavily on the assumption that the beam travels more or less on axis throughout the beam line. In some types of machine (for example, a non-scaling FFAG such as EMMA) this is not a good assumption. In this paper, a method for calculating fringe fields based on analytical expressions is presented, which allows fringe field effects to be included at the start of an accelerator design project. The magnetostatic Maxwell equations are solved analytically and a solution that fits all orders of multipoles derived. Quadrupole fringe fields are considered in detail as these are the ones that give the strongest effects. Two examples of quadrupole fringe fields are presented. The first example is a magnet in the LHC inner triplet, which consists of a set of four quadrupoles providing the final focus to the beam, just before the interaction point. Quadrupoles in EMMA provide the second example. In both examples, the analytical expressions derived in this paper for quadrupole fringe fields provide a good approximation to the field maps obtained from a numerical magnet modelling code.


Introduction
Fringe fields represent regions that lie at the edges of a magnet where there is a transition from the nominal field to zero field, or to the nominal field in an adjacent magnet. In multipole magnets, the nominal field has no longitudinal component. However, in the fringe field region where the fields vary with longitudinal position, Maxwell's equations require the presence of a non-zero longitudinal field component. For dipoles, the nominal field has only a single component. In dipole fringe fields, therefore, the field has only two components, and (assuming that the fields are independent of the horizontal transverse co-ordinate) analytical expressions for the field can be obtained by solving the two-dimensional Laplace equation. For quadrupoles and higher-order multipoles, however, the fields in the fringe region have three components, and analytical expressions for these fields must be obtained by solving the three-dimensional Laplace equation. Fringe fields can impact the motion of particles passing through magnets in a number of ways. For example, they can introduce nonlinearities in the equation of motion, or they can make a substantial contribution to the desired effects of the nominal field. The latter situation is the case in EMMA, for example, where the large aperture of the quadrupoles compared to their lengths means that the fringe fields dominate the focusing effects of these magnets. Nonlinearities from fringe fields can become important if the transverse size of the beam is large, or if the beam traverses a multipole at an angle and some distance from the magnetic axis: this is often the case in final focus quadrupoles in colliders. The implementation of fringe fields is also important in some tracking codes which include effects like space charge, as is the case in GPT (General Particle Tracer) [1], for example. This requires that all fields be continuous so that there are smooth regions where the fields transit from their maximal value to zero and vice versa.
There are several models available for the study of fringe fields for multipoles, see for example [2,3] and references therein. However, existing models are usually limited to on-axis and mid-plane approximations, meaning that the field in the full space of the multipole typically has to be computed with elliptic integrals. Therefore, for simple but accurate particle tracking, a significant amount of effort and computing power goes into just the calculation of the fields that particles see. The results obtained in this paper make possible an alternative method, based on analytical expressions for the fields as functions of position, that provide exact solutions to the static Maxwell equations in three dimensions. This allows for arbitrarily smooth fields to be constructed and used in tracking codes as well as the possible creation of transfer maps. The fringe fields considered here have associated scalar and vector potentials. From the scalar potential it is possible to inspect the shape of the pole-face for an iron-dominated magnet generating the given field. Using the vector potential it is possible to perform symplectic integration of the equations of motion for a particle in the field, leading to the construction of transfer maps for the fringe region. Further, it may be possible to improve the efficiency of the magnet design process by making some initial assumptions based on the formulae presented in this paper. However, a full three-dimensional numerical field map will ultimately be needed (obtained from a numerical magnet modelling code) in order to achieve the accuracy that is needed for validating the design of many accelerators. This paper develops the mathematical framework that was initially presented in [4], gives further details and results, and presents two examples. Following the Introduction, fringe fields for dipoles are briefly reviewed. The formalism used for dipoles is extended to fully three-dimensional fields in Section 3. A complete solution to the static Maxwell equations, in a form suitable for application to multipole fringe fields, is derived and presented. Expressions for fringe fields in multipoles of arbitrary order are then given in Section 4. A particular case of a quadrupole fringe field, together with a of the simple fall-off in the form of an Enge function [5], is then presented in Section 5. All the salient properties are described in order to demonstrate that the field behaves in the way expected of a quadrupole, first inside the magnet, then in the fringe field region, and finally at a large distance from the magnet (so that the field effectively falls to zero). In Section 6, the scalar and vector potentials for fringe fields in the case of a multipole of arbitrary order are discussed. In Section 7, two quadrupole examples are presented. The first is the HL-LHC inner triplet, where the beam size and trajectory in the quadrupoles providing strong focusing close to the interaction point make fringe field effects significant. The second example is based on the quadrupole magnets in the non-scaling FFAG, EMMA. The large aperture of these magnets compared to their length means that the fringe fields make a dominant contribution to the focusing effects. Conclusions are given in Section 8.

Fringe fields for dipoles
The goal is to derive expressions for multipole fringe fields that satisfy Maxwell's equations. To ensure the validity of the solution and the corresponding assumptions, it is important to write all equations explicitly. For static fields in the absence of any electric current, the equations for the magnetic field B are: For dipole magnets, it is sufficient to consider a two dimensional version of the equations. Taking B x = 0, we are left with: together with: which excludes all dependence on x. Further, we seek fringe fields which have a possible fall-off along the axis of the magnet given by the six parameter Enge function [5,6]: 1 1 + e E(z) , with E(z) given by: All coefficients a n are constants determined by modelling and/or experiment. D is the full aperture of the dipole. The main advantage of the Enge function is that it is analytic and can be made to tend to asymptotic values arbitrarily fast, if required. The main disadvantage of this function is that, if all coefficients a n with n = 0, . . . , 5 are included, varying any one coefficient changes the effective length of the magnet (in this case, a dipole). Other functions which decay sufficiently rapidly may be used instead of the Enge function [7,8]. For simplicity, we consider in this paper only the case where a 1 = 0 and all other coefficients are set to zero. This has the further advantage that varying a 1 only changes the fringe field "hardness" without altering its length. Additionally, because we only have one non-zero coefficient, we can normalise it to 1 without loss of generality. Maxwell's equations (1) imply: where ∆ y,z = ∂ 2 y + ∂ 2 z . Both wave equations (for B y and B z ) can easily be solved: Requiring that equations (1) be solved as well, we end up with: If we further restrict ourselves to real magnetic fields, we obtain: B y and B z are given by twice the real and imaginary parts of the function e(z + iy), respectively. Traditionally, dipole fringe fields are described by: where B 0 is the nominal strength of the dipole field, and g is a parameter (related to the aperture of the magnet) affecting the precise shape of the fringe field. Using the result: B z may be re-written as: whence e(z + iy) has the form: Then, using equation (3) we find B y to be given by: which may be converted into the same form as the B y component of the magnetic field used in GPT (General Particle Tracer) [9], for example. So, in summary, given a complex function e(z + iy), Maxwell's equations may be satisfied automatically by setting B y = 2Re(e(z + iy)) and B z = 2Im(e(z + iy)).
A possibility for having a magnetic field whose B y component falls off on-axis is given by the six parameter Enge function: This would force B z to have the form: for some complex function E(z + iy). If we consider the simple case E(z + iy) = z + iy (with aperture having unit diameter) then equations (5) and (6) simplify to: B y = (1 + e z cos y) 1 + 2e z cos y + e 2z , B z = −e z sin y 1 + 2e z cos y + e 2z .
This may be extended to include as many parameters of the Enge function as desired, the only restriction being that E = E(z + iy). So, if E(z + iy) = b 1 (z + iy) + b 2 (z + iy) 2 , with b 1 and b 2 arbitrary constants, we have: where f = b 1 z + b 2 (z 2 − y 2 ) and h = y(b 1 + 2b 2 z). This can be generalized to: with b j constants if an N parameter Enge function is desired. Indeed, any function can be used, provided it is a function of z + iy.
We can plot the fringe fields (5) and (6) in the simplest case (E(z + iy) = z + iy), as shown in Fig. 1. Each field component has singularities, of which two are visible in the plots. From the denominator of (5) and (6), it can be seen that the singularities appear when (re-introducing the arbitrary aperture diameter D) 1 + 2e z/D cos(y/D) + e 2z/D = 0 or cos(y/D) = −1 → y/D = ±nπ where n is an integer. This does not cause any problems for modelling the field as it is always possible to arrange that the singularities are located outside the region of interest. It is not possible to avoid the appearance of singularities since any solution to Laplace's equation either has singularities or is constant. In practical terms, the singularities can be associated with the current in the windings of the dipole.
Note that the singularities have a tan-like behaviour and that, at the exact point where z = 0, the value of the field is precisely half of the maximum and everything is smooth. This can be verified by setting z = 0 in the expressions for the dipole fringe fields above and using l'Hopital's rule.
Locally, it is possible to find potentials for the dipole fringe fields. The usual vector potential A is related to the field by B = ∇ × A. For magnetostatic fields in the absence of currents, the field can be derived from a scalar potential ϕ, by B = ∇ϕ. The existence of the vector potential is ensured by the Maxwell equation ∇ · B = 0, while the scalar potential exists because ∇ × B = 0. For the simple case discussed above, the scalar potential is given by: with the only gauge freedom being given by the constant. From the scalar potential, it is possible to obtain a description of the pole-face geometry, since this is given by surfaces where ϕ is constant. This is shown in Fig. 2 where two profiles of the pole-face can be seen. It should be remembered that there is a scale invariance in the expressions depending on the dimensions of the dipole. Figure 2 also shows that the pole-face profiles encompass the singularities: this is consistent with the assertion made earlier that the singularities are associated with the current in the coils of the dipole. The vector potential has extensive gauge freedom and, in one of its simplest forms, can be given as: with the other components of A set to zero.

General three dimensional solution in a fringe field
General expressions for the fringe fields in a dipole followed from writing Maxwell's equations in the form (1). In a dipole, we only needed to consider two field components (B y and B z ) as functions of two co-ordinates (y and z). For higher order multipoles, it is necessary to consider the dependence of all three field components on all three co-ordinates. One approach might be to look for solutions to the three dimensional Laplacian. A formal solution (due to Whittaker [10,11,12]) is known, and can be expressed as: where: ∆ x,y,z ϕ ≡ ∂ 2 x ϕ + ∂ 2 y ϕ + ∂ 2 z ϕ = 0. However, it is difficult to use this result to find expressions that are useful in practice. The only well-known practical solution is f = (z + ix cos ϑ + iy sin ϑ) −1 , which gives the standard solution ϕ = 2π/r with r = x 2 + y 2 + z 2 .
An alternative approach, which we develop here, is to define new variables in terms of which Maxwell's equations can be written for a general three-dimensional magnetostatic field in a form very similar to (1). This raises the possibility of finding expressions for fringe fields in higher order multipoles by applying similar techniques to those described in Section 2 for dipoles. To write Maxwell's equations for three-dimensional fields in a form similar to (1), we define new variables: We express the magnetic field in terms of components: In terms of the new variables, Maxwell's equations can be written: From (7) and (8), one can see immediately that, in the absence of any fringe fields, the general solution of Maxwell's equations for any magnet, acting transversely only and without fringe (B ζ = 0) is given by: for any functions f (v) and h(u). The case of a multipole of order n (n = 1 for a quadrupole, n = 2 for a sextupole, and so on) is given by: For example, a quadrupole is described by B u = iv, B v = −iu and B ζ = 0.
We now assume a very general form that we expect multipole fringe fields should take: The form of this field is based on a generalisation of the fringe fields for the dipole case. Substituting (11)-(13) into Maxwell's equations (7)-(10) gives a set of constraints on the possible forms of the functions f 1 , f 2 etc. To obtain useful expressions for multipole fringe fields, we need to find solutions satisfying the various constraints: this is the task that we address in the remainder of this section.
Essentially, there are only two types of derivative that we need to consider. These are: where A = 1 + 2f 3 e ζ + e 2ζ . For the remaining derivatives, we simply implement the following changes sequentially: As all equations (7)-(10) are equal to zero, we can take out a factor of A 2 to give: (∂ x f 6 + ∂ x f 7 e ζ )(1 + 2f 3 e ζ + e 2ζ ) − (∂ ζ f 1 + ∂ ζ f 2 e ζ + f 2 e ζ )(1 + 2f 3 e ζ + e 2ζ ) −2(f 6 + f 7 e ζ )e ζ ∂ x f 3 + 2(f 1 + f 2 e ζ )(f 3 e ζ + e 2ζ ) = 0, We can now equate coefficients of e ζ giving: Note that we have not included all the steps and the above equations represent the original set with all possible simplifications taking into account the set itself. Note also that equations (14)- (17) and (26)-(29) may be solved independently of the rest and can therefore be dealt with later. From equation (21), using (18) and (22), we see that: Had we looked at equations (20) and (24) instead, using (19) and (23), we would have had: with the same result from equation (25). Now, f 6 and f 7 cannot both be zero as this would mean B ζ = 0. Therefore, we must have: The general solution for f 3 , using the method of characteristics [10], is: where b and c are constant. Without loss of generality and for simplicity, we set c = 0. Substituting this back into (18)-(25) gives the relations: The equations reduce to just two independent equations that may be written as: Using f 1 = b 2 f 4 and equations (26) and (27) we see that we require which can again be solved by the method of characteristics to give f 4 = f 4 (h, z). Using this with equations (28) and (29) we see that (14) and (15) and, subsequently (16) and (17) gives and f 7 = f 7 (h, z). This leaves six equations to be satisfied from the original system (14)-(29), namely (31) and (32) together with: After cross-differentiation, equations (35) and (36) give Now, we can re-express the partial derivatives in u and v in terms of h only. The equations simplify to: We introduce the co-ordinates: Note that this operation is equivalent to complex conjugation in the ζ co-ordinate only and the function h is untouched. Therefore, we have the solutions: Substituting this back into (35) and (36), we see that the solutions are further constrained to: with k constant, from which we can get f 4 via f 4 = 1 b 2 f 1 . Subsequently, we can get f 2 from (32) and hence f 5 via f 5 = 1 b 2 f 2 and f 7 from (31). The general result, in terms of p 6 and q 6 may be summarised as follows (with h = 1 b u + bv): We are left with equations (33) and (34) to be solved. Upon substitution of (37)-(43), these are actually seen to be trivially satisfied with no further constraints on any of the functions f i . In fact, the results constitute a Darboux transformation, where, given a solution to Maxwell's equations expressed by (26)-(29), a new solution may be created, given by (11)-(13), provided (14)-(25) are satisfied. Further, the results can be seen to imply the following solution to Maxwell's equations: with h defined by (30). This solution could have been arrived at by a shorter method, but the above workings show that, for the kind of functions we are interested in, it is the only solution that fits. We may have been forced to consider more complicated functions by introducing non-linearities in the coordinate dependance, for example. Note that the constant b is purely a scaling constant of the coordinates, as well as giving a proportionality between B u and B v . For the purposes of the remainder of this article, we shall refer to the solution (44), (45) and (46) as an elementary solution of Maxwell's equations in three dimensions. In the original (Cartesian) co-ordinates ( 2) the elementary solution reads: where: and h is now expressed as: For some applications, it is useful to have expressions for the scalar potential ϕ and vector potential A, from which the field can be derived by: or: It can be readily verified by substitution into (50) that an expression for the scalar potential for the elementary solution can be written: Setting the lower limit of each integral to zero fixes the gauge so that ϕ = 0 at the origin; a change of gauge can be made simply by adding a constant to ϕ. It can be verified by substitution into (51) that the vector potential for the elementary solution may be written: Here, we have chosen a gauge in which A x = 0. Note that any function of h ± iζ may also be expressed as a function of ζ ∓ ih: we shall use both forms in what follows in order to simplify the expressions as appropriate.

Multipole magnets with general fringe fields
The relations f 2 = b 2 f 5 and f 1 = b 2 f 4 found earlier imply that B u ∝ B v : this means that no physical magnetic fields can be represented this way. However, because of the linearity of Maxwell's equations, any linear combination of elementary solutions gives a solution to Maxwell's equations. Physical solutions corresponding to multipole fields can be constructed by taking appropriate combinations of elementary solutions. Further (physical) constraints are that the field decays as z → ∞ and that the field matches the nominal multipole field inside the magnet. In this section, we describe in detail the procedure for constructing an expression for the fringe field in a quadrupole, and then generalise our results to multipoles of any order. To obtain the correct behaviour of the field as a function of z, the functions f (h + iζ) and g(h − iζ) are each written as a product of two factors. The first factor, (h + iζ) n and (h − iζ) n respectively, represents the (nominal) multipole for some constant n. The second factor, F (h + iζ) and G(h − iζ) respectively, are multiplicative functions which are chosen to give the desired decay as z becomes large (z → ∞), and to give the correct field in the body of the magnet. Thus, we write: Following the form of (44), (45) and (46), let B u , B v and B ζ be given by linear combinations of the elementary solution: where h j = 1 bj u + b j v, and m is a constant determined by the number of copies of the elementary solution needed to construct a physical multipole of the desired order. The constant n determines the order of the multipole. Within the body of the magnet, B u = iv n , B v = −iu n , and B ζ = 0. We also expect to have to satisfy several relations between the constants b j and c j .

Quadrupole magnets
A quadrupole field is obtained by putting n = 1 in equations (52) and (53). Therefore, within the body of the magnet (far from any fringe field) F (h j + iζ) = G(h j − iζ) = ±1 for any j. The relative signs of F (h j + iζ) and G(h j − iζ) are to be determined by the required gradient inside the quadrupole and the fact that there should be no dependence on z at the centre of the magnet. The required behaviour of the fields in the body of the magnet can be obtained by taking m = 2 (larger values of m can be used, but are not required). Then, putting B u = iv, B v = −iu and B ζ = 0 gives: Thus, we have three equations to satisfy: Therefore, inside the magnet, the coefficients b j and c j must satisfy the following constraints: We are left with the freedom of choosing one constant, which we take to be b 2 . We shall consider the significance of this constant in more detail later, but for now we note that to avoid b 1 , c 1 or c 2 becoming singular, b 2 must not be equal to 0 or 1. Also, the field is unchanged if we replace b 2 by 1/b 2 . Therefore, before implementing any simplifications, we have the following: Physically, it is necessary that our choices of constants should give a real field in the end. Now, since the two co-ordinates are related via u =v (where the bar denotes complex conjugation) h 2 =h 1 to within a factor ±1. In detail, the first constraint above implies the following two cases: Now, in either case, the second constraint (c 1 = −c 2 ) is the same and we factor the resulting constant out for clarity so, in terms of magnetic field components we are left with (for b 1 = +1/b 2 ): where we have put h 1 = h. On the other hand, if we take b 1 = −1/b 2 , then the field components read: So we can see that, for both choices, B z is real, however, this is not the case for the other two components. In the first case (b 1 = +1/b 2 ), it is not possible to say anything about B x or B y unless f = g, in which case both are real. However, in the second case (b 1 = −1/b 2 ), all components are automatically real, irrespective of the choices of f and g. Note that allowing the constant b 1 to take complex values only increases the complexity of the equations without bringing in any useful additional parameters; therefore, without loss of generality, we consider only cases where b 1 takes real values. Therefore, and so as to be able to keep f and g independent of each other, we choose The results for fields in quadrupole magnets can be extended to higher order multipoles in a straightforward way. However, the higher the order, the more copies of the elementary solution are needed to construct fields with the required properties. In the following subsections, we shall consider explicitly the cases of sextupole and octupole magnets, to establish a pattern from which the results for a general multipole magnet can be written down.

Sextupole magnets
A sextupole field is obtained by putting n = 2 in equations (52) and (53). To obtain the constraints on the various coefficients c j and b j , we can proceed in close analogy with the case of the quadrupole magnet. In particular, within the body of the magnet (far from any fringe field) F (h j + iζ) = G(h j − iζ) = ±1 for any j, and the field is given by: The corresponding constraints are: There are four equations to satisfy, namely: Finally, the constraints on the coefficients b j and c j may be written: Note that, for sextupoles and higher order multipoles, it is unlikely that a particular choice of the associated constants results in the field components being automatically real. This is because of the additional number of constants involved. However, this should not result in any problem because the real part can just be taken at the end of the explicit construction given that any complex solution to Maxwell's equations implies that, both the real and imaginary parts, are also solutions of the same equations.

Octupole magnets
An octupole field is obtained by putting n = 3 in equations (52) and (53). The field is given by: Following the same procedure as for the quadrupole and the sextupole, the constraints on the coefficients c j and b j are: There are five equations to satisfy, namely: The constraints on the coefficients b j and c j can be written: , ,

General multipole magnets
The derivation of the constraints on the coefficients c j and b j in the case of quadrupoles, sextupoles and octupoles can be extended to any order of multipole. In the body of the multipole, the components of the field are: Substituting from (54), (55) and (56) gives: There are n + 2 equations to satisfy. Two of the equations can be written: The remaining equations have the form: where, for convenience, we have defined α = n − 2p with p a positive integer. The number of equations is determined by the order of the multipole and the restriction that α ≥ 0. Solving the equations leads to the following constraints on the coefficients b j and c j : Ultimately, we shall see that the constant b, for quadrupoles, can be used as a parameter to fit the field behaviour with distance from the axis, in the same way that the Enge coefficients can be used as parameters to fit the field behaviour with distance along the axis. The constant b and the Enge coefficients are ultimately dependent on the geometry of the magnet.

Significance of parameters and scaling laws
The number of constants increases as (n − 2)/2 with increasing order n of the multipole. For a quadrupole, there are effectively two transverse directions in which the rate of the fall-off of the field in the fringe region can vary. The rate of fall-off in one direction can be adjusted by an overall re-scaling; to vary the rate of fall-off in the second direction (independently of the first direction) requires a separate parameter. Similarly, in a sextupole, there are three independent transverse directions in which the rate of fall-off of the field can vary, which means that there are two independent parameters to control the transverse behaviour. For octupoles, with four independent transverse directions, three parameters are required, and so on.

Multipole magnets with Enge-type fringe fields
Having constructed general expressions for fringe fields in multipole magnets, it is worth investigating these expressions further, to show that the field has the appropriate behaviour. For a multipole magnet of arbitrary order, the expressions for the field can be rather complicated: therefore, we consider in detail only the case of a quadrupole. In this section, we shall discuss the behaviour of the fields in a quadrupole, first for the elementary solution, and then for the full solution (in which the fields are given by real numbers) obtained by adding two versions of the elementary solution. In order to plot the field, we need to make some assumption for the form of the roll-off of the field along the axis of the magnet: we shall consider the case that the roll-off is described by an Enge function. At the end of the section, we shall generalise the expressions for a full quadrupole solution with Enge roll-off of the gradient to higher order multipole magnets.
Of particular interest is the appearance of singularities in the magnetic field. Singularities are expected from the properties of Laplace's equation: in two dimensions, solutions to Laplace's equation are either constant everywhere, or have singularities somewhere. Our expressions for multipole fringe fields have been obtained by extending the two dimensional case to three dimensions. Mathematically, it is no surprise that singularities appear, but if the expressions we have derived for the fields are to be applied in physical situations, we should understand the position and nature of the singularities.

Elementary solution for a quadrupole with Enge-type fringe field
We take n = 1 for the elementary solution and examine its behaviour. For simplicity we assume that the fall-off functions are related by F 1 = −G 1 . To simplify the notation we drop the subscript j = 1, so the field components (54)-(56) become: We shall see when we consider the full quadrupole solution in Section 5.2 that an Enge-type fall-off in the quadrupole gradient occurs when we choose the function G(ζ) as follows: The components of the field (in a Cartesian basis) can be written in this case: Along the axis of the magnet (h = 0), the transverse components of the field vanish and the longitudinal component of the field is purely imaginary. If we identify the real parts of the field with the physical magnetic field, then the longitudinal component of the field also vanishes along the axis. Singularities in the field occur when: where is any odd integer. In terms of the Cartesian co-ordinates, the singularities occur at: Note, however, that the singularities in the different terms in the transverse components of the field cancel out when z = 0. The behaviour of the field can be seen in Fig. 4, which shows the real parts of the B x and B z components of the magnetic field as functions of the transverse co-ordinates at various z locations. For the elementary quadrupole-like solution, the constant b amounts to a re-scaling in the coordinates and the field; therefore, we show plots for only a single value b = 0.1. Also, since no new features occur in the behaviour of B y compared to the behaviour of B x , we only show plots for B x and B z . The plots in Fig. 4 show that the field has the expected behaviour for a quadrupole: in particular, within the body of the magnet, B x is linear in the co-ordinate y, and independent of x. At increasing values of z, the slope of B x versus y decreases (the gradient falls off); at z = 0, the quadrupole gradient is half the value at large negative z. As z increases further, the gradient (and the field itself) falls to zero. The singularities in the field have the behaviour expected from (59). In particular, we see that at z = 0, the singularities in the transverse field disappear completely.
Note that it was not necessary to demand that F = −G in the fall-off functions above. This was only done for simplicity and it could well be the case that a non-symmetric fall-off is desired and this would require F and G to be completely different. As they are just functions that dictate the behaviour the fringe fall-off, nothing changes inside the quadrupole and the magneto-static Maxwell equations are still satisfied.

Full solution for a quadrupole with Enge-type fringe field
As shown above, we can construct physical (real) fields in a quadrupole by adding two versions of the elementary solution with n = 1: where: h j = d j x + ie j y, and: together with the definitions of b j and c j given in (57) and (58), keeping in mind we took the negative solution in equation (57). Working in cylindrical polar co-ordinates r, θ, z, with the x axis corresponding to the line z = θ = 0, the radial field component can be expressed as a Taylor series in r: where G is the derivative of G with respect to its argument. Equation (63) describes the behaviour we would expect from a magnet in which the quadrupole gradient g(z) varies along the axis as: A conventional model for the fringe field in a quadrupole magnet describes the gradient as an Enge function [13]. Again to keep the analysis as simple as possible, we consider the case that the gradient varies as an Enge function with a single parameter: Writing the gradient in this way determines z = 0 as the location along the axis at which the quadrupole gradient falls to half its nominal value within the body of the magnet. Integrating (64) gives: where the constant of integration has been chosen so that G(ζ) remains finite as ζ → 0. The components of the field can now be written (in a Cartesian basis): The quadrupole gradient g(ζ) and the corresponding function G(ζ) are plotted as a function of distance along the axis of the quadrupole magnet in Fig. 5. If we plot the two constants d j and e j as functions of b j (as shown in Fig. 3) we see that the field components B x and B y will decay at different rates. The transverse field components only show the same behaviour when b 2 = 0; this value of b 2 is not allowed, since as b 2 → 0, e 2 → ±∞. Therefore, the symmetry B x ↔ B y under x ↔ y has to be imposed. This can be done by adding equivalent expressions to B x , B y and B z , with x and y interchanged but keeping b 2 the same, and dividing the result by 2. Therefore, in what follows, the full quadrupole solution refers to the fully symmetric quadrupole solution. There may be cases where exact quadrupole symmetry is not required; however, we do not consider such cases here. Figure 6 shows the behaviour of the field components B x , B y and B z for the full quadrupole solution, with the fixed value of the parameter b 2 = 0.1, and with the imposition of the symmetry B x ↔ B y under x ↔ y. Within the body of the quadrupole (large negative z) we see that, as expected, B x ∝ y, B y ∝ x and B z ≈ 0. At z = 0 the quadrupole gradient (the constant of proportionality between B x and y, or between B y and x) falls to half of its value within the body of the magnet. At a large distance from the quadrupole (large positive z), the field approaches zero.

Full solution for a multipole magnet with Enge-type fringe field
The full solution for a multipole of order n can be written: The function G 0 (ζ) is an Enge function (with one coefficient): Functions G n (ζ) for n > 0 are obtained by induction: It is then found that expanding the radial field component B r as a Taylor series in r gives: where O(r n+1 ) represents terms in r of order n + 1. Thus, the field has the behaviour expected of a multipole magnet, with gradient that decays as an Enge function in the fringe field region. The functions G n (ζ) can also be written in the form: where Li n (ζ) is the polylogarithm (or Jonquière function [14]) of order n.

Potentials for multipole magnet fringe fields
For some applications it may be useful to have analytical expressions for scalar and vector potentials from which the fringe fields in multipole magnets can be derived. For example, in an iron dominated magnet the pole faces correspond to surfaces of constant scalar potential; this makes it possible to inspect the geometry of a magnet that would have a fringe field of a given form. The vector potential can be useful for symplectic tracking of particles through the field (see, for example, [15]). General expressions for the scalar and vector potentials for the elementary solution were given in Section 3. In this section, we consider in particular the case that the multipole gradient g n (ζ) has a roll-off in the fringe field that can be represented by an Enge function.

Scalar potential and pole-face geometry
Consider a fringe field of the form: In the case that: where the functions G n (ζ) are defined in (65) and (66), this field represents the full solution for the fringe field in a multipole magnet of order n. For the moment, we do not assume any particular form for the function g n (ζ). The field (67)-(69) may be derived from a scalar potential given by: where:g This expression for the scalar potential may be verified by substituting into B = ∇ϕ. If it is possible to integrate the function g n (ζ) describing the fringe field, then it is possible to write down a scalar potential for the field. In the case of the fringe field of a multipole magnet of order n, with roll-off given by an Enge function,g n (ζ) is given by: Hence, in a multipole magnet of order n with fringe field falling off as an Enge function, the scalar potential can be written: In an iron dominated magnet, the pole faces can be identified with surfaces of constant scalar potential: an example of such a surface (for an inner triplet quadrupole for the HL-LHC) is shown in Fig. 11. The surface has the geometry expected of the poles in an iron dominated quadrupole. Within the body of the magnet, the intersection of the pole face with a plane z =constant forms a hyperbola. As a function of position along the axis, the pole face has little dependence on z for large negative z; but in the fringe field region the pole face abruptly "flattens" off to lie close to the plane z = 0.

Vector potential
It can be useful to know the vector potential for a given field for symplectic integration of the equations of motion for a charged particle moving through the field: see, for example, [15]. It is possible to write down an expression for a vector potential A from which the field (67)-(69) may be derived. In a gauge with A x = 0, the components A y and A z are given by: This form for the vector potential may be verified by substitution 1 into B = ∇ × A.
In the case that the field (67)-(69) represents the fringe field in a multipole magnet: with G 0 (ζ) an Enge function, it is possible to express the vector potential in terms of the functions G n (ζ) (defined in (65) and (66)) using:g

Examples
To illustrate the application of the methods and results described in the previous sections, we consider two examples: a quadrupole in a final focus or "inner triplet" region of the high luminosity upgrade of the Large Hadron Collider (HL-LHC), and a quadrupole in the EMMA non-scaling fixed-field alternating gradient accelerator. These magnets are chosen to provide two contrasting cases in terms of magnet technology and parameter regime. The HL-LHC inner triplet quadrupole is a superconducting magnet with large aperture and gradient (respectively, 150 mm diameter and 140 T/m). Design studies are still in progress; some field maps are available, but the impact of the fringe fields on the beam dynamics are still under investigation. The EMMA quadrupoles are normal-conducting electromagnets with more conventional gradient. Here, we consider one of the two types of EMMA magnet, namely the F quadrupoles [16]. These magnets have large apertures (74 mm diameter) given the length of the iron poles (73 mm); as a result, the fringe fields make a dominant contribution to the focusing effects. The specified integrated gradient is 0.387 T. Although it is possible to represent an EMMA quadrupole by 1 The result d 2 j = e 2 j + 2 may be useful for verifying the expressions for the vector potential. a hard-edge model, beam dynamics studies (supported by experimental results) [17] indicate the need for a more realistic representation in order to give an accurate description of the longitudinal and transverse dynamics in the machine. We should emphasise that the purpose of considering these illustrative cases is not to demonstrate close agreement between the numerical and analytical field models, but simply to show that it is possible, using analytical formulae from the previous sections, to construct a model of the fringe field in each case that satisfies Maxwell's equations, and is closer to reality than a simple hard-edge model.
For the HL-LHC inner triplet quadrupole and for the EMMA F quadrupole we use the full solution for a quadrupole with Enge-type fringe field, presented in Section 5.2. The numerical field data in each case are fitted using a roll-off function such that the quadrupole gradient has the form of an Enge function: .
The fit is optimised by varying the parameters a 0 , a 1 and a 2 . Following the same procedure as in Section 5.2, we find that the radial component of the magnetic field can be written to first order in the radial co-ordinate r: If the radial component of the field is known as a function of longitudinal position for some θ and (small) r, then a 0 , a 1 and a 2 can be obtained from a nonlinear fit to the data. With the gradient described by an Enge function as in (71), the function G(ζ) in (60)-(62) is given by: All field components at any point in the magnet (including the fringe field) can then be obtained using (60)-(62).

HL-LHC inner triplet quadrupoles
As a first example, we show the results of a fit to the field in an HL-LHC inner triplet quadrupole. Field data are obtained from a magnetic model, and a fit to the field based on a function of the form (73) is then obtained using the method described above. The parameter values obtained by fitting the radial field component along a line at θ = π/4 and various values of r are shown in Table 1. If the chosen function provides a good description of the field, then we expect little variation in the parameters obtained by fitting along lines with different θ and r, at least for small r (since the fit is based on the linear term in a series expansion of the field). For the HL-LHC inner triplet quadrupole, inspection of the values shown in Table 1 shows changes in the fit parameters of less than 0.7%, comparing values obtained by fitting the field along lines with r = r max /10 and r = r max /6 (and fixed θ = π/4).
Having obtained values for the fit parameters, we can make direct comparisons of the field given by equations (60)-(62) with the field obtained by the numerical magnetic model. An example of such a comparison, for the radial field component, is shown in Fig. 7. Although the analytical model does not match the numerical model exactly, the general behaviour of the field is reproduced quite closely. The impact of the residuals from the fit on the beam dynamics still needs to be studied; but the analytical model does include features of the field that would be completely omitted in a hard-edged magnet model. Figure 8 shows the radial field component as a function of position along lines of fixed radius and for different values of the polar co-ordinate θ. Again, there is reasonable agreement between the numerical field data (black lines) and the analytical fit (red line). Finally, Fig. 9 shows the longitudinal field component as a function of position along lines of given radius and with fixed θ = π/4. Here, it appears that there are more significant discrepancies between the numerical field data and the analytical model; however, there is general agreement in the main features, especially for the region close to the axis of the magnet. The more detailed structure that appears at larger distances from the axis cannot be reproduced by the relatively simple (Enge) function that is used to describe the fall-off of the quadrupole gradient in the fringe field.
It is worth considering the dependence of the field on the parameter b 1 (= 1/b 2 ). Changing the value of b 1 has an effect on the way that the fringe field varies with distance from the magnetic axis. This can be seen by comparing the plots in Fig. 10, which show the radial field component as a function of position along lines with fixed polar co-ordinate θ = π/4 and different distance r from the axis of the magnet. The plot on the left in Fig. 10 shows a comparison between the numerical field map and the analytical model with b 1 = 1.5; the plot on the right compares the numerical field map and the analytical  Table 1 with radius of fit r max /10, and b 1 = 2.5.  Table 1 with radius of fit r max /10, and b 1 = 2.5.  Table 1 with radius of fit r max /10, and b 1 = 2.5.   Table 1. model with b 1 = 3.5. In regions close to the axis of the magnet, changes in the value of b 1 in the range 1.5 to 3.5 have little effect on the field; but differences are apparent at larger distances from the axis. The parameter b 1 therefore allows control over the field behaviour at large r, after the parameters a 0 , a 1 and a 2 have been chosen to fit the field behaviour close to the axis of the magnet.
Although the HL-LHC inner triplet quadrupoles are superconducting magnets, we can nevertheless inspect the shape of the pole face that would be needed to give the same field. Assuming poles with infinite magnetic permeability, the shape of a pole face corresponds to a surface of constant magnetic scalar potential ϕ; in the case of a quadrupole with gradient falling off as an Enge function, the scalar potential is given by (70). Figure 11 shows an equipotential surface for the analytical field model fitted to the numerical HL-LHC inner triplet quadrupole field data. The surface has been chosen so that the magnetic scalar potential has the (arbitrary) value ϕ = 0.25 T m. We see that the equipotential surface has the shape that might be expected of an iron-dominated quadrupole, with the curved surace following an hyperbola in a plane of constant negative z, and the end of the pole being reasonably flat (close to a plane of constant z ≈ 0).

EMMA quadrupoles
The quadrupole magnets in EMMA are iron-dominated, normal conducting magnets. The unusual feature of the EMMA quadrupoles is that the diameter of the aperture is comparable to the length of the magnet; this allows the accelerator to have a large transverse acceptance in a lattice consisting of magnets packed very close together. However, the relatively large aperture of the quadrupoles means that the gradient falls off rapidly from the centre of the magnet. There is no appreciable  Figure 12: Radial component of the magnetic field in an EMMA F quadrupole as a function of position along lines of given radius and with cylindrical polar co-ordinate θ = π/4. Each pair of black and red lines shows the field at a different radius, from r max /5 to 4r max /5, with r max = 36 mm. The black lines show the field obtained from the (numerical) magnetic model; the red lines show the analytical model (60)-(62), with fit parameters given in Table 2 with radius of fit r max /10, and b 1 = 1.8. The centre of the quadrupole is at the far left of the plot, z = −53.6 mm.
distance along the axis for which the gradient is constant, and a realistic model of the field must include some representation of the fringe fields.
Given numerical field data from a magnetic model of an EMMA quadrupole, we can repeat the analysis used for the HL-LHC inner triplet quadrupole in Section 7.1. We again use an Enge function of the form (71) to represent the fall-off of the gradient along the axis of the magnet. The field data cover an area closely approaching the pole tip (the field values are given on a rectangular grid, with transverse co-ordinates extending to 36 mm). The results of fitting the parameters in the Enge function to field data along a line parallel to the axis (at different distances from the axis and fixed polar angle θ = π/4) are shown in Table 2. We see that there is much larger variation in the parameters if the fit is performed at different distances from the axis, compared to the case of the HL-LHC inner triplet quadrupole. This suggests that the quality of the fit will not be as good. Making a direct comparison of the field based on the analytical formula with the numerical field data confirms that this is the case: see Figs. 12, 13 and 14. Inspecting Fig. 12 suggests the reason for the poor quality of the fit. The gradient initially falls off quite rapidly along the axis from the centre of the magnet, but there is a long "tail" as the gradient approaches zero: this asymmetric behaviour cannot be represented accurately using an Enge function with a small number of coefficients. Using a larger number of Enge coefficients improves the quality of the fit for the field at a given radius (and polar angle θ), but then performing the integral of the quadrupole gradient g(z) to find the function G(ζ) becomes difficult.
It is of course possible to plot an equipotential surface to represent the shape of the pole in an EMMA quadrupole, in the same way that we did for a "normal conducting equivalent" HL-LHC inner triplet quadrupole. However, we find that the shape of the pole is much as expected for a normal conducting quadrupole, i.e. the plot is qualitatively very similar to that shown in Fig. 11.

Conclusions
A closed form analytic expression was presented for fringe fields in multipole magnets. For quadrupoles, the field described by the analytic expression was shown to have the expected properties. The expression can be extended to describe multipoles of any order. For ease of explanation and illustration, we looked in particular at fringe fields in which the multipole gradient had a roll-off along the axis of the magnet described by an Enge function with only a single parameter in the exponent. However, the technique can be applied to any function with the appropriate dependence on the co-ordinates (i.e. any function that depends on the co-ordinates combined in the form √ 2z ± ih). Examples of other (non-Enge) functions that may be suitable for describing fringe fields may be found in [7,8].
Expressions were also given for scalar and vector potentials from which the multipole fringe fields presented here could  Figure 13: Radial component of the magnetic field in an EMMA F quadrupole as a function of position along lines of fixed distance r = r max /4. Each pair of black and red lines shows the field at a different value of the cylindrical polar co-ordinate θ, from π/16 to π/4 (in steps of π/16). The black lines show the field obtained from the (numerical) magnetic model; the red lines show the analytical model (60)-(62), with fit parameters given in Table 2 with radius of fit r max /10, and b 1 = 1.8. The centre of the quadrupole is at the far left of the plot, z = −53.6 mm.  Figure 14: Longitudinal component of the magnetic field in an EMMA F quadrupole as a function of position along lines of given radius and with cylindrical polar co-ordinate θ = π/4. Each pair of black and red lines shows the field at an increasing radius, from r max /5 to 4r max /5, with r max = 36 mm. The black lines show the field obtained from the (numerical) magnetic model; the red lines show the analytical model (60)-(62), with fit parameters given in Table 2 with radius of fit r max /10, and b 1 = 1.8. The centre of the quadrupole is at the far left of the plot, z = −53.6 mm. be derived. Again, the expressions for the potentials can be extended to apply to multipole magnets of any order. The scalar potential is of interest since, in iron-dominated magnets, the pole faces form surfaces of constant scalar potential. This provides a connection between studies of the dynamics of particles moving through the fringe fields a particular magnet, and design studies of the magnet geometry. It is hoped that by having access to realistic analytical descriptions of fringe fields at an early stage in the design of an accelerator beamline, the design process (typically involving many iterations between beam dynamics studies and magnet design work) may be made more efficient.
The vector potential is of interest for particle tracking. In particular, some techniques for symplectic integration of the equations of motion for particles moving in magnetic fields are based on analytical expressions for the vector potential (see, for example, [15]). Again, it is hoped that there will be benefits in being able to perform symplectic tracking through realistic fringe field models at an early stage in the design of an accelerator.
In some types of magnet, such as those used in non-scaling FFAGs, fringe fields dominate the effects of the magnet. In such cases, being able to study the impact of fringe fields at an early stage of the accelerator design is essential for making efficient progress with the design. It should be possible to implement the methods presented here in standard accelerator tracking codes; this will allow accurate modelling of fringe field effects in multipole magnets of arbitrary order, and enhance the range of tools available for accelerator design and simulation.