Scaling for faster macroparticle simulation in longitudinal multiparticle dynamics

Macroparticle tracking is a direct and attractive approach to following the evolution of a longitudinal phase space distribution. When the particles inter-act through short range wake (cid:12)elds or when the inter-particle force is included, calculations of this kind require a large number of macroparticles. However, it is possible to reduce both the number of macroparticles required and the number of tracking steps per unit simulated time by employing a simple scaling. It is demonstrated that the Vlasov equation is unchanged by introduction of the scaled quantities. It is further shown that under rather general conditions the speed of calculation improves with the fourth power of the scaling constant. Two examples comparing scaled and original cases illustrate the e(cid:11)ectiveness of the approach. Limitations to the amount of re-scaling are discussed.


I. INTRODUCTION
Multiparticle tracking programs have a long history and established utility for modeling the evolution of the longitudinal phase space distributions for particles in accelerators as the particles respond to the rf in acceleration or bunch manipulation. The macroparticle distribution can be used to approximate the evolution of the beam current distribution or fourier spectrum throughout the process being modeled. Passing from single particle dynamics to multiparticle dynamics by calculating the beam current every time step and including the contribution of the fields induced by it to the single particle motion is a direct extension of the technique which makes a wide range of interesting problems accessible.
However, the question of the number of macroparticles needed or the relevant bandwidth for quantities calculated in frequency domain requires careful attention, and it is very easy to generate spectacular spurious instabilities by excessive time steps or an insufficient number macroparticles. This question is the source and focus of what follows; recent studies of high brightness injectors and the so-called factory accelerators have given it currency. It will be shown that modeling which might appear to require of order 10 6 macroparticles and days of time on a fast computer can be scaled to the desktop and normal patience.

II. SCALING CONCEPT
The objective of the scaling sought is to allow more rapid calculation of the time evolution of the energy and rf phase of particles in a synchrotron or storage ring taking into account the contribution of beam image current, direct interparticle force, etc. In a multiparticle tracking calculation, the positions of representative macroparticles in an energy-phase plane are found by iterating a pair of first order difference equations with an iteration step almost always taken as one period of circulation of particles around the ring. The iteration step can be fractional or multiple turns, however. The interaction of the particles is calculated each iteration by an energy increment for each particle arising from the fields generated by the total beam current.
A useful approximation for the single turn map is [1] where ϕ is a phase difference between a particle and the synchronous phase φ s , likewise ε represents is the energy difference between a particle and the synchronous energy E s , i labels particles, and m labels turns. All symbols are defined in Table I.
By inspection of the map parameters it is plausible that the phase space motion can be accelerated by scaling the phase slip factor η and the potential up by the same factor.
Notice the potential is not necessarily just that from the rf system; the collective potential produced by the beam current enters identically. A trial comparing the evolution of some distribution mapped according to eqs. 1 and 2 with that obtained from the map with η and V multiplied by a scaling factor λ reveals that the distributions, i. e. plots of macroparticle locations on the energy-phase plane, are identical when the plot for time t in the un-scaled calculation is compared to time t/λ in the scaled calculation.
The obvious gain is a factor λ −1 in the computing time by speeding up the clock in the scaled calculation. However, the scaling up of the time means that frequencies associated with the motion, like the rf frequency for example, have also been scaled up. A consequence of the frequency scaling is that a resonant potential in the problem, like a higher order cavity mode for example, must be entered into the calculation with its resonant frequency scaled up by the same factor. If this is done, the result of scaled and original mappings may be indistinguishable from one another for a reasonable choice of λ. When broadband impedance or the direct interparticle forces (variously called space charge or perfectly conducting wall forces) are included, the consequences of the frequency scaling are different and very advantageous. If it has been concluded that the effects of beam pipe impedance must be covered over some particular frequency bandwidth, or, more or less equivalently wake fields down to some minimum range, it will be necessary to bin the charge distribution into bins sufficiently narrow to reflect detail in the distribution to that scale. It is easiest to formulate the requirement in the frequency domain: the bandwidth of interest must be spanned by fourier components of the beam current. If beam circulation frequency is f • and the required bandwidth is W , the charge distribution must be divided into 2W/f • bins for the finite transform. However, in the scaled system f • is λ times higher; the number of bins can be reduced by a factor λ −1 . The number of harmonics is reduced, but the frequency range covered is the same. What has been sacrificed is sensitivity to features in the frequency dependence of the longitudinal impedance Z on the scale of f • .
For example, the space charge term in frequency domain is where n is the harmonic number,

III. FORMAL ANALYSIS
The discussion based on ad hoc scaling of the mapping was not offered as proof that the scaling is valid, rather that it is plausible. Here it is demonstrated that the Vlasov equation for the original problem is unchanged by introducing the scaled variables and parameters.
Symbols used in the following are defined in Table I.
Consider a bunch for which the small amplitude synchrotron oscillation frequency ω s = ω • Q s . The steady state collective voltage experienced by a particle with phase φ is where the sum over k is taken over all particles. The single particle equations of motion in canonical energy-time variables are The rf potential is sinusoidal with amplitude V rf . The energy increment eV rf sin φ s keeps a particle with phase φ s synchronous as the guide field is changed or energy is lost to the real part of the longitudinal impedance. The total hamiltonian for the motion of all of the particles is The time evolution of the phase space distribution for the multiparticle system is described by the Vlasov equation where Ψ(ε, τ ; t) is the particle distribution function which depends on t only implicitly through ε and τ . It expresses Liouville's theorem on conservation phase space density and is valid in the absence of diffusion or external damping of the motion.
The scaling introduced into the map in Sec. II creates an apparently new physical system which will be denoted as the primed system; it is related to the original system by However, if the two systems are physically equivalent, then Ψ | λt = Ψ| t . By direct substitution one finds The Vlasov equation in the primed system is which leads immediately to Because the scaling constant λ = 0 appears to the same power in every term of the equation, the solution is independent of the scaling.

IV. COMPARISON OF SCALED AND UN-SCALED TRACKING
In Sec. II there is a comparison of a scaled and an un-scaled tracking with the collective potential generated by an h=3 HOM and 4 · 10 13 protons in the ring. A more realistic case is to add the space charge force to this example. In the earlier example 8000 macroparticles were quite sufficient because the HOM resonance is at a rather low frequency. When the space charge force is evaluated, the gradient of the charge distribution is calculated; there must be enough macroparticles to produce a smooth distribution. An insufficient number was adjusted to avoid spurious breakup. The plot in Fig. 4 was generated from tracking with the scaling λ = 2 and used 2.5 · 10 5 macroparticles; it required 1.9 hours on the same computer, just a little better than λ 4 times faster. It may not always be possible to scale this much, although in some cases even more may be possible.

V. UTILITY
If the scaling factor λ can be significantly larger than 1 without significantly reducing the accuracy of the result, much time can be saved in macroparticle modeling. One limitation on the choice of λ has been discussed above, viz., that the spacing of harmonics of the circulation angular frequency ω • should not be greater than the scale of important features of the frequency dependence of the longitudinal impedance Z (ω). Returning to the difference equations eqs. 1 and 2, one sees that, in the absence of any collective potential, the scaling by integral or 1/integer λ is equivalent to a multi-turn or fractional-turn map respectively.
This indicates that an additional feature of the scaling with λ > 1 is the introduction of an artificial increase in the effective time step. This artificial discretization error leads to distortion of the macroparticle trajectories that can be very apparent in extreme cases. The parameter to control is the synchrotron oscillation phase advance between energy increments.
When this is below 0.01π or so, the calculation should be at least qualitatively reasonable.
Another condition, which is generally the same as the limitation on synchrotron oscillation phase advance per iteration is that rf phase slip per iteration should be small with respect to a bunch length, and energy increment per iteration should be small with respect to the beam energy spread. The scaling described is so simple to implement and so innocuous in typical applications that it seems reasonable to employ it when calculation time for the un-scaled case is more than a few minutes. The precaution of considering the granularity of the sampling of Z (ω) and the adjustment of the frequency of any narrow resonances requires forethought, but the test for the appropriateness of the iteration step can easily be an automatic check in the modeling code.
The conceptual discussion above has used mostly a frequency domain description for simplicity. The two examples shown, however, are pure time domain calculations using a slightly modified version of the ESME tracking code. [4] The developmental version of ESME now incorporates λ scaling as an option, but its compatability with all existing features, especially frequency domain options, has not been verified. A public version is planned when experience has established internal consistency, robustness, and ease of use.  Note that the time scale is t = t/λ, one half of that in Fig. 3.