Nonlinear perturbative particle simulation studies of the electron-proton two-stream instability in high intensity proton beams

014401-1 Two-stream instabilities in intense charged particle beams, described self-consistently by the nonlinear Vlasov-Maxwell equations, are studied using a 3D multispecies perturbative particle simulation method. The recently developed Beam Equilibrium, Stability and Transport code is used to simulate the linear and nonlinear properties of the electron-proton (e-p) two-stream instability observed in the Proton Storage Ring (PSR) experiment for a long, coasting beam. Simulations in a parameter regime characteristic of the PSR experiment show that the e-p instability has a dipole-mode structure, and that the growth rate is an increasing function of beam intensity, but a decreasing function of the longitudinal momentum spread. It is also shown that the instability threshold decreases with increasing fractional charge neutralization and increases with increasing axial momentum spread of the beam particles. In the nonlinear phase, the simulations show that the proton density perturbation first saturates at a relatively low level and subsequently grows to a higher level. Finally, the nonlinear spacecharge-induced transverse tune spread, which introduces a major growth-rate reduction effect on the e-p instability, is studied for self-consistent equilibrium populations of protons and electrons.


I. INTRODUCTION
In periodic focusing accelerators and transport systems [1][2][3][4], when a second charge component is present, it has been recognized for many years, both in theoretical studies and in experimental observations [5][6][7][8][9][10][11][12][13][14][15][16][17][18], that the relative streaming motion of the high intensity beam particles through a background charge species can provide the free energy to drive the classical two-stream instability, appropriately modified to include the effects of dc space charge, relativistic kinematics, presence of a conducting wall, etc.A background population of electrons can result by secondary emission when energetic particles strike the chamber wall or through ionization of background neutral gas by the beam ions.A welldocumented example is the electron-proton (e-p) twostream instability observed in the Proton Storage Ring (PSR) [12 -14], although a similar instability also exists for other ion species, including (for example) electronion interactions in electron storage rings [15][16][17][18].
At the high beam currents and charge densities of practical interest, it is increasingly important to develop an improved theoretical understanding of the influence of the intense self-fields using a kinetic model based on the nonlinear Vlasov-Maxwell equations [1,[19][20][21][22].Recently, the f formalism, a low-noise, nonlinear perturbative particle simulation technique for solving the Vlasov-Maxwell equations, has been developed for intense beam applications [23,24].The 3D multispecies nonlinear f formalism has been implemented in the Beam Equilibrium, Stability and Transport (BEST) code [25], which has been applied to a wide range of important collective processes in intense beams [1,[25][26][27].In this paper, we study the electron-proton two-stream instabil-ity numerically using the nonlinear f method, with particular emphasis on the parameter regime characteristic of the PSR experiment [12 -14] for a long, coasting beam.Following a brief description of the nonlinear f method in Sec.II, we present the simulation results in Sec.III.In particular, the dependences of the instability growth rate on beam intensity, fractional charge neutralization, and longitudinal momentum spread are investigated in detail.The nonlinear phase of the instability is studied as well.Finally, in Sec.IV, we present a detailed description of the nonlinear space-charge-induced transverse tune spread, which provides an effective damping mechanism for the instability that can significantly reduce the growth rates of the unstable modes.

II. NONLINEAR f FORMALISM FOR THE VLASOV-MAXWELL SYSTEM
The theoretical model employed here is based on the nonlinear Vlasov-Maxwell equations [1].We consider a thin, continuous, high intensity ion beam j b, with characteristic radius r b propagating in the z direction through background electrons j e, with each component described by a distribution function f j x; p; t [1,5,20,22] in the phase space x; p.The charge components j b; e propagate in the z direction with characteristic axial momentum j m j j c.While the nonlinear f formalism outlined here is readily adapted to the case of a periodic applied focusing field [28], for present purposes we make use of a smooth-focusing model in which the applied focusing force is described by F foc j ÿ j m j ! 2 j x ?, where x ?x ê e x y ê e y is the transverse displacement from the beam axis, and !j const is the effective applied betatron frequency for transverse oscillations.In the electrostatic and magnetostatic approximations, we represent the self-electric and self-magnetic fields as E s ÿrx; t and B s r A z x; t ê e z .The nonlinear Vlasov-Maxwell equations can be approximated by [1,5] @ @t v @ @x ÿ j m j ! 2 j x ?e j r ÿ v z c r ?A z @ @p f j x; p; t 0; (1) and r 2 ÿ4 X j e j Z d 3 pf j x; p; t; X j e j Z d 3 pv z f j x; p; t: (2) Detailed descriptions of the theoretical model can be found in Refs.[1,25].In the nonlinear f formalism [23][24][25][26][27], we divide the total distribution function into two parts, f j f j0 f j , where f j0 is a known equilibrium solution @=@t 0 to the nonlinear Vlasov-Maxwell equations ( 1) and (2), and the numerical simulation is carried out to determine the detailed nonlinear evolution of the perturbed distribution function f j .This is accomplished by advancing the weight function defined by w j f j =f j , together with the particles' positions and momenta.The equations of motion for the particles, obtained from the characteristics of the nonlinear Vlasov equation (1), are given by dx ?ji dt j m j ÿ1 p ?ji ; dz ji dt v zji j c ÿ3 j m ÿ1 j p zji ÿ j m j j c; dp ji dt ÿ j m j ! 2 j x ?ji ÿ e j r ÿ v zji c r ?A z : ( Here the subscript ''ji'' labels the ith simulation particle of the jth species.Furthermore, the dynamical equations for w ji are [23,[25][26][27] where ÿ 0 and A z A z ÿ A z0 .Here, the equilibrium solutions ( 0 , A z0 , f j0 ) solve the steady-state (@=@t 0) Vlasov-Maxwell equations ( 1) and (2).Awide variety of axisymmetric equilibrium solutions to Eqs. ( 1) and (2) has been investigated in the literature.The perturbed distribution f j is obtained through the weighted Klimontovich representation [1] where N j is the total number of actual jth species particles, and N sj is the total number of simulation particles for the jth species.Maxwell's equations are also expressed in terms of the perturbed fields and perturbed density according to r 2 ÿ4 X j e j n j ; where Here, Sx ÿ x ji is a shape function distributing particles on the grids in configuration space [25].
The nonlinear particle simulations are carried out by iteratively advancing the particle motions, including the weights they carry, according to Eqs. ( 3) and ( 4), and updating the fields by solving the perturbed Maxwell's equations (6) with appropriate boundary conditions at the cylindrical, perfectly conducting wall.Even though it is a perturbative approach, the f method is fully nonlinear and simulates completely the original nonlinear Vlasov-Maxwell equations.Compared with conventional particle-in-cell simulations, the noise level in f simulations is significantly reduced.The f method can also be used to study detailed linear stability properties, provided the factor 1 ÿ w ji in Eq. ( 4) is approximated by unity, and the forcing terms in Eq. ( 3) are replaced by the unperturbed force.Implementation of the 3D multispecies nonlinear f simulation method described above is embodied in the BEST code [25][26][27] developed at the Princeton Plasma Physics Laboratory.On the IBM SP supercomputer at the National Energy Research Scientific Computing Center, the BEST code typically advances 4:2 10 11 particles time steps when simulating the e-p two-stream instability in the Proton Storage Ring experiment.

III. SIMULATION OF THE TWO-STREAM INSTABILITY
In high intensity accelerators and storage rings, there exist many discrete collective eigenmodes (excitations) of the ion beam.Among them, the dipole surface mode can exhibit instability when a background electron population is present [1,[5][6][7][8][9][10][11][12][13][14].This instability is basically of the two-stream type and is strongest when the ions are relatively cold in the propagation direction.The directed velocity difference, V b ÿ V e , between the beam ions and the background electrons provides the free energy for the collective modes to grow.The instability observed in the Proton Storage Ring [12 -14] is believed to have this two-stream characteristic.We have simulated the two-stream instability using the BEST code, which self-consistently solves the nonlinear Vlasov-Maxwell equations using the f method.The simulation results presented here are for system parameters typical of the PSR experiment for a long, coasting beam.The background distribution functions f j0 r; p under quasi-steady-state conditions @=@t 0 are assumed to be bi-Maxwellian with temperature T j? const in the x-y plane, and temperature T jk const in the z direction, i.e., Here, n n j is the particle number density on axis r 0 of the jth species, and 0 and A z0 are equilibrium self-field potentials, which are coupled with Eq. ( 8) through the nonlinear Maxwell equations In the simulations, we take is the on-axis r 0 relativistic proton plasma frequency.The equilibrium self-field potentials 0 and A z0 as functions of the radial coordinate r, obtained by numerically solving Eqs. ( 8) and ( 9), are plotted in Fig. 1.
The equilibrium density profiles n 0 j r j b; e for the protons and the electrons ; can be calculated self-consistently from the solutions for 0 r and A z0 r: Previous numerical studies for the baseline case [25,27] have shown that the e-p two-stream instability has a dipole-mode structure with azimuthal mode number l 1.It is important to emphasize that the simulations are based on first principles -the nonlinear Vlasov-Maxwell equations.Unlike an instability analysis based on the beam centroid model where only the dipole mode is allowed, all possible mode excitations are allowed in the simulations.Simulations using typical operating parameters in the PSR experiment [12 -14] for a long, coasting beam indicate that the l 1 dipole mode is indeed the most unstable mode.Detailed properties of the linear phase of the instability for the baseline case, such as the dependence of the instability growth rate on the electron density profile and the normalized axial wave number k z V b =! b (k z 2n=L n=R, where R is the ring radius, and n is the mode number), have been numerically investigated and reported in Refs.[25,27].The numerical results were found to be in good agreement with the experimental observation in PSR and qualitatively consistent with the analytical results obtained for uniform-density beams [5,6].
In the present study, we have carried out detailed numerical investigations of the e-p instability for a wide range of beam intensities and fractional charge neutralization.The space-charge intensity varies from moderate to strong, corresponding to 0:008 , but a decreasing function of the longitudinal momentum spread, which qualitatively agrees with previous analytical results [7].A larger longitudinal momentum spread induces stronger Landau damping by parallel kinetic effects and therefore reduces the growth rate of the instability, whereas higher beam intensity provides more free energy to drive a stronger instability.
In addition to the effects of longitudinal Landau damping by the beam ions, we are also able to simulate the stabilizing effects due to the nonlinear space-charge-induced tune spread, because the simulations are carried out for realistic beams with radially varying equilibrium density profiles.As a result of the presence of several important damping and growth-rate reduction mechanisms, an instability threshold is observed in the simulations.Plotted in Fig. 3 is the instability threshold in terms of the normalized beam density n n b = n n b0 as a function of momentum spread p bk =p bk for different values of fractional charge neutralization f.Evidently, larger momentum spread and smaller fractional charge neutralization imply a higher density threshold for the instability to occur.For a specified value of f, if p bk =p bk ; n n b = n n b0 fall below the curves in Fig. 3, then there is no two-stream instability.
Finally, in Fig. 4, we simulate an unstable case to its fully nonlinear phase.This case corresponds to n b 0:079, f 0:1, and p bk 0 p ek at t 0: In Fig. 4, the time history of the density perturbations at a fixed spatial location is shown for both species.There are basically two phases for the evolution of the instability.The first phase is the linear stage where the density perturbations for both species grow exponentially.However, due to the large mass ratio between the protons and the electrons, the density perturbation amplitude for the electrons is much larger than that for the protons.When the linear growth saturates, the saturation level for the electron density perturbation is therefore much larger.The saturation level for the electron density perturbation shown in Fig. 4 is about 8%, whereas the saturation level for the proton density is very small (less than 0:1%).
The second phase of the instability is the nonlinear phase, starting approximately at t 500=! b , during which the electron density perturbation level stays nearly constant around the 8% level.In other words, in the nonlinear phase the electron density perturbation shows no extra dynamical behavior other than the initial nonlinear saturation.This is probably because the initial saturation level for the election density perturbation is already quite large compared with the proton density perturbation.On the other hand, in the nonlinear phase the proton density perturbation grows first slowly and then very fast after t 1400=!b to a high level, considerably larger than that of the electron density perturbation.This simulation result suggests that the late-time growth of the e-p instability observed experimentally in PSR has likely passed the initial linear growth and saturation phase, and entered the second stage of strong nonlinear growth evident in Fig. 4. It also points to the possible physical mechanism proposed by Channell [29] that due to the large mass ratio, the electron density perturbation quickly saturates long before the proton density perturbation becomes sizable, and the large electron density fluctuation level then provides a newly developed background force that drives the proton density perturbations to a large level on a longer time scale [29].

IV. NONLINEAR SPACE-CHARGE-INDUCED TUNE SPREADS AND THEIR EFFECT ON GROWTH-RATE REDUCTION
It has long been recognized that transverse tune spread can have an important damping effect on various types of beam instabilities [9,10].In many cases, instabilities can be completely suppressed by the transverse tune spread.In accelerator experiments, tune spreads are often introduced through the machine chromaticity as a result of the beam momentum spread in the longitudinal direction.However, the tune spread induced by nonlinear spacecharge effects is often neglected.As the beam intensity of contemporary accelerators increases, the transverse tune spread due to the nonlinear space-charge force becomes increasingly large and in some cases can dominate the tune spread induced by the longitudinal momentum spread.As mentioned earlier, the growth rates of the e-p two-stream instability observed in the PSR experiment and in the numerical simulations are generally significantly smaller than theoretical predications based on a centroid model or a kinetic model using the Kapchinskij-Vladimirskij (KV) distribution [5].This difference exists for arbitrary momentum spread when the beam intensity is above the instability threshold.One of the physical effects that contributes to this difference is the spacecharge-induced tune spread which is not incorporated in a centroid model or a kinetic KV model with a flattop density profile.When the beam density profile is bell shaped rather than flattop, particles experience different local transverse oscillation frequencies at different radial locations due to the nonlinear space-charge potential.The corresponding spread in transverse oscillation frequency has been shown to reduce the instability growth rate in Sec.V of Ref. [5].It is clear from the analysis in Ref. [5] that any tune spread associated with nonlinearity in the space-charge field will provide a growth-rate reduction mechanism.However, a non-self-consistent model equilibrium close to a step function is adopted in Ref. [5] in order to analytically carry out the derivation.For the self-consistent bell-shape beam density profile used in the present numerical simulations, it is not possible to carry out a similar quantitative analysis because there are no analytical expressions for the eigenmode structure and corresponding dispersion relation.Nevertheless, the qualitative features of the growth-rate reduction mechanism are expected to be the same, and it is still useful to obtain a measure of the space-charge-induced tune spread.In this section, we estimate the space-chargeinduced tune spread for self-consistent equilibrium populations of protons and electrons, in a parameter regime typical of the PSR experiment.
To simplify the analysis, we use cylindrical polar coordinate r; in the transverse plane.This is because the equilibrium beam profiles are axisymmetric @=@ 0 and the canonical angular momentum is conserved.Therefore, we describe the transverse tunes (oscillation frequencies) in terms of the tunes for the r; motion.For the axisymmetric equilibrium described by Eqs. ( 8) and ( 9), the transverse equations of motion for a single particle of species j j e; b are j m j r r ÿ r _ 2 j m j ! 2 j r e j @ 0 r @r ÿ j @A z0 r @r 0; (10) where 0 r and A z0 r are the equilibrium self-field potentials determined from Eqs. ( 8) and (9).Examples of self-consistent equilibrium solutions are plotted in Figs. 1.For present purposes, we consider the class of particle orbits that pass through the beam axis r 0 with P 0 and _ v =r 0. In this case, the radial orbit equation in Eq. (10)  where !!j r is the local depressed transverse betatron frequency defined by ! !j ! 2 j e j j m j r @ @r 0 r ÿ j @ @r A z0 r 1=2 : (13) Note from Eq. ( 13) that ! !2 j r is a known function of r once 0 r and A z0 r are determined from Eqs. ( 8) and (9).Plotted in Fig. 5  The distribution of particles as a function of !!j can be related to the equilibrium density profile by which gives Here r r !!j is the inverse of the function !! j!!j r, and N j 2 R r w 0 rn 0 j rdr is the axial line density of the jth species.Plotted in Fig. 6  From the plots, we note that a larger beam intensity corresponds to a larger spread in !!j for both species, and the distribution functions F j !!j vary significantly for different beam intensities in the range of 0:1 n n b = n n b0 2:0.After determining the distribution F j !!j , the average transverse oscillation frequency !!j can be determined from In Figs. 7 and 8, the average oscillation frequency and rms frequency spread are plotted versus normalized beam intensity for the case where f n n e = n n b 0:1.As expected, higher beam intensity results in a larger tune depression and frequency spread for the protons and a larger oscillation frequency for the electrons.While the frequency spread for the electrons is much larger than that of the protons, it stays relatively constant for different beam intensities.In the baseline case, the tune depression and frequency spread for the protons are !! e 24%.In summary, at the high beam intensities in the PSR experiment, the space-charge-induced frequency spread in the transverse particle orbits can be large, particularly for the electrons.During the linear growth phase of the instability, it is expected that this frequency spread accounts for the much lower growth rates of the e-p instability observed in the simulation studies and in the PSR experiment relative to those calculated from a centroid model or a kinetic KV model with a flattop density profile.

V. CONCLUSIONS
In conclusion, a 3D multispecies nonlinear perturbative particle simulation method has been developed to study the electron-ion two-stream instability described self-consistently by the Vlasov-Maxwell equations.Important properties of this instability were investigated numerically and are found to be in qualitative agreement with theoretical predictions [5][6][7] and the PSR experiment [12 -14].Numerically, the instability threshold was found to decrease with increasing fractional charge neutralization and increase with increasing axial momentum spread of the beam particles.In the nonlinear phase, the simulation results showed that the instability first saturates at a relatively low level and subsequently grows to a much larger level.The results suggest that due to the large mass ratio, the electron density perturbation quickly saturates long before the proton density perturbation becomes sizable, and the large electron density fluctuation level then provides a newly developed background force that drives the proton density perturbation to a large level on a longer time scale [29].Finally, the nonlinear spacecharge-induced transverse tune spread, which provides an important mechanism for reducing the growth rate of the instability, was analyzed for self-consistent equilibrium populations of the protons and the electrons.Investigations of the long-time nonlinear evolution of the e-p instability, together with the inclusion of electron production mechanisms and the effects of finite bunch length, are important extensions of the present simulation studies, carried out for a long, coasting beam.We will continue our investigations in these areas and report new results in future publications.
b 1:85, m e =m b 1=1836, !b 40 MHz, r w 5 cm, V e 0, and !e 0 (corresponding to axially stationary electrons).The system parameters for the coasting-beam ''baseline'' case in the present study are taken to be ! !2 pb =2 2 b ! 2 b 0:079, T b? = b m b V 2 b 3:61 10 ÿ6 , T e? = b m b V 2 b 5:86 10 ÿ7 , and fractional charge neutralization f n n e = n n b 0:1, corresponding to an average proton current of 35 A in the PSR experiment.Here, !! p 4 n n b e 2 b = b m b 1=2 FIG.1.Plots of the normalized equilibrium self-field potential profiles.

FIG. 2 .
FIG. 2. Maximum growth rate versus normalized beam density for different values of initial axial momentum spread of the beam ions.
are !!b r=! b and !! e r=! b versus r=r w for several values of normalized beam intensity n n b = n n b0 .Here n n b0 9:41 10 8 cm ÿ3 is the density for the baseline case corresponding to an average current of 35 A in the PSR experiment.The fractional charge neutralization f n n e = n n b is taken to be 10% in Fig. 5. Evidently, the local transverse oscillation frequency for the protons is a minimum at the beam center r 0 and decreases with increasing beam intensity.Because e 0 and !e 0 are assumed, the local transverse oscillation frequency of the electrons is due to the spacecharge force only, and maximizes at the beam center, and increases with increasing beam intensity.For the baseline case with n n b = n n b0 1:0, the maximum local tune depression for the protons is 1 ÿ !!b =! b 0:4% at r 0, where the local transverse oscillation frequency of the electrons also reaches its maximum value of !! e =! b 28:5.
are F b !!b and F e !! e for different values of normalized beam intensity.
FIG. 6. Plots of the distributions function F j !!j for the (a) protons and (b) electrons.