Longitudinal solitons in bunched beams

044402-1 Stable, coherent, longitudinal oscillations have been observed in several accelerators. Within the context of perturbation theory, the beam parameters and machine impedance often suggest these oscillations should be Landau damped. When nonlinear effects are included, long-lived, stable oscillations become possible for low intensity beams. In this paper we report observations of stable humps in bunched beams and present a theoretical framework for their description. Implications for bunched beam stochastic cooling are considered.


I. INTRODUCTION
Solitary waves in the form of holes or humps have been observed in plasmas and coasting beams and the theory has been developed [1][2][3][4][5][6]. These solitons are stationary solutions of a coupled Vlasov-Poisson system.
For bunched beams below transition, the space charge force reduces the incoherent synchrotron frequency and, as such, is generally referred to as a defocusing impedance. Above transition space charge is focusing, leading to an increase in the incoherent synchrotron frequency. The change from focusing to defocusing at transition is an example of the mass conjugation theorem [4,7], and we will generally refer to focusing and defocusing impedances so that our results apply equally well above and below transition. Since we consider relativistic beams, the wall impedance is also included.
Solitons in bunched beams were suggested by observations of ''rf activity'' during bunched beam stochastic cooling studies in the SPS [8] and the Tevatron [9,10]. The rf activity was strong and extended to very high frequency, making the design of a cooling system difficult. A bunched beam stochastic cooling system for RHIC is in the design stage [11] and the possibility of a halo cooling system for the LHC has been considered [12]. It is therefore both timely and prudent to acquire a working knowledge of these phenomena. Before proceeding with the new material we review the coasting beam case.
Let denote machine azimuth and ! 0 be the angular revolution frequency. We use the comoving phase ÿ ! 0 t as the longitudinal coordinate. Assume that the beam current I; t varies slowly in the comoving frame so that the longitudinal voltage per turn is given by where W is the longitudinal wake potential, and we have assumed a short range wakefield to remove multiturn contributions to the voltage. This wakefield includes both the space charge piece as well as the contributions from bellows, transitions, etc. We have made the usual approximation that the frequency spread in the beam is small enough so that the locations of the impedance producing objects in the ring are irrelevant and only the sum of all the individual contributions matters. The equation of motion for a single particle is where 1= 2 T ÿ 1= 2 is the frequency slip factor, q is the charge per particle, E is the central energy, and v=c. Set p d=dt and introduce the phase space density f; p; t, where fdpd is the number of particles in dpd. The Vlasov equation for f is then @f @t p @f @ ÿ ! 2 0 q 2 2 E V; t @f @p 0; where Soliton solutions are independent of time and nonlinear. When W W 0 2 sgn expÿjj; where > 0 and W 0 are constants, the voltage is related to the current via For this case analytical progress is possible [4]. Certain other wakefields are tractable as well [2 -4]. For the present case we will set 0 so that the coasting beam Hamiltonian is approximated as The constant is given by where L is the net inductance. The net inductance satisfies j! 0 L Z=n, where Z is the sum of the broadband wall impedance and the longitudinal space charge impedance. For a round Gaussian beam with rms radius , in a round pipe of radius b , the longitudinal space charge impedance is well approximated by (1) where 1= 1 ÿ 2 p is the Lorentz factor, Z 0 377, and the derivation is given in the appendix. Figure 1 shows a model due to Sacherer [13], in the frame comoving with the soliton. The piecewise constant, phase space density is either 0, f 0 , or f 0 f 1 ; and the distribution is independent of time. Since the phase space density is constant on contours of constant H, one obtains the equations H 1 ; p 1 H 2 ; 0 and H 1 ; p 2 H 2 ; p 0 . These may be combined to give p 2 2 ÿ p 2 0 p 2 1 2kf 0 p 0 ÿ p 2 ÿ f 1 p 1 : Taking p 0 , f 0 , and p 1 as input parameters yields Taking p 1 p 0 and expanding the square root gives If the correction term kf 0 =p 0 is set to zero, the condition is identical to that for a phase space density of f 1 to self bunch. The change in the line density due to the soliton is ÿp 2 1 =k. For an inductive impedance above transition k > 0. This results in a deficit of phase space density and one observes a hole in the line density. When Þ 0 the problem is more difficult. In a forthcoming publication [14], in which the delta effect is taken into account, we shall investigate the extent to which these predictions of the Sacherer model can be maintained. We go on to present some bunched beam data for the RHIC.

II. DATA
Long-lived coherence has been observed in the SPS [8,12], the Tevatron [9,15], and now the RHIC [16]. Figure 2 shows a mountain range plot of the line density for freshly injected protons with 25:9. The amplitude of the coherent oscillation increased steadily, and  flattop with 107. In all the cases shown, only the 28 MHz accelerating cavities were operating and the total acquisition time was 4000 turns 50 ms. RHIC's transition energy is T 23:8, so all the data are above transition. All the data show a coherent oscillation which corresponds to a region of overdensity, or hump, in the longitudinal phase space. This behavior is commonplace in RHIC and we have never observed a stable hole. Measurements of RHIC's broadband impedance [17] give Z=n j3 1 for the inductive wall contribution. The longitudinal space charge impedance at 25:9 is Z=n ÿj1:3 and the space charge impedance becomes negligible at store. For all the data shown the coherent force leads to a reduction in the incoherent synchrotron frequency and as such it is defocusing, in the usual sense. We have also obtained data showing humps in deuteron beams with 10:7 < T and Z=n ÿ5j. Again, the deuterons show humps for impedances that are usually thought of as defocusing. Therefore, we always see humps with a defocusing impedance, which is just the reverse of what one expects for coasting beams. We go on to consider the theory of solitons in bunched beams.

III. BUNCHED BEAM MODEL
For bunched beams we simply include the rf voltage in the coasting beam equations of motion. To simplify notation let denote the position in the bunch measured in units of rf radians, ! s;0 denote the small amplitude angular synchrotron frequency, and use s ! s;0 t as the evolution variable. Let ; s be the normalized line density of the particles Z ÿ d; s 1: The rf voltage is Vt V rf sin and ! rf is the angular rf frequency. The total charge in the bunch is Q which is assumed positive, and the effect of transition energy is included in the sign of the rf voltage with V rf > 0 below transition and V rf < 0 above transition. Since solitons are long lived we neglect the effect of any synchronous phase. The equation of motion for is To simplify notation set Space charge dominated beams below transition have ' > 0. For a steady state, matched bunch, a positive value of ' defocuses the beam and leads to an incoherent synchrotron frequency that is less than the synchrotron frequency for ' 0. Set p d=ds so that dp ds ÿ sin ÿ ' @; s @ : A drift-kick computer code has been written to simulate these equations of motion. Figure 6 shows simulation results for a short, roughly matched bunch that is displaced in phase from the stable fixed point. As is clear from the figure, a positive value of ' produces solitons of high density. Solitons in the RHIC accelerator always have ' > 0. However, for some ranges of the simulation parameters a hump can form with a focusing impedance. Next we present an analytic model of this phenomena.

IV. EQUATIONS OF MOTION
Taking p as the momentum coordinate and as the position coordinate, Eq. (3) follows from the Hamiltonian H; p; s p 2 =2 1 ÿ cos '; s: Introduce the distribution function f; p; s, where f; p; sddp is the fraction of particles in ddp. The line density is then Consider a canonical transformation of Goldstein's [18] first type with a generator given by These are rotating coordinates so both J and are ''slow'' variables as long as the bunch is not too long. Since the transformation is canonical dpd dJd . The Hamiltonian in the new variables is Assume that f varies slowly compared with the sinusoids which is equivalent to assuming that the coherent frequency of the soliton is close to the small amplitude synchrotron frequency. Define the phase average of a generic function hJ; ; s as hhJ; ; si s 1 2 with the slow variables held constant during the average. For the first line in Eq. (7) ÿJ hJcos 2 si s 1 ÿ hcos 2J where J 0 x is the Bessel function, and we have introduced J both for convenience and to allow for other forms of rf voltage. Averaging the Hamiltonian in (6) gives To perform the phase average of the delta function notice 2J where is a phase that depends on J, J 1 , , and 1 . And The last steps follow since and R are slow, and the delta function in the integrand on the right is nonzero only when s =2. The Hamiltonian becomes FIG. 6. Scatter plots of p versus from longitudinal simulations with a sinusoidal rf voltage and a broadband impedance. The initial conditions (top panels) are identical and the only difference during the evolution is the sign of the broadband impedance. Effects of a defocusing impedance are shown on the left and those for a focusing impedance are shown on the right. Plots at 400 synchrotron periods (middle panels) and 800 synchrotron periods (bottom panels) clearly show that a defocusing impedance can lead to solitons.
The equations of motion for and J are The corresponding Vlasov equation is We seek uniformly rotating solutions (i.e., the solitons), f ; J; s g r rs; J; wherer r is a constant. Inserting this functional form in the Vlasov equation gives r where g g ; J is a function of two variables. For some purposes it is useful to transform back to Cartesian phase space coordinates.
This is a canonical transformation so the density gA; B satisfies @K @B @g @A ÿ @K @A @g @B 0; with the new Hamiltonian where the single particle motion is generated by and collective effects are due to

V. APPROXIMATE SOLUTIONS
For physical solutions of Eq. (12), local extrema of g occur at stable fixed points of K. Additionally, contours of constant values of g and K coincide. Suppose we have a situation like that shown in the bottom left panel of Fig. 6. Choose the A axis to pass through the center of the soliton so it is centered at A A 0 , B 0. We will consider solutions for which this is a fixed point of both R and V. In the vicinity of the fixed point the leading order approximation to the single particle generator is To this order @R=@A 2A ÿ A 0 and @R=@B 0. In this approximation the curvature of the soliton apparent on the left side of Fig. 6 will be absent, but for solitons of small angular extent, the results should be good.
If one could find a class of functions gA; B which were constant on elliptical contours, and yielded linear coherent forces within the beam, the Vlasov equation would reduce to an algebraic equation among the various parameters. Much of the relevant mathematics is the same as for transverse space charge [19] but will be reproduced here for completeness. First, notice that Eq. (15) is the electrostatic potential of a two dimensional charge sheet in the plane containing the sheet. Let us add a third dimension z.
The generalization of Eq. (15) is which satisfies the equation VA; B; z ÿ4'gA; B; z: (17) Consider distributions of the form, and define The field gradients are given by Substituting Eqs. (20) - (22) in (17), and integrating by parts proves that the solution works up to the boundary condition terms. To prove the boundary conditions assume that nT 0 for T > T 0 , and that nT n 0 for all T. Assume B 2 =b 2 > T 0 , then which goes to zero as B ! 1. Letting A ! 1 yields the same limit and proves the boundary conditions.
Next consider the potential in the z 0 plane VA;B;0 ' As c ! 0 we obtain a charge sheet with the two dimension potential where the last integral is independent of c. For distributions where nT is a constant for T < T 0 and zero otherwise, the gradients (20) - (22) are linear functions in the region where n Þ 0. Integrating over the z dimension yields the class for functions needed for the solution of the Vlasov equation where a and b are parameters, and we have shifted the center of the distribution to the stable fixed point.
In the region where g Þ 0 the coherent forces are With Eqs. (25) -(27) the Vlasov equation (12) is equivalent to where we have introduced the notation V AA @ 2 V=@A 2 A A 0 ; B 0. The potential VA; B is ''rounder'' than the distribution gA; B. More explicitly, if a < b then We have verified this via numerical calculation of the various integrals [20]. Intuitively, notice that VA; B is the convolution of gA; B with a cylindrically symmetric function.
Generally, one has V AA ÿ'=a 3 G a a=b and V BB ÿ'=a 3 G b a=b. Define r a=b G a r 3r 3 2 The summations are combinations of various elliptic integrals [21]. Equation (28) reduces to ÿ 2a 3 ' G b r=r 2 ÿ G a r ÿ0:77r lnr=1 r; where the approximate expression is accurate with a few percent for 0 < r < 10. The right side of this equation is plotted as a function of r in Fig. 7. The maximum value of 0.212 24 occurs for r 0:269 82. Figure 8 shows the individual G a and G b . For r > 0:4, r 2 G a =G b r 1=2 , in accord with inequality (29). Notice that any r > 0 is permissible. For r < 1 the soliton is due to a defocusing impedance while a focusing impedance is required for r > 1.
When a background beam is present, as shown in Figs. 2 -5, analytic solutions of Eq. (12) have not been found. In the next section we present some numerical results.

VI. SELF-CONSISTENT NUMERICAL SOLUTIONS
Equation (12) implies that contours of constant g coincide with those of constant K. However, this does not imply that there always exists a single function GK such that gA; B GKA; B. To see this consider Fig. 9 which shows the Hamiltonian through a line containing both the origin and the center of the soliton for both positive (solid red line) and negative (blue dashed line) values of '. For the moment consider the defocusing (solid red line) case, which has ' > 0. There are four values of A for which K 0:01 but the density gA; B need not have the same value at all four points. The key is the separatrix, on which K K s , which is shown in Fig. 10.
In each of the three regions bounded by the separatrix one is free to choose a different functional form g G i K. There are limitless possibilities but we take the fairly simple prescription gA;B C 4 8 < : in region 1; where RA; B is the single particle generator defined in Eq. (14). We obtain self-consistent solutions for ' > 0 using the following procedure:  (1) Chooser r. Set n 0. Choose initial parameters C 0 0 , C 1 , and C 2 . Set up an initial line density g 0 A; B which is symmetric around B 0.
(2) Use g n A; B to calculate K n A; B through (34). Find the separatrix on which K n A; B K n s . In region 1, K n A; B is maximum at the stable fixed point (A A n 0 ; B 0) and minimum at the unstable fixed point. Both fixed points always lie on the line B 0.
Define the electrostatic amplitude n K n O ÿ K n X where O and X denote the stable and unstable fixed points, respectively.
(3) If n > 0 change C 0 through the formula where is an input parameter. If n 0 use C 1 0 C 0 0 . (4) Calculate g n1 using Eq. (33) with K n A; B and C n 0 . Calculate the global normalization constant C 4 and increase n.
The variation in C 0 , defined through Eq. (35), greatly increases the initial parameter space leading to interesting solutions; typically, 1. Figure 11 shows the evolution of the density as the iteration proceeds.
Depending on the initial parameters, a large range of stable solutions can be obtained. An outcome of this numerical procedure is that there is a definite relationship betweenr r and the converged quantities C 0 and , denoted by C 0 ;r r shown in Fig. 12. It can be considered as a nonlinear dispersion relation (as in the standard soliton cases [1][2][3][4]), since it relates the soliton amplitude and the trapping parameter C 0 with the phase speedr r. As seen from (33), C 0 can be expressed by where g0 [gX] is the distribution function at the stable (unstable) fixed point. It is clear that C 0 is positive for a hump in region 1 of Fig. 10, zero for a flattopped and negative for a notch, and hence C 0 reflects the status of trapped particles. We learn from Fig. 12 that C 0 is strictly positive, corresponding to a hump (hot spot), and that there are two different values for (for a givenr r and C 0 ). For the larger the majority of particles is residing in the hot spot (as in Fig. 13), whereas for the smaller value of the particles are mostly contained in the background. Figure 13 shows the phase space density and accompanying line density for a bunch with most of the beam trapped in the soliton, while Fig. 14 shows a small soliton. Both solutions have ' 0:01 and A 0 0:025 showing the broad range of stable solutions that are possible for the same bunch charge and machine impedance. This is qualitatively consistent with Figs. 2 -5. In the next section we will compare the results of the iterative solutions with simulations and with the approximate analytic solution.

VII. COMPARISON WITH SIMULATIONS
Longitudinal simulations were carried out with the update equations: p n1 p n ÿ sin n F c n ; (37) n1 n p n1 ; where 2=1000 is the update in s. The coherent force F c was obtained by (i) using linear interpolation to allocate the particle charges to a fine grid, (ii) smoothing the finely gridded data, (iii) applying a centered difference numerical differentiation, and (iv) using linear interpolation to get the forces.
For some of the simulations all the test particles had the same charge and for others the values of the test charge differed with the particle [22]. All particles had the same charge to mass ratio. In all cases the force of particle i on particle j was equal and opposite to the force of particle j on particle i. Since the simulation is in one dimension  there were no problems associated with thermalization effects so that having a range of particle charges produced no secular evolution. Figure 15 compares the simulations to the analytic formula (32). The initial conditions were set up using Eq. (32) although a small amount of curvature was included by using the mapping A cosB=A, p A sinB=A. The simulation parameters are shown in Table I. Simulations 3 and 5 have ' < 0 while the others have ' > 0. For ' > 0 one has a < b and all the simulations but 7 show no filamentation. Simulations 3 and 5 both show some filamentation. Examining the blue, dashed curve in Fig. 9 suggests an explanation. The curves in Fig. 9 were obtained from Eq. (13) using VA; B numerically integrated from a solution of Eq. (32). The solutions of (32) had ab p 0:08 and ' 8 10 ÿ4 . For ' > 0, a=b 0:833; for ' < 0, a=b 1:33. Consider the potential difference between the edge of the soliton and the separatrix. For ' > 0 this value was 5.1 times larger than its value in the ' < 0 case. In other words, with all else equal, ' > 0 solitons make deeper buckets than ' < 0 solitons. These simulations confirm the accuracy of Eq. (32) for isolated solitons. Next we consider the utility of this equation for the case when a background beam is present. Figure 16 shows simulation results for a soliton with a background beam. To obtain the initial conditions we used a modified form of Eq. (32): where was the fraction of the total bunch charge in the hump, and ' is calculated using to total charge in the bunch. The parameters a, r, and were calculated using  the soliton parameters. In this case we used initial conditions with the macroparticles evenly spaced on a rectangular grid, and each macroparticle was assigned a charge proportional to the value of f; p at the 1:1 10 6 macroparticle locations. The coherent force was calculated using 600 fine bins and smoothed to obtain 50 statistically independent bins within the bunch. There was no sign of thermalization, probably because the simulation was in one dimension. From Fig. 16 it is clear that Eq. (39) does not yield an exact match, but the two mountain range plots are quite similar even for as small as 1%. Larger values of produced better agreement between the simulations and Eq. (39). This is reminiscent of the coasting beam result in the introduction, where the soliton parameters were close to the self-bunching condition for the density difference. Comparable simulations have been done using the solutions of the iterative code as initial conditions. The iterative code produced solitons that were well matched, providing an additional check to the correctness of the theory. The theory, of course, considers a smooth phase space distribution but actual beams are composed of particles. We go on to consider the long time evolution of these solitons when the discrete nature of the distribution is included.

VIII. COLLISIONAL EFFECTS
The proton bunches in Figs. 2 -5 show well-defined solitons. RHIC also accelerates and stores gold ions, but no long-lived solitons have been observed in gold beams. The present section considers a possible reason for this discrepancy. The Coulomb collision rate for a single species plasma is [23] where n is the number of ions per unit volume, r p is the classical proton radius, Z is the atomic number of the ion, A is the atomic mass, c is the speed of light, v is the rms velocity, and ln is the Coulomb logarithm. Assuming fixed bunch charge, spatial size, and velocity distribution, the collision rate scales as Z 3 =A 2 , which is 13 for gold and 1 for protons.
A more sophisticated approach employs the Fokker-Planck equation to describe the evolution of the particle distribution due to this intrabeam scattering [24]. The code assumes that the transverse distributions are Gaussian with equal emittances, and that the longitudinal distribution is a function of action alone. The initial longitudinal distribution was set to be a smooth background with a density enhanced core and a Fokker-Planck solver was used to evolve the system. Figure 17 shows the fraction of the bunch charge remaining in the density enhanced core as a function of time for both gold and protons and for a range of intensities. As one can see from the curves, the gold bunch always smooths out more  quickly than the proton bunch, and the factor of 13 obtained earlier seems reasonable. Along with the difference in collision rates the values of ' were usually quite different for protons and gold, with the proton values being larger by about an order of magnitude. This was predominantly due to the low rf voltage required for protons, which were injected just above the transition energy. Hence, for a given value of , the dimensions of a soliton in a proton bunch would be about twice as large as a soliton with the same in a gold bunch. If we take a simple diffusion model, the rate at which a soliton dissipates will scale as coll =a 2 , with a being a typical dimension of the soliton. With this extra factor, the dissipation rate for solitons in a gold bunch is about 50 times faster than the dissipation rate for a comparable soliton in a proton bunch. Deuterons have not been studied as carefully as protons, but the data are roughly consistent with this picture.

IX. CONCLUSIONS
A theoretical framework for understanding stable solitons in bunched beams has been presented. A key result is that a defocusing impedance can lead to humps, while such an impedance can only give holes in coasting beams. Data from RHIC show that humps with a defocusing impedance are present in actual bunches. An analytic formula describing isolated solitons has been presented and verified using numerical simulations. Iterative solutions of the Vlasov equation also agree with simulations so there are three different theoretical approaches that are in agreement with the data.

ACKNOWLEDGMENTS
We thank J. M. Brennan, S. Koscielniak, and S. Peggs for insightful comments, and F. Severino for developing  The plot on the left shows the evolution for the first synchrotron oscillation and the plot on the right shows the evolution over the last of 100 synchrotron oscillations. There are 9 10 5 macroparticles of varying charge and the smoothing length for the coherent force calculation is 1=50th of the bunch length.