Fully 3D Multiple Beam Dynamics Processes Simulation for the Tevatron

We present validation and results from a simulation of the Fermilab Tevatron including multiple beam dynamics effects. The essential features of the simulation include a fully 3D strong-strong beam-beam particle-in-cell Poisson solver, interactions among multiple bunches and both head-on and long-range beam-beam collisions, coupled linear optics and helical trajectory consistent with beam orbit measurements, chromaticity and resistive wall impedance. We validate individual physical processes against measured data where possible, and analytic calculations elsewhere. Finally, we present simulations of the effects of increasing beam intensity with single and multiple bunches, and study the combined effect of long-range beam-beam interactions and transverse impedance. The results of the simulations were successfully used in Tevatron operations to support a change of chromaticity during the transition to collider mode optics, leading to a factor of two decrease in proton losses, and thus improved reliability of collider operations.


I. MOTIVATION
The Fermilab Tevatron [1] is a p-p collider operating at a center-of-mass energy of 1.96 TeV and peak luminosity reaching 3.53×10 32 cm −2 s −1 . The colliding beams consist of 36 bunches moving in a common vacuum pipe. For high-energy physics operations, the beams collide head-on at two interation points (IPs) occupied by particle detectors. In the intervening arcs the beams are separated by means of electrostatic separators; long-range (also referred to as parasitic) collisions occur at 136 other locations. Effects arising from both head-on and long-range beam-beam interactions impose serious limitations on machine performance, hence constant efforts are being exerted to better understand the beam dynamics. Due to the extreme complexity of the problem a numerical simulation appears to be one of the most reliable ways to study the performance of the system.
Studies of beam-beam interactions in the Tevatron Run II mainly concentrated on the incoherent effects, which were the major source of particle losses and emittance growth. This approach was justified by the fact that the available antiproton intensity was a factor of 10 to 5 less than the proton intensity with approximately equal transverse emittances. Several simulation codes were developed and used for the optimization of the collider performance [2,3].
With the commissioning of electron cooling in the Recycler, the number of antiprotons available to the collider substantionally increased. During the 2007 and 2008 runs the initial proton and antiproton intensities differed by only a factor of 3. Moreover, the electron cooling produces much smaller transverse emittance of the antiproton beam (≃ 4π mm mrad 95% normalized vs. ≃ 20π mm mrad for protons), leading to the head-on beam-beam tune shifts of the two beams being essentially equal. The maximum attained total beam-beam parameter for protons and antiprotons is 0.028.
Under these circumstances coherent beam-beam effects may become an issue. A number of theoretical works exist predicting the loss of stability of coherent dipole oscillations when the ratio of beam-beam parameters is greater than ≃ 0.6 due to the suppression of Landau damping [4]. Also, the combined effect of the machine impedance and beam-beam interactions in extended length bunches couples longitudinal motion to transverse degrees of freedom and may produce a dipole or quadrupole mode instability [5].
Understanding the interplay between all these effects requires a comprehensive simulation. This paper presents a macroparticle simulation that includes the main features essential for studying the coherent motion of bunches in a collider: a self-consistent 3D Poisson solver for beam-beam force computation, multiple bunch tracking with the complete account of sequence and location of long-range and head-on collision points, and a machine model including our measurement based understanding of the coupled linear optics, chromaticity, and impedance.
In Sections II-V we describe the simulation subcomponents and their validation against observed effects and analytic calculations. Section VI shows results from simulation runs which present studies of increasing the beam intensity. Finally, in Section VII we study the coherent stability limits for the case of combined resistive wall impedance and long-range beam-beam interactions.

II. BEAMBEAM3D CODE
The Poisson solver in the BeamBeam3d code is described in Ref. [6]. Two beams are simulated with macroparticles generated with a random distribution in phase space. The accelerator ring is conceptually divided into arcs with potential interaction points at the ends of the arcs. The optics of each arc is modeled with a 6 × 6 linear map that transforms the phase space {x, x ′ , y, y ′ , z, δ} coordinates of each macroparticle from one end of the arc to the other. There is significant coupling between the horizontal and vertical transverse coordinates in the Tevatron. For our Tevatron simulations, the maps were calculated using coupled lattice functions [7] obtained by fitting a model [8] of beam element configuration to beam position measurements. The longitudinal portion of the map produces synchrotron motion among the longitudinal coordinates with the frequency of the synchrotron tune. Chromaticity results in an additional momentum-dependent phase advance δµ x(y) = µ 0 C x(y) ∆p/p where C x(y) is the normalized chromaticity for x (or y) and µ 0 is the design phase advance for the arc. This is a generalization of the definition of chromaticity to apply to an arc, and reduces to the normalized chromaticity (∆ν/ν)/(∆p/p) when the arc encompasses the whole ring.
The additional phase advance is applied to each particle in the decoupled coordinate basis so that symplecticity is preserved.
The Tevatron includes electrostatic separators to generate a helical trajectory for the oppositely charged beams. The mean beam offset at the IP is included in the Poisson field solver calculation. Different particle bunches are individually tracked through the accelerator. They interact with each other with the pattern and locations that they would have in the actual collider.
The impedance model applies a momentum kick to the particles generated by the dipole component of resistive wall wakefields [9]. Each beam bunch is divided longitudinally into slices containing approximately equal numbers of particles. As each bunch is transported through an arc, particles in slice i receive a transverse kick from the wake field induced by the dipole moment of the particles in forward slice j: The length of the arc is L, N j is the number of particles in slice j, r 0 is the classical electromagnetic radius of the beam particle e 2 /4πǫ 0 m 0 c 2 , z ij is the longitudinal distance between the particle in slice i that suffers the wakefield kick and slice j that induces the wake. r j is the mean transverse position of particles in slice j, b is the pipe radius, c is the speed of light, σ is the conductivity of the beam pipe and βγ are Lorentz factors of the beam. Quantities with units are specified in the MKSA system.

III. SYNCHROBETATRON COMPARISONS
We will assess the validity of the beam-beam calculation by comparing simulated synchrobetatron mode tunes with a measurement performed at the VEPP-2M 500 MeV e + e − collider and described in Ref. [12]. These modes are an unambiguous marker of beam-beam interactions and provide a sensitive tool for evaluating calculational models. These modes arise in a colliding beam accelerator where the longitudinal bunch length and the transverse beta function are of comparable size. Particles at different z positions within a bunch are coupled through the electromagnetic interaction with the opposing beam leading to the development of coherent synchrobetatron modes. The tune shifts for different modes have a characteristic evolution with beam-beam parameter ξ = Nr 0 /4πγǫ, in which N is the number of particles, r 0 is the classical electromagnetic radius of the beam particle, and ǫ is the unnormalized one-sigma beam emittance.
There are two coherent transverse modes in the case of simple beam-beam collisions between equal intensity beams without synchrotron motion: the σ mode where the two beams oscillate with the same phase, and the π mode where the two beams oscillate with opposite phases [10]. Without synchrotron motion, the σ mode mode has the same tune as unperturbed betatron motion while the π mode frequency is offset by Kξ, where the parameter K is approximately equal to and greater than 1 and depends on the transverse shape of the beams [11]. The presence of synchrotron motion introduces a more complicated spectrum of modes whose spectroscopy is outlined in Fig. 1 in Ref. [12].
We simulated the VEPP-2M collider using Courant-Snyder uncoupled maps. The horizontal emittance in the VEPP-2M beam is much larger than the vertical emittance. The bunch length (4 cm) is comparable to β * y = 6 cm so we expect to see synchrobetatron modes. In order to excite synchrobetatron modes, we set an initial y offset of one beam sigma approximately matching the experimental conditions. Simulation runs with a range of beam intensities corresponding to beam-beam parameters of up to 0.015 were performed, in effect mimicking the experimental procedure described in Ref. [12]. For each simulation run, mode peaks were extracted from the Fourier transform of the mean bunch vertical position. An example of the spectrum from such a run is shown in Fig. 1 with three mode peaks indicated. In Fig. 2, we plot the mode peaks from the BeamBeam3d simulation as a function of ξ as red diamonds overlaid on experimental data from Ref. [12] and a model using linearized coupled modes referred to as the matrix model described in that reference and Ref. [13,14]. As can be seen, there is good agreement between the observation and simulation giving us confidence in the beam-beam calculation.

IV. IMPEDANCE TESTS
Wakefields or, equivalently, impedance in an accelerator with a conducting vacuum pipe gives rise to well known instabilities. Our aim in this section is to demonstrate that the wakefield model in BeamBeam3d quantitatively reproduced these theoretically and experimentally well understood phenomena. The strong head-tail instability examined by Chao [9]  arises in extended length bunches in the presence of wakefields. For any particular accelerator optical and geometric parameters, there is an an intensity threshold above which the beam becomes unstable.
The resistive wall impedance model applies an additional impulse kick in addition to the application of the map derived from beam optics. The tune spectrum is computed from the Fourier transform of the beam bunch positions sampled at the end of each arc. In order for the calculation to be a good approximation of the wakefield effect, the impedance kick should be much smaller than the x ′ or y ′ change due to regular beam transport so we divide the ring into multiple arcs. which brings up the question is how many is sufficient. The difference in calculated impedance tune shift for a 12 arc division of the ring or a 24 arc division is only 2 × 10 −4 , which is less than 3% of the synchrotron tune (0.007 in this study), the relevant scale in these simulations. We perform the calculation with 12 arcs for calculational efficiency.
In the absence of impedance, we would expect to see the tune spectrum peak at 20.574, the betatron tune of the lattice. With a pipe radius of 3 cm and a bunch length of 20 cm, resistive wall impedance produces the spectrum shown in Fig. 3 for a bunch of 4 × 10 12 protons at 150 GeV [22]. In this simulation, the base tune ν β is 20.574 and the synchrotron tune is 0.007. Three mode peaks are clearly evident corresponding to synchrobetatron modes with frequencies ν β − ν s shifted up by the wakefield (point A), ν β shifted down (point B), and ν β + ν s shifted upward (point C) as would be expected in Ref. [15].
In Fig. 4, we show the evolution of the two modes as a function of beam intensity. With the tune and beam environment parameters of this simulation, Chao's two particle model predicts instability development at intensities of about 9 × 10 12 particles, which is close to where the upper and lower modes meet. We show two sets of curves for two slice and six slice wakefield calculations. The difference between the two slice and six slice simulations is accounted for by the effective slice separation,ẑ, that enters Eq. 1. With two slices, the effectiveẑ is larger than than the six slice effectiveẑ, resulting in a smaller W 0 . With the smaller wake strength, a larger number of protons is necessary to drive the two modes together as is seen in Fig. 4. the threshold of strong head-tail instability has a parabolic dependence on beam intensity.
The wakefield calculation reproduces this feature, as shown in Fig. 5. The growth rate is slowly increasing up to the instability threshold at 5.42×10 12 , after which it has the explicitly quadratic dependence on beam intensity (I) of growth rate = −0.100 + 0.0304I − 0.00207I 2 .
Chromaticity interacts with impedance to cause a different head-tail instability. We simulated a range of beam intensities and chromaticity values. The two particle model and the more general Vlasov equation calculation [9] indicate that the growth rate scales by the head-tail phase χ = 2πCν βẑ /cη, where η is the slip factor of the machine andẑ is roughly the bunch length. The head-tail phase gives the size of betatron phase variation due to chromatic effects over the length of the bunch.
Some discussion of the meaning of the slip factor in the context of a simulation is necessary. In a real accelerator, the slip factor has an unambiguous meaning: The momentum compaction parameter α c is determined by the lattice and γ is the Lorentz factor. We simulate longitudinal motion by applying maps to the particle coordinates z and δ in discrete steps. The simulation parameters specifying longitudinal transport are the longitudinal beta function β z and synchrotron tune ν s . Note that these parameters do not make reference to path length travelled by a particle. However, path length enters into the impedance calculation because wake forces are proportional to path length. In addition, analytic calculations of the effect of wake forces depend on the evolution of the longitudinal particle position which in turn depend explicitly on the slip factor. For our comparisons with analytic results to be meaningful, we need to use a slip factor that is consistent with the longitudinal maps and the path lengths that enter the wake force calculations. The relationship between the slip factor η and the simulation parameters is β z = ηL /2πν s , where L is the length of the accelerator and β z = σ z /σ δ is the longitudinal beta function [16,17] which may be derived by identifying corresponding terms in the solution to the differential equations of longitudinal motion and a one term linear map.
When the growth rate is normalized by Nr 0 W 0 /2πβγν β , which includes the beam intensity and geometric factors, we expect a universal dependence of normalized growth rate versus head-tail phase that begins linearly with head-tail phase [18] and peaks around -1 [23]. The phenomenon of bunch dependent emittance growth is observed experimentally [20].
The beam-beam simulation with 36-on-36 bunches shows similar effects. We ran a simulation of 36 proton on 36 antiproton bunches for 50000 turns with the nominal helical orbit.
The proton bunches had 8.8 × 10 11 particles (roughly four times the usual to enhance the effect) and the proton emittance was the typical 20π mm mrad. The antiproton bunch intensity and emittance were both half the corresponding proton bunch parameter. The initial  The observed bunch emittance variation is similar to the simulation results. Another beambeam simulation with the beam separation at the closest head-on IP expanded 100 times its nominal value resulted in curve (b) of Figure 8 showing a much reduced bunch-to-bunch variation. We conclude that the beam-beam effect at the long-range IPs is the origin of the bunch variation observed in the running machine and that our simulation of the beam helix is correct.

A. Single bunch features
We looked at the tune spectrum with increasing intensity for equal intensity p andp beams containing one bunch each. As the intensity increases, the beam-beam parameter ξ increases. Fig. 9 shows the spectrum of the sum and difference of the two beam centroids for ξ = 0.01, 0.02, 0.04, corresponding to beam bunches containing 2.2 × 10 11 , 4.4 × 10 11 and 8.8 × 10 11 protons. The abscissa is shifted so the base tune is at 0 and normalized in units of the beam-beam parameter at a beam intensity of 2.2 × 10 11 . The coherent σ and π mode peaks are expected to be present in the spectra of the sum and difference signals of the two beam centroids. The coherent σ modes are evident at 0, while the coherent π modes should slightly greater than 1, 2, and 4 respectively. Increasing intensity also causes larger induced wake fields which broaden the mode peaks, especially the π mode, as shown in Fig. 9.
The 4D emittances at higher intensities show significant growth over 20000 turns as shown in Fig. 10. The kurtosis excess of the two beams remains slightly positive for the nominal intensity, but shows a slow increase at higher intensities indicating the the beam core is being concentrated as shown in Fig. 11. Concentration of the bunch core while emittance is growing indicates the development of filamentation and halo.

B. Simulation of bunch length, synchrotron motion and beam-beam interactions
Synchrotron motion in extended length bunches modifies the effects of the beam-beam interaction by shifting and suppressing the coherent modes. The plots in Fig. 12   simulated spectra for sum and difference signals of the beam centroid offsets for one-onone bunch collisions in a ring with Tevatron-like optics, with both short and long bunches, at three different synchrotron tunes. The sum signal will contain the σ mode while the difference signal will contain the π mode. In this Tevatron simulation, the beam strength is set so that the beam-beam parameter is 0.01, the base tune in the vertical plane is 0.576, and β y is approximately 30 cm. Subplots a and b of Fig. 12 show that with small synchrotron tune both the σ and π mode peaks are evident with short and long bunches. The σ mode peak is at the proper place, with the π mode peak shifted upwards by the expected amount, larger than the beam-beam splitting (subplots e and f ), short bunches still exhibit strong coherent modes, but with long bunches the coherent modes are significantly diluted. In the case of long bunches, the σ mode has been shifted upwards to 0.580, and the π mode is not clearly distinguishable from the continuum. At ν s of 0.01 and 0.02, the synchrobetatron side bands are clearly evident. with runs of two-on-two bunches and six-on-six bunches before moving on to investigate the situation with the full Tevatron bunch fill of 36-on-36 bunches. Two-on-two bunches will demonstrate the bunches coupling amongst each other, but will not be enough to demonstrate the end bunch versus interior bunch behavior that characterizes the Tevatron. For that, we will look at six-on-six bunch runs.
In these studies, we are only filling the ring with at most six bunches in a beam. Referring to Fig. 7, we see that only the head-on location at B0 is within the green shaded region where beam-beam collisions may occur with six bunches in each beam. Because of the beam-beam collisions, each bunch is weakly coupled to every other bunch which gives rise to multi-bunch collective modes.
We began the investigation of these effects with a simulation of beams with two bunches each. The bunches are separated by 21 RF buckets as they are are in normal Tevatron operations. Collisions occur at the head-on location and at parasitic locations 10.5 RF buckets distant on either side of the head-on location. To make any excited modes visible, we ran with 2.2 × 10 11 particles, which gives a single bunch beam-beam parameter of 0.01.
There are four bunches in this problem. We label bunch 1 and 3 in beam 1 (proton) and bunch 2 and 4 in beam 2 (antiproton) with mean y positions of the bunches y 1 , . . . y 4 .
By diagonalizing the covariance matrix of the turn-by-turn bunch centroid deviations, we determine four modes, shown in Fig. 13. Fig. 13(a) shows the splitting of the σ mode. The coefficients of the two modes indicate that this mode is primarily composed of the sum of corresponding beam bunches (1 with 2, 3 with 4) similar to the σ mode in the one-on-one bunch case. The other two modes in Fig 13(b) have the character and location in tune space of the π mode, from their coefficients and also their reduced strength compared to the σ mode.
With six on six bunches, features emerge that are clearly bunch position specific. Fig. 14(a) shows the turn-by-turn evolution of 4D emittance and (b) y kurtosis for each of the six proton bunches. It is striking that bunch 1, the first bunch in the sequence, has a lower emittance growth than all the other bunches. Emittance growth increases faster with increasing bunch number from bunches 2-5, but bunch 6 has a lower emittance growth than even bunch 4. The kurtosis of bunch 1 changes much less than that of any of the other bunches, but bunches 2-5 have a very similar evolution, while bunch 6 is markedly closer to bunch 1. One difference between the outside bunches (1 and 6) and the inside bunches (2)(3)(4)(5) is that they have only one beam-beam interaction at the parasitic IP closest to the head-on collision, while the inside bunches have one collision before the head-on IP, and one after it. The two parasitic collision points closest to the head-on collision point have the smallest separation of any of the parasitics, so interaction there would be expected to disrupt the beam more than interactions at other parasitic locations.
To test this hypothosis, we did two additional runs. In the first, the beam separation at the parasitic IP immediately downstream of the head-on IP was artificially increased in the simulation so as to have essentially no effect. The effect of this is that the first proton bunch will not have any beam-beam collisions at an IP close to the head-on IP, while all the other bunches will have one collision at a near-head-on IP. The corresponding plots of emittance and kurtosis are shown in Fig. 15. The kurtosis data shows that bunches 2-5 which all suffer Figure (a) shows the two modes that are most like σ modes. σ mode 1 is 0.53y 1 + 0.53y 2 + 0.59y 3 − 0.31y 4 , σ mode 2 is 0.39y 1 + 0.49y 2 − 0.46y 3 − 0.63y 4 . Figure (b) shows the two π-like modes. π mode 1 is 0.74y 1 − 0.66y 2 − 0.08y 3 , π mode 2 is 0.12y 1 + 0.20y 2 − 0.66y 3 + 0.31y 4 . The absolute scale of the Fourier power is arbitrary, but the relative scales of plots (a) and (b) are the same.
one beam-beam collision at a close parasitic IP are all together while bunch 1 which does not have a close IP collision is separated from the others.
Emittance and kurtosis growth in simulations where the beam separation at the closest upstream and downstream parasitic IPs was increased is shown in Fig. 16. In this configuration no bunch suffers a strong beam-beam collision at a parasitic IP close to the head-on location so the kurtosis of all the bunches evolves similarly.

VII. LOWER CHROMATICITY THRESHOLD
During the Tevatron operation in 2009 the limit for increasing the initial luminosity was determined by particle losses in the so-called squeeze phase [21]. At this stage the beams are separated in the main interaction points (not colliding head-on), and the machine optics is gradually changed to decrease the beta-function at these locations from 1.5 m to 0.28 m. With proton bunch intensities currently approaching 3.2×10 11 particles, the chromaticity of the Tevatron has to be managed carefully to avoid the development of a head-tail instability. It was determined experimentally that after the head-on collisions are initiated, the Landau damping introduced by beam-beam interaction is strong enough to maintain beam stability at chromaticity of +2 units (in Tevatron operations, chromaticity is ∆ν/(∆p/p).) At the earlier stages of the collider cycle, when beam-beam effects are limited to long-range interactions the chromaticity was kept as high as 15 units since the concern was that the Landau damping is insufficient to suppress the instability. At the same time, high chromaticity causes particle losses which are often large enough to quench the superconducting magnets, and hence it is desireable to keep it at a reasonable minimum.
Our multi-physics simulation was used to determine the safe lower limit for chromaticity.
The simulations were performed with starting beam parameters listed in Table I. With chromaticity set to -2 units, and no beam-beam effect, the beams are clearly unstable as seen in Fig. 17. With beams separated, turning on the beam-beam effect prevents rapid oscillation growth during the simulation as shown in Fig. 18. The bursts of increased amplitude is sometimes indicative of the onset of instability, but it is not obvious within the limited duration of this run. The RMS size of the beam also does not exhibit any obvious unstable tendencies as shown in Fig. 19.
Based on these findings the chromaticity in the squeeze was lowered by a factor of two, and presently is kept at 8-9 units. This resulted in a significant decrease of the observed particle loss rates (see, e.g., Fig. 5 in [21]).

VIII. SUMMARY
The key features of the developed simulation include fully three-dimentional strong-strong multi-bunch beam-beam interactions with multiple interaction points, transverse resistive wall impedance, and chromaticity. The beam-beam interaction model has been shown to reproduce the location and evolution of synchrobetatron modes characteristic of the 3D strong-strong beam-beam interaction observed in experimental data from the VEPP-2M collider. The impedance calculation with macroparticles excites both the strong and weak head-tail instabilities with thresholds and growth rates that are consistent with expectations from a simple two-particle model and Vlasov calculation. Simulation of the interplay between the helical beam-orbit, long range beam-beam interactions and the collision pattern qualitatively matches observed patterns of emittance growth.
The new program is a valuable tool for evaluation of the interplay between the beam-