Towards Quantitative Simulations of High Power Proton Cyclotrons
Abstract
We describe a large scale simulation effort using OPAL (Object Oriented Parallel Accelerator Library), that leads to a better quantitative understanding of the existing PSI high power proton cyclotron facility. The 1.3 MW of beam power on target poses stringent constraints on the controlled and uncontrolled beam losses. We present initial conditions for the Ring simulation, obtained from the new timestructure measurement and the many profile monitors of the 72 MeV transfer line. A trim coil model is developed, needed to avoid the dangerous resonance. By properly selecting the injection position and angle (eccentric injection), the flattop voltage and phase, very good agreement between simulations and measurements at the radial probe RRE4 is obtained. We report on orders of magnitude in dynamic range when comparing simulations with measurements. The relation between beam intensity, rms beam size, and accelerating voltage is studied and compared with measurement. The demonstrated capabilities are mandatory in the design and operation of the next generation high power proton drivers. In an outlook we discuss our future plans to include more physics into the model, which eventually leads to an even larger dynamic range in the simulation.
pacs:
29.20.dg;29.27.Bd;41.85.EwI Introduction
PSI operates a cyclotron based high intensity proton accelerator routinely at an average beam power of 1.3 MW. With this power the facility is at the worldwide forefront of high intensity proton accelerators. An upgrade program is under way to ensure high operational reliability and push the intensity to even higher levels. The beam current is limited in practice by losses at extraction and the resulting activation of accelerator components. Further intensity upgrades and new projects aiming at an even higher average beam power, are only possible if the relative losses can be lowered in proportion, thus keeping absolute losses at a constant level.
Maintaining beam losses at levels allowing handson maintenance is a primary challenge in any high power proton machine design and operation. For a MW beam in the PSI Ring cyclotron this corresponds to a transmission of taking controlled and uncontrolled losses into account. In a 10 MW class machine we require the losses to be on the same level which is a challenging task and is asking for precise beam dynamics calculation. In consequence, predicting beam halo at these levels is a great challenge and will be addressed in this paper.
High power hadron drivers have being used in many disciplines of science and, a growing interest in cyclotron technology for high power hadron drivers has be shown very recently. Two very recent papers demonstrate this fact: 1) The search for violation in the Neutrino sector Conrad and Shaevitz (2010) calls ultimately for three machines in the megawatt range at an energy of 800 MeV. 2) in Abderrahim et al. (2010), a white paper on Accelerator and Target Technology for Accelerator Driven Transmutation and Energy Production the cyclotron technology is advertised, quote: On the whole, the development status of accelerators is well advanced, and beam powers of up to 10 MW for cyclotrons and 100 MW for linacs now appear to be feasible …..
This report will briefly introduce OPAL, a tool for precise beam dynamics simulations including 3D space charge. One of OPAL’s ”flavors” (OPALcycl) is dedicated to high power cyclotron modeling and is explained in greater detail. We then explain how to obtain initial conditions for our PSI Ring cyclotron which still delivers the world record in beam power of 1.3 MW in continuous wave (cw) operation. Several crucial steps are explained, necessary to be able to predict tails at the level of in the PSI Ring cyclotron. We compare our results at the extraction with measurements, obtained with a MW cw production beam. Based on measurement data, we develop a simple linear model to predict beam sizes of the extracted beam as a function of intensity and confirm the model with simulations. A conclusion and discussions to include more physics into the model, which eventually leads to a even larger dynamic range in the simulation, closes the paper.
Ii Basic Equations and Physical Model
ii.1 A Brief Look at Opal
OPAL (Object Oriented Parallel Accelerator Library) is a tool for chargedparticle optic calculations in large accelerator structures and beam lines including 3D space charge. OPAL is built from first principles as a parallel application, it admits simulations of any scale: on the laptop and up to the largest High Performance Computing (HPC) clusters available today. Simulations, in particular HPC simulations, form the third pillar of science, complementing theory and experiment. OPAL includes various beam line element descriptions and methods for single particle optics, namely maps up to arbitrary order, symplectic integration schemes and lastly time integration Adelmann et al. (2008). OPAL is based on IPPL (Independent Parallel Particle Layer) Adelmann (2009) which adds parallel capabilities. Main functions inherited from IPPL are: structured rectangular grids, fields and parallel FFT and particles with the respective interpolation operators. Recently we added a powerful iterative solver to OPAL taking into account complicated boundary conditions Adelmann et al. (2010). More details on cyclotron modeling which are direct relevant to this article can be found in Yang et al. (2010). Several flavors of OPAL are available. For details we refer to the User Manual Adelmann et al. (2008). In this paper we use OPALt for the tracking of 72 MeV beam line, connecting two cyclotrons, the Injector 2 and the Ring Cyclotron. The other OPAL flavor  OPALcycl  is designed specially for cyclotron beam dynamics and, is explained in the next section.
ii.2 The Beam Dynamics Model of Opalcycl
In the cyclotrons and beam lines under consideration, the collisions between beam particles can be neglected because the typical bunch densities are low. In time domain, the general equations of motion of a charged particle in electromagnetic fields can be expressed by
where are rest mass, charge and the relativistic factor. With we denote the momentum of a particle, is the speed of light, and is the normalized velocity vector. In general the time () and position () dependent electric and magnetic vector fields are written in abbreviated form as .
If is normalized by , Eq. (II.2) can be written in Cartesian coordinates as
(1)  
The evolution of the beam’s distribution function can be expressed by a collisionless Vlasov equation:
(2) 
where and include both external applied fields, and space charge fields
(3) 
In order to model a cyclotron, the external electromagnetic fields are given by measurements or by numerical calculations.
The space charge fields can be obtained by a quasistatic approximation. In this approach, the relative motion of the particles is nonrelativistic in the beam rest frame, so the selfinduced magnetic field is practically absent and the electric field can be computed by solving Poisson’s equation
(4) 
where and are the electrostatic potential and the spatial charge density in the beam rest frame. The electric field can then be calculated by
(5) 
and back transformed to yield both the electric and the magnetic fields, in the lab frame, required in Eq. (II.2) by means of a Lorentz transformation. Because of the large vertical gap in our cyclotron, the contribution of image charges and currents are minor effects compared to space charges Baartman (1995), and hence it is a good approximation to use open boundary conditions.
The combination of Eq. (2) and Eq. (4) constitutes the VlasovPoisson system. In the following, the method of how to solve these equations in cyclotrons using PIC methods is described in detail.
Considering that particles propagate spirally outwards in cyclotrons, and the longitudinal orientation changes continuously, three righthanded Cartesian coordinate systems are defined, as shown in Fig. 1. The first coordinate system is the fixed laboratory frame , in which the external field data is stored and the particles are tracked.
Its origin is the center of the cyclotron and its plane is the median plane and the positive direction of axis points vertically upwards.
The second coordinate system is the local instantaneous frame , which is a temporal auxiliary frame for the space charge solver. Its origin is the mass center of the beam and the orientation of the axis is coincident with the average longitudinal direction and the positive orientation of the axis points vertically upwards.
The third coordinate system is the beam rest frame , which is comoving with the centroid of the beam. It has the same orientation and origin as , but the length in longitudinal direction is scaled by due to relativistic effects.
At each time step, the frames and are redefined according to the current 6D phase space distribution, and all particles are transformed from to , then a Lorentz transformation is performed to transform all particles to . The Poisson equation is then solved in the frame . In a 3D Cartesian frame, the solution of the Poisson equation at point can be expressed by
(6) 
with the 3D Green function
(7) 
assuming open boundary conditions. Details of the space charge field calculation can be found in Hockney and Eastwood (1988).
The model of the external magnetic field is based on midplane field measurements with excited trim coils. In consequence we have a vertical field, , measured on the median plane () as a function of azimuthal position (). Since the magnetic field outside the median plane is required to compute trajectories with , the field needs to be expanded in the direction. According to the approach given by Gordon and Taivassalo Gordon and Taivassalo (1985), by using a magnetic potential and the measured on the median plane at the point in cylindrical polar coordinates, the 3 order field can be written as
(8)  
where and
(9)  
All the partial differential coefficients are computed on the median plane data by interpolation, using Lagrange’s 5point formula.
Finally both the external fields and space charge fields are used to track particles for one time step using a 4 order RungeKutta (RK) integrator, in which the fields are evaluated for four times in each time step. Space charge fields are assumed to be constant during one time step, because their variation is typically much slower than that of external fields. More details and unique features can be found in Yang et al. (2010).
Iii Obtaining Initial Conditions for the Ring Cyclotron
At the extraction region of the Injector 2 we only have a very limited number of measurement data, however in the injecting beam line connecting the two cyclotrons we have vertical and horizontal beam profile monitors available for high intensity operation. Three timestructure measurements, one at the last turn of the Injector 2, one 27 meters downstream and one at the first turn of the PSI Ring cyclotron give important information on the longitudinal beam size Dölling (2010). In an overview (Fig 2) the starting point of the simulations, and some of the diagnostics are shown. We note that from a beam dynamics point of view, the particles travel in the order of km, from the marked start of the simulation to the RRE4, the probe covering the last 9 turns of the PSI Ring cyclotron.
We start the OPALt simulations from the middle of the last valley before extraction from the Injector 2 and perform a full 3D simulation until the magnetic injection channel (MIC) of the PSI Ring cyclotron. In Tab. 1 the initial values for the simulation of important beam parameters are shown. At the MIC we resample the distribution and switch to OPALcycl for the PSI Ring cyclotron simulation. The new distribution is sampled using the moments obtained from the transfer line simulation (at MIC).
Distribution  (mmmrad)  (mmmrad)  (mm)  (mm)  (mm)  (%)  

3D Gaussian  2.22  0.43  2.9  0.5  6.2  0.06  0.14  0.07  0.92 
Figure 3 shows the comparison of the beam width between OPALt and the measurements. The model predicts very well the evolution of the envelope from the beginning to the end of the transfer line.
Figure 4 shows the comparison of the predicted bunch length by the model and the measurements using the timestructure probes. The large error bar at the
injection to the Ring is because of the large background. The longitudinal initial conditions, in the center of the valley between sector magnet 3 and 4 (see Fig: 2), are derived from the time structure measurement, which is m upstream. In the center of the valley, the major axis of the bunch ellipse is along the longitudinal direction, and hence we obtain the transformed initial condition by a simple rotation in the longitudinal and radial plane.
Iv Towards Realistic High Power Cyclotron Simulations
The beam losses during the operation of the cyclotron usually limits the intensity that can be extracted. The PSI 590 MeV Ring routinely delivers 2.2 mA of cw beam, having a very low integrated loss rate, of the order of 0.02%. This thight margin avoids excessive activation of accelerator components and hence keeps the radiation dose imposed on the personnel involved in maintenance at acceptable levels. Furthermore, about % of the losses are located at injection and extraction.
Therefore, the understanding of the beam dynamics and, the knowledge of the detailed beam distribution especially at the extraction region, is one of the key points to be addressed especially if power levels increase in future projects Abderrahim et al. (2010).
Several important effects which need to be carefully modeled to keep extraction losses in the order of 0.02% are:

the turn separation at the position of the extraction septum must be made as large as possible,

the radial beam size at the extraction region must be smaller than the turn separation,

the halo, especially at the extraction, has to be minimized,

in case of the PSI Ring cyclotron, a long ”pencil” beam is used and hence the linear space charge effects must be effectively compensated to avoid the formation of a Sshaped beam which apparently increases the effective radial beam size.
We now discuss these issues related to the PSI Ring cyclotron which however can be considered universal for high power cyclotrons, and hence are certainly important for future high intensity related projects Abderrahim et al. (2010).
iv.1 The Flattop Phase
Although a compact beam is observed at the extraction of the Injector 2 cyclotron, the bunch length increases from about to about at injection into the Ring after passing through the almost 60 m long (72 MeV) transfer line. For such a long ”pencil” beam, a flattop cavity is needed to compensate the energy difference from the main cavity and avoid the formation of the Sshape beam caused by space charge effects.
When there is no space charge effect, the ideal flattop makes the total energy gain of any particle almost the same independent of the RF phase. Considering a high current beam, the flattop phase must be shifted such that the tail particles gain more energy than the head particles. This compensates exactly the linear part of the space charge force. Therefore the phase of the flattop is adjusted intensitydependent and, there exists an optimum flattop phase for a given intensity.
For our simulation we use % of the sum of the main cavity voltages as the flattop cavity voltage (as set in the control room) and adjust the phase to obtain the same phase difference between main and flattop cavity as set in the control room.
iv.2 The Effect of the Trim Coil TC15
A small manufacturing error produces a slight deviation in the average field profile and a corresponding shift in the tunes and . This requires a strong excitation of the trim coil TC15. Without this correction, the coupling resonance would be crossed four times at energies of 490, 525, 535 and 585 MeV, respectively. A large horizontal oscillation would be transformed into a large vertical one at the coupling resonance which can lead to large vertical beam losses. An analytic model was developed which mimics the field due to real trim coil characteristics Adam and Joho (1974). It is described by Eq. (10).
(10) 
where R is the radius, is the maximum magnetic field, the constants , , , , , , for our case. It provides an additional magnet field and field gradient in the radial direction as shown in Fig. 5.
The trim coil gives a maximum magnetic field of 14 Gauss, and furthermore has a long tail towards smaller radii in order to make the integrated strength of the trim coil over the radius to zero. The radial and vertical tune shifts caused by TC15 are given by,
(11) 
where is the orbit radius, is the hill field, is the average field gradient in radial direction. Careful beam dynamics studies have shown the meaningfulness of such detailed modeling in order to obtain a complete and precise pictures of the beam dynamics in the PSI Ring cyclotron. The modified tune diagram by TC15 is shown in figure. 6.
Without TC15, in simulations and in the operation of the Ring, we observe severe vertical beam losses and can not obtain the required extraction efficiency.
iv.3 The Injection Position and Angle
The PSI Ring cyclotron has a single turn extraction, hence a large radial turn separation between the last two turns is required. The turn separation for a centered beam is defined as
(12) 
where is the field index. For the PSI Ring cyclotron this gives about mm (Fig. 7 upper part) at the extraction region, which is not enough for high current operation and would result in large losses.
To increase the turn separation, a noncentered injection into the PSI Ring cyclotron is used. Since at extraction, adjusting the injection position and angle, results in the betatron amplitude being almost equal to the increase in radius per turn. The formation of the turn pattern under this condition, for the last nine turns, is shown in Fig. 7 (lower part).
This is a special turn pattern because the last turn is well separated from the overlapping second, third and fourth last turns. In this case, the turn separation at the extraction turn is as large as mm, hence it allows the extraction of a high intensity beam with very low losses.
iv.4 Comparing the Radial Intensity Profile at Extraction with Measurements
Up to now we have described the most important steps in setting up a precise beam dynamics simulation of the PSI Ring cyclotron. We now compare simulations with measurements from a radial probe (RRE4) covering the last turns of the PSI Ring cyclotron. The probe is located cm upstream from the m thick electrostatic extraction septum and, hence gives a very good picture of the beam distribution at the septum.
This probe is able to measure at the full intensity of the MW cw beam. In order to compare the simulations with measurements, not only a radial probe is implemented in OPALcycl but also all other parameters, described in the previous sections, can be entered into the simulation.
The flattop phase and the injection position and angle are optimized to get the largest turn separation and smallest beam size at the extraction region, in both the simulation and operation of the PSI Ring cyclotron.
The effect of the trim coil TC15 on the turn pattern is shown in Fig. 8, black denotes the measurement and the colors distinguishing simulations with and without TC15 . For a fixed energy the shift is given by . For turn 180 this shift is: , hence the center of turn 180 moves to the exact position of the measurement when considering the effect of the TC15.
In the PSI cyclotron facility, the beam is heavily collimated during the early stage of acceleration in the Injector 2 and in the beam transfer line to the PSI Ring cyclotron. As a consequence, the beam profiles do not follow a Gaussian distribution. In Fig. 9 we show again the intensity pattern at the last 9 turns but for different starting distributions at MIC. Using the properties of a binomial distribution Joho (1980), we can vary a single parameter from 0 to and cover a wide range of distributions: from KV to Gaussian. We find that a parabolic distribution () matches best the measurement at the crucial point, the septum. On the extreme side, as expected, the Gaussian distribution (truncated at ) with its tails would fill up the intensity dip and hence would increase the losses at the septum. This indicated that there are indeed sharp edges in the real distribution, which still differs from the assumed idealized distribution, as suggested by the remaining differences between simulation and measurement (Fig 9, e.g at ).
Nevertheless, this remarkable agreement is obtained after km of tracking the millions of macro particles through external and self fields. This is only possible because of the parallel nature of OPAL which allows such simulations to be carried out on large high performance computing clusters. The statistical error of the measurement is indicated at mm in Fig. 8 and 9. At other radii the error bars are significantly smaller and hence not shown in the figures. The statistical errors of the simulation are smaller than those of the measurement errors due to the large number of simulation particles.
iv.5 Scaling Law of Beam Size With Respect to Current
The energy spread caused by the linear longitudinal space charge force after revolutions is given by Joho Joho (1981) using the sector model:
(13)  
where (form factor), and is the average current, is the phase width, the turn number and , where is the final velocity of the beam.
The linear energy spread can be compensated with a tilted flattop voltage. This reduces the energy of leading particles and increases the energy of trailing particles. There remains a non linear part of the energy spread , which can not be compensated. Let’s define
(14) 
and note that depends strongly on the beam distribution and is in our case an open parameter in the range of .
According to Eq. 12, the space charge induced energy spread leads to a radial spread () that results in an increase in beam size. In Fig. 10 we compare beam sizes at the extraction for beam currents from 10 to 2.2 mA with simulations.
Even though the measurements where done over a time span of 4 years with very different machine configurations we obtain a good agreement between the simulations (theory) and the measurements. Hence we can predict very well the extracted beam size as function of the intensity. This simple model and the precise simulations shown in the previous paragraph constitutes a benchmarked model for the prediction of the most delicate parameters in high intensity cyclotrons.
The relation between the average energy gain and the turn number is:
(15) 
where is the final energy and is the initial energy. For single turn extraction the loss on the septum is limited by the ratio
leading to the condition
(16) 
We obtain empirically, for a centered beam a value for which is , where as for an eccentric beam, the turn separation is enhanced, and hence can be as high as .
Putting (13), (14), (15) and (16) together, we get for the current limit from longitudinal space charge forces
(17) 
where and . Since the turn number is inverse proportional to the cavity voltage , we see the big advantage of a large cavity voltage
(18) 
as experimentally demonstrated by the historical current development in the PSI Ring cyclotron shown in Fig 11.
This relation was first predicted by Joho Joho (1981) and confirmed during almost 36 years of operation of the PSI Ring cyclotron.
The prediction of the current limit in the PSI Ring cyclotron gives with (17) , using the following parameters: MV, MV, , turns ( MeV), and estimating and (eccentric injection).
This is remarkably close to the present current limit of mA (2010), given the crude assumptions:

no turn structure inside the charge sheet (see Joho (1981) figure 3),

non relativistic approximation,

no radial boundary condition and

uncertainties in and .
We note that for the PSI Injector 2 (17) is not applicable due to phase mixing Adam (1985) and hence would give a pessimistic value for the current limit.
V Conclusions and Discussions
In this paper, we present novel precise simulations for the beam dynamics in high intensity cyclotrons. For the first time we are able to obtain a realistic and detailed understanding of the beam dynamics in the very complex PSI Ring cyclotron by means of 3D particle simulations. By a rough estimation of the initial distribution, according to measurements of beam profile monitors, and the timestructure of the beam, realistic simulations of the PSI Ring cyclotron are presented and compared to measurements.
Very good agreement for the radial probe between the simulation and measured data is obtained by adjusting the injection position, angle, flattop voltage, and the trim coil TC15. These parameters are all in agreement with settings obtained from the control room.
The presented results with a level of detail large enough to predict limiting tails on the extraction septum (at beam width levels of ), and can be seamless extrapolated to future high power cyclotrons and enable the precise prediction of crucial parameters, such as losses, based on an existing cw megawatt facility experiences. However a crucial part is the knowledge of the initial distribution from the evaluation of measurements of beam profiles and time structure.
Primary modeling limitations include an accurate knowledge of the initial particle distribution in the full 6D phase space, and the lack of particlematter interaction in our model. ParticleMatter interaction models and resulting struggled primary particles and electrons will play an important role when intensity levels increase while at the same time, the losses must be held at present levels. We plan to include these effects in future studies, preliminary results on the particlematter interaction model are reported in Bi et al. (2010) and ideas for secondary electron creation and field emission can be found in Wang et al. (2010).
Vi Acknowledgments
The authors thank the Accelerator Modeling and Advanced Simulation (AMAS) group members C. Kraus, Y. Ineichen and J. J. Yang for many discussions regarding programming and T. Schietinger for providing the postprocessing tool H5PartRoot. We also thank H. Zhang for providing information of the 72 MeV injection line and the PSI Ring cyclotron. This work was partly performed on the felsim cluster at the Paul Scherrer Institut and on the Cray XT5 at Swiss National Supercomputing Center (CSCS) within the “Horizon” collaboration.
References
 Conrad and Shaevitz (2010) J. M. Conrad and M. H. Shaevitz, Phys. Rev. Lett. 104, 141802 (2010).
 Abderrahim et al. (2010) H. A. Abderrahim et al., Tech. Rep., DOE (2010), Accelerator and Target Technology for Accelerator Driven Transmutation and Energy Production, URL http://www.science.doe.gov/hep/files/pdfs/ADSWhitePaperFinal.%pdf.
 Adelmann et al. (2008) A. Adelmann, Y.Bi, C. Kraus, Y. Ineichen, S. Russel, and J. Yang, Tech. Rep. PSIPR0802, Paul Scherrer Institut (2008).
 Adelmann (2009) A. Adelmann, Tech. Rep. PSIPR0905, Paul Scherrer Institut (2009).
 Adelmann et al. (2010) A. Adelmann, P. Arbenz, and Y. Ineichen, Journal of Computational Physics 229, 4554 (2010), ISSN 00219991, URL http://www.sciencedirect.com/science/article/B6WHY4YHP08T1/%2/41309c23eb7fa1b4af95d9401a21da39.
 Yang et al. (2010) J. J. Yang, A. Adelmann, M. Humbel, M. Seidel, and T. J. Zhang, Phys. Rev. ST Accel. Beams 13, 064201 (2010).
 Baartman (1995) R. Baartman, in Proc. 14th Int. Conf. on Cyclotrons and their Applications (Capetown, 1995), p. 440.
 Hockney and Eastwood (1988) R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (Hilger, New York, 1988).
 Gordon and Taivassalo (1985) M. M. Gordon and V. Taivassalo, IEEE Trans. Nucl. Sci. 32, 2447 (1985).
 Dölling (2010) R. Dölling, in Proc. of HB2010 (Morschach, Switzerland, 2010), ”MOPD62”.
 Adam and Joho (1974) S. Adam and W. Joho, Tech. Report TM1113, PSI (1974).
 Joho (1980) W. Joho, Tech. Report TM1114, PSI (1980).
 Joho (1981) W. Joho, in 9th Int. Conf. on Cyclotrons, p. 337 (Caen, 1981).
 Adam (1985) S. Adam, IEEE Trans. on Nuclear Science 32, 2507 (1985).
 Bi et al. (2010) Y. J. Bi, A. Adelmann, R. Dölling, W. Joho, M. Seidel, C. X. Tang, and T. J. Zhang, in Proc. of HB2010 (Morschach, Switzerland, 2010), ”TUO2A03”.
 Wang et al. (2010) C. Wang, A. Adelmann, and Y. Ineichen, in Proc. of HB2010 (Morschach, Switzerland, 2010), ”MOPD55”.