- Email: [email protected]
Scattering in a magnetic field
ARTICLE IN PRESS
Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
Scattering in a magnetic field David C. Carey* Fermi National Accelerator Laboratory, B.D. Special Projects, MS 221, P.O. Box 500, Batavia, IL 60563, USA
Abstract A scheme of ionization cooling proposed at Fermilab involves energy loss by passing a muon beam through liquid hydrogen. Focussing is done by solenoids. A particle travelling in a straight line, will, after being scattered, follow a helical trajectory. After a full circle is made, the trajectory will be downstream at a point where some of the previously unscattered particles will now scatter. The portions of phase space occupied in the two scatterings now overlap and will occupy less space than if there were no overlap. Results of detailed calculations and approximations made will be discussed. r 2003 Published by Elsevier B.V.
1. Introduction The fixed target program at Fermilab has come to an end. New projects are in the planning stage. Among them is a muon storage ring. Up to the present, all storage rings in high-energy physics have carried stable particles, namely the electron and proton and their antiparticles. The muon is unstable and decays with a mean lifetime of 2:0 106 sec: Two types of cooling have been used in the past. One is stochastic cooling where an electrode is used to detect the positions of the particles and send a signal to another position across the ring. Through successive applications of this technique, the phase space is ultimately greatly reduced and beams can be made to collide with a useful event rate. The second type of cooling is electron cooling. Here protons and electrons are made to travel together for a short distance. Equipartition causes transfer of transverse energy of the protons to that of the electrons. *Corresponding author. Tel.: +1-630-840-3639; fax: +1630-840-6039. E-mail address: [email protected] (D.C. Carey). 0168-9002/$ - see front matter r 2003 Published by Elsevier B.V. doi:10.1016/j.nima.2003.11.186
Neither of these methods is fast enough to allow acceleration of a sufficient number of muons up to maximum energy before they decay. A new method known as ionization cooling has been proposed [1]. The muons are cooled by passing them through a container of liquid hydrogen. The energy loss reduces both transverse and longitudinal momentum. The longitudinal momentum is restored with RF cavities. The net result is to maintain the longitudinal momentum while cooling the transverse momentum. To minimize the total travel distance of the muons the liquid hydrogen is placed inside the focusing solenoids. The question arises as to whether the presence of the solenoids influences the phase space occupied by the muons. After the muon scatters it has transverse momentum. In a constant longitudinal magnetic field the trajectory wraps around the field lines and coincides in momentum and position with a particle which scatters one cycle later. Here we calculate the change in emittance for both a drift space and a solenoid. We find that the presence of the solenoid does cause a reduction in phase space. Shown below are both a derivation of
ARTICLE IN PRESS 462
D.C. Carey / Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
Fig. 1. Scattering of a muon by a slice of liquid hydrogen.
the behavior of the muon phase space and a plot showing the strength of the effect described (Fig. 1).
2. Multiple scattering The particles proceed to the right from a plane marked ‘‘begin.’’ They are scattered by a centrally located slice of material. From the slice of material, the particles continue to the right, finally reaching the plane at the right end of the diagram. The section spanning the distance from the initial plane to the centrally located slice of material has the transfer matrix Ra : The section spanning the distance from the centrally located slice of material to the end on the right has the transfer matrix Rb : The rms scattering angle of a charged particle passing through a section of material medium, in a single transverse plane, is given by the equation [2]: rffiffiffiffiffiffi 13:6 MeV x x y0 ¼ z 1 þ 0:038 log : ð1Þ bcp X0 X0 Here p is the muon momentum and bc is the velocity. The letter z is the charge number of the particle. For a muon, this number is, of course, 1.0. The thickness of the slab of material is x; and X0 is the radiation length of the material. To answer our question about the wrapping up of phase space we shall make a few approximations. Some of these approximations are necessary. However, their effect is small enough so as not to seriously change the answer to the question.
Instead, by simplifying some part of the derivation we add greater clarity to the question about the reduction of phase space. First we ignore any energy loss. The distance a muon can travel in material will be much longer than a wavelength of magnetic focusing in a solenoid. Also, it is difficult to calculate analytically the trajectory of a particle in a solenoid when the magnitude of the momentum is continuously changing. We also ignore any energy dependence of the rms width of the scattering on momentum. Secondly, the value of y0 derived is from a fit to a Moliere distribution. It is not the value that would be obtained from a direct fit to a Gaussian. We delete the term from Eq. (1) containing the logarithm. Then the effect of two scatters would be to add in quadrature, making the mathematics much simpler. Some numerical justification can be obtained from the fact that the radiation length of liquid hydrogen is approximately 9 m: The term containing the logarithm is then about 10% of the remainder of the expression. Finally, we assume that the magnetic field in the solenoid is uniform. As stated above, the equations of motion are simplified, since the projection of the motion on a plane perpendicular to the axis is a circle. If we restrict ourselves to first order, there are no singularities in the end fields. We can integrate right through the ends without any sudden discontinuities.
3. The effect of scattering on the beam matrix From Eq. (1) we can determine that the effect on the beam (sigma) matrix at the point of scatter is given by Eq. (2). The Dz indicates that the change in the transverse second moments is linear in the longitudinal coordinate z: The subscripts x and y represent the coordinates in the transverse configuration space. In order to keep with prevailing usage, we have changed the notation from the previous section. In the previous section, the longitudial coordinate was x: Here z is longitudinal and x and y are transverse. Dsxx ðsÞ ¼ Dsyy ðsÞ ¼ K Dz:
ð2Þ
ARTICLE IN PRESS D.C. Carey / Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
Here we use the argument ‘‘s’’ to indicate the point of scatter. We shall also use the argument ‘‘b’’ and ‘‘e’’ of the sigma matrix to identify the beginning and end of the system respectively. The indices on the transfer matrix R and beam matrix s correspond to a six-dimension vector of kinematic quantities [3]. The vector is ðx; x0 ; y; y0 ; c; dÞ: In our particular application, only the first four components are significant. The symbols c and d indicate the two longitudinal coordinates. They are longitudinal separation and fractional momentum deviation respectively. Transforming the sigma matrix to the end of the system, we get Dsij ðeÞ ¼ K Ri2 ðbÞRj2 ðbÞDz þ KRi4 ðbÞRj4 ðbÞDz: ð3Þ We break the total system up into two steps. The first step (a) is from the beginning of the system to the point of scatter. The second step (b) is from the point of scatter to the end of the system. The total transfer matrix is then RðtÞ ¼ RðbÞRðaÞ
ð4Þ
or RðbÞ ¼ RðtÞRðaÞ1 :
ð5Þ
Substituting for RðbÞ we get X X 1 Dsij ðeÞ ¼ K Rik ðtÞR1 k2 ðaÞRjc ðtÞRc2 ðaÞDz k
þK
c
X X k
Rik ðtÞR1 k4 ðaÞRjc ðtÞ
c
R1 c4 ðaÞDz:
ð6Þ
Integrating Eq. (6) gives the result X X Dsij ðeÞ ¼ Rik ðtÞRjc ðtÞ k
þ
Z Z
c
1 R1 k4 ðaÞRc4 ðaÞ
K dz:
the right side of the equation. Z 1 R1 Dsij ðbÞ ¼ k2 ðaÞRc2 ðaÞ Z 1 1 þ Rk4 ðaÞRc4 ðaÞ K dz:
ð7Þ
This equation gives the components of the beam matrix at the end of the system. The equation for the beam matrix at the beginning of the system is obtained by omitting the two summations and the two factors of RðtÞ on
ð8Þ
Since we will eventually be taking determinants, this omission will not change the answer.
4. Factorized transfer matrices Taking the inverse or determinant of a 4 4 matrix is a procedure fraught with peril, and guaranteed to produce more errors than correct calculations. Fortunately in this case the total transfer matrix factors into a product of matrices. Each of these matrices is much simpler than the product matrix. From the analytic form of the product of matrices, we can see that one matrix represents a geometric rotation of two planes [4,5]. The rotating frame is known as the Larmor frame. The dynamics all take place in the rotating plane. The motions in the two transverse planes are independent. They consist of harmonic oscillations in each of the two planes. If the magnetic field is not uniform the dynamic motion in the rotated plane is a general matrix with cosinelike and sinelike functions executing a harmonic motion of changing period. Since the Larmor frame is a geometric rotation, the components of the transformation matrix can be given directly in terms of trigonometric functions. Non-uniform field: 0
1 R1 k2 ðaÞRc2 ðaÞ
463
cu B c0 B R ¼B u @0
sy s0u
0 0
0
cv
1 0 0C C C sv A
0 0
0
c0v
s0v
cos b B 0 B B @ sin b 0
0 cos b
sin b 0
0 sin b
cos b 0
1 0 sin b C C C: 0 A cos b
ð9Þ
ARTICLE IN PRESS D.C. Carey / Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
464
Uniform field: 0 1 sin kz B cos kz k B B k sin kz cos kz B R ¼B B B 0 0 @ 0 0 0 cos kz 0 B 0 cos kz B B @ sin kz 0 sin kz
0
Inverse matrix: 0 cos kz 0 B 0 cos kz B R1 ¼ B @ sin kz 0 0 sin kz 0
For a single transverse plane 0
0
1
C C 0 0 C C C 1 C sin kz C cos kz A k k sin kz cos kz 1 sin kz 0 0 sin kz C C C: ð10Þ cos kz 0 A cos kz
0
sin kz 0 cos kz 0
1 0 sin kz C C C 0 A cos kz
1 0 B cos kz k sin kz B B k sin kz cos kz 0 B B B B 0 0 cos kz @ 0 0 k sin kz
0
1 4 L : Determinant ¼ 12
ð14Þ
6. Application to solenoids Evaluating the matrix for a solenoid takes a little more effort. Here we let the program MATHEMATICA [6] do the manipulation. We then find the result 1 CS B 2k2 L k B B S2 B B 2 Z B 2k sum ¼ B B B 0 B B B @ 0 0
S2 2k2 1 CS Lþ 2 k
0 0
0 0 1 CS L 2k2 k
S2 2k2
0
1
C C C C C 0 C C; C S2 C 2 C C 2k C 1 CS A Lþ 2 k
ð15Þ 1
C C C 0 C C: 1 C sin kz C A k cos kz
ð11Þ
where C ¼ cos kL S ¼ sin kL qB k¼ : 2p For a single transverse plane S2 2 2 Determinant ¼ 1=4k L 2 : k
ð16Þ
ð17Þ
5. The integrated product of inverse matrices 7. Ratio of emittances 5.1. Drift spaces If we integrate the equation for a drift space, we find the results: 0 1 3 1 12 L2 0 0 3L Z B 1 L2 L 0 0 C B C sum ¼ B 2 ð12Þ 1 3 1 2C @ 0 0 2 L A 3L 0 1 Determinant ¼ 144 L8 :
0
12 L2
L ð13Þ
We drop the factor K as we shall be taking the ratio of determinants and factors of K will cancel.
The horizontal and vertical emittances in terms of elements of the beam matrix are given by the relations e2H ¼ s11 s22 s221
ð18Þ
e2V ¼ s33 s44 s243 : From this we can derive the ratio of emittances in a single transverse plane: eS for a solenoid, and eD for a drift (Fig. 2). sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi eS 3 S2 ¼ 1 Ratio ¼ : ð19Þ eD f2 f2
ARTICLE IN PRESS D.C. Carey / Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
465
Acknowledgements RATIO OF EMITTANCES
1.00
The author wishes to express appreciation for suggestions and advice to Valeri Belbekov, Mike Syphers, Michael Kriss, Norman Gelfand, and Paul LeBrun.
0.75
0.50
References 0.25
0.00 0
2
4
6
8
10
PHI
Fig. 2. The ratio of horizontal emittances for solenoids and drift spaces.
The Larmor phase f is given by f ¼ kL: and S is given in (16).
ð20Þ
[1] Ankenbrandt et al., Ionization cooling research and development program for a high luminosity muon collider, Proposal, April 1998. [2] Eur. Phy. J. 15 (2000) 1. [3] D.C. Carey, The Optics of Charged Particle Beams, Harwood Academic Publishers, New York, 1987. [4] M. Reiser, Theory and Design of Charged Particle Beams, Wiley, New York, 1994. [5] K. Halbach, private communication. [6] Wolfram, Stephen, The MATHEMATICA BOOK, Wolfram Media and Cambridge University Press, Cambridge, 1991.
Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
Scattering in a magnetic field David C. Carey* Fermi National Accelerator Laboratory, B.D. Special Projects, MS 221, P.O. Box 500, Batavia, IL 60563, USA
Abstract A scheme of ionization cooling proposed at Fermilab involves energy loss by passing a muon beam through liquid hydrogen. Focussing is done by solenoids. A particle travelling in a straight line, will, after being scattered, follow a helical trajectory. After a full circle is made, the trajectory will be downstream at a point where some of the previously unscattered particles will now scatter. The portions of phase space occupied in the two scatterings now overlap and will occupy less space than if there were no overlap. Results of detailed calculations and approximations made will be discussed. r 2003 Published by Elsevier B.V.
1. Introduction The fixed target program at Fermilab has come to an end. New projects are in the planning stage. Among them is a muon storage ring. Up to the present, all storage rings in high-energy physics have carried stable particles, namely the electron and proton and their antiparticles. The muon is unstable and decays with a mean lifetime of 2:0 106 sec: Two types of cooling have been used in the past. One is stochastic cooling where an electrode is used to detect the positions of the particles and send a signal to another position across the ring. Through successive applications of this technique, the phase space is ultimately greatly reduced and beams can be made to collide with a useful event rate. The second type of cooling is electron cooling. Here protons and electrons are made to travel together for a short distance. Equipartition causes transfer of transverse energy of the protons to that of the electrons. *Corresponding author. Tel.: +1-630-840-3639; fax: +1630-840-6039. E-mail address: [email protected] (D.C. Carey). 0168-9002/$ - see front matter r 2003 Published by Elsevier B.V. doi:10.1016/j.nima.2003.11.186
Neither of these methods is fast enough to allow acceleration of a sufficient number of muons up to maximum energy before they decay. A new method known as ionization cooling has been proposed [1]. The muons are cooled by passing them through a container of liquid hydrogen. The energy loss reduces both transverse and longitudinal momentum. The longitudinal momentum is restored with RF cavities. The net result is to maintain the longitudinal momentum while cooling the transverse momentum. To minimize the total travel distance of the muons the liquid hydrogen is placed inside the focusing solenoids. The question arises as to whether the presence of the solenoids influences the phase space occupied by the muons. After the muon scatters it has transverse momentum. In a constant longitudinal magnetic field the trajectory wraps around the field lines and coincides in momentum and position with a particle which scatters one cycle later. Here we calculate the change in emittance for both a drift space and a solenoid. We find that the presence of the solenoid does cause a reduction in phase space. Shown below are both a derivation of
ARTICLE IN PRESS 462
D.C. Carey / Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
Fig. 1. Scattering of a muon by a slice of liquid hydrogen.
the behavior of the muon phase space and a plot showing the strength of the effect described (Fig. 1).
2. Multiple scattering The particles proceed to the right from a plane marked ‘‘begin.’’ They are scattered by a centrally located slice of material. From the slice of material, the particles continue to the right, finally reaching the plane at the right end of the diagram. The section spanning the distance from the initial plane to the centrally located slice of material has the transfer matrix Ra : The section spanning the distance from the centrally located slice of material to the end on the right has the transfer matrix Rb : The rms scattering angle of a charged particle passing through a section of material medium, in a single transverse plane, is given by the equation [2]: rffiffiffiffiffiffi 13:6 MeV x x y0 ¼ z 1 þ 0:038 log : ð1Þ bcp X0 X0 Here p is the muon momentum and bc is the velocity. The letter z is the charge number of the particle. For a muon, this number is, of course, 1.0. The thickness of the slab of material is x; and X0 is the radiation length of the material. To answer our question about the wrapping up of phase space we shall make a few approximations. Some of these approximations are necessary. However, their effect is small enough so as not to seriously change the answer to the question.
Instead, by simplifying some part of the derivation we add greater clarity to the question about the reduction of phase space. First we ignore any energy loss. The distance a muon can travel in material will be much longer than a wavelength of magnetic focusing in a solenoid. Also, it is difficult to calculate analytically the trajectory of a particle in a solenoid when the magnitude of the momentum is continuously changing. We also ignore any energy dependence of the rms width of the scattering on momentum. Secondly, the value of y0 derived is from a fit to a Moliere distribution. It is not the value that would be obtained from a direct fit to a Gaussian. We delete the term from Eq. (1) containing the logarithm. Then the effect of two scatters would be to add in quadrature, making the mathematics much simpler. Some numerical justification can be obtained from the fact that the radiation length of liquid hydrogen is approximately 9 m: The term containing the logarithm is then about 10% of the remainder of the expression. Finally, we assume that the magnetic field in the solenoid is uniform. As stated above, the equations of motion are simplified, since the projection of the motion on a plane perpendicular to the axis is a circle. If we restrict ourselves to first order, there are no singularities in the end fields. We can integrate right through the ends without any sudden discontinuities.
3. The effect of scattering on the beam matrix From Eq. (1) we can determine that the effect on the beam (sigma) matrix at the point of scatter is given by Eq. (2). The Dz indicates that the change in the transverse second moments is linear in the longitudinal coordinate z: The subscripts x and y represent the coordinates in the transverse configuration space. In order to keep with prevailing usage, we have changed the notation from the previous section. In the previous section, the longitudial coordinate was x: Here z is longitudinal and x and y are transverse. Dsxx ðsÞ ¼ Dsyy ðsÞ ¼ K Dz:
ð2Þ
ARTICLE IN PRESS D.C. Carey / Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
Here we use the argument ‘‘s’’ to indicate the point of scatter. We shall also use the argument ‘‘b’’ and ‘‘e’’ of the sigma matrix to identify the beginning and end of the system respectively. The indices on the transfer matrix R and beam matrix s correspond to a six-dimension vector of kinematic quantities [3]. The vector is ðx; x0 ; y; y0 ; c; dÞ: In our particular application, only the first four components are significant. The symbols c and d indicate the two longitudinal coordinates. They are longitudinal separation and fractional momentum deviation respectively. Transforming the sigma matrix to the end of the system, we get Dsij ðeÞ ¼ K Ri2 ðbÞRj2 ðbÞDz þ KRi4 ðbÞRj4 ðbÞDz: ð3Þ We break the total system up into two steps. The first step (a) is from the beginning of the system to the point of scatter. The second step (b) is from the point of scatter to the end of the system. The total transfer matrix is then RðtÞ ¼ RðbÞRðaÞ
ð4Þ
or RðbÞ ¼ RðtÞRðaÞ1 :
ð5Þ
Substituting for RðbÞ we get X X 1 Dsij ðeÞ ¼ K Rik ðtÞR1 k2 ðaÞRjc ðtÞRc2 ðaÞDz k
þK
c
X X k
Rik ðtÞR1 k4 ðaÞRjc ðtÞ
c
R1 c4 ðaÞDz:
ð6Þ
Integrating Eq. (6) gives the result X X Dsij ðeÞ ¼ Rik ðtÞRjc ðtÞ k
þ
Z Z
c
1 R1 k4 ðaÞRc4 ðaÞ
K dz:
the right side of the equation. Z 1 R1 Dsij ðbÞ ¼ k2 ðaÞRc2 ðaÞ Z 1 1 þ Rk4 ðaÞRc4 ðaÞ K dz:
ð7Þ
This equation gives the components of the beam matrix at the end of the system. The equation for the beam matrix at the beginning of the system is obtained by omitting the two summations and the two factors of RðtÞ on
ð8Þ
Since we will eventually be taking determinants, this omission will not change the answer.
4. Factorized transfer matrices Taking the inverse or determinant of a 4 4 matrix is a procedure fraught with peril, and guaranteed to produce more errors than correct calculations. Fortunately in this case the total transfer matrix factors into a product of matrices. Each of these matrices is much simpler than the product matrix. From the analytic form of the product of matrices, we can see that one matrix represents a geometric rotation of two planes [4,5]. The rotating frame is known as the Larmor frame. The dynamics all take place in the rotating plane. The motions in the two transverse planes are independent. They consist of harmonic oscillations in each of the two planes. If the magnetic field is not uniform the dynamic motion in the rotated plane is a general matrix with cosinelike and sinelike functions executing a harmonic motion of changing period. Since the Larmor frame is a geometric rotation, the components of the transformation matrix can be given directly in terms of trigonometric functions. Non-uniform field: 0
1 R1 k2 ðaÞRc2 ðaÞ
463
cu B c0 B R ¼B u @0
sy s0u
0 0
0
cv
1 0 0C C C sv A
0 0
0
c0v
s0v
cos b B 0 B B @ sin b 0
0 cos b
sin b 0
0 sin b
cos b 0
1 0 sin b C C C: 0 A cos b
ð9Þ
ARTICLE IN PRESS D.C. Carey / Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
464
Uniform field: 0 1 sin kz B cos kz k B B k sin kz cos kz B R ¼B B B 0 0 @ 0 0 0 cos kz 0 B 0 cos kz B B @ sin kz 0 sin kz
0
Inverse matrix: 0 cos kz 0 B 0 cos kz B R1 ¼ B @ sin kz 0 0 sin kz 0
For a single transverse plane 0
0
1
C C 0 0 C C C 1 C sin kz C cos kz A k k sin kz cos kz 1 sin kz 0 0 sin kz C C C: ð10Þ cos kz 0 A cos kz
0
sin kz 0 cos kz 0
1 0 sin kz C C C 0 A cos kz
1 0 B cos kz k sin kz B B k sin kz cos kz 0 B B B B 0 0 cos kz @ 0 0 k sin kz
0
1 4 L : Determinant ¼ 12
ð14Þ
6. Application to solenoids Evaluating the matrix for a solenoid takes a little more effort. Here we let the program MATHEMATICA [6] do the manipulation. We then find the result 1 CS B 2k2 L k B B S2 B B 2 Z B 2k sum ¼ B B B 0 B B B @ 0 0
S2 2k2 1 CS Lþ 2 k
0 0
0 0 1 CS L 2k2 k
S2 2k2
0
1
C C C C C 0 C C; C S2 C 2 C C 2k C 1 CS A Lþ 2 k
ð15Þ 1
C C C 0 C C: 1 C sin kz C A k cos kz
ð11Þ
where C ¼ cos kL S ¼ sin kL qB k¼ : 2p For a single transverse plane S2 2 2 Determinant ¼ 1=4k L 2 : k
ð16Þ
ð17Þ
5. The integrated product of inverse matrices 7. Ratio of emittances 5.1. Drift spaces If we integrate the equation for a drift space, we find the results: 0 1 3 1 12 L2 0 0 3L Z B 1 L2 L 0 0 C B C sum ¼ B 2 ð12Þ 1 3 1 2C @ 0 0 2 L A 3L 0 1 Determinant ¼ 144 L8 :
0
12 L2
L ð13Þ
We drop the factor K as we shall be taking the ratio of determinants and factors of K will cancel.
The horizontal and vertical emittances in terms of elements of the beam matrix are given by the relations e2H ¼ s11 s22 s221
ð18Þ
e2V ¼ s33 s44 s243 : From this we can derive the ratio of emittances in a single transverse plane: eS for a solenoid, and eD for a drift (Fig. 2). sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi eS 3 S2 ¼ 1 Ratio ¼ : ð19Þ eD f2 f2
ARTICLE IN PRESS D.C. Carey / Nuclear Instruments and Methods in Physics Research A 519 (2004) 461–465
465
Acknowledgements RATIO OF EMITTANCES
1.00
The author wishes to express appreciation for suggestions and advice to Valeri Belbekov, Mike Syphers, Michael Kriss, Norman Gelfand, and Paul LeBrun.
0.75
0.50
References 0.25
0.00 0
2
4
6
8
10
PHI
Fig. 2. The ratio of horizontal emittances for solenoids and drift spaces.
The Larmor phase f is given by f ¼ kL: and S is given in (16).
ð20Þ
[1] Ankenbrandt et al., Ionization cooling research and development program for a high luminosity muon collider, Proposal, April 1998. [2] Eur. Phy. J. 15 (2000) 1. [3] D.C. Carey, The Optics of Charged Particle Beams, Harwood Academic Publishers, New York, 1987. [4] M. Reiser, Theory and Design of Charged Particle Beams, Wiley, New York, 1994. [5] K. Halbach, private communication. [6] Wolfram, Stephen, The MATHEMATICA BOOK, Wolfram Media and Cambridge University Press, Cambridge, 1991.