Three-dimensional multispecies nonlinear perturbative particle simulations of collective processes in intense particle beams

Collective processes in intense charged particle beams described self-consistently by the VlasovMaxwell equations are studied using a 3D multispecies nonlinear perturbative particle simulation method. The newly developed beam equilibrium, stability, and transport (BEST) code is used to simulate the nonlinear stability properties of intense beam propagation, surface eigenmodes in a high-intensity beam, and the electron-proton (e-p) two-stream instability observed in the Proton Storage Ring (PSR) experiment. Detailed simulations in a parameter regime characteristic of the PSR experiment show that the dipolemode two-stream instability is stabilized by a modest spread (about 0.1%) in axial momentum of the beam particles.


I. INTRODUCTION
Periodic focusing accelerators and transport systems [1 -4] have a wide range of applications ranging from basic scientific research to applications such as spallation neutron sources, heavy ion fusion, and nuclear waste transmutation.As 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 produced by the beam space charge and current on detailed equilibrium, stability, and transport properties.For a one-component highintensity beam, considerable progress has been made in describing the self-consistent evolution of the beam distribution function f b ͑x, p, t͒ and the self-generated electric and magnetic fields in kinetic analysis [5 -11] based on the Vlasov-Maxwell equations.In many practical accelerator applications, however, an (unwanted) second charge component is present.For example, a background population of electrons can result by secondary emission when energetic beam ions strike the chamber wall, or through ionization of background neutral gas by the beam ions.When a second charge component is present, it has been recognized for many years, both in theoretical studies and in experimental observations [12][13][14][15][16][17][18][19][20][21], that the relative streaming motion of the high-intensity beam particles through the background charge species provides the free energy to drive the classical two-stream instability [22], appropriately modified to include the effects of dc space charge, relativistic kinematics, presence of a conducting wall, etc.A well-documented example is the electron-proton (e-p) instability observed in the Proton Storage Ring (PSR) [17,18], although a similar instability also exists for other ion species including (for example) electron-ion interactions in electron storage rings [19 -21].
Recently, the df formalism, a low-noise, nonlinear perturbative particle simulation technique, has been developed for intense beam applications and applied to matched-beam propagation in a periodic focusing field [23,24] and other related studies.The present paper reports recent advances in applying the df formalism to investigate nonlinear collective processes in intense charged particle beams.The BEST ( beam equilibrium, stability, and transport) code [25] described here is a newly developed 3D multispecies nonlinear perturbative particle simulation code, which can be applied to a wide range of important collective processes in intense beams, such as the electron-ion two-stream instability [12][13][14][15][16][17][18] and the periodically focused beam propagation [11].While several features of the 3D BEST code are summarized in Ref. [25], together with initial simulation results, the present paper is the first comprehensive description of the simulation capability and its application to the nonlinear dynamics of the two-stream instability.The 3D and multispecies capability of the simulation code are required by the physics of the collective processes under investigation.In general, collective processes have 3D spatial structures, and a second or third species of charged particle is introduced into the system either intentionally or unintentionally.
Following a description of the nonlinear df formalism (Sec.II), this paper presents detailed simulation results (Sec.III) for the electron-ion two-stream instability, with particular emphasis on the parameter regime characteristic of the PSR experiment [17,18].Most importantly, the simulations show that the instability can be stabilized by a modest spread in axial momentum of the beam particles.

II. NONLINEAR df FORMALISM
The theoretical model employed here that incorporates collective effects is based on the nonlinear Vlasov-Maxwell equations.We consider a thin, continuous, highintensity ion beam ͑ j b͒, with characteristic radius r b propagating in the z direction through background electron and ion components ͑ j e, i͒, each of which is described by a distribution function f j ͑x, p, t͒ [10][11][12].The charge components ͑ j b, e, i͒ propagate in the z direction with characteristic axial momentum g j m j b j c, where V j b j c is the average directed axial velocity, g j ͑1 2 b 2 j ͒ 21͞2 is the relativistic mass factor, e j and m j are the charge and rest mass, respectively, of a jth species particle, and c is the speed of light in vacuo.While the nonlinear df formalism outlined here is readily adapted to the case of a periodic applied focusing field, for present purposes we make use of a smooth-focusing model in which the applied focusing force is described by where x Ќ x êx 1 y êy is the transverse displacement from the beam axis, and v bj const is the effective applied betatron frequency for transverse oscillations.Furthermore, in a frame of reference moving with axial velocity b j c, the motion of a jth species particle is assumed to be nonrelativistic.The space-charge intensity is allowed to be arbitrarily large, subject only to transverse confinement of the beam ions by the applied focusing force, and the background electrons are confined in the transverse plane by the space-charge potential f͑x, t͒ produced by the excess ion charge.In the electrostatic and magnetostatic approximation, we represent the selfelectric and self-magnetic fields as E s 2=f͑x, t͒ and B s = 3 ͓A z ͑x, t͒ê z ͔.
The nonlinear Vlasov-Maxwell equations in the six-dimensional phase space ͑x, p͒ can be approximated by [10][11][12]] and In the nonlinear df formalism, we divide the total distribution function into two parts, f j f j0 1 df j , where f j0 is a known solution to the nonlinear Vlasov-Maxwell equations [(2) and (3)], and the numerical simulation is carried out to determine only the detailed nonlinear evolution of the perturbed distribution function df j .This is accomplished by advancing the weight function defined by w j ϵ df 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 (2), are given by dx ji dt ͑g j m j ͒ 21 p ji , Here the subscript "ji" labels the ith simulation particle of the jth species.The weight functions w j , as functions of phase space variables, are carried by the simulation particles, and the dynamical equations for w j are easily derived from the definition of w j and the nonlinear Vlasov equation (2).Following the algebra in Refs.[23,24] , we obtain where df f 2 f 0 and dA z A z 2 A z0 .Here, the equilibrium solutions (f 0 , A z0 , f j0 ) solve the steady-state (≠͞≠t 0) Vlasov-Maxwell equations ( 2) and ( 3) with ≠͞≠z 0 and ≠͞≠u 0. A wide variety of axisymmetric equilibrium solutions to Eqs. ( 2) and ( 3) have been investigated in the literature.The perturbed distribution df 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 where Here, S͑x 2 x ji ͒ represents the method of distributing particles on the grids in configuration space.
The nonlinear particle simulations are carried out by iteratively advancing the particle motions, including the weights they carry, according to Eqs. ( 4) and ( 5), and updating the fields by solving the perturbed Maxwell's equations (7) with appropriate boundary conditions at the cylindrical, perfectly conducting wall.Even though it is a perturbative approach, the df method is fully nonlinear and simulates completely the original nonlinear 084401-2 084401-2 Vlasov-Maxwell equations.Compared with conventional particle-in-cell (PIC) simulations, the noise level in df simulations is significantly reduced.The dominant numerical noise mechanisms in particle simulations, such as numerical collisions, are statistical.The df method reduces the noise level of the simulations because the statistical noise, which is of order O͑N 21͞2 s ͒ for the total distribution function in the conventional PIC method, is associated only with the perturbed distribution function in the df method.If the same number of simulation particles is used in the two approaches, then the noise level in the df method is reduced by a factor of f͞df relative to the PIC method.To achieve the same accuracy for the perturbed fields, the number of simulation particles used in the df method is reduced by a factor of ͑f͞df͒ 2 .For the electron-proton two-stream instability in the Proton Storage Ring experiment studied in this paper, we obtain satisfactory results with about 10 5 simulation particles using the df method.If the conventional PIC method is used, for a nonlinear saturation level of 1%, about 10 4 times more simulation particles would be needed to achieve the same accuracy.In addition, the df method can be used to study linear stability properties, provided the factor ͑1 2 w ji ͒ in Eq. ( 5) is approximated by unity and the forcing term in Eq. ( 4) is replaced by the unperturbed force, which is equivalent to integrating along unperturbed particle orbits for the linearized system.
Implementation of the 3D multispecies nonlinear df simulation method described above is embodied in the BEST code [25] developed at the Princeton Plasma Physics Laboratory.The code advances the particle motions using a leapfrog method and solves Maxwell's equations in cylindrical geometry.For those fast particle motions which require much larger sampling frequency 1͞Dt than the frequency of the mode being studied, the code uses an adiabatic field pusher to advance the particles many time steps without solving for the perturbed fields.The upper limit for Dt, the time step to advance the particle's phase space position, is normally determined by the Courant condition.For the electron-proton two-stream instability, the electron's transverse motion requires the smallest Dt, and the mode frequency is comparable to the electron bounce frequency in the transverse direction.If the number of grids in the radial direction is 100, then the frequency for advancing the electron motion should be 100 times larger than the mode frequency.We can therefore update the electrons' phase space positions more often than the field.As in many particle simulation schemes, an implicit scheme was not chosen to be used due to the memory and CPU time considerations for large-scale simulations.
The code has achieved an average speed of 40 ms͞ ͑ particle 3 step͒ on a DEC alpha personal workstation 500au computer and has been parallelized using both the OpenMP and MPI parallel schemes.Written in Fortran 90/ 95, the code extensively utilizes the object-oriented features provided by the computer language.The NetCDF scientific data format is implemented for large-scale diagnostics and visualization.

III. SIMULATION RESULTS
We first present application of the code to a singlespecies thermal equilibrium ion beam in a constant focusing field.It is assumed that the beam is centered inside a cylindrical chamber with perfectly conducting wall located at r r w , and the equilibrium is one-dimensional, depending only on the radial coordinate r ͑x 2 1 y 2 ͒ 1͞2 .The isotropic thermal equilibrium distribution function in the phase space ͑r, p͒ is given by where nb is the number density of beam particles at r 0 and T b const is the temperature of the beam ions in energy units.The equilibrium self-field potentials f 0 and A z0 can be determined numerically from the nonlinear Maxwell's equations in Eq. (3).As an example, we examine the nonlinear propagation properties of a heavy ion beam with g b 1.08, mass number A 133, and normalized space-charge intensity v2 pb ͞2g 2 b v 2 bb 0.95.Here, v2 pb 4p n2 b e 2 b ͞g b m b is the relativistic plasma frequency squared on axis ͑r 0͒.A random threedimensional initial perturbation is introduced into the system, and the beam is propagated from t 0 to t 1200t b , where t b ϵ v 21  bb .The simulation results show that the perturbations do not grow and the beam propagates quiescently over large distance, which agrees with the nonlinear stability theorem [9] for the choice of monotonically decreasing equilibrium distribution function in Eq. (9).Shown in Fig. 1 is a plot of the density perturbation at one spatial location versus normalized time v bb t, for perturbations about the thermal equilibrium distribution in Eq. ( 9).The amplitudes of the initial random perturbation in weights in Fig. 1 are of order 10 25 , which leads to a very small density perturbation.It is evident from Fig. 1 that the perturbations remain extremely small, and the beam propagates quiescently over very large distances, as expected.
As a second example, we study the linear surface mode for perturbations about a thermal equilibrium ion beam in the space-charge-dominated regime, with flattop density profile.These modes are of practical interest because they can be destabilized by a two-stream electron-ion interaction when background electrons are present [12 -18].The BEST code, operating in its linear stability mode, has recovered very well-defined eigenmodes with mode structures and eigenfrequencies which agree well with theoretical predications [1,12].For the dipole mode with azimuthal mode number l 1, the dispersion relation for these modes is given by [12] v where r b is the radius of the beam edge and r w is the location of the conducting wall.In Eq. ( 10), v2 pb 4p nb e 2 b ͞g b m b is the ion plasma frequency-squared and vpb ͞ p 2 g b Ӎ v bb has been assumed in the spacecharge-dominated limit.Shown in Fig. 2(a) is a comparison between plots of the eigenfrequency versus r w ͞r b obtained from the simulations (diamonds and triangles) and that predicted by Eq. ( 10) (solid curves).The system parameters for this simulation are chosen close to the space-charge limit, and the perturbation has normalized axial wave number k z V b ͞v bb 2p.It is clear from Fig. 2 that the simulation results agree very well with theoretical predictions.
In a high-intensity ion beam, the surface mode described above can be destabilized by the presence of a background electron population [12][13][14][15][16][17][18].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 2 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 [17,18] is believed to have this two-stream characteristic.For the case of a flattop density profile for both ions and electrons, assuming V e 0, v be 0, and negligible spread in axial momentum, we obtain the analytical dispersion relation for dipole perturbations with azimuthal mode number l 1 [12], where f ϵ ne ͞ nb is the fractional charge neutralization and v 2 e and v 2 b are defined by Equation ( 11) supports an unstable eigenmode with eigenfrequency v closely tuned to the electron collective oscillation frequency v Ӎ v e and axial wave number k z V b Ӎ v e 1 v b .For a near-Gaussian density profile, there is no known analytical description of the eigenmode.However, the simulation results indicate that the basic characteristics of the unstable mode are qualitatively similar to the case of a flattop density profile.For example, the most unstable mode is a l 1 dipole mode and is localized in the region where the ion beam density gradient is largest.We present here simulation results for the electronproton two-stream instability with moderate space-charge .86 3 10 27 , and f ϵ ne ͞ nb 0.1, V e 0, and v be 0 (stationary electrons).These system parameters correspond to the typical operating parameters in the PSR experiment [17,18].To achieve a good theoretical understanding of the electron-proton twostream instability, it is necessary to use a fully 3D, kinetic, low-noise simulation method.This is because, first of all, the instability has a 3D mode structure which depends on (x, y, z).In the r direction, the mode peaks at the location with maximum equilibrium density gradient.In the u direction, the mode has an l 1 dipole structure.The mode structure varies in the z direction as well, because the resonance condition of the mode corresponds to a certain wavelength band in the z direction.Full kinetic effects should be included in the study because they dominate the stabilization process and the nonlinear saturation of the instability-questions of significant importance for experiments and therefore for theoretical studies.Because of the large mass ratio between the ions and the electrons, and the fact that the growth rate of the instability is much smaller than the real frequency of the eigenmode, it takes a relatively long time to simulate the instability.A low-noise method, such as the df method used here, is therefore highly desirable.In the simulations of the e-p instability, we take the background distribution functions f j0 ͑r, p͒ to be the bi-Maxwellian generalization of Eq. ( 9), with temperature T jЌ const in the x-y plane and temperature T jk const in the z direction.Because the e-p instability is strongest when the beam ions are cold in the parallel direction [12] (no Landau damping by parallel kinetic effects), we take T bk 0 and T ek 0 in the simulations presented in Figs.3-7 effects at increasing values of T bk ͞T bЌ is illustrated in Fig. 8. Illustrated in Fig. 3 is a typical unstable case, where the x-y projection (at fixed value of z) of the perturbed spacecharge potential df͑x, y, t͒ grows exponentially with time during the linear phase of the instability.Clearly, the unstable mode is a dipole mode with azimuthal mode number l 1 [12].It is important to emphasize that the simulations are based on first principles-the nonlinear Vlasov-Maxwell equations.All possible mode excitations are allowed in the simulations.Simulations using typical operating parameters in the PSR experiment [17,18] indicate that the l 1 dipole mode is the most unstable mode.For this dominant mode in Fig. 3, the real part of the eigenfrequency is Rev 25.13v bb , and the normalized wavelength in the longitudinal direction is k z V b ͞v bb 26.17.These results are in good agreement with those measured in the PSR experiments [17,18].
Plots of the instability growth rate g Imv versus T eЌ ͞T bЌ and k z V b ͞v bb , with other parameters kept constant, are shown in Fig. 4. The k z V b ͞v bb dependence of the growth rate is qualitatively consistent with the analytical results obtained for uniform-density beams [12].The important physics here is that only for a certain range of k z V b ͞v bb can the collective mode of the beam ions effectively resonate with the electrons and produce instability.From Fig. 4 has an important effect on the growth rate.In order to maximize their energy exchange with the beam ions, the electrons must spatially overlap the region where the eigenmode of the beam ions is localized (approximately the region with the largest transverse gradient in ion density), which requires sufficiently large T eЌ ͞T bЌ .The electrons are radially confined by the space-charge potential of the beam ions, and the perpendicular electron temperature determines the radial extent of the electron density profile.The growth rate is therefore strongly dependent on T eЌ ͞T bЌ .Shown in Figs. 5 and 6 is a comparison between two cases with identical parameters except for the values of T eЌ ͞T bЌ .When T eЌ ͞T bЌ 0.018 (Fig. 5), the electrons are relatively cold and localized at the beam center and no instability develops over 370v 21  bb .When T eЌ ͞T bЌ 0.130 (Fig. 6), however, the electrons are sufficiently hot that the electron density profile overlaps that of the beam ions and the onset of a strong e-p instability is observed.
Finally, for T eЌ ͞T bЌ 0.130, the simulation results for the linear and nonlinear phases of the instability are shown in Fig. 7. nonlinearly saturates at t ϳ 400v 21 bb at a normalized amplitude of dn b ͞n b ϳ 0.3%.
In the simulation results for the e-p instability presented above, we have assumed cold beam ions in the longitudinal direction (T bk 0) to maximize the growth rate of the instability.In general, when the longitudinal temperature of the beam ions is finite, Landau damping by parallel ion kinetic effects provides a mechanism that reduces the growth rate.The decrease in the linear growth rate due to Landau damping of the unstable modes can be estimated to be FIG.8.The maximum linear growth rate ͑Imv͒ max of the electron-proton instability decreases as the longitudinal temperature of the beam ions increases.of order k z y Tbk , where y Tbk ͑2T bk ͞g b m b ͒ 1͞2 .Shown in Fig. 8 is a plot of the maximum linear growth rate ͑Imv͒ max versus T bk ͞T bЌ and k z y Tbk obtained in numerical simulations of the e-p instability using the BEST code.As is evident from the figure, the growth rate decreases dramatically as T bk ͞T bЌ and k z y Tbk increase.When T bk is high enough that k z y Tbk is comparable to the linear growth rate for the T bk 0 case, the mode is stabilized by longitudinal Landau damping by the beam ions.Because the phase velocity of the mode in the longitudinal direction is far removed from the electron velocity distribution j v͞k z j ¿ V e 1 y Tek , we do not expect the longitudinal electron temperature to affect significantly the growth rate of the instability.

IV. CONCLUSIONS
In conclusion, a 3D multispecies nonlinear perturbative particle simulation method has been developed to study collective processes in intense charged particle beams described self-consistently by the Vlasov-Maxwell equations.The simulation results show that an isotropic thermal equilibrium ion beam in a constant focusing field is nonlinearly stable and can propagate quiescently over hundreds of lattice periods.For the surface eigenmodes excited in a uniform-density beam, the simulation results agree well with analytical results [12].Introducing a background component of electrons, a strong electronproton two-stream instability is observed in the simulations when T bk is sufficiently small.Several properties of this instability have been investigated numerically and are found to be in qualitative agreement with theoretical predictions [12].Simulation results for the typical operating parameters in the PSR experiment [17,18] are in good agreement with the experimental results in terms of the dipole mode structure, wavelength, and eigenfrequency.Most importantly, in the simulations, a small spread in axial momentum of the beam particles is found to be effective in stabilizing the two-stream instability.The newly developed BEST code has been tested and applied in different parameter regimes.As a 3D multispecies perturbative particle simulation code, it provides several unique capabilities.Since the simulation particles are used to simulate only the perturbed distribution functions and the perturbed self-fields, the simulation noise is reduced significantly.The perturbative approach also enables the code to investigate different physics effects separately as well as simultaneously.The code can be easily switched between linear and nonlinear operation and used to study both linear stability properties and nonlinear beam dynamics.These features, combined with 3D and multispecies capabilities, provide an effective tool to investigate the electron-ion two-stream instability, periodically focused solutions in alternating gradient field configurations, halo formation, and many other important problems in nonlinear beam dynamics and accelerator physics.Finally, the BEST code is readily adapted to the case where the applied focusing force F foc j corresponds to a periodic focusing quadrupole field or solenoidal field.Results of these studies will be reported in future publications.

FIG. 1 .
FIG. 1.Time history of dn b ͞ nb for small-amplitude perturbations about a thermal equilibrium ion beam.

FIG. 2 .
FIG. 2. Excitation of the l 1 surface mode in a uniformdensity ion beam, including (a) a plot of the normalized oscillation frequency v͞v bb versus r w ͞r b and (b) a plot of the frequency spectrum for r w ͞r b 2.2.

v2 pb ͞2g 2 b v 2
bb 0.074, g b 1.85, and m e ͞m b 1͞1836.The equilibrium distribution functions f j0 are chosen to be thermal equilibrium distributions for both species with T bЌ ͞g b m b V 2 b 3.61 3 10 26 , T eЌ ͞g b m b V 2 b 5

FIG. 3 .
FIG.3.The x-y projection (at fixed value of z) of the perturbed electrostatic potential df͑x, y, t͒ for the electron-proton two-stream instability growing from a small initial perturbation, shown at (a) t 0 and ( b) v bb t 200.
FIG. 4. Plots of the linear growth rate g Imv, including (a) g͞v bb versus k z V b ͞v bb and (b) g͞v bb versus T eЌ ͞T bЌ .
FIG.5.Plots of (a) equilibrium proton and electron density profiles and (b) perturbation time history for cold electrons with T eЌ ͞T bЌ 0.018.

FIG. 7 .
FIG.6.Plots of (a) equilibrium proton and electron density profiles and (b) perturbation time history for warm electrons with T eЌ ͞T bЌ 0.130.