Kinetic description of electron-proton instability in high-intensity proton linacs and storage rings based on the Vlasov-Maxwell equations

The present analysis makes use of the Vlasov-Maxwell equations to develop a fully kinetic descrip of the electrostatic, electron-ion two-stream instability driven by the directed axial motion of a hig intensity ion beam propagating in the z direction with average axial momentumgbmbbbc through a stationary population of background electrons. The ion beam has characteristic radius rb and is treated as continuous in the z direction, and the applied transverse focusing force on the beam ions modeled byF foc ­ 2gbmbv 02 bbx' in the smooth-focusing approximation. Here, v 0 bb ­ const is the effective betatron frequency associated with the applied focusing field, x' is the transverse displacement from the beam axis,sgb 2 1dmbc is the ion kinetic energy, andVb ­ bbc is the average axial velocity, wheregb ­ s1 2 b bd. Furthermore, the ion motion in the beam frame is assumed to b nonrelativistic, and the electron motion in the laboratory frame is assumed to be nonrelativistic. T ion charge and number density are denoted by 1Zbe and nb , and the electron charge and number density by2e and ne. For Zbnb . ne, the electrons are electrostatically confined in the transvers direction by the space-charge potential f produced by the excess ion charge. The equilibrium and stability analysis retains the effects of finite radial geometry transverse to the beam propaga direction, including the presence of a perfectly conducting cylindrical wall located at radius r ­ rw. In addition, the analysis assumes perturbations with long axial wavelength, k2 z r 2 b ø 1, and sufficiently high frequency that jvykz j ¿ yTez andjvykz 2 Vb j ¿ yTbz , whereyTez andyTbz are the characteristic axial thermal speeds of the background electrons and beam ions. In this regime, Landau dampin axial velocity spaceyz) by resonant ions and electrons is negligibly small. We introduce the ion plasm frequency squared defined bŷ v pb ­ 4pn̂bZ 2 beygbmb, and the fractional charge neutralization defined by f ­ n̂eyZbn̂b , wheren̂b and n̂e are the characteristic ion and electron densities. The equilibrium and stability analysis is carried out for arbitrary normalized beam intensitŷ v pbyv 2 bb, and arbitrary fractional charge neutralization f, consistent with radial confinement of the beam particles. For th moderately high beam intensities envisioned in the proton linacs and storage rings for the Acceler for Production of Tritium and the Spallation Neutron Source, the normalized beam intensity is typica v̂ pbyv 2 bb & 0.1. For heavy ion fusion applications, however, the transverse beam emittance is v small, and the space-charge-dominated beam intensity is much larger, with v̂ pbyv 2 bb & 2g 2 b. The stability analysis shows that the instability growth rate Im v increases with increasing normalized beam intensity v̂ pbyv 2 bb and increasing fractional charge neutralization f. In addition, the instability is strongest (largest growth rate) for perturbations with azimuthal mode number , ­ 1, corresponding to a simple (dipole) transverse displacement of the beam ions and the background electrons. For the ca overlapping step-function density profiles for the beam ions and background electrons, correspondin monoenergetic ions and electrons, a key result is that there is no threshold in beam intensity v̂ pbyv 2


I. INTRODUCTION
Periodic focusing accelerators and transport systems [1][2][3][4] have a wide range of applications ranging from high-intensity 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 E s ͑x, t͒ and B s ͑x, t͒ in kinetic analyses [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] based on the nonlinear Vlasov-Maxwell equations.For example, in a recent calculation [23,24], a three-dimensional, kinetic stability theorem based on the Vlasov-Maxwell equations has been developed for a high-intensity ion beam (or charge bunch) in the smoothfocusing approximation.It is found that a beam equilibrium f 0 b ͑x 0 , p 0 ͒ that is a monotonically decreasing function of total particle energy H 0 b in the beam frame is nonlinearly stable to perturbations with arbitrary amplitude and polarization.The analysis [23,24] is valid for arbitrary beam intensity consistent with transverse confinement of the beam particles by the focusing field and includes the effects of a perfectly conducting cylindrical wall located at radius r r w .
In many practical accelerator applications, however, an (unwanted) second charge component is present.For example, a background population of electrons can result locally when an H 2 beam is injected through a stripper foil into a proton storage ring or when energetic ions strike the chamber wall.When a second charge component is present, it has been recognized for many years, both in theoretical studies and in experimental observations [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42], 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 [43 -45], appropriately modified to include the effects of dc space charge, relativistic kinematics, presence of a conducting wall, etc.For electrons interacting with a proton beam, as in the Proton Storage Ring (PSR), this instability is usually referred to as the electron-proton ͑e-p͒ instability [31][32][33][34][35], although a similar instability also exists for other ion species, including (for example) electron-ion interactions in electron storage rings [36][37][38][39][40][41].Moreover, a related instability (known as the "ion-resonance" instability), driven by the relative average motion of ion and electron components, also exists in electron-rich non-neutral plasmas [26,29,30] and in collective acceleration schemes such as the electron ring accelerator [42].
Theoretical treatments of the e-p instability are traditionally based on models [46] that analyze the center-ofmass motion of the ion and electron charge components.Such models, while treating accurately several bulk features of the instability, are limited in scope and difficult to generalize to include the dependence of stability behavior on the detailed phase-space properties of the ion and electron distribution functions.Therefore, in the present analysis, we develop and apply a theoretical formalism based on the Vlasov-Maxwell equations [1,47] that describe the self-consistent interaction of the ion and electron distribution functions, f b ͑x, p, t͒ and f e ͑x, p, t͒, with the applied field and the self-generated electric and mag-netic fields.Furthermore, in integrating the linearized Vlasov-Maxwell equations, we make use of the method of characteristics [44,47] to integrate along the particle trajectories in the equilibrium field configuration.Such an approach has proven to be a powerful technique for describing the stability properties of spatially nonuniform, non-neutral plasmas and intense beam systems [1,47], as well as spatially nonuniform, electrically neutral plasmas [44,48].As noted below, the present analysis does include the (stabilizing) influence of a perfectly conducting cylindrical wall located at radius r r w const.However, as a consequence, the analysis does not include other important instabilities, such as the resistive-wall instability [49] or the beam breakup instability [50], that can result from finite wall resistivity or the presence of structures attached to the chamber wall.
To briefly summarize, the present analysis makes use of the Vlasov-Maxwell equations to develop a fully kinetic description of the electrostatic, electron-ion two-stream instability driven by the directed axial motion of a highintensity ion beam propagating in the z direction with average axial momentum g b m b b b c through a stationary population of background electrons.The ion beam has characteristic radius r b and is treated as continuous in the z direction, and the applied transverse focusing force on the beam ions is modeled by F b foc 2g b m b v 0 2 bb x Ќ in the smooth-focusing approximation.Here, v 0 bb const is the effective betatron frequency associated with the applied focusing field, x Ќ is the transverse displacement from the beam axis, ͑g b 2 1͒m b c 2 is the ion kinetic energy, and V b b b c is the average axial velocity, where g b ͑1 2 b 2 b ͒ 21͞2 .Furthermore, the ion motion in the beam frame is assumed to be nonrelativistic, and the electron motion in the laboratory frame is assumed to be nonrelativistic.The analysis generally allows for an applied transverse focusing force on the electrons modeled by F e foc 2m e v 0 2 be x Ќ , where v 0 be const.We denote the ion charge and number density by 1Z b e and n b , and the electron charge and number density by 2e and n e .For Z b n b .n e , the electrons can be electrostatically confined in the transverse direction by the space-charge potential f produced by the excess ion charge, even when v 0 be 0. The present equilibrium and stability analysis retains the effects of finite radial geometry transverse to the beam propagation direction, including the presence of a perfectly conducting cylindrical wall located at radius r r w .In the stability analysis, the z and t dependences of perturbed quantities are assumed to be of the form exp͑ik z z 2 ivt͒, where k z is the axial wave number and v is the complex oscillation frequency, with Imv .0 corresponding to instability (temporal growth).The present analysis assumes perturbations with long axial wavelength, k 2 z r 2 b ø 1, and sufficiently high frequency that jv͞k z j ¿ y Tez and jv͞k z 2 V b j ¿ y Tbz , where y Tez and y Tbz are the characteristic axial thermal speeds of the background electrons and beam ions.In 054401-2 054401-2 PRST-AB 2 KINETIC DESCRIPTION OF ELECTRON-PROTON … 054401 (1999) this regime, Landau damping (in axial velocity space y z ) by resonant ions and electrons is negligibly small.For step-function density profiles, we introduce the ion plasma frequency-squared defined by v2 pb 4p nb Z 2 b e 2 ͞g b m b , and the fractional charge neutralization defined by f ne ͞Z b nb , where nb and ne are the ion and electron densities.The present analysis based on the Vlasov-Maxwell equations is carried out for arbitrarily normalized beam intensity v2 pb ͞v 0 2 bb , and arbitrary fractional charge neutralization f, consistent with radial confinement of the beam particles.For the moderately high beam intensities envisioned in the proton linacs and storage rings for the spallation neutron sources and tritium production, the normalized beam intensity is typically v2 pb ͞v 0 2 bb & 0.1 [51].For heavy ion fusion applications [6,7], however, the transverse beam emittance is very small, and the space-charge-dominated beam intensity is much larger, with v2 pb ͞v 0 2 bb & 2g 2 b .The present stability analysis shows that the instability growth rate Imv increases with increasing normalized beam intensity v2 pb ͞v 0 2 bb , and increasing fractional charge neutralization f.In addition, the instability is strongest (largest growth rate) for perturbations with azimuthal mode number ᐉ 1, corresponding to a simple (dipole) transverse displacement of the beam ions and the background electrons.
The organization of this paper is the following.In Sec.II, we summarize the basic assumptions (Sec.II A) and describe the theoretical model based on the Vlasov-Maxwell equations (Sec.II B).Examples of self-consistent equilibrium solutions ͑≠͞≠t 0͒ to the Vlasov-Maxwell equations are then presented for the case of an intense, continuous ion beam propagating through a stationary background population of electrons (Sec.II C).In Sec.III, we formally integrate the linearized Vlasov-Maxwell equations using the method of characteristics (Sec.III A) and discuss properties of the ion and electron orbits in the applied field plus equilibrium self-field configuration.The orbit equations are analyzed both for the case of step-function ion and electron density profiles (Sec.III B), corresponding to monoenergetic beam ions and monoenergetic electrons, and for the case where the equilibrium density profiles have a continuous variation with radius r, corresponding to a spread in (depressed) betatron frequencies (Sec.III C).In Sec.IV, the necessary orbit integrals are evaluated in closed analytical form for the case of step-function ion and electron density profiles, leading to a kinetic dispersion relation which is valid for arbitrary normalized beam intensity v2 pb ͞v 0 2 bb , fractional charge neutralization f, and azimuthal mode number ᐉ (Sec.IV A).The resulting dispersion relation is analyzed in detail for the case of mode number ᐉ 1, corresponding to a simple transverse displacement of the beam ions and electrons (Sec.IV B), and a brief discussion of stability behavior for quadrupole perturbations with mode number ᐉ 2 is presented (Sec.IV C).For monoenergetic ions and electrons and the corresponding step-function density profiles considered in Sec.IV, a key result is that there is no threshold in beam intensity v2 pb ͞v 0 2 bb or fractional charge neutralization f for the onset of instability.Finally, for the case of continuously varying density profiles with parabolic profile shape, we make a semiquantitative estimate in Sec.V of the effects of the corresponding spread in (depressed) betatron frequency on stability behavior, including an estimate of the instability threshold for the case of weak density nonuniformity.

II. THEORETICAL MODEL AND ASSUMPTIONS
In this section, we summarize the basic assumptions (Sec.II A) made in the present analysis and describe the theoretical model based on the Vlasov-Maxwell equations (Sec.II B).Finally, examples of self-consistent equilibrium solutions ͑≠͞≠t 0͒ of the Vlasov-Maxwell equations are presented (Sec.II C) for the case of an intense ion beam propagating through a background population of electrons.

A. Basic assumptions
We consider a thin, high-intensity ion beam with distribution function f b ͑x, p, t͒, characteristic radius r b , and axial momentum g b m b b b c propagating in the z direction through a background population of electrons with distribution function f e ͑x, p, t͒.While the ions have high directed axial velocity V b b b c in the z direction, the background electrons are assumed to be nonrelativistic and stationary with R d 3 p p z f e Ӎ 0 in the laboratory frame.In the context of the smooth-beam approximation, the ion beam is assumed to be continuous in the z direction, and the applied transverse focusing force on a beam ion is modeled by where x Ќ x êx 1 y êy is the transverse displacement from the beam axis, ͑g b 2 1͒m b c 2 is the characteristic ion kinetic energy, m b is the ion rest mass, c is the speed of light in vacuo, and v 0 bb const is the effective betatron frequency for transverse ion motion in the applied focusing field.The focusing force in Eq. ( 1) would correspond to the transverse electric force produced by a uniformly distributed, fixed charge background with charge density r foc 2g b m b v 0 2 bb ͞2pZ b e const, where 1Z b e is the charge of a beam ion.Such a model is often used to describe the average focusing properties of an alternating-gradient lattice of magnetic or electric quadrupoles.For the background electrons, to the extent that the beam ion density exceeds the background electron density, the space-charge force on an electron, F s e e=f, provides transverse confinement of the background electrons by the electrostatic potential f͑x, t͒.
However, for completeness, the present analysis also incorporates the effects of an applied transverse focusing force on the electrons modeled by F e foc 2m e v 0 2 be x Ќ , where m e is the electron rest mass and v 0 be const is the effective betatron frequency for transverse electron motion in the applied focusing field.
It is further assumed that the ion motion in the beam frame is nonrelativistic and that the transverse momentum components of a beam ion, p x and p y , and the characteristic spread in axial momentum, dp z p z 2 g b mb b c, are small compared with the directed axial momentum; i.e., While the space-charge intensity in the present analysis is allowed to be arbitrarily large, subject only to transverse confinement of the beam ions by the focusing force in Eq. ( 1), it is assumed that where n B is Budker's parameter, N b ϵ R dx dy n b is the number of beam ions per unit axial length, and n b ͑x, t͒ R d 3 p f b ͑x, p, t͒ is the number density of the beam ions.Equation (3) assures that the space-charge potential is sufficiently weak that jZ b ef͞g b m b c 2 j ø 1; however, the electrostatic potential energy of a beam ion, Z b ef, is allowed to be comparable to or larger than the transverse kinetic energy ͑p 2 x 1 p 2 y ͒͞2g b m b of a beam ion.In addition, the present analysis is carried out in the electrostatic approximation, where the self-generated electric field produced by space-charge effects is and the electrostatic potential f͑x, y, z, t͒ is determined self-consistently from Poisson's equation Here, n e ͑x, t͒ R d 3 p f e ͑x, p, t͒ is the electron number density.Furthermore, to determine the self-generated magnetic field produced by the axial ion current, it is assumed that the axial velocity profile V zb ͑x, t͒ Ӎ b b c is approximately uniform over the beam cross section.In this case, in the magnetostatic approximation, the z component of vector potential A z ͑x, y, z, t͒ is determined self-consistently from where use has been made of the assumption that the electrons carry zero axial current in the laboratory frame, i.e., n e V ze R d 3 p͑p z ͞m e ͒f e Ӎ 0. Finally, under equilibrium conditions ͑≠͞≠t 0͒, the present analysis assumes that ion and electron properties are spatially uniform in the z direction with ≠͞≠z 0. However, in the stability analysis (Secs.III and IV), we assume small-amplitude perturbations with z and t variations proportional to exp͑ik z z 2 ivt͒, where k z 2pn͞L is the axial wave number and v is the (complex) oscillation frequency, with Imv .0 corresponding to instability.Here, n is an integer, L is the fundamental axial periodicity length of perturbed quantities in straight (e.g., linac) geometry, and L 2pR for the case of a storage ring with (large) radius R ¿ r b .As noted earlier, the electron motion is assumed to be nonrelativistic and the axial momentum spread of the ions is assumed to be small [see Eq. ( 2)].For our purposes here, the present stability analysis assumes electrostatic perturbations with sufficiently long axial wavelength l z 2p͞k z and sufficiently high frequency v that Here, y Tbz ͑2T bz ͞g b m b ͒ 1͞2 and y Tez ͑2T ez ͞m e ͒ 1͞2 are the characteristic axial thermal speeds of the beam ions and the background electrons, respectively.The inequalities in Eq. ( 8) lead to several simplifications in the Vlasov-Maxwell equations (Sec.II B).For example, because 5) and (7) for f͑x, t͒ and A z ͑x, t͒ can be approximated by the perpendicular Laplacian = 2 Ќ ≠ 2 ͞≠x 2 1 ≠ 2 ͞≠y 2 .Furthermore, because of the inequalities in Eq. ( 8), the perturbed axial forces on the electrons and ions [e.g., dF e e͑≠͞≠z͒df êz and dF b 2Z b e͑≠͞≠z͒df êz ] are treated as negligibly small.The subsequent analysis therefore neglects the effects of Landau damping (by resonant particles) due to the axial momentum spread [44] of the beam ions and background electrons.

B. Model based on nonlinear Vlasov-Maxwell equations
We now make use of the assumptions delineated in Sec.II A to simplify the theoretical model of ion beam interaction with the background electrons based on the Vlasov-Maxwell equations [1,47].First, for narrow axial momentum spread, we introduce the reduced distribution functions F b ͑x, p Ќ , t͒ F b ͑x, y, z, p x , p y , t͒ and F e ͑x, p Ќ , t͒ F e ͑x, y, z, p x , p y , t͒ defined by where integrations are over axial momentum p z .Because R dp z p z f e Ӎ 0 for the electrons, and axial forces are treated as negligibly small, the nonlinear Here, 2e is the electron charge, = Ќ ϵ êx ≠͞≠x 1 êy ≠͞≠y is the perpendicular gradient, and use has been made of the assumptions enumerated in Sec.II A. The ions, however, have large directed axial velocity V b Ӎ b b c.We therefore approximate v ?≠͞≠x Ӎ ͑p Ќ ͞g b m b ͒ ?≠͞≠x Ќ 1 V b ≠͞≠z, in the Vlasov equation for the reduced ion distribution function F b ͑x, p Ќ , t͒, and the perpendicular self-field force on an ion is approximated by F bЌ Z b e͓2= Ќ f 1 b b êz 3 ͑= Ќ A z 3 êz ͔͒, where f and A z are determined self-consistently from Eqs. ( 5) and (7).Consistent with the assumptions in Sec.II A, the Vlasov equation for the reduced ion distribution function F b ͑x, p Ќ , t͒ therefore becomes Here, 1Z b e is the ion charge and c͑x, t͒ is the combined potential defined by where the electrostatic potential f͑x, t͒ and the combined potential c͑x, t͒ solve In obtaining Eqs. ( 13) and ( 14), use has been made of Eqs. ( 6), (7), and (12); and n e ͑x, t͒ R d 2 p F e ͑x, p Ќ , z͒ are the ion and electron number densities, respectively; and we have approximated Ќ ≠ 2 ͞≠x 2 1 ≠ 2 ͞≠y 2 by virtue of the thinbeam approximation and the inequality k 2 z r 2 b ø 1 assumed in Eq. (8).
Equations (10), (11), (13), and ( 14) constitute a complete nonlinear description of the collective interaction of the beam ions with the background electrons based on the nonlinear Vlasov-Maxwell equations, consistent with assumptions enumerated in Sec.II A. In the subsequent analysis, we further assume that the ion beam propagates axially through a perfectly conducting cylindrical pipe with radius r r w , where r ͑x 2 1 y 2 ͒ 1͞2 is the radial distance from the beam axis.Enforcing ͓E s u ͔ rr w ͓E s z ͔ rr w ͓B s r ͔ rr w 0 at the conducting wall readily gives the boundary conditions f͑r r w , u, z, t͒ 0, c͑r r w , u, z, t͒ 0 .
Here, we have introduced cylindrical polar coordinates x r cosu and y r sinu, and the constant values of the potentials f and c at r r w have been taken equal to zero without loss of generality.
in Eqs.(10), (11), (13), and ( 14), we readily conclude that equilibrium distribution functions ͑≠͞≠t 0͒ for the beam ions and background electrons (denoted by F 0 b and F 0 e ) of the general form [47] exactly solve the nonlinear Vlasov-Maxwell equations, where H Ќb and H Ќe are the single-particle Hamiltonians defined by Here, for ≠͞≠u 0 ≠͞≠z, H Ќb and H Ќe are exact single-particle constants of the motion, and the constants ĉ0 ϵ c 0 ͑r 0͒ and f0 ϵ f 0 ͑r 0͒ are the on-axis ͑r 0͒ values of c 0 ͑r͒ and f 0 ͑r͒.
There is clearly very wide latitude in specifying the functional forms of the equilibrium distribution functions F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒.Once F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒ are specified, however, the equilibrium self-field potentials c 0 ͑r͒ and f 0 ͑r͒ can be calculated self-consistently from Eqs. ( 13) and ( 14) with ≠͞≠u 0 ≠͞≠z, i.e., Here, the equilibrium ion and electron density profiles, n 0 b ͑r͒ and n 0 e ͑r͒, are defined by A simple class of equilibrium distribution functions F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒ [8,9], which correspond to overlapping step-function density profiles for the beam ions and background electrons, is given by [30,47] Here, nb and ne ϵ fZ b nb are positive constants corresponding to the ion and electron densities, f const is the fractional charge neutralization, and TЌb and TЌe are constants corresponding to the on-axis ͑r 0͒ values of the transverse ion and electron temperatures, respectively.Without presenting details, some algebraic manipulation that makes use of Eqs. ( 19)-( 21) gives the step-function density profiles [30,47] and For overlapping density profiles with r e r b (Fig. 1), the equilibrium potential profiles calculated from Eqs. ( 19), (22), and ( 23) are given by e ͑r͒ [Eqs.(22) and (23)] for the choice of ion and electron distribution functions in Eq. (21).For the purpose of illustration, we have taken Z b 1 and f 0.5 in the figure . is related to TЌb and TЌe and other system parameters by the equilibrium constraint conditions " In Eq. ( 26), we have introduced the ion plasma frequencysquared defined by where N b nb pr 2 b .Note from Eqs. ( 24) and ( 25) that the constants f0 ϵ f 0 ͑r 0͒ and ĉ0 ϵ c 0 ͑r 0͒ are readily determined by enforcing the boundary conditions in Eq. ( 15 a wide range of equilibrium system parameters.This is illustrated in Fig. 2 for the case where v 0 be 0, corresponding to zero applied focusing force on the electrons.In this case, the electrons are radially confined by the electrostatic potential of the beam ions.Shown in Fig. 2 is the allowed region in ͑f, v2 pb ͞v 0 2 bb ͒ parameter space corresponding to radial confinement of the ions and electrons consistent with Eq. (26).Here, the solid bounding curves in Fig. 2  pb .The specific choices of distribution functions in Eq. ( 21) lead to the particularly simple forms for the equilibrium density and potential profiles in Eqs. ( 22)- (25).Most notably, the particle trajectories in the equilibrium field configuration consistent with Eqs. ( 22)-( 25) can be calculated in closed analytical form, allowing detailed stability properties to be determined from the linearized Vlasov-Maxwell equations (Secs.III and IV).
There are clearly many possible choices for the equilibrium distribution functions F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒.Another example would correspond to the thermal equilibrium distributions [47] Here, H Ќb and H Ќe are defined in Eq. ( 18), T Ќb and T Ќe are positive constants corresponding to ion and electron temperatures (energy units), and nb n 0 b ͑r 0͒ and ne n 0 e ͑r 0͒ are the on-axis ion and electron densities.Substituting Eq. ( 28) into Eq.( 20) gives for the equilibrium density profiles , When Eq. ( 29) is substituted into Eq.( 19), we note that the resulting coupled equations for the potentials f 0 ͑r͒ and f 0 ͑r͒ are highly nonlinear and must generally be solved numerically.Requiring radial confinement of the ions and electrons with n 0 b ͑r !`͒ 0 and n 0 e ͑r !`͒ 0 generally imposes restrictions on the allowed range of system parameters v 0 2 bb , v2 pb 4p nb Z 2 b e 2 ͞g b m b , ne ͞Z b nb , etc.For example, in the absence of electrons ͑ ne 0͒, careful examination FIG. 2. The area to the right of the solid curve is the allowed region of ͑f, v2 pb ͞v 0 2 bb ͒ parameter space corresponding to radial confinement of the ions and electrons consistent with Eq. ( 26) and v 0 be 0, TЌe $ 0, and TЌb $ 0. For the purpose of illustration, g b p 2 1.414 is assumed in the figure . of Eqs. ( 19) and (29) shows that the ion density profile n 0 b ͑r͒ is bell shaped [20], assuming a maximum value nb at r 0, and decreasing monotonically with increasing r, provided the inequality v2 pb ͞2g 2 b , v 0 2 bb is satisfied.This is simply a statement that (repulsive) space-charge forces must be weaker than the (applied) transverse focusing force.Typical numerical solutions to Eqs. (19) and (29) are illustrated in Fig. 3, where n 0 b ͑r͒ and n 0 e ͑r͒ are plotted versus radius r for the choice of system parameters v 0 be 0, g b 1.85, v2 pb ͞v 0 2 bb 0.1, and on axis fractional charge neutralization f ne ͞Z b nb 0.2.

III. LINEARIZED VLASOV-MAXWELL EQUATIONS
In preparation for the stability analysis in Sec.IV, we formally integrate the linearized Vlasov-Maxwell equations using the method of characteristics (Sec.III A) and discuss properties of the ion and electron orbits in the applied field plus equilibrium self-field configuration (Secs.III B and III C).The orbit equations are analyzed both for the case of the step-function density profiles and corresponding potential profiles in Eqs. ( 22)-( 25) (Sec.III B) and for the case where the equilibrium density profiles have a continuous variation with radius r corresponding to a spread in (depressed) betatron frequencies.
Equations ( 40) -(43) represent the final system of eigenvalue equations derived for small-amplitude perturbations about general equilibrium distribution functions F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒, consistent with the kinetic model and assumptions summarized in Sec.II.Equations ( 40)-(43) have a wide range of applicability and can be used to determine the complex oscillation frequency v and detailed stability properties for a wide range of system parameters and choices of distribution functions F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒.The principal challenge in analyzing Eqs. ( 40)-( 43) is twofold.First, depending on the equilibrium profiles, the transverse orbits ͑r 0 , u 0 ͒ or ͑x 0 , y 0 ͒ are often difficult to calculate in closed analytical form 054401-9 054401-9 PRST-AB 2 RONALD C. DAVIDSON et al. 054401 (1999) from Eqs. ( 37) and (38).Second, once the orbits in the equilibrium fields are determined, the integrations over t 0 in Eqs. ( 40) and ( 41) are challenging because the r 0 orbits occur explicitly in the arguments of the (yet unknown) eigenfunction amplitudes d ĉᐉ ͑r 0 ͒ and d fᐉ ͑r 0 ͒.For future reference in Secs.III B, III C, and IV, we express the ion and electron orbit Eqs.(37) and (38) in the convenient forms and Here, n 2 b ͑r͒ and n 2 e ͑r͒ are the (depressed) betatron frequencies-squared including applied-field plus self-field effects, defined by for the ions, and for the electrons.For specified equilibrium distribution functions F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒, and corresponding equilibrium density profiles n 0 b ͑r͒ and n 0 e ͑r͒, the equilibrium self-field contributions in Eqs. ( 46) and ( 47) can be expressed as and Here, use has been made of Eq. ( 19).Whenever the ion charge density exceeds the electron charge density with Z b n 0 b ͑r͒ .n 0 e ͑r͒, we note from Eqs. ( 45), (47), and (49) that the equilibrium self-field force on the electrons is always focusing, even when the applied betatron frequency v 0 be 0.

B. Particle orbits for step-function density profiles
The orbit equations ( 44) and ( 45) simplify considerably for the case of the step-function density profiles in Eqs. ( 22) and (23), which correspond to the choice of equilibrium distribution functions F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒ in Eq. ( 21).Substituting Eqs. ( 22) and ( 23) into Eqs.( 46)-( 49) readily gives in the beam interior where n2 b and n2 e are constants (independent of radius r).Here, f ϵ ne ͞Z b nb is the fractional charge neutralization and v2 pb 4p nb Z 2 b e 2 ͞g b m b is the ion plasma frequency squared.Because nb and ne are constants, Eqs. ( 44) and ( 45) are readily integrated with respect to t 0 .Enforcing the boundary conditions x 0 ͑t 0 t͒ x, ͑dx 0 ͞dt 0 ͒ t 0 t p x ͞g b m b , etc. readily gives for the ions, and for the electrons.Here, t t 0 2 t, nb and ne are defined in Eq. ( 50), and Eqs. ( 51) and ( 52) are valid in the beam interior where 0 # r 0 ͑x 0 2 1 y 0 2 ͒ 1͞2 , r b .Equations ( 51) and ( 52) constitute a Cartesian representation of the orbits in the equilibrium field configuration.The transverse orbits can also be expressed in a cylindrical coordinate representation ͑r 0 , u 0 ͒, where x 0 r 0 cosu 0 and y 0 r 0 sinu 0 .Introducing p x p Ќ cosw and p y p Ќ sinw, where w is the azimuthal momentum phase, it is readily shown from Eq. ( 51) that r 0 2 ͑t 0 ͒ x 0 2 ͑t 0 ͒ 1 y 0 2 ͑t 0 ͒ can be expressed as where ͑x, y͒ ͑r cosu, r sinu͒.Moreover, from Eq. ( 51), the azimuthal orbit determined from tanu 0 y 0 ͞x 0 is given by Note from Eq. ( 53) that the motion of r 0 2 ͑t 0 ͒ corresponds to that of a displaced simple harmonic oscillator, oscillating at constant frequency 2 nb .Also note from Eqs. ( 53) and (54) that ͑r 0 , u 0 ͒ ͑r, u͒ at t t 0 2 t 0, as required.Finally, from Eq. ( 52), the transverse orbits r 0 ͑t 0 ͒ and u 0 ͑t 0 ͒ for the electrons are identical in form to Eqs. ( 53) and (54) provided we make the replacements nb !ne and g b m b !m e .The simple form of the transverse orbits in Eqs. ( 51)-(54) will permit an exact evaluation in Sec.IV of the orbit integrals in Eqs. ( 40) and ( 41) for the choice of equilibrium distribution functions in Eqs. ( 21) and the corresponding step-function density profiles in Eqs. ( 22) and (23).

C. Particle orbits for continuously varying equilibrium profiles
We now examine the ion and electron orbit equations ( 44) and ( 45) for the case of continuously varying equilibrium density profiles n 0 b ͑r͒ and n 0 e ͑r͒.In this case, from Eqs. ( 48) and (49), the (depressed) betatron frequencies n 2 b ͑r͒ and n 2 e ͑r͒ defined in Eqs. ( 46) and (47) generally vary with radial coordinate r.Because the equilibrium is axisymmetric ͑≠͞≠u 0͒, an exact consequence of Eqs. ( 44) and ( 45) is the conservation of angular momentum; i.e., for both the electrons and the ions.
Here, P u xp y 2 yp x rp Ќ sin͑w 2 u͒ is the angular momentum expressed in phase-space variables, and use has been made of ͑x, y͒ ͑r cosu, r sinu͒ and ͑p x , p y ͒ ͑p Ќ cosw, p Ќ sinw͒.Another exact consequence of Eqs. ( 44) and ( 45) is the conservation of particle energy, i.e., H 0 Ќb ͑t 0 ͒ H Ќb const for the ions, and H 0 Ќe ͑t 0 ͒ H Ќe const for the electrons, where H Ќb ͑r, p Ќ ͒ and H Ќe ͑r, p Ќ ͒ are defined in Eq. (18).For example, conservation of energy for the electrons in cylindrical coordinates can be expressed as where P u const is the angular momentum.In obtaining Eq. ( 56), use has been made of m e r 0 2 ͑t 0 ͒du 0 ͑t 0 ͒͞dt 0 P u to express ͑1͞2͒m e r 0 2 ͑du 0 ͞dt 0 ͒ 2 P 2 u ͞2m e r 0 2 , and the electrostatic potential f 0 ͑r 0 ͒ is determined selfconsistently from Eq. (19).
For the ions, the equation of motion for r 0 ͑t 0 ͒ is identical in form to Eq. ( 56) provided we make the replacements m e !g b m b , v 0 2 be !v 0 2 bb , 2e !Z b e, f 0 ͑r 0 ͒ !c 0 ͑r 0 ͒, and H Ќe !H Ќb in Eq. (56).For present purposes, we therefore focus on an examination of Eq. (56).
It is convenient to rewrite Eq. ( 56) as an equation for r 0 2 ͑t 0 ͒; i.e., For specified potential profile f 0 ͑r͒, Eq. ( 57) is a convenient form for direct integration and determination of r 0 2 ͑t 0 ͒.Alternatively, taking the derivative of Eq. ( 57) with respect to t 0 gives Equations ( 57) and (58) are very useful representations of the radial orbit equation for r 0 2 ͑t 0 ͒.For example, for the step-function density profiles in Eqs. ( 22) and ( 23), it follows from Eq. ( 24) that e͓f 0 ͑r 0 ͒ 2 f0 ͔ 2͑1 2 f͒p nb Z b e 2 r 0 2 for 0 # r 0 , r b , and the coefficient of r 0 2 in Eq. (58 , where ne is defined in Eq. (50).Therefore, as expected from Eq. (58), the motion of r 0 2 ͑t 0 ͒ corresponds to a displaced oscillator, oscillating at frequency 2 ne const.Finally, for general equilibrium profiles, once the radial orbit r 0 2 ͑t 0 ͒ has been determined from Eq. (57) or Eq. ( 58), the azimuthal orbit u 0 ͑t 0 ͒ can be determined from m e r 0 2 ͑t 0 ͒du 0 ͑t 0 ͒͞dt 0 P u const, where P u is the angular momentum.Integrating with respect to t 0 and enforcing u͑t 0 t͒ u gives In circumstances where the equilibrium density profiles n 0 b ͑r͒ and n 0 e ͑r͒ vary continuously with r, the coefficient ͕• • •͖ of r 0 2 in Eq. (58) will depend on r 0 2 , corresponding to a spread in the (depressed) betatron-frequency due to self-field effects.To illustrate this effect, we consider a simple example where n parabolic form (Fig. 4) Here, e is a positive constant in the range 0 # e # 1.For e 0, Eq. ( 60) has the step-function form in Eq. ( 22), whereas for e 1, n 0 b ͑r͒ decreases monotonically to zero at the beam radius r r b .For simplicity, we further assume that the electron density profile n 0 e ͑r͒ has an identical shape to Eq. (60) with n 0 e ͑r͒ fZ b n 0 b ͑r͒, where f const is the fractional charge neutralization.Substituting into the equilibrium Poisson equation ( 19) for f 0 ͑r͒ then gives Here, v2 pb 4p nb Z 2 b e 2 ͞g b m b is the on-axis ͑r 0͒ ion plasma frequency-squared.We substitute Eq. (61) into Eq.(57) or Eq. ( 58) to determine the r 0 2 motion.For example, Eq. (58) becomes where n2 e is defined in Eq. ( 50) and the dimensionless coefficient h e is defined by For v 0 be 0, we note that Eq. (63) reduces to h e ͑3͞4͒e.Note also that the nonuniform density variation ͑h e fi 0͒ leads to a nonlinear frequency shift in Eq. (62).Substituting Eq. (61) into Eq.(57), or, alternatively, integrating Eq. ( 62) once with respect to t 0 gives where P u const is the angular momentum.The radial equation of motion for the ions in the equilibrium potential c 0 ͑r͒ calculated from Eqs. ( 19), (59), and ( 60) is identical in form to Eq. ( 64) provided we make the replacements ne !nb , H Ќe !H Ќb , m e !g b m b , and h e !h b in Eq. ( 64), where h b is defined by For step-function density profiles with e 0, and therefore h e 0 h b , the electron and ion orbits for r 0 2 ͑t 0 ͒ are identical to those calculated in Sec.III B, with oscillatory components at the constant frequencies 2 ne and 2 nb .For e fi 0 and h j fi 0, however, the spatial nonuniformity in the equilibrium density profiles produces a spread in (depressed) betatron frequencies for the particle orbits.For example, for specified values of H Ќe and P u , electrons that are confined within the beam ͑r 0 2 , r 2 b ͒ still exhibit periodic motion as a function of t t 0 2 t, but the period t e ͑H Ќe , P u ͒ for the r 0 2 ͑t 0 ͒ motion depends on the energy H Ќe and the angular momentum P u and is no longer equal to the constant value t 0 e 2p͞2 ne obtained for e 0 (see the Appendix).
It is also important to recognize that spatial nonuniformity in the equilibrium density profiles n 0 b ͑r͒ and n 0 e ͑r͒ can have a larger influence on the electron motion than the ion motion in many applications of practical interest.To illustrate this point, we consider the interesting case where v 0 be 0 and v 0 bb fi 0, and Eq. ( 63) reduces to Even for the high-intensity proton beams envisioned for the next-generation linacs and storage rings for spallation neutron sources and tritium production, the beam intensity is such that v2 pb ͞v 0 2 bb & 0.1.Therefore, from Eqs. (65) and (66), h b ø e for the protons, whereas h e 3e͞4 for the electrons, and the nonlinear effects and frequency spread are correspondingly larger for the electron orbits.
We now return to the electron orbit equation for r 0 2 ͑t 0 ͒ in Eq. (64), keeping in mind that the ion equation of motion is similar in form.As shown in the Appendix, for arbitrary inhomogeneity strength 0 # e # 1, Eq. ( 64) can be solved exactly in terms of elliptic integrals of the first kind.For present purposes, and future reference in Sec.V, we note here that the approximate orbit for r 0  correct to order h e .For e 0 h e , we note that Eq. (A33) reduces exactly to the electron analog of Eq. ( 53), as expected.On the other hand, for small h e fi 0, the nonuniformity in the equilibrium density profiles introduces a frequency spread in the r 0 2 ͑t 0 ͒ motion that depends on the energy H Ќe according to Eq. ( 67).
In Sec.V, we will make use of a simple model to obtain semiquantitative estimates of the influence of a nonuniformity-induced frequency spread on stability behavior.

IV. KINETIC STABILITY PROPERTIES FOR STEP-FUNCTION DENSITY PROFILES
We now return to the linearized Vlasov-Maxwell equations ( 40)-( 43), specializing to the case where the equilibrium ion and electron density profiles have the simple step-function forms in Eqs. ( 22), ( 23), and Fig. 1, and the corresponding monoenergetic equilibrium dis-tributions for F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒ are specified by Eq. (21).In this case, the particle motion in the beam interior ͑0 # r 0 # r b ͒ is that described in Sec.III B. In this section, the necessary orbit integrals for the ions and electrons are evaluated in closed analytical form, leading to a kinetic dispersion relation which is valid for arbitrary beam intensity v2 pb , fractional charge neutralization f, and azimuthal mode number ᐉ (Sec.IV A).The resulting dispersion relation is then analyzed in detail for the case of azimuthal mode number ᐉ 1, corresponding to a simple transverse displacement of the beam ions and background electrons (Sec.IV B).Finally, a brief discussion of stability behavior for quadrupole perturbations ͑ᐉ 2͒ is presented (Sec.IV C).

A. Kinetic dispersion relation
In Eqs. ( 42) and (43) for the potential amplitudes d ĉᐉ ͑r͒ and d fᐉ ͑r͒, the integrations over transverse momentum are expressed as , and w is the azimuthal phase of p Ќ in the transverse plane.Because the single-particle Hamiltonians, H Ќb ͑r, p Ќ ͒ and H Ќe ͑r, p Ќ ͒, are independent of phase w [see Eq. ( 18)] it follows that ≠F 0 b ͑H Ќb ͒͞≠H Ќb and ≠F 0 e ͑H Ќe ͒͞≠H Ќe are also independent of w.Therefore, when calculating the perturbed ion and electron charge densities from Eqs. ( 40) and (41), what is required are the phase-averaged orbit integrals, I ᐉ b ͑r, p Ќ ͒ and I ᐉ e ͑r, p Ќ ͒, defined by Here, t t 0 2 t, and u 0 ͑t 0 ͒ and r 0 ͑t 0 ͒ are the transverse particle orbits in the equilibrium fields that pass through u and r at time t 0 t.For the step-function density profiles consistent with Eqs. ( 21)-( 23), we will subsequently find that Eqs. ( 40 Here, ĉᐉ and fᐉ are constant amplitudes.Because the ion and electron orbits for x 0 ͑t 0 ͒ and y 0 ͑t 0 ͒ have similar functional forms [compare Eqs. ( 51) and ( 52)], we focus here on an evaluation of the ion orbit integral I ᐉ b ͑r, p Ќ ͒ defined in Eq. (68).
It is now straightforward to evaluate the perturbed ion and electron densities, 74).For the choice of equilibrium distribution functions in Eq. ( 21), it can be shown that [30,47] 2p Here, n2 b and n2 e are defined in Eq. ( 50), ne fZ b nb is the electron density, where f is the fractional charge neutralization, and the beam radius r b is related self-consistently to other equilibrium parameters by the force-balance constraints in Eq. (26).
We make use of Eqs. ( 74 We note that the perturbed charge and current densities on the right-hand side of Eqs. ( 78) and (79) correspond to surface-charge perturbations localized to the beam surface at r r b (see Fig. 1).Therefore, the exact solutions to Eqs. ( 78) and (79) can be expressed as and In Eqs. ( 80) and ( 81 and " where e !0 1 .Note that Eqs. ( 82) and ( 83) relate the discontinuities in the perturbed radial electric field, 2͑≠͞≠r͒d fᐉ ͑r͒, and the perturbed azimuthal magnetic field, 2͑≠͞≠r͒d Âᐉ z ͑r͒ [where to the perturbed surface charge and current densities at r r b .Substituting Eqs. ( 80) and (81) into Eqs.( 82) and (83), we obtain two coupled equations relating the potential amplitudes ĉᐉ and fᐉ ; i.e., " The condition for a nontrivial solution to Eq. ( 84) with nonzero ĉᐉ and fᐉ is that the two-by-two determinant of the coefficients of ĉᐉ and fᐉ vanish.This gives Here, the ion and electron susceptibilities, e ͑v͒, are defined in Eq. ( 75) and the (depressed) betatron frequencies, nb and ne , are defined in Eq. (50).Equation ( 85) is the final form of the fully kinetic dispersion relation, derived from the linearized Vlasov-Maxwell equations for small-amplitude perturbations about the equilibrium distribution functions, F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒, in Eq. ( 21), and the corresponding stepfunction density profiles in Eqs. ( 22) and ( 23).As such, Eq. ( 85) can be used to determine the complex oscillation frequency v over a wide range of system parameters, including beam intensity ͑ v2 pb ͒, fractional charge neu-tralization ͑f ne ͞Z b nb ͒, focusing field strength ͑v 0 2 bb ͒, azimuthal mode number ͑ᐉ͒, axial wave number ͑k z ͒, etc., subject only to the simplifying assumptions summarized in Sec.II.
In the absence of electrons ͑ ne 0 v2 pe ͒, the dispersion relation (85) reduces to the simple result When background electrons are present ͑ ne fi 0͒, however, Eq. ( 85) supports unstable solutions ͑Imv .0͒ with instability resulting from the axial streaming ͑V b fi 0͒ of the beam ions through the background electrons.

B. Electron-ion instability for azimuthal mode number ᐉ ᐉ ᐉ 5 1
We defer an analysis of the dispersion relation (85) for general azimuthal mode number ᐉ to a subsequent investigation and focus the present analysis on dipole perturbations with azimuthal mode number ᐉ 1, corresponding to a simple transverse displacement of the beam ions and background electrons.A brief discussion of the dispersion relation for quadrupole perturbations with mode number ᐉ 2 is given in Sec.IV C.
For azimuthal mode number ᐉ 1, it follows from the definitions of the electron and ion susceptibilities in Eq. ( 75) that Here, v2 pb 4p nb Z 2 b e 2 ͞g b m b and v2 pe 4p ne e 2 ͞m e are the ion and electron plasma frequency-squared, respectively, and the (depressed) betatron frequencies, nb and ne , are defined in Eq. ( 50).Substituting Eq. ( 87) into the dispersion relation (85) gives for ᐉ 1.We define the electron and ion collective oscillation frequencies, v e and v b , by and where v2 pe has been expressed as v2 pe ͑g b m b ͞ Z b m e ͒f v2 pb .Substituting Eqs. ( 89) and (90) into Eq.( 88) and rearranging terms, the ᐉ 1 dispersion relation (88) can be expressed in the equivalent (compact) form where v f is defined by In the absence of background electrons (f 0 and v f 0), Eq. ( 91) gives stable collective oscillations of the ion beam with frequency v 2 k z V b 6v b , where v b is defined in Eq. ( 90).For f fi 0, however, the ion and electron terms on the left-hand side of Eq. ( 91) are coupled by the v 4 f term on the right-hand side, leading to one unstable solution with Imv .0. The instability is two-stream in nature and results from the directed ion motion with axial velocity V b through the (stationary) background electrons.Equation ( 91) is a fourth-order algebraic equation for the complex oscillation frequency v.Some straightforward analysis shows that there are two stable solutions to Eq. ( 91) with purely real v and two complex solutions that are complex conjugates (one is growing with Imv .0 and the other is damped with Imv , 0) for certain ranges of axial wave number k z .Equation ( 91) can of course be solved numerically for v over a wide range of system parameters, including fractional charge ͑f͒, normalized beam intensity ͑ v2 pb ͞v 0 2 bb ͒, proximity of the conducting wall ͑r b ͞r w ͒, etc.
In the interesting application to proton linacs and storage rings, the natural oscillation frequencies, v e and v b , occurring in Eq. (91) tend to be large in comparison with v f , even for fractional charge neutralization f 1 and moderately large normalized beam intensity v2 pb ͞v 0 2 bb & 0.5.In this case, it is found that the unstable branch in Eq. ( 91) has real frequency and axial wave number ͑v, k z ͒ very closely tuned to the values ͑v 0 , k z0 ͒ defined by Here, we consider the electron branch with positive frequency, v Ӎ 1v e , and the down-shifted ion branch (in the beam frame) with v 2 k z V b Ӎ 2v b , which are (unstably) coupled by the v 4 f term in Eq. ( 91).Expressing v v 0 1 dv and k z k z0 1 dk z , and assuming jdvj ø 2v e and jdv 2 dk z V b j ø 2v b , it is straightforward to show that Eq. ( 91) can be approximated by the quadratic dispersion relation  94) supports an unstable solution with Imdv .0 for dk z in the interval 22G 0 , dk z V b , 12G 0 .In this range of dk z , the unstable solution to Eq. ( 94) is given by Redv As evident from Eq. ( 95) and illustrated in Fig. 5, the growth rate Imdv is a symmetric function of dk z , achieving a maximum value of G 0 for dk z 0 and decreasing to zero for dk z V b 62G 0 .From Eqs. ( 89), (90), and (92), the maximum growth rate, ͑Imdv͒ max G 0 , is given explicitly by Equation ( 96) can be applied over a wide range of system parameters subject to the assumptions jdvj ø 2v b , 2v e , or, equivalently, Several points are noteworthy from Eq. ( 96).First, the growth rate increases with increasing fractional charge neutralization ͑f͒ and increasing normalized ion beam intensity ͑ v2 pb ͞v 0 2 bb ͒.Second, the growth rate decreases when the conducting wall is in close proximity to the beam (larger values of r b ͞r w ).Third, there is no threshold value of f for the onset of instability, which is likely a consequence of the fact that there is no spread in (depressed) betatron frequencies for the step-function density profiles in Eqs. ( 22) and ( 23) and the corresponding monoenergetic ion and electron distributions in Eq. ( 21) considered in this section.Finally, the linearized Vlasov equations ( 30) and (31), and therefore the resulting dispersion relation (88), neglect the effects of Landau damping [44] associated with an axial velocity spread in the beam ions or background electron population.The expression for ͑Imdv͒ max in Eq. (96) simplifies further for the case of negligible applied focusing force on the electrons ͑v 0 be 0͒ and large conducting wall radius ͑r w ¿ r b ͒.
Setting v 0 be 0 and r b ͞r w !0 in Eq. ( 96) gives For example, for a proton beam ͑Z b 1, m b ͞m e 1836͒ with relativistic mass factor g b 1.85, a moderate value of normalized beam intensity v2 pb ͞v 0 2 bb 0.1, and fractional charge neutralization f 0.1, Eq. ( 98) gives ͑Imdv͒ max 0.127v 0 bb , corresponding to a particularly virulent growth rate for the e-p instability.For this choice of system parameters, the central oscillation frequency and wave number calculated from Eq. ( 93) are v 0 13.03v 0 bb and k z0 V b 14.03v 0 bb .For completeness, shown in Fig. 6 is a plot of normalized growth rate ͑Imdv͒ max ͞v 0 bb versus normalized beam intensity v2 pb ͞v 0 2 bb calculated from Eq. (96) for the choice of system parameters v 0 be 0, g b 1.85, m b ͞m e 1836, r b ͞r w 0.5, and several values of fractional charge neutralization corresponding to f 0.1, 0.3, 0.6, and 0.9.Note that ͑Imdv͒ max increases with increasing beam intensity and increasing charge neutralization, as expected.
While the beam distribution function in the PSR experiment [31][32][33][34][35] has a (small) energy spread, and is not modeled well by the monoenergetic distribution in Eq. ( 21), it is nonetheless instructive to apply the stability results obtained in Eqs. ( 94) and (96) to characteristic parameters in PSR.An instability, believed to be caused by trapped electrons in the proton beam, has been observed in PSR [31][32][33][34][35] for coasting beams at currents exceeding I b 2. pb ͞v 02 bb 0.025.Substituting into Eq.( 96) gives a predicted maximum growth rate ͑Imdv͒ max 6.9 3 10 5 s 21 for step-function density profiles.On the other hand, the observed growth rate in PSR at a beam current of I b 2.23 A is in the range of 3 3 10 4 s 21 .The discrepancy between theory and experiment is likely due to the reduction in growth rate caused by a spread in depressed betatron frequencies [see, e.g., Ref. [28] and Sec.V], an effect which is not incorporated in the monoenergetic distributions in Eq. ( 21).
The quadratic approximation to the dispersion relation given in Eq. ( 94) is valid for the moderately high beam intensities ͑ v2 pb ͞v 0 2 bb ͒ envisioned in the proton linacs and storage rings for tritium production and spallation neutron sources [51].For heavy ion fusion applications [6,7], however, the beam emittance, which is proportional to TЌb in Eq. ( 26), is very low and the normalized beam intensity is such that v2 pb ͞2g 2 b v 0 2 bb can approach unity in the absence of background electrons ͑f 0͒.This is evident from the equilibrium force constraint for the ions in Eq. (26), where, for f 0 and 2 TЌb ͞g b m b r 2 b v 0 2 bb ø 1, it follows that v2 pb ͞2g 2 b v 0 2 bb ! 1.At such high beam intensities, it follows that it is necessary to solve the full quartic dispersion relation (91) for the complex oscillation frequency v. Typical numerical results obtained from Eq. (91) are illustrated in Fig. 7. Here, ͑Imv͒͞v 0 bb and ͑Rev 2 v e ͒͞v 0 bb are plotted versus ͑k z 2 k z0 ͒V b ͞v 0 bb for several values of v2 pb ͞v 0 2 bb ranging from 0.1 to 2.0.Other system parameters in Fig. 7 correspond to Z b 1, mass number A m b ͞m p 200, ͑g b 2 1͒m b c 2 10 GeV, r b ͞r w 0.5, f 0.1, and v 0 be 0.For sufficiently small values of v2 pb ͞v 0 2 bb , the numerical results obtained in Fig. 7 from the full quartic dispersion relation (91) are in excellent agreement with the approximate quadratic dispersion relation in Eq. (94).On the other hand, at very high beam intensity with v2 pb ͞v 0 2 bb 2, say, it is evident from Fig. 7 that the growth rate Imv͞v 0 bb has very large bandwidth and becomes significantly skewed about k z k z0 [in contrast with the symmetric results obtained from the quadratic approximation in Eq. ( 95)].It is also striking from Fig. 7 that the instability growth rate can be large for the very high beam intensities of interest for heavy ion fusion.For example, ͑Imv͒ max 2.15v 0 bb for v2 pb ͞v 0 2 bb 2 in Fig. 7.In concluding Sec.IV B, it is important to recognize that the dispersion relation (85) has been derived under the assumptions [Eq.( 8)] of long-wavelength, high-frequency perturbations satisfying k 2 z r 2 b ø 1, jv͞k z 2 V b j ¿ y T bz , and jv͞k z j ¿ y T ez , which allowed us to neglect kinetic effects (such as Landau damping) in the z direction [44].
Here, y T bz ͑2T bz ͞g b m b ͒ 1͞2 and y T ez ͑2T ez ͞m e ͒ 1͞2 are the characteristic axial thermal speeds of the beam ions and background electrons, respectively.For mode number ᐉ 1, we estimate v Ӎ v 0 v e and k z Ӎ k z0 V 21  b ͑v e 1 v b ͒ [see Eq. ( 93)] and make use of Eqs. ( 26) and ( 50 where nb , v b , and v e are defined in Eqs. ( 50), (89), and (90).Equation (99) clearly requires that the directed axial velocity V b be large in comparison with the thermal speeds y T bЌ , y T bz , and y T ez .Moreover, the first two inequalities in Eq. ( 99) are the most difficult to satisfy because v e is typically larger than v b and nb .

C. Dispersion relation for azimuthal mode number ᐉ ᐉ ᐉ 5 2
Detailed analysis of the kinetic dispersion relation (85) for perturbations with azimuthal mode numbers ᐉ $ 2 will be the subject of a future investigation.For present purposes, we briefly summarize here properties of the dispersion for quadrupole perturbations with ᐉ 2. For ᐉ 2, Eq. ( 75) gives for the ion and electron susceptibilities Substituting Eq. (100) into Eq.( 85), we obtain for ᐉ 2 Equation ( 101) is similar in form to Eq. (88) (for ᐉ 1), but the particle resonances in the denominators in Eq. ( 101 Comparing Eq. (102) with Eq. (91), it is clear that the ᐉ 2 dispersion relation (102) is identical in form to the ᐉ 1 dispersion relation (91), with simple redefinitions of the frequencies v b , v e , and v f .Therefore, for f fi 0, Eq. (102) always has one unstable solution with Imv .0 for certain ranges of k z , and the analysis developed in Sec.IV B for ᐉ 1 can also be applied to Eq. ( 102) for azimuthal mode number ᐉ 2. Without presenting algebraic details, for specified values of f, v2 pb ͞v 0 2 bb , g b , etc., the ᐉ 2 growth rate calculated from Eq. ( 102) is smaller than the ᐉ 1 growth rate calculated from Eq. (91).For example, the quadratic dispersion relation estimate, ͑Imdv͒ max v 2 f ͞2͑v e v b ͒ 1͞2 , gives a smaller growth rate for ᐉ 2, using the new definitions for v b , v e , and v f .In this regard, the ᐉ 1 mode is the most "dangerous" mode because it has the largest growth rate.

V. EFFECTS OF A SPREAD IN BETATRON FREQUENCIES
The general kinetic eigenvalue equations ( 40)-(43) developed in Sec.III can be applied to electrostatic pertur-bations about a wide range of nonmonoenergetic equilibrium distribution functions, F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒, and corresponding self-consistent equilibrium density profiles, n 0 b ͑r͒ and n 0 e ͑r͒, that vary continuously with radial coordinate r.A detailed, self-consistent stability analysis based on Eqs. ( 40)-(43) for continuously varying equilibrium profiles is beyond the scope of the present article and will be the subject of a future investigation.For present purposes, based on the insights gained in Secs.III C and IV and the Appendix, we summarize the results of a simple model that illustrates semiquantitatively the stabilizing influence [28] that a (weak) density nonuniformity and the corresponding spread in betatron frequencies can have on stability behavior.
The model assumes overlapping ion and electron density profiles with parabolic profile shape specified by Eq. (60) for the ion density profile n 0 b ͑r͒, and electron density profile specified by n 0 e ͑r͒ Z b fn 0 b ͑r͒, where f const is the fractional charge neutralization.The model further makes the ansatz that the dispersion relation (considered here for azimuthal mode number ᐉ 1) has similar form to Eq. ( 85), obtained by making the susceptibility replacements [see Eq. ( 87 In Eq. ( 104), h b and h e are both proportional to the inhomogeneity parameter e and are defined in Eqs. ( 65) and (63), respectively.Moreover, nb and ne are the (constant) orbital oscillation frequencies defined in Eq. ( 50) for e 0. It follows trivially for e 0 that Eq. ( 103) reduces exactly to the susceptibility expressions in Eq. ( 87), obtained for step-function density profiles.For e fi 0, however, Eq. ( 103) incorporates the effects of spatial nonuniformity in the equilibrium density profiles and the corresponding spreads in betatron oscillation frequencies, at least in the context of the present simple model.Some straightforward algebraic manipulation shows that the r integrations in Eq. ( 103) can be carried out exactly for e in the interval 0 # e # 1 to give , For present purposes, we limit the analysis of Eq. ( 105) and the corresponding dispersion relation to the case of weak spatial nonuniformity with e, h b , and h e ø 1.
Taylor expanding Eq. (105) gives to first order in e x e ͑v͒ 2 v2 where a b and a e are defined by Here, h b ͞e and h e ͞e are defined in terms of v2 pb , v 0 2 bb , etc., in Eqs. ( 65) and (63), respectively.As expected, Eq. (106) reduces exactly to the susceptibility expressions in Eq. (87) in the limit of step-function density profiles with e 0.
We substitute Eqs. ( 103) and (106) into Eq.(85) for azimuthal mode number ᐉ 1. Rearranging terms, this gives the dispersion relation where the coefficients a 0 b and a 0 e are defined by a 0 b ϵ a b ͑v 0 2 k z0 V b ͒ and a 0 e ϵ a e ͑v 0 ͒, i.e., The quadratic dispersion relation (110) of course reduces exactly to Eq. ( 94) for the case of uniform step-function density profiles with e 0. In this case, as discussed in Sec.IV B, there is no threshold in v2 pb and f for the onset of instability.The situation is different when e fi 0 and the (small) frequency spreads Dv 0 2 b and Dv 0 2 e defined in Eq. ( 111) are nonzero.A straightforward analysis of Eq. ( 110) for dk z 0 shows that the necessary and sufficient condition for instability (existence of a solution with Imdv .0) is given by 1 4 Here, we have approximated ͑1 2 ea 0 e ͒ ͑1 2 ea 0 b ͒ Ӎ 1 to good accuracy on the left-hand side of Eq. ( 113).Whenever the frequency spread is sufficiently large that the inequality in Eq. ( 113) is violated, the system is stable with Imdv 0. Making use of the definitions of v 4 f , Dv 0 2 b , and Dv 0 2 e in Eqs. ( 92) and ( 111), the necessary and sufficient condition for instability in Eq. ( 113) can be expressed in the equivalent form where v e , v b , a 0 e , and a 0 b are defined in Eqs. ( 89), (90), and (112).
Equation ( 114) leads to a threshold in beam intensity for the onset of the two-stream instability.To illustrate this point, we consider the case where v 0 be 0 and r w ͞r b !`.Equations (89), (90), and (112) then give Making use of Eq. (115) in Eq. ( 114), the necessary and sufficient condition for the onset of instability can be expressed as The ion term proportional to a 0 b ͑m e ͞m b ͒ 1͞2 on the righthand side of Eq. ( 116) is small in comparison with the electron term ͑fa 0 e ͒ in the parameter regimes of practical interest.Therefore, making use of fa 0 e 3͞16 1 ͑5͞16͒f, Eq. (116) reduces to The derivation of the quadratic dispersion relation (110) and therefore the instability threshold condition in Eq. ( 117) have assumed weak spatial nonuniformity with e ø 1 and moderate values of beam intensity with v2 pb ͞v 0 2 bb & 0.5, say (see discussion in Sec.IV B).
As a simple example, we consider Eq. (117) for inhomogeneity parameter e 0.075 and protons with g b 1.85, Z b 1, and m b ͞m e 1836.Equation (117) then reduces to If, for example, the fractional charge neutralization by the background electrons is f 0.
To briefly summarize, the present simple model provides semiquantitative evidence that a spread in betatron frequencies (particularly for the electrons) induced by spatial nonuniformity in the density profiles has a stabilizing influence on the two-stream instability, leading to a threshold for the onset of instability.It should be reiterated, however, that a fully self-consistent treatment of electrostatic stability properties for nonmonoenergetic distributions with continuously varying equilibrium profiles should be based on the kinetic eigenvalue equations ( 40) -(43) derived in Sec.III.

VI. CONCLUSIONS
The present analysis made use of the Vlasov-Maxwell equations to develop a fully kinetic description of the electrostatic, electron-ion two-stream instability driven by the directed axial motion of a high-intensity ion beam propagating through a stationary population of background electrons.The basic assumptions, theoretical model, and examples of self-consistent equilibrium solutions for a continuous ion beam-electron plasma system were discussed in Sec.II.In Sec.III, the linearized Vlasov-Maxwell equations were formally integrated using the method of characteristics and the properties of the ion and electron orbits in the applied field plus equilibrium self-field configuration were discussed, both for the case of overlapping step-function ion and electron density profiles, corresponding to monoenergetic beam ions and monoenergetic electrons, and for the case where the equilibrium density profiles have a continuous variation with radius r, corresponding to a spread in (depressed) betatron frequencies.The necessary orbit integrals were evaluated in closed analytical form in Sec.IV for the case of stepfunction ion and electron density profiles, leading to a kinetic dispersion relation which is valid for arbitrary normalized beam intensity v2 pb ͞v 0 2 bb , fractional charge neutralization f, and azimuthal mode number ᐉ.The resulting dispersion relation was analyzed in detail for the case of azimuthal mode number ᐉ 1, which corresponds to the strongest instability (largest growth rate).As a general remark, the instability growth rate is found to increase with increasing beam intensity v2 pb ͞v 0 2 bb and increasing fractional charge neutralization f, and decrease with increasing proximity of the conducting wall r b ͞r w .For monoenergetic ions and electrons and the corresponding step-function density profiles considered in Sec.IV, a key result is that there is no threshold in beam intensity or fractional charge neutralization for the onset of instability.Finally, for the case of continuously varying density profiles with parabolic profile shape, in Sec.V we made use of a simple model to obtain a semiquantitative estimate of the effects of the corresponding spread in (depressed) betatron frequency on stability behavior, including an estimate of the instability threshold for the case of weak density nonuniformity.As expected, it is the spread in the electron oscillation frequency that has the largest stabilizing influence.
In conclusion, we reiterate that the kinetic eigenvalue equations ( 40) -( 43) can be applied to electrostatic perturbations about a wide range of nonmonoenergetic equilibrium distribution functions, F 0 b ͑H Ќb ͒ and F 0 e ͑H Ќe ͒, and corresponding equilibrium density profiles, n 0 b ͑r͒ and n 0 e ͑r͒, that vary continuously with radial coordinate r.In future investigations, we will make use of the kinetic eigenvalue equations ( 40)-(43) to determine self-consistently the influence of spreads in electron energy H Ќe and ion energy H Ќb on stability behavior for continuously varying equilibrium profiles.

ACKNOWLEDGMENTS
This research was supported by the U.S. Department of Energy, and by the APT Project, the Short Pulse Spallation Source Project, and the LANSCE Division of Los Alamos National Laboratory.It is a pleasure to acknowledge the benefit of useful discussions with Thomas P. Wangler and Han S. Uhm.

APPENDIX: PARTICLE MOTION IN PARABOLIC DENSITY PROFILES
To analyze the radial orbit equation (64), we introduce the dimensionless quantities R 0 2 ͑t 0 ͒, t 0 , h, and p 2 defined by where the particle motion in the equilibrium field configuration is restricted to the region 0 # R 0 2 , 1 (the beam interior).Substituting Eq. (A1) into Eq.(64) gives in dimensionless variables Uniform beam density ͑e 0 h e ͒: For h e 0, Eq. (A2) reduces to where Motion is allowed by Eq. (A3) for R 0 2 ͑t 0 ͒ in the interval where R 0 2 ͑t 0 t͒ R 2 ϵ r 2 ͞r 2 b and use has been made of t 0 2 ne ͑t 0 2 t͒.Note that for e 0 the oscillations in R 0 2 ͑t 0 ͒ are exactly at frequency 2 ne const, as expected [compare with Eq. ( 53)].
Cartesian orbits x 0 ͑t 0 ͒ and y 0 ͑t 0 ͒ for weak inhomogeneity ͑e ø 1͒: In the previous section of this Appendix, we showed for weak spatial nonuniformity with e ø 1 and h e ø 1 that the radial orbit for r 0 2 ͑t 0 ͒ x 0 2 ͑t 0 ͒ 1 y 0 2 ͑t 0 ͒ oscillates at one distinct frequency 2n e ͑H Ќe ͒, where n e ͑H Ќe ͒ is defined in Eq. (A30) [see Eq. (A33)].This is true both for P u 0 and for P u fi 0. Whenever P u fi 0, however, analysis of the complete Cartesian orbit equation for x 0 Ќ ͑t 0 ͒ x 0 ͑t 0 ͒ êx 1 y 0 ͑t 0 ͒ê y in the beam A careful analysis of Eq. (A35) for h e ø 1 shows that the transverse orbits for x 0 ͑t 0 ͒ and y 0 ͑t 0 ͒ are given to leading order by x 0 ͑t 0 ͒ A 1 cos͓n 1 e ͑t 0 2 t͒ 2 c 1 ͔ 1 A 2 cos͓n 2 e ͑t 0 2 t͒ 2 c 2 ͔ , y 0 ͑t 0 ͒ 2A 1 sin͓n 1 e ͑t 0 2 t͒ 2 c 1 ͔ (A38) 1 A 2 sin͓n 2 e ͑t 0 2 t͒ 2 c 2 ͔ , Here, the constant amplitudes ͑A 2 , A 1 ͒ and phases ͑c 2 , c 1 ͒ are defined in terms of the phase-space variables ͑x, y, p x , p y ͒ by To the level of accuracy of Eq. (A38), a careful examination of the transverse orbits for x 0 ͑t 0 ͒ and y 0 ͑t 0 ͒ shows that the radial orbit for r 0 2 ͑t 0 ͒ x 0 2 ͑t 0 ͒ 1 y 0 2 ͑t 0 ͒ oscillates at a single distinct frequency given by which is identical to 2n e ͑H 0 Ќe ͒ defined in Eq. (A30).Whenever P u fi 0, however, the individual x 0 ͑t 0 ͒ and y 0 ͑t 0 ͒ motions have frequency components at the two distinct frequencies, n 1 e ͑H 0 Ќe , P u ͒ and n 2 e ͑H 0 Ќe , P u ͒, defined in Eq. (A36), which are separated by an amount Finally, it should be pointed out that the orbit analysis in this Appendix can also be applied to the ion motion in parabolic density profiles by making the obvious replacements, m e !g b m b , 2e !Z b e, h e !h b , etc.

Here, f ϵ
FIG. 1. Equilibrium step-function density profiles n 0 b ͑r͒ and n 0 e ͑r͒ [Eqs.(22) and (23)] for the choice of ion and electron distribution functions in Eq.(21).For the purpose of illustration, we have taken Z b 1 and f 0.5 in the figure.

FIG. 3 .
FIG. 3. Plots versus radius r of (a) the ion density profile n 0 b ͑r͒͞ nb and ( b) the electron density profile n 0 e ͑r͒͞ ne obtained numerically from Eqs. (19) and (29) for the choice of system parameters v 0 be 0, g b 1.85, v2 pb ͞v 0 2 bb 0.1, Z b 1, m b ͞m e 1836, and f ne ͞ nb 0.2.

FIG. 4 .
FIG. 4. Plots versus normalized radius r͞r b of the equilibrium ion density profile n 0 b ͑r͒ defined in Eq. (60) and the electron density profile n 0 e ͑r͒ Z b fn 0 b ͑r͒.For the purpose of illustration, we have taken e 0.5, f 0.5, and Z b 1 in the figure.
), ĉᐉ and fᐉ are constant amplitudes and we have enforced the boundary conditions, d ĉᐉ ͑r r w ͒ 0 d fᐉ ͑r r w ͒ at the perfectly conducting wall at r r w .We have also enforced continuity of d fᐉ ͑r͒ and d ĉᐉ ͑r͒ at the surface of the beam ͑r r b ͒.The remaining boundary conditions are obtained by integrating Eqs.(78) and (79) from r b ͑1 2 e͒ to r b ͑1 1 e͒ across the surface of the beam at r r b and taking the limit e !0 1 .Operating on Eqs.(78) and (79) with R r b ͑11e͒ r b ͑12e͒ dr r • • • readily gives "
23 A. If we take illustrative parameters for PSR to be I b p nb r 2 b eb b c 2.23 A, g b 1.85, b b 0.84, f 0.02, r b 2 cm, r w 5 cm, v 0 be 0, and v 0 bb 4.05 3 10 7 s 21 , then nb 4.4 3 10 7 cm 23 and v2

FIG. 7 .
FIG. 7. Plots of (a) normalized growth rate Imv͞v 0 bb and (b) normalized real frequency ͑Rev 2 v e ͒͞v 0 bb versus shifted axial wave number ͑k z 2 k z0 ͒V b ͞v 0 bb obtained numerically from the full dispersion relation (91) for the unstable branch with positive real frequency.System parameters correspond to Z b 1, m b ͞m p 200, ͑g b 2 1͒m b c 2 10 GeV, r b ͞r w 0.5, f 0.1, and v 0 be 0. Curves are shown for several values of normalized beam intensity v2 pb ͞v 0 2 bb ranging from 0.1 to 2.0.

) (for ᐉ 2 )
occur at the second harmonics of nb and ne , i.e., at v 2 k z V b 62 nb and v 62 ne .Expressing v2 pe ͑g b m b ͞Z b m e ͒f v2 pb , where f ne ͞Z b nb , the ᐉ 2 dispersion relation (101) can be rewritten in the form ( assumes, for weak spatial nonuniformity with e ø 1, that the eigenfunctions d ĉᐉ ͑r͒ and d fᐉ ͑r͒ have (approximately) the same radial dependence [Eqs.(80) and (81)] as obtained for the case of step-function density profiles in Sec.IV.In Eq. (103), v 2 pj ͑r͒ v2 pj ͑1 2 er 2 ͞r 2 b ͒ for j b, e, and we make use of the radial orbit equation (62) to express

4 f
ea b ͒ , (108) where v 2 e , v 2 b , and v 4 f are defined exactly as in Eqs.(89), (90), and (92) in Sec.IV, and the inhomogeneity-induced frequency spreads, Dv 2 b and Dv 2 e , are defined by pb 4p nb Z 2 b e 2 ͞g b m b , and use has been made of v2 pe ͑g b m b ͞Z b m e ͒f v2 pb .In Eq. (108), the frequency spreads Dv 2 b and Dv 2 e are treated as small for e ø 1 and, similar to Sec.IV B, the dispersion relation is analyzed for moderate values of normalized beam intensity v2 pb ͞v 0 2 bb .In this case, a quadratic approximation to Eq. (108) is valid in which the frequency and wave number ͑v, k z ͒ of the unstable branch are very closely tuned to the values v 0 v e and v 0 2 k z0 V b 2v b defined in Eq. (93).Expressing v v 0 1 dv and k z k z0 1 dk z , and paralleling the analysis in Sec.IV B, Eq. (108) reduces to the approximate form √ 4v e v b ͑1 2 ea 0 e ͒ ͑1 2 ea 0 b ͒ .
n e ͑H Ќe , P u ͒ ϵ ne r # r b and arbitrary e in the interval 0 # e # 1.

2 e
is linearly proportional to the angular momentum ͑P u ͒ and the strength of the density nonuniformity ͑h e ͒.