Extended 1D Method for Coherent Synchrotron Radiation including Shielding

Coherent Synchrotron Radiation can severely limit the performance of accelerators designed for high brightness and short bunch length. Examples include light sources based on ERLs or FELs, and bunch compressors for linear colliders. In order to better simulate Coherent Synchrotron Radiation, the established 1-dimensional formalism is extended to work at lower energies, at shorter bunch lengths, and for an arbitrary configuration of multiple bends. Wide vacuum chambers are simulated by means of vertical image charges. This formalism has been implemented in the general beam dynamics code"Bmad"and its results are here compared to analytical approximations, to numerical solutions of the Maxwell equations, and to the simulation code"elegant".


I. INTRODUCTION
It is envisioned that future accelerators will call for shorter beams of higher intensity. A possible limiting factor in these efforts is an increase in energy spread and transverse emittance, as well as a micro-bunching instability, due to Coherent Synchrotron Radiation (CSR).
The first CSR calculations were performed by Schwinger in 1945. Using a Green's function method, he arrives at the power spectrum of a single charge bending in free space as well as between infinite conducting plates, and thereby computes the coherent power radiated by a collection of charges [1]. Warnock extends this work to include the longitudinal impedance on a bunched beam [2]. Many papers covering the history and importance of CSR forces can be found in [3].
This paper uses an approach to calculate the CSR wake-field originating with Saldin et al. [4] and generalized by Sagan [5]. Here we calculate the CSR force between two charges traveling on the same curve, and integrate over a longitudinal bunch distribution to give a longitudinal wake-field. Transverse particle coordinates and transverse force components are neglected. The formalism developed here is generalized to include arbitrary lattice configurations of bends and drifts, including, for example, radiation from one bend entering another, and bend radiation extending into drift regions.
Simulating CSR effects is the subject of a number of codes. The method here is implemented in the particle tracking code Bmad [6]. Our simulation results are compared with two of the codes described in [7] and with approximate analytic formulas.

II. TWO PARTICLE INTERACTION
In order to compare the force acting on one particle from the radiation that is emitted by another, the analysis starts by considering two particles of charge e following the same trajectory as shown in Fig. 1. The Liénard-Wiechert formula [8] gives the electric field E(P) at the position of the kicked particle at point P and time t due FIG. 1: A particle at point P ′ kicks a particle at point P.
to the source particle at point P ′ and retarded time t ′ E(P) = e 4πǫ 0 1 γ 2 (L − Lβn ′ ) + 1 c 2 L × [(L − Lβn ′ ) × a ′ ] (L − L · βn ′ ) 3 (1) It will be assumed that both particles have the same constant speed β = v/c, n ′ and n are the unit velocity vectors for the source and kicked particles respectively, and L = P−P ′ is the vector from the source point to the kick point. The retarded time t ′ is related to t via t−t ′ = L/c. At time t, the source particle has a longitudinal position z ′ with respect to the bunch center and the longitudinal position of the kicked particle is z. The distance ζ ≡ z−z ′ between the particles at constant time can be computed via the equation where L s is the path length from P ′ to P. Generally, the relativistic approximation β = 1 will be made. However, some terms in 1 − β ≃ 1/2 γ 2 will need to be retained. The first term on the right hand side of Eq. (1) has a 1/ζ 2 singularity at small distances. Following Saldin et al. [4], this singularity is dealt with by dividing the electric field into two parts. The space charge component E SC , which contains the singularity, is the field that would result if the particles where moving without acceleration along a straight line. The CSR term, E CSR , is what is left after subtracting off the space charge term The rate K ≡ dE/ds at which the kicked particle is changing energy due to the field of the source particle is Following Saldin et al. [4], the transverse extent of the beam will be ignored in the calculation of K CSR . However, the inclusion of the finite beam size will be needed to remove the singularity in the calculation of K SC as discussed in Section IV.

III. CSR CALCULATION
The source point P ′ and the kick point P will, in general, not be within the same lattice element. Because the transverse extent of the beam is being ignored, all elements will be considered to be either bends or drifts.
In Fig. 1, R is the bending radius and g = 1/R is the bending strength of the element that contains the source point P ′ . The magnitude of the acceleration is a ′ ≃ c 2 /R. This element ends at point O. The angle between P ′ and O is φ, and d = R φ is the path length between P ′ and O.
Between point O and the kick point P, d i is the path length within the i th element, i = 1, . . . , N , where N is the number of elements in this region. For the last element, d N is the distance from the start of the element to point P. For the i th element, φ i is the bend angle, R i is the bend radius, and g i = 1/R i is the bend strength. For a drift φ i = g i = 0.
In Fig. 1, (v, w) are the coordinates of point P with respect to point O with the v-axis parallel to the orbit's longitudinal s-axis at point O and the w-axis pointing upwards towards the inside of the element containing the point P ′ .
With this notation, the difference in v and w from the beginning of an element to the end is where ψ i is the orientation angle at the entrance end of the element The above formulas are able to handle negative bends (beam rotating clockwise). For a negative bend R i , g i and φ i are negative while d i = R i φ i is always positive.
With the assumption that all bend angles are small, v and w can be approximated by where and the small angles have been retained to second order. The angle θ of the vector n with respect to the v-axis is θ = N i=1 g i d i . In terms of v and w, the components of the vector L are Again angles are retained to second order. The path length is simply This, With Eq. (2), gives where terms to second order in combinations of angles and 1/γ are retained. Substituting these expressions into Eq. (1), and defining the individual terms in Eq. (1) read as Putting these together yields While we have used S.I. units, the classical radius r c and the mass m has been used to make the formula independent of the unit system. From Eq. (14), K • CSR , which is K CSR restricted to the special case where points P and P ′ are within the same bend, reduces to Eq. (32) of Saldin et al. [4], whereφ ≡ γ φ, α = Rφ 2 /2, κ =φ, and τ = Rφ. This equation is valid forφ > 0; forφ < 0, K CSR is, to a very good approximation, zero.
In the limit of small ζ, K • CSR has a limiting value of At large values of ζ, ζ is cubic in φ so that φ ≃ (24 ζ/R) 1/3 . With this, Eq. (15) becomes which corresponds to Eq.  The fact that K CSR is highly peaked in amplitude near ζ = 0 can be problematic for simulations at ultrarelativistic energies because the characteristic longitudinal distance between particles or mesh points needs to be less than R/γ 3 . One way of dealing with the peaked nature of K CSR is to first consider the kick from a line of particles of density λ(z) and then to integrate by parts where I • CSR , which is I CSR for P and P ′ in the same bend, is plotted in Fig. 3. The peaked nature of K CSR has been smoothed over at the cost of having to deal with a derivative of λ. For ζ ≫ R/γ 3 , the approximation of Eq. (17) can be used to calculate an explicit ultrarelativistic equation for I • CSR as in [9], Equation (21) is also plotted in Fig. 3. While, in general, it is helpful to have explicit formulas, for the purposes of evaluation within a simulation program this is not needed. The alternative is to use an exact implicit solution. Because Eq. (11) and Eq. (14) are rational functions, Eq. (20) can be integrated. The last term that compensates K SC can be integrated to r c mc 2 /(γ 2 ζ). The other terms can be written with as With ∂τ /∂d = γ, ∂α/∂d = γgτ , and ∂κ/∂d = γg, one can further simplify to This can be integrated over ζ to yield It can be shown that, while quantities like d and g are discontinuous across element boundaries, τ , α, and κ are continuous and hence I CSR is a continuous function as it should be. Equation (25) is the main result of this paper. Using Eq. (25), the integration of Eq. (19) in a simulation program can be done via interpolation of Eq. (11). Equation (25) has several advantages over equations like Eq. (21). It is applicable at lower values of γ 3 ζ, that is, at lower energies and/or smaller length scales. Additionally, Eq. (25) has no singularity at small ζ, and it can be used to handle any combination of elements between the source and kick points.

IV. SPACE CHARGE CALCULATION
The singularity at small ζ in the space charge term E SC is removed by considering the finite transverse beam size. This term is equivalent to the problem of calculating the field given a static distribution of charges. It will be assumed that at any longitudinal position the transverse profile of the beam is Gaussian. Thus, a longitudinal slice of the beam will produce an energy change for a particle at longitudinal z and transverse offset (x, y) from the slice center of A heuristic solution for Eq. (26) in the region of interest (x 3 σ x and y 3 σ y ) is where sign(z) is 1 for positive and −1 for negative z, and ζ = z − z ′ . Equation (28) is exact in the limit z = 0 and z → ∞, and is an excellent approximation in the region in between. This is illustrated in Fig. 4, which shows K SC as a function of z as computed from an integration of Eq. (26) and from the approximate Eq. (28). The particular parameters chosen for the computation are given in the figure. Two cases were considered. One where the kicked particle is on-axis, and the other where the kicked particle is displaced by σ x off-axis. As can be seen, Eq. (28) gives an excellent approximation to the longitudinal space charge kick.
At high energies, the CSR energy kick is independent of the beam energy as indicated by Eq. (21). On the other hand, the factors of γ in the denominator of Eq. (28) insures the K SC will be a decreasing function of γ. Past some point in an accelerator, the increasing γ will make the effect of K SC small compared to the effect of K CSR . To estimate that point, consider the maximum energy kick in a Gaussian bunch of N e particles with The space charge kick is maximum at x, y = 0 and using Eq. (28) gives The dominant term in the denominator in the integrand is either σ x σ y for small ζ, or γ 2 ζ 2 for large ζ.
As an approximation, the middle term in the denominator will therefore be ignored. The approximation sinh zζ/σ 2 z ≈ zζ/σ 2 z will also be made. This approximation is justified since, as shown below, the region of interest for z is around z ≈ σ z and for realistic beam parameters, σ x σ y ≪ γ 2 σ 2 z and therefore the only significant contribution to the integral will come in the region ζ < σ z . With these approximations, Eq. (30) becomes The maximum CSR kick at large energies is, from Eq. (14), The ratio is The condition for K SC begin small can be written as The condition for this can be well approximated by noting that the quantity in square brackets on the LHS of Eq. (35) is slowly varying and never extremely large. Thus the required condition is where M is a number of order unity. In this case, σ z must be interpreted as the characteristic longitudinal distance over which the bunch density is changing.

V. CSR IN BMAD
The above algorithm for simulating CSR and Longitudinal Space Charge (SC) has been implemented as part of the Bmad [6] subroutine library for relativistic chargedparticle simulations. Bmad simulates a beam as a set of particles. The beam is tracked through a lattice element FIG. 5: Bmad implementation of the CSR algorithm. The beam of particles is divided up into a number of bins. The contribution of a particle to a bin's total charge is determined by the overlap of the particle's triangular charge distribution and the bin.
by dividing the element into a number of slices. Tracking through a slice involves first propagating the particles independently from each other and then applying the CSR and SC energy kicks. To calculate the energy kick, the beam is divided longitudinally into N b bins as shown in Fig. 5. For computing the charge in each bin, each beam particle is considered to have a triangular charge distribution. The overlap of the triangular charge distribution with a bin determines that particle's contribution to the total charge in that bin. The width of the particle's triangular charge distribution and the number of bins are set by the user. The bin width is dynamically adjusted at each time step so that the bins will span the bunch length. Increasing the particle width smooths the distribution at the cost of resolution.
The charge density λ i at the center of the i th bin is taken to be λ i = ρ i /∆z b where ρ i is the total charge within the bin and ∆z b is the bin width. The charge density is assumed to vary linearly in between the bin centers. The CSR energy kick for a particle at the center of the j th bin after traveling a distance ds slice according to Eq. (19) is then where Evaluation of I CSR,j involves inversion of Eq. (11) to obtain d. Because z is a monotonic function of d, Newton's method [10] is used to find numbers d 1 and d 2 which bracket the root and then Ridders' Method [10] is used to quickly find d.
In deriving Eq. (37), the approximation has been used. Generally this is an excellent approximation, except when j = 0 and ∆z b ≫ R/γ 3 , as shown in Fig. 3. Here, however, the integral can be done exactly assuming that the source and kick points lie within the same element (40) Once the energy kick at the centers of the bins is calculated, the energy kick applied to a particle is calculated via interpolation assuming a linear variation of the kick between bin centers.
In calculating the energy kick, The computational time for calculating the charge in the bins charge scales as N p , the number of particles in the simulation. The computational time for calculating the energy kick at the bin centers scales as N 2 b , and the time for calculating the energy kick of the particles scales as N p .

A. Chamber Walls
The simulation incorporates the shielding of the top and bottom chamber walls by using image currents. Appendix A explains why neglecting the width of a chamber is a good approximation when it's width is larger than its hight. Here, because the image current is well separated from the actual beam, there are no singularities to deal with, K SC does not have to be subtracted, and a straight forward integration is done using Eq. (18), and Eq. (14).
where q i = λ i ∆z b is the charge in a bin, h is the chamber height, and k indexes the image currents at vertical displacement y = ±k h. The number of image charges N i needs to be chosen large enough so that the neglected image currents do not have a significant effect on the simulation results. Because the relevant angles are not small, the image charge kick K must be calculated without the small angle approximation, as in Eq. (14).

VI. AGOH AND YOKOYA CSR CALCULATION
Agoh and Yokoya (A&Y) have developed a code to calculate CSR wake fields by directly integrating Maxwell's equations on a mesh representing a rectangular beam chamber [11]. The approach depends on the paraxial approximation, a rigid Gaussian bunch density, small chamber dimensions relative to the bending radius, and ultra-relativistic particles. Using a co-moving coordinate system (x, y, s) in Fourier space, they are able to reduce the problem to a tractable two-dimensional differential equation, where E ⊥ is a complex 2-dimensional vector related to the perpendicular electric field, ρ 0 is the charge density, k is the wave number, and R is the magnet bending radius. It is solved using a finite-differencing method.

VII. CSR IN ELEGANT
The particle tracking code elegant (version 17.2.2) uses Eq. (43) to compute CSR kicks without shielding by a vacuum chamber [12]. The charge distribution λ(z) and its derivative dλ/dz are calculated by binning the macroparticles and then employing a smoothing filter.

A. Transient Effects at Magnet Edges
Using retarded fields, Saldin et al. [4] derive, in the ultra-relativistic limit, a formula for the wake-field due to a bunch entering from a drift region into a bend where φ is the angle traveled into the magnet by the bunch center. Eq. (43) reduces to the free space steadystate wake-field of Eq. (21) in the limit s L → ∞. As worked out by Emma and Stupakov [9], synchrotron radiation will continue to propagate and affect the bunch beyond the end of a bending magnet. For a finite magnet of length L m , an ultra-relativistic bunch at a distance x from the end of this magnet experiences the free space exit wake-field

B. Steady State CSR in a Bend
CSR effects in a vacuum chamber have been computed by the Green's function of grounded parallel plates [1,2]. These formulas are difficult to compute numerically, due to the presence of high order Bessel functions, so we will use an excellent approximation developed by A&Y [11]. The impedance for the steady-state in a dipole with horizontal plates separated by a distance h is where Z 0 = cµ 0 is the free space impedance, k is the wave number, and Ai and Bi are Airy functions. The parallel plate wake-field due to a bunch with longitudinal density λ(z) is obtained by Fourier transform:

IX. COMPARISON BETWEEN BMAD, THE AGOH AND YOKOYA, AND ELEGANT
In order to validate our method, we compare simulations from Bmad to those using the A&Y code and elegant. For ease of reading, all magnet and bunch parameters used are enumerated in Tab. I, and will be referred to by a letter.
Note that all electric fields in the graphs are normalized by which approximately describes the maximum amplitude of the CSR-Wake dE/ds. All simulations use a bunch charge of 1nC, and an energy of 1 GeV, unless otherwise noted. We used Bmad with the following parameters: number of bins N b = 800, number of macro-particles N p = 4 · 10 5 , number of image charges N i = 32, and tracking step ds step = 1mm. The triangular bin width, as in Fig. 5, is 32 bins. Figure 6 shows that this choice of N b and N p is reasonable by varying N p and N b , using parameter set A. While reducing the number of particles leads to a visibly less smooth wake field. The number of bins has been increased, until the wakefield starts becoming less smooth. The other parameters were similarly varied to find the applied parameters. Figure 7 shows the steady-state CSR kick in a bend as a function of z for various values of the chamber height. The parameters used correspond to set A of Table I and are the same as used for Fig. 1 of A&Y [11]. Figure 7 shows excellent agreement among the CSRmesh code, Bmad, and the CSR-Wake formula Eq. (49). Note that Bmad computes the wake by tracking a bunch, and therefore shows the result only in the range of the length of the bunch. This is not problematic, because the wake only influences particles in the region of the bunch.

B. Transient Case
Parameter sets B, C, and D are used to explore the transient case where some of the kick is generated from particles in the drift region before the bend. Set B corresponds to values used in Fig. 3 of Agoh and Yokoya [11], which has a chamber of relatively large cross section. Figure 8 (left) shows agreement between Bmad, the A&Y code, and the CSR-Wake formula Eq. (43). However, for parameter set C in Fig. 8 (right), where the chamber is smaller as in Fig. 4 of [11], a discrepancy appears in the front of the bunch. This discrepancy is noted in [11] and explained to be due to the finite chamber size.
We use parameter set D to analyze whether the chamber width is responsible for this difference, and one sees in Fig. 9 that the effect remains and thus appears to be due to the reduced chamber height. It also is larger for longer magnets. Varying the number of Fourier coefficients used in the A&Y calculation in Fig. 10 does not change this result either, verifying that we chose a reasonable number of Fourier coefficients for solving the Maxwell equations numerically. Reversing the set D height and width, as in Fig. 11, somewhat reduces the discrepancy. Thus, again indicating that the effect is apparently due to the reduced chamber height. Practically speaking, the difference in these wake-fields only becomes appreciable far in front of the bunch, where there are few particles to affect. To see how strong the deviation becomes for especially small chamber heights, we repeat parameter sets B and C with h = 2cm in Fig. 12, and the disagreement between Bmad and the A&Y code is again not very large but significant.

C. Realistic Magnets
To evaluate how significant the differences are in realistic magnets, we use parameter sets E, F, and G which correspond to the JLab TH2, CESR Analyzer magnet, and Cornell ERL CE magnets, respectively. Wake-fields for the first two are plotted in Fig. 13, showing good agreement between the A&Y code and Bmad. The free space steady-state wake-field Eq. (43) and results from elegant are plotted for reference.
The principal detrimental effects of the CSR-Wake in an accelerator are energy loss and increase in energy spread of a bunch. The transverse bunch distribution can also be damaged and is mostly influenced by increases of the energy spread which, through dispersive orbits, couples to transverse motion. To visualize CSR- driven energy loss and energy spread in a single magnet, we plot in Fig. 14 the average and RMS electric field across the bunch distribution as a function of distance into the magnet. This is done using parameter set G, where we examine how well one can ignore the chamber width, as done in the Bmad calculation. One sees that, in this shielded case, the chamber width w in the A&Y code begins to change the wake field when it is comparable or less than 4cm, the height of the chamber. This effect is heuristically explained in appendix A. It again shows that ignoring the chamber width, as in Bmad is a reasonable approximation when the chamber is wider than high.
It is apparent that shielding by a vacuum chamber reduces the power emitted by CSR very effectively. It does not reduce the energy spread nearly as much, and is therefore not as effective for preserving bunch properties as one might have concluded from the reduced radiation power. Interestingly, the 4cm wide chamber even produces larger RMS E s than wider chambers.

D. Exit Wake
The method in this paper can correctly account for the CSR-Wake in a drift section following a magnet. Using parameter set A, Fig. 15 shows this wake as a function of the distance d from the end of the magnet in free space and with shielding. The free space case shows excellent agreement between Bmad and the CSR-Wake formula Eq. (44).
The other codes do not allow simulations with shielding in this regime. We therefore compare Bmad with a numerical solution using retarded fields. In order to avoid the complication of computing retarded times and positions of the bunch distribution and its image charges, we do not use Liénard-Wiechert fields, but rather Jefimenko's equations [8]. In general, the electric field due to a 1-dimensional charge density ρ(s, t) and current density J(s, t) at a position s, time t, and height h is where L is the vector from position s ′ to position s at height h, L is its magnitude, and the brackets [ ] are evaluated at the retarded time t ′ = t − L/c. As compared with integrating over retarded fields at s ′ , using Jefimenko's equations has the advantage that one never has to solve for the retarded time in equations of the type t ′ = t − |r − r(t ′ )|/s. The total electric field due to alternating image charges is then Applying this to the geometry of a bend followed by

E. Coherent Energy Loss
Our final test of Bmad compares the total coherent energy lost for various particle energies and chamber heights to the integration of the power spectrum. In general, for N particles traveling on the same curve at different phases φ n , the N particle power spectrum is where dP (1) /dω is the single particle power spectrum, and λ(z) is the longitudinal particle distribution. The first term in Eq. (55) is the incoherent power spectrum, while the second is the coherent power spectrum. In the presence of conducting parallel plates, the single particle power spectrum is given in Eq. (47) of [1]. Using this, Eq. (55) can be integrated numerically, and in Fig. 16  resulting coherent part shows excellent agreement with Bmad for energies down to 5MeV and chamber heights down to 2mm. At smaller heights, the number of image layers used in the simulation (N i = 64 here) is not sufficient to correctly model the CSR.
In absence of shielding plates, Eq. (55) can be integrated exactly for a Gaussian distribution using the wellknown free space single particle power spectrum. With a standard deviation σ z , the total power lost by N particles is where is the power lost by a single particle, and This result agrees well with Bmad in Fig. 16. The function T (a) can be expanded asymptotically, giving the leading order coherent energy change in a length L as This is consistent in the scaling and magnitude of E 0 in Eq. (51).

X. CONCLUSION
A general implicit formula for the longitudinal kick due to the coherent synchrotron radiation has been developed for particles on a common orbit. This formalism will handle any geometry of bends and drifts. For simulations, this formula is to be preferred over the explicit ultrarelativistic formula because the implicit formula does not have a singularity at z = 0 and is applicable at lower particle energies and smaller length scales.
Additionally, a heuristic formula for the longitudinal space charge kick has been presented which takes into account transverse displacements of the kicked particles. This formalism has been implemented in Bmad. We show that the longitudinal wake field compares well with the A&Y code of A&Y [11] and the CSR-Wake formula of Warnock [2] for the steady state, with and without CSR shielding by parallel plates. In the transient case where the A&Y code often does not follow the CSR-Wake formula of [13] exactly, Bmad does agree well with that formula.

XI. ACKNOWLEDGMENT
The authors would like to acknowledge help from Tsukasa Miyajima and Ivan V. Bazarov. And we thank Michael Borland for useful discussions. This work has been supported by NSF cooperative agreement PHY-0202078.

APPENDIX A: HEURISTIC SHIELDING ARGUMENT
FIG. 17: Numerical results using the A&Y code with varying dimensions of a rectangular chamber. A 1m, 7.5 degree bend with a 2ps long bunch of 0.82nC charge was used. Contours represent the energy change induced due to the CSR wake field in 2000eV increments. The shielding effect thus not only due to the vertical dimension, but primarily due to the smaller of the two beam-pipe dimensions.
It initially might seem surprising that vacuum chambers of many centimeter width can have a shielding effect on synchrotron radiation with much smaller wavelength. We therefore add a heuristic explanation here. Starting with Schwinger [1], shielding by a vacuum chamber has often be considered by studying infinite horizontal, conducting plates. But the following heuristic argument indicates why it is both the vertical and the horizontal boundary that determines shielding of CSR. Numerically solving Maxwell's equation in the vacuum chamber supports this argument as shown in Fig. 17.
A highly relativistic particle emits synchrotron radiation within a narrow cone. For the radiation's component of wavelength λ, the opening angle of this cone is approximately ∆θ = ( λ R ) 1 3 in the horizontal and vertical. The opening angle in the vertical determines the vertical width of the radiation load on the vacuum chamber wall, and the horizontal divergence can be observed by shining through a pinhole.
The radiation field builds up within a radiation buildup time of ∆t = R c ( λ R ) 1 3 . During this time, the radiation fields produced by the electron coherently add up to form the full radiation power. If the radiation does not interfere with an obstacle, for example the vacuum pipe, within this time, the electron looses as much energy as it would without any vacuum pipe.
The width and hight of the radiation cone that builds up during the radiation buildup time is therefore given by w r ≈ h r ≈ c∆t∆θ = R( λ R ) 2 3 . Vacuum chambers that have smaller dimensions interfere with the radiation process and shield the part of radiation for which Wavelengths are therefore shielded when they are above a length that is much smaller than the chamber dimensions.
While we have used a very approximate heuristic argument, Fig. 17 computed by the A&Y code, indeed shows that both dimensions can lead to shielding, and that to first approximation only the smaller of the two dimensions is relevant.