Beam quality and beam loss predictions with space charge for SIS100

The SIS100 synchrotron as a part of the new FAIR accelerator facility at GSI should be operated at the “space charge limit” for light and heavy-ion beams. Losses due to space charge induced resonance crossing should not exceed a few percent during a full cycle. Detailed magnet field measurements are now available for 72 out of the total 108 main SIS100 dipole magnets. Particle tracking studies including nonlinear field errors up to 7th order in the main magnets together with different space charge models are performed. Because of the long time scales reduced space charge models are employed for tune scans. First comparisons with simulations using a self-consistent space charge solver are discussed as well as potential measures to further improve the options in tune space for the reference intensities and beyond.

chosen as a reference case in this study. After acceleration (up to 2.7 GeV/u for U 28+ ) the beam is either extracted slowly over one or more seconds or re-bunched and compressed before fast extraction.
The budget for space charge induced heavy-ion beam loss is determined by the effectiveness of the vacuum system to stabilise the residual gas pressure against beam loss induced desorption. In addition, and for all ions, beam loss induced component damage and activation have to remain below their threshold values. Dedicated collimation systems for stripping and space charge induced 'halo' losses are foreseen in SIS100 [5]. Overall, space charge induced loss should remain below a few percent to allow for additional loss mechanisms due to stripping or at transition energy, for example.
The budget for space charge induced transverse emittance growth in SIS100 is determined by the requirements at the production targets. The in-flight separation of radioactive ions requires a small beam spot size on the target, for example. In the longitudinal plane there is a tight budget for emittance growth because of the short bunch length required for the subsequent storage of the produced secondaries.
At injection energy heavy-ion beams fill approximatelly half of the available SIS100 aperature. Besides direct space charge forces also image forces from the elliptical beam pipe can modify the incoherent tune shifts and should be included in simulation models. In the present contribution we focus on the effect of direct space charge.
Besides the 3D space charge forces the particle tracking simulation have to resolve the local magnet errors, extracted from magnetic field measurements of individual magnets, if available. Because of the demanding time scales (up to 1 second or 10 5 turns) non self-consistent space charge models have been used so far for the 2D tune scans, leading to the identification of suitable working point areas. With the improved performance of our self-consistent Particle-In-Cell (PIC) schemes we can now use a sufficiently large number of macroparticles for noise reduction and validate the identified areas, within accceptable computing times.
A very first goal of our studies is to identify working point areas for low loss operation for the design space charge tune shifts and for the present linear and nonlinear SIS100 magnet error model, based on the ongoing magnet measurements. A second goal is to identify optimisation knobs to be used in the later high-intensity operation, to further reduce the losses and also to determine the actual 'space charge limit' of SIS100 according to our computer model. In particular, we look at the dependence on the beta-beat correction and on the bunch profile.
This contribution is organised as follows: first we present the SIS100 magnet error model, second beam intensities and the resulting space charge tune shifts in SIS100 are discussed, the simulation models and codes are presented together with their underlying space charge models, simulation results are presented and interpreted, finally we conclude with an outlook.

Main magnet field errors
The SIS100 magnet system [6] consists of 108 main dipole magnets, 166 main quadrupole magnets, and various correction magnets. The magnets are presently under production and testing.
The dipole magnets [6,7] are superferric superconducting magnets. The iron yoke is bent, the length is 3 m and the curvature radius is ρ = 52.632 m. A magnet testing facility with a rotating -2 -coil system [8] has been established at GSI Darmstadt. Presently, the field error data for 72 series dipole magnets is available.
The magnetic field representation using the harmonic coefficients B n , A n is where r 0 is the reference radius, B n are the normal components, A n are the skew components. B 1 is the dipole component, B 2 is the quadrupole component, etc. The relative field error components for a dipole magnet are defined as integrated over the magnet length. For a quadrupole magnets, the definition is corresponding, with the leading component B 2 which is used in the denominator. Once the field data for all of the magnets is available, the (b n , a n )-sets for each individual magnet will be used for numerical studies. Presently, we use a model for the SIS100 dipole magnets, based on the data from the available magnets. It is assumed that each error coefficient is described by a "systematic" component µ n , which is the mean of the distribution, and by a "random" component σ n , which is the standard deviation of the Gaussian distribution, The random component is extracted from the available magnets as the sample standard deviation. For the geometrically allowed harmonics (for a dipole magnet B 3 , B 5 , B 7 ), the systematic components are given by the average of the available magnets. For the rest of the harmonics (non-allowed), the systematic components have been assumed zero. The present status of the resulting model for the SIS100 dipole magnets is shown in figure 1, left (r 0 = 30 mm). A unit implies ×10 −4 for a value. The dots correspond to the systematic components, the error bars show ±2σ of the random components.
The model for the SIS100 main quadrupole magnets (figure 1, right) is based on measurements of prototype magnets at JINR, Dubna (internal project reports). The assumption about the systematic components applies to the only allowed component B 6 . The notation is equivalent to the dipole magnets case, the reference radius here is r 0 = 40 mm. The assumptions about the random components correspond to the "worst-case" scenario. It is expected that in the series production the random errors, especially for b 4 , b 6 , can be smaller, but in this study we use these assumptions.

Space charge in SIS100 beams
In this study we focus on a SIS100 reference U 28+ beam at injection energy, because it is characterized by strong space charge tune shift and relatively large beam size, see table 1. The large beam sizes compared to the physical aperture makes beam loss with space charge especially important, as any emittance growth easily results in beam loss. Lighter ions and protons, with similar space charge tune shifts, will have smaller transverse beam sizes due to higher injection energies.   Table 1. Parameters of the nominal heavy-ion bunch in SIS100 during the accumulation. E 0 is the kinetic injection energy, N b the number of ions per bunch. h is the harmonic number, Q s is the synchrotron tune, B f is the bunching factor (for an rf bucket). ε x,y are the unnormalized transverse emittances (4 times the rms emittances). σ z is the rms bunch length, δp/p is the rms momentum spread, δQ ξ is the rms tune spread due to chromaticity. ∆Q sc is the maximum incoherent space charge tune shift for the horizontal/vertical plane.
The space charge tune shifts for the reference heavy-ion scenario in SIS100 are illustrated in figure 2. The parameter ∆Q sc here is the maximal tune shift for the particles with small transverse amplitudes in the longitudinal bunch center, where R = C/(2π) is the effective ring radius, β and γ are the relativistic parameters, r p = -4 -q 2 ion /(4π 0 mc 2 ) is the classical particle radius, λ 0 is the peak line density (at the bunch center), and ε ⊥ is the total transverse emittance. The factor "2" is because of the Gaussian transverse profile. Due to the elliptical transverse cross section the maximal vertical tune shifts are stronger than that in the horizontal plane. For this, the parameter ε ⊥ in eq. (3.1) should be replaced by here for the vertical plane, for the horizontal plane it is corresponding. The parameters ε y , ε x are the vertical, horizontal rms emittances, Q 0y , Q 0x are the machine set (or "bare") tunes. There are bunch manipulations after acceleration, for the generation of a single short bunch, which are not shown in figure 2. The space charge tune shifts are generally much smaller at top energy. Only during a very short (several turns) period the space charge tune shift of the fast rotated bunch reaches up to −0.7 [3]. The bunch compression in SIS100 is not covered in this report.
Particle tracking simulations for beam loss predictions in this report are focused on the reference heavy-ion bunch during the accumulation plateau. The corresponding incoherent tune footprint covers several potential resonance lines, see figure 3. Resonances up to octupole order (b 4 , a 4 ) are shown.
Due to the variation of the transverse space charge force along the bunch, the particles cross the resonances as they perform synchrotron oscillations. The black dot shows an exemplary machine tune (close to the nominal tunes), in the simulation scans we investigate tunes in the quadrant 18.5 < Q x,y < 19.

Particle tracking tools used for SIS100
The tracking code "elegant" [9,10] is employed to model the time evolution of a nominal U 28+ bunch in SIS100 during the accumulation plateau at the injection energy. A 3D Gaussian distribution (cut -5 -2020 JINST 15 P07020 at 2.5 σ) is assumed, the corresponding frozen space charge force model is available in the "elegant" code [11]. 6D symplectic particle tracking is performed in the complete SIS100 lattice, with the realistic magnet geometries and with the systematic/random magnet field errors (figure 1) taken into account. For longitudinal dynamics, SIS100 cavities with the nonlinearities in the RF bucket are included. We complement the "elegant" code results with simulations using two further tracking codes, viz. the "MAD-X" code and the "SixTrackLib" code. The symplectic thin-lens tracking module of MAD-X [12] with the frozen space charge module [13] provides reference results here. The recent SixTrackLib engine is a flexible and portable reimplementation of the "SixTrack" physics in an architecture independent way [14,15]. SixTrackLib simulations can be run on (multi-core) CPUs and GPUs (graphics processing units) and thus profit from tremendous speed-ups when symplectically tracking 6D phase space. The option to track full machine lattices such as SIS100 at up to O(1000) times faster on the latest high-end GPU models (such as NVIDIA V100) than on a single CPU core render SixTrackLib an attractive alternative to the aforementioned codes. On the other hand, the flexible interface of SixTrackLib enables the integration with other codes. In the framework of this study, the GPU accelerated space charge modules of PyHEADTAIL [16,17] are added via space charge kicks into the non-linear lattice tracking of SixTrackLib.
For this study, the above described non-linear magnet field error model of SIS100 has been implemented in SixTrackLib on top of the full lattice.
The single-particle beam dynamics in SixTrackLib are confirmed to be identical to MAD-X: tracking the nominal SIS100 lattice over 20'000 turns, SixTrackLib deviates from MAD-X within O(10 −8 ) relative error. With RF cavities switched off, the trajectory deviation remains even at a -6 -constant O(10 −12 ) relative error. This indicates a minor systematic difference in the RF cavity tracking between both codes, while tracking only the dipole and quadrupole magnets with exact drifts in between yields identical results up to numerical precision. When including the non-linear field imperfections up to 6th order with RF on, trajectories are equivalent within a relative error of 10%. Example trajectories from both codes are shown in figure 4 starting from the same initial condition and tracking single particle trajectories for 20'000 turns. When including the SIS100 aperture model covering all the magnets as well as the septum, tracked particles are found to be lost at identical elements in the same turn when comparing MAD-X and SixTrackLib simulations. Correspondingly, beam loss studies for the SIS100 scenario find identical results for MAD-X and SixTrackLib without space charge.

Space charge descriptions for particle tracking
Direct space charge as the effect of the beam self-fields on the macro-particles can be modelled in a variety of ways. In the course of this study, we investigate and compare fast approximative space charge models with the computationally more demanding self-consistent particle-in-cell (PIC) model. The most realistic simulation of space charge reconstructs the beam self-fields from the macroparticle distribution itself. The standard approach to this self-consistent modelling is to employ PIC algorithms. The PyHEADTAIL 2.5D PIC model used in the present study bins the macro-particles into longitudinal slices, discretises the macro-particles onto a transverse regular grid for each slice and then solves the transverse Poisson equation on these grids based on open boundary conditions with the FFT algorithm using Hockney's cyclic expansion. The thus obtained beam fields on the grids are interpolated back to the macro-particles. Such PIC simulations are able to self-consistently take into account local field changes due to deformations of the particle distribution e.g. under the influence of resonances. Realistic simulations with PIC for 3D bunches typically require a very large number of macro-particles on the order of O(10 6 ) in order to minimise numerical artifacts from grid heating [18], while properly resolving the beam self-fields with a sufficient number of -7 -grid cells (in this study we consider 64 slices with 128 × 128 transverse grid cells). Therefore, PIC simulations are generally known to be computationally expensive.
In order to obtain much faster space charge trackers one often imposes assumptions on the shape of the space charge fields, while the simulated results should ideally resemble the self-consistent results. Such models are referred to as frozen space charge. For SIS100 one assumes 6D Gaussian distributed bunches, which in the particular case of U 28+ fill a considerable fraction of the aperture. Significant changes in the RMS beam sizes (induced e.g. by resonances) would then immediately be followed by beam loss. The assumption of approximating the beam shape to remain Gaussian throughout the entirety of the simulation should be valid for the low-loss tunes of interest. On the other hand, working points leading to non-negligible beam loss would render the absolute beam loss figure non-realistic. However, the present beam loss study aims to identify and delimit areas in tune space with negligible beam loss from the lossy working point regions.
For long bunches, the transverse fields of the Gaussian distribution can be computed with the Bassetti-Erskine formula [19]. In the longitudinal plane one then readily applies a Gaussian line charge density profile. Three typical variants of this frozen Gaussian space charge model may be qualitatively distinguished: 1. fixed frozen: the above Gaussian field map is computed before the simulation with the initial RMS values for the transverse and longitudinal planes and it is centred on the computed orbit. During the simulation, these parameters remain constant and the particles of the tracked distribution are treated as test particles probing this distribution.
2. adaptive frozen: during the simulation, the RMS values for the three planes in the field map are regularly updated from the tracked macro-particle distribution. While the shape of the distribution remains frozen, the resonance conditions for particles within the implied tune spread can move with the RMS beam size changes.
3. matched frozen: in addition to the adaptive frozen model, also the centroid of the tracked macro-particle distribution is subtracted before the forces are computed in order to centre the field map on the bunch. The space charge model in the "elegant" code is of the matched frozen type. In the present simulation set-up, the RMS values are updated every turn at every space charge kick. MAD-X offers space charge in the fixed frozen and the adaptive frozen variants. SixTrackLib by itself supports the fixed frozen space charge model. In combination with PyHEADTAIL it can be extended to all three frozen models besides the aforementioned self-consistent PIC setup.
To suppress noise and inherently affiliated emittance growth from the macro-particle statistics computations, the adaptive and matched frozen variants require more macro-particles (around 5'000 to 10'000) than the fixed frozen space charge model. We found the latter to yield representative beam loss results with 1'000 macro-particles. Although the evaluation of the Bassetti-Erskine formula with the computationally costly Faddeeva or complex error function takes an order of magnitude more time per macro-particle than PIC, the required lower numbers of macro-particles significantly reduce the simulation run time. As a practical result, we found O(min) for 20'000 turns with the fixed frozen space charge model on an NVIDIA V100 GPU with SixTrackLib, while the equivalent PIC scenario takes O(day).

Space charge model comparison
In order to compare the discussed three frozen space charge variants with PIC we set up a simple and exactly reproducible case based on the linear SIS100 lattice without random field imperfections using SixTrackLib with PyHEADTAIL space charge kicks. At a resonance-free working point (Q x = 18.95, Q y = 18.93), the tune foot prints of all four space charge models are confirmed to cover the expected ∆Q sc y = −0.3 (within 3%). For the comparison, the horizontal tune is fixed at Q x = 18.86 far enough from any horizontal resonances. The vertical tune Q y approaches the half-integer resonance Q y = 18.5 from above. Two of the quadrupole magnets ("1" for S52QD11 and "2" for S52QD12) in the extraction region are powered at slightly different integral field values compared to the focusing ("f") and defocusing ("d") quadrupole families: These two quadrupole magnets, further referred to as the "warm magnets", act as a local source for beta-beat increasing the half-integer stop-band width. During normal operation the resulting beta-beat would be reduced by local correctors.
The simulations track macro-particles over 200 turns, apertures are not considered here in order to observe only RMS emittance growth. The matched frozen space charge model is chosen to update the RMS beam sizes of the field map at every kick to follow the evolution of the tracked macro-particle distribution closely. The adaptive frozen space charge model is chosen to be updated only "slowly", i.e. every 5 turns at each kick with a cumulatively averaged RMS beam size. This effectively low-pass filters envelope oscillations.
All four space charge models entail constant horizontal emittance for all considered tunes (the relative difference between end and initial values remain below at most 3 parts in 1000). The increasing vertical emittance growth towards the half-integer resonance is illustrated in figure 5. The PIC results have been exactly reproduced when increasing the macro-particle number by factor four. The matched frozen results have been confirmed at a factor hundred increased number of macro-particles. The adaptive and fixed frozen simulations use the same resolution as the matched frozen setup.
The fixed frozen space charge model practically always implies the least emittance growth. Due to the absent update of the field map, the resonance condition does not change and thus no new particles can be fed into the resonance after depletion at the resonant amplitudes. The PIC results are expected to represent the most realistic results. Below a tune of Q y ≤ 18.62, they indicate the largest vertical emittance growth of all models -demonstrating the well-known mechanism of the beam escaping the resonance by weakening of space charge through increasing (vertical) RMS size. Above Q y > 18.63, the matched frozen and the adaptive model feature the largest relative emittance increase. They exceed the PIC induced emittance growth figures and thus lead to overly weakened space charge. Figure 6 shows the evolution of the vertical emittance during the 200 turns for two chosen working points in both conditions, fairly strong half-integer influence at Q y = 18.6 and very weak influence at Q y = 18.7. While matched and fixed frozen space charge develop similarly at Q y = 18.6, at Q y = 18.7 the fixed frozen space charge shows no significant influence of the -9 -2020 JINST 15 P07020   half-integer resonance any more as opposed to all the other space charge models. In contrast, the matched and adaptive frozen model show an equivalent behaviour with a much faster initial rise compared to the realistic PIC case and then they increase the vertical emittance at a slightly higher rate than PIC. At this point it is interesting to note that the overly increased emittance growth for the matched frozen model at Q y = 18.7 is not due to lacking resolution and corresponding noise effects -much on the contrary, the emittance growth curve remains exactly the same for 10 and even 100 times more macro-particles. This rules out noise effects and rather points to the enforced Gaussian shape of the field map as a reason for the "overshoot", which is an inherent feature of the updating Gaussian field map models (no update as in the fixed frozen case prevents such reflexive dynamics by construction). In the PIC case, on the other hand, the particles are free to redistribute locally to counteract and escape the resonance yielding a transverse profile deviating from a Gaussian.
In this limited setup with a strong and fast acting half-integer resonance we can already identify different dynamics implied by the various space charge models. Other longer-term effects such as incoherent resonance crossing and trapping, which depend on the synchrotron motion, are excluded -10 -here due to the relatively short simulation time below a synchrotron period. Nevertheless, such incoherent effects are known to be well captured also by the frozen models in comparison to self-consistent PIC, as many contributions in literature such as [20] prove.
In conclusion, we differentiate between both areas: under strong resonance influence, the matched frozen and the fixed frozen model give similar results which are too optimistic in terms of the RMS emittance growth in comparison to self-consistent PIC. Under weak resonance influence, fixed frozen space charge leads to underestimation of the effect while the matched frozen space charge (and adaptive model) tend to overestimate the resonance impact.

Simulation results with the full nonlinear SIS100 lattice
At present we rely on results obtained with the matched frozen space charge model in "elegant", keeping in mind that differences with PIC still require more dedicated studies. Results of twodimensional tune scans are presented in figures 7, 8. The color scale shows the beam loss after 20000 turns, i.e. 0.1 means that 90% particles remained in the bunch. The results for the reference U 28+ bunch (figure 7, right) indicate a "good" tune area visible by the violet color. Long simulations for 160000 turns (corresponding to the accumulation time in SIS100) for the tunes from this area showed only very slow further losses. The overall beam losses below 5% would be then possible, which fulfills the low-loss requirement. The expected closed-orbit distortions will be corrected by the SIS100 steerers. For the simulations, we assume a residual closed orbit distortion of ∆x rms = 2 mm, ∆y rms = 1 mm. 18  A comparison to the simulation scan performed without space charge (figure 7) demonstrates the important role of space charge for beam loss in SIS100.
The expected variation of the integrated gradient in the quadrupole magnets corresponds to an rms beta beat of 5%. In order to analyse the effect of the beta beat on the losses, simulations with a lower variation of the quadrupole gradient are performed. figure 8,  of the losses. Simulations for a lattice without nonlinear field errors (for n > 2, all b n , a n are zero, figure 8, right) show a broad stopband close to the half-integer resonance and along the coupling line. Space charge induced nonlinear error resonances, caused by the beta beat and nonlinear space charge, are not visible. Presently we observe such resonances in other simulation models and this is still a topic of ongoing investigations. Although this case might be academic and represents an ideal magnetic system with only linear errors, it is interesting for space charge model comparisons and to idenfify possible space charge resonances in the full scenario. Finally, it is important to point out that within the identified low loss tune area also the emittance growth is very small, the beam profiles do not distort, and the beam quality remains within the requirement.

Mitigation measures
For the present SIS100 magnet error model and for the heavy-ion reference bunches we identified a low loss area in tune space, potentially suitable for high-intensity operation. However, more safety margin and further flexibility in tune space would be useful, also if one wants to operate at even higher intensities. The following measures are under investigation and could allow for lower losses and enlarged low loss areas in tune space: 1. quadrupole field quality: our present model for the nonlinear field components of the SIS100 quadrupole magnets is rather conservative and the actual magnets might show a better field quality, especially regarding the random component. As soon as the field measurements of the first quadrupole units will be available we will update our models and support the magnet quality control during the series production..
2. beta beat correction: the further correction of the beta beat should reduce the half-integer error stopbands and also space charge induced resonances.
-12 -3. dual rf buckets: during accumulation the bunch profiles can be flattend using dual harmonic buckets. First simulation results show an improvement of the losses, corresponding to the higher bunching factor.
4. pulsed electron lenses: for the existing SIS18 pulsed electron lenses for partial compensation of the space charge tune spread are under consideration [21]. Studies for the SIS100 are ongoing.

Conclusions and outlook
Tracking studies using the present SIS100 magnet error model together with different space charge models were performed. Two-dimensional tune scan were performed with the code "elegant" with a frozen space charge model. A low loss area potentially suitable for high-intensity operation has been identified.
One of the applications of this study is the magnet field quality assessment for the SIS100 magnets during production. The magnetic field data from 72 (out of total 108) series dipole magnets has been used in our tracking simulation. For the main quadrupole magnets, a very conservative model based on the first prototype magnet has been applied. Conservative assumptions for other imperfections (beta beat, closed-orbit distortions, etc.) have been applied. The present beam loss predictions allow for the conclusion that the SIS100 dipole magnets show the high field quality required for high-intensity, low loss operation.
First comparisons with the newly developed tracking suite SixtrackLib/PyHEADTAIL and self-consistent space charge solvers indicated some discrepancies, to be understood more in detail in the future. Ongoing work also focuses on different optimisation measures towards a further understanding and increase of the "space charge limit" in SIS100.