- Email: [email protected]
Tracking protons in an FFAG and a synchrotron with tiny arcs
Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
Contents lists available at SciVerse ScienceDirect
Nuclear Instruments and Methods in Physics Research A journal homepage: www.elsevier.com/locate/nima
Tracking protons in an FFAG and a synchrotron with tiny arcs$ K.M. Hock a,b,n, C.S. Edmonds a,b a b
Department of Physics, University of Liverpool, Liverpool L69 7ZE, United Kingdom Cockcroft Institute, Daresbury Science and Innovation Campus, Daresbury, Warrington WA4 4AD, United Kingdom
a r t i c l e i n f o
a b s t r a c t
Article history: Received 1 August 2011 Accepted 26 September 2011 Available online 9 November 2011
We study a ray tracing method that uses tiny arcs to track a proton through a magnetic field. We show how to apply this to an FFAG and a synchrotron that are designed for cancer therapy. Assuming that the field traversed by a proton is uniform over each time step, the required equations and a simple scheme for treating a lattice are developed. We compare the computation time with the code Zgoubi for a FODO lattice. To test its usefulness, it is applied to the spiral FFAG, designed in the RACCAM project, to find the horizontal phase space trajectories. It is also applied to a compact synchrotron, designed in Japan, to determine the phase space trajectories during slow extraction of the proton beam. & 2011 Elsevier B.V. All rights reserved.
Keywords: FFAG Synchrotron Proton therapy Tracking Ray tracing
1. Introduction This work is motivated by the need for a simple and efficient way to track protons in a cancer therapy machine. Tracking of charged particles is an important part of the design of an accelerator. It involves calculating the path of the particles to see if they are still stable after many turns. There are many methods for doing so. These can be broadly classified into two maingroups—transfer maps and ray tracing. A transfer map is a way to quickly compute the position and velocity at some point along a beamline, given the co-ordinates at a starting point. This is done without the need to compute the positions of the trajectory in between. This map could make use of a matrix, a power series, or Lie algebra [1]. It would normally involve some approximations, such as assuming that the map is linear, or truncating the power series after a few orders. These make the calculations very fast, but usually mean that the results are accurate only for small displacements from the reference orbit. Many codes have been written using these methods [2]. They are particularly useful for high energy accelerators, where the particles are expected to stay close to the reference orbit. At low energies, or in more complex elements of the accelerator, trajectories could either have large curvatures or move sideways by large distances. For accurate tracking, ray tracing has to be used. This involves integrating the equation of motion of the particle in an electromagnetic field in some way. This means that
$
Work supported by the Science and Technology Facilities Council, UK. Corresponding author at: Department of Physics, University of Liverpool, Liverpool L69 7ZE, United Kingdom. E-mail address: [email protected] (K.M. Hock). n
0168-9002/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2011.09.066
unlike transfer maps, the co-ordinates between any two points on a trajectory must be computed. The methods used include the Runge–Kutta method [3], power series expansion of the equations of motion [4], and tiny arcs [5]. The time or distance steps can be reduced until sufficient accuracy is achieved. For long trajectories, very small steps are required, which could be time consuming. In this paper, we develop on the method that uses tiny arcs. The field traversed by a proton in each time step is assumed to be uniform. The trajectory in this time step is therefore an arc. This method has been developed by another group for tracking in cyclotrons [6]. We modify this to make it more efficient for accelerators that are made up of lumped elements with large drift spaces in between. In Section 2, we develop the algorithm and test the code on a FODO lattice. In Section 3, we apply this to an FFAG and a synchrotron. In Section 4, we summarise the results and suggest further work.
2. The arc method A charged particle travelling through a magnetic field experiences a force given by F ¼ qv B
ð1Þ
where q is the charge of the particle, v is the velocity and B the magnetic field. To determine the path through a nonuniform field, numerical methods are needed. The code Zgoubi [7], for example, uses Taylor expansions up to a few orders to approximate the path in each time step. In these notes, I consider using arcs instead. As I intend to develop this method for the study of the proton accelerator, I shall take proton as the charge particle.
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
Consider a proton travelling through a nonuniform field. In a very small time step, the field over which the proton would move would be approximately constant. If the velocity is perpendicular to the field, the path would be a perfectly circular arc. If the velocity has a component parallel to the field, the arc also gains a component in the direction of the field. The arc becomes a spiral arc. Assuming that the velocity is approximately constant in this short time interval, we would still be able to calculate the endpoint of this spiral arc. The idea of the arc method is to start with the new field at the endpoint of this arc, and find the next arc. This is repeated, and the whole trajectory is obtained by joining many arcs together. After that, it is important to repeat the calculations for a number of smaller time steps to check that the answer has converged. From the starting position x0 , we need a formula to calculate the endpoint x1 of the arc after a time interval dt, in terms of B and v0 , the velocity at the starting position. We first consider only the perpendicular component v? and find the circular arc s, as shown in Fig. 1. Then we add the displacement due to the parallel component vJ . To do these, we first obtain the following. The constant magnetic field should ideally be taken as the field at the midpoint of the arc. However, we do not know where this midpoint is before finding the arc. A simple approximation would be to go in a straight line from the starting position for half a time step to the point x ¼ x0 þv0 dt=2
ð2Þ
and take the magnetic field here as the constant field B. The unit vector in the direction of the field is B B^ ¼ : B
ð3Þ
The parallel component of the velocity is ^ B: ^ vJ ¼ ðv0 :BÞ
ð4Þ
The perpendicular component is then v? ¼ v0 vJ :
ð5Þ
The radius of the circular arc is given by r¼
gmv?
ð6Þ
qB
where m is the rest mass and g is the Lorentz factor. We first consider the component of the motion that is perpendicular to the field. The circular arc starts off from x0 , goes in the direction v? , and curves towards F —since F points to the centre of the arc. After time dt, the angle y that the arc subtends at the centre is
y¼
v? d t ¼ odt r
ð7Þ
where we have also defined o. To express the endpoint of the arc in terms of these, we use the following unit vectors. The unit
141
vector in the direction of v? is given by v^ ? ¼
v? v?
ð8Þ
and the unit vector in the direction of F is given by F F^ ¼ : F
ð9Þ
Therefore, the endpoint of the circular arc is x1? ¼ x0 þðrr cos yÞF^ þ r sin yv^ ? :
ð10Þ
Finally, we add the displacement parallel to the field to get the endpoint of the spiral arc, x1 ¼ x0 þ ðrr cos yÞF^ þ ðr sin yÞv^ ? þvJ dt:
ð11Þ
Differentiating this with respect to time, we obtain the velocity at the endpoint: v1 ¼ ðr o sin yÞF^ þ ðr o cos yÞv^ ? þvJ :
ð12Þ
We shall call the last two equations the arc formula. To use the arc formula to develop a computer code for tracking a proton in an accelerator with lumped elements, we assume the following: 1. The beamline is made up of elements of a finite length, containing magnetic fields. 2. A proton would usually enter from one side of an element, and exit from the opposite side. 3. In between elements are drift spaces, with no field. Elements with electric fields are not considered in this paper. With the above assumptions, we can always define two imaginary planes for each element, one on each side. For example, the planes could be the entrance and exit faces of a quadrupole or dipole. These planes could be roughly perpendicular to the beamline. A proton would enter one plane, go through the element, and come out through the other plane. In between the two planes, there would be a magnetic field. The field can be of any distribution. The two planes for each element can be used to define the boundaries of the element in the code. The arc formula can then be used to compute the trajectory of the proton within these boundaries. There is one problem. The arc formula used a fixed time step dt. When the trajectory formed by the arcs crosses the plane, the field could drop abruptly to zero just outside the plane. This happens if we use the hard edge models for a dipole or a quadrupole. Since the final arc in the element would in general not terminate exactly on the exit plane, the arc is partly inside the field and partly outside. This is different from the original assumption that the field is approximately constant over the extent of the arc, and would lead to inaccuracy. To avoid this problem, we can make the final arc terminate exactly on the exit plane by applying the Newton– Raphson method. The formula is given by t1 ¼ t0
f ðt 0 Þ 0 f ðt 0 Þ
ð13Þ
where f(t) is a function that crosses zero for some t. By applying the formula iteratively—computing t1, replacing t0 by this, computing t1 again—t1 would get closer and closer to the root. We can apply this by taking f(t) as the distance from a point on the final arc to the exit plane, and renaming dt as t. To find f(t), rewrite the arc formula in this form: x1 ¼ x0 þ ðrr cos otÞF^ þ ðr sin otÞv^ ? þ vJ t: Fig. 1. Trajectory in a uniform field.
ð14Þ
Next, define the exit plane by a point on the plane p2 , and the unit vector normal to the plane n2 . The perpendicular displacement of
142
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
a point on the arc from the plane is then f ðtÞ ¼ ðx1 p2 Þ:n2 :
ð15Þ
Differentiating this gives the derivative: 0 f ðtÞ ¼ ½ðv? sin otÞF^ þ ðv? cos otÞv^ ? þ vJ :n2 :
ð16Þ
Since dt tends to be small, the start point of the final arc is already quite close to the exit plane. So a few iterations is often enough to achieve good accuracy. In our simulations, we stop when the change in the next iteration is within 10 5 of dt. The value of t obtained is then used in place of dt in the arc formula for the final arc before the exit plane. A C code has been written based on the above descriptions. In the first test, we compute the path of a proton through a 1 T uniform field. The result is shown in Fig. 2. This is a trivial case, since the arc formula is exact for the uniform field, regardless of the dt. For this figure, I choose 1000 time steps, just to make the curve look smooth. It is only to check that the code is working. The spiral path is produced as expected. In the next test, we use a FODO lattice with 100 cells. The same 10 MeV proton is used. Each quadrupole is 10 cm long, and has a field gradient of 10 T/m. The distance between adjacent quadrupoles is 0.2 m. The trajectory computed using the arc formula is shown in Fig. 3. We check for convergence by looking at the y coordinate of the endpoint of the final arc, which falls on the exit face of the final quadrupole in the lattice. We start with a step size
of a few centimetres, carry out the tracking, record the y value and the computational time. The step size is then halved, and the tracking repeated. Each time the step size is halved, the difference between the current and previous y values is calculated. Assuming that the calculation converges, this difference would decrease as the step size decreases. We define the tolerance to be the fraction of this difference over the current y value. When the tolerance falls below a value that is small enough for our purpose, we can say that the y value, and therefore the tracking, has converged. Fig. 4 shows the tolerance, plotted against step size, for the arc method and for Zgoubi. Notice that to achieve the same tolerance, the step size required for the arc method is smaller than the step size required for Zgoubi. Fig. 5 shows the computation time plotted against tolerance. Our computer has a Pentium 3 GHz processor. The trends for the two methods are quite clear. The times would be different on a different computer, but the relative time between the two methods should remain the same, as this depends mainly on the amount of computations needed. For higher tolerance, the computation time for the arc method is shorter. There is a tolerance below which the time for Zgoubi becomes shorter.
101 100 10-1
Zgoubi Arc
x (m)
Tolerance
10-2 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4
10-3 10-4 10-5 10-6 10-7
0
0
2
4
z (m)
6
10-8
0.1
-0.1 -0.2 -0.3 -0.4 m) -0.5 y( -0.6 10
8
10-9 10-3
10-2
10-1
100
101
Step size (cm) Fig. 4. The fractional change in the endpoint y coordinate each time the step size is halved and is plotted against step size.
Fig. 2. Trajectory in a uniform field.
0.02
100 Zgoubi
0.015
Arc 10-1 Compute time (s)
0.01
y (m)
0.005 0 -0.005 -0.01
10-2
10-3
-0.015 -0.02 0
10
20
30
40
50
60
z (m) Fig. 3. Trajectory through 100 FODO cells. It crosses the quadrupoles at the thickened parts of the curve.
10-4 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 101 Tolerance Fig. 5. The computation time for tracking through 100 FODO cell is plotted against the tolerance required.
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
Tolerance is the fractional change in the result from one step size to the next. It is not the same as accuracy, which should be a measure of the error of the result compared to the correct answer. Fig. 6 shows the y coordinate of the endpoint of the trajectory plotted against the tolerance for both methods. The results of both methods appear to approach the same limiting value. Since we do not know the correct answer, we assume that this is the average of the values obtained using the smallest step sizes in the two methods. We then define the error for each step size, as the difference between the y value at the endpoint calculated by each method, and this nominal answer. Fig. 7 shows a graph of the computation time against the absolute error of each method. Note that the graphs would be less reliable for very small errors. Since our nominal answer is only a guess, the graphs are likely to be less accurate around the smallest step sizes. The graphs for the 100 FODO cells show that the arc method is faster if the required error is above 1 mm. Below that, Zgoubi is faster. Further study would be needed to determine the usefulness of the arc method in accelerator design, particularly when tracking for many turns round a ring is required. An issue that is often a concern in tracking methods is symplecticity [8]. Symplecticity is a condition on the transfer map that follows from a Hamiltonian system. Making a tracking
Zgoubi
-1.2986
143
method symplectic removes one source of error. One example is to approximate each tracking step using a drift, a kick, and a drift [8], so that each step is symplectic. Another example is to use Lie algebra [9] to create symplectic maps. A symplectic tracking method is still an approximation, of course. There could be other sources of errors that may not be small and have to be carefully checked. In the arc method, we see that within each time step, the tracking is symplectic. This is because the arc is the exact solution, once we assume that the magnetic field is constant over this step. So the joining of tiny arcs is in essence a concatenation of symplectic maps, which must also be symplectic. Although this is interesting to know, in practice, the arc is still an approximation to the real path and the arc length has to be made small enough to give good accuracy. Another issue that could be of interest is time reversibility. After tracking through a lattice, if we reverse the particle velocity, reverse all the magnetic fields, and track using the same method, would the particle arrive back at the same starting point? The answer is no. The reason is that at each time step, the field at the beginning of the step would in general be different from the field at the end of the step. If time is reversed, the end would become the start. Then over the same step, the field could be slightly different. This difference could be made small by reducing the time step. For fairly symmetrical magnetic fields, we may expect that the error in the time reversed tracking is comparable to the error in the forward tracking.
Arc Endpoint y (cm)
3. Results and discussion
-1.2988 3.1. Spiral FFAG
-1.299
We start with the tracking of protons in the spiral FFAG. This spiral FFAG is designed in the RACCAM project [10–12] for proton therapy. The lattice we simulated is shown in Fig. 8. It follows the parameters in Ref. [11]. Those relevant to our simulation are summarised in Table 1. The number of cells refer to the number of spiral magnets, including the associated drift spaces. The cells are identical.
-1.2992
-1.2994 10-7
10-6
10-5 10-4 Tolerance
10-3
4
10-2
3
Fig. 6. The y coordinate at the end of the tracking is plotted against the tolerance.
magnet
180 MeV 60 MeV
proton
15 MeV 2
100
(X,Y) Zgoubi
1
r
marker
Compute time (s)
10
Y (m)
Arc -1
0
β r1
-1
10-2
θ
r2
A
-2 10-3
-3 entrance plane exit plane
fringe edges 10-4 10-7
-4 10-6
10-5
10-4 Error (cm)
10-3
10-2
10-1
Fig. 7. The computation time is plotted against the absolute error for each method.
-4
-3
-2
-1
0 X (m)
1
2
Fig. 8. Lattice and closed orbits in the spiral FFAG.
3
4
144
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
2
Table 1 Parameters for the spiral FFAG [11]. N p k
z r1 r2 B2
10 0.34 5 53.71 2.7 m 3.5 m 1.8 T
Imagine drawing a circle that is centred at the centre of the accelerator, and that passes through the magnets. The packing factor is the fraction of the circle that is covered by the magnets. This fraction is the same for different radii, as long as the circle passes through the magnets. The field index k is the variation of the magnetic field with the radius r. This variation is given by
1.5
Bz (T)
Number of cells Packing factor Field index Spiral angle Minimum radius Maximum radius Maximum field
2.7 m 3.1 m 3.5 m
1
0.5
0 -30
Bz ¼ br
k
where Bz is the vertical component in the mid-plane in the magnet region. Ideally, the horizontal component at the midplane is zero. This is also what we assume in the simulation. The constant of proportion b may be obtained using the parameters in Table 1 that Bz is 1.8 T when r is 3.5 m. Using the equation and the parameters in this way gives us a simplified model of the measured field data in Ref. [12]. A sophisticated model is given in Ref. [11]. The maximum and minimum radii in Table 1 indicate the region where we expect the protons to travel. The actual size of the magnet is larger. The spiral angle is the angle between a radius and a point on a spiral edge of a magnet. Each magnet has four curved edges. Two of these are circular arcs. These are the one nearest to the centre of the accelerator, and the one furthest away. The other two are spiral arcs. The radius refer to a straight line from the centre of the accelerator. Imagine drawing a radius to a point on any one of the spiral arcs. Then this straight line would always make an angle of 53.71 with the tangent to the spiral arc, at the point where they meet. This is used as a condition to define the exact shape of the spiral arc. From this condition, the equation of the spiral arc can be derived as r ¼ r 0 exp
y
tan z
ð18Þ
where z is the 53.71 spiral angle. The parameter r0 would be different for each of the spiral edges. A prototype of the magnet has been constructed and the field measured [12]. A significant part of the field is fringe field. We use a simple model for this field. Imagine again a circle through the magnets. If we move into the field of a magnet along this circle, the field increases, stay the same, then decreases again. The increase and decrease are the fringe fields. In our simulation, we simplify the field using a linear variation for the fringe field along the circle. To describe this variation, we may define an angular function FðyÞ. The overall field for a magnet may then be given by k
-25
ð17Þ
Bz ¼ br Fðyy0 Þ:
ð19Þ
The resulting field is shown in Fig. 9 for the three radial distances 2.7 m, 3.1 m and 3.5 m. These correspond roughly to the positions of the closed orbits at 15 MeV, 60 MeV and 180 MeV respectively—see Fig. 8. To understand what y0 is, consider the magnet just below the ‘marker’ label in Fig. 8. Imagine a spiral curve that is mid way between the two spiral edges, and that is described by the same Eq. (18) as the edges. The point on this curve closest to the centre of the ring is at distance r1 from the centre. Suppose y is zero
-20
-15
-10 θ (deg)
-5
0
5
10
Fig. 9. Field model for the spiral FFAG magnet.
there. The curve would be given by y r ¼ r 1 exp : tan z
ð20Þ
The y in this equation is then equal to y0 for each value of r. In this way, the angular function Fðyy0 Þ is made to follow the spiral. Using the same angular function for different r is possible because the length of the circular arc through the magnet is always equal to A. This is true for any r, since one spiral edge can be obtained from the other by a rigid rotation about the centre of the ring. However, the actual extent of the fringe fields change at different r, so using the same angular function is an approximation. The field for each of the N magnets is then given by k
Bz ¼ br Fðyy0 nbÞ
ð21Þ
where n ¼ 0; 1, . . . ,N1. This completes our simple model of the magnetic field in the whole lattice. We should emphasise that we make this model just to demonstrate the arc method. The accurate model in Ref. [11] can be used instead, as can any field map. Next, we need to specify the entrance and exit planes for each magnet. These are shown by the two straight lines next to each magnet in Fig. 8. These planes have to be just beyond the fringe fields, so that the full extent of the field would be sampled during tracking of the proton. The planes should not be too far from the fringe fields, so as to give more drift distance between the exit plane of one magnet and the entrance plane of the next. This is because tracking over this distance requires little computation. The most direct way to determine the positions of the entrance and exit planes is by accurately plotting both their positions and the boundaries of the fringe fields on the same drawing, and checking them by eye. For our simulation, we plotted these as lines and curves on a graph. The curves labelled ‘fringe edges’ in Fig. 8 are the boundaries of the fringe fields. We typed in the positions of the planes by hand, replotted the graphs, and tried this a few times until we got a reasonable position for each plane. If the positions of the planes are correct, we would be able to track a proton though the ring and get a trajectory that looks like one of the orbits shown in Fig. 8. In practice, we find it useful to make the code track the proton through one magnet at a time, each time plotting out the trajectory with the magnet outlines, as in Fig. 8, to check. We start by finding the closed orbit at 60 MeV, using a method described in Ref. [7]. We start near the middle of the plane labelled ‘marker’ in Fig. 8, with a velocity in the direction of the tangent to
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
3.2. Compact synchrotron A compact synchrotron is designed in Japan for cancer therapy [13–15]. The lattice is outlined in Fig. 11. It makes use of six combined function magnets. The fields in these magnets contain dipole and quadrupole components. Each magnet is an F–D–F triplet. Since these are linear, tracking with transfer maps are usually sufficient for trajectories close to the reference orbit. We consider here the extraction region of the orbit, where large displacements from the reference are expected. The relevant parameters for the synchrotron lattice are listed in Table 2. Each F–D–F triplet forms one of the curved shape in Fig. 11. The ideal path of a proton through this sector has a radius of curvature of 1.9 m. This path is a circular arc that subtends an angle of 601 at the centre of curvature. The F sector on either end takes up 151 of the arc. The 301 at the centre makes up the D sector. The field in the F sector contains a dipole component and a quadrupole component. The dipole component is 1.28 T. This is the strength of the constant field that a 250 MeV proton following the ideal path would see. The quadrupole component has an n value of 5.855. The n value is defined by n¼
r @B B @r
ð22Þ
4
F1
D F2
30°
sextupole 2 2
sector
2m 1.9 m ρ
15°
θ
15°
proton
centres of curvature Y (m)
the circle through that point. After tracking through a few magnets, the proton missed the next magnet and literally went off at a tangent. The problem is that the tangent of the circle is not close to the direction for the closed orbit. Nor does the closed orbit cross near the middle of the marker plane. Fortunately, with some trial and error that involves shifting the start position along the marker plane and steering the initial direction of the proton, we manage to persuade the proton to move through the next magnet, followed by the next, and so on. Once it is able to do 10–20 turns, the proton continues for 100 turns without getting lost. In order to find the closed orbit, we track the proton for one turn round the ring and record the position and velocity when it intersects the marker plane again. We repeat this for 100 turns. The radial distances and radial velocities are then plotted. The result is one of the loops labelled 60 MeV in the phase space (Fig. 10). This loop is also called a trajectory. We find the average position (X, Y) and the average of the velocity component vX of these points. The component vY can be computed given the energy. Then we use these as the new starting conditions, and repeat the above steps. After each iteration, the trajectory in the phase space would shrink to a much smaller loop. When the size of loop falls below an acceptable level of tolerance, we consider that the result has converged. We can then take one turn of the trajectory as the closed orbit. Using the above method, we find the closed orbits at the energies of 15 MeV, 60 MeV and 180 MeV. The results are shown in Fig. 8. For each of the energies, a few phase space trajectories on the marker plane are shown in Fig. 10. To obtain these, we go the point where the closed orbit at, say, 60 MeV intersects the marker plane. We then move a few mm away horizontally along the plane and use the new position as the start point for the proton. Tracking this for 100 turns gives one of the phase space trajectories in Fig. 10. To get another trajectory, we go back to the marker plane, move a few mm further, and track for 50 turns again to get another phase space trajectory. We repeat this until the oscillations are so large that the particle gets lost. After that, we select a few and plotted them in Fig. 10. Then we repeat it for 15 MeV and 180 MeV. The largest for each energy is close to the stability limit, beyond which the proton gets lost—typically, it leaves a magnet at such a large angle that it does not enter the next magnet. It is interesting to compare Fig. 10 to the corresponding results in Ref. [11]. There are similarities in shapes, sizes and positions. The quantitative differences may be attributed to our simplification of the magnet field, and our choice of the marker plane position which is not given in Ref. [11].
145
entrance plane
0
exit plane
sextupole 1
-2 marker -4
extraction
0.2 60 MeV
180 MeV
0.15
-4
0.1
r’ (rad)
0 X (m)
2
Fig. 11. Lattice and closed orbit in compact synchrotron.
15 MeV
0.05
-2
0 Table 2 Parameters for the compact synchrotron.
-0.05 -0.1 -0.15 -0.2 -0.25 2.6
2.7
2.8
2.9
3
3.1 3.2 r (m)
3.3
3.4
Fig. 10. Phase space trajectories in the spiral FFAG.
3.5
3.6
Number of sectors Bending angle of magnet Radius of curvature Drift distance Angle of D sector Angle of F sector n value of D sector n value of F sector Energy Dipole field component
6 601 1.9 m 2m 301 151 6.164 5.855 250 MeV 1.28 T
4
146
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
where r is the distance from the centre of curvature of the ideal path, and B is the magnetic field. The right hand expression is to be evaluated at r ¼ 1:9 m, the radius of curvature of the ideal path. The quadrupole gradient @B=@r can then be obtained. Denote this by kF. From the measured data in Refs. [14,15], we make a simple model of the magnetic field. We only consider motion in the mid plane here, and assume that the field in this plane is vertical. Let the dipole component be B0. The quadrupole component is given by B1 ¼ k F x
ð23Þ
where x is the transverse, horizontal displacement from the ideal path. We model the fringe field with straight edges again. If a proton moves along the ideal path through the F sector, it sees a linear increase in the field from zero, followed by a constant field in the magnet, and then a linear decrease as it leaves. Denote this angular dependance by F 1 ðyÞ. The overall field in the F sector is then Bz ¼ ðB0 þ kF xÞF 1 ðyÞ:
ð24Þ
Denote the angular dependances for the other F and the D sectors by F 2 ðyÞ and DðyÞ respectively, and the gradient in the D sector by kD. The field for the complete F–D–F triplet is then Bz ðr, yÞ ¼ ðB0 þ kF xÞF 1 ðyÞ þðB0 þ kD xÞDðyÞ þðB0 þkF xÞF 2 ðyÞ:
ð25Þ
The angular functions are shown in Fig. 12. Notice how the fringe fields are approximated by the straight lines. We simulate the trajectory for the extraction of the protons. This synchrotron is designed for slow extraction [13]. Slow extraction is normally used in proton therapy because a small continuous current is needed. It makes use of sextupoles. These reduce the dynamic aperture of the beam. Even for a stable beam, there would be some betatron oscillation of a proton about the closed orbit. If only quadrupoles are present, the beam would become unstable when the amplitude is larger than some value (at which point the amplitude grows with time). The amplitude at which this starts to happen is the dynamic aperture. By increasing the sextupole strengths, the range of stable amplitudes would decrease, that is the dynamic aperture shrinks. If it shrinks below the amplitude of the oscillating proton, the proton becomes unstable and its amplitude starts to grow. It is possible to arrange the positions of the sextupoles such that when the proton reaches the place where it needs to be extracted, it has a large displacement towards the septum magnet. The proton would then come under the influence of the septum field and be guided into the extraction beamline [16].
Slow extraction can be done in other ways. One of these is to keep the sextupole strengths fixed, and use a magnet at some point along the ring to steer the beam. Steering the beam increases the transverse velocity, which in turn increases the amplitude of the betatron oscillation. When this amplitude exceeds the dynamics aperture imposed by the sextupoles, the proton again becomes unstable and can be extracted in the same way as before. We can observe this by simulation. The sextupoles used for extraction are at the positions labelled in Fig. 11, and the extraction location [15] is labelled ‘marker’. We first find the closed orbit when the sextupole strengths are zero. We launch a 250 MeV proton from a point on the marker plane, with velocity directed along þX. This is the obvious start condition to try since it is along in reference trajectory, which is known. The closed orbit is slightly different from this reference path because of the presence of the fringe fields. Even so, with the above start conditions, the proton remains close to the reference path after we track it for many turns. This means that the closed orbit, shown in Fig. 11, also remains near the reference path. This stability is not surprising, since there are no nonlinear field yet. Things change when we turn on the sextupoles. We now set the strength of one sextupole so that @2 Bz ¼ 600 T=m2 : @x2
ð26Þ
The other sextupole has the same strength, but with the field in the opposite direction. So in the mid plane of the magnets, the ‘Sextupole 1’ in Fig. 11 has its field downwards: Bz ¼ 300x2 T
ð27Þ
(x is in metres) and ‘Sextupole 2’ has its field upwards: Bz ¼ þ 300x2 T:
ð28Þ
In the simulation, we steer the beam by changing the initial direction of the proton at the marker plane. In practice, this would have to be done elsewhere, since the marker is at the extraction location. The proton is then tracked around the ring. The position and velocity at the plane after each turn are recorded. The Y components of the positions and velocities are then plotted in Fig. 13 for a range of steering angles. Note the characteristic triangular shapes [16]. The long tail extending to the left means large displacements in the Y direction. Looking at the marker position in Fig. 11, such a displacement makes it possible to extract the proton if a septum magnet is installed just below the
0.03 2 F1 D F2
1.5
0.02
Y’ (rad)
Angular Function
0.01 1 0.5
0 -0.01
0
-0.02 -0.5
-0.03 -1
-3.7 -10
0
10
20
30
40
50
60
70
-3.68 -3.66 -3.64 -3.62
-3.6
-3.58 -3.56
Y (m)
Angle (deg) Fig. 12. Angular functions for sectors in the compact synchrotron.
Fig. 13. Phase space trajectories near the septum of the compact synchrotron during low extraction.
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
marker in the figure. This is also where the extraction beamline is designed to be in Ref. [14].
4. Conclusion In the small FFAG community in the UK [17–19], Zgoubi [7] is often used in tracking simulations—for the EMMA accelerator at Daresbury and the PAMELA accelerator being designed at Oxford, just to name two examples. As the code is sophisticated and the learning curve steep, one of the authors (K.M.H.) hoped to have a simple tracking method so that we can easily write the code ourselves. The idea of using tiny arcs seems simple enough. The codes used for the simulations in this paper have indeed been straightforward to write, and the performance compared well with Zgoubi for the simple test case we have used. When we started writing this paper, a wider literature search revealed that the same method has been developed by PSI in 1987, in a code called TRACK [6]. We are not able to find in the literature mathematical details of the algorithm used, and hope that our contribution in this paper fills this small gap. TRACK has been mainly applied to cyclotrons [5] that are used for cancer therapy. In cyclotrons, the magnetic field fills much of the space covered by the proton trajectory. In FFAGs and synchrotrons, the proton tends to spend more of its time drifting in spaces with zero field. The fields are concentrated in lumped elements of magnets. It would reduce computational time significantly if we simply draw a straight line in between magnets, and only compute the tiny arcs inside each magnet. To provide this feature in our method, we define entrance and exit planes for each magnet. We use the Newton–Raphson method to accurately determine the position and velocity where a trajectory crosses the exit plane. Since the trajectory moves in a straight line to the next magnet, the position and velocity where it arrives at the next entrance plane can be computed quickly. This procedure works well when we apply it to the spiral FFAG and the compact synchrotron. To benchmark our code, we tested it on a FODO lattice. We compare this with a simulation using Zgoubi on the same lattice. In each case, the steps sizes are reduced until the results converged to specified tolerances. The trajectories computed using both methods agree well. For this lattice, the arc method is faster when the required accuracy is lower, whereas Zgoubi is faster when higher accuracy is needed. Whether the arc method is an acceptable solution would therefore depend on the type of lattice and tracking studies required. For our simulations on the spiral FFAG, we have used a very simple model to approximate the field. We have computed the closed orbits and the phase space trajectories for a few energies. The results are quite close to those in Ref. [11] which uses the detailed field maps obtained from measurements. For the compact synchrotron, we computed the phase space trajectories at the extraction location, when the extraction sextupoles are on and the beam is steered. The results are the triangular shaped distributions that are expected in the phase space when slow extraction is executed [16]. As far as we know, this phase space result has not been previously published by the group [14] that designed this synchrotron.
147
The above results demonstrate that the arc method works for accelerators other than cyclotrons, and that it can be applied in a flexible way to simulate different beam behaviours. In each case, the length of the arcs have to be reduced until the trajectories converged. The subsequent calculations for each phase space trajectory in the above examples usually takes less than a few minutes. This computing time could be reduced if we use longer step size where the field varies more gently, a feature that is employed in TRACK [6]. This is one area for future work. Another is to develop the method for tracking in electric fields, particularly the time varying ones in the rf cavities.
Acknowledgment We would like to thank Andy Wolski for helpful discussions on symplecticity.
References [1] D.R. Douglas, Particle Accelerators 19 (1986) 119. [2] J.C. Smith, et al., Comparison of tracking codes for the international linear collider, in: Proceedings of the Particle Accelerator Conference, Albuquerque, USA, 2007, pp. 3020–3022. [3] V.N. Aseev, et al., TRACK: the new beam dynamics code, in: Proceedings of the Particle Accelerator Conference, Knoxville, USA, 2005, pp. 2053–2055. [4] F. Me´ot, Nuclear Instruments and Methods in Physics Research Section A 427 (1999) 353. [5] J.M. Schippers, D.C. George, V. Vrankovic, Results of 3D beam dynamic studies in distorted fields of a 250 MeV superconducting cyclotron, in: Cyclotrons 2004, Tokyo, 18–22 October 2004. [6] V. Vrankovı´c, D.C. George, Particle ray-tracing program TRACK applied to accelerators, in: PSI Scientific and Technical Report 2001, vol. 6: Large Research Facilities, pp. 27–28. [7] F. Me´ot, F. Lemuet, Zgoubi’s Users’ Guide. CEA DAPNIA SEA-97-13 /http:// sourceforge.net/projects/zgoubi/S, 1997. [8] J.R. Rees, Symplecticity in beam dynamics: an introduction, in: SLAC-PUB9939, 2003, pp. 44–46 /http://www.slac.stanford.edu/pubs/slacpubs/9000/ slac-pub-9939.htmlS. [9] A.J. Dragt, Nuclear Instruments and Methods in Physics Research Section A 258 (1987) 339. [10] J. Fourrier, et al., Nuclear Instruments and Methods in Physics Research Section A 589 (2008) 133. [11] S. Antoine, Nuclear Instruments and Methods in Physics Research Section A 602 (2009) 293. [12] M.-J. Leray, et al., Magnetic measurements of the RACCAM prototype FFAG dipole, in: Proceedings of the Particle Accelerator Conference, Vancouver, BC, Canada, 2009, pp. 5032–5034. [13] A. Noda, et al., Design study of compact proton synchrotron for biomedical use, in: Proceedings of the European Particle Accelerator Conference, Berlin, Germany, 1992, pp. 1663–1665. [14] A. Morita, et al., A compact synchrotron with combined-function lattice dedicated for cancer therapy, in: Proceedings of the Particle Accelerator Conference, New York, USA, 1999, pp. 2528–2530. [15] A. Morita, et al., Field measurement of combined function magnet for medical proton synchrotron, in: Proceedings of the European Particle Accelerator Conference, Vienna, Austria, 2000, pp. 2110–2112. [16] G.H. Rees, Extraction, in: Proceedings of the CERN Accelerator School, CERN 85-19, Paris, France, 1984, pp. 346–357. [17] Y. Giboudot, et al., Particle tracking studies using dynamical map created from finite element solution of the EMMA cell, in: Proceedings of the Particle Accelerator Conference, Vancouver, BC, Canada, 2009, pp. 3290–3292. [18] R. Fenning, et al., An FFAG transport line for the PAMELA project, in: Proceedings of the Particle Accelerator Conference, Vancouver, BC, Canada, 2009, pp. 3741–3743. [19] S.L. Sheehy, The EMMA accelerator and implications for hadron non-scaling FFAGs, in: Proceedings of HB2010, Morschach, Switzerland, 2010, pp. 82–85.
Contents lists available at SciVerse ScienceDirect
Nuclear Instruments and Methods in Physics Research A journal homepage: www.elsevier.com/locate/nima
Tracking protons in an FFAG and a synchrotron with tiny arcs$ K.M. Hock a,b,n, C.S. Edmonds a,b a b
Department of Physics, University of Liverpool, Liverpool L69 7ZE, United Kingdom Cockcroft Institute, Daresbury Science and Innovation Campus, Daresbury, Warrington WA4 4AD, United Kingdom
a r t i c l e i n f o
a b s t r a c t
Article history: Received 1 August 2011 Accepted 26 September 2011 Available online 9 November 2011
We study a ray tracing method that uses tiny arcs to track a proton through a magnetic field. We show how to apply this to an FFAG and a synchrotron that are designed for cancer therapy. Assuming that the field traversed by a proton is uniform over each time step, the required equations and a simple scheme for treating a lattice are developed. We compare the computation time with the code Zgoubi for a FODO lattice. To test its usefulness, it is applied to the spiral FFAG, designed in the RACCAM project, to find the horizontal phase space trajectories. It is also applied to a compact synchrotron, designed in Japan, to determine the phase space trajectories during slow extraction of the proton beam. & 2011 Elsevier B.V. All rights reserved.
Keywords: FFAG Synchrotron Proton therapy Tracking Ray tracing
1. Introduction This work is motivated by the need for a simple and efficient way to track protons in a cancer therapy machine. Tracking of charged particles is an important part of the design of an accelerator. It involves calculating the path of the particles to see if they are still stable after many turns. There are many methods for doing so. These can be broadly classified into two maingroups—transfer maps and ray tracing. A transfer map is a way to quickly compute the position and velocity at some point along a beamline, given the co-ordinates at a starting point. This is done without the need to compute the positions of the trajectory in between. This map could make use of a matrix, a power series, or Lie algebra [1]. It would normally involve some approximations, such as assuming that the map is linear, or truncating the power series after a few orders. These make the calculations very fast, but usually mean that the results are accurate only for small displacements from the reference orbit. Many codes have been written using these methods [2]. They are particularly useful for high energy accelerators, where the particles are expected to stay close to the reference orbit. At low energies, or in more complex elements of the accelerator, trajectories could either have large curvatures or move sideways by large distances. For accurate tracking, ray tracing has to be used. This involves integrating the equation of motion of the particle in an electromagnetic field in some way. This means that
$
Work supported by the Science and Technology Facilities Council, UK. Corresponding author at: Department of Physics, University of Liverpool, Liverpool L69 7ZE, United Kingdom. E-mail address: [email protected] (K.M. Hock). n
0168-9002/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2011.09.066
unlike transfer maps, the co-ordinates between any two points on a trajectory must be computed. The methods used include the Runge–Kutta method [3], power series expansion of the equations of motion [4], and tiny arcs [5]. The time or distance steps can be reduced until sufficient accuracy is achieved. For long trajectories, very small steps are required, which could be time consuming. In this paper, we develop on the method that uses tiny arcs. The field traversed by a proton in each time step is assumed to be uniform. The trajectory in this time step is therefore an arc. This method has been developed by another group for tracking in cyclotrons [6]. We modify this to make it more efficient for accelerators that are made up of lumped elements with large drift spaces in between. In Section 2, we develop the algorithm and test the code on a FODO lattice. In Section 3, we apply this to an FFAG and a synchrotron. In Section 4, we summarise the results and suggest further work.
2. The arc method A charged particle travelling through a magnetic field experiences a force given by F ¼ qv B
ð1Þ
where q is the charge of the particle, v is the velocity and B the magnetic field. To determine the path through a nonuniform field, numerical methods are needed. The code Zgoubi [7], for example, uses Taylor expansions up to a few orders to approximate the path in each time step. In these notes, I consider using arcs instead. As I intend to develop this method for the study of the proton accelerator, I shall take proton as the charge particle.
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
Consider a proton travelling through a nonuniform field. In a very small time step, the field over which the proton would move would be approximately constant. If the velocity is perpendicular to the field, the path would be a perfectly circular arc. If the velocity has a component parallel to the field, the arc also gains a component in the direction of the field. The arc becomes a spiral arc. Assuming that the velocity is approximately constant in this short time interval, we would still be able to calculate the endpoint of this spiral arc. The idea of the arc method is to start with the new field at the endpoint of this arc, and find the next arc. This is repeated, and the whole trajectory is obtained by joining many arcs together. After that, it is important to repeat the calculations for a number of smaller time steps to check that the answer has converged. From the starting position x0 , we need a formula to calculate the endpoint x1 of the arc after a time interval dt, in terms of B and v0 , the velocity at the starting position. We first consider only the perpendicular component v? and find the circular arc s, as shown in Fig. 1. Then we add the displacement due to the parallel component vJ . To do these, we first obtain the following. The constant magnetic field should ideally be taken as the field at the midpoint of the arc. However, we do not know where this midpoint is before finding the arc. A simple approximation would be to go in a straight line from the starting position for half a time step to the point x ¼ x0 þv0 dt=2
ð2Þ
and take the magnetic field here as the constant field B. The unit vector in the direction of the field is B B^ ¼ : B
ð3Þ
The parallel component of the velocity is ^ B: ^ vJ ¼ ðv0 :BÞ
ð4Þ
The perpendicular component is then v? ¼ v0 vJ :
ð5Þ
The radius of the circular arc is given by r¼
gmv?
ð6Þ
qB
where m is the rest mass and g is the Lorentz factor. We first consider the component of the motion that is perpendicular to the field. The circular arc starts off from x0 , goes in the direction v? , and curves towards F —since F points to the centre of the arc. After time dt, the angle y that the arc subtends at the centre is
y¼
v? d t ¼ odt r
ð7Þ
where we have also defined o. To express the endpoint of the arc in terms of these, we use the following unit vectors. The unit
141
vector in the direction of v? is given by v^ ? ¼
v? v?
ð8Þ
and the unit vector in the direction of F is given by F F^ ¼ : F
ð9Þ
Therefore, the endpoint of the circular arc is x1? ¼ x0 þðrr cos yÞF^ þ r sin yv^ ? :
ð10Þ
Finally, we add the displacement parallel to the field to get the endpoint of the spiral arc, x1 ¼ x0 þ ðrr cos yÞF^ þ ðr sin yÞv^ ? þvJ dt:
ð11Þ
Differentiating this with respect to time, we obtain the velocity at the endpoint: v1 ¼ ðr o sin yÞF^ þ ðr o cos yÞv^ ? þvJ :
ð12Þ
We shall call the last two equations the arc formula. To use the arc formula to develop a computer code for tracking a proton in an accelerator with lumped elements, we assume the following: 1. The beamline is made up of elements of a finite length, containing magnetic fields. 2. A proton would usually enter from one side of an element, and exit from the opposite side. 3. In between elements are drift spaces, with no field. Elements with electric fields are not considered in this paper. With the above assumptions, we can always define two imaginary planes for each element, one on each side. For example, the planes could be the entrance and exit faces of a quadrupole or dipole. These planes could be roughly perpendicular to the beamline. A proton would enter one plane, go through the element, and come out through the other plane. In between the two planes, there would be a magnetic field. The field can be of any distribution. The two planes for each element can be used to define the boundaries of the element in the code. The arc formula can then be used to compute the trajectory of the proton within these boundaries. There is one problem. The arc formula used a fixed time step dt. When the trajectory formed by the arcs crosses the plane, the field could drop abruptly to zero just outside the plane. This happens if we use the hard edge models for a dipole or a quadrupole. Since the final arc in the element would in general not terminate exactly on the exit plane, the arc is partly inside the field and partly outside. This is different from the original assumption that the field is approximately constant over the extent of the arc, and would lead to inaccuracy. To avoid this problem, we can make the final arc terminate exactly on the exit plane by applying the Newton– Raphson method. The formula is given by t1 ¼ t0
f ðt 0 Þ 0 f ðt 0 Þ
ð13Þ
where f(t) is a function that crosses zero for some t. By applying the formula iteratively—computing t1, replacing t0 by this, computing t1 again—t1 would get closer and closer to the root. We can apply this by taking f(t) as the distance from a point on the final arc to the exit plane, and renaming dt as t. To find f(t), rewrite the arc formula in this form: x1 ¼ x0 þ ðrr cos otÞF^ þ ðr sin otÞv^ ? þ vJ t: Fig. 1. Trajectory in a uniform field.
ð14Þ
Next, define the exit plane by a point on the plane p2 , and the unit vector normal to the plane n2 . The perpendicular displacement of
142
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
a point on the arc from the plane is then f ðtÞ ¼ ðx1 p2 Þ:n2 :
ð15Þ
Differentiating this gives the derivative: 0 f ðtÞ ¼ ½ðv? sin otÞF^ þ ðv? cos otÞv^ ? þ vJ :n2 :
ð16Þ
Since dt tends to be small, the start point of the final arc is already quite close to the exit plane. So a few iterations is often enough to achieve good accuracy. In our simulations, we stop when the change in the next iteration is within 10 5 of dt. The value of t obtained is then used in place of dt in the arc formula for the final arc before the exit plane. A C code has been written based on the above descriptions. In the first test, we compute the path of a proton through a 1 T uniform field. The result is shown in Fig. 2. This is a trivial case, since the arc formula is exact for the uniform field, regardless of the dt. For this figure, I choose 1000 time steps, just to make the curve look smooth. It is only to check that the code is working. The spiral path is produced as expected. In the next test, we use a FODO lattice with 100 cells. The same 10 MeV proton is used. Each quadrupole is 10 cm long, and has a field gradient of 10 T/m. The distance between adjacent quadrupoles is 0.2 m. The trajectory computed using the arc formula is shown in Fig. 3. We check for convergence by looking at the y coordinate of the endpoint of the final arc, which falls on the exit face of the final quadrupole in the lattice. We start with a step size
of a few centimetres, carry out the tracking, record the y value and the computational time. The step size is then halved, and the tracking repeated. Each time the step size is halved, the difference between the current and previous y values is calculated. Assuming that the calculation converges, this difference would decrease as the step size decreases. We define the tolerance to be the fraction of this difference over the current y value. When the tolerance falls below a value that is small enough for our purpose, we can say that the y value, and therefore the tracking, has converged. Fig. 4 shows the tolerance, plotted against step size, for the arc method and for Zgoubi. Notice that to achieve the same tolerance, the step size required for the arc method is smaller than the step size required for Zgoubi. Fig. 5 shows the computation time plotted against tolerance. Our computer has a Pentium 3 GHz processor. The trends for the two methods are quite clear. The times would be different on a different computer, but the relative time between the two methods should remain the same, as this depends mainly on the amount of computations needed. For higher tolerance, the computation time for the arc method is shorter. There is a tolerance below which the time for Zgoubi becomes shorter.
101 100 10-1
Zgoubi Arc
x (m)
Tolerance
10-2 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4
10-3 10-4 10-5 10-6 10-7
0
0
2
4
z (m)
6
10-8
0.1
-0.1 -0.2 -0.3 -0.4 m) -0.5 y( -0.6 10
8
10-9 10-3
10-2
10-1
100
101
Step size (cm) Fig. 4. The fractional change in the endpoint y coordinate each time the step size is halved and is plotted against step size.
Fig. 2. Trajectory in a uniform field.
0.02
100 Zgoubi
0.015
Arc 10-1 Compute time (s)
0.01
y (m)
0.005 0 -0.005 -0.01
10-2
10-3
-0.015 -0.02 0
10
20
30
40
50
60
z (m) Fig. 3. Trajectory through 100 FODO cells. It crosses the quadrupoles at the thickened parts of the curve.
10-4 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 101 Tolerance Fig. 5. The computation time for tracking through 100 FODO cell is plotted against the tolerance required.
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
Tolerance is the fractional change in the result from one step size to the next. It is not the same as accuracy, which should be a measure of the error of the result compared to the correct answer. Fig. 6 shows the y coordinate of the endpoint of the trajectory plotted against the tolerance for both methods. The results of both methods appear to approach the same limiting value. Since we do not know the correct answer, we assume that this is the average of the values obtained using the smallest step sizes in the two methods. We then define the error for each step size, as the difference between the y value at the endpoint calculated by each method, and this nominal answer. Fig. 7 shows a graph of the computation time against the absolute error of each method. Note that the graphs would be less reliable for very small errors. Since our nominal answer is only a guess, the graphs are likely to be less accurate around the smallest step sizes. The graphs for the 100 FODO cells show that the arc method is faster if the required error is above 1 mm. Below that, Zgoubi is faster. Further study would be needed to determine the usefulness of the arc method in accelerator design, particularly when tracking for many turns round a ring is required. An issue that is often a concern in tracking methods is symplecticity [8]. Symplecticity is a condition on the transfer map that follows from a Hamiltonian system. Making a tracking
Zgoubi
-1.2986
143
method symplectic removes one source of error. One example is to approximate each tracking step using a drift, a kick, and a drift [8], so that each step is symplectic. Another example is to use Lie algebra [9] to create symplectic maps. A symplectic tracking method is still an approximation, of course. There could be other sources of errors that may not be small and have to be carefully checked. In the arc method, we see that within each time step, the tracking is symplectic. This is because the arc is the exact solution, once we assume that the magnetic field is constant over this step. So the joining of tiny arcs is in essence a concatenation of symplectic maps, which must also be symplectic. Although this is interesting to know, in practice, the arc is still an approximation to the real path and the arc length has to be made small enough to give good accuracy. Another issue that could be of interest is time reversibility. After tracking through a lattice, if we reverse the particle velocity, reverse all the magnetic fields, and track using the same method, would the particle arrive back at the same starting point? The answer is no. The reason is that at each time step, the field at the beginning of the step would in general be different from the field at the end of the step. If time is reversed, the end would become the start. Then over the same step, the field could be slightly different. This difference could be made small by reducing the time step. For fairly symmetrical magnetic fields, we may expect that the error in the time reversed tracking is comparable to the error in the forward tracking.
Arc Endpoint y (cm)
3. Results and discussion
-1.2988 3.1. Spiral FFAG
-1.299
We start with the tracking of protons in the spiral FFAG. This spiral FFAG is designed in the RACCAM project [10–12] for proton therapy. The lattice we simulated is shown in Fig. 8. It follows the parameters in Ref. [11]. Those relevant to our simulation are summarised in Table 1. The number of cells refer to the number of spiral magnets, including the associated drift spaces. The cells are identical.
-1.2992
-1.2994 10-7
10-6
10-5 10-4 Tolerance
10-3
4
10-2
3
Fig. 6. The y coordinate at the end of the tracking is plotted against the tolerance.
magnet
180 MeV 60 MeV
proton
15 MeV 2
100
(X,Y) Zgoubi
1
r
marker
Compute time (s)
10
Y (m)
Arc -1
0
β r1
-1
10-2
θ
r2
A
-2 10-3
-3 entrance plane exit plane
fringe edges 10-4 10-7
-4 10-6
10-5
10-4 Error (cm)
10-3
10-2
10-1
Fig. 7. The computation time is plotted against the absolute error for each method.
-4
-3
-2
-1
0 X (m)
1
2
Fig. 8. Lattice and closed orbits in the spiral FFAG.
3
4
144
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
2
Table 1 Parameters for the spiral FFAG [11]. N p k
z r1 r2 B2
10 0.34 5 53.71 2.7 m 3.5 m 1.8 T
Imagine drawing a circle that is centred at the centre of the accelerator, and that passes through the magnets. The packing factor is the fraction of the circle that is covered by the magnets. This fraction is the same for different radii, as long as the circle passes through the magnets. The field index k is the variation of the magnetic field with the radius r. This variation is given by
1.5
Bz (T)
Number of cells Packing factor Field index Spiral angle Minimum radius Maximum radius Maximum field
2.7 m 3.1 m 3.5 m
1
0.5
0 -30
Bz ¼ br
k
where Bz is the vertical component in the mid-plane in the magnet region. Ideally, the horizontal component at the midplane is zero. This is also what we assume in the simulation. The constant of proportion b may be obtained using the parameters in Table 1 that Bz is 1.8 T when r is 3.5 m. Using the equation and the parameters in this way gives us a simplified model of the measured field data in Ref. [12]. A sophisticated model is given in Ref. [11]. The maximum and minimum radii in Table 1 indicate the region where we expect the protons to travel. The actual size of the magnet is larger. The spiral angle is the angle between a radius and a point on a spiral edge of a magnet. Each magnet has four curved edges. Two of these are circular arcs. These are the one nearest to the centre of the accelerator, and the one furthest away. The other two are spiral arcs. The radius refer to a straight line from the centre of the accelerator. Imagine drawing a radius to a point on any one of the spiral arcs. Then this straight line would always make an angle of 53.71 with the tangent to the spiral arc, at the point where they meet. This is used as a condition to define the exact shape of the spiral arc. From this condition, the equation of the spiral arc can be derived as r ¼ r 0 exp
y
tan z
ð18Þ
where z is the 53.71 spiral angle. The parameter r0 would be different for each of the spiral edges. A prototype of the magnet has been constructed and the field measured [12]. A significant part of the field is fringe field. We use a simple model for this field. Imagine again a circle through the magnets. If we move into the field of a magnet along this circle, the field increases, stay the same, then decreases again. The increase and decrease are the fringe fields. In our simulation, we simplify the field using a linear variation for the fringe field along the circle. To describe this variation, we may define an angular function FðyÞ. The overall field for a magnet may then be given by k
-25
ð17Þ
Bz ¼ br Fðyy0 Þ:
ð19Þ
The resulting field is shown in Fig. 9 for the three radial distances 2.7 m, 3.1 m and 3.5 m. These correspond roughly to the positions of the closed orbits at 15 MeV, 60 MeV and 180 MeV respectively—see Fig. 8. To understand what y0 is, consider the magnet just below the ‘marker’ label in Fig. 8. Imagine a spiral curve that is mid way between the two spiral edges, and that is described by the same Eq. (18) as the edges. The point on this curve closest to the centre of the ring is at distance r1 from the centre. Suppose y is zero
-20
-15
-10 θ (deg)
-5
0
5
10
Fig. 9. Field model for the spiral FFAG magnet.
there. The curve would be given by y r ¼ r 1 exp : tan z
ð20Þ
The y in this equation is then equal to y0 for each value of r. In this way, the angular function Fðyy0 Þ is made to follow the spiral. Using the same angular function for different r is possible because the length of the circular arc through the magnet is always equal to A. This is true for any r, since one spiral edge can be obtained from the other by a rigid rotation about the centre of the ring. However, the actual extent of the fringe fields change at different r, so using the same angular function is an approximation. The field for each of the N magnets is then given by k
Bz ¼ br Fðyy0 nbÞ
ð21Þ
where n ¼ 0; 1, . . . ,N1. This completes our simple model of the magnetic field in the whole lattice. We should emphasise that we make this model just to demonstrate the arc method. The accurate model in Ref. [11] can be used instead, as can any field map. Next, we need to specify the entrance and exit planes for each magnet. These are shown by the two straight lines next to each magnet in Fig. 8. These planes have to be just beyond the fringe fields, so that the full extent of the field would be sampled during tracking of the proton. The planes should not be too far from the fringe fields, so as to give more drift distance between the exit plane of one magnet and the entrance plane of the next. This is because tracking over this distance requires little computation. The most direct way to determine the positions of the entrance and exit planes is by accurately plotting both their positions and the boundaries of the fringe fields on the same drawing, and checking them by eye. For our simulation, we plotted these as lines and curves on a graph. The curves labelled ‘fringe edges’ in Fig. 8 are the boundaries of the fringe fields. We typed in the positions of the planes by hand, replotted the graphs, and tried this a few times until we got a reasonable position for each plane. If the positions of the planes are correct, we would be able to track a proton though the ring and get a trajectory that looks like one of the orbits shown in Fig. 8. In practice, we find it useful to make the code track the proton through one magnet at a time, each time plotting out the trajectory with the magnet outlines, as in Fig. 8, to check. We start by finding the closed orbit at 60 MeV, using a method described in Ref. [7]. We start near the middle of the plane labelled ‘marker’ in Fig. 8, with a velocity in the direction of the tangent to
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
3.2. Compact synchrotron A compact synchrotron is designed in Japan for cancer therapy [13–15]. The lattice is outlined in Fig. 11. It makes use of six combined function magnets. The fields in these magnets contain dipole and quadrupole components. Each magnet is an F–D–F triplet. Since these are linear, tracking with transfer maps are usually sufficient for trajectories close to the reference orbit. We consider here the extraction region of the orbit, where large displacements from the reference are expected. The relevant parameters for the synchrotron lattice are listed in Table 2. Each F–D–F triplet forms one of the curved shape in Fig. 11. The ideal path of a proton through this sector has a radius of curvature of 1.9 m. This path is a circular arc that subtends an angle of 601 at the centre of curvature. The F sector on either end takes up 151 of the arc. The 301 at the centre makes up the D sector. The field in the F sector contains a dipole component and a quadrupole component. The dipole component is 1.28 T. This is the strength of the constant field that a 250 MeV proton following the ideal path would see. The quadrupole component has an n value of 5.855. The n value is defined by n¼
r @B B @r
ð22Þ
4
F1
D F2
30°
sextupole 2 2
sector
2m 1.9 m ρ
15°
θ
15°
proton
centres of curvature Y (m)
the circle through that point. After tracking through a few magnets, the proton missed the next magnet and literally went off at a tangent. The problem is that the tangent of the circle is not close to the direction for the closed orbit. Nor does the closed orbit cross near the middle of the marker plane. Fortunately, with some trial and error that involves shifting the start position along the marker plane and steering the initial direction of the proton, we manage to persuade the proton to move through the next magnet, followed by the next, and so on. Once it is able to do 10–20 turns, the proton continues for 100 turns without getting lost. In order to find the closed orbit, we track the proton for one turn round the ring and record the position and velocity when it intersects the marker plane again. We repeat this for 100 turns. The radial distances and radial velocities are then plotted. The result is one of the loops labelled 60 MeV in the phase space (Fig. 10). This loop is also called a trajectory. We find the average position (X, Y) and the average of the velocity component vX of these points. The component vY can be computed given the energy. Then we use these as the new starting conditions, and repeat the above steps. After each iteration, the trajectory in the phase space would shrink to a much smaller loop. When the size of loop falls below an acceptable level of tolerance, we consider that the result has converged. We can then take one turn of the trajectory as the closed orbit. Using the above method, we find the closed orbits at the energies of 15 MeV, 60 MeV and 180 MeV. The results are shown in Fig. 8. For each of the energies, a few phase space trajectories on the marker plane are shown in Fig. 10. To obtain these, we go the point where the closed orbit at, say, 60 MeV intersects the marker plane. We then move a few mm away horizontally along the plane and use the new position as the start point for the proton. Tracking this for 100 turns gives one of the phase space trajectories in Fig. 10. To get another trajectory, we go back to the marker plane, move a few mm further, and track for 50 turns again to get another phase space trajectory. We repeat this until the oscillations are so large that the particle gets lost. After that, we select a few and plotted them in Fig. 10. Then we repeat it for 15 MeV and 180 MeV. The largest for each energy is close to the stability limit, beyond which the proton gets lost—typically, it leaves a magnet at such a large angle that it does not enter the next magnet. It is interesting to compare Fig. 10 to the corresponding results in Ref. [11]. There are similarities in shapes, sizes and positions. The quantitative differences may be attributed to our simplification of the magnet field, and our choice of the marker plane position which is not given in Ref. [11].
145
entrance plane
0
exit plane
sextupole 1
-2 marker -4
extraction
0.2 60 MeV
180 MeV
0.15
-4
0.1
r’ (rad)
0 X (m)
2
Fig. 11. Lattice and closed orbit in compact synchrotron.
15 MeV
0.05
-2
0 Table 2 Parameters for the compact synchrotron.
-0.05 -0.1 -0.15 -0.2 -0.25 2.6
2.7
2.8
2.9
3
3.1 3.2 r (m)
3.3
3.4
Fig. 10. Phase space trajectories in the spiral FFAG.
3.5
3.6
Number of sectors Bending angle of magnet Radius of curvature Drift distance Angle of D sector Angle of F sector n value of D sector n value of F sector Energy Dipole field component
6 601 1.9 m 2m 301 151 6.164 5.855 250 MeV 1.28 T
4
146
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
where r is the distance from the centre of curvature of the ideal path, and B is the magnetic field. The right hand expression is to be evaluated at r ¼ 1:9 m, the radius of curvature of the ideal path. The quadrupole gradient @B=@r can then be obtained. Denote this by kF. From the measured data in Refs. [14,15], we make a simple model of the magnetic field. We only consider motion in the mid plane here, and assume that the field in this plane is vertical. Let the dipole component be B0. The quadrupole component is given by B1 ¼ k F x
ð23Þ
where x is the transverse, horizontal displacement from the ideal path. We model the fringe field with straight edges again. If a proton moves along the ideal path through the F sector, it sees a linear increase in the field from zero, followed by a constant field in the magnet, and then a linear decrease as it leaves. Denote this angular dependance by F 1 ðyÞ. The overall field in the F sector is then Bz ¼ ðB0 þ kF xÞF 1 ðyÞ:
ð24Þ
Denote the angular dependances for the other F and the D sectors by F 2 ðyÞ and DðyÞ respectively, and the gradient in the D sector by kD. The field for the complete F–D–F triplet is then Bz ðr, yÞ ¼ ðB0 þ kF xÞF 1 ðyÞ þðB0 þ kD xÞDðyÞ þðB0 þkF xÞF 2 ðyÞ:
ð25Þ
The angular functions are shown in Fig. 12. Notice how the fringe fields are approximated by the straight lines. We simulate the trajectory for the extraction of the protons. This synchrotron is designed for slow extraction [13]. Slow extraction is normally used in proton therapy because a small continuous current is needed. It makes use of sextupoles. These reduce the dynamic aperture of the beam. Even for a stable beam, there would be some betatron oscillation of a proton about the closed orbit. If only quadrupoles are present, the beam would become unstable when the amplitude is larger than some value (at which point the amplitude grows with time). The amplitude at which this starts to happen is the dynamic aperture. By increasing the sextupole strengths, the range of stable amplitudes would decrease, that is the dynamic aperture shrinks. If it shrinks below the amplitude of the oscillating proton, the proton becomes unstable and its amplitude starts to grow. It is possible to arrange the positions of the sextupoles such that when the proton reaches the place where it needs to be extracted, it has a large displacement towards the septum magnet. The proton would then come under the influence of the septum field and be guided into the extraction beamline [16].
Slow extraction can be done in other ways. One of these is to keep the sextupole strengths fixed, and use a magnet at some point along the ring to steer the beam. Steering the beam increases the transverse velocity, which in turn increases the amplitude of the betatron oscillation. When this amplitude exceeds the dynamics aperture imposed by the sextupoles, the proton again becomes unstable and can be extracted in the same way as before. We can observe this by simulation. The sextupoles used for extraction are at the positions labelled in Fig. 11, and the extraction location [15] is labelled ‘marker’. We first find the closed orbit when the sextupole strengths are zero. We launch a 250 MeV proton from a point on the marker plane, with velocity directed along þX. This is the obvious start condition to try since it is along in reference trajectory, which is known. The closed orbit is slightly different from this reference path because of the presence of the fringe fields. Even so, with the above start conditions, the proton remains close to the reference path after we track it for many turns. This means that the closed orbit, shown in Fig. 11, also remains near the reference path. This stability is not surprising, since there are no nonlinear field yet. Things change when we turn on the sextupoles. We now set the strength of one sextupole so that @2 Bz ¼ 600 T=m2 : @x2
ð26Þ
The other sextupole has the same strength, but with the field in the opposite direction. So in the mid plane of the magnets, the ‘Sextupole 1’ in Fig. 11 has its field downwards: Bz ¼ 300x2 T
ð27Þ
(x is in metres) and ‘Sextupole 2’ has its field upwards: Bz ¼ þ 300x2 T:
ð28Þ
In the simulation, we steer the beam by changing the initial direction of the proton at the marker plane. In practice, this would have to be done elsewhere, since the marker is at the extraction location. The proton is then tracked around the ring. The position and velocity at the plane after each turn are recorded. The Y components of the positions and velocities are then plotted in Fig. 13 for a range of steering angles. Note the characteristic triangular shapes [16]. The long tail extending to the left means large displacements in the Y direction. Looking at the marker position in Fig. 11, such a displacement makes it possible to extract the proton if a septum magnet is installed just below the
0.03 2 F1 D F2
1.5
0.02
Y’ (rad)
Angular Function
0.01 1 0.5
0 -0.01
0
-0.02 -0.5
-0.03 -1
-3.7 -10
0
10
20
30
40
50
60
70
-3.68 -3.66 -3.64 -3.62
-3.6
-3.58 -3.56
Y (m)
Angle (deg) Fig. 12. Angular functions for sectors in the compact synchrotron.
Fig. 13. Phase space trajectories near the septum of the compact synchrotron during low extraction.
K.M. Hock, C.S. Edmonds / Nuclear Instruments and Methods in Physics Research A 664 (2012) 140–147
marker in the figure. This is also where the extraction beamline is designed to be in Ref. [14].
4. Conclusion In the small FFAG community in the UK [17–19], Zgoubi [7] is often used in tracking simulations—for the EMMA accelerator at Daresbury and the PAMELA accelerator being designed at Oxford, just to name two examples. As the code is sophisticated and the learning curve steep, one of the authors (K.M.H.) hoped to have a simple tracking method so that we can easily write the code ourselves. The idea of using tiny arcs seems simple enough. The codes used for the simulations in this paper have indeed been straightforward to write, and the performance compared well with Zgoubi for the simple test case we have used. When we started writing this paper, a wider literature search revealed that the same method has been developed by PSI in 1987, in a code called TRACK [6]. We are not able to find in the literature mathematical details of the algorithm used, and hope that our contribution in this paper fills this small gap. TRACK has been mainly applied to cyclotrons [5] that are used for cancer therapy. In cyclotrons, the magnetic field fills much of the space covered by the proton trajectory. In FFAGs and synchrotrons, the proton tends to spend more of its time drifting in spaces with zero field. The fields are concentrated in lumped elements of magnets. It would reduce computational time significantly if we simply draw a straight line in between magnets, and only compute the tiny arcs inside each magnet. To provide this feature in our method, we define entrance and exit planes for each magnet. We use the Newton–Raphson method to accurately determine the position and velocity where a trajectory crosses the exit plane. Since the trajectory moves in a straight line to the next magnet, the position and velocity where it arrives at the next entrance plane can be computed quickly. This procedure works well when we apply it to the spiral FFAG and the compact synchrotron. To benchmark our code, we tested it on a FODO lattice. We compare this with a simulation using Zgoubi on the same lattice. In each case, the steps sizes are reduced until the results converged to specified tolerances. The trajectories computed using both methods agree well. For this lattice, the arc method is faster when the required accuracy is lower, whereas Zgoubi is faster when higher accuracy is needed. Whether the arc method is an acceptable solution would therefore depend on the type of lattice and tracking studies required. For our simulations on the spiral FFAG, we have used a very simple model to approximate the field. We have computed the closed orbits and the phase space trajectories for a few energies. The results are quite close to those in Ref. [11] which uses the detailed field maps obtained from measurements. For the compact synchrotron, we computed the phase space trajectories at the extraction location, when the extraction sextupoles are on and the beam is steered. The results are the triangular shaped distributions that are expected in the phase space when slow extraction is executed [16]. As far as we know, this phase space result has not been previously published by the group [14] that designed this synchrotron.
147
The above results demonstrate that the arc method works for accelerators other than cyclotrons, and that it can be applied in a flexible way to simulate different beam behaviours. In each case, the length of the arcs have to be reduced until the trajectories converged. The subsequent calculations for each phase space trajectory in the above examples usually takes less than a few minutes. This computing time could be reduced if we use longer step size where the field varies more gently, a feature that is employed in TRACK [6]. This is one area for future work. Another is to develop the method for tracking in electric fields, particularly the time varying ones in the rf cavities.
Acknowledgment We would like to thank Andy Wolski for helpful discussions on symplecticity.
References [1] D.R. Douglas, Particle Accelerators 19 (1986) 119. [2] J.C. Smith, et al., Comparison of tracking codes for the international linear collider, in: Proceedings of the Particle Accelerator Conference, Albuquerque, USA, 2007, pp. 3020–3022. [3] V.N. Aseev, et al., TRACK: the new beam dynamics code, in: Proceedings of the Particle Accelerator Conference, Knoxville, USA, 2005, pp. 2053–2055. [4] F. Me´ot, Nuclear Instruments and Methods in Physics Research Section A 427 (1999) 353. [5] J.M. Schippers, D.C. George, V. Vrankovic, Results of 3D beam dynamic studies in distorted fields of a 250 MeV superconducting cyclotron, in: Cyclotrons 2004, Tokyo, 18–22 October 2004. [6] V. Vrankovı´c, D.C. George, Particle ray-tracing program TRACK applied to accelerators, in: PSI Scientific and Technical Report 2001, vol. 6: Large Research Facilities, pp. 27–28. [7] F. Me´ot, F. Lemuet, Zgoubi’s Users’ Guide. CEA DAPNIA SEA-97-13 /http:// sourceforge.net/projects/zgoubi/S, 1997. [8] J.R. Rees, Symplecticity in beam dynamics: an introduction, in: SLAC-PUB9939, 2003, pp. 44–46 /http://www.slac.stanford.edu/pubs/slacpubs/9000/ slac-pub-9939.htmlS. [9] A.J. Dragt, Nuclear Instruments and Methods in Physics Research Section A 258 (1987) 339. [10] J. Fourrier, et al., Nuclear Instruments and Methods in Physics Research Section A 589 (2008) 133. [11] S. Antoine, Nuclear Instruments and Methods in Physics Research Section A 602 (2009) 293. [12] M.-J. Leray, et al., Magnetic measurements of the RACCAM prototype FFAG dipole, in: Proceedings of the Particle Accelerator Conference, Vancouver, BC, Canada, 2009, pp. 5032–5034. [13] A. Noda, et al., Design study of compact proton synchrotron for biomedical use, in: Proceedings of the European Particle Accelerator Conference, Berlin, Germany, 1992, pp. 1663–1665. [14] A. Morita, et al., A compact synchrotron with combined-function lattice dedicated for cancer therapy, in: Proceedings of the Particle Accelerator Conference, New York, USA, 1999, pp. 2528–2530. [15] A. Morita, et al., Field measurement of combined function magnet for medical proton synchrotron, in: Proceedings of the European Particle Accelerator Conference, Vienna, Austria, 2000, pp. 2110–2112. [16] G.H. Rees, Extraction, in: Proceedings of the CERN Accelerator School, CERN 85-19, Paris, France, 1984, pp. 346–357. [17] Y. Giboudot, et al., Particle tracking studies using dynamical map created from finite element solution of the EMMA cell, in: Proceedings of the Particle Accelerator Conference, Vancouver, BC, Canada, 2009, pp. 3290–3292. [18] R. Fenning, et al., An FFAG transport line for the PAMELA project, in: Proceedings of the Particle Accelerator Conference, Vancouver, BC, Canada, 2009, pp. 3741–3743. [19] S.L. Sheehy, The EMMA accelerator and implications for hadron non-scaling FFAGs, in: Proceedings of HB2010, Morschach, Switzerland, 2010, pp. 82–85.