Electron cloud instabilities in the Proton Storage Ring and Spallation Neutron Source

Electron cloud instabilities in the Los Alamos Proton Storage Ring (PSR) and those forseen for the Oak Ridge Spallation Neutron Source (SNS) are examined theoretically, numerically, and experimentally.

A transverse, bunched beam instability has been seen in the PSR. There are several curious features, described in section 2, which strongly argue that the instability is due to coupled electron-proton oscillations. Electron cloud driven instabilities were observed in the ISR [5,6] and are known or suspected in several positron rings. These include the KEK photon factory [7], CESR [8], KEKB [9,10], and PEP-II [11]. For these machines the bunch length is relatively short, and much theoretical and numerical work has been done to explain the observations [7,[12][13][14][15][16]. They have also been seen in the PS [17] and SPS [18] with proton beams. For the PSR and SNS the bunches are long and a natural starting point is the application of coasting beam dispersion relations [5,6,[19][20][21]. These dispersion relations have been applied to the PSR [1,22,23] and, recently, to the proposed J-PARC project [24]. Simulations of both coasting [29,30] and bunched beams [24][25][26] have been started. Studies based on beam break-up models [27] as well as studies of behavior well beyond threshold [28] have also been done. This paper aims to provide an understanding of and insights into the electron cloud instability in the PSR that will be useful for a timely estimate of the impact on SNS. Parameters for the two machines are shown in * Electronic address: mmb@bnl.gov  Table I.

II. OVERVIEW OF EXPERIMENTAL DATA
The experimental characteristics of the PSR instability may be summarized as follows: 1. The central frequency of the instability f c increases with intensity as f c ∝ √ I [1].
2. The PSR instability is controlled, in practice, by applying a sufficiently high RF voltage V rf . For fixed bunch length τ b the maximum number of stored protons N b scales linearly with V rf [3]. 3. The threshold value of the RF voltage for a given intensity increases when some unchopped beam is injected into the gap [3]. 4. A broad band tranverse resistance of order 1 MΩ/m is needed to match the observed growth rates [1]. 5. For a fixed RF voltage the maximum number of stored protons depends only weakly on bunch length [1,31] . 6. Sustained, coherent oscillations are observed below the loss threshold for a well-conditioned ring 7. At moderate to high intensities an intense electron flux at the wall is observed as the bunch passes [32]. 8. Over accumulated operating time, for a fixed bunch intensity, the threshold RF voltage decreases [33]. 9. Increases in vacuum pressure and losses have a marginal effect on the stability threshold.
The first two items are difficult to reconcile with an impedance-driven instability since a fixed impedance should drive a given range of frequencies and given a fixed impedance the threshold intensity should scale linearly with momentum spread (or synchrotron frequency) and hence as V rf Item 3 would be relevant to an impedance driven instability if the beam in the gap was adequate to keep the offending resonator driven at a sufficient level. Data show that the threshold voltage doubles when ∼ 3% of the injected turns are unchopped. In this case the ratio of the instantaneous current in the middle of the gap to the peak current is less than 2%.
Item 4 is difficult to reconcile with a machine circumference of C = 90 m and an average pipe radius of 5 cm.
Item 5 is at odds with both narrow and broad band impedance driven stability models, since reducing the bunch length reduces the momentum spread in the beam.
Item 6 suggests that there is a difference between the linear threshold of the instability and the intensity required for beam loss.
Items 7 and 8 [23,32,34] are consistent with the model that has evolved to explain the previous 6, namely that the instability is driven by electrons [1][2][3][4]. Fig. 1 shows the threshold RF voltage as a function of beam current for different dates and a variety of injected bunch lengths. The injected bunch length is set by the width of the chopper pulse and, since the linac beam has non-zero momentum spread, is a lower limit on the bunch length in the ring. The curves are nearly straight (item 2) and independent of bunch length (item 5). The historical threshold voltage lies above the curves obtained later in the run (item 8). The tune of the machine is quite stable, and the reduction in threshold voltage over time might be due to a reduction in secondary yield due to conditioning by electrons.
Evolution of the instability is shown in Fig. 2. Notice that the electron signal peaks after the instability has risen, suggesting that the electron cloud responsible for starting the instability is not present near the detector. Other electron detectors show an observable signal before the instability begins. The electron signal represents the electron flux into the wall and gives an indication of the cloud density. This will be addressed in Section III. Fig. 3 shows the beam current and a mountain range plot of transverse position along the bunch. The position was obtained by integrating the sum and difference pickup signals and taking the ratio. The signal grows from noise but the amplitude saturates. These data are for a full bunch length of 275 ns and V rf = 15.3 kV. Data have also been taken for shorter bunches. Persistent oscillations are present as well. The oscillations set in when V rf is about 10% above the fast loss threshold for bunch lengths between 200 ns and 275 ns. along a PSR bunch (red) and the beam current(blue). The total time was 490 µs, increasing vertically, with every 30th turn plotted. The RF voltage was 10% above the threshold for fast losses.

III. FORMATION OF THE ELECTRON CLOUD
The properties of the electron cloud are fundamental to the question of beam stability. An inital electron population is created by beam loss, residual gas ionization, and various processes related to the stripping foil [1,2,35]. The inital population interacts with the beam and, via secondary emission processes, is amplified. For short bunches the problem has been addressed by several authors [12,[36][37][38]. For PSR and SNS the bunches are long and an electron trapped by the beam performs > ∼ 25 transverse oscillations during a single bunch passage [39]. For electrons that exist within the beam pipe before the bunch arrives, and for those created by residual gas ionization, the electron transverse amplitude remains small during the bunch passage and multiplication via secondary emission is unlikely. On the other hand, free electrons created by losses have a transverse amplitude equal to the beam pipe radius when the instantaneous beam current is large and strike the beam pipe with appreciable energy. These processes are illustrated in Fig.  4 for typical PSR parameters.
The dynamics illustrated in Fig. 4 are amenable to both analytic estimates and simulations. First we estimate the energy with which an electron can strike the wall. The electron equation of motion is approximated by where y is the vertical position of the electron, I(t) is the instantaneous proton beam current, β = v/c for the protons, Z 0 = 377 Ω, and σ is the rms radius of the proton beam. Here we neglect the electric field due to the other cloud electrons. We consider motion with an amplitude comparable to the beam pipe radius b and introduce an effective frequency ω e (t) so that the oscillator equation y + ω 2 e (t)y = 0 has the same frequency as equation (1). To estimate ω e take ω e (t) = eZ 0 I(t) 2πβm e (σ 2 + 2b 2 /π) , which interpolates between small and large amplitude motion, and agrees with the numerical integration of equation (1) to < 5%. For the harmonic oscillator the adiabatically conserved action is J = E/ω e (t) where E is the total energy. Now suppose a situation similar to the red electron line in Fig.  4. At t ≈ 0.2 µs the electron nearly hits the lower wall at y = −b. Half an oscillation period later the electron bounce frequency has changed by ∆ω e and the electron would reach an amplitude of y = b(1 − ∆ω e /2ω e ) if the upper wall was absent (∆ω e < 0 for trailing edge electrons). The increase in amplitude yields the energy with which the electron strikes the wall. Substituting units, the electron strikes the wall with an energy The strike energy is positive on the trailing edge of the bunch (ω e < 0) and equation (3) is valid only when the electron frequency does not have a large fractional change per period, |ω e | ω 2 e . Equations (2) and (3) predict E strike = 55 eV for the strike energy. From the PSR simulation the electron has an energy of 45 eV for the first wall strike. Notice that a small error in vertical steering will cause the grazing and strike to occur on the same side of the beam pipe. For this case the strike energy is ∼ 100 eV hence amplification of the electron cloud by secondary emission is even more likely. In comparison, typical numbers for ISIS give E strike < ∼ 10 eV for a centered beam hence no secondary amplification is expected. Secondary emission involves the collisional liberation of electrons from within the first few tens of nanometers of the surface. Other processes such as backscatter and rediffusion are also present. A useful experimental measure is the secondary emission yield δ(E). For a normally incident electron with kinetic energy E, δ(E) is the average number of electrons leaving the surface due to all processes. Fig. 5 shows δ(E) for titanium nitridecoated stainless steel. The samples were manufactured by BNL vacuum group under the direction of H.C. Hseuh. The measurements were performed by B. Henrist of the CERN vacuum group. Notice that δ(0) ≈ 0.8 and is due to backscattering. The peak yield varies by ±10% and the peak value of δ(E) for unbaked, unconditioned, uncoated stainless steel is around 2.5.
To model the effect of a given curve the seconday yield is parameterized as The terms proportional to R 0 and R inf in equation (4) approximate the contribution from reflected electrons while the last is due to true secondaries. For non-normal incidence the entire secondary yield is multiplied by where cos θ = 1 for normal incidence and α θ ≈ 0.5. Electrons emitted off the surface have a range of energies and angles. For true secondaries the energy distribution is modeled by Elastically reflected electrons have an energy equal to the incident energy and rediffused electrons have an energy distribution which is uniform between zero and the incident energy. The conditional probability for rediffusion P r,c is an input parameter.
This model hase been implemented in the code CSEC (Cylindrically Symmetric Electron Cloud). With the wall conditions specified the simulation proceeds as follows. Initially, electron macro-particles are created at the wall. These particles are actually parallel filaments which are long compared to the pipe radius and parallel to the pipe axis. Only transverse motion is considered. At each time step the electric fields are calculated assuming constant density in the longitudinal direction. The macroparticles are produced at regular time intervals (hundreds to thousands) along the bunch and the charge on each macro-particle is proportional to to the instantaneous beam intensity. The average charge liberated per meter per second is an input parameter. The beam and pipe are round and the cloud field is calculated assuming cylindrical symmetry. However, individual electron macroparticles have both radial and azimuthal velocities to include the effects of angular momentum. A drift-kick algorithm is used to update the macro-particles. When a macroparticle hits the wall the reflection probability, , is calculated and a uniform random deviate r1 between 0 and 1 is chosen. If r1 < P r the macroparticle is reflected and its charge is unchanged. Another random deviate is calculated to choose between elastic and rediffused reflection and the velocity off the wall is obtained.
When r1 > P r the macroparticle charge is multiplied by The macroparticle charge is compared with upper (Q max ) and lower (Q min )values, which are input parameters. If the charge is less that Q min the macro-particle is dropped. If the macro-particle charge is larger than Q max the macroparticle is split into smaller macroparticles so that each has a charge smaller than Q max . For each of the macroparticles a random deviate r2 between 0 and 1 is chosen and the secondary energy is calculated by  The velocity of the macroparticle is set to the electron velocity corresponding to E. Given the speed of the particles off the wall, the angle between the wall normal and the velocity is chosen according to cos θ = √ r3 where r3 is a uniform deviate between 0 and 1, corresponding to an angular distribution dN/dΩ ∝ cos θ.
This simulation has been benchmarked against the code POSINST [40][41][42] for a PSR test case. Fig. 6 shows the total line density for both electrons and protons. The blue line is the output of CSEC. The red line is POSINST output using the full model for the energy spectrum of the rediffused electrons. Fig. 7 shows the electron current into the wall for CSEC and POSINST. For some parameters a 50% discrepancy is apparent but in the following we will show that this is small within the present context. Table II shows the results of nonlinear least squares fitting of equation (4) to the curves in Fig. 5. For all but one case the best fit asymptotic reflection probability was negative, so we set R inf ≡ 0. The worst fit, with an rms error of 0.025, is shown in Fig. 5. The other SNS simulation parameters are an rms beam radius of   Figure 8 is the most relevant for stability analysis. To see this consider a round, uniform beam with line density λ p and radius a p , and a uniform electron cloud with line density λ e and radius a e . Let y be the (small) offset between the beam and cloud centroids. The force per unit length on the proton beam is .
The last term on the right of (5) is the electron line density within the beam and, since a p is nearly constant   Figure 5. 24  700  43  1  36  350  39  2  61  41  10  3  24  700  45  5  49  110  22  6  56  78  17 along the beam, a linear response model should use the electron line density within the beam. Fig. 9 shows the electron current into the wall. When integrated with respect to time one gets the electron dose per unit area. Table III gives the total electron dose as well as the average electron strike energy and the dose with strike energy greater than 100 eV.
Over time, electrons striking the wall reduce the secondary yield. This process, known as conditioning, is a standard tool in RF and microwave engineering. Conditioning using the beam has proved beneficial in B factories, light sources, the SPS, and is important for the LHC. Figure 1 shows the benefit of conditioning for the PSR. Figure 10 shows data obtained from the articles by Henrist et al [43] and Kijima et al [44]. The Henrist data are for clean copper. The Kijima data are for a variety of materials. In particular, the curve for the OFHC surface was rinsed with water and then electron conditioned with no bake in between.
For all curves in Figure 10 an electron dose of 0.1 C/cm 2 reduces the peak secondary yield below 1.8. For all the SEY curves in Figure 5 the peak secondary yield is greater than 1.8 and one can reasonably expect that reducing the peak SEY to 1.8 would yield smaller densities than any shown in Figure 8. time required for such conditioning consider the dose with E strike > 100 eV in Table III. Typical values are > ∼ 15 pC/cm 2 /turn. Simulations show that this value varies by about 50% over the last 400 turns of the 1200 turn SNS cycle, so the average dose per SNS cycle is 4.5 nC/cm 2 /cycle. Dividing this into 0.1 C/cm 2 and taking a 60-Hz rep rate gives a conditoning time of 106 hours.
The SNS beam must be stable at reasonably high intensity for such conditioning to take place. In the following sections we consider what is known about the PSR and use various theoretical tools to extrapolate to the SNS parameter regime.

IV. LINEAR STABILITY THEORY
The stability eigenvalue problem for positron bunches in the KEKB-LER and protons in the CERN-SPS has been considered in [15,16]. A key parameter is the number of electron oscillations during the bunch passage, which is fairly small for both machines. In this case one can get good results using traditional basis expansion techniques [45] while maintaining a managable matrix size.
For the PSR and SNS consider the beam position monitor traces in Fig. 2. The peak oscillation frequency is about f max = 200 MHz and the bunch length is about τ b = 200 ns. Using Perevedentsev's notation [16], the smallest matrix which would allow for this fidelity involves a sum over all and k ≥ 0 with | | + 2k ≤ 2f max τ b = 80. This yields a 3321 × 3321 matrix. The eigenvalues would need to be calculated and the matrix size increased until convergence was found. Note that these matrices are dense with no special symmetry properties and a generic routine is needed [47].
An alternate technique is to model the longitudinal dynamics using a square well potential [46,48,49]. Particles coast longitudinally within the bunch and are reflected at both the head and tail. The accompanying line density is a boxcar distribution and the longitudinal velocity distribution is independent of position within the bunch. With this special distribution we will find that the dimension of the eigenvalue problem scales linearly with the number of electron oscillations within the bunch, which keeps the matrices manageable.
For zero chromaticity, the vertical equation of motion for a single proton is In equation (6) y p is the proton coordinate, θ is the machine azimuth which will be used as the timelike variable, Q 0 is the bare betatron tune, ∆Q sc is the incoherent space-charge tune shift, andȳ e is the centroid position of the electron cloud which depends on both θ and longitudinal position within the bunch. The linear approximation for the space charge force is equivalent to taking the Kapchinskij-Vladimirskij distribution for the transverse phase space density, and is discussed in the Appendix. The strength of the electron cloud interaction is determined by the parameter where λ e is the magnitude of electron line density, ω 0 is the angular revolution frequency, a e is the radius of the electron cloud, a p is the radius of the proton beam, γ is the relativistic factor of the beam and m p is the proton mass. Notice that Q p is the betatron tune the protons would have in the absence of other focussing forces. Let the longitudinal coordinate within the bunch be φ = ω 0 τ with 0 < τ < τ b . As the bunch passes the electrons remain at a fixed value of θ and the electron centroid obeys where the boundary conditions areȳ e (0, θ) = 0, and ∂ φȳe (0, θ) = 0. In equation (8) the transverse electron oscillation frequency is given by in analogy to (7), and the spread in electron oscillation frequency is parameterized by α ≡ Q e /2Q r with Q r being the effective quality factor of the electron cloud. Factors contributing to this spread include variations in a p , via variations in the lattice functions, and the dependence of the electron frequency on the electron's amplitude. Taking only lattice variations we estimate Q r ≈ 2 for both PSR and SNS. Integrating equation (8) one obtains the electron centroid.
Given equations (6) and (10)  The equation for y p is found by making the substitution in equation (6), where U (φ) is the longitudinal potential associated with the square well. For parameters relevant to a synchrotron, the betatron tune shifts will be small compared to the betatron tune hence we can define y p (φ, v, θ) = Y (φ, v, θ) exp(−iQ 0 θ) and neglect second derivatives of Y . One obtains whereȲ with ρ(v) the normalized velocity density, Also, the approximation Q 2 0 ≈ Q 2 0 + Q 2 p has been used in calculating the coefficients on the right hand side of (11).

A. Action-Angle approach
The Hamiltonian for longitudinal motion is H = v 2 /2 + U (φ) where U (φ) is a square well with walls at φ = 0 and φ =φ [50]. As a proton executes a single longitudinal oscillation it encloses a phase space area of A = 2|v|φ. Define the action-angle variables I and ψ. Since I is constant and ψ increases by 2π during each oscillation, A = 2πI. In the action-angle variables the longitudinal Hamiltonian is H = (πI/φ) 2 /2. With this Hamiltonian dψ/dθ ∝ |v| so we expect φ =φŝ(ψ)/π where the sawtooth function is given byŝ(ψ) = |ψ| for |ψ| < π and for other valuesŝ(ψ) =ŝ(ψ + 2π). A canonical transformation of Goldstein's third type may be used to verify the coordinate change [51]. (∆Q + ∆Q sc )Y (ψ, The total wake potential is given by (13) Next expand Y (ψ, I) as Use Fourier orthogonality to isolate Y n (I) and definê Since φ depends only on ψ and not I the second line of (12) depends only on the values ofŶ n . Isolate the values ofŶ n on the first line of (12) and define where δ k,0 is the Kronecker delta. This yields the final dispersion relation, The dispersion integral is given by with ν = ∆Q + ∆Q sc and the impedance matrix is where This has the form of a coasting beam dispersion relation. Consider a wave with frequency ω e . The value of k which creates modulations at this frequency is k = ω e τ b /π. Substituting this value and generalizing to an arbitrary coherent tune shift ∆Q 0 yields With ∆Q 0 = ∆Q sc this is the dispersion relation for space charge waves of frequency Q e on a coasting beam with current equal to the peak current in the bunched beam.
For definiteness take ∆Q 0 = ∆Q sc and a parabolic momentum distribution with ρ(v) = (3/4v)(1−v 2 /v 2 ) for |v| <v. Equation (18) predicts a threshold for coherent oscillations which is given by Consider the PSR threshold for N p = 4×10 13 with V rf = 15 kV and τ b = 220 ns. Substituting these values in the last equation one obtains F = 1.52. Using design parameters for the SNS at 2 MW one obtains F = 0.42.
The space-charge threshold also scales differently from the observed behavior in PSR. Assuming that the bunch shape remains constant (ie, parabolic ) as the RF voltage and bunch length are changed, the space charge threshold is given by where K depends one beam energy, betatron tune and other machine constants. The maximum number of stored protons scales linearly with RF voltage but the scaling with bunch length does not agree with the observations. Up to this point the calculation is equivalent to that presented in [49]. We go on to include the effects of the electron forces. Toward this end consider velocity distributions of the form [52] where α j+1 > α j > 0, and the C j s are obtained using simple matrix techniques. The shape of the distribution depends on the values of α j chosen, with distributions approaching Gaussians being easy to construct. One such formula is For M > 10 the fine cancellation implied by equation (20) appears to require more than 16 bytes of numerical precision and we will generally use M = 5. The dispersion integrals (15) are given by To compare these results to those of the Gaussian distribution consider the coasting beam dispersion relations ∆Q 0 = 1/D(ν, 1) whereφ = π and ν varies over the real numbers. Fig. 11 shows threshold curves for M = 5, parabolic, and Gaussian distributions with σ(v) = 1.

Substituting expression (21) for the dispersion integrals into the eigenvalue problem (14) yields
To proceed set A m,j . Insert this expression into the previous matrix and demand equality for each value of j, without the sum. This results in For practical applications the infinite sum over p needs to be truncated and a value of M chosen. We will take M = 5 and test convergence with a simple model. Set ω 0 τ b = 2π/35 and ω e τ b = 5π. The beam current and electron line density are chosen to give a cold, coasting beam tune shift of Plots of Im(∆Q) for the most unstable mode versus the rms tune spread in the cental line (Q e σ(v)) for various truncation values as well as those obtained from equation (18) are shown in Figure 12. Notice that the coasting beam dispersion relation gives a reasonable threshold estimate for the eigenvalue problem using the larger number of modes, 0 ≤ p ≤ 105. Also notice that the red line takes a dip in the vicinity of Q e σ(v) = 0.015. This dip is quite important since an accurate solution of equation (11) must approach the beam breakup limit as σ(v) → 0. For the beam breakup limit the amplitude of the oscillation grows as ln |y| ∝ √ θ and as the beam breakup limit is approached the growth rates of all modes go to zero [46, Fig. 12].
For parameter regimes appropriate to the PSR and SNS, ω e τ b /π 1 and we will sum for |p − ω e τ b /π| < N . Also, for comparison to a real beam, appropriate bunch lengths and momentum spreads need to be obtained. For this we equate rms bunch length the rms momentum spread in the real and modeled beams. For the PSR  a typical, real bunch length is 250 ns. When modeled by a square bunch we set τ b = 177 ns. To avoid confusion we will refer to the bunch length of the square-modeled bunch as the bunch length. Fig. 13 shows results for a bunch charge of 6.4 µC and a range of PSR parameters. The threshold voltage increases with increasing electron density, and also increases as the bunch gets shorter. The latter observation is consistent with equation (19). We note that Im(∆Q) = 0.015 gives an equivalent transverse resistance of about 1 MΩ/m. Fig. 14 shows the threshold RF voltage versus space charge tune shift for a bunch length of 177 ns and a bunch charge of 7 µC. The plot was made by varying ∆Q sc in equation (13) while leaving the second term in the wake- field fixed. Space charge clearly increases the threshold voltage which also is consistent with (19). This result is different from the increase in threshold current with space charge tune shift reported in [46]. The set of wake potentials used in [46] led to coupling between low lying modes, even with large space charge tune shift. The wake potentials used here change sign several times allowing for coupling between high order modes. Since space charge reduces the distance between high order modes, space charge reduces the wake induced tune shift needed for coupling. Figure 15 shows thresholds for various electron line densities and compares the matrix analysis with the coasting beam estimate. The good agreement between  (18), and a parabolic momentum distribution with the same rms width as for Fig. 15. The threshold curve for historical PSR data is shown for comparison. For the data, the momentum spread from the LINAC is present even with zero RF voltage.
the two techniques suggests that a coasting beam estimate based on a more accurate momentum distribution (ρ(v)) would yield a better estimate of behavior in the PSR. Fig. 16 shows threshold estimates for a parabolic distribution. While the threshold estimates for small beam current are comparable in Figures 15 and 16, the threshold voltage increases more rapidly with beam current in Fig. 16. This is a consequence of the dispersion diagrams shown in Fig. 11. The M = 5 dispersion curve has extended "wings" while the dispersion curve for the parabolic distribution goes to zero before a real tune shift of 2σ. This implies that the parabolic threshold is more sensitive to the space charge tune shift, so a fair approximation to the threshold voltage can be obtained from equation (19). Also notice that equation (19) implies that threshold voltage increases as bunch length decreases. The factor of 3 discrepancy between theory and experiment shown in Figure 16 becomes worse as the bunch length is reduced. Also notice that the calculated threshold voltage is a weak function of the average electron line density, which is at odds with the conditioning effects implied by Figures 1 and 10.

B. Nonlinear space charge forces
For the PSR the transverse beam profile is roughly Gaussian. SNS target requirements dictate a transverse density that is nearly constant within the beam, and a sharp beam edge. This implies that the transverse amplitude dependence of the space charge tune shift in the PSR is significantly larger than in the SNS. The impact of this difference on instability threshold estimates for the PSR are the subject of the present section.
Consider a coasting beam instability and take the evolution variable to be time (t). Consider only one transverse dimension (y, v =ẏ). Take the fractional momentum deviation, δ = (p − p 0 )/p 0 , to be the longitudinal momentum-like coordinate and θ to be the longitudinal position coordinate. The equation for the electron centroid is approximated as For no momentum spread the proton centroid obeys Setȳ p =ŷ p exp(ik(ω 0 t − θ) + iω c t) andȳ e = y e exp(ik(ω 0 t−θ)+iω c t). Assume that |kω 0 −ω c −ω e | ω e /Q r so thatŷ e = −iQ rŷp . Additionally assume that |ω c − ω β | ω β so that coupling to the other betatron sideband can be ignored. Then ω c − ω β = iω 2 p Q r /2ω β is the coherent frequency shift of the protons. Next consider momentum spread in the absence of collective Threshold diagrams for various values of ∆Qsc,max/W and highly nonlinear space charge with q/κ = 3/5. Tune shifts corresponding to points below a given curve are stable for that value of ∆Qsc,max/W . forces. For a particle with momentum offset δ its kth betatron sideband occurs at a frequency ξ is the normalized chromaticity, and ω β is the betatron frequency for an on-momentum particle in the absence of collective effects [57]. The frequency spread in the beam is given by the momentum distribution, ρ(δ), and the relation δω β (δ) = η(kω 0 − ω β )δ + ξω β δ ≈ −ηω e δ.
We may now use the formalism in [52] to estimate the effects of nonlinear space charge on stability thresholds. Considering the results of the previous section we take a parabolic momentum distribution with half width at base δ max , and numerically calculate equation (33) in [52]. Figure 18 shows stability diagrams for a generic coherent tune shift and various amounts of space charge. The parameter W = δω β (δ max )/ω 0 is the normalized frequency spread in the sideband due to momentum spread, and the curves are labeled according to the central space charge tune shift via ∆Q sc,max /W . These curves assume a round beam with transverse density proportional to (r 2 max − r 2 ) 2 with nonlinearity parameter q/κ = 3/5 [52]. This regime corresponds to a soft upper limit for space charge tune spread in real beams and yields a tune shift with betatron amplitude (ŷ) that is given by For our purposes the width of the tune distribution is W = δ max |η|Q e . By any measure, ∆Q sc Q r Q 2 p /2Q β , so the coherent tune shift is ∆Q 0 = iQ r Q 2 p /2Q β . Figure 19 shows threshold voltage versus beam current for long and short bunchs with maximal space charge tune spread.
For both bunch lengths the threshold voltage with nonlinear space charge is no larger than the threshold assuming linear space charge forces. Also notice that the threshold voltage for the short bunch is significantly reduced when the electron line density is reduced. The value of 0.25 nC/m was in fact calculated using CSEC with a 7 µC bunch and wall parameters consistent with partially conditioned stainless steel. The large reduction in threshold voltage due to reduced λ e is in sharp contrast to the results in Figure 16.
The possibility that threshold voltage is a strong function of electron survival during the gap has been suggested before, and the conjecture that it influences the dependence of threshold voltage on bunch length has been made [1,31]. There are other experimental facts that have not been included in the calculations leading to the Figures. The electron line density surviving the gap is a strong function of bunch current [33]. Also, the threshold scaling like τ 0 b in Fig. 1 is different from earlier observations, which showed something more like τ −1 b [1,2]. That is to say, the fact that the curves for different bunch lengths are essentially identical in Fig. 1 does not indicate a fundamental symmetry of the physics. Finally, from Fig. 2 it is clear that the electron cloud responsible for the onset of the instability does not always occupy the entire ring.
By including these sort of effects it might be possible to fit both the intensity and bunch length scaling in PSR. Instead, we will remember that the model used here is fairly rough and attack the problem using simulations.

A. Description of algorithms
The parameter regimes appropriate to the PSR and SNS make a direct particle in cell calculation of the electron cloud instability difficult. In this section we develop some phenomenological equations of motion which contain the various interactions in simplified form that allow for greater computational speed. Consider a single proton macro-particle. The continuum version of its equations of motion are taken to be The longitudinal coordinate is arrival time, τ b , and the motion is simple harmonic with synchrotron tune Q s . The vertical motion has bare tune Q y , with transverse, linear space charge forces and a vertical force on the proton due to the electron cloud F y,e . The assumption of linear transverse space charge increases computational speed. From Figures 18 and 19 this approximation could significantly underestimate threshold current for PSR. The horizontal motion is unperturbed with betatron tune Q x . Neglecting the horizontal collective forces roughly halves the simulation time though they could easily be included. The code actually allows for non-zero chromaticity, but its effect is negligible for PSR and SNS parameter regimes. Electron macro-particles are assumed to have transverse motion only. When the macroparticle is inside the pipe d 2 y e dτ 2 = F y,p (y e , θ, τ ) + κF y,e (y e , θ, τ ) The vertical force on an electron macroparticle is due to both the protons F y,p and the electron cloud κF y,e . When nonzero, the parameter κ accounts for the mass ratios as well as the fact that the proton time-like variable is taken to be azimuth (θ), while the electrons evolve in real time during the passage of the proton bunch. Horizontal electron motion is analogous to the vertical for motion in a drift (I d = 1). For motion in a dipole magnet the horizontal/longitudinal Larmor motion is neglected and I d = 0. Wall interactions will be discussed later. The continuum version of the collective force on the electrons due to the protons is taken to be, .
(31) The actual beam radius a(θ) is allowed to vary with azimuth as appropriate to a strong focusing machine. This means that x p is actually the normalized position coordinate and θ is the phase advance divided by the tune. The proton beam is assumed to be round which simplifies the electron update equations. The value a(θ) is an input function, not derived from the beam characteristics, with a(θ) = a(θ + 2π). The range of values of a(θ) reproduce the variation in electron bounce frequency for the actual lattice calculated using variations in both the horizontal and vertical beam dimensions.
To get the collective force on the protons due to the electrons the average and mean square positions of the cloud are obtained for each azimuth and each time step along the bunch. Let the macroparticles be denoted by index j, then y e (θ, τ ) = j λ e,j y e,j j λ e,j The average and variance are then multiplied by a 0 /a(θ) and [a 0 /a(θ)] 2 , respectively, to account for the normalized proton coordinates. The macroparticles can have different charges (actually line densities) due to interactions with the walls. The force on the protons due to the electrons is The initial conditions of the electron cloud and secondary emission are considered as follows. For each azimuth at which the beam is updated, usually 10 times per betatron oscillation, the electron cloud is generated by taking N e electron macroparticles with random positions within the pipe and zero velocity. Each macroparticle has the same initial line density and their sum is an input parameter, λ e . A given macro-particle is evolved by equations (29) and (30) until it strikes the wall with energy E. Upon striking the wall equation (4) is evaluated and that macroparticles charge is multiplied by the secondary emission yield δ(E). After the charge is updated the macroparticle remains in the same location, but with zero velocity. This neglects the complications associated with the secondary energy distribution and should have a small effect on the proton dynamics.
To discretize the equations of motion we take ∼ 10 Q y equally spaced thin lenses in machine azimuth to implement the collective forces. Between collective kicks the transverse and longitudinal proton motions are approximated by rotation matrices. Consider the implementation of the collective forces at a given thin lens. Choose a longitudinal smoothing length τ e and a longitudinal bin size δτ = T rev /N g . Generally τ e /δτ > ∼ 10 and ω e τ e < 1, where ω e is the maximum value of the electron bounce frequency within the bunch. Use linear interpolation to calculate estimates of λ p (τ ) and λ p (τ )ȳ p (τ ) on the N g grid points. Smooth these arrays using S(t) = (1+4|t|/τ e ) exp(−4|t|/τ e )/τ e which involves one exponetiation and O(N g ) additions and multiplications when the convolution is expressed as an autoregressive filter. The smoothed values λ p (τ ) andȳ p (τ ) (obtained by division) are used to drive the electron cloud via equations (29), (30) and (31). As the proton bunch passes equations (32), (33) and the electron line density λ e (τ ) are stored in arrays. Then equation (34) is used to get the transverse kicks on each of the proton macroparticles.

B. Results
The simulation code has 5 purely numerical parameters and 14 physical parameters. Fiducial values for these are given in Table IV. To characterize the evolution during the simulation consider the average value of the coherent amplitude on turn n, Y (n) . Definē p p (θ, τ ) = (1/Q y )(dy p /dθ) ,  in direct analogy toȳ p (θ, τ ). For turn n .
(35) Figure 20 shows the evolution of Y for the SNS and PSR parameters in Table 2. The PSR is unstable and losses on the tail of the bunch are apparent in Figure 21.
There is no simulation flag associated with a proton hitting the pipe wall so behavior after the beam gets outside the pipe is unphysical. Figure 22 shows the PSR current pulse and the dipole density, I(τ )ȳ p (τ ). The rms amplitude of 3 mm would produce an easily measureable  Table 2.  Table 2.
instability signal. The same parameters for SNS are plotted in Figure 23. The dipole density corresponds to the Y ∼ 0.5 mm amplitude displayed by SNS throughout the simulation. There is no sign of instability in this case. The PSR parameters in Table 2 correspond to a marginally stable beam in the actual PSR, while the simulation predicts a strong instability. Simulations of the PSR for other intensities, bunch lengths and RF voltages have been done. In general, the onset of the instability roughly corresponds to the coasting beam estimates for parameters in the middle of the bunch. This is not the same as equating rms quantities, since the longitudinal profile of the PSR beam is typically quite peaked. The simulations also show some evidence of nonlinear saturation, which is not surprising given equation (34). However, for such nonlinear saturation to play a fundamental role in the PSR it would be necessary to observe the linear threshold of the instability at intensities well below those required for beam losses. This would be especially true for short bunches and it is not observed in the actual machine.
For SNS, the rms bunch length is a good indicator of peak current. The simulations agree fairly well with the thresholds show in Figure 17.

VI. CONCLUSIONS
Electron cloud instabilities in the PSR and SNS have been explored. Estimates of the SNS cloud density have been made using measured secondary yield data. For 2 × 10 14 protons per bunch we expect less than 5 nC/m of electrons to survive the gap. Similar simulations have been done for the PSR and agree well with experiments. Linear stability theory has been applied to the PSR and SNS. For the PSR, the linear model tends to predict instability for lower currents than are actually observed. This may in part be due to our conservative approximations. For the SNS with 2 × 10 14 protons per bunch, the linear theory predicts that an harmonic-one RF voltage of 15 kV should be adequate to stabilize the beam for an electron density of 5 nC/m. The harmonic one design voltage for SNS is 40 kV. Simulations of the electron cloud instability have been performed. For SNS, the simulations are in fair agreement with the predictions of the linear stability analysis. With 2 × 10 14 protons per bunch and a 60 Hz repetition rate, conditioning rates for electrons with more than 100 eV of kinetic energy are (0.1 → 0.5) C/cm 2 /week for the unconditioned surfaces considered here. A dose of 0.1 C/cm 2 should reduce the peak secondary yield to 1.8 or less. To motivate the linear space charge term in equation (6) consider the coupled Vlasov-Maxwell equations for a beam in a straight channel with uniform, linear focusing. We use (x, y, z) and (p x , p y , p z ) as the phase space coordinates, and clock time t as the evolution variable; with F (x, y, z, p x , p y , p z , t)dxdydzdp x dp y dp z denoting the total charge in the phase space volume dxdydzdp x dp y dp z . The particle velocity is dr/dt = v = p/γm and q is the charge per particle.
The Vlasov equation is given by Multiplying equation (A1) by x ⊥ and using integration by parts yields The key to obtaining (A16) is to notice that the integral of all terms proportional to ∂F/∂p vanish, since the coefficients of these terms do not depend on p. Multiplying equation (A1) by p and proceeding similarly gives −qxk x D x (z, δ, t) − qŷk y D y (z, δ, t) +q (xE x (z, t) +ŷE y (z, t)) Ψ 0 (z, δ) (A17) The first term on the right of (A17) is due to space charge, (A18) The terms on the second and third lines are variants of well known expressions, see e.g. [56]. However, we have been somewhat sloppy with the δ dependence in the terms proportional to k x and k y . To include these terms to leading order, as well as the effects of closed orbit curvature and to change the timelike variable to θ, the reader only needs to check that the head-tail phase shifts are appropriate in the final equations. The point here is to deal with T sc ; the aforementioned subtleties will be neglected.
To continue, notice that the integration with respect to p ⊥ in equation (A18) only affects F so we may consider F (x ⊥ ,z, δ, t) = d 2 p ⊥ F (x ⊥ , p ⊥ ,z, δ, t).
We will now employ first order pertubation theory witĥ F = R 0 (x ⊥ ,z)Ψ 0 (z, δ) +F 1 (x ⊥ ,z, δ, t). Splitting the unperturbed distribution into this product form assumes that the transverse motion does not influence the longitudinal motion. Substituting the pertubation in equation A18 and using A11 where κ = q/2π 0 γ 2 . The zeroth order terms are neglected, but vanish identically when the full equation is considered. The second order terms have been dropped. The next assumption is to take a uniform, elliptical density for the unperturbed distribution.
Strictly speaking a x and a y should vary withz due to the variation in space charge defocussing along the bunch. Substituting this expression into (A19) and noticing that F 1 is zero for x 2 /a 2 x + y 2 /a 2 y > 1 gives T sc (z, δ, t) =