Towards Quantitative Simulations of High Power Proton Cyclotrons

PSI operates a cyclotron based high intensity proton accelerator routinely at an average beam power of 1.3MW. With this power the facility is at the worldwide forefront of high intensity proton accelerators. The beam current is practically limited 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 hands-on maintenance is a primary challenge in any high power proton machine design and operation. In consequence, predicting beam halo at these levels is a great challenge and will be addressed in this paper. High power hadron driver have being used in many disciplines of science and, a growing interest in the cyclotron technology for high power hadron drivers are being observed very recently. This report will briefly introduce OPAL, a tool for precise beam dynamics simulations including 3D space charge. One of OPAL's flavors (OPAL-cycl) 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 continuous wave (cw). Several crucial steps are explained necessary to be able to predict tails at the level of 3\sigma ... 4\sigma in the PSI Ring cyclotron. We compare our results at the extraction with measurements, obtained with a 1.18 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 intensities and confirm the model with simulations.

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 hands-on maintenance is a primary challenge in any high power proton machine design and operation. For a 1.3 MW beam in the PSI Ring cyclotron this corresponds to a transmission of 99.97% 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 CP violation in the Neutrino sector [1] calls ultimately for three machines in the megawatt range at an energy of 800 MeV. 2) in [2], 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" (OPAL-cycl) 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 3σ . . . 4σ in the PSI Ring cyclotron. We compare our results at the extraction with measurements, obtained with a 1.18 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.

A. A BRIEF LOOK AT OPAL
OPAL (Object Oriented Parallel Accelerator Library) is a tool for charged-particle 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 [3]. OPAL is based on IPPL (Independent Parallel Particle Layer) [4] 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 [5]. More details on cyclotron modeling which are direct relevant to this article can be found in [6]. Several flavors of OPAL are available. For details we refer to the User Manual [3]. In this paper we use OPAL-t for the tracking of 72 MeV beam line, connecting two cyclotrons, the Injector 2 and the Ring Cyclotron. The other OPAL flavor -OPAL-cycl -is designed specially for cyclotron beam dynamics and, is explained in the next section.

B. THE BEAM DYNAMICS MODEL OF OPAL-cycl
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 m 0 , q, γ are rest mass, charge and the relativistic factor. With p = m 0 cγβ we denote the momentum of a particle, c is the speed of light, and β = (β x , β y , β z ) is the normalized velocity vector. In general the time (t) and position (x) dependent electric and magnetic vector fields are written in abbreviated form as B and E.
If p is normalized by m 0 c, Eq. (1) can be written in Cartesian coordinates as The evolution of the beam's distribution function f (x, cβ, t) can be expressed by a collisionless Vlasov equation: where E and B include both external applied fields, and space charge fields 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 quasi-static approximation. In this approach, the relative motion of the particles is non-relativistic in the beam rest frame, so the self-induced magnetic field is practically absent and the electric field can be computed by solving Poisson's equation where φ and ρ are the electrostatic potential and the spatial charge density in the beam rest frame.
The electric field can then be calculated by and back transformed to yield both the electric and the magnetic fields, in the lab frame, required in Eq. (4) 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 [7], and hence it is a good approximation to use open boundary conditions. The combination of Eq. (3) and Eq. (5) constitutes the Vlasov-Poisson 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 right-handed Cartesian coordinate systems are defined, as shown in Fig. 1. The first coordinate system is the fixed laboratory frame S lab , in which the external field data is stored and the particles are tracked.
Its origin is the center of the cyclotron and its X − Y plane is the median plane and the positive direction of Z axis points vertically upwards. The second coordinate system is the local instantaneous frame S local , which is a temporal auxiliary frame for the space charge solver. Its origin O is the mass center of the beam and the orientation of the Y axis is coincident with the average longitudinal direction and the positive orientation of the Z axis points vertically upwards. The third coordinate system is the beam rest frame S beam , which is co-moving with the centroid of the beam. It has the same orientation and origin as S local , but the length in longitudinal direction is scaled by 1/γ due to relativistic effects.
At each time step, the frames S local and S beam are redefined according to the current 6D phase space distribution, and all particles are transformed from S lab to S local , then a Lorentz transfor-mation is performed to transform all particles to S beam . The Poisson equation is then solved in the frame S beam . In a 3D Cartesian frame, the solution of the Poisson equation at point (x, y, z) can be expressed by with G the 3D Green function assuming open boundary conditions. Details of the space charge field calculation can be found in [8].
The model of the external magnetic field is based on mid-plane field measurements with excited trim coils. In consequence we have a vertical field, B z , measured on the median plane (z = 0) as a function of azimuthal position (θ). Since the magnetic field outside the median plane is required to compute trajectories with z = 0, the field needs to be expanded in the Z direction. According to the approach given by Gordon and Taivassalo [9], by using a magnetic potential and the measured B z on the median plane at the point (r, θ, z) in cylindrical polar coordinates, the 3rd order field can be written as where B z ≡ B z (r, θ, 0) and All the partial differential coefficients are computed on the median plane data by interpolation, using Lagrange's 5-point formula. 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 [6]. 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 4 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 OPAL-t simulations from the middle of the last valley before extraction from the

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 90% 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 [2].
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 S-shaped 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 [2].

A. The Flattop Phase
Although a compact beam is observed at the extraction of the Injector 2 cyclotron, the bunch length increases from about σ = 6 mm to about σ = 31 mm 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 S-shape 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 intensity-dependent and, there exists an optimum flattop phase for a given intensity.
For our simulation we use 11.5% of the sum of the main cavity voltages as the flat-top cavity voltage (as set in the control room) and adjust the phase to obtain the same phase difference between main and flat-top cavity as set in the control room.

B. 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 ν r and ν z . This requires a strong excitation of the trim coil TC15.
Without this correction, the coupling resonance ν r = 2ν z 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 [11]. It is described by Eq. (11).
where R is the radius, B a is the maximum magnetic field, the constants A 1 = −1.08, A 2 = 1.08, 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, where R is the orbit radius, B is the hill field, dB dR is the average field gradient in radial direction. 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.

C. 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 where k is the field index. For the PSI Ring cyclotron this gives about 6.0 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 non-centered injection into the PSI Ring cyclotron is used.
Since ν r ≈ 1.7 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 16 mm, hence it allows the extraction of a high intensity beam with very low losses.

D. 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 9 turns of the PSI Ring cyclotron. The probe is located 30 cm upstream from the 50 µ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 1.3 MW cw beam. In order to compare the simulations with measurements, not only a radial probe is implemented in OPAL-cycl 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  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 [12], we can vary a single parameter m from 0 to ∞ and cover a wide range of distributions: from KV to Gaussian. We find that a parabolic distribution (m = 2) matches best the measurement at the crucial point, the septum. On the extreme side, as expected, the Gaussian distribution (truncated at 3σ) 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 R = 4434). The energy spread ∆E sc (linear) = e∆U sc (linear) caused by the linear longitudinal space charge force after n revolutions is given by Joho [13] using the sector model: where g 1c ≈ 1.4 (form factor), Z 0 = 1/4πε 0 c = 30Ω and I is the average current, ∆φ is the phase width, n the turn number and and note that f n depends strongly on the beam distribution and is in our case an open parameter in the range of 0.1 . . . 0.5.
According to Eq. 13, the space charge induced energy spread leads to a radial spread (∆R sc ) that results in an increase in beam size. In Fig. 10 we compare beam sizes at the extraction for beam currents from 10 µA to 2.2 mA with simulations. The relation between the average energy gain dE/dn and the turn number n is: where E f is the final energy and E i is the initial energy. For single turn extraction the loss on the septum is limited by the ratio leading to the condition ∆E sc (nonlinear) < µ n dE dn .
We obtain empirically, for a centered beam a value for µ n which is µ n ≈ 1/3, where as for an eccentric beam, the turn separation is enhanced, and hence µ n can be as high as µ n ≈ 1.
Putting (14), (15), (16) and (17) together, we get for the current limit from longitudinal space charge forces where U f = E f /e and U i = E i /e. Since the turn number n is inverse proportional to the cavity voltage V cav , we see the big advantage of a large cavity voltage  This is remarkably close to the present current limit of 2.3mA (2010), given the crude assumptions: • no turn structure inside the charge sheet (see [13] figure 3), • non relativistic approximation, • no radial boundary condition and • uncertainties in f n and µ n .
We note that for the PSI Injector 2 (18) is not applicable due to phase mixing [14] 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 time-structure 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 3σ . . . 4σ), 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 6-D phase space, and the lack of particle-matter interaction in our model. Particle-Matter 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 particle-matter interaction model are reported in [15] and ideas for secondary electron creation and field emission can be found in [16].