Calculation of coherent synchrotron radiation using mesh

054403-1 We develop a new method to simulate coherent synchrotron radiation numerically. It is based on the mesh calculation of the electromagnetic field in the frequency domain.We make an approximation in the Maxwell equation which allows a mesh size much larger than the relevant wavelength so that the computing time is tolerable. Using the equation, we can perform a mesh calculation of coherent synchrotron radiation in transient states with shielding effects by the vacuum chamber. The simulation results obtained by this method are compared with analytic solutions. Though, for the comparison with theories, we adopt simplifications such as longitudinal Gaussian distribution, zero-width transverse distribution, horizontal uniform bend, and a vacuum chamber with rectangular cross section, the method is applicable to general cases.


I. INTRODUCTION
Electrons traveling along curved trajectories emit synchrotron radiation.The component of radiation whose wavelength is longer than the electron bunch length is emitted coherently.This component is called coherent synchrotron radiation (CSR).As is well known, its transverse and longitudinal forces cause serious influences on the particles such as emittance growth and increase of energy spread [1][2][3].Recent accelerators tend to adopt short bunches for various demands.For example, freeelectron lasers require a high peak current with a given bunch charge, energy recovery linacs often require a short duration of radiation, and B factories and linear colliders demand a short bunch for higher luminosity.In the design of such machines, CSR is an important matter of concern.
Until now, a number of analytic and numerical studies on CSR have been made.Analytic solutions have been found under some conditions, namely, steady states in vacuum [4], steady states between parallel plates [5], and transient states in vacuum [6].These analytic formulas describe the essential properties of CSR.However, if transient effects are substantial, and suppression by vacuum chamber is not negligible, then these analytic formulas cannot be applied.Furthermore, the transverse beam size is neglected and the bunch shape has to be rigid in these formalisms.For these reasons, one has to rely on numerical approaches in the practical research.
Several numerical analyses have been studied on CSR effects.Borland [7] has developed a semianalytic tracking code where CSR is given by the analytic solution with a line charge distribution.It simulates the energy change by CSR.Emma et al. has developed a tracking code [8] in which CSR is also given by an analytic formula.Although they neglect transverse forces by CSR, emittance growth via energy change is considered.It works fast because of the simple 1D model without shielding.Hajima has estimated emittance growth due to shielded CSR [3].In his study CSR is given by the analytic formula with infinite parallel plates.Considering more realistic cases, Li [9] has developed a tracking code in which the bunch has 2D charge distribution (no vertical size).A 3D tracking code TRAFIC 4 [10] has been developed by Dohlus, Kabel, and Limberg.These two codes are based on integration of the retarded potential, considering the shielding effects by infinite parallel plates.Assuming a perfectly conductive metal, the parallel plates can be substituted with image charges.In the case of strong shielding, the simulation takes a long time due to the enormous number of image charges.A different approach has recently been developed by Stupakov and Kotelnikov [11].It calculates the eigenmode of synchrotron radiation in a toroidal chamber.The different point is that it solves approximate Maxwell equations in the frequency domain.Considering steady states, they have solved the eigenvalue problem with a troidal boundary.The chamber cross section is rectangular or round.This method cannot take into account the transient effects.
In this paper, we present a numerical method based on the mesh calculation of the electromagnetic field in the frequency domain.We assume ultrarelativistic electron bunches as source charges 1, where is the Lorentz factor.We assume, as in [11], that the size of the chamber cross section a is much smaller than the radius of curvature, namely, a= p 1.In general, in order to carry out the correct simulations, the mesh size should be much smaller than the wavelength of the field.In our formalism, however, the mesh size is allowed to be larger than the radiation wavelength, owing to the neglect of the radiation propagating in the opposite direction to the beam.
We begin with Maxwell equations in vacuum and Fourier transform them into a frequency domain.Then we approximate these equations on the assumptions as above and solve them numerically by a finite difference.A reference orbit which gives the coordinates, is arbitrary PHYSICAL REVIEW SPECIAL TOPICS -ACCELERATORS AND BEAMS, VOLUME 7, 054403 (2004) 054403-1 1098-4402=04=7(5)=054403(9)$22.50  2004 The American Physical Society 054403-1 in the vacuum chamber.The source bunch can be set to any 3D charge distribution around the reference orbit, though it should be sufficiently smooth to require few Fourier modes.The evolution term in the basic equation enables us to simulate CSR in transient states.The cross section of the vacuum chamber can be any shape and vary along the reference orbit, although we present in this paper a rectangular vacuum chamber with a uniform size.In this method we solve Maxwell equations by using meshes with boundary conditions.We present the basic idea and derive the equations in Sec.II.The algorithm in the numerical calculation is given in Sec.III.In Sec.IV we show the results obtained by using this method with analytic formulas, and then we conclude in Sec.V.

II. THEORY A. Assumptions and notation
To derive the basic equations, we adopt some approximations summarized below.
(a) The size of the chamber cross section a is much smaller than the bending radius , namely, a= p 1. (b) The bunch consists of ultrarelativistic electrons 1.Under this assumption, the longitudinal component of the electric field is much smaller than the transverse one.
(c) The radiation components propagating at large angles with respect to the beam are ignored (paraxial approximation).In particular, this assumption excludes a vacuum chamber having a projection from the wall, which would cause a wave propagating along the opposite direction.
(d) The bunch shape does not change.This assumption can be relaxed so as to include ''predictable'' changes (those estimated by a simple optics calculation).The dynamic change of the bunch shape due to the CSR itself cannot be included.
We employ a coordinate system x; y; s that is often used in accelerator physics, namely, s being the length along the reference orbit, and x and y are perpendicular to s.For simplicity we assume a planar orbit on the x; s plane with curvature radius s, though it is easy to extend to the case of the general 3D orbit.Thus, the unit vectors e x ; e y ; e s satisfy de x ds e s s ; de s ds ÿ e x s : (1) We mainly work in the frequency domain.The variables in the time domain have a symbol t, and others are variables in the frequency domain.Let k be the wave number.We define Fourier transformation of a function f as follows: where t ÿ s.This means that fk is a normal Fourier transform of f t with respect to t, but the factor e iks is taken off.We put c 1 (velocity of light) throughout the present paper.

B. Equations
In this subsection we reduce Maxwell equations to a simpler form for numerical calculation.We denote the electric and magnetic field by E, B, the charge density and current density by J 0 , J.
Maxwell equations in the accelerator coordinate are where g 1 x=s.We suppose that the current flows along the reference orbit J t J t s e s .Let us go to the frequency domain.Note that the field variables are functions not only of (x; y) but also of s.Fourier transformation of Eqs. ( 3) and (4) with higher order terms O 3 ignored gives Neglecting higher order terms O 2 , Eq. ( 5) can be expressed with the electric field: B s can be obtained from E x ; E y as follows: Thus, the longitudinal components can be expressed in terms of the transverse components.For these two expressions only the leading terms (B x ÿE y , B y E x ) are needed in ( 6) and (7).The Lorentz force in the transverse direction up to O 2 is given by the following equations: PRST-AB 7 CALCULATION OF COHERENT SYNCHROTRON . . .

(2004)
054403-2 054403-2 These are obtained by using the full expressions ( 6) and ( 7), but they will not be discussed in the present paper where the energy gain/loss is the main concern.Now we derive an equation which describes the evolution of the transverse fields.From Maxwell equations, we obtain We Fourier transform Eq. ( 12) into the frequency domain and neglect the higher order terms O 2 : where The last term in Eq. ( 16) might be large where the curvature changes rapidly, in particular, at the edge of the bending magnets.Integrating the term over the edge, we can roughly estimate the jump of E x approximately as We can neglect the last term in Eq. ( 16).
Similarly, the first terms in the square brackets in Eqs. ( 16) and ( 17) might be large.Integration of the term over the edge gives a rough estimation of the jump of E x , E y as Thus the first terms in Eq. ( 16) and (17) are negligible.
In the case of the hard edge magnet, the second term in the square brackets in Eqs. ( 16) and (17) can be the square of the function because E x;y itself can be step function like.The jump of E x at the edge is where l edge is the edge length.If the bunch length is comparable or shorter than the edge length of the magnet s 1=k & l edge , the jump E x is O 2 so that we can neglect the term.Thus the terms C x , C y are negligible.
Finally, we make the most important approximation.We defined the Fourier transform with the factor e iks taken off.If the wave propagates nearly parallel to the beam (paraxial approximation), the remaining s dependence of Ex; y; s will be weak.Thus, we neglect the terms of the second derivative with respect to s in Eqs. ( 13) and ( 14).This approximation means that we neglect the wave which has a phase ks t; that is, we ignore the radiation which propagates in the opposite direction to the beam.If the vacuum chamber has projections on the wall, we cannot use this approximation because of the reflection waves.As a result, all the terms on the right-hand sides (RHS) in Eqs. ( 13) and ( 14) are small.Thus finally we obtain where r ?@ x ; @ y and r 2 ?@ 2 x @ 2 y .The RHS of this equation is identical to the fundamental equation in [11] except that we included the source term and assumed 1.Our equation describes the evolution of the field whereas the equation in [11] describes the steady state only.
Let us see Eq. ( 21) from a different point of view.If we remove the source term from Eq. ( 21), it is identical to the equation of laser propagation in the nonuniform medium in terms of the Cartesian coordinates x; y; z: Here n is the refraction index of the medium.It is a function of the coordinate in case of a nonuniform medium.The laser beam is bent due to a nonuniform medium.Comparing Eq. ( 21) with Eq. ( 22), we find that the term 2k 2 x= in Eq. ( 21) plays a role of the refraction index nx 1 ÿ x= varying as x.This means that the radiation follows a bent path in vacuum.The reason is that we use curvilinear coordinates, so that the radiation propagating straight deviates from the curved axis.
We also find Eq.( 21) is a Schro ¨dinger-like equation.We can in fact obtain the Schro ¨dinger equation from the Klein-Gordon equation by using the same mathematical procedure, i.e., to take out the time dependence e ÿimt (m: rest mass) from the wave function and to ignore the remaining second derivative with respect to t.
Equation ( 21) can be solved as an initial value problem with respect to s under the boundary condition of the vanishing tangential electric field on the wall.The Maxwell equation is basically elliptic with respect to the spatial coordinate.Nevertheless, our equation can be PRST-AB 7 T. AGOH AND K. YOKOYA 054403 (2004) 054403-3 054403-3 solved step by step with s.The reason is that there is no backward wave in our case.The charge distribution J 0 x; y; k is arbitrary in our method.We discuss here, however, only the case where the transverse and longitudinal dependences are factorizable: where k is a line charge density in the frequency domain and 'x; y is the charge distribution in the transverse dimension.We adopt a Gaussian distribution for t .The line charge distribution in the frequency domain k is also Gaussian: where s is the rms bunch length.
Before giving an algorithm of solving Eq. ( 21) numerically in the next section, we discuss the possible singularity coming from the transverse distribution 'x; y.When 'x; y is smooth enough, one can apply the algorithm in the next section directly to Eq. ( 21).In some applications as in the present paper, however, one may not be interested in the transverse charge distribution and the transverse force.In such cases one can use a line charge distribution,

'x; y xy:
(26) (We assume the charge is located at the origin x y 0.) In this case E x (E y ) has a strong singularity like xy=x 2 y 2 near the origin.To relax the singularity coming from 'x; y, we separate the field into two parts where E r is the radiation field, and E b is the beam field which is defined as the electric field along the straight trajectory (with or without the beam-pipe boundary).

E b ? E b
x ; E b y satisfies the equation: Substitution of Eqs. ( 27) and (28) into Eq.( 21) eliminates the source term Similarly, Eq. ( 8) becomes x @E r y @y : (30) Note that E b s 0 for E b satisfying Eq. ( 28) for ultrarelativistic beams.Obviously, we have to use E r E b in imposing the boundary conditions.
In the case of the line charge (26), E b x and E b y are given without the beam-pipe boundary condition by (31) (32) where N is the number of electrons in the bunch and e is the elementary charge.This expression is enough for the present purpose.However, there is still a logarithmic singularity which gives rise to the so-called Talman force.It is absent in E s but must correctly be treated when the transverse force is needed.In such a case one possible way of relaxing the singularity of E r is to include the (known) logarithmic term in E b without imposing the condition (28).Another way to avoid the singularity is to introduce the finite transverse distribution instead of the line charge (26).One then need not introduce E b .Even in such cases, however, it is better to separate E b when the beam transverse dimension is small compared with the characteristic scale length introduced in the next section, because otherwise a smaller mesh size would be needed to resolve the small transverse distribution.For such purposes a distribution for which an analytic expression of E b ? is known is better (like a 2-dimensional Gaussian distribution).These prescriptions are necessary only when the transverse force is needed.In the present paper we adopt the line charge distribution.

A. Finite difference
We here assume, for simplicity, the chamber is rectangular with a constant size.We introduce a finite difference to Eqs. ( 29) and (30).x i ix, i 0 m, y j jy, j 0 n, s ' 's, ' 0 u, where i, j, and ' are the indices of x, y, and s, respectively.For the convenience of using the central difference, we define the field components at s s ' as E x;i;j E x x iÿ 1 = 2 ; y j ; 1 i m; 0 j n; E y;i;j E y x i ; y jÿ 1 = 2 ; 0 i m; 1 j n; E s;i;j E s x i ; y j ; 0 i m; 0 j n; (34) where x iÿ 1 = 2 x iÿ1 x i =2, etc.We apply the central difference for the second derivative with respect to x; y: As for the derivative with respect to s, we employ the From now on, we omit the symbol of radiation r and put the index ' for s.
The finite difference equation of ( 29) is then given by E ' Similarly, Eq. ( 30) becomes Assuming perfectly conducting walls, the tangent components of the field are zero on the walls.That is, the x component on the horizontal walls and the y component on the vertical walls are zero.Then they are given by E x and E y are decoupled because of the rectangular cross section.
Another boundary condition, valid for constant cross section, is E s 0 on the wall.

B. Initial condition and tracking
Let us give the initial condition at the entrance of a magnet.It is reasonable to assume that the field is in the steady state in the chamber with 1.We put @ s E x;y 0 in Eq. ( 29).Then we have to solve G x;yi;j H x;yi;j I x;yi;j 0 with the boundary condition.The leapfrog method requires two initial conditions at s s 0 and at s ÿ1 , but they are the same.
After the initial conditions are given, repeat the following procedure 1 to 4 in each step of s ' ' 1; 2; . . .; n.
(3) Calculate E s inside the chamber using Eq. ( 43). ( 4) Give the boundary values to E s .The whole procedure including initial conditions is repeated for different wave numbers k p for Fourier transformation.(This is unnecessary if only the impedance at one wave number is needed.)As is explained in the next subsection, the mesh size depends on the wave number k p .Hence we have to interpolate the field to get the fields on a common mesh for Fourier transformation using Eq. ( 2).

C. Conditions of mesh size
To execute the simulation, the most important factor is the mesh size.It should be determined taking into account two aspects.One is the accuracy, and the other is the stability of iteration.
The typical scale of the field in (x; y; s) direction is obtained by the dimensional analysis of Eq. ( 29).The mesh size must be considerably smaller than the corresponding scale length: x; y x typ ; y typ =k 2 1=3 ; (47) Note that these scale lengths are much larger than the wavelength 1=k.
On the other hand, the mesh size must satisfy another condition for stable iteration.For the analysis of the stability let us ignore the curvature term in Eq. (21): Let us consider a solution of the form From Eqs. ( 49) and (50), the dispersion relation of the leapfrog iteration is given by To keep the iteration stable, must be real for any real number and .Accordingly, the mesh size must satisfy the following relation: In numerical calculations, the mesh size must satisfy Eqs. ( 47), (48), and (52).Since we repeat the tracking for different wave numbers k p p 1; 2; . .., these conditions, assigned on the mesh size, change with k p .

IV. COMPARISON WITH ANALYTIC THEORIES
We show the results obtained by this method.The results are compared with analytic solutions derived by (A) Derbenev et al. [4], (B) Warnock [5], and (C) Saldin et al. [6].They are derived with different conditions as tabulated in Table I.We suppose that the line charge density has a Gaussian longitudinal distribution given by Eq. ( 23).We compare the longitudinal electric field E s obtained by our simulation with the analytic expressions.
In order to compare the numerical results with cases (A) and (C), the cross section of the chamber has to be large enough to avoid shielding effects.In the case of (B), only the side walls has to be far enough.We simulate the infinite space by choosing the chamber size as large as possible within the limitation w= p , h= p 1, where w is a full width, and h is a full height of the chamber cross section.Since the analytic formulas (A) and (B) describe the field of CSR in the steady state, we have to track the particles over a long distance, L m 23 2 s 1=3 .If the tracking length is too long, however, the radiation reflected on the wall can reach the observation points and disable the comparison.Therefore, we have to choose an appropriate tracking distance.
In the following plots the horizontal axis is the position in the bunch in units of the rms bunch length (bunch head to the right), and the vertical axis is the longitudinal electric field E s normalized by the factor E 0 2Nr e mc 2 = 2 p 3 2 4 s 1=3 .We used the rms bunch length s 0:3 mm and the curvature radius 10 m in all the simulations.
Unshielded CSR in steady states (A) is shown in Fig. 1.One finds the larger cross section gives closer results to the analytic formula.In the case where the chamber has a 34 cm 28 cm cross section (plotted with blue dots), the result agrees well with the analytic solution.The tiny residual difference behind the bunch is due to the reflection on the wall.
Shielded CSR in steady states (B) is shown in Fig. 2 for various sizes of the chamber height.In all cases the simulation results agree well to the analytic formula (see the Appendix about Warnock's formula).The residual small error can be attributed to the chamber of finite width in the simulation.
Unshielded CSR in transient states (C) is shown in Fig. 3 for various distance s from the entrance of the magnets.The size of the chamber cross section is w h 34 cm 28 cm in the simulation.The results again agree well with the analytic formula in transient states, though the case of s 100 cm shows a small difference.
To see the effects of the finite distance of the wall we simulated the case of the smaller chamber 10 10 cm and plotted the results in Fig. 4. Comparing the two cases, one can recognize that the residual difference from the analytic formula is in fact due to the finite distance to the wall.Thus, in all cases our method gives very good agreement with theories.Since existing theories demand infinite space in either transverse or longitudinal dimension, the computing time for the comparison is somewhat long (more than several hours).However, in an actual situation, where the 3D dimension is limited, the computing time is much less.In fact, 1 min is often enough in the case of strong shielding such as the parameters for storage rings.

V. CONCLUSION
We have shown that the CSR can be calculated by using mesh in a moderate computing time.The most important point is the paraxial approximation which enables the adoption of the mesh size larger than the relevant wavelength.We compared the simulation results with existing analytic formulas and found excellent agreements.Although we tried only simple cases for which an analytic formula is available, the method is very flexible and can be extended to more general cases such as the following: (i) finite transverse beam size, (ii) chamber cross section other than rectangular, (iii) chamber size changing gradually, (iv) finite beam energy (but still large ), (v) changing beam profile (but not due to the CSR itself).
It should also be possible to include the finite resistivity of the chamber wall by imposing the resistive boundary condition.
We have calculated only the energy change due to E s but it is also possible to compute the transverse force.This issue will be discussed somewhere else.Another issue to be discussed is the range of validity of the paraxial approximation.

TABLE I .
Conditions of the analytic solutions.
FIG. 1. (Color)The longitudinal electric field E s in steady state without shielding.The Derbenev's formula is plotted with a solid line.The dots are the simulation results for various sizes of the vacuum chamber corresponding to w h 6 6 cm (magenta), 8 8 cm (green), and 34 28 cm (blue), respectively.Length of the magnet is 2.1 m.E s is calculated at the center x y 0 of the chamber.FIG.2.(Color)The longitudinal electric field E s in steady state at x y 0 between infinite parallel plates.Warnock's formula is plotted with solid lines.The dots are the simulation results.The different colors show the gap between the two horizontal parallel plates, h 1:5 cm (black), 2.5 cm (red), 4 cm (magenta), 6 cm (green), and 8 cm (blue), respectively.The width of the chamber is w 50 cm and the length of the magnet is 3 m.FIG.3.(Color)The longitudinal electric field E s at x y 0 in transient state without shielding.The analytic solution is plotted with solid lines and the dots are the simulation results.The different colors indicate the tracking length from the entrance of the magnet: s 20 cm (black), 40 cm (red), 60 cm (magenta), 80 cm (green), and 100 cm (blue), respectively.The size of the chamber cross section is w h 34 cm 28 cm.