Space charge dynamics in high intensity rings

Uncontrolled beam losses due to space-charge-induced halo generation are a concern in high intensity rings, which are characterized by high beam intensities and low uncontrolled beam loss requirements. It is therefore important to investigate the dynamics of space charge in high intensity rings. We report here the results of extensive calculations using a particle-tracking approach with a self-consistent particle-in-cell model and alternatively with a particle core model. We ﬁnd that the inclusion of space charge forces provides agreement between calculated and experimentally observed beam proﬁle shapes in the high intensity proton storage ring. We also conﬁrm computationally the extension to rings of the accepted dynamics of halo generation with rms beam mismatch exciting the parametric resonance. In addition, we propose a new two-stage mechanism for halo production in rings in which space-charge-driven lattice resonances generate beam mismatch that excites the parametric resonance. Because of its dependence on lattice resonances, this mechanism is peculiar to rings and is capable of generating halo even from initially matched beams. It is also very sensitive to the operating point in tune space, as we show in the results of a vertical tune scan simulating injection into the Spallation Neutron Source accumulator ring. Our results extend and enhance the understanding of fundamental space charge physics, which has been developed for linear accelerators, to rings.


I. INTRODUCTION
Beam dynamics in high intensity rings has become very relevant due to a number of new machines under consideration, including the Spallation Neutron Source (SNS), European Spallation Source (ESS), Japan Hadron Facility (JHF), m-Collider Driver, and others. These machines are characterized by large beam currents and by very stringent uncontrolled beam loss requirements. For example, the accelerator system of the SNS [1] will deliver a 1 GeV pulsed proton beam to a liquid Hg target at 60 Hz. The accumulator ring is being designed to support 2 MW of beam, which implies that it must be capable of holding more than 2 3 10 4 protons in each pulse. In order to expedite hands-on maintenance, the requirement for uncontrolled losses is set to about one part in 10 6 per meter. Because of the necessity of high beam intensity and low uncontrolled beam loss, space charge contributions to beam loss through halo formation are an essential concern for the SNS project.
Space charge effects have been studied for many years, primarily for linear accelerators. Kapchinskij and Vladimirskij included space charge forces to derive an envelope equation and a self-consistent equilibrium solution of the Vlasov equation, called the KV distribution [2]. Since this pioneering work, much progress has been made through the development and generalization of such concepts as rms emittance, the envelope equation, etc. [3]. Numerical simulations, mostly for linacs, and the particle core model indicate that the principal cause of space-charged-induced halo is the 2:1 parametric resonance excited by the envelope oscillations of mismatched beams [3,4].
The extension of this work to detailed numerical calculations of space charge dynamics in rings is the subject of this paper. We report here the results of extensive calculations using a particle-tracking approach with a self-consistent particle-in-cell (PIC) model and alternatively with a particle core model (PCM). Section II describes the physics models employed in our calculations. The discussion includes not only the physics that forms the basis of our PIC model, but also the generalization of the envelope equations to include bending effects in order to provide a PCM for rings. Section III discusses the numerical methods and convergence of our implementation of these models in our particle-tracking computer code, ORBIT [5]. We alternatively utilize a modified version of the particle tracking and injection computer code, ACCSIM [6], in our calculations. ORBIT has been developed in C11 as an object-oriented tracking code, and it contains most of the features in ACCSIM and has been benchmarked to other computer codes, including ACCSIM and SIMPSONS [7]. Section IV presents a detailed comparison of PIC code calculations with experimental results from the proton storage ring (PSR) [8]. Excellent agreement is obtained when space charge forces are included, but not when space charge terms are omitted from the calculations. These results suggest that space charge forces are essential to understanding the formation of experimentally observed beam profile shapes under typical conditions. In Sec. V, we present results, from both PIC and PCM models, that confirm computationally the extension to rings of the accepted dynamics of halo generation, with rms beam mismatch exciting the parametric resonance. Section VI proposes and demonstrates a two-stage mechanism for halo production in rings in which space-charge-driven lattice resonances generate beam mismatch that excites the parametric resonance. Because of its dependence on lattice resonances, this mechanism is peculiar to rings and is capable of generating halo even from initially matched beams. It is also very sensitive to the operating point in tune space. Finally, in Sec. VII we discuss the need for detailed parametric studies with space charge in the design of high intensity rings and we present our conclusions, namely, that the results of the present studies both extend and enhance the understanding of fundamental space charge physics, which has been developed for linear accelerators, to rings.

II. PHYSICS MODELS
High intensity rings are characterized by the separation of longitudinal and transverse scales. In SNS, for example, the longitudinal bunch length is on the order of 100 m, compared with transverse beam dimensions of a few cm; and the longitudinal tune is about 10 23 , compared with transverse tunes of about 5.8. For this reason, it is possible, with good approximation, to separate the longitudinal and transverse dynamics in high intensity rings and, for the study of space charge effects and halo, to concentrate on the transverse dynamics. For particletracking calculations in six dimensional phase space, the rapidly varying transverse canonical coordinates couple into the longitudinal phase equation. However, the only coupling of the slowly varying longitudinal dynamics back into the transverse equations is through the momentum deviation d ϵ Dp͞p 0 and through the modulation of the space charge distribution along the bunch. For this reason, the discussion presented in this section will address space charge in transverse phase space.
Particle tracking is simply the solution of Newton's law for the tracked particles. In transverse phase space this takes the form of a second order differential equation for the particle coordinates q x and y, where prime denotes differentiation with respect to the longitudinal coordinate s and the forces have been separated into K lin q 3 q, the linear magnet force, F sc q , the space charge force, and F nl q , the nonlinear magnet force. The linear focusing force and nonlinear magnet force are external, involving no interaction between the beam particles. They depend on the lattice and are evaluated independently for each particle. The linear terms can be evaluated using first order transport matrices, as generated by either of the computer codes MAD [9] or DIMAD [10]. Although the second order transport matrices from either of these codes could be used to evaluate the nonlinear magnet force terms, the resulting integration would not be symplectic. It is more appropriate to treat these terms symplectically using, for example, thin lens kicks or Lie maps. The space charge force is collective, involving the interaction of all the particle, and it may be written as where the summation is over the other N 0 particles and P is the generalized perveance of the beam [11]. The selfconsistent solution of Eqs. (1) and (2) forms the basis of the PIC calculations that comprise our main approach to the study of space charge in rings. However, the PCM provides an alternative model for the evaluation of the space charge term, which is useful for theoretical studies and for some tracking applications. The analysis and understanding of space charge effects for particle beams in linear accelerators has been greatly facilitated by the use of particle core models [3]. PCMs represent the collective dynamics of the beam by envelope equations that contain the effects of the lattice focusing forces, the beam emittances, and the space charge. PCMs can also be employed to provide space charge forces for particle-tracking calculations to study both collective and individual particle dynamics. One advantage of PCMs is the ability to track specified particles individually or in small groups. PIC calculations require the entire particle distribution in order to evaluate the space charge forces for tracking.
Until recently, PCMs have been applied primarily to straight channels with strong space charge forces. To extend their application to rings, the representation of the beam via the envelope equations must be generalized to include the effects of dispersion. Two similar treatments of this dispersion problem have been carried out [12,13], and we adopt the approach in Ref. [13]. We now present a brief discussion of the derivation and properties of the model, given fully in Refs. [13,14].
The model is derived using an accelerator ordering scheme, based on the separation of transverse and longitudinal dynamics and the small energy spread of the beam, valid for high intensity rings. This ordering is applied to derive a single particle Hamiltonian that decouples the pure betatron motion at each particle energy from the longitudinal motion and from the dispersion. The transverse canonical coordinates and momenta represent pure betatron oscillations about the closed orbit at the precise particle energy. As a by-product of the derivation, an equation for a dispersion function including space charge corrections is obtained. Next, a statistical approach, using the kinetic distribution function in the canonical phase space of the separated Hamiltonian, is applied to perform a moments analysis of the beam at fixed energy values. Effective emittances, conserved to leading order, are defined for the transverse canonical coordinates, and rms envelope equations are derived for the azimuthal variation of the deviations of the second moments. These emittances and envelope equations are decoupled from dispersion in the sense that they pertain to pure betatron motion about the closed orbit at the appropriate energy. The effect of dispersion appears only in determining the spatial charge distribution for the evaluation of the space charge forces. This is found by integrating the charge distribution implied by the fixed-energy envelope equations over the energy distribution of the beam. Thus, the result depends on the energy distribution of the beam. In any case, closure requires the resulting space charge distribution to be consistent with that appearing in the dispersion and envelope equations. To ensure this, an iterative process may be required. We now present a very simple version of this PCM with dispersion for rings, based on the assumption of a generalized KV distribution. By generalized KV distribution we mean a particle distribution with energy spread and having a given KV distribution for its transverse distribution at each energy.
Let q be one of the transverse coordinates and let p be its conjugate momentum. The effective emittance is given by´q where ͗ ͘ denotes an average over the beam distribution function taken at fixed energy. In principle,´q is a function of energy, but we take it to be a constant here. The quantities s p and s q are the rms standard deviations of p and q, respectively, also at fixed energy. The emittances are conserved to leading order in the adopted ordering scheme. The rms envelope equations for x and y can be written as where a 2s x , b 2s y , K lin a 3 a and K lin b 3 b are the linear focusing forces of the lattice magnets, and P is the perveance. These envelope equations pertain to betatron oscillations about the appropriate closed orbit at each energy in the beam. A similar equation is obtained for a generalized dispersion function with space charge, where r is the local bending radius and K lin D 3 D x is the linear focusing force of the lattice magnets. The quantities A and B are analogous to the elliptical axes in the usual PCM, and are found to satisfy where ͗d 2 ͘ tot is the standard deviation of the momentum spread Dp͞p 0 . Thus, A is determined statistically to contain separate contributions from pure betatron oscillations and from the energy spread of the beam. Because Eqs. (4)-(6) are interdependent, they must be solved consistently as a system. For particle tracking, the transverse space charge forces are calculated in this model by assuming a generalized KV distribution with axes A and B, so that the effect of dispersion and energy spreading is contained in Eq. (6) and those of space charge are seen directly in Eqs. (4) and (5) and indirectly through D x in Eq. (6). We note that, for particle tracking, the forces used in this PCM model are linear for particles inside the uniform charge distribution, but drop off nonlinearly for particles outside the distribution; i.e., the finite size of the distribution is considered in evaluating the forces. The PIC and PCM approaches to space charge dynamics are complementary. The object of PIC calculations, embodied in Eqs. (1) and (2), is the simulation of the full space charge dynamics of the beam. For example, the particle distribution can manifest a complete set of moments that, through the space charge potential, can excite a variety of lattice resonances. Also, the particle distribution responds in a self-consistent fashion to its own space charge force, so that the relaxation of the beam due to space charge effects can be observed. Because of the necessity of evaluating Eq. (2), the PIC approach requires the simultaneous tracking of an entire ensemble of particles, and this ensemble must be large enough to provide a good representation of the beam with respect to the dynamics under study. Because it consists of a simple, tractable, closed set of equations, the PCM has been used in theoretical studies of space charge dynamics and has been instrumental in the identification of the parametric resonance as the cause of halo generation in linacs [4]. Computationally, the PCM in rings involves the solution of Eqs. (4)- (6). The range of dynamic phenomena that is accessible using this approach is dictated by the limitations on the beam symmetry implicit in a description in terms of second moments. The phenomena include envelope oscillations, which are simply mismatched or aperiodic (with respect to the lattice periodicity) solutions of the equations, and, when augmented by the coordinates of the beam centroid, the excitation of integer and half-integer lattice resonances. Because the envelope oscillations provide the accepted explanation of space-charge-induced halo generation through the excitation of the parametric resonance, they have been the main topic of theoretical and computational study using the PCM. Although the solutions of the PCM equations do show response to integer and half-integer lattice resonances, they neither excite nor respond to the wide range of resonances seen with the PIC model. Particle tracking can also be carried out using the PCM by calculating the space charge force obtained from the parameters A and B and from the assumed charge distribution. In this case, the space charge force can be regarded as external, and Eq. (1) can be used to track test particles individually without the need for an entire ensemble. This approach is useful for visualizing the systematic effects of the envelope oscillations on the topology of the transverse phase spaces and on x-y coupling. However, in the sense that the tracked particle motions do not couple back into the space charge potential, this approach is not self-consistent. In the work that follows, we use primarily the PIC model to perform self-consistent calculations of beam evolution, but we also occasionally carry out PCM calculations for systematic studies or to observe coupling effects on individual particle orbits.

III. NUMERICAL METHODS AND CONVERGENCE
An essential consideration in particle-tracking calculations is the choice of integration scheme. For Hamiltonian systems, Liouville's theorem and the conservation of phase space require this scheme to be symplectic. For calculations involving the long-term stability of beams in storage rings, symplectic schemes accurate to high order in the integration step size are used. These schemes employ Lie algebraic [15] and truncated power series algebraic [16] methods to derive symplectic nonlinear maps for an entire turn. The advantage of these one-turn maps is that they can be iterated rapidly to track particles accurately over many turns. The use of such maps is restricted, however, to tracking under the influence of external forces. For intense beams, space charge forces cannot be neglected, and numerical stability conditions typically limit integration step sizes to a small fraction of a turn. Although symplecticity remains essential in particle tracking with space charge, the limitation to small step sizes removes the necessity of a high order scheme. Consequently, the integration scheme we have chosen is explicit second order symplectic, with a transport matrix representation of all linear external magnet forces and the inclusion of all nonlinear forces, including space charge, as symplectic thin lens kicks. We include two alternatives for the treatment of space charge: (i) a particle-in-cell model [17] with fast Fourier transforms (FFTs) to evaluate the forces, and (ii) the generalized KV PCM described in Sec. II. The incorporation of nonlinear magnet forces into the model is recent, and we do not utilize this feature in the present results. Consequently, the space charge potential provides the only nonlinearity in the calculations presented here. We now discuss the adopted integration scheme for the PIC model of space charge forces. The only difference when using the PCM is that the envelope equations (4)-(6) are solved along with the particle tracking, and the resulting parameters, rather than the PIC FFT scheme, are used to evaluate the space charge forces.
The numerical model solves the dynamic equations using an explicit, second order, symplectic integration scheme: (i) Divide the lattice into N az elements (all or parts of drifts, bends, quadrupoles, solenoids, and higher multipoles). (ii) Transport macroparticles linearly through these elements using a transport matrix approach, including dispersion. (iii) Treat nonlinear forces as kicks, applied at the ends of linear elements with strengths appropriate to a second order symplectic scheme (force 3 1 2 length of associated element). (iv) External nonlinear forces are applied independently to each macroparticle as symplectic thin lens kicks immediately before and after each linear transport. (v) Space charge forces involve interaction of all beam macroparticles, and require special treatment. They are applied as kicks between the external nonlinear kicks with a strength (force 3 1 2 length of previous plus following elements).
The space charge forces are evaluated as nonlinear kicks using a PIC model and FFTs: (i) At each azimuth s, select a regular rectangular ͑x, y͒ mesh of N FFT 3 N FFT points, centered on the beam. Because of the periodicity of the FFT, the mesh must be at least twice the extent of the beam in each direction. (ii) Obtain the particle density on the mesh by bilinear distribution of the macroparticle charges to adjacent mesh points. (iii) Obtain the FFT of space charge forces at mesh points as the convolution of the FFT of the particle density on the mesh with that of the force due to a unit charge. The inverse FFT of this quantity gives the space charge forces on the mesh points. (iv) Obtain the space charge force on each macroparticle as a bilinear interpolation of the forces at the mesh points to the location of the macroparticle.
The study of space charge using the explicit second order PIC model involves three numerical parameters: N az , the number of azimuthal integration steps; N, the number of tracked particles, or macroparticles; and N FFT , the spatial resolution (grid parameter for the FFT algorithm). The advantage of using the FFT procedure is one of speed, with the number of operations scaling as N for distributing particles to the FFT nodes and for interpolating forces from the nodes to the particles and as N FFT ln͑N FFT ͒ for the transforms. Convergence studies have shown that the accuracy of the solution depends both on having a sufficient number of macroparticles, N, and enough FFT grid points, N FFT 3 N FFT , but that the numerical stability of the integration scheme depends only on having a sufficiently fine azimuthal grid, as determined by N az . The parameters N and N FFT are not completely independent because a sensible charge density representation at the FFT nodes is obtained only when there is a sufficient number of macroparticles per FFT cell. For the calculation of transverse space charge forces this relationship implies a scaling of N~N 2 FFT . In typical calculations, we currently use N 10 5 and N FFT 128, which evaluates to six particles per FFT cell. However, because the extent of the FFT grid is twice that of the beam in each direction, the effective number of particles per cell in the beam is approximately 25 in this case.
In many numerical calculations for collective phenomena in systems of pairwise interacting particles, a smoothing parameter´is introduced to eliminate noise due to strong binary interactions between close macroparticles. Choice of the smoothing parameter´is perhaps one of the more subjective aspects of PIC calculations. Too large a value smooths the charge distribution to the point of erasing the dynamics. It has been argued that too small a value obscures the collective dynamics in noise due to binary interactions between macroparticles. This argument is valid when a direct summation over macroparticle pairs is used to evaluate the space charge force, but the FFT algorithm used here contains its own smoothing effect. This effect is strongest in smoothing out the pairwise forces from nearby particles. Thus, even in the limit of no explicit smoothing, the FFT algorithm smooths the strong interactions with nearby particles. We find that when the number of macroparticles per FFT cell is sufficiently greater than one, no explicit smoothing is necessary. Thus, we typically set the smoothing parameter´to zero. We now study the convergence properties of an example calculation with respect to the parameters N az , N, and N FFT .
We study the effect of the numerical parameters on the convergence of the model for an example case, namely, halo generation by the parametric resonance [4] in a doublet lattice. We consider a doublet lattice with fourfold symmetry similar to the SNS FODO lattice, length 220.668 m, with linear focusing only (48 quadrupoles, 32 sector bends, and 80 drifts), and bare tunes of n x 5.85 and n y 5.70. We consider a coasting beam with a KV distribution, energy and energy spread E 0 1 GeV, DE max 9.4 MeV, DE rms 4.7 MeV, x and y emit-tances´x ,y 100p mm mrad, and number of particles N 3.08 3 10 14 . The initial beam is rms mismatched, resulting in envelope oscillations of about 10% around the matched values. We follow the subsequent evolution for 1200 turns.
The evolution is observed to be as follows: Because the beam is mismatched, envelope oscillations occur, as shown in Fig. 1(a), which plots the beam averages of ͗͑x 2 ͗x͒͘ 2 ͘ and ͗͑y 2 ͗y͒͘ 2 ͘ [͗͑Dx͒ 2 ͘ and ͗͑Dy͒ 2 ͘]. Because the y tune is less than the x tune, the focusing is slightly weaker in y than in x, and through coupling the oscillation energy is transferred from the x direction into the y direction. As the beam relaxes, some particles cross the separatrix of the parametric resonance, driven by the oscillations in ͗͑Dy͒ 2 ͘, becoming halo particles and leading to a growth in the rms y emittance [ Fig. 1(b)]. All emittances presented in this paper will be rms emittances calculated by applying Eq. (3) to the actual beam distribution. This halo-generating process removes energy from the envelope oscillations, which then diminish. As this driving force is removed, the instability saturates leaving the beam in a new steady state with halo. Figure 1(c) plots the final beam cross section in the y-y 0 phase plane.
That the observed evolution is a result of physics, and not a numerical instability, was tested by carrying out the same calculation with the only difference being the method for distributing the charge to the FFT grid. In this latter method, each particle is regarded as a representative of a population of macroparticles, all having the same action values ͑J x , J y ͒ but distributed uniformly over betatron oscillation phases ͑f x , f y ͒. In this "phase averaged" model, envelope oscillations are suppressed so that the parametric resonance is not driven, but no changes are made to numerical methods or parameters. The result of the phase-averaged calculation is emittance conservation and no halo generation, as expected for a stable numerical scheme. This indicates that the observed instability is a result of the dynamics.
We have performed convergence studies in the three numerical parameters, N az , N, and N FFT . The convergence study in the number of azimuthal integration points, N az , resulted in the following: in the presence of space charge, N az must exceed a minimum value to assure numerical stability of the integration scheme, but once the scheme is stable the results are quite insensitive to the number of integration points. In the cases examined here, the tracking with the representation of the lattice by 292 azimuths was numerically unstable. However, both the cases having 488 and 880 azimuths were numerically stable and gave nearly identical results. As a guideline derived from these calculations, we require approximately 20 integration points per lattice cell or 40 integration points per envelope oscillation period.
In the variation of the transverse phase space parameters, the main effect of increasing N and N FFT is a delay of the onset of the instability coupled with a faster rise of the instability following onset. The initial and final configurations are basically the same as are the main features of the evolution. Apparently, the passage to higher resolution decreases the numerical viscosity, thus sharpening the observed evolution.
In order to illustrate the effects of the numerical parameter variations, we now compare the beam evolution for the "base case" (with N 7680 macroparticles, N FFT 3 N FFT 32 3 32 FFT grid cells, and N az 488 steps͞turn to that of the big case (with N 245 760 macroparticles, N FFT 3 N FFT 128 3 128 FFT grid cells, and N az 880 steps͞turn). In both cases the main features of the beam evolution are the same. In the big case the onset of instability is delayed, the growth of the instability is faster, and the beam halo at saturation contains a slightly higher fraction of the beam. Figure 2(a) shows the beam averages of ͗͑Dx͒ 2 ͘ and ͗͑Dy͒ 2 ͘ for the two cases, plotted versus turn number. The most noticeable difference is in the decrease of the vertical fluctuations driving the instability. In the big case this decrease is both delayed and sharper relative to the base case. Figure 2(b) plots the evolution of the rms y emittances for the two cases, and it is again seen that the rise in the emittance for the big case is delayed but sharper in comparison with the base case. The difference in saturation value of the two emittances is very small, with that of the big case being slightly larger. Thus, the base case provides an adequate description of the gross features of the evolution, but higher resolution is required to represent the details of the transient behavior of the instability. Comparison with intermediate cases indicates that, for this process, the big case has converged. To summarize, convergence studies have been carried out for the calculation of halo generation via the parametric resonance, excited by rms mismatch. In these studies, the number of macroparticles per FFT grid cell varied from a minimum of 7.5 to a maximum of 30, although the effective number in the occupied cells is at least 4 times higher since the FFT grid is twice the extent of the beam in both directions. The number of integration points per quadrupole pair is 20 in the case of 488 integration steps and 37 for 880 steps. It was found that the base case of 7680 macroparticles resulted in a reasonable description of the overall dynamics of the beam. Increased resolution revealed a delay in the onset of the instability together with a faster rise time. However, in the case studied here, more than 1% of the beam macroparticles migrated to the halo and the vertical rms emittance growth was more than 3%. In high intensity rings, such as SNS, the limit on uncontrolled losses is 10 24 , or 2 orders of magnitude below the observed losses here. To accurately calculate the dynamics of cases with such small losses, high resolution will be a requirement.

IV. EFFECTS OF SPACE CHARGE: BENCHMARK WITH PSR EXPERIMENT
As a benchmark of the ORBIT code, we now present comparisons of simulated and experimentally measured transverse beam profiles at the end of injection in the Los Alamos Neutron Science Center PSR. This work is described in detail in Ref. [18]. Experimental beam profiles were generated using H 2 foil injection from an 800 MeV linac. Injection was carried out for a duration of 825 ms, or 2305 turns, for two painting schemes and three beam intensities. The horizontal and vertical bare tunes for the ring were set to n x 3.17 and n y 2.14, respectively. Subsequently, the beam was extracted in a single turn and transported to a wire scanner beam profile diagnostic in the extraction beam line.
In the benchmark calculations, care was taken to match the full PSR lattice parameters, injection scenarios, and beam intensities [18]. Because the absolute transverse location of the wire scanner was not precisely known, the centers of the calculated and experimental distributions were matched to facilitate the comparison. The required shift proved to be independent of beam intensity. Also, because the calculated and measured distributions were obtained in different arbitrary units, the vertical scales of the distributions were adjusted so that the calculated height matched that from the wire scanner data at the center of the distribution. We emphasize that the width and shape of the calculated distribution were not adjusted.
Comparisons were made for two injection scenarios, both with and without a vertical closed orbit bump at the foil. The PSR does not have bump magnets to paint in the horizontal direction, but a horizontal spread is provided through the injected beam energy distribution since there is a dispersion of 1.4 m at the foil location.
The overall results reveal good agreement between the measured and calculated transverse profiles. For both injection scenarios, there is a noticeable broadening of the vertical distribution at high beam intensity in both the experimental measurements and calculated profiles. The horizontal profiles are much less sensitive to beam intensity. Also, the horizontal profiles are less affected by space charge than are the vertical profiles. A likely explanation of the vertical broadening at high beam intensity is that the space charge tune shift is sufficient to cause the beam to encounter the resonance at n y 2. The calculations reveal strong coherent second moment (envelope) oscillations in the vertical direction, and these oscillations occur with a tune close to 4.0. The injected beam has a rectangular shape, indicating strong quadrupole and octupole moments. Figure 3 shows the horizontal and vertical quadrupole moments of the beam normalized to remove the lattice effects, ͗͑x 2x͒ 2 ͘͞b x and ͗͑y 2ȳ͒ 2 ͘͞b y , calculated over the last three injection turns for the high intensity case with 3 3 10 13 injected particles and no vertical closed orbit bump. A well-defined coherent vertical envelope oscillation of the beam at a frequency of 4 oscillations͞turn is evident. The horizontal quadrupole moment behavior is much less regular, with variations resulting from the interplay of dispersion and space charge. A plot of the octupole moments of the beam shows very similar behavior. The same calculations without space charge do not develop these moment oscillations. It is important that the number of vertical oscillations shown in Fig. 3 is very close to 4 per turn. This resonance between the coherent vertical second moment  oscillations and the lattice will enhance the parametric resonance and contribute to beam broadening. The octupole moments of the beam may also contribute to the redistribution of space charge. One possible mechanism for this is through the resonance 2 3 ͑n x 2 n y ͒ 2. Because there are a number of other beam moments with coherent oscillations, it is possible that the redistribution of the beam by space charge is carried out through a combination of mechanisms. Further studies are underway to elucidate the details of the dynamics.
Our results suggest that the good agreement between experiment and calculations is very dependent on including the space charge force in the calculations. This is demonstrated in Fig. 4, which shows the experimental and calculated horizontal and vertical beam profiles at full intensity, calculated both with and without the space charge force. Although space charge does not strongly affect the horizontal beam profile, the vertical profiles with and without space charge are very different. For this case the injection is offset vertically with no bumps. This is reflected by the "horns" on the vertical profile calculated without the space charge force, which occur at the location of injection and have widths corresponding to the linac beam size. In both the experiment and the calculation with space charge, space charge forces spread these horns, causing the distribution to fill in and to broaden. The transverse space charge model predicts most of the observed beam broadening. Although it is possible that other nonlinear effects that are neglected in the calculations, such as magnet errors, could be responsible for the observed beam profiles, the dependence of the spreading on beam intensity, together with the observed moment oscillations at high intensity, is more consistent with space charge for its explanation.
It is interesting that for this case the space charge tune shifts are approximately 0.15, giving a 5%-8% tune depression in the vertical direction. In linacs this value of tune depression would lead to the characterization of the beam as emittance dominated, and space charge forces would be deemed unimportant. In rings, however, the beam effectively traverses a much longer path than in linacs, and it interacts with the periodicity of the lattice. For these reasons, space charge forces that give tune depressions of a few percent, which are modest in the linac regime, can be essential to describe the correct physics in rings. The present benchmark comparison suggests that the space charge calculations are both correct and important, particularly considering that the intensity of the PSR beam at high currents is comparable to that in proposed spallation neutron sources.

V. HALO GENERATION IN RINGS VIA THE PARAMETRIC RESONANCE
In a ring, the beam is rms matched when the rms beam radii vary periodically with the ring periodicity. The rms parameters of matched beams are closed solutions of the rms envelope equations having the periodicity of the ring. rms mismatch manifests itself as envelope oscillations at frequencies slightly less than twice the bare tune in the corresponding plane [4]. These oscillations drive the parametric resonance, which is characterized by approximately two envelope oscillations for each particle betatron oscillation. Because of this relationship between the betatron tunes and the envelope oscillation frequencies, the resonance typically occurs as a phase space separatrix in an annular region outside the beam core. For rms mismatched beams, particles that encounter this separatrix track around its edge, and some can be diffused or bumped across the unstable fixed points to become halo particles.
The dynamics of halo formation due to the parametric resonance for a mismatched beam in a doublet ring was described briefly in Fig. 1 and the related discussion. Returning to that case, we show in Fig. 5(a) the fraction of the beam in the halo versus the number of turns, where a particle is defined to be in the halo when its emittance exceeds the initial maximum emittance by at least 50% in either plane. By comparison of Figs. 5 and 1, the growth in the number of halo particles is seen to closely track the emittance growth in the vertical plane, and more than 1% of the beam particles ultimately enter the halo. A check of the halo particles for this case shows that all have vertical emittances and none have horizontal emittances exceeding 150p mm mrad, again demonstrating that the instability lies in the vertical plane. Further confirmation of this is given in Figs. 5(b) and 5(c), which plot the vertical versus horizontal actions and tunes, respectively, of the final beam distribution. The tunes were calculated by tracking the phase advances of each particle over a single turn. In action space the initial generalized KV distribution is a uniformly populated line from 0.5 on the vertical intercept to 0.5 on the horizontal intercept. The interaction of the particles as the beam evolves causes this initial line to broaden into a band. However, Fig. 5(b) shows a sizable population of halo particles having large emittances near the vertical axis. The halo can be seen in tune space as a cloud of particles with vertical tunes above those of the main beam core, Fig. 5(c). Because they spend more time at greater distance from the beam, halo particles have large actions and are less tune shifted by space charge than are the core particles. Finally, Fig. 5(d) is a vertical phase space plot of the beam distribution at 600 turns, just at the outset of halo formation. This figure clearly shows the topology of the parametric resonance, including the particles crossing the unstable fixed points of the separatrix to enter the outside halo region.
In addition to simulating the basic halo formation process via the parametric resonance, it is also possible to   6. (Color) Vertical emittance evolution for four beam intensities. Increased beam intensity yields increased emittance growth due to halo formation. use particle tracking to examine the systematic behavior as various physical parameters are changed. We now present the results of a series of parameter scans using particle tracking that describe more fully the systematic behavior of the halo formation process. The parametric resonance is driven by the space charge force and by rms beam mismatch, so the first variations we consider are of those quantities. For these studies we have used the fourfold symmetric FODO lattice for the SNS accumulator ring, with bare tunes of n x 5.85 and n y 5.70. We consider a coasting beam with a generalized KV distribution, energy and energy spread E 0 1 GeV, DE max 9.4 MeV, DE rms 4.7 MeV, x and y emittances´x ,y 100p mm mrad. Although we use a coasting beam and uniformly populated longitudinal phase space distribution, the stated number of beam particles will be normalized to 65% of the ring length of 220.668 m. The calculations were carried out using 20 000 macroparticles. Figure 6 shows the vertical emittance versus turn number for four cases with 10% radial rms beam mismatch, identical except for the number of beam particles assumed in calculating the space charge force. It is seen that as the strength of the space charge force, i.e., the beam intensity, increases, so does the saturated value of the emittance, while the rise time decreases. Thus, increasing the beam intensity leads to increased beam halo for comparably mismatched beams. Note that with no space charge force the emittance is constant.
The effect of beam intensity on individual particles can be conveniently studied using the PCM. In these calculations, the rms envelope equations are solved according to the prescription in Eqs. (3)-(6) using the desired space charge beam intensity, and the resulting envelope parameters are used to calculate the force on individual particles under the assumption of a generalized KV distribution. We again note that the forces are evaluated nonlinearly for particles outside the distribution. Figure 7 presents Poincaré plots of individual particle trajectories, taken once per envelope oscillation period, in vertical phase space with different colors corresponding to different particles. Because the period of the envelope oscillations is incommensurate with that of the lattice, it was necessary to divide y and to multiply y 0 by the square root of the vertical beta function, p b y . The particles tracked in Fig. 7 were initialized with zero action in the horizontal plane and, because of the symmetry of the lattice and the generalized KV distribution, the transverse motion remains vertical. The plot in Fig. 7(a) is taken for a beam intensity corresponding to 1 3 10 14 particles, and the plot in Fig. 7(b) is for a beam intensity corresponding to 4 3 10 14 particles, both with radial rms mismatch of 10%. It is seen that increasing the beam intensity does not substantially increase the width of  the parametric resonance island structure. Rather, the island shifts inward toward the beam core with increasing intensity, thus becoming more accessible to beam core particles and allowing a greater fraction of the beam to cross the separatrix and become halo.
To examine the effects of rms mismatch on halo formation, we now fix the number of beam particles to 2 3 10 14 and vary the radial mismatch. Figure 8, generated using the PIC model, plots the evolution of the rms emittance in the vertical plane for a rms matched beam and for 10%, 20%, and 30% radially mismatched beams. There is a significant increase in both the saturation value and the growth rate of the emittance with increasing mismatch. Comparison with Fig. 6 shows a stronger dependence of halo formation on mismatch than on beam intensity for this case. Figure 9 shows the effects of increased rms mismatch on the parametric resonance obtained by tracking individual particles in vertical phase space using the PCM. The plot in Fig. 9(a) was made with 10% radial mismatch, and the plot in Fig. 9(b) was made with 40% radial mismatch, both with 1 3 10 14 particles. The resonance island structure becomes significantly larger and moves away from the core as the mismatch increases. Although the island is farther from the beam core, the increased excitation due to the greater mismatch couples with the larger island size to give the large observed emittance increases.
Unlike linacs, rings display the effects of dispersion because of the presence of bending magnets. We examine the coupling of these effects with the particle energy distribution using the PIC model by varying the particle energy distribution for otherwise identical cases. For example, consider again the case of the mismatched doublet lattice shown in Figs. 1 and 5. In the initial case, the energy distribution satisfied E 0 1 GeV, DE max 9.4 MeV, and DE rms 4.7 MeV. In that case, the vertical emittance saturated at about a 3% increase. We now consider two cases identical except for the energy distribution. One case has zero energy spread, with DE max DE rms 0, and the other case has a large energy spread of DE max 18.8 MeV, and DE rms 9.4 MeV. Figure 10 shows the evolution of the horizontal and vertical emittances for these two cases. The result is that the emittance growth is much stronger for the case with zero energy spread than for the initial case, which is in turn stronger than that in the case with large energy spread. We have found this result consistently over many cases: for a given lattice with bending, increasing the energy spread tends to stabilize the parametric resonance. One mechanism that contributes to this is described by Eq. (6): dispersion broadens the beam with increasing energy spread and this broadening decreases the space charge force felt by the beam particles. Also, this broadening is determined by the lattice periodicity and is thus incommensurate with the frequency of the parametric resonance. We note that these results are for varying energy spread and fixed lattice configuration; we have In the parametric resonance studies presented thus far, the bare tunes have been n x 5.85 and n y 5.70, and the emittance growth has been in the vertical plane. It has been stated that halo formation prefers the direction of lesser tune. In order to demonstrate this, Fig. 11 presents the results of the dispersionless calculation of Fig. 10 together with those from an identical calculation, except that the tunes have been flipped to n x 5.70 and n y 5.85. Although the behavior is less dramatic when the tunes are reversed, it is seen in this case that the horizontal emittance, rather than the vertical emittance, grows. In other cases we have observed this same tendency for the halo to form primarily in the plane corresponding to the lower tune value. We interpret this as a result of weaker focusing forces in the lower tune direction.
The PCM can be used with particle tracking to understand some aspects the complicated nature of single particle motion in the parametric resonance. One example is the dependence of the motion on x-y coupling. Figure 12 shows, for the SNS FODO lattice case with 2 3 10 14 particles and 10% rms radial mismatch, Poincaré plots of horizontal phase space. The particles tracked in Fig. 12(a) have pure horizontal motion lying on the vertical midplane, while those in Fig. 12( b) were initialized with equal horizontal and vertical actions. The outer surfaces display some distortion due to the coupling, but the main effect is that the parametric resonance island moves inward radially and displays a different detailed structure for different particles. Thus, from the single particle standpoint, the details of the structure of the parametric resonance depend on the coupled x-y motions of the individual particles. Figure 13 further elaborates the effects of coupling by plotting, for a single particle with initially equal horizontal and vertical actions, the horizontal and vertical phase space Poincaré plots and the horizontal and vertical coordinate values versus turn number. The coupling effect is clear in the beating of the coordinate amplitudes as the motion progresses. This beating explains the complicated appearance of the island structure, which changes with the relative amplitudes and phases of the horizontal and vertical motion.
Before leaving the subject of the parametric resonance, we make a number of comments in summary. The parametric resonance provides a mechanism for halo generation in rings. The parametric resonance is excited by the space charge force and by envelope oscillations caused by rms beam mismatch. In the presence of dispersion, energy spread tends to damp the parametric resonance and retard halo formation. Coupling effects complicate attempts to visualize the dynamics of the parametric resonance in terms of single particle motion. The results shown here indicate, for conditions appropriate to high intensity rings, that space charge tune shifts are on the order of 10 21 , or a few percent of the bare tunes, so that the beam is emittance dominated. Thus, rings present a different regime for space charge dynamics than linacs, in which space charge tune shifts are a large fraction of the bare tunes. In linacs, halo formation occurs within a few envelope oscillations. However, according to the cases presented here, halo formation in rings requires hundreds of turns, or thousands of envelope oscillations. Although this process is slow by linac standards, in a ring the beam undergoes many turns, providing substantially more time for the parametric resonance to act than in linear accelerators. The point is that there is time for halo formation to occur. In the next section we consider the interaction between space charge, lattice resonances, and the parametric resonance in generating halo.

VI. EXCITATION OF BETATRON RESONANCES BY SPACE CHARGE AND A TWO-STAGE MECHANISM FOR HALO GENERATION IN RINGS
In addition to exciting the parametric resonance, space charge can also excite betatron resonances [19]. We now present two examples of the excitation of resonances by space charge, one using PIC code results and the other using the PCM. Figure 14(a) plots the fraction of the beam in the halo as a function of time for a doublet lattice with 10% rms mismatch and 2 3 10 14 particles. The three cases differ only in the value of the vertical tune, which is n y 5.70 (red), 5.79 (blue), and 5.83 (green) for a horizontal tune n x 5.85. The figure shows that as the vertical tune moves closer to the horizontal tune, the generation of halo increases. This is a manifestation of the Montague resonance, which occurs at 2 3 n x 2 3 n y and which will be discussed in more detail below. Figure 14( b) was obtained by solving the PCM with n x 6.035, n y 5.885, and 2 3 10 14 particles and plotting once each turn the ratio of the calculated horizontal and vertical envelope amplitudes to those of the matched beam solution. Under these conditions the horizontal envelope amplitude frequency is very close to 12, which is a structure resonance at 3 times the lattice periodicity. Even though the envelope amplitudes were initialized to those of a numerically matched single turn solution, the horizontal amplitude develops large oscillations within a few turns, and the vertical amplitude, driven by the horizontal oscillations, also develops oscillations of about 10% around the matched solution. It has been pointed out by Baartman [20] that space-charge-driven betatron resonances in rings are collective in nature and occur at the coherent, rather than bare tunes. This is illustrated in Fig. 14(b) where the bare tune is 6.035, but the x envelope oscillation has a tune very close to 12. It is clear that the excitation of betatron resonances by space charge is very sensitive to the tunes, which is not the case for the excitation of the parametric resonance through beam mismatch. However, when betatron resonances are excited, by space charge or otherwise, one of the consequences is a resulting mismatch of the beam. This mismatch can drive the parametric resonance, leading to halo formation, even for initially matched beams. Thus, although it is necessary to have well-matched beams to avoid halo formation, this alone may not be sufficient. We now show that, depending on the operating point in tune space, space charge forces can drive resonances that lead to beam mismatch, even for initially rms-matched beams, and that this mismatch can lead to halo formation through the parametric resonance. This mechanism for halo formation was first pointed out by Jeon et al. [21].
Consider three cases in the fourfold symmetric SNS ring lattice having initially matched beams with generalized KV distributions. For these cases, the initial transverse rms beam emittances are chosen to be 120p mm mrad, the proton beam energy is 1 GeV, the beam intensity is set so that the initial tune shift is about 0.08 giving a tune depression of 1.4%, and the initial horizontal bare tune is n x 5.82. The three cases are identical except for varying the choice of vertical bare tune n y , which is taken to be n y 5.67, 5.82, and 5.77.
Before proceeding, we note that an ideal KV distribution supports only even terms in the space charge potential. The numerical generation of such distributions for PIC codes is approximate, involving the random statistical population of phase space by macroparticles according to the mathematical distribution function. Accordingly, all terms will be present if there are asymmetric nonuniformities in the actual numerically generated beam distribution. However, we expect the even mode terms, which are consistent with the beam symmetry, to be dominant. If the operating point in tune space is sufficiently close to a resonance that can be excited by an even mode potential, it is possible for the space charge to excite this resonance and destabilize the beam. As shown below, this will apply even for the case of an rms matched beam.
Considering first the case n y 5.67, the operating point is such that no betatron resonances are excited. In a calculation tracking the rms matched beam for 1250 turns, the transverse rms emittances are observed to remain constant, envelope fluctuations remain small, and the beam phase space plots show no halo. To observe halo generation for this operating point, it is necessary to mismatch the initial envelope radii by more than about 5%, beyond which the extent and population of the halo increases with increasing mismatch.
We next consider the case of n x n y 5.82. This point in tune space lies atop the Montague difference resonance 2n x 2 2n y 0. Difference resonances are not usually regarded as dangerous to the operation of accel-erators because of the conservation of J x 1 J y , where J x͑ y͒ is the x͑y͒ action. However, it has been shown [22] using self-consistent Vlasov-Poisson equations that the space charge potential supports an even mode contributioñ a 0 x 4 1 a 2 x 2 y 2 1 a 4 y 4 capable of exciting the Montague resonance. Given such a driving term, the Montague resonance facilitates halo generation in two ways: (i) by increasing the strength of the 2:1 parametric resonance through the coupling of the horizontal and vertical envelope oscillations; and (ii) by helping particles, through coupling, to cross the separatrices of the 2:1 parametric resonance. Figure 15 shows the evolution of the transverse rms emittances and envelope fluctuations for the case of a well-matched beam at n x n y 5.82. The small initial oscillations in the second moments reflect the quality of the match. Even though a 0 x 4 1 a 2 x 2 y 2 1 a 4 y 4 is not present in the ideal KV distribution, and is initially very small in the numerically generated generalized KV distribution, this term can grow to become significant and excite the Montague resonance. Evidence of the resulting x-y coupling is clear in the anticorrelation of the FIG. 15. (Color) Evolution of (a) rms emittances, (b) second moments ͗Dx 2 ͘ and ͗Dy 2 ͘, and (c) fourth moments ͗Dx 2 Dy 2 ͘ for n x n y 5.82. Although not plotted here, ͗Dx 4 ͘ and ͗Dy 4 ͘ are similar to ͗Dx 2 Dy 2 ͘ in appearance. horizontal and vertical emittances in Fig. 15(a). This leads to rms mismatch, as shown by the increasing oscillation amplitudes of the moments shown in Figs. 15(b) and 15(c), and ultimately to the strong excitation of the 2:1 parametric resonance. Little halo formation is observed for this case for over 2000 turns, but ultimately substantial halo formation accompanies a large emittance increase at around 2300 turns. The behavior of the fourth moments, which reflect the coupling term in the space charge potential, is shown in Fig. 15(c). This behavior is similar to that of the second moments, and it supports the involvement of the Montague resonance in the beam destabilization. The strength of this process is due to the enhancement of the parametric resonance by the Montague resonance, which couples the horizontal and vertical envelope oscillations and aids the particles in crossing the separatrix through x-y coupling.
Let us now consider the case of n x 5.82 and n y 5.77. In this case, the fourth order resonance 4n y 23 is excited by the fourth order even mode term, a 0 x 4 1 a 2 x 2 y 2 1 a 4 y 4 , in the space charge potential. Even for the initially well-matched beam, this induces a mismatch that leads to significant halo formation. Unlike odd mode terms in the space charge potential, which are skew terms in the Hamiltonian that are usually small in lattices with midplane symmetry, even mode terms can be significant. Consequently, when the space charge potential is the driving force in such lattices, the 4n y n resonance associated with the fourth order even mode potential term is significantly more excited than the 3n y n resonance associated with the third order odd mode term, a 1 x 2 y 1 a 3 y 3 [22].
When the coherent tune is near a betatron resonance [20] associated with an even mode space charge potential, the excited resonance can generate a mismatch, even though the initial beam is well matched. This mismatch can drive the 2:1 parametric resonance and expedite halo formation. Figure 16 illustrates this process for the case with n x 5.82 and n y 5.77. The moderate y rms emittance increase around 200 turns in Fig. 16(a) is due to the excitation of the fourth order 4n y 23 resonance. This can be seen in Fig. 16(b), which is a snapshot of the particle distribution in y phase space taken at 250 turns. An rms beam mismatch is generated by this resonance as can be seen in the plot of the beam moments in Fig. 16(c). This mismatch leads to the excitation of the 2:1 parametric resonance and to the subsequent large increase in the y emittance. The anticorrelation of the x and y emittances indicates that coupling is also playing a role in the destabilization. The y phase space plot in Fig. 16(d) superposes the beam distributions at 1000 turns and 1250 turns, and shows the 1:2 parametric resonance structure and halo. Clearly, it is advantageous to avoid nonlinear resonances associated with even mode space charge potentials when choosing an operating point for intense beam circular accelerators. This is especially true when the horizontal and vertical tunes are close.
The main conclusion of this section is that, depending on the operating point in tune space, space charge forces can excite resonances that lead to beam mismatch. This mismatch can then excite the parametric resonance, thus leading to substantial halo formation, even for initially matched beams. Furthermore, if the x and y tunes are close together, the resulting coupling can contribute to the destabilization of the parametric resonance and to the resulting halo formation.

VII. FUTURE WORK AND CONCLUSIONS
We have carried out a number of injection studies for the SNS ring, including optimization of the injection scheme for minimum halo formation [23]. The most recent injection studies include a survey of vertical tune space for fixed horizontal tune, n x 5.82, and a comparison of doublet and FODO lattices having the same global parameters, including length, beam size, and bare tunes. These studies track the entire injection process, including such features as incident linac beam distribution, foil stripping model, and bumping schemes, for 1158 turns. More details of the injection model are presented in Ref. [23]. Figure 17 shows the results of the vertical tune survey by plotting the percentage of particles with x or y emittance exceeding 180p mm mrad as a function of bare tune at the end of injection. The resolution of this tune scan is coarse, with cases separated by Dn y 0.05 in the range from 4 , n y # 6, but it is clear that halo formation is sensitive to the operating point in tune space. It is interesting that, although bare tunes of exactly five or six give low halo formation, the worst operating points flank these integral values. The divergence of the results as n y approaches four is due to the structure resonance with the fourfold symmetric SNS lattice. We also note that there are a number of operating points with no particles exceeding 180p mm mrad, indicating that it may be possible to avoid producing halo through proper choice of operating point.
Calculations have also been carried out comparing halo formation for doublet and FODO lattices having the same global parameters. For both lattices, the results are sensitive to the tune values, but in some cases the doublet lattice produces less halo and in other cases the FODO lattice produces less.
Practical design studies with space charge using realistic injection scenarios and beam distributions to explore the operating space of high intensity rings have been and are being carried out, but these studies are of an exploratory, rather than comprehensive nature. It remains to explore operating space in much greater detail than that shown in Fig. 17, and to understand the reasons for both the good and bad operating points with respect to halo formation. Further, it remains to perform calculations including higher order (nonlinear) magnet forces and errors and to understand their effect on the physics of space charge, both for idealized beam distributions and for realistic injection scenarios. In spite of the work yet to be done, we have presented here a number of detailed studies and new results regarding the dynamics of space charge in rings.
We have described two space charge force models, a particle core model for rings that includes the effect of dispersion and a PIC model, that are both incorporated into a particle tracking code. We have described the convergence of these models in terms of the numerical parameters. We have demonstrated by benchmark comparison with experiment that the inclusion of space charge forces is an important component in any beamdynamic model of high intensity rings. We have studied in detail halo generation for mismatched beams by the parametric resonance in rings. We have also shown that nonlinear space charge forces can excite various lattice resonances that lead to rms beam mismatch. This can trigger a two-stage process in which the mismatched beam in resonance with the lattice excites the parametric resonance, resulting in substantial halo formation. These results, along with realistic injection studies surveying the tune dependence of halo formation, emphasize the sensitivity of beam stability to the operating point in tune space. Comparison of halo formation in FODO and doublet lattices with the same global parameters shows that results are sensitive to the tune, but reveals no performance bias for one or the other. Future work will extend these parametric space charge studies and consider space charge dynamics in the presence of nonlinear magnetic forces, such as errors, fringe fields, and higher multipoles, which have recently been incorporated into the ORBIT code.