Anisotropy-driven collective instability in intense charged particle beams

1098-4402= The classical electrostatic Harris instability is generalized to the case of a one-component intense charged particle beam with anisotropic temperature including the important effects of finite transverse geometry and beam space charge. For a long, coasting beam, the eigenmode code bEASt have been used to determine detailed 3D stability properties over a wide range of temperature anisotropy and beam intensity. A simple theoretical model is developed which describes the essential features of the linear stage of the instability. Both the simulations and the analytical theory clearly show that moderately intense beams are linearly unstable to short-wavelength perturbations provided the ratio of the longitudinal temperature to the transverse temperature is smaller than some threshold value. The delta-f particle-in-cell code BEST has been used to study the detailed nonlinear evolution and saturation of the instability.


I. INTRODUCTION
Periodic focusing accelerators, transport systems and storage rings [1][2][3] have a wide range of applications ranging from basic scientific research in high energy and nuclear physics, to applications such as heavy ion fusion, spallation neutron sources, tritium production and nuclear waste transmutation, to mention a few examples.Of particular importance at the high beam currents and charge densities of practical interest, are the effects of the intense self fields produced by the beam space charge and current on determining the detailed equilibrium, stability and transport properties.Charged particle beams confined by external focusing fields represent an example of nonneutral plasma [4].A characteristic feature of such plasmas is the nonuniformity of the equilibrium density profiles and the nonlinearity of the self fields, which makes detailed analytical investigation very difficult.The development and application of numerical tools such as eigenmode codes [5,6] and Monte-Carlo particle simulation methods [5,[7][8][9][10][11] is often the only tractable approach to understand the underlying physics.It is therefore particularly important to develop simplified models of different instabilities familiar in electrically neutral plasmas which include the important effects of strong space-charge and finite geometry, and which may cause a degradation in beam quality.One such instability is the electrostatic Harris instability [12] driven by a large temperature anisotropy.Anisotropic temperatures T jjb =T ?b 1 in the beam frame, where subscript jj? denotes parallel (perpendicular) to the beam propagation direction, develop naturally in accelerators.Indeed, due to the conservation of energy for particles with charge e b and mass m b accelerated by a voltage V, the energy spread of particles in the beam does not change, and (nonrelativistically) E bi m b v 2  bi =2 E bf m b V b v bf , where V b 2e b V=m b 1=2 is the average beam velocity after acceleration.Therefore, the velocity spread-squared, or equivalently the temperature, changes according to T jjbf T 2 jjbi =2e b V (for a nonrelativistic beam).In addition, the transverse temperature may increase due to nonlinearities in the applied focusing field, intense self-field forces, nonstationary beam profiles, and beam mismatch.It is well known that in electrically neutral plasmas with uniform magnetic field and strongly anisotropic electron distribution T jje =T ?e 1, where subscript jj? denotes parallel (perpendicular) to the magnetic field, an electrostatic (Harris-like) collective instability [12] may develop if the plasma is sufficiently dense that !pe > !ce , where !pe 4e 2 n=m 1=2 is the electron plasma frequency, and !ce eB=mc is the electron cyclotron frequency.For the case of charged particle beams in accelerators, the cyclotron oscillations in the applied magnetic field are replaced by the betatron oscillations of the beam particles in the combined applied and self-generated fields [5][6][7][8][13][14][15][16][17][18][19][20].Heavy ion experiments for high energy density physics and fusion applications require transporting high-current beams when the average depressed betatron frequency of the beam particles is much smaller than the average plasma frequency of the beam particles.Figure 1 shows a comparison between the conditions for the Harris instability to occur in electrically neutral plasma confined by an axial magnetic field B 0 , and for a nonneutral ion beam confined transversely in a focusing channel.In both cases the electric field of the density perturbation couples to the transverse oscillatory particle motion.The resulting instability is somewhat similar to the two-stream instability, with the source of free energy to drive the instability associated with the transverse particle motion.The instability may lead to a deterioration of the beam quality and an increase in the longitudinal velocity spread, which would make focusing the beam difficult and impose a limit on the minimum spot size achievable in focusing the intense beam onto a target.
The organization of this paper is the following.In Sec.II, we discuss the physical mechanism for the instability and present a simplified kinetic model for the most unstable dipole mode.The linear numerical results obtained with the eigenmode code bEASt are also described in Sec.II and compared with the model.Many of the linear numerical results in Sec.II are based on an invited paper presented at the 2005 Particle Accelerator Conference [20].The nonlinear f simulation method [21] is briefly described in Sec.III, and detailed simulation results of nonlinear stage and saturation of the instability [20] obtained with the nonlinear f beam equilibrium, stability and transport (BEST) code [9] are presented.

II. SIMPLIFIED MODEL OF ELECTROSTATIC HARRIS-TYPE INSTABILITY
For simplicity, the subsequent analysis is carried out in the beam frame (V b 0).In what follows, it is convenient to introduce the effective depressed betatron frequency !?defined by [7] where T ?b is the transverse beam temperature, r b is the root-mean-square beam radius, m b is the mass of a beam particle, and is the average beam plasma frequency squared.The normalized tune depression = 0 is defined by where !f const is the transverse frequency associated with the applied focusing field in the smooth-focusing approximation.
We now briefly illustrate the physical mechanism for the electrostatic Harris instability in intense particle beams with a Kapchinskij-Vladimirskij (KV) distribution.As shown in previous studies [5][6][7][8], the dipole mode has the highest growth rate, and for T jjb 0 the growth rate is an increasing function of k z r b and approaches a maximum value for k 2 z r 2 b 1.Therefore, we consider dipole-mode perturbations with k 2 z r 2 b 1, which in lowest order correspond to a displacement of the beam charge mainly along the beam propagation direction, arranged as a dipole-mode perturbation as shown in Fig. 2. One can distinguish three possibilities, !pb !?, ! pb !?, and !pb !? .If !pb !?, the charge oscillation will be mostly along the beam propagation direction due to the electrostatic restoring force, and as the result, the mode frequency will be close to the plasma frequency !!pb .In the opposite case, when ! pb !?, the charge perturbations will oscillate with the frequency !?, mainly in the direction transverse to the beam propagation direction due to the restoring betatron force.In this case, the mode frequency will be close to the average betatron frequency !? .Finally, if ! pb !?, the charge perturbation, moved longitudinally by the electrostatic restoring force, will at the same time traverse the beam transversely and will enhance the already existing perturbation as shown in Fig. 2(c).In this case, the arrangement of the charge perturbation does not change and the mode is purely growing with Re! 0 and Im! !? .One can also see the effects that longitudinal temperature has on the instability.If, during the period of transverse oscillation 2=! ?, a particle with average velocity v th jjb travels a distance larger than the wavelength of the perturbation z , then the perturbation enhancement shown in Fig. 2 will not occur and the instability is absent.This provides the threshold condition for the onset of instability, i.e., where in the inequality equation we have used Eq. ( 1).Since we have assumed k 2 z r 2 b 1, Eq. ( 4) implies We can now summarize the necessary conditions for instability as follows.The instability exists for sufficiently intense beams !pb !?with large temperature anisotropy T jjb =T ?b 1. Unstable modes have short wavelengths with k 2 z r 2 b 1.In what follows, we give a more quantitative description of the instability mechanism following the physical picture described above [6,20].Since a KV beam has a uniform density profile, the perturbed electrostatic potential for a dipole-mode perturbation inside the beam has the form where is the perturbation amplitude, and ! and k z are the perturbation frequency and the longitudinal wave number, respectively, and the perturbed electric field is given approximately by E ÿik z êz .Here x is the transverse distance from the beam axis in x direction (for example).
Next, we consider a beam particle oscillating longitudinally in the perturbed electric field and at the same time performing transverse betatron oscillations.For general distribution function, the equilibrium self-electric field is nonlinear, and the transverse betatron oscillation of the beam particles will contain many harmonics of the betatron frequency, which generally depend on the particle energy and angular momentum [5][6][7].For purposes of illustrating the physical mechanism, we consider here the simplified model of an equivalent KV beam where all of the particles oscillate with the same frequency, equal to the average depressed betatron frequency !?defined in Eq. (1), i.e., where 0 is the oscillation phase at t 0, x 2H x =m b p =! ? is the oscillation amplitude, and H x is the transverse energy.
Making use of Eqs. ( 5) and ( 6), the longitudinal equation of motion for a beam particle becomes Integrating Eq. ( 7) with respect to time t, we obtain where 0 !?t.To calculate the average displacement hzix; z; t in the z direction we average over all particles with the same transverse position x and longitudinal position z z 0 v z t at time t.This gives hzix; z; t where E z ÿik z [see Eq. ( 5)].Note that even though the individual particle motion [Eq.( 8)] has two characteristic frequencies, !ÿ k z v z ÿ !?and !ÿ k z v z !?, the average quantity hzi oscillates at the perturbation frequency ![Eq.( 9)].The beam displacement creates a restoring electric field Combining Eqs. ( 9) and ( 10), we obtain the simple dispersion relation [6,20] where has been made use of the average value of the plasma frequency introduced in Eq. ( 2) to take into account the beam density profile shape in a lowest-order sense.For a longitudinally cold beam, using the definition of the depressed tune [Eqs.( 1) and ( 3)], we can rewrite Eq. ( 11) as where n = 0 is the normalized depressed tune, and !n !=! f is the normalized mode frequency.Equation ( 12) is easily solved to give From Eq. ( 13), we determine that the mode with the lower sign in Eq. ( 13) is unstable and purely growing for 0:39 occurring for max = 0 1=3 p 0:58.Also, for very intense beams with = 0 !0, the normalized growth rate becomes Im!=! f ' = 0 .Figs. 3(a) and 3(b) show the normalized growth rate Im!=! f and real oscillation frequency Re!=! f at maximum growth, plotted as functions of the normalized tune depression = 0 obtained numerically using the eigenmode code bEASt [6].The thick solid curve corresponds to the simple estimate in Eq. (13).Despite the approximations made in the present simplified model, the agreement is reasonably good.
To illustrate the influence of the finite longitudinal temperature on the instability we assume that longitudinal velocity distribution in the beam frame is Maxwellian with Here, T jjb m b v 2 jjb .Substituting Eq. ( 14) into Eq.( 11) and integrating, we obtain the dispersion relation where Zs is the plasma dispersion function [22].Using Eq. ( 15) we can determine the finite bandwidth introduced by the finial longitudinal temperature.Indeed, for k z r b FIG. 3. Plots of the normalized growth rate Im! max =! f and real frequency Re! max =! f at maximum growth versus normalized tune depression = 0 for T jjb =T ?b 0 and azimuthal mode number m 0 (the two thin solid curves correspond to two different unstable modes) and m 1 (dotted curve).Results are obtained using the eigenmode code bEASt [6].The thick solid curve corresponds to the simple estimate in Eq. (13).
k z r b max the growth rate and the real frequency are Im! 0 and Re! 0 and we can rewrite Eq. ( 15) as Figure 5 shows the threshold value of the anisotropy T jjb =T ?b as a function of the normalized tune depression = 0 obtained using the eigenmode code bEASt [6].Note from Fig. 5 that the maximum threshold value, T th kb =T ?b 0:11, is achieved for moderately intense beams with = 0 0:4.

III. NONLINEAR SIMULATIONS OF THE TEMPERATURE ANISOTROPY INSTABILITY WITH BEAM EQUILIBRIUM STABILITY AND TRANSPORT (BEST) CODE
To investigate the nonlinear stage of instability, we make use of the nonlinear f method [21] described below, as implemented in the beam equilibrium, stability and transport (BEST) code [1,9,11].Here we give a short description of the f method.For a more detailed description of the BEST code the reader is referred to the extensive discussions in Refs.[1,5,7,9,11].
In the f approach, the solutions to the nonlinear Vlasov-Poisson equations are expressed as f b f 0 b f b , 0 , where f 0 b ; 0 are known equi- librium solutions (@=@t 0).In the electrostatic approximation the self-electric-field is represented as E s ÿrx; t.The perturbed potential satisfy the equation [9] where e b is the particle charge, and f b x; p; t is given by the weighted Klimontovich representation, Here, N sb is total number of beam simulation particles, N b is total number of actual beam particles, and the weight function is defined by w b f b =f b .The nonlinear particle simulations are carried out by iteratively advancing the particle motion, including the weights they carry, according to [9] and updating the fields by solving the perturbed Poisson's equation with appropriate boundary conditions at the cylindrical, perfectly conducting wall at radius r w .The f approach is fully equivalent to the original nonlinear Vlasov-Poisson equations, but the noise associated with representation of the background distribution f 0 b in conventional particle-in-cell (PIC) simulations is removed.In the f approach, the simulation particles are used to represent only a small part of the entire distribution f b f b ÿ f 0 b , and therefore the statistical error in the simulations is proportional to f w bi = N sb p , whereas the error in PIC simulations is proportional to pic 1= N sb p .Therefore, the typical gain in accuracy in f simulations compared to PIC simulations with the same number of particles is f = pic w bi [9].This fact allows much more accurate simulations of the nonlinear dynamics and instability thresholds when j w bi j 1.In addition, the f method can be used to study linear stability properties, provided all nonlinear terms in the dynamical equations ( 20)-( 23) are neglected [9][10][11].This corresponds to replacing the term 1 ÿ w bi with 1 in Eq. ( 22) for the weights, and moving particles along the trajectories calculated in the unperturbed potential 0 .
The f method described above has been implemented in the three-dimensional electromagnetostatic particle-in-cell code (BEST) in cylindrical geometry with a perfectly conducting cylindrical wall at radius r w .Poisson's equation ( 18) is solved using fast Fourier transform (FFT) techniques in the longitudinal and azimuthal directions.The particle positions [Eqs.(20) and ( 21)] and weights [Eq.( 22)] are advanced using a second-order predictorcorrector algorithm.The BEST code is parallelized using the message passing interface (MPI) with domain decomposition in the direction of beam propagation.The NetCDF data format is used for large-scale diagnostics and visualization.Typical simulation runs consist of 10 7 simulation particles and are performed on the IBM SP/RS 6000 at NERSC.
Figures 6 -11, obtained using the BEST simulation code, illustrate the nonlinear stage of the instability [20] for a beam with initial temperature ratio T jjb =T ?b 10 ÿ4 .Figure 6 shows the time history of the two azimuthal components of the density perturbation corresponding to m 0 (red), and m 1 (black) for an intense beam with = 0 0:6.As can be seen from Fig. 6, the m 1 mode dominates during the linear stage of the instability, whereas the m 0 mode becomes dominant after the instability has nonlinearly saturated.During the linear stage of the instability, the two modes have frequencies !0 for m 1, and !0:55! f for m 0. During the nonlinear stage, the frequency of the dominant m 0 mode is almost doubled to ! ' !f .Figure 7 shows the Fourier spectrum of the density perturbation n b plotted versus k z r w at two different times corresponding to the end of the linear stage (! f t 50), and the fully developed nonlinear stage (! f t 150) for an intense beam with = 0 0:6.Initially, the instability is excited over a wide range of longitudinal wave numbers k z r b > 1 (see also Fig. 4).During the nonlinear stage, the spectra shifts into the long-wavelength region with k z r b ' 2. Figure 8 shows the behavior of the effective longitudinal temperature T jjb m b hv 2 jj i as a function of time for a beam with = 0 0:6.Comparing with Fig. 6, note that the longitudinal temperature stays almost constant during the entire linear stage of the instability, growing superexponentially during a short time period 50 < !f t < 70, and then slowly increasing until reaching a constant level T jjb =T ?b ' 0:088 at time !f t 150.The temperature behavior is better illustrated by examining the longitudinal phase space of the beam.Figure 9 shows the contour plots in longitudinal phase space for a beam with = 0 0:6 at times (a) !f t 65, (b) !f t 75, (c) !f t 100, and (d) !f t 150.Since the linear stage of instability is dominated by the m 1 mode with zero real frequency, and a wide range of longitudinal wave numbers, at time !f t 65 the entire distribution is trapped inside a region of velocity space which corresponds to a spread in particle longitudinal energies Here, is the linear growth rate for the m 1 mode, and max denotes maximum value over the entire unstable spectrum of the m 1 mode for T jjb 0, and for a particular value of normalized depressed tune = 0 .During the time interval 50 < !f t < 70, particles are accelerated longitudinally and fill the trapped region of phase space.During the acceleration phase, the particle velocities are randomly scattered by the many waves excited, which lead to a quasilinear mixing in phase space and longitudinal thermalization with T jjb =T ?b ' 0:075 [Fig.9(b)] by time !f t 75.This longitudinal temperature can be estimated from Eq. ( 24) as T jjb ' E jj .During the remaining time, the nonlinear interactions lead to more mixing [Fig.9(c)] and small longitudinal heating to T jjb =T ?b ' 0:088 [Fig.8], and an eventual shift of the oscillation spectrum into the m 0 mode and into longer wavelength modes with k z r b 2 [Fig.9(d)].The final distribution has the clear signature of an m 0 mode with almost single-mode excitation at wavelength ' r w 3r b and frequency !' !f .The nonlinear evolution leads to a doubling of the mode frequency, and as a result to a shift of the wave phase velocity into the region of the phase space where the wave is stable.
The shape of the final longitudinal velocity distribution and the properties of the nonlinear mode at saturation depend on the beam intensity.Figure 10 shows the wings of the longitudinal phase space in the vicinity of the velocity resonant with the wave for two values of the normalized depressed tune, = 0 0:6 [Fig.10(a)] and = 0 0:7 [Fig.10(b)].Because of the reduction of the m 1 mode linear growth rate [Fig.4], and the increase of the longitudinal wave number where =k z is maximized, the longitudinal width in the latter case is reduced according to Eq. (24).Also the amount of particles trapped inside the m 0 mode structure is larger for the beam with depressed tune = 0 0:7.Evidently, the shape of the longitudinal velocity distribution is not Gaussian, but has a triangular shape for = 0 0:6 [Fig.11

IV. CONCLUSIONS
To summarize, we have generalized the analysis of the classical Harris instability to the case of a one-component intense charged particle beam with anisotropic temperature.The instability is kinetic in nature and is due to the coupling of the particles' transverse betatron motion with the longitudinal plasma oscillations excited by the pertur-bation.For a long, coasting beam, the f particle-in-cell code BEST and the eigenmode code bEASt have been used to determine detailed 3D stability properties over a wide range of temperature anisotropy and beam intensity.We have also used the f particle-in-cell code BEST to study the nonlinear evolution and saturation of the instability.The nonlinear saturation is governed by longitudinal particle trapping by a spectrum (not just a single wave) of fastgrowing waves with a broad band of longitudinal wave numbers and zero oscillation frequency.The presence of many waves leads to the nesting and overlapping of particle resonances in longitudinal phase space, and as a consequence, to the fast randomization of the trapped-particle distribution and longitudinal heating of the beam particles.The nonlinear interactions lead to the shift of the waves spectrum into the long-wavelength region (m 1 !m 0 and k z r b !2).This can be explained qualitatively as follows.The instability is self-quenching, since as the instability grows, the plasma heats up longitudinally due to particle trapping and phase-mixing.The large wave number modes saturate first, but as the smaller wave number modes grow and saturate they increase the longitudinal temperature even further [see Eq. ( 24)], which in turn begins to damp the saturated modes at large wave numbers.As a result, the wave energy spectrum transfers to the smaller wave number region.This is qualitatively consistent with the behavior of the mode spectrum observed in the simulations.The final longitudinal velocity distribution is not Maxwellian and can be characterized by a remnant temperature anisotropy T jjb =T ?b ' 10%, where T jjb m b hv 2 jj i.A very interesting feature of the nonlinear saturation is the formation of a stable, longitudinal, nonlinear BGK-like wave structure with a negligible number of trapped particles.

FIG. 1 .
FIG. 1.Comparison between Harris instability in (a) a neutral plasma n 0e n 0 i confined by an axial magnetic field, and (b) for a nonneutral ion beam n 0 b Þ 0; n 0 e 0 confined transversely in a focusing channel.

FIG. 2 .
FIG. 2. Physical mechanism for the electrostatic Harris instability in intense particle beams for the case of a dipole-mode perturbation.Three possibilities are illustrated, (a) for !pb !? the mode frequency is !!pb ; (b) for !pb !? the mode frequency is !!?; and (c) for !pb !? the mode is purely growing Re! 0 and Im! !? .

dt expt 2 ÿ 1 ;Figures 4
Figures 4(a) and 4(b) show the results obtained from the eigenmode code bEASt [6] for the normalized growth rate Im!=! f plotted as a function of the normalized longitudinal wave number k z r b for = 0 0:3 and for several values of the temperature ratio T jjb =T ?b 0; 0:01 and 0:05.Here, Figs.4(a) and 4(b) correspond to azimuthal mode number m 0 and m 1, respectively.From Fig. 4(b), the instability bandwidth is k z r b max ' 11; 4, whereas Eq. (17) predicts values k z r b max ' 10; 4:5 for the temperature ratios T jjb =T ?b 0:01; 0:05 respectively.Despite the approximations made in the present simplified model, Eq. (17) gives a reasonably accurate prediction of the instability bandwidth.As the temperature anisotropy is reduced T jjb =T ?b ! 1 the bandwidth gets smaller, until k 2 z r 2 b max ' 1, and the instability disappears.The threshold value of the anisotropy T jjb =T ?b required for the onset of instability cannot be determined from the simplified model because of the assumption k 2 z r 2 b 1.Figure5shows the threshold value of the anisotropy T jjb =T ?b as a function of the normalized tune depression = 0 obtained using the eigenmode code bEASt[6].Note from Fig.5that the maximum threshold value, T th kb =T ?b 0:11, is achieved for moderately intense beams with = 0 0:4.

FIG. 4 .
FIG. 4. Plots of the normalized growth rate Im!=! f versus k z r b for = 0 0:3 and several values of the temperature ratio T jjb =T ?b 0; 0:01; 0:05.Here, parts (a) and (b) correspond to azimuthal mode numbers m 0 and m 1, respectively.Results are obtained using the eigenmode code bEASt[6].

FIG. 6 .
FIG. 6. (Color)The time history of the two azimuthal components of the density perturbation is shown for m 0 (red), and m 1 (black) for a beam with = 0 0:6.
Figures 11(a) and 11(b) show the z-averaged longitudinal velocity distribution for two values of normalized tune depression, = 0 0:6 and = 0 0:7, respectively, plotted at different times during the nonlinear evolution of the instability.The average longitudinal velocity distribution in both case becomes almost stationary after time !f t > 120, notwithstanding the obvious mode structure that is still clearly evident in the phase-space contour plots [Figs.10(a) and 10(b)].

FIG. 7 . 7 FIG. 10
FIG. 8. Plot of effective longitudinal temperature T jjb m b hv 2 jj i as a function of time for a beam with = 0 0:6.