Fixed field phase shifters for a multipass recirculating superconducting proton linac

The multi-pass recirculating proton linac can significantly improve the usage efficiency of rf superconducting cavities by accelerating the proton beam through the same cavity multiple times. However, in order to achieve acceleration for different energy proton beams, longitudinal phase synchronous conditions among these beams have to be satisfied. In this paper, we propose a fixed field phase shifter system based on two superconducting bending magnets with an optimized field distribution to meet the synchronous conditions. Such a phase shifter provides both longitudinal phase delay and transverse focusing of multiple energy beams. Self-consistent macroparticle simulations including space-charge effects through the system show that the proton beam emittances can be well preserved after the phase shifter for multiple energy beams.


I. INTRODUCTION
Superconducting proton linacs have been constructed or proposed as a driver for high energy neutrino physics study [1,2], for production of spallation neutron sources [3,4] and tritium [5], and for nuclear waste transmutation and energy production [6][7][8][9]. However, those linacs are expensive due to the inefficient usage of superconducting rf cavities with single pass of the proton beam. Recently, we proposed a concept of recirculating superconducting proton linac that has the potential to substantially save accelerator project costs by reusing the same superconducting rf cavities multiple times to accelerate the beam to multiple GeVs [10,11]. In that proposed example, 17 superconducting cavities with 10.3 MeV energy gain per cavity was assumed to accelerate the proton beam from 150 MeV to 500 MeV in two passes, 39 cavities with 9.6 MeV energy gain per cavity to accelerate beam from 500 MeV to 2 GeV in 4 passes, and 50 cavities with 20 MeV energy gain per cavity to accelerate the proton beam from 2 GeV to 8 GeV. This new concept can drastically reduce the number of superconducting cavities needed to reach the final 8 GeV beam energy from near 500 to around 100. It also avoids the potential intensity and repetition rate limits of the rapid cycling synchrotron [12].
Beam dynamics design studies [13,14] have been done for the first section of proposed multi-GeV recirculating linac using the IMPACT code [15,16]. The separation between superconducting cavities is adjusted to satisfy the longitudinal synchronous condition so that proton beams with two different energies can be accelerated using the same rf cavity. However, for multiple passes (>2), multiple longitudinal synchronous conditions among these passes need to be satisfied for fixed design phase acceleration through rf cavities. The phase slippage can no longer be compensated by a single control variable of space between cavities and requires a new approach. An energy dependent path length change during recirculation is needed.
In the previous study, a chicanelike phase shifter concept using three bending magnets was proposed for multipass recirculating proton linacs [17]. However, that concept requires discontinuous change of magnetic field distribution, which is not practically achievable. Moreover, using three bending magnets between two rf cavities is not cost effective. An additional challenge for the design of phase shifting magnets comes from continuous wave (CW) operation. In such a mode it is not possible to ramp the magnetic field fast enough, requiring that the magnets have a momentum acceptance large enough to cover the full range of energies in the recirculating section with fixed magnetic field. In this study, we propose a fixed field achromatic superconducting magnetic system that consists of two superconducting bending magnets providing longitudinal phase shifts for different energy beams. A schematic plot of this phase shifter is shown in Fig. 1. Here, a proton beam passing through the rf cavity will enter the bending magnet field region at 45 degrees with respect to the field's entrance edge. In the middle plane of the magnet, there is only a vertical component of the field that varies only as a function of x, the distance perpendicular to the entrance edge. With sufficiently long magnet, the magnetic field is only a function of x and y, and there is no z component of magnetic field in the bending magnet. Given this magnetic field distribution, a proton beam entering the bending magnet with a given angle will exit from the magnet with the same angle. In this case, with 45 degree entry angle, the proton beam will exit the magnet at another 45 degrees with respect to the field entrance edge. The entire magnet provides 90 degree bending of the proton beam. A proton beam with different energies after the rf cavity entering the first bending magnet from the same path will be separated out inside the magnet and exit the magnet with different horizontal separations. Using the symmetry as shown in the dashed line of the above plot, those separated beams entering the second bending magnet, that has the same magnetic field distribution as the first bending magnet, will be bent for another 90 degrees and merge into a single straight path at the exit of the magnet fringe field edge due to the symmetric condition. These two 90 degree bending magnets form an achromatic system for different energy beams. This approach using the symmetry of bending magnets to achieve large momentum acceptance has also been used for electron accelerators [18], and has recently been considered for superconducting gantries for proton therapy [19]. By appropriately designing the field distribution along the direction perpendicular to the entrance edge, i.e., B y ðxÞ, the time differences of multiple proton beam passes through this phase shifter can be adjusted to satisfy the longitudinal synchronous conditions so that those proton beams with different energies can be accelerated by the same rf cavity for multiple times.
The organization of this paper is as follows: after the Introduction, we discuss the computational model used in this study in Sec. II; we then give several beam dynamics design examples of phase shifters in Sec. III; and draw conclusions in Sec. IV.

II. COMPUTATIONAL MODEL
In the multi-pass recirculating superconducting proton linac, the proton beam passing through the same cavities separated by a phase shifter has to satisfy a phase synchronous condition in order to be accelerated during the multiple passes. Let t im denote the time elapsed between two rf cavities i and i þ 1 during the mth beam pass of the accelerator, t in the time elapse between the two cavities  FIG. 3. Horizontal centroid (hxi) trajectories of four energy passes through the phase shifter using the field distribution from the Table I. during the nth beam pass. If the difference of the two time elapses satisfies the following synchronous condition: where T rf is the oscillation period of the rf field inside the cavity, the proton beam will see the same design phase of the cavity during multiple passes of the accelerator. If one would like to have a different design phase on each pass, the above longitudinal synchronous condition can be adjusted correspondingly for each pass to attain the desired phase.
In addition to satisfying the above longitudinal synchronous condition, the proton beam also needs to be transversely focused to compensate the large vertical defocusing effect of the magnet edge due to entering at a 45 degree angle, and beam emittances should be preserved after transporting through the phase shifter. In order to meet both longitudinal and transverse requirements, the magnetic field distribution inside the bending magnet has to be carefully designed. In this study, we assume that the middle plane vertical magnetic field B y ðxÞ can be approximated using multiple Gaussian wavelets as: The coefficients a i , b i , and c i will be adjusted during the optimization process to obtain a field distribution that yields both correct longitudinal phase shifts among multiple energies and appropriate transverse focusing of the beam. From the middle plane magnetic field B y ðxÞ, the general off-plane expression of the magnetic fields inside the bending magnet can be written as [20]: where ðx; y; zÞ refers to the local coordinate system inside the bending magnet as shown in Fig. 1.
Extensive numerical optimizations were employed to search for appropriate magnetic field distribution B y ðxÞ. The trajectory of a charged particle inside the bending  (Table I)  magnet can be obtained from the solution of the equations of motion: where r ¼ ðx; y; zÞ T is the spatial position vector, v is the velocity vector, z, c is the speed of light in vacuum, m is the rest mass of the particle, q is the charge of the particle. The single proton trajectory (i.e., reference particle trajectory) inside the middle plane can be obtained using a leap-froglike method to solve the above equations of motion. The advances of particle coordinates and velocities after one time step τ are In the above scheme, the positions are updated twice in one time step, and the velocities are updated once in one time step. We update only one velocity component directly and make use of energy conservation of a charged particle inside the magnetic field (neglecting synchrotron radiation and wakefields) to update the other component. This results in an exact energy conservative numerical integrator.
In order to monitor the proton beam transverse sizes during the bending magnet field distribution optimization, we use a linear beam dynamics model to compute the beam root mean square (RMS) size evolution through the phase shifter quickly. Given the middle plane bending magnet  (Table I)  field distribution B y ðxÞ, we track a single particle using this field following the above numerical scheme. From the trajectory of this particle, one can obtain the distance L z along the z axis that the particle travels from the entrance to the exit of the bending magnetic field. Next, we calculate the linear transfer matrix through the bending magnetic field along this reference trajectory using a truncated power series algebra (TPSA) method [21]. Here, we use z in bending magnet coordinate as shown in Fig. 1 as the independent variable and define a TPSA particle with coordinates ðX; V x ; Y; V y ; T; V z Þ. Each coordinate of the TPSA particle denotes a TPSA variable that contains the variable value and its first, second, and higher order derivatives with respect to the independent variable z. For example, for a variable that is a one-dimensional function fðzÞ, its TPSA representation at z 0 can be written as: where superscript denotes the derivative with respect to the variable z. For a constant c, its TPSA representation is Dc ¼ ½c; 0; 0; …; 0, for a variable z, Dz ¼ ½z; 1; 0; …; 0.
With appropriate definition of addition rule and multiplication rule of the TPSA variable, the calculation of derivatives becomes algebra operations without explicit function derivatives. This significantly speeds up the numerical calculation of those derivatives on a computer by using an available library [22]. The linear and high order transfer matrix between two locations along z can be obtained from those derivatives. We applied the similar leap-frog method for approximating the solution of equations of motion. Following the numerical solution of the equations of motion of a particle inside the magnetic field, the new coordinates and velocities of the TPSA particle after one position step τ z are: The linear transfer matrix through the bending magnet is obtained after tracking this TPSA particle through the bending magnetic field. Using the linear transfer matrices through the drift and the bending magnet, the evolution of the transverse RMS sizes of the proton beam between two locations z 1 and z 2 can be calculated as: where Σ is the moment matrix of the particle relative coordinates (with respect to the reference particle tra- R is the linear transfer matrix between location 1 and 2, and superscript T denotes the transpose of the matrix.
Using this linear dynamics model, the proton beam sizes can be quickly computed at several locations of the phase shifter.
The above beam dynamics tracking model is integrated into an optimization package for the longitudinal phase optimization. This optimization package uses a combination of a global stochastic differential evolution method and a local deterministic Simplex optimization method. The differential evolution method uses a group of potential solutions for global search of optimal solution. After a number of generations, when the solution is near the optimum, we switch to the Simplex optimizer from the numerical recipe [23] for fast convergence to the solution. Using the combination method greatly improves the speed of finding the solution that satisfies the longitudinal synchronous conditions.

III. APPLICATION EXAMPLES
We applied the above computational model to the beam dynamics design of two multipass phase shifters in a superconducting recirculating proton linac. In this application, we assume that a proton beam enters the recirculating linac with an initial kinetic energy of 150 MeV, and exits from the linac with a final energy of 600 MeV. The linac consists of 10 superconducting rf cavities operating at 650 MHz with an average energy gain of 11.25 MeV per cavity. The proton beam passing through the same rf cavity four times with an energy separation of 112.5 MeV per pass before exiting the linac. Such a linac can be used as an injector of the accelerator driven system study resembling the one in MYRRHA [7].
A. Phase shifter at the entrance of the linac In the first example, we would like to find the middle plane vertical magnetic field distribution of the bending magnet in a phase shifter that can provide appropriate longitudinal phase shift and transverse focusing for four proton beam passes with 150 MeV, 262.5 MeV, 375 MeV, and 487.5 MeV kinetic energies. Such a constant energy gain per pass is based on the assumption of average energy gain and can have somewhat variation in a real accelerator design. This shifter corresponds to the phase shifter at the entrance of the linac. During the optimization of the middle plane magnetic field distribution, we set a constraint of maximum 5 mm on the proton beam transverse RMS size at the exit of the first bending magnet for the first three passes and after the second bending magnet for the fourth pass. These constraints are based on the assumption that there would be sufficient drift space between the two bending magnets for the first three passes so that quadrupole magnets can be inserted to control the transverse beam sizes inside and after the second magnet. For the fourth pass, the separation between the two bending magnets could be too close to insert quadrupoles between them. The time differences among those four passes are optimized with respect to the parameters used to define the field distribution in Eq. (2) in order to meet the synchronous condition Eq. (1). Here, we used four Gaussian wavelets to approximate the middle plane vertical magnetic field distribution B y ðxÞ. The analytical expression of this field is given as: where the coefficients in these Gaussian functions after the optimization are given in Table I. Figure 2 shows the middle plane vertical magnetic field distribution from the optimization (Table I) after the longitudinal synchronous conditions were attained for all four passes. The peak magnetic field is about 4 Tesla in this case. Given the middle plane vertical field, we can calculate the threedimensional magnetic field components using Eq. (5).
Using these fields, we carried out multiparticle tracking simulation using a parallel particle tracking code, IMPACT-T. [16]. We assumed an initial Waterbag proton beam distribution with 1 mm RMS size and 1 mm mrad RMS  (Table II). emittance in all three dimensions. Figure 3 shows horizontal centroid (hxi) trajectories through the phase shifter for the four energy passes. Here, we have used a local coordinate system for each pass. The simulation starts with L 1 ¼ 0.6 meter drift before the bending magnet. This choice of drift length is not critical as shown in another example. The distance between the first dipole entrance point and the center symmetry line L 2 ¼ 2.39 meter is obtained from the longitudinal optimization. At the entrance of the bending magnet field (edge of the fringe field), the horizontal plane X − Z coordinates are rotated 45 degrees into the bending magnet coordinates as shown in Fig. 1. At the exit of this magnet field, the X − Z coordinates are rotated another 45 degrees into the horizontal drift coordinates between the two bending magnets. At the entrance of the second bending magnet field, the coordinates are rotated another 45 degrees into the second bending magnet coordinate system. At the exit of the second bending magnet field, the coordinates are rotated another 45 degrees into the vertical down drift coordinates. The difference of the coordinate system between the initial up drift and final down drift is 180 degrees. From Fig. 3, it is seen that the proton beams with different energies are separated into different trajectories inside the bending magnets. The separation between the two trajectories in the drift between two bending magnets is about half a meter. These separations in the drifts between the two bending magnets allow one to insert separate quadrupoles for each beam energy between the bending magnets to provide extra transverse focusing for each energy pass. In this example, we inserted four quadrupoles in the drift between the two bending magnets for the 150 MeV energy pass and for the 262.5 MeV energy pass separately. The strength of those four quadrupoles are adjusted to avoid the final transverse emittance growth while keeping reasonable transverse beam sizes at the end of the phase shifter. Figure 4 shows the proton beam RMS emittance evolutions in all three dimensions through the phase shifter with 150 MeV, 262.5 MeV, 375 MeV, and 487.5 MeV kinetic energies. It is seen that the initial emittances in three dimensions for all four proton beam energies are well preserved after the phase shifter. There is little emittance growth after the phase shifter. Figure (Table II)  (rotation in x-z plane) between different section of elements through the phase shifter as shown in Fig. 1.
In the above simulation, we assumed zero current and the space-charge effects were not included in the simulation. The nonlinear space-charge effects could cause beam emittance growth and even beam blow up. In order to evaluate the potential space-charge effects on the beam quality, we reran the simulation with an average 40 mA current proton beam at 650 MHz frequency. This corresponds to an average 10 mA current proton beam at 162.5 MHz that was proposed in a high intensity superconducting linac for the accelerator driven system application [8]. Figures 6 and 7 show the proton beam transverse and longitudinal emittance and RMS size evolutions through the phase shifter including the space-charge effects of 40 mA beam current. It is seen that the impact of the space-charge effects on the proton beam quality is small. There is no significant extra emittance growth or RMS size increase due to the space-charge effects.  (Table II) (Table III). In the above simulations, we assumed an initial Waterbag distribution. In order to check the impact due to different initial distributions on emittance growth, we reran the simulation using an initial Gaussian distribution with the same RMS values. Figures 8 and 9 show the proton beam transverse and longitudinal emittance and RMS size evolutions through the phase shifter with an initial Gaussian distribution including the space-charge effects of 40 mA current. The impact from different initial distributions is small. There is little extra emittance growth or RMS size increase due to the use of the initial Gaussian distribution.
In the above example, we assumed an initial 0.6 meter drift before entering the first bending magnetic field region. The use of this drift is not critical. In the following, we searched for a new solution that results in a different initial drift. In this example, the initial drift is 1.06 meters. We used three Gaussian functions to approximate the vertical magnetic field distribution in the middle plane. The analytical expression of this field is given as: The coefficients in these Gaussian functions were adjusted so that the time differences of four energy passes satisfy the longitudinal synchronous conditions (1). The resultant coefficients are given in Table II. The vertical magnetic field distribution in the middle plane is given in Fig. 10. This field has a similar distribution to the preceding example with a 0.6 meter initial drift. Given the middle plane vertical field, we can calculate the three-dimensional magnetic field components following the same Eq. (5). Using these fields and the middle drift length L 2 ¼ 1.9, we carried out self-consistent macroparticle tracking including the space-charge effects of 40 mA beam current. Figures 11  and 12 show the proton beam transverse and longitudinal emittance and rms size evolutions through this phase shifter. There is little emittance growth after the phase shifter for all four passes of the proton beam similar to the preceding example. The transverse RMS beam sizes are also reasonably controlled within a few millimeters after the phase shifter.

B. Phase shifter at the exit of the linac
In the above example, we did beam dynamics optimization of the first phase shifter for four proton beam passes FIG. 14. Proton beam transverse and longitudinal RMS emittance evolution through the phase shifter (Table III)  starting with the entrance 150 MeV energy. In the following example, using the above computational model, we obtained a solution for the four passes starting with an energy entering the last phase shifter. The proton beam passes through this phase shifter with a kinetic energy of 251.25 MeV, 363.75 MeV, 476.25 MeV, and 588.75 MeV. We used six Gaussian functions to approximate the vertical magnetic field distribution in the middle plane of the bending magnet: where the coefficients a i , b i , and c i are adjusted so that the time differences among the four passes satisfy the longitudinal synchronous conditions. We also put the same constraints on transverse beam size at the exit of the first bending magnets for the first three passes and the second bending magnet after the fourth pass. The resultant coefficients are given in Table III. The vertical magnetic field distribution in the middle plane of the bending magnet is also shown in Fig. 13. The peak of this field is around 10 Tesla. Such a high peak might result from the higher proton beam energy through this phase shifter. A higher magnetic field also helps keep the phase shifter smaller. Using this magnetic field, we carried out self-consistent macroparticle tracking simulation through the phase shifter for the four proton beam passes including the space-charge effects of 40 mA beam current. Figures 14 and 15

IV. CONCLUSIONS
In this study, we proposed a fixed field, two bending magnet phase shifter concept that can provide appropriate longitudinal phase delay and transverse focusing of different energy proton beams for a multi-pass recirculating FIG. 15. Proton beam transverse and longitudinal RMS size evolution through the phase shifter (Table III)  superconducting proton linac. This concept is illustrated with beam dynamics optimization design of two four energy pass phase shifters that yields realistic magnetic field distribution in the middle plane of the bending magnet. Detailed engineering design of those magnets can be done using an inverse Biot-Savart optimization method [24] and will be reported in a future publication. Self consistent macroparticle tracking simulations including the spacecharge effects of 40 mA current are carried out for proton beams transporting through these phase shifters. The proton beam emittances with four kinetic energies are well preserved after the phase shifters. Such a fixed field phase shifter for multiple energy beams forms an important beam line component in an efficient recirculating superconducting proton linac. In this paper, we have not included effects such as rf cavity failure and rf cavity coupler limit in the study. These effects and the detailed beam dynamics design of a recirculating linac will be investigated in the future study and reported in another publication.

ACKNOWLEDGMENTS
We would like to thank Drs. Z. Liu and Y. Hao for the TPSA library used in this study. This work was supported by the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and used computer resources at the National Energy Research Scientific Computing Center.

APPENDIX: UNIFIED DIFFERENTIAL EVOLUTION OPTIMIZATION METHOD
The differential evolution is a simple but powerful method for global parameter optimization [25][26][27]. Compared with the other evolutionary algorithms such as the genetic algorithm, the differential evolution algorithm makes use of the differences of parent solutions to attain gradient information. This helps improve the convergence speed of the algorithm in comparison to the genetic algorithm. Meanwhile, compared with the particle swarm method, the differential evolution algorithm has a crossover stage to enhance the diversity of solutions. This helps the differential evolution algorithm to avoid converging to a premature solution.
The differential evolution algorithm generates new offspring using two operations: mutation and crossover. During the mutation stage, for each population member (target vector) ⃗ x i , i ¼ 1; 2; 3; …; NP at generation G, a new mutant vector ⃗ v i is generated by following a mutation strategy. A number of mutation strategies have been proposed for the conventional standard differential evolution algorithm. The presence of multiple mutation strategies in the literature complicates the use of the differential evolution algorithm. Recently, we proposed a single mutation expression that can unify most conventional mutation strategies used by the differential evolution algorithm [28]. This single unified mutation strategy can be written as: Here, the second term on the right-hand side of equation (A1) denotes the contribution from the best solution found in the current generation, the third term denotes the rotationally invariant contribution from the random solution, and the fourth term and the fifth term are the same terms as those used in the original differential evolution algorithm to account for the contribution from the difference of parent solutions. Those last two terms divert the mutated solution away from the best solution and help to improve the algorithm's exploration of the decision parameter space. The four parameters F 1 , F 2 , F 3 , and F 4 are the weights from each contribution. This unified expression represents a combination of exploitation (using the best found solution) and exploration (using randomly chosen solutions) when generating the new mutant solution. The use of four control parameters in the above unified differential evolution scheme is to increase the diversity of the algorithm in search of the global optimal solution. On the other hand, this might also slow down the speed of convergence to the final optimal solution, which will be improved by switching to a faster deterministic Simplex optimizer after a number of generations.
Using the Eq. (A1), the multiple mutation strategies of the standard differential evolution algorithm can be included in a single expression. For example, the standard differential evolution algorithm such as DE/rand/1 and DE/rand/2 can be attained by setting F 1 ¼ 0, and F 2 ¼ 1. This new expression provides an opportunity to explore more broadly the space of mutation operators. Using a different set of parameters F 1 , F 2 , F 3 , F 4 , a new mutation strategy can be achieved. Moreover, by adjusting these parameters during the evolution, the multiple mutation strategies and their combinations can be used during different stages of optimization. Thus, the unified mutation expression has the virtue of mathematical simplicity and also provides users with flexibility for broader exploration of different mutation strategies.
A crossover operation between the new generated mutant vector ⃗ v i and the target vector ⃗ x i is used to further increase the diversity of the new candidate solution. This operation combines the two vectors into a new trial vector ⃗ U i ; i ¼ 1; 2; …; NP, where the components of the trial vector are obtained from the components of ⃗ v i or ⃗ x i according to a crossover probability Cr. In the binomial crossover scheme for a D dimensional control parameter space, the new trial vector ⃗ U i , i ¼ 1; 2; …; NP is generated using the following rule: where rand j is a randomly chosen real number in the interval [0, 1], and the index mbr i is a randomly chosen integer in the range ½1; D. This ensures that the new trial vector contains at least one component from the new mutant vector. During the mutation stage and the crossover stage at generation, each individual solution ⃗ x i , i ¼ 1; 2; …; NP has a set of control parameters F 1;i ; F 2;i ; F 3;i ; F 4;i and Cr i associated with it. Before generating a new mutant solution using the unified differential evolution expression (A1), a set of control parameters F 1;i ; F 2;i ; F 3;i ; F 4;i and Cr i are calculated as: Cr i ¼ Cr min þ r i ðCr max − Cr min Þ: ðA5Þ For each generation, where r ji ; r i ; j ¼ 1, 2, 3, 4 are uniform random values in the interval [0, 1], F j min and F j max for j ¼ 1, 2, 3, 4 are the minimum and the maximum allowed values of those control parameters, Cr min and Cr max are the minimum and the maximum crossover probability, and τ j ; j ¼ 1, 2, 3, 4, 5 represents the probability to use a new value or to keep the old value for the j th control parameter. The values of τ j are normally kept small so that better control parameters associated with surviving solutions will be reused to generate the new trial solution. In this study, we set τ j ¼ 0.1 following the Ref. [29]. The values of F j min and F j max are set to 0 and 1 respectively in this study. We also set Cr min ¼ 0.0 and Cr max ¼ 1. The selection of these values is based on the consideration that the various conventional differential evolution mutation strategies can be covered by the settings of those parameters, and in the literature, F 3 and F 4 are rarely greater than one.