Simple method for particle tracking with coherent synchrotron radiation

Coherent synchrotron radiation (CSR) is of great interest to those designing accelerators as dri free-electron lasers (FELs). Although experimental evidence is incomplete, CSR is predicted t potentially severe effects on the emittance of high-brightness electron beams. The performanc FEL depends critically on the emittance, current, and energy spread of the beam. Attempts to in the current through magnetic bunch compression can lead to increased emittance and energy sp to CSR in the dipoles of such a compressor. The code ELEGANT was used for design and simulation the bunch compressor [M. Borland et al., in Proceedings of the 2000 Linear Accelerator Conference, Monterey, CA (SLAC, Menlo Park, CA, 2001), p. 863] for the low-energy undulator test line (LEUT FEL [S. V. Milton et al., Phys. Rev. Let.85, 988 (1999)] at the Advanced Photon Source (APS). order to facilitate this design, a fast algorithm was developed based on the 1D formalism of and co-workers [E. L. Saldin, E. A. Schneidmiller, and M. V. Yurkov, Nucl. Instrum. Methods P Res., Sect. A398, 373 (1997)]. In addition, a method of including CSR effects in drift spaces lowing the chicane magnets was developed and implemented. The algorithm is fast enough to running hundreds of tolerance simulations including CSR for 50 000 particles. This article describ details of the implementation and shows results for the APS bunch compressor.


I. INTRODUCTION
It has long been known that electrons traveling through a bending magnet emit radiation.Only recently has it become widely accepted that, when very short electron bunches travel through a dipole magnet, they can emit coherently at wavelengths that are comparable to the bunch length and can propagate in the vacuum chamber.A number of authors (e.g., [1]) have delved in detail into the effects of coherent synchrotron radiation (CSR) on the electron beam.Here, we describe in simple terms why CSR has an impact on the electron beam.
Because electrons traverse the dipole on an arc of a circle, radiation from one part of the arc can catch up with electrons on another part of the arc.This, together with the fact that the emission is coherent and at a wavelength comparable to the bunch length, results in a modulation of the energy along the bunch.This is similar to a wakefield but, unlike wakefields, CSR affects the particles ahead of the emitting particles rather than behind.(For convenience, in this paper the phrase "CSR wake" is used in spite of this difference.) Since this is happening inside a dipole magnet, the energy modulation results in a modulation of the transverse slopes in the bending plane, which increases the projected transverse emittance in the bending plane.The mechanism behind the emittance growth per se is the same as that for emittance growth due to quantum excitation, i.e., variation among particles of the loss of energy in a dispersive region.Of course, quantum excitation and CSR are very different physical phenomena.One important difference is that quantum excitation is random and hence increases the slice emittance as well as the projected emittance.
The CSR model used by ELEGANT [2] is based on an equation [3] for the energy change of an arbitrary linecharge distribution as a function of the position in the bunch and in a dipole magnet.The effects of changes in the longitudinal distribution within a dipole are included, with the following limitation: ELEGANT computes the CSR effects at a given location in the dipole based on the longitudinal distribution at that location, meaning that retardation effects are not properly included.In addition, because Saldin's result is for a line-charge distribution, this model does not include the effect of the transverse distribution of the beam on the CSR nor does it include the variation of the CSR across the beam.The model does include the effect of the transverse beam distribution on the amount of emittance growth due to the energy modulation.
CSR in drift spaces is included by propagating the terminal CSR wake in each bend through the drifts with the beam.This method has not been derived in a rigorous fashion, but arguments supporting its plausibility are given below.
Most prior work on simulation of beam transport with CSR has involved detailed and computer-intensive computation of the electric field including retardation effects.See, for example, the work of Dohlus and Limberg [4] and Li [5].An effort similar to the present work but using a simplified formulation and an assumption of Gaussian longitudinal distributions was reported for PARMELA by Dowell and O'Shea [6].This latter work did not attempt to include CSR in the drift spaces following dipoles.

A. Theoretical basis
We begin this section with a recapitulation of the relevant results from Saldin et al.The rate of change of energy in the bunch is given by where R is the bend radius, f is the angle into the bend (i.e., the average electron position), s is the position within the bunch, and t is time.T 1 is responsible for most of the CSR effect and transitions into the steady-state result.It is given by where l͑z͒ is the linear charge density, s L Rf 3 24 is the slippage length, and ͑3R 2 ͒ 1͞3 .The slippage length is the difference between the path length from the dipole entrance traveled by an electron (on a curved path) and the straight-line distance traveled by radiation emitted at the entrance.The equation for T 1 reflects the fact that a test electron can only be affected by radiation from electrons less than a slippage length behind it.Radiation from electrons farther behind has not had time to reach the test electron.As the beam travels farther into the dipole, the slippage length becomes longer.Eventually, if the dipole is long enough, the slippage length will become long enough that each electron is affected by radiation from all the electrons behind it.
The other term, T 2 , is an entrance transient that eventually dies out in a sufficiently long dipole.It is given by Again, for a sufficiently long dipole, s L will eventually become large enough that T 2 is zero.Stated another way, the component of the wake due to T 2 gradually moves ahead of the beam and eventually does not affect it.To facilitate testing, ELEGANT supports a "steady-state" mode in which, effectively, s L ! `.This means that T 2 !0 and the lower limit on the integral for T 1 goes to 2`.

B. Implementation
The implementation of these results in ELEGANT relies on splitting each dipole into a user-specified number of pieces.Tests show that 60 pieces is adequate, though typically 100 are used.For each piece, the following steps are performed: (i) Propagate the entire beam through the piece using a second-order or fourth-order canonical integrator [7].The order is user selectable, but fourth-order is pre-ferred.(ii) Compute the CSR wake.(iii) Apply the CSR energy kicks.
Steps (ii) and (iii) are detailed below.A more accurate implementation would include the CSR energy kicks in the canonical integration.However, this would have complicated the code and is probably unnecessary in any case.Because the canonical integration through the dipole fields includes the energy dependence of the transport to all orders, the dispersive effects on the transverse coordinates of the particles are computed automatically.
Computation of the CSR wake is performed as follows: (a) Particle arrival times at the end of the dipole piece are binned; (b) this density histogram is smoothed using fast-Fourier transform (FFT) convolution with a Savitzky-Golay filter [8]; (c) the same filter is used to take the derivative of the smoothed density distribution; (d) the T 1 and T 2 functions are computed for each bin.Simple sums are used for the integral required by T 1 ; (e) each particle's energy is changed by dE cdt Ds for the bin it occupies, where Ds is the central path length of the dipole piece.
One difficulty in these computations comes from noise in the linear density histogram due to the use of a finite number of particles and a large number of bins.This is a particular problem when taking the derivative of l.Smoothing is used to overcome this problem, at the expense of some loss of information (e.g., real spikes in the current will be reduced).Tests show that using at least 40 000 particles and 3072 bins is required.
For relatively smooth input distributions and modest bunch compression, a first-order Savitzky-Golay smoothing filter is used with a half-width of 50 bins.This reduces noise effects and reproduces the results of runs with much larger numbers of particles.
However, smoothing must be used cautiously when spikes are produced in the longitudinal distribution, which happens when employing high compression factors or with unsmooth input distributions, as discussed below.In these cases, the most damaging effects occur due to the CSR produced by the spikes and it is important not to smooth the distribution excessively.In such cases, spurious emittance growth due to noise is negligible compared to the emittance growth due to the CSR from the spikes.We have found that for full-compression and overcompression studies of the Advanced Photon Source (APS) bunch compressor, use of about 40 000 particles and a Savitzky-Golay filter half-width of six bins is adequate to reproduce results of 400 000 particle runs with a filter half-width of one bin.Users of the code are cautioned that it is necessary to vary both the number of particles and the amount of smoothing until convergence is obtained.The parameters used here serve only as a rough guide.a Gaussian input beam with 100 000 particles, charge of 1 nC, and an rms bunch length of 50 mm.The bending radius was 1.5 m.This matches the parameters used in the example in Saldin et al.

C. Tests of the implementation
First, we compared the dE cdt functions computed by ELEGANT and by numerical integration.This comparison shows excellent agreement, as seen in Fig. 1.Differences appear to be due to deviation of the distribution from Gaussian in ELEGANT, which hints at the importance of using a detailed beam distribution in the computations rather than analytical results based on an assumption of a Gaussian line density. The where L b is the length of the dipole and s s is the rms bunch length.ELEGANT reproduces these results very accurately, as seen in Figs.2-5.

III. SIMULATION OF CSR IN DRIFT SPACES
The simulation of CSR in drift spaces is as important as simulation in dipoles, yet has garnered less attention in the literature.The approximate method used in ELEGANT is based on the assumption that the terminal CSR wake from the dipole propagates with the beam through the elements following the dipole.The code assumes that the wake does not change shape and that it moves at the speed of light, while the beam moves at a lesser speed.Support for these assumptions is given by the work of Dohlus and Limberg [9].In particular, their Fig. 3 shows the results of detailed computations of CSR in a drift space.One sees that the wakes at each subsequent location in the drift space have nearly the same shape but diminished intensity.
In practice, drifts following dipoles are split into pieces of about 1 cm in length in order to integrate the effect of CSR.If nondrift elements (e.g., quadrupoles) are in the beam line, the effect of CSR in those elements is deferred until the beginning of the next drift space.This was done to simplify the implementation.
Given that it would be unphysical for the effect of the wake to continue indefinitely, ELEGANT supports several methods for modifying the intensity of the wake.
Two of these are of interest here: use of the "overtaking length" as the exponential attenuation length and use of analytical expressions given by Saldin, Schneidmiller, and Yurkov [3].
The overtaking length [9] is essentially the length required for radiation from trailing particles to overtake and interact with leading particles.It is given by the expression L o ͑24s s R 2 ͒ 1͞3 , where s s is the rms bunch length and R is the bending radius.In Fig. 3 of Ref. [9], the peak intensity of the wake falls off in an approximately exponential fashion with a decay constant within 10% of the overtaking length.There is some reason to think that the overtaking length is applicable to the exit of the dipole, because it is simply a geometric factor related to the bunch length and curvature of the arc in the dipole.That is, in addition to giving the length over which CSR overtakes the head of the bunch, it gives the length over which changes in the conditions of the tail (e.g., whether the tail is inside the dipole) catch up with the head.
A second, more rigorous method provided in ELEGANT for modifying the intensity of the wake is based on Eqs. ( 53) and (54) in Ref. [3], which give the normalized energy gain for a particle inside a uniform rectangular bunch as a function of the normalized position ŝ sg 3 ͞R within the bunch and the normalized distance x xg͞R from the exit of the dipole.(Here s is the position in the bunch relative to the tail and R is the bending radius.)ELEGANT uses these equations solely for the purpose of determining the dependence of the energy gain on distance from the exit of the dipole.To do so, it averages the function over a series of values of the position s within the bunch.It was found that averaging over 20 values was more than sufficient to produce a converged result.The values of s are chosen to span the central 68% of the bunch (corresponding to 6s s for a Gaussian bunch).
The resulting table of values is normalized to a peak value of unity.When 20 values of s are averaged, the peak value is found to occur at x 0, i.e., at the exit of the dipole, which is consistent with the physical concept used here; namely, that the CSR wake in the drift is simply the terminal CSR wake in the dipole multiplied by some decay factor.To avoid having to make a very large table, the values of x for the table are spaced exponentially, which is suitable given the rapid decay of the function.Although a table of as few as ten values would probably be adequate, 100 values are typically used to be conservative.The code adjusts this value if necessary so that the initial point spacing does not exceed the user's specification of the interval for CSR kicks.The default initial point spacing in the table is 0.01% of the maximum value of x, given by Saldin et al. as s͑͞1 2 b͒, where s is the total bunch length and b y͞c.
In the next section, a comparison of the effect of choosing different models for the CSR in the drifts is made for the APS bunch compressor.

IV. APPLICATION TO THE APS BUNCH COMPRESSOR
ELEGANT was used to design and explore the performance of the APS bunch compressor [10], which features a four-magnet chicane with movable dipoles.The system starts with a photoinjector consisting of an S-band photocathode gun followed by an S-band capture section that accelerates the beam to 40 MeV.Following this, the beam is accelerated to (typically) 140 MeV by a set of four S-band accelerating waveguides (the "precompressor" linac), after which it enters the bunch compressor chicane.
The program PARMELA [11] was used to simulate the photoinjector [12], from which point ELEGANT was used.Figures 6-8 show the shape of the CSR wake at several points in the fourth chicane dipole.One sees that it is dramatically different from what would be expected from a Gaussian bunch.This is due to the appearance of a highcurrent spike at the head of the bunch.In this example, a 0.4 nC bunch is compressed to yield an average current of about 600 A, which is typical of what we plan to run at the APS facility in the near future.CSR is included in drift spaces using the second model discussed in the previous section.(This is true of all simulations unless otherwise noted.) To understand the appearance of the spikes at the head of the bunch, note that the predicted longitudinal distribution from the photoinjector shows a nonlinear dependence of momentum on position in the bunch, which manifests itself as one or more spikes in the momentum distribution entering the precompressor accelerator.Because of the relatively large energy chirp imparted to the beam by the precompressor accelerating section and the nearly Gaussian temporal distribution, the bunch at the entrance to the chicane has nearly Gaussian distributions in time and momentum.However, the "rotation" of the longitudinal phase space that occurs during bunch compression reveals the original nonlinear time-momentum relationship as spikes in the temporal distribution.Other contributing factors are  wakefields in the accelerator and nonlinearities in the beam transport through the compressor.Because of such effects, ELEGANT's use of the actual beam distribution (as opposed to a Gaussian approximation) in computing the CSR wake is expected to be important.
The appearance of the CSR wakes in these figures can be understood qualitatively by referring to Eqs. ( 1)-( 3).Referring to Eq. ( 2), one sees that T 1 ͑s͒ depends strongly on the derivative of the longitudinal density distribution just behind s since the 1͑͞s 2 z͒ 1͞3 term in the integrand gives emphasis to the region near s.Hence, when there is a spike in the longitudinal distribution, a significant component of the CSR wake will resemble the derivative of that spike, which is what is seen in the figures.In addition, referring to Eq. ( 3), one sees that T 2 may also be similar to a derivative of the longitudinal density distribution, given that it has the form l͑z͒ 2 l͑z 2 Dz͒.In the case of T 2 , however, the particles near s are subjected to a wake that resembles the derivative of l at s 2 2.5 ‫ء‬ s L , provided s L is small compared to the bunch length.In the APS compressor, this generally holds in the first quarter of the dipole.
It is interesting to look at the predicted final momentum distribution for this compressor as a function of the degree of compression, which is varied by changing the phase of the precompressor accelerator (in order to change the momentum chirp).The accelerator gradient is varied to maintain a fixed average beam momentum at the entrance to the chicane.Results from the simulation are shown in Fig. 9.
One sees that the momentum distribution has a spike at the upper end for f 90 ± (beam on crest in the precompressor linac); this is present in the results of the photoinjector simulation, as discussed above.For full compression, which occurs at f 66 ± , the momentum distribution is dramatically affected by CSR.Less dramatic effects are seen for partial compression and in the overcompression region (i.e., f , 66 ± ).These results are similar to those seen from other codes and in experiments at APS and elsewhere [13,14].Detailed inspection of the results from the simulations shows that the prominent spike for f 66 ± in Fig. 9 is generated primarily in the fourth dipole and following drifts, with some effect from the third dipole and following drifts.This momentum spike results from a large temporal spike at the head of the bunch, similar to the one seen in Fig. 8.The spike in the momentum distribution for f 71 ± is generated by the fourth dipole, whereas the distortion for f 56 ± is generated by the second dipole.For this latter case, the energy chirp is approximately double the energy chirp required for full compression, so that the bunch is already fully compressed at the exit of the second dipole and begins to lengthen in the third dipole.
A qualitative understanding of the appearance of the CSR-driven spikes in the momentum distribution can be reached by reviewing the remarks above about the dependence of the CSR wake on the derivative of the density distribution.Assume that the compression process creates a spike in the temporal distribution inside a given dipole.As a result of CSR, particles on the leading edge of the spike will be accelerated, while those on the trailing edge will be decelerated.For compression the leading particles must have lower energy than the trailing particles.Hence, the lower-energy particles near the temporal spike are accelerated while the higher-energy particles are decelerated, which produces a concentration, or spike, in the momentum distribution.Multiple momentum spikes may be produced by the action of CSR in several dipoles, by the occurrence of multiple temporal spikes as the bunch folds over itself, and by the complexity of the CSR wake as exhibited in the T 1 and T 2 terms.
It is worthwhile to compare the results of using the two models discussed above for CSR in the drifts.Figures 10  and 11 show, respectively, the momentum spread and horizontal emittance as a function of precompressor phase for three different situations: (i) no CSR in drift spaces, (ii) CSR in drift spaces using the overtaking length as the attenuation length, determined by analysis of the bunch longitudinal distribution, and (iii) CSR in drift spaces using Eqs.( 53) and (54) from Saldin, Schneidmiller, and Yurkov [3].
The contribution to the rms momentum spread from CSR in the drift spaces is small compared to the contribution from the chirp imposed for bunch compression, making it difficult to see any difference between the models.A clearer indication of the impact of CSR comes from the emittance.The contribution to the emittance due to CSR in the drift spaces is very large.The appearance of two peaks in the emittance curve is not unexpected, given the momentum spectra in Fig. 9.One peak occurs under the condition of full compression, near a phase of 66 ± ; in this case, a short bunch emerges from the third dipole and is corrupted by CSR primarily in the drift space between the third and fourth dipoles.The second peak occurs under the condition of maximum compression at the exit of the second dipole; in this situation, the compressed beam travels with a CSR wake through a drift section in a location of high dispersion.Given the approximate nature of the models, the agreement is good.In addition, as will be reported in upcoming publications, experiments with the APS bunch compressor show some agreement with these simulations.A clear peak in the measured emittance is observed as a function of phase.The measured peak emittance is about 50% higher than in simulation.We also see indications of a secondary emittance peak in the overcompression region, although it appears smaller than predicted.The differences may be due FIG. 11. (Color) Normalized rms horizontal emittance as a function of phase of the precompressor accelerator without CSR in drifts (1), with CSR in drifts using the overtaking length as the exponential decay length (2), and with CSR in drifts using Saldin's Eqs. ( 53) and (54) (3).
to the lack of detailed knowledge of the incoming temporal and momentum distribution, both of which clearly play an important role.
In summary, the corruption of the emittance is a result of intense CSR that results when spikes occur in the temporal distribution inside a dipole.Since the spikes in the temporal distribution originate in the nonlinear relationship between momentum and time in the beam at the input to the compressor, the appearance of prominent spikes in the momentum distribution after compression has its origin in the details of the distribution from the photoinjector and nonlinearities in the beam transport.This underscores the importance of integrated simulations of the photoinjector, beam transport, acceleration, and compression.With the addition of CSR computations using the method discussed in this paper, ELEGANT is a valuable component of such integrated simulations.

2 PRST-AB 4 SIMPLE
FIG. 1. (Color) Comparison of ELEGANT simulation and numerical integration, showing the change in g in successive pieces of a dipole.Successive frames are separated by an angle of 20 mrad.

FIG. 6 .
FIG.6.(Color) Bunch distribution and CSR wakes near the beginning of the fourth dipole.In this and following figures, Ds is the total length of the bunch.

FIG. 7 .
FIG. 7. (Color) Bunch distribution and CSR wakes at the center of the fourth dipole.

FIG. 8 .
FIG. 8. (Color) Bunch distribution and CSR wakes at the end of the fourth dipole.
second test made use of the steady-state CSR mode in ELEGANT, allowing comparison with analytical calculations for Gaussian bunches.From results given FIG. 2. (Color) Mean fractional momentum offset vs bending magnet radius.by Saldin et al., one can deduce the mean energy offset and the standard deviation of the energy offset.Letting d ͑p 2 p o ͒͞p o , where p o is the central momentum, one obtains