Anomalous beam diffusion near beam-beam synchro-betatron resonances

The diffusion process near low order synchro-betatron resonances driven by beam-beam interactions at a crossing angle is investigated. Macroscopic observables such as beam emittance, lifetime and beam profiles are calculated. These are followed with detailed studies of microscopic quantities such as the evolution of the variance at several initial transverse amplitudes and single particle probability distribution functions. We present evidence to show that the observed diffusion is anomalous and the dynamics follows a non-Markovian continuous time random walk process. We derive a modified master equation to replace the Chapman-Kolmogorov equation in action-angle space and a fractional diffusion equation to describe the density evolution for this class of processes.


Introduction
Diffusion of particle beams due to nonlinear fields is often a major source of emittance growth and beam loss in an accelerator. Measurements of diffusion coefficients have been reported from several hadron accelerators [1,2,3]. The diffusion equation was also used to explain the change in beam lifetime following the failure of a separator during a Tevatron store [4]. In collision mode the beam-beam interactions are usually the dominant nonlinearity. Diffusion coefficients in the absence of low order resonances have been calculated for head-on interactions [5] and for long-range interactions [6]. Diffusion due to nonlinear resonances is more complex and the study of this phenomenon has a long history, see e.g [7,8,9,10,11]. Resonances when modulated, either by dynamical effects such as synchrobetatron coupling or due to ripple in magnet currents, can sweep across phase space and transport particles to large amplitudes [12].
In this article we will study the nature of the diffusion process due to synchro-betatron resonances driven by beam-beam interactions with a crossing angle. This was first investigated at the DORIS collider [13] and has since been observed at other colliders. Our aim is to establish the correct statistical mechanical model that describes the evolution of the beam density. We examine the possibility that the diffusion process is anomalous with detailed tracking simulations and derive a master equation and a related fractional diffusion equation that may describe the transport process. A preliminary version of this study was reported in [14]. An example of anomalous diffusion observed in particle beams as a consequence of rf phase modulation was reported in [15]. Anomalous diffusion processes have been reported in several areas of physics including plasma turbulence [16], and in the motion of laser cooled atoms on a lattice [17].

Synchro-betatron resonances due to crossing angles
Synchro-betatron resonances (SBRs) due to beam-beam interactions at a crossing angle are convenient to study resonantly driven amplitude growth for several reasons. At large amplitudes, the non-linear force vanishes, hence particle excursions do not go to arbitrarily large amplitudes which is not the case for resonances due to multipole nonlinearities. This removes numerical instabilities and also allows the entire beam to be probed for the particle dynamics. Another advantage is that the resonances can be studied in one transverse plane since these resonances are driven by energy pumped from the longitudinal plane to the transverse plane with very little impact on the longitudinal dynamics.
When beams collide at an angle, the transverse distance of a test particle from the center of the opposing bunch depends on the longitudinal position of the particle. Consequently synchrotron oscillations of the particle couple to the transverse beam-beam force leading to excitation of synchro-betatron resonances. Since the beam-beam force goes to zero at large transverse separations, the effects of these resonances are experienced by particles only within a certain range of transverse amplitudes.
For simplicity, we choose the resonances to be in only one transverse plane, here the horizontal plane. In order to observe effects over relatively short computation times, we choose low order resonances. The tunes we choose are unrealistic for operating colliders but it is likely that the dynamics near high order resonances is similar but occurs over a longer time scale.
Linear motion and the beam-beam interactions can be described by the equations of motion resulting from the Hamiltonian where (ν x , ν y , ν s ) are the tunes, and (J x , J y , J s ) are the actions. U (x, y, s) is the beam-beam potential, δ P is the periodic delta function, φ is the azimuthal coordinate and the sum extends over the number N IP of interaction points. Assuming Gaussian distributions in all three planes, crossing angles of (2φ x , 2φ y ) in the horizontal and vertical planes respectively, the beam-beam potential for colliding proton bunches can be written as where N b is the bunch intensity of the opposing bunch, r p is the classical proton radius, γ p is the proton energy in units of its rest mass and σ x , σ y are the rms beams sizes of the opposing beam at the interaction point (IP). The potential can be expanded as a Fourier series This potential can excite synchro-betatron resonances given by the resonance condition m x ν x + m y ν y + m s ν s = p where (m x , m y , m s , p) are integers. It can be shown from the structure of the Fourier harmonics U m x ,m y ,m s that they are non-zero only when the sum m x + m y + m s is even. The Fourier harmonics can also be used to calculate the tune shifts with amplitude and the resonance driving terms, as was done in [18]. As one example, we write down the zero transverse amplitude tune shift for round beams. This tune shift now depends on the longitudinal oscillation amplitude a s σ s as and a similar expression for ∆ν y . Here ξ = N b r p /(4πε N ) is the usual beam-beam parameter, (a x σ x , a y σ y ) are the transverse amplitudes of the particle, I 0 , I 1 are modified Bessel functions and the other dimensionless parameters are As a consequence, only those zero transverse amplitude particles with zero longitudinal amplitude a s experience the full beam-beam tune shift ξ . Particles with non-zero amplitude a s experience a smaller tune shift.
Since the LHC employs crossing angles in its collision scheme, we will use the LHC beam parameters in the simulations reported here. As in the LHC, the crossing angle is in the horizontal plane at one IP and in the vertical plane at the second IP. We consider resonances excited in the horizontal plane only, so they are of the form m x ν x + m s ν s = p with m x +m s even. In our model the only sources of tune spread are the beam-beam interactions. These interactions between protons lowers the betatron tunes at small amplitudes. We choose the large amplitude tunes, i.e. the tunes with only the linear lattice, to satisfy one of the SBR resonance conditions. Having chosen a particular resonance m x ν x + mu s ν s = p to be satisfied by the bare lattice tunes, the tunes inside the bunch are determined by the beambeam parameter ξ , the synchrotron tune ν s and the amplitudes (a x , a y , a s ) of the particle. The nominal LHC horizontal tune is 0.31 at collision, so we searched among the following resonances: 3ν x ± ν s = 1, 2(3ν x ± 2ν s ) = 2 as well as 2(4ν x ± ν s ) = 2 and 4ν x ± 2ν s = 1 to find those that cause large growth of the emittance and beam tails. Given that the betatron tune spread from head-on beam-beam interactions is about 0.007 and the small amplitude synchrotron tune is ∼0.002, the choices 2(3ν x −2ν s ) = 2 and 4ν x −2ν s = 1 had the greatest impact on the beam. With these choices, low amplitude particles are resonant with the third and fourth order betatron resonances respectively, and the synchrotron oscillations modulate these resonances leading to large amplitude growth. The other resonances are resonant at larger amplitudes and consequently have a smaller impact on the bunch. The bare lattice (which become the large amplitude) betatron tunes corresponding to these resonances are shown in Table 1. Some of these parameters may be slightly different from the present LHC design values, e.g the LHC design value of the crossing angle is 285µrad.

Simulations of beam variables
In this section we will describe multi-particle simulation results. These will include the emittance growth, evolution of beam profiles, amplitude growth at different initial amplitudes, and also the growth of the variance in action at these initial amplitudes. This will allow us to probe both the macroscopic and microscopic beam behaviour.
The simulations were performed with a simple numerical model consisting of six dimensional linear transport between the two collision points, a sinusoidal longitudinal map through an rf cavity and weak-strong beam-beam interactions at the two IPs. The beambeam interactions occur with a horizontal crossing angle at one IP and a vertical crossing angle at the second IP. The strong beam was assumed to have a Gaussian distribution in all three planes. Magnetic nonlinearities are not included, both to keep the model as simple as possible and also to avoid particle amplitudes from growing exponentially fast far from the beam core. Limiting amplitude growth to finite values allows us to keep all particles in the distribution and hence study the growth of the beam tails with good statistics.

Emittance growth and lifetimes
Emittance growth was calculated by evolving ensembles of N particles (5000 ≤ N ≤ 20000) starting with Gaussian distributions in all planes. Typically 10,000 particles sufficed to obtain results that did not change much with a larger number of particles. The calculated emittance was the rms emittance, e.g. ε Figure 1 shows the emittance growth with 20,000 particles on the two resonances. We find that the growth follows a simple power law, the fits are also shown in the figure. We observe that the horizontal emittance growth after 10 6 turns is more than 2.5 times larger on the 2(3ν x − 2ν s ) = 2 resonance than on the 4ν x − 2ν s = 1 resonance. The vertical emittance growth is much smaller than the horizontal, about a factor of five smaller for the first resonance and it is practically zero for the second resonance.
By imposing a finite aperture restriction, we can find the escape time needed by particles to reach this aperture. This has been calculated for several different apertures and for both resonances. Apertures were placed from 5σ to 10σ at intervals of 1σ . On the 2(3ν x − 2ν s ) = 2 resonance, we find that about 7% of particles reach 8σ , a handful reach 9σ and none reach 10σ . On the 4ν x − 2ν s = 1 resonance, about 4% of particles reach 6σ , a few reach 7σ and none reach 8σ . The amplitude distribution of the particles reaching 8σ on the first resonance and of the particles reaching 6σ on the second resonance are shown in the top plots of Fig 2. The initial distribution in each case was a Gaussian with 40,000 particles. On the 2(3ν x −2ν s ) = 2 resonance, the maximum of the amplitude distribution occurs close to 1.5σ -an amplitude close to the lower edge of the resonance islands, shown later in Fig  10. The minimum amplitude that reaches the aperture is 0.25σ . On the 4ν x − 2ν s = 1 resonance, the corresponding peak in the amplitude distribution is close to 1.8σ , also at the lower edge of the resonance islands seen in Fig. 11. The minimum amplitude that reaches the aperture on this resonance is 0.9σ . cape time with 20,000 particles yielded similar values showing that these numbers have converged to stable values. The average escape time increases by an order of magnitude or more for each increase in aperture by 1σ . The average escape time at 8σ on the 2(3ν x − 2ν s ) = 2 resonance is about the same as at 6σ on the 4ν x − 2ν s = 1 resonance. At a fixed aperture, the differences in escape times between the two resonances increases by about two orders of magnitude at 5 and 6 σ and three orders of magnitude at 7σ . One would expect this trend of increasing lifetimes to continue with higher order resonances.

Beam profiles
The beam profiles were found for the same distributions and resonances. The left plot in Fig 3 shows a mountain range view of the horizontal beam profiles (i.e. distribution function of the horizontal position), initially and then at other intervals up to 10 6 turns with tunes on the resonance 2(3ν x − 2ν s ) = 2. After the initial time, the subsequent horizontal profiles develop long non-Gaussian tails which extend out to±8σ compared to the initial Gaussian distribution which was limited to ±3.5σ . The vertical beam profiles (not shown here) however stayed Gaussian and close to the initial distribution. The right plot in this figure shows the horizontal profiles but with tunes on resonance 4ν x − 2ν s = 1. We observe that in this case as well that the tails are non-Gaussian and extend out to about ±6σ , not quite as far as on the first resonance. Again there is very little change in the vertical profile. We observe that the beam tails do not appear to change very much after the first 100,000 turns or so. It is most likely that the regions of enhanced diffusion are depleted within these turns. The particles in the vicinity of the resonance islands are transported to larger amplitudes quickly and are detuned from the resonance. The amplitudes to which they move have much smaller diffusion, so the beam tails do not change much. As we will see in the next subsection, the evolution in the beam core shows growth even after several hundred thousand turns. However these particles do not migrate to the tails during the time duration followed. Thus we continue to observe emittance growth.
In order to find distributions that can best fit the non-Gaussian tails, we first look to the Central Limit Theorem (CLT) which explains the ubiquity of the Gaussian distribution. This powerful theorem states that the distribution of a sum of a sequence of random, identically distributed and independent variable with finite mean and second moment tends to a Gaussian distribution in the limit that the number in the sequence approaches infinity. Generalizing the CLT by dropping the requirement of a finite second moment leads to the family of Levy stable distributions [20]. For applications in beam dynamics, these distributions will still have a finite second moment because they do not extend to infinity but are truncated at the beam pipe or the closest physical apertures.
Levy stable distribution functions are defined by an inverse Fourier transform of a stretched exponentially decaying function in Fourier space There is no known closed form expression for arbitrary values of α. Special cases include: the Lorentz distribution L 1 (z) while L 2 (z) is the Gaussian distribution. There are more general asymmetric versions of the Levy stable distribution with additional parameters but we shall not need them here. Some basic properties of these functions are [21] • These functions are normalized : ∞ −∞ dzL α (z) = 1 • They are even functions : , which increases rapidly when α → 0.
• At large values of z, the distributions decay as We find that the non-Gaussian horizontal profiles can be fit by these Levy stable distributions L α . The left plot in Fig. 4 shows the fit of the final horizontal profile for the resonance 2(3ν x − 2ν s ) = 2 with a Levy stable distribution with parameter α = 0.95. This profile is narrower than a Lorentzian and decays at large x as |x| −1.95 . The right plot in this figure shows the final distribution on the resonance 4ν x − 2ν s = 1 can also be fit by a Levy stable distribution with a larger central width and corresponding to α = 1.3. This profile is wider than a Lorentzian and decays at large x as |x| −2.3 . The Levy stable distributions were generated with a Mathematica package [19].
It is known [22] that the Levy stable distributions serve as Green's functions to fractional diffusion equations for a density ρ(x,t) of the type where χ is a constant diffusion coefficient and −∞ D α x is the Riemann-Liouville fractional space derivative of order α given by, The solution of the fractional diffusion equation above is where ρ 0 (x) is the initial density. Levy stable distributions have also been shown to be solutions of other fractional diffusion equations [23]. There is no reason to believe that either respectively. Initially 4000 particles each were placed at horizontal amplitudes from 0.5 to 4σ in steps of 0.5σ . The vertical amplitude was kept constant at 0.1σ . In the plots the vertical scale has been truncated to 400 in order to show clearly the particle numbers in the final distribution at large amplitudes.
Equation (6) or of the type in reference [23] are appropriate for our problem. However the fact that the long time beam profiles are described by these Levy distributions is our first indication that the amplitude growth process may be described by an appropriate fractional diffusion equation rather than the regular diffusion equation. In Appendix A we derive a different fractional diffusion equation that may describe the dynamics observed here.

Growth at individual amplitudes
We now take a closer look inside the beam distribution to determine how the amplitude growth changes with amplitude. Instead of a Gaussian distribution in phase space, we consider delta function distributions in action. We select a discrete number of horizontal actions and at each action we place 4000 particles uniformly distributed in angle. The vertical amplitude was kept constant at 0.1σ for all particles. The initial distribution in transverse action angle space can be written as where J 0.1 is the action at an amplitude of 0.1σ , P(θ x ) is a uniform distribution in the horizontal angles etc. The initial longitudinal variables were chosen to be the same for all particles: z = 1σ s , δ p/p = 1σ p in these simulations. We let these distributions evolve and record the final distribution in amplitude after 10 6 turns. The left plot in Fig 5 shows the initial (red) and final (blue) distributions for resonance I. We observe that particles at 0.5σ stay close to their initial amplitude. At 1σ , many particles have moved to larger amplitudes but a sizable fraction stay in their original neighbourhood. This shows a large variation in final amplitude depending on their initial angle or sensitivity to their initial conditions. It suggests that motion in the neighbourhood of 1σ could correspond to bounded chaos.  showing that diffusion at these amplitudes is weaker than in the first resonance.

Variance of the action and diffusion type
We now examine the diffusion from individual amplitudes. In regular diffusion the variance of the diffusing quantity, here the action, grows linearly with time which allows one to define time independent diffusion coefficients D(J) = (∆J) 2 /∆t. We check the validity of this assumption for the beam-beam driven SBRs. using the same initial distributions as used in Fig 5. Variances are calculated over particles at the same initial action. Figure 6 shows the growth in the variance of the horizontal action at several initial actions for both resonances. The vertical amplitude was constant at y = 0.1σ . Initially the variance is zero at all actions but then grows at different rates depending on the action. The growth in the variance is not linear at any action. In most cases there is a sharp initial transient growth which is followed by a slower long term growth. This long term growth can be modeled (again in most cases) by a power law behavior of the form where the coefficients (C x ,C y ) and the powers (p x , p y ) depend on the initial action. Exponents less than 1 indicate sub-diffusive behavior while exponents greater than 1 imply super-diffusive motion Figure 6 also shows the fits with this power law. Growth of the variance in the vertical action can also be fit by a single power law with small values of (C y , p y ) showing that there is no appreciable diffusion in that plane. On the resonance 2(3ν x − 2ν s ) = 2, there is significant growth in the action at amplitudes of 2 and 2.5σ compared to neighboring actions both lower and higher. The exception to the single power law fit occurs at x = 1σ where the variance stays nearly constant after the initial transient and then after about 400,000 turns grows by an order of magnitude over the next 600,000 turns but with oscillations in the variance. These oscillations occur because of the large sensitivity to the initial angle at this amplitude. The oscillations decrease significantly when the number of particles at the same initial action is increased from 4000 to 20000 particles, which results in a more complete sampling of the initial angle. Simulations show that this greater sensitivity to the initial angle is also present at amplitudes in the range 1.0σ ≤ x ≤ 1.3σ with y = 0.1σ . The fits to a power law in this zone are applied after the variance starts to grow rapidly but with 20000 particles. The average action with initial |x| = 1.0σ grows about 10% after 10 6 turns while the average action with initial |x| = 1.5σ grows by about a factor of two over this time. So the narrow zone around |x| = 1.0σ corresponds to a zone of bounded chaos.
At the resonance 4ν x − 2ν s = 1, the growth in variance is largest in the range x = 2.5 − 3σ and drops for both smaller and larger initial actions. The large oscillations in the variance occur in a range around x = 2.0σ and again these oscillations are reduced when the number of particles is increased from 4000 to 20,000. For this resonance, the zone around |x| = 2.0σ is a zone of bounded chaos. Similar behaviour is seen at other values of y but the width of the zone of bounded chaos changes.
The exponents in the power laws were calculated for several values of the horizontal amplitude and for different vertical initial amplitudes. Fit 7 shows the exponents for both resonances. On the resonance 2(3ν x − 2ν s ) = 2 there is a spike in the exponent to values well above 1 in the regions of bounded chaos for y = 0.1, 0.5σ suggesting super-diffusive behavior. Above the zone of bounded chaos, the exponent falls well below 1 suggesting sub-diffusive behavior. At y = 1σ the exponent stays well below 1 for all x showing that zones of bounded chaos have disappeared. On the 4ν x − 2ν s = 1 resonance, the exponent rises above 1only in a narrow zone around x = 2σ at y = 0.1σ . At y = 0.5σ the exponents stay well below 1 at all x with a small spike at x = 2σ . The motion is sub-diffusive at all x values studied when y = 1σ . Since the super-diffusive regions are narrow, it is possible that they may appear for |y| ≥ 1σ when the motion is studied with a finer resolution or even when the longitudinal variables are changed. We remark that we have observed here three different signatures of bounded chaos: large variations in final amplitude when starting from the same initial amplitude (seen in Fig 5), large oscillations in the action variance over time (seen in Fig. 6 and a spike in the power law for the growth of the variance (seen in Fig. 7). These signatures apply to an ensemble of particles at the same amplitude but different initial angle as opposed to the Lyapunov exponent criterion which is applied to a pair of particles that are initially infinitesimally close.
The picture that emerges is that near synchro-betatron resonances, phase space is divided into several zones. At small amplitudes there is no diffusion. At larger amplitudes there is a zone of bounded chaos with super-diffusive motion. The next zone outward in phase space is wider with sub-diffusive motion. Finally at even larger amplitudes, the motion becomes linear again and consequently there is no diffusion. Fig. 8 shows a qualitative sketch of these different zones. The width of the super-diffusive zone with bounded chaos depends on the resonance, on the amplitude of the orthogonal transverse amplitude and on The fact that the sub-diffusive regions seem to be dominant in this perturbed Hamiltonian system should not be unexpected due to the existence of hyperbolic fixed points and the existence of perturbed KAM tori. These fixed points and tori lead to orbits which stay in their vicinity for long time periods and consequently to slower growth. Similar phenomena have been reported for the standard map by Balescu [24].

Statistics of single particle behavior
We saw in the previous section that in most regions of phase space, the variance grows slower than linearly with time. If we define an instantaneous or 'running' diffusion coefficient [24] as D J x = (1/2)∂ ∆J 2 x /∂t, then this coefficient would be time dependent and would vanish in the very long time limit. Near both resonances we did not observe any zone of regular diffusion with constant diffusion coefficients. We also saw that the beam profile was given by a Levy stable distribution which is known to be the solution of a fractional diffusion equation. These suggest that the dynamics near these resonances cannot be described by the regular diffusion equation but instead that the diffusion is anomalous which needs a different diffusion equation. In order to test this possibility in more detail, we will examine the validity of the assumptions behind the regular diffusion equation.

Continuous Time Random Walks
The regular diffusion equation arises after assuming that the particle dynamics can be modeled as a classical random walk following a Markov process. This implies that particle jumps occur at regular time intervals and there is a well defined time scale such that events A well known alternative to the standard random walk picture is the Continuous Time Random Walk (CTRW) model introduced by Montroll and Weiss [25] to consider processes where both the times at which jumps occur as well as the sizes of the jumps in space are random functions. A review of CTRW and connections to fractional diffusion equations may be found in [26].
A general dynamical process may not have a characteristic time scale. In those cases a Markov description may not be applicable. The CTRW model introduces the concepts of a probability distribution w for the waiting times before a jump occurs and a probability distribution Ψ for the size of a jump In beam dynamics there is no diffusion when the motion is linear and the usual Courant-Snyder actions are conserved. Consequently it makes sense to define the jumps in action space when the motion is nonlinear and diffusive. Hence we define w(t, J)∆t to be the probability that a particle waits for a time between t and t + ∆t at action J before making a jump. and define Ψ(∆J; J,t)∆J to be the probability of making a jump by ∆J at the action J at time t. These distributions are normalized, i.e.
The concept of a waiting time endows the system with memory. The CTRW model reduces to the classical random walk model on which the regular diffusion equation is based, when the waiting time follows an exponential behavior in time e −t/τ with a characteristic time scale τ.
These waiting time and jump size distributions can be used in many cases to determine the evolution followed by the density distribution ρ(J,t). The canonical CTRW model assumes a power law waiting time distribution, a Gaussian for the jump size distribution and a constant diffusion coefficient. These lead to a fractional diffusion equation for the density [26]. In our case the dynamics near resonances is sufficiently complicated that we need to establish the evolution equation for the density from first principles. We therefore need to determine the forms of the jump size distribution and the waiting time distributions from the dynamics. Simulations discussed in the rest of this section are used to extract these distributions.
A check of the whether the CTRW model may be applicable here can be done by examining the time series of single particles. Fig 9 shows one example of a time series of the amplitude 2β x J x for a single particle on the resonance 2(3ν x = 2ν s ) = 2. The left plot shows that a particle may perform small amplitude quasi-periodic oscillations for a while before a major qualitative change occurs. The middle and the right plots show that step sizes can be large (several σ ), of varying amplitude, and there are intermittent sequences of varying duration where there are smaller steps. The time dependent behaviour of this sequence and the non-locality of the changes establish that this is a process with a distribution of waiting times and a distribution of action step sizes, the key ingredients of the CTRW model.

Jump size distributions
We now calculate the jump size distributions by following a single particle for 10 6 turns and find the changes ∆x, ∆J x in position and action per turn. Fig. 10 shows the phase space, and jump distributions of ∆x, ∆J x on the resonance 2(3ν x − 2ν s ) = 2 with initial values of x = (0.2, 2, 8)σ . At the smallest initial position x 0 = 0.2σ , the phase space is a distorted ellipse with no trace of the resonance island; motion here is quasi-linear. The plot for the distribution function of ∆x also has the distribution function for a periodic function shown in dotted lines. When the argument of a periodic function like sine or cosine is sampled from a random distribution, the distribution function for the periodic function f has the

Number
Step size in relative action form The distribution function has local maxima wherever the function itself becomes stationary, so that many more points are sampled from the neighbourhood of these stationary points.
Since the motion at small and large amplitudes is quasi-periodic in our model, it is to be expected that the distribution in ∆x is close to that of a periodic function. The distribution function for ∆J x is plotted on a semi-log scale and shown as discrete points, for greater clarity. At x 0 = 0.2σ , the distribution for ∆J x lies on a single curve but not given by any simple expression. As the particle's initial position increases to 2 σ , the nonlinearity of the beam-beam force manifests and we see resonance islands in phase space and large excursions. The distribution function for ∆x undergoes a qualitative change to resembling a parabolic curve but with a dip in the center and with peaks close to the center. The distribution function for ∆J x now falls on two separate curves. Similar distributions for ∆x, ∆J x are seen for initial particle amplitudes in the range 1.5σ ≤ |x 0 | ≤ 6.5σ . At x 0 = 8σ , the phase space returns to a distorted ellipse with considerable smear, and the distribution functions also resemble those seen at x 0 = 0.2σ .

Number
Step size in relative action

Waiting time distributions
The waiting time distribution is the important distribution that determines the nature of the diffusion process. As remarked earlier, a waiting time distribution that follows an exponential law reduces to a Markov process, otherwise the process is non-Markovian. The waiting time for each initial amplitude is found here by tracking a particle at that amplitude for 10 6 turns. The phase space region in action angle coordinates that is visited by the particle is divided into different zones and the time that the particle stays in the zone before leaving is one instance of the waiting time. The choice of the width of the zone is somewhat arbitrary since there is no dynamics dependent action scale which is applicable to all of phase space. For example, the resonance width is not relevant at small or large amplitudes and if there were multiple resonances, there would be multiple widths. We therefore calculated the waiting time distribution twice, once with a chosen width such that there was enough statistics in each zone and the second time with twice the width. In most cases we found  that the parameters of the distribution change by less than 10%; we take this to be a sign of convergence of the distribution. We find that the exponential function is not a good fit to the distribution for either resonance. The results for a fit to a power law distribution are shown in Figure 12. The distributions are plotted on a log-log scale for several initial amplitudes where there is significant amplitude growth. On the 2(3ν x − 2ν s ) = 2 resonance, most of the points (with the exception of the single occurrence events with long waiting times) lie on straight lines showing that a power law is a reasonable fit. The power law exponents for the different amplitudes are close. For the amplitudes shown in this figure, the waiting law distributions are On the 4ν x − 2ν s = 1 resonance, the waiting time distribution can also be fit by a power law distribution but the range of variation in the exponent α is larger: 1.4 ≤ α ≤ 2.7. The greater variability in the exponent is expected to have an impact of the dependence of the diffusion rate at different amplitudes on this resonance.

Fractional diffusion equation
Since the waiting time distribution suggests that the transport near resonances is non-Markovian, we need to establish an alternative to the regular diffusion equation. For a Markov process, the regular diffusion equation is obtained from the Chapman-Kolmogorov master equation, a derivation is sketched in Appendix A. In Appendix A we also derive a different master equation using general jump size and waiting time distributions for a CTRW process in action space following a method outlined in [27]. The master equation for the density in action angle space that we obtain is where L t is an integral operator given by Here L −1 is an inverse Laplace transform, τ is a time parameter in the waiting time distribution w(t, J),ŵ(s; J),ρ(s; J) are the Laplace transforms of the waiting time distribution and the density respectively. In the Appendix we then show that expanding this master equation in a Taylor series in the same manner as is done for the Chapman-Kolmogorov equation, the following fractional diffusion equation is obtained for a power law waiting Here the exponent α depends on the action J which will be true in general and D kl are action dependent diffusion coefficients, defined in the appendix. It remains to be verified that this fractional diffusion equation describes the dynamics near resonances, as seen in the particle tracking simulations. However, this diffusion equation has been derived under general considerations of a CTRW process which the dynamics near the SBR resonance appears to follow. Given the large variations in the diffusion coefficients, the solution of this diffusion equation will likely require a special purpose numerical algorithm. The density can then be used to calculate the beam lifetime and various moments such as the emittance .

Discussion
We have studied the detailed transport process near two low order horizontal synchrobetatron resonances driven by beam-beam interactions at a crossing angle. We found that the horizontal beam profiles develop long beam tails. The horizontal beam distribution evolves from an initially Gaussian distribution to a Levy stable distributions on both resonances. The Levy stable distributions are solutions of simple fractional diffusion equations which describe some anomalous diffusion processes. The evolution of the variance in action at several initial values characterizes the nature of the diffusion in phase space. At small amplitudes there is no diffusion, then there is a narrow region where the motion is super-diffusive (the variance grows faster than linearly with time), followed by a broad region where the motion is sub-diffusive (the variance grows slower than linearly with time) and finally no diffusion at large amplitudes. The width and the location of the superdiffusive region depends on the resonance, the width is narrower for the weaker resonance. This super-diffusive region is also marked by signatures of bounded chaos and particles do not experience large amplitude growth. For both resonances, this region is located at the lower edge of the resonance islands. The broad sub-diffusive region abuts the superdiffusive region and continues until about 5-6σ depending on the resonance. Here particles do migrate to larger amplitudes. We do not observe regular diffusion anywhere in phase space on either resonance with the particle distributions we used.
The jump size distribution and the waiting time distribution, key ingredients of a continuous time random walk process, were found by analysis of single particle tracking data. The jump size distributions for both resonances were similar -in the linear regions of phase space, the distributions in ∆x are close to the arcsine distribution while in the nonlinear regions they have a more complex shape. The similarity of these distributions for the two resonances suggests that these may be universal features near such resonances. When the waiting time distributions follows an exponential law, the stochastic process is Markovian. We find that the waiting time distribution follows instead a power law, again for both resonances. Since the process is non-Markovian, the regular diffusion equation cannot be used to describe the evolution of the density. For a general CTRW process, we derived a master equation in action-angle space which is applicable to processes with arbitrary jump size and waiting time distributions. A fractional diffusion equation was derived from this master equation. Numerical solutions of this diffusion equation will allow computations of beam observables such as lifetimes and emittance growth.
This model can be tested against beam observations when anomalous diffusion is suspected. Comparison of beam profiles with Levy stable distributions would be a first check. Another indicator would be if the emittance of pencil beams grow nonlinearly with time. This could then be followed by measurements of diffusion coefficients at different amplitudes, using them in the fractional diffusion equation and comparing the numerically calculated emittance growth and beam lifetime with the measured values.
In this article we considered low order synchro-betatron resonances so as to observe effects on a short time scale. Based on comparisons of the two resonances studied here, we expect that the physics at high order resonances (and hence more applicable to operational accelerators) will be similar but on longer time scales. When multiple such resonances are present simultaneously, the diffusion is likely to be anomalous but the phase space dynamics will be more complicated. It is possible that the physics near space charge driven resonances may be similar to that obtained here but that remains to be investigated.
We now consider a more general master equation for a CTRW process in action angle space with arbitrary jump size and waiting time distributions. We use a method outlined in [27]. It uses two basic balance conditions: the first states that a change of density arises from the difference in the incoming flux Γ + (J, θ ,t) and the outgoing flux Γ − (J, θ ,t). The outflux at (J, θ ,t) has contributions from particles that were present initially but left after waiting for time t and those that arrived later before leaving We assume that the same relation as in Eq. (A.4) between the drift and diffusion coefficients holds. Then we have as the modified diffusion equation In cases where (1/τ)L t ρ = ρ, this is the regular diffusion equation.