Space charge and octupole driven resonance trapping observed at the CERN Proton Synchrotron

124201-1 The combined effect of space charge and nonlinear resonance on beam loss and emittance was measured in a benchmarking experiment over a 1.2 s long flat bottom at 1.4 GeV kinetic energy in the presence of a single controllable octupole. By lowering the working point towards the resonance, a gradual transition from a loss-free core emittance blowup to a regime dominated by continuous loss was found. We compare the observation with 3D simulations based on a new analytical space charge model and obtain good agreement in the emittance blowup regime. Our explanation is in terms of the synchrotron oscillation, which causes a periodic tune modulation due to space charge, and leads to trapping and detrapping on the resonance islands. For working points very close to the resonance this induces a beam halo with large radius. The underlying dynamics is studied in detail, and it is claimed that the predicted halo in conjunction with a reduced dynamic aperture for the real machine lattice is the source of the loss observed in the experiment.


I. INTRODUCTION
A detailed theoretical understanding as well as confidence in simulation modeling of the long-term (10 5 -10 6 turns) effects of high intensity or high phase space density is crucial for the SIS100 of the GSI future project [1,2], where it is necessary to hold the high-intensity bunches between injections over typically 1 s at a loss level not exceeding 1%, likewise for the optimum performance of the CERN Proton Synchrotron (PS) for high-intensity beams.A major focus of such studies is the combined effect of space charge and nonlinear resonances and its impact on halo formation and/or beam loss.
Up to now comparison of simulation with experimental work for second or higher order resonances dominated by space charge has been successfully carried out in the millisecond time frame [3,4].In the realm of long-term behavior, instead, where self-consistent 3D simulation is beyond current computer capabilities, the question of adequate approximations in space charge calculation is a challenging issue.Moreover, an explanation of the proper mechanisms describing the combined effect of nonlinearity, space charge, and synchrotron oscillation is necessary to gain sufficient confidence in the simulation modeling.In the simplified simulation model of Ref. [5] it was recently suggested that the single-particle tune modulation caused by space charge and synchrotron motion may lead to trapping and detrapping on the resonance islands, which are moving in and out.This process is related to the trapping during a single passage through a higher order nonlinear resonance as a result of a changing tune (see Refs. [6,7]).Experimental studies of the effect of machine tune modulation on island trapping or detrap-ping, in the absence of space charge, have been carried out at the CERN Super Proton Synchrotron [8,9], where the modulation was found to lead to diffusion effects, at the Fermilab Tevatron [10], and at the Indiana University Cyclotron Facility cooler ring [11].The latter work gave first evidence of the single-particle motion near islands by taking advantage of the small emittances obtained by electron cooling.The main point of the present work is to study these phenomena in the context of the intrinsic tune modulation caused by space charge and synchrotron motion, which has not been addressed so far.This introduces additional complexity due to the following features: (i) the tune of an individual particle is strongly modulated, by an amount comparable with the incoherent tune shift, depending on two parameters, the synchrotron phase as well as the betatron amplitude; (ii) the latter is itself evolving in time, depending on the preceding trapping and detrapping events; and (iii) there is a selfconsistent time evolution due to the global changes of space charge.These circumstances do not allow experimental verification of single-particle behavior.We use 3D computer simulation with space charge to demonstrate the long-term balance of trapping and detrapping as competing mechanisms and to compare the global predictions on the rms emittance and halo growth with the measurements.

II. MEASUREMENTS
The measurements were carried out as part of a highintensity machine development time at the PS in October 2002.The number of protons in the single bunch (200 ns long at 4) was 1:1 10 12 (parameters, see A vertical maximum space charge tune shift of 0.12 and a horizontal one of 0.075 (for minimum amplitude particles) were achieved with relatively small emittances of x 9 mm mrad and y 4:5 mm mrad (unnormalized at 2).The space charge tune shift was chosen significantly below the maximum possible in the PS to avoid overlap with other resonances.The bunch profiles measured 10 ms after injections were found to be Gaussian in all directions in the absence of the octupole.The vertical machine tune was set to Q y 6:12, and the horizontal one was varied in the interval 6:25 < Q x < 6:32.The chromaticity is close to the natural one, hence the small momentum spread of 3 10 ÿ3 (at 2) allows one to ignore the chromatic effects.The kinetic energy was kept at the injection value of 1.4 GeV with a measurement window of 1.2 s (4:4 10 5 turns) over which the bunch intensity was monitored with a current transformer.The calibrated octupole (here K 3 1:215 I m ÿ3 ) was powered to 40 A at 110 ms after injection to excite the resonance 4Q x 25.We used the transverse profiles measured with the flying wire (20 m=s), fitted them with a Gaussian profile, and determined the corresponding rms emittances.Initial and, in most cases, final profiles were actually found quite close to Gaussian.
In Fig. 1 results of final measurements 1.2 s after injection are plotted as a function of the machine working point.
Our main finding is the existence of two regimes: an emittance growth dominated regime for Q x sufficiently above the resonance (in our example Q x > 6:28) and a loss dominated regime for Q x < 6:28.It is noted that for the working point of maximum loss (Q x 6:27) the emittance also shrinks, since large amplitude particles are predominantly lost.The time evolution of the bunch intensity for Q x 6:27 is shown in Fig. 2. Note the continuous loss at a nearly constant rate after an initially enhanced loss (the intensity drop at 1200 ms is caused by the kicker event).

III. SIMULATION
Interpretation of these measurements must rely on adequate computer simulation.The comparison is also a necessary basis for code benchmarking on such a long time scale, which has not been undertaken so far to our knowledge.To initiate such a process we have carried out a series of simulation runs in 2D and 3D.

A. The model
The main factors to explain the observed behavior are space charge, the octupole, and synchrotron motion.While the octupole is the driving mechanism for the resonance, the detuning effect, as an amplitude limiting mechanism, is not only due to the octupole and other lattice nonlinearities, but mainly due to space charge.Here we are dealing with the transverse amplitude dependence of space charge as well as the strong periodic modulation of the transverse tune shift during a synchrotron period, which is the main source of tune modulation in our experiment.Chromaticity dependent tune modulation is more than an order of magnitude lower and ignored in our simulations.
The variation of space charge between beam core and halo is best seen on a plot of single-particle tunes in a tune diagram for the initial ensemble considered in Sec.III D as shown in Fig. 3.
The ''tune footprint'' reflects particles in the bunch center with small synchrotron and betatron amplitudes (left lower tip) as well as particles at the bunch ends (right upper tip).In the example shown (machine tune at Q x 6:26) the resonance line intercepts with such particles near the bunch ends.Transverse halo particles are absent in the initial footprint, since the initial transverse distribution has only particles inside a sharp edge at 3. We are referring to this region as the ''beam core,'' whereas particles beyond 3 (where they get to as a result of the resonance) are called ''halo.''Note that the larger horizontal emittance causes an asymmetry in the footprint.
We have replaced, for simplicity, the linear PS focusing lattice by constant focusing and ignored lattice nonlinearities besides the contribution from the localized octupole.We have also ignored the smaller vertical beam  emittance of the experiment and assumed a circular cross section for space charge calculation in the analytical 3D space charge model, which is based on a rotational ellipsoid.The horizontal emittance has been redefined accordingly such as to reproduce accurately the maximum horizontal space charge tune shift extracted from the measurement, which we believe is the crucial issue here since we are not dealing with a coupling resonance.

B. Dynamic aperture
The loss observed in the experiment must be related to the shrinking of the dynamic aperture, since the beam was too small to hit a physical aperture.To roughly explore the effect on the dynamic aperture of the octu-pole alone we have carried out a numerical test by searching the maximum stable radius of test particles placed into 20 different directions in the upper half of the x-y plane and ignoring space charge at this point.We have found that the nominal octupole (40 A) leads to a dynamic aperture (10 5 turns) of about 5 near Q x 6:25, where is the horizontal rms beam size of the injected beam.This value is not small enough to explain the observed loss of particles as will be seen in the subsequent simulations.Hence, a more complete knowledge of machine nonlinearities at the working points used here may be needed to explain the observed loss, which has to be the subject of future measurements.Assuming 200 A octupole current, we have calculated that the dynamic aperture shrinks to a radius of 2:5 near Q x 6:25 for 10 3 turns and about 2:2 for 10 5 turns.

C. 2D simulation
We first attempt a comparison with the fully self-consistent 2D particle-in-cell (PIC) version of the MICROMAP code [12] with 10 5 simulation particles, mainly to demonstrate that synchrotron motion and 3D are indispensable.We employ a Gaussian distribution function and a 64 64 grid filling a rectangular boundary of 70 70 mm size.We find no loss for 40 A: the rms emittance growth remains below 2%.In the absence of space charge it is 15% at Q x 6:25, independent of the octupole strength, which only determines the time it takes to reach saturation (about inversely proportional to the octupole strength).This low saturation level of 2% is a result of the well-known mechanism of nonlinear detuning [6].The effect is enhanced here by the large space charge detuning relative to the natural detuning effect of the octupole.It is crucial, for this detuning to be effective, that the single-particle tune (including space charge) is static in the 2D simulation due to the absence of synchrotron  124201-3 124201-3 oscillation.A more pronounced effect, with a loss regime and an emittance growth regime analogous to the experiment, though quantitatively still much weaker, is obtained assuming a 5 times stronger octupole.This is shown in Fig. 4 after 1000 turns, where the effect is practically saturated.The saturated maximum emittance growth is now about 8%, which reflects the reduced space charge detuning relative to the octupole strength.The loss regime confirms the shrinking of the dynamic aperture for Q x !6:25.Note that a comparison of this self-consistent round beam simulation with a case, where the horizontal tune shift is unchanged, but the vertical emittance is reduced to the experimental value, has shown that the emittance asymmetry has a weak effect only.
In order to explore the evolution over significantly longer times we have replaced the fully self-consistent space charge calculation by an analytical calculation for the ''frozen'' initial profile.While such analytical space charge models ignore the dynamically changing space charge force, they have the advantage of being much faster and eliminating completely the inherent noise of a PIC simulation.As a first example we have taken the above 2D coasting beam case with exactly the same transverse rms values and Gaussian distribution.The round beam approximation allows us to use the Gauss law for space charge calculation, which we have applied to the initial profile.The result after 10 3 turns-without updating the space charge electric field-is found to deviate by not more than 10% from the self-consistent simulation.Beyond this, the analytical space charge model shows practically no change of rms emittance between 10 3 and 10 5 turns.This finding gives us confidence that self-consistency might not be an important factor as long as emittance growth or loss is small enough.

D. 3D simulation with synchrotron motion
The failure of 2D simulation to describe the experiment justifies the need for including the longitudinal motion in a 3D bunch model, which induces a tune modulation primarily via space charge.In Appendix A we show that for certain classes of density profiles we can perform an exact integration of Poisson's equation using a fully analytical method.This method generalizes a technique used in Ref. [13] for 2D and parabolic profiles to the 3D case as well as more general (relatively arbitrary) density profiles.It uses only integrals of algebraic expressions in x, y, and z and is therefore sufficiently fast and noise free.We employ a density distribution following a higher order parabola of the kind 1 ÿ x 2 =a 2 3 (in x and similar in y and z), which is analyzed in detail in Appendix B. Its core is sufficiently close to that of a Gaussian, but it has a finite beam edge at 3. Using 2000 test particles we generate the corresponding (initially) consistent distribution in 6D phase space, assuming the same bunch length (200 ns at 4) and a synchrotron period comparable with that of the experiment.
The dependence on the working point is seen (Fig. 5) to have a similar trend as in the experiment for Q x > 6:28, but no loss for smaller tunes.For better comparison with the experiment we have applied a Gaussian fit to the simulation data and determined the rms emittance from this fit, which puts the emphasis on the core emittance rather than rms emittances.Note that the relatively large emittance growth without accompanying loss reflects the large physical aperture in both experiment and simulation, if compared with the initial beam size.
The resulting maximum halo increases for Q x !6:25 (Fig. 6).This is due to the fact that larger betatron amplitudes must be adopted to compensate the increasing space charge, while the particle moves longitudinally to the bunch center and trapping on the resonance island is maintained.Note that for Q x !6:32, where the resonance loses its effect, the halo shrinks to the initial beam edge of 3.
This halo formation by island trapping is the reason for beam loss in the experiment, where apparently further 0.5 nonlinearities cause a smaller dynamic aperture than in the simulation and lead to halo extraction.An example of a simulation transverse beam density cut with pronounced halo determined after 5 10 5 turns is shown in Fig. 7 on a logarithmic vertical scale.The total number of particles in the halo at this instant is about 1%.This quantity does not include the particles, which have been temporarily in the halo at some earlier time, and which would determine the total loss, if a physical aperture would intercept the halo.We will return to this issue of integrated halo in Sec.III F.

E. Interpretation
Our interpretation of the significant difference between 2D and 3D emittance growth relates this to synchrotron motion: in 2D particles have static tunes and are on resonance for basically one value of the betatron amplitude.Once on resonance they get easily detuned again after a small amplitude increase due to the dominant space charge detuning; in 3D the synchrotron motion makes the particles oscillate between high and low space charge, which induces an efficient periodic tune modulation.Eventually such particles are locked to the resonance islands, which implies that they may be carried to a larger transverse amplitude to compensate the enhanced space charge when moving back to the bunch center.A relatively large number of particles, those with sufficiently large synchrotron amplitude to reach the island tunes, is thus able to periodically cross the resonance, until trapping occurs.As was shown in Ref. [5], such trapping may be followed by detrapping after some time unless the particle hits the aperture before.This process is shown in Fig. 8, where we plot the time evolution for a single test particle of the simulation in Sec.III D with tune Q x 6:257.x / max σ The particle was arbitrarily chosen with maximum synchrotron amplitude, hence initiated at the bunch end.Units on the abscissa are horizontal single-particle emittances relative to the initial (transverse) beam edge emittance corresponding to 3. The initial relative single-particle emittance is chosen to be 0.8.The tune of Q x 6:257 is sufficiently close to the resonance so that a large halo radius may be expected, which is confirmed by the simulation.The growth of the single-particle emittance up to 5.5 times the beam edge emittance is consistent with Fig. 6.Noting that a full synchrotron period takes 625 turns, the main structure of this plot is seen to be closely linked to half synchrotron periods.The first major spike near turn 1750 indicates island trapping for a shorter fraction of the synchrotron period: after trapping the synchrotron motion carries the particle from one side of the bunch through the center (where it adopts maximum emittance) to the other side; there it gets detrapped at small emittance and remains so for about another half synchrotron period, followed by another trapping event.
The large excursion at about turn 3350 leads to immediate detrapping and the particle remains at large betatron amplitude for half a synchrotron period when it gets trapped again, apparently in the phase where it longitudinally goes through the bunch center, and is lost from the island at a somewhat smaller betatron amplitude, and so on.This whole process follows a chaotic pattern as is seen from the extended plot on the bottom of Fig. 8.
Over a larger number of synchrotron periods the particle may stay at increased amplitude -likewise at any level of amplitude -continuously over a certain interval of time, because it was detrapped after a preceding trapping.On a long run it thus frequently jumps between halo and core, which implies also a large jump in the periodic tune modulation due to the strong dependence of the betatron tune on transverse distance.This process of trapping and detrapping by an island is, in principle, closely related to previous studies of the effect of machine tune modulation in the absence of space charge [7,8,10,11,14].The main difference lies in the fact that the space charge tune modulation is exposed to these large jumps following trapping or detrapping events, hence the relative position of islands changes in time and from particle to particle.
To illustrate the process of trapping and detrapping in more quantitative terms we have analyzed the singleparticle emittance of the specific test particle of Fig. 8.In Fig. 9 we plot, for a total history over as many as 10 7 turns, the probability that the single-particle emittance exceeds a certain ''reference emittance.''The latter we have chosen as twice the emittance corresponding to the unperturbed beam edge at 3. For every crossing of the reference emittance we determine the number of turns, and the particle continuously exceeds this value.Referring to the bottom diagram of Fig. 8 we note that, for example, there is an event following turn 10 800, where x = x;edge > 2 for 2000 turns.This particular event yields a probability of n 2000=10 7 at bin 2000 on the abscissa, if such an event occurs n times over the full length of 10 7 turns (we are ignoring events lasting less than 50 turns).As expected from Fig. 8 these probabilities are clustered around multiples of half synchrotron periods.Note that the probability for the particle to be continuously beyond the reference emittance for more than six synchrotron periods is only at the 10 ÿ4 level.For less periods there is almost a linear dependence on the logarithmic scale.The integrated effect is such that this particle spends about 50% of its time outside the reference emittance.

F. Instantaneous and integrated halo intensity
The 50% total halo residence probability for the above described test particle is specific to this particle, and for a sufficiently large number of turns to get good statistics.Particles with different initial betatron or synchrotron amplitudes may take a different time to get trapped.It is therefore useful to compare, at a given time, the instantaneous halo intensity with the integrated halo intensity.The latter is defined as number of particles which have been ''kicked'' into the halo ( x = x;edge > 1) at least once up to this time, no matter for how long.The actual loss on a physical aperture limitation is given by this integrated halo intensity.The result is shown in Fig. 10 for the case with Q x 6:26, where it is seen that the instantaneous halo intensity levels off at about 1% relatively early (compare with the corresponding final density cut in Fig. 7).The integrated halo intensity, instead,

G. The rms self-consistency
The simulation rms emittance evolution for Q x 6:29 using the fully frozen analytical space charge is shown in Fig. 11.The time profile compares well with the measured data, but saturates somewhat below the experimental values.Similar to the instantaneous halo intensity in Sec.III F, this saturation is due to the fact that an equilibrium occurs between trapping and detrapping among the particles with a sufficiently large synchrotron amplitude (depending on Q x ) to cross the resonance.We compare this result with a modification, where the growing rms emittance is used to update the horizontal rms size.As a result of this rms self-consistency space charge gets weaker with progressing emittance growth.This allows even particles to cross the resonance, which have initially had too small of synchrotron amplitudes, correspondingly large tune shifts, to be able to reach the resonance condition.With this increasing number of potentially resonant particles, further emittance growth takes place, and better agreement with the measurement is achieved.While this simple modification works well in the emittance growth regime, it does not help to improve the agreement in the loss regime, where better knowledge of the dynamic aperture is needed.

IV. CONCLUSION
Synchrotron oscillations in a 3D bunch have been shown to enhance significantly the response on the octupole, if compared with the coasting beam 2D case.In the emittance growth regime quite good agreement is achieved with the measurements over half a million turns, which supports our 3D space charge model and the interpretation of the observed phenomena in terms of tune modulation caused by space charge.We predict the formation of a halo increasing in radius for Q x !6:25 and claim this is the source of the measured loss.Future measurements should consider weaker octupoles, where the predicted halo might be entirely inside the dynamic aperture and observable by scrapers.Obviously more refined measurements are desirable to further deepen the understanding of the underlying complex nonlinear dynamics processes.Theoretical efforts to explain the detailed processes have to be advanced accordingly.Cross-checks of the simulations presented here with fully self-consistent 3D simulation may be feasible over some 10 4 turns, which requires pushing the octupole excitation to the limits in the experiment (planned for the near future).

APPENDIX A: ANALYTICAL SOLUTION OF POISSON'S EQUATION
It is well known that Poisson's equation can be integrated explicitly for the special class of ellipsoidal symmetry density profiles, where the real density (normalized to 1) can be written as where is a normalization constant and F an arbitrary function.The coordinates are normalized by using the single-particle emittances as invariants, hence x x x 2 x x 02 , y ỹ y 2 ỹ y 02 , and z z z 2 z z 02 , with maximum values x , y , and z .Using these invariants to construct the distribution function is based on the approximation that the underlying motion is linear.Note that for the considered weak nonlinearities such a matching of the initial distribution consistent with the linear part in the focusing forces is common practice.
The electric field for N protons in the bunch is then given by (see Refs. [15,16]) For simplicity we limit the discussion to rotationally symmetric bunches, hence a b and T r 2 =a 2 z 2 =c 2 , with r x 2 y 2 p .We assume a power series expansion: Note that such a general ansatz allows us to consider sufficiently realistic density profiles.In previous work in connection with the 2D coasting beam case only the linear term was considered, which creates a parabolic density profile (see Ref. [13]).While we present the general formalism here, we will, later on, keep the linear as well as quadratic terms.Using Eq. (A6), we obtain By inserting this into Eq.(A5) we find the general expansion of the electric field as The starting integral is This results in the following expressions for the lowest orders: .

APPENDIX B: EXAMPLES
These expressions will be used to calculate explicitly the electric field for three examples, where the E y is obtained from E x replacing x by y.Note that a full Gaussian requires extending the series expansion of n nt to infinity.We confine the evaluation up to the second order terms, which yield sufficiently smooth density profiles.

Zeroth order
The normalized density zero order term leading to uniform real space density is taken as In this case the electric field can be written as E x 3eN 8 0 xI 1;0 ; E z 3eN 8 0 zI 0;1 : (B2)

First order
Using the first order term (similar to the 2D treatment of Ref. [15]

FIG. 1 .
FIG. 1. (Color) Experimental results on final rms emittance (of Gaussian fit) and beam current relative to initial values.

FIG. 7 .
FIG.7.Particle density profile normalized to the bunch center density in simulation for Q x 6:26 (radial units in initial ).Note the evolution of tails beyond 3.

FIG. 8 .
FIG.8.Trapping and detrapping of a single test particle and Q x 6:257.Shown are transverse relative single-particle emittances as functions of turns (top frame zoom of bottom frame abscissa).

FIG. 9 .
FIG. 9. Statistical evaluation of trapping and detrapping of a single test particle over 10 7 turns and Q x 6:257.Shown is the probability with which a given particle stays continuously outside the ''reference emittance'' for a certain number of turns.

TABLE I .
Experimental parameters.
in 2D) we have This leads to a transverse density profile nx 15=16a1 ÿ x=a 2 2 (similar in y, z), in contrast with the parabolic density obtained from the same order in 2D.The equivalent distribution function is F / 1 ÿ t ÿ1=2 .In this case the electric field can be written as This is the case used for the present paper.At the beam edge it is sufficiently smooth (with vanishing second order derivatives of density) representing a better match with a Gaussian density than the previous one.The n nt is based on the quadratic term: xI 1;0 ÿ 2r 2 I 2;0 ÿ 2z 2 I 1;1 z 4 I 1;2 2r 2 z 2 I 2;1 r 4 I 3;0 ; zI 0;1 ÿ 2r 2 I 1;1 ÿ 2z 2 I 0;2 z 4 I 0;3 2r 2 z 2 I 1;2 r 4 I 2;1 :(B6)