High intensity space charge effects on slip stacked beam in the Fermilab Recycler

In order to achieve high intensity beams for Fermilab ’ s neutrino program, slip stacking is used in the Recycler in which bunches overlap and slip through other bunches multiple times. As the bunch intensity is increased in the future, space charge effects during the overlap periods become more detrimental. In order to investigate the size of these space charge effects, the beam simulation program S YNERGIA has been extended to allow copropagation of bunch trains at different momenta.


I. INTRODUCTION
The Fermilab accelerator complex has recently been reconfigured and upgraded to meet the needs of the U.S. high energy physics research community as defined by the Particle Physics Project Prioritization Panel (P5) report of 2014 [1]. Under the proton improvement plan (PIP), several steps were taken to increase the delivered intensity from the complex, achieving 700 kW on target to the neutrino program.
Proton improvement plan-II (PIP-II) [2] is Fermilab's plan for further increasing the beam power to the laboratory's experiments. The increased beam power will position Fermilab as the leading laboratory in the world for accelerator-based neutrino experiments. PIP-II will see the addition of an 800-MeV superconducting linear accelerator called SCL, which capitalizes on the lab's expertise in superconducting radio-frequency cavity technologies. Along with improvements to Fermilab's existing Main Injector and Recycler accelerators, PIP-II will aim to provide a 1.2 MW proton beam that is needed for the DUNE experiment [3][4][5][6]. The key parameter differences between PIP and PIP-II are shown in Table I. For PIP-II, the Recycler will be required to accumulate 50% more beam than it does under current operating conditions. With increasing intensity come increasing space charge effects which may cause emittance growth. Of particular concern is the effect on our intensity-doubling process known as slip stacking [7,8] during which bunches partially overlap. Slip stacking is used in the Recycler to double the intensity of the bunches sent to the Main Injector. The Recycler has space for seven booster batches, however, six batches are injected operationally to allow a gap for the kickers to fire. Slip stacking is used to allow the injection of 12 batches which are stacked and combined to form six double intensity batches. A schematic of slip stacking is shown in Fig. 1. Initially, six batches of 81 bunches known as the A group are injected into the Recycler at its central momentum (a) and then decelerated by 1260 Hz (b). This frequency separation is determined by the repetition rate of the booster as the batches must slip the equivalent of a single batch between injections such that the injection gap is always clear. A second rf cavity is then turned on at the original frequency and a further six batches known as the B group are injected (c). Due to the lower momenta of the A group compared to the B group, the A group will slip with respect to the B group. Near the end of the cycle, both groups are accelerated to optimize the aperture by 630 Hz. When the A and B groups line up (d), the double intensity batches are extracted to the Main Injector where they are accelerated to 120 GeV.
The purpose of this paper is to investigate the size of these space charge effects at intensities for PIP-II and beyond and to determine if they will be detrimental to beam operation. In order to do this, the collective beam effects code SYNERGIA [9,10] was upgraded to deal with two copropagating beams.

II. SYNERGIA
SYNERGIA is a particle-in-cell (PIC) based framework for self-consistent, high fidelity modeling of charged particle beam transport in accelerators or storage rings in the presence of collective effects such as space charge and wakefields. SYNERGIA tracks macroparticles contained within one or more beam bunches through the accelerator. For purposes of calculating collective effects, each macroparticle represents the charge of 100-10 000 real beam particles, but responds to fields as a single particle. One or more bunches of macroparticles are tracked through a set of common magnetic elements and rf cavities. Detailed aperture descriptions are associated with each beam line element for simulating the source of beam losses. Apertures may be rectangular or elliptical masks of arbitrary size and transverse offset or an arbitrary polygon. Macroparticles which strike an aperture are removed from the simulation.
Macroparticles are organized within one or more beam bunches which exist within one or two bunch trains. Bunches within a train have a fixed separation. However, bunch trains may have different momenta with the consequence that they slip longitudinally relative to each other. The relative position of the bunches is taken into account in the space charge solver. The bunch initial particle distribution may be matched to the machine using nonlinear matching techniques or specified by an arbitrary input distribution.
The particle time coordinates are maintained relative to a master rf frequency which is set to be h=T where h is the harmonic number and T is the orbit time of the reference particle. SYNERGIA also allows additional rf cavities that may have frequencies displaced from the master frequency. In this case, the phase of those cavities slips relative to the master clock during each revolution of the bunch train by hTΔω, where Δω is the angular frequency difference. Bunch trains traveling at the momentum of the reference particle are synchronized and have bunch stability with the master frequency. The bunch average does not receive any additional longitudinal kick from the master rf cavities. Any kicks that they might receive from off-frequency cavities average to 0, and thus those bunches remain synchronized with the reference particle. Bunch trains traveling at a displaced momentum synchronized with the off-frequency rf cavities similarly maintain longitudinal bunch stability at their frequency, but their longitudinal position slips relative to the primary trains synchronized with the master clock.
Collective effects are combined with particle tracking using the split-operator method [11]. For a Hamiltonian separable into single-particle and collective pieces, The time evolution mapping for a time step τ can be approximated using the split operator method: SYNERGIA provides a variety of space charge solvers. The simulations described in this paper use a 2.5D open boundary solver. Participating particles are binned in a twodimensional rectangular grid to obtain the collapsed charge density distribution. The Poisson equation is solved with fast Fourier transforms as described in Hockney and Eastwood [12] to obtain the potentials which are interpolated at each macroparticle position to determine the transverse Electromagnetic field. The field at each particle and the size of the momentum kick is modulated by the longitudinal charge density at the particle's location. In the case of two bunch trains, the space charge solver has to superimpose the fields generated by overlapping bunches. Figure 2 shows an example of the fields used in the space charge solver for a bunch in the high momentum train.
Within a train, bunches exist inside of rf buckets determined by the rf frequency. At the point in the lattice where space charge is applied, the buckets of both bunch trains are matched according to their timing offset to identify overlapping buckets. Because opposite bunch trains may be traveling at a different velocity, each bucket may overlap with two buckets in the opposite train. For each bucket, the solver calculates the potential produced by the particles contained within. That potential produces a kick that is applied to particles within that bucket, and also to particles contained within overlapping buckets. The procedure is repeated for each bunch train, resulting in the application of the space charge kick to both bunch trains. The phase space coordinates of all macroparticles are saved each turn during the simulation for later analysis. Macroparticles that strike apertures during the simulation are removed and their location at the time of loss is saved.

III. RECYCLER SIMULATIONS
We performed simulations to investigate the effect of multiple bunch overlaps during the slip stacking process under operating conditions. The first condition investigated is referred to as the base condition. Additional simulations changed the initial tunes, initial longitudinal distribution, transverse emittance, and chromaticity. Further simulations look at how the PIP-II bunch structure modifies the emittance evolution. A schematic of the evolution of the simulation is shown in Fig. 3. For computational tractability, the two bunch trains have eight bunches each. The lower momentum train is offset in the negative Z direction such that at the beginning of the simulation, the eighth lower momentum bunch (LMB 8) is one bucket away from the first bunch in the higher momentum train (HMB 1). Due to the momentum offset, the lower momentum bunches slip with respect to the higher momentum bunches. At the end of the simulation, HMB 1 will be overlapped with LMB 1 and so on for each subsequent bunch.
The positions and strengths of the Recycler magnetic elements are specified by a MAD8 [13] file. The optical functions produced from this file by MAD8 and SYNERGIA are consistent. The simulation input contains measured machine imperfections up to fifth order. Each element and drift segment has an associated aperture with the correct dimensions and shape to enable realistic modeling of beam loss magnitude and location. The simulation includes the correct beam steering through the Lambertson magnets [14] with their aperture restrictions.

A. Base
The initial case to be investigated is referred to as the base case. It is set up using the parameters shown in Table II. Input transverse and longitudinal bunch distributions for simulations were prepared using the following procedure. The transverse distribution is a matched Gaussian with normalized 95% emittance of 15π mm mrad which is consistent with measurement observations. The longitudinal phase space is populated with a measured tomography distribution which has been allowed to undergo filamentation for 4600 turns in a 80 kV bucket. The resulting distribution is used as the input for higher momentum bunches is shown in Fig. 4 (left). A measured distribution was not used directly as many of the overlaps occur with bunches that have been in the machine for a significant number of turns as well as to avoid confusing any effects from the initial longitudinal motion. To determine the distribution for lower momentum bunches, the higher momentum bunch distribution is decelerated for a number of turns equivalent to 50 ms until the rf frequency reaches the desired separation of 1260 Hz from the original frequency.  This distribution is shown in Fig. 4 (right). Particles that fall outside of the bucket during filamentation and deceleration are removed.
With the discussed input distributions, the simulations are performed as described in the previous section. During the run, the six-dimensional phase space of all the macroparticles is saved at the same location every turn, a sample of which is used to calculate the tune footprint. The footprint at the start and middle of the simulation is shown in Fig. 5. The top plot shows the combined tune footprint for HMB 1 and LMB 1 at the beginning of the simulation. The histograms show a one-dimensional projection of the footprint for each bunch with the higher momentum bunch in blue and the lower momentum bunch in green. The bottom plot shows the same, but at turn 282, which is where HMB 1 is overlapped with the LMB 5 bunch. Focusing on the tune footprint, it is apparent that the footprint for HMB 1 has increased due to its interaction with LMB 5, whereas the footprint for LMB 1 is unchanged as it is yet to overlap any higher momentum bunches.
To show this effect further, Fig. 6 shows how the mean of the tune distribution in the horizontal plane changes during the simulation. For the HMB 1 and LMB 8, which both undergo eight overlaps, the mean shifts more than 0.02 during an overlap. For LMB 1 and HMB 8, there is no increase until the end when these two bunches finally overlap with each other. The incoherent tune shift is dependent on the bunching factor. Therefore, as the bunches oscillated in the bucket, this longitudinal motion will result in a modulation of the tune shift. The different frequencies that can be seen are related to the 1260 Hz separation and the synchrotron frequency.
The momentum offset between bunch trains needed for slip stacking is determined by the booster cycle frequency. The Recycler holds seven booster batches. For slip stacking to work, the off-momentum booster batch should slip 1=7 of the Recycler circumference during one booster cycle. Bunch trains with fractional momentum offset Δδ are synchronous with rf cavities with a frequency offset: where f rev is the revolution frequency, h is the harmonic number of the Recycler and η is the slip factor. For PIP, the momentum separation is −0.0027, while it is −0.0037 for PIP-II to meet the larger slipping rate required for the booster running at 20 Hz. The nonzero chromaticity ξ in the Recycler acting on the momentum separation results in a tune shift ΔQ given by The on-momentum beam and off-momentum beam will be separated by ΔQ for a given set machine tune, restricting the available tune space. The predicted shift, which to first order is simply given by the chromaticity During the course of slip stacking, any particular bunch will temporarily overlap with a succession of bunches in the opposite bunch train during which the peak of which the local density will be double the nonoverlapped density. In order to understand the effect of these overlaps, the evolution of the 4D emittance (in the absence of coupling, the 4D emittance is the product of the horizontal and vertical emittance) is shown for each bunch in Fig. 7. As the bunch overlaps another bunch, a step in emittance growth is observed. HMB 1 undergoes the most overlaps; after overlapping eight bunches, the 4D emittance is 15% larger. For HMB 8, which does not overlap at all, the final 4D emittance growth is just ∼4%. The cause of the growth is believed to come from resonant stopbands. The lower momentum bunches follow a similar pattern, except that the size of the growth is smaller for these bunches. This can be explained due to the chromaticity-induced tune shift discussed in the paragraph above. Thus, if the 4D emittance growth is caused by hitting a resonance line, the positive tune shift caused by chromaticity on the lower momentum bunch will move it further away from the stopband.

B. Tune effect
The expected cause of the 4D emittance growth is due to the space charge tune shift causing particles to cross resonance lines. In this case, the most likely lines would be third order lines. To investigate this, the simulations were then repeated with different set tunes. The horizontal tune was lowered and increased by 0.02 and the vertical was decreased by 0.02 and increased by 0.03. As expected, lowering the tune results in a larger 4D emittance growth as more particles will cross the third order resonance line. Figure 8 shows the results of these scans on HMB 1 and LMB 8, the bunches that undergo the most overlaps.

C. Phase offset
Operationally, the longitudinal emittance can be diluted and the bunch length effectively increased by introducing an injection phase offset. To simulate this, the initial tomography distribution was offset by 1.5 ns and the filamentation procedure was repeated to produce the input distributions that are shown in Fig. 9. The bunch length for this distribution is 20% bigger than without the phase offset.
The 4D emittance growth seen in Fig. 10 shows that adding a phase offset appears to cut the 4D emittance growth in half. This is expected as the bigger bunch length from the phase offset will lead to a smaller space charge tune shift.  Lowering the tune leads to increased 4D emittance growth and raising the tune leads to a decrease in 4D emittance growth which is consistent with the source of the growth coming from third order resonance lines.

D. High chromaticity
High chromaticity operation was also investigated by setting both horizontal and vertical chromaticity to −18. The Recycler operates below transition energy and therefore is operated at negative chromaticity for stability, thus, high chromaticity in this case means high absolute chromaticity. The 4D emittance growth for the high chromaticity cases are shown in Fig. 11. For HMB 1, there was a small increase in the 4D emittance in both cases which is most likely due to the increased tune spread from chromaticity. The lower momentum bunches see a different effect from chromaticity due to their momentum offset. As discussed earlier, the larger chromaticity causes the footprint to shift to higher tunes. The result is a smaller 4D emittance growth when the vertical chromaticity is increased as the footprint is pushed away from the third order resonance line. However, when the horizontal chromaticity is increased, the 4D emittance increases. This is because the horizontal set tune is much higher, so raising the chromaticity pushes the footprint closer to the half integer lines resulting in a larger 4D emittance growth. Thus, to optimize the tune space used, the set tune cannot simply be set as high as possible to avoid the third integer lines but must also be low enough to avoid the lower momentum bunches touching the half integer resonances.

E. Transverse emittance
In order to see the effect of increasing the transverse emittance, the simulations were repeated where the input transverse distribution was matched to 18π mm mrad and 20π mm mrad. Figure 12 shows that increasing the transverse emittance results in less overall 4D emittance growth. This is to be expected as space charge is inversely proportional to the transverse emittance.

F. Base case simulation conclusions
A series of simulations were performed studying space charge effects when two trains of bunches overlap. A base case was investigated looking at its sensitivity to different parameters such as tune, chromaticity, phase offset and transverse emittance. The dominant effect observed in these simulations was emittance growth when two bunches overlap, doubling the space charge tune shift causing some particles to cross resonant lines. The effect can be reduced  FIG. 11. The effect of high chromaticity on the 4D emittance growth. High chromaticity has a smaller negative effect on the higher momentum bunches. The larger tune offset caused by chromaticity can move away from the bunch from resonant lines at lower tune but potential push into higher tune lines such as the half-integer resonance.
FIG. 12. The 4D emittance evolution for when the input distributions have a different starting transverse emittance. The larger the initial transverse emittance, the smaller the 4D emittance growth. by using a higher set tune and implementing a phase injection offset which results in a larger bunch length. Running with a set tune of 25.46 horizontally and 24.42 vertically saw the 4D emittance growth reduced to 3% for the high momentum bunches.

G. PIP-II
The cases of PIP and PIP-II are compared with the base case. The PIP case is the same as the base case except the intensity is a factor of 2 smaller. PIP-II has the same intensity as the base case, however, the input distributions are filamented in 140 kV buckets and the lower momentum bunches are decelerated to a larger separation of 1680 Hz as the booster repetition rate is increased to 20 Hz. The PIP-II simulations are also run for fewer turns as the larger momentum offset means the time between overlaps is ∼53 turns rather than ∼71. The 4D emittance growth for the three cases is shown in Fig. 13.
The 4D emittance growth for PIP is considerably smaller compared to the other cases. For the PIP-II case, the 4D emittance growth is a little bit larger. This is most probably caused by an increased space charge effect due to the 140 kV bucket decreasing the bunch length resulting in higher local charge densities. Furthermore, for PIP-II, the sensitivity to tunes and an injection phase offset was found to be similar as for the base case. Figure 14 shows the base case run with higher tunes (Q h ¼ 25.46, Q v ¼ 24.42) and PIP-II as the same high tunes. Both show significant improvement by raising the working point. Table III shows the bunch length and momentum spread for each input for the base case and PIP-II cases. The amount of uncaptured beam was also calculated as the number of particles outside of the eight bucket windows at the end of the simulation. For the base case, adding a phase offset roughly doubles the amount of beam that is lost from the buckets. Moving to the larger frequency separation for PIP-II reduces the amount of beam lost by around a factor of 3. Adding a phase offset for PIP-II does not seem to result in more uncaptured beam.

IV. RESONANCE SCANS AND COMPENSATION
It is believed that the primary cause of the 4D emittance growth is due to the large space charge tune shifts causing particles to enter resonances. While it has been shown that raising the working point can help reduce this effect, operationally, this is not always the case so the possibility to run at lower tunes would be useful. Tune scans were performed in both planes without space charge to identify the strongest resonances driven by just the lattice. While the tune in one plane was varied, the tune in the other plane was left at its base value. At each tune point, a bunch of 65 536 particles was tracked for 565 turns. In the horizontal plane, the two most prominent resonances are the 3Q x ¼ 76 and the 2Q x þ Q y ¼ 75. Small 4D emittance growth can also be observed at Q x ¼ 0.39 due to the linear coupling resonance. In the vertical plane, the two most prominent resonances are the 3Q y ¼ 73 and the 2Q y þ Q x ¼ 74. Given these results, we attempted to compensate the 3Q x ¼ 76 resonance by adjusting sextupoles to eliminate the third order resonance driving force. The generating term corresponding to the 3Q x ¼ 76 resonance is given by [15,16]   where h 3000 is the Hamiltonian coefficient which depends on the MAD8 field strength K 2 , l is the length of the magnet, β x is the beta function and Δϕ x is the phase advance. In order to compensate this resonance, sextupole magnets can be used to make the numerator of (3) zero. In the case of the simulation, the lattice is defined so the calculation of f 3000 is simple. As f 3000 is complex, two sextupoles with compensation currents I 1 and I 2 are required to achieve full compensation. When the two sextupoles are orthogonal, one sextupole can compensate the real part and one can compensate the imaginary part: where g is the gradient of the sextupoles per current and Bρ is the magnetic rigidity. Powering two sextupoles to compensate a resonance will also result in a chromaticity change. Therefore, two more sextupoles are needed to account for this which are placed at 3Δϕ x ¼ 180°from the other sextupoles and ran at the same currents but with opposite polarities. The linear system above needs to be modified for four magnets and two more constraints are added such that the current of the sextupoles pairs sum to zero. Applying this compensation scheme to the simulations, the tune scans were repeated with the results shown in Fig. 15. A large improvement in transmission is observed at the 3Q x ¼ 76 stopband. It is believed that higher order terms which are not compensated by the first order calculation of f 3000 become dominant close to the resonance. The same scheme was then applied to the slip stacking simulations presented earlier for the base case and the results are shown in Fig. 16. The compensation scheme shows positive effects, reducing the final 4D emittance growth by more than a factor of 2. When the resonance is compensated, saturation is observed within eight overlaps.

V. MEASUREMENTS IN THE RECYCLER
Dynamic tune scans [17] were performed in the Recycler to identify resonances in order to compare with the simulations. A two-dimensional tune space scan was performed in the range 0.32-0.46 horizontally and vertically. The numerical derivative of beam transmission as a function of tune allows the identification of resonances in the ring. This was then repeated for different vertical tune settings and then repeated ramping the vertical tune at fixed horizontal tunes. This allowed a two-dimensional scan of the resonances in the ring which is shown in Fig. 17. The 3Q x ¼ 76 resonance was found to be the strongest resonance just as found in simulations.
Compensation of this resonance was then attempted using sextupoles already installed that are used to control chromaticity during current operations. There are 24 sextupoles at horizontal locations and 36 at vertical locations installed for chromaticity correction where the sextupoles are powered in pairs so a single sextupole cannot be powered independently. Thus, each pair was considered to be a single sextupole for the point of compensating the resonance. To compensate the 3Q x ¼ 76 resonance, two pairs of sextupoles are required which satisfy 3Δϕ x ¼ 90°and another two pairs are needed which satisfy 3Δϕ x ¼ 180°from the other two pairs to compensate any chromaticity change resulting from the other sextupole pairs. Four pairs of sextupoles were identified, s202/204, s128/130 and s116/118, s412/414. In order to test compensation, all the sextupoles were set to zero current and the tune was ramped from 0.42 to 0.32 to cross the third order resonance. This was then repeated with different currents in the sextupoles mentioned above. Compensation of the resonance was observed as shown in Fig. 18. Transmission was improved from 3% to 90%.

VI. SUMMARY
A number of simulations were performed focusing on the space charge effects during slip stacking at high intensity. At PIP-II intensities and beyond, the large space charge tune shifts caused when bunches overlap result in significant 4D emittance growth as the individual particles tunes cross the third order resonance. A compensation scheme was developed and tested in simulations and shown to reduce the 4D emittance growth. A demonstration of this technique was then applied to the actual machine and improved transmission while crossing the third order resonance from 3% to 90%. A dedicated scheme is now in the process of being developed.
The simulations presented in this paper were limited to eight bunches per train which is much less the 81 bunches used for operation. However, in the case of a suitable resonance compensation scheme, after some initial emittance growth, the space charge effect is no longer strong enough to deteriorate the emittance further. Given that this effect was observed within eight overlaps, simulations with more bunches are not required and PIP-II is expected to be feasible in terms of space charge if the 3Q x ¼ 76 resonance is suitably compensated.