Efficient computation of coherent synchrotron radiation in a rectangular chamber

We study coherent synchrotron radiation (CSR) in a perfectly conducting vacuum chamber of rectangular cross section, in a formalism allowing an arbitrary sequence of bends and straight sections. We apply the paraxial method in the frequency domain, with a Fourier development in the vertical coordinate but with no other mode expansions. A line charge source is handled numerically by a new method that rids the equations of singularities through a change of dependent variable. The resulting algorithm is fast compared to earlier methods, works for short bunches with complicated structure, and yields all six field components at any space-time point. As an example we compute the tangential magnetic field at the walls. From that one can make a perturbative treatment of the Poynting flux to estimate the energy deposited in resistive walls. The calculation was motivated by a design issue for LCLS-II, the question of how much wall heating from CSR occurs in the last bend of a bunch compressor and the following straight section. Working with a realistic longitudinal bunch form of r.m.s. length 10 . 4 μ m and a charge of 100 pC we conclude that the radiated power is quite small (28 W at a 1 MHz repetition rate), and all radiated energy is absorbed in the walls within 7 m along the straight section.


I. INTRODUCTION
Effects of coherent synchrotron radiation (CSR) need to be studied in almost all advanced accelerator projects. In single-pass systems, for instance free electron lasers, the effects are usually deleterious, for instance in causing transverse emittance degradation in bunch compressors. In electron storage rings CSR may cause unwanted bunch instabilities, but may also be useful in providing intense radiation in the THz domain. In spite of ambitious efforts to provide models and computational tools to describe CSR, there is still a lot of room for improvements in models and algorithms. There are several distinct aspects to the problem: (i) modeling the charge/current distribution in the bunch; (ii) modeling the vacuum chamber; (iii) computing fields of very short bunches which may have microstructures; (iv) providing all the field components necessary to describe the experimental situation.
We believe that the paraxial method in the frequency domain offers opportunities for improvements in most of these aspects. It was introduced for the CSR problem in 2003-2004 by Stupakov and Kotelnikov [1] and Agoh and Yokoya [2,3], and pursued since then by Stupakov and Kotelnikov [4], Zhou [5], Zhou et al. [6], and Bizzozero et al. [7][8][9]. An outstanding feature of the method is that it works better as the frequency, equivalent to wave number k, is increased. This has roots in a relation to the ray picture of optics, and provides a chance to study the fields of very short bunches with possible interior microbunching. Of course, one has to demonstrate the practicality of discretizing in k over a wide range, and taking the inverse Fourier transform to construct fields in space-time. In the present paper these steps are shown to be entirely practical for bunch parameters close to the current state of the art, namely for a 10 μm bunch with a realistic longitudinal profile for a section of LCLS-II. We could go to even smaller bunch lengths and account for microbunching as well. Thus we have a response to item (iii) above.
For item (ii), the vacuum chamber, most of the cited papers assume a rectangular chamber with perfectly conducting walls. An exception is Ref. [1], which treats a round chamber. On the other hand, codes that aspire to the realistic simulation of beams usually model the vacuum chamber by infinite parallel plates, if a vacuum chamber is included at all. This might be because the paraxial codes with the rectangular model are thought to be too slow to apply in tracking of macroparticles. We show that the solution of the paraxial equations with the rectangular model can be speeded up greatly by a more efficient discretization in the transverse coordinates ðx; yÞ. Instead of a finite difference or finite element representation of derivatives in ðx; yÞ-space [2,3,5,6,8], or a two-dimensional mode expansion [4], we make a Fourier mode development in the vertical coordinate y, treating only x by finite differences. The motivation for this goes back to early work [10,11], in which it was noticed that the Fourier series in y converges quickly and also affords a way to enforce boundary conditions on the horizontal walls. In fact, bringing in the y-expansion establishes a fruitful connection to the soluble model of a complete circular torus and its whispering gallery modes. This connection shows the proper way to enforce low frequency cutoffs at modedependent "shielding thresholds." Because we avoid a mode expansion in x it is possible to treat a chamber of varying width in the x-direction, provided that variations are not too extreme. That treatment, motivated by our study of a flared chamber at the Canadian Light Source (CLS) [12], will be covered in a later report.
Another requirement for the CLS study was to find transverse electric field components at the location of a detector far off the beam axis. Off-axis fields are also needed in the present work to treat resistive wall heating. Thus we have examples of the requirement (iv) above, to find various field components at any point in the chamber. This requirement is met very easily in the present framework, because the Fourier transforms of all field components are expressed in terms of those for the y-components of E and H. After these two field components are computed by solving independent paraxial wave equations it costs almost nothing to find the other four. In particular, the transverse Lorentz force on the beam is readily available.
Item (i), the description of the source, is a crucial step in any scheme and deserves the closest attention. We look forward to a self-consistent scheme in which the bunch is modified by the field it produces, ideally through the Vlasov equation integrated by the method of local characteristics. In a less costly approach a macroparticle method replaces the Vlasov description. In any event we want at least a two dimensional representation of the bunch in ðs; xÞ space, where s is the longitudinal coordinate, in order to study bunch compression. As in the cited papers [1][2][3][4][5][6][7] we here consider only a fixed, one-dimensional charge/current, but in such a way as to clear the path toward higher dimensions. At the end of the paper we give our conception a full three-dimensional, self-consistent scheme.
Our charge density has the factored form qλðs − βctÞδðxÞHðyÞ, where the vertical distribution HðyÞ is fixed and has a finite width. Other authors put HðyÞ ¼ δðyÞ, which implies an infinite field at x ¼ y ¼ 0, and so are led to special arguments to extract the relevant finite part. It turns out that the finite width of HðyÞ is essential in constructing field components other than the longitudinal electric. This is a new story, which must be understood prior to attempting a full theory without factorization of the charge distribution.
In another important departure from earlier work we found an efficient way to deal numerically with δðxÞ in the source for the wave equation. By a simple change of dependent variable, the wave equation acquires a new effective source in which θðxÞ, the unit step function, appears instead of δðxÞ. Successive transformations can make the effective source arbitrarily smooth. We find that two transformations, for an effective source proportional to xθðxÞ, produces good results. There is a wider scope for this idea, since it can be generalized to smooth and broaden an x-distribution which is narrow but not a delta function. This would seem to be the proper approach for a self-consistent scheme describing a low emittance beam, better than the obvious idea of devoting more mesh points to a region near x ¼ 0.
It seems clear that our methods could provide a relatively fast self-consistent scheme with macroparticles. The cost of a competing method, which computes fields in space-time through retarded potentials and uses the parallel plate model, is strongly dominated by field evaluations rather than charge/current density construction and particle pushing [13]. At each point of observation of the wake field a separate integral over histories is performed, which is very costly. Our method provides a markedly faster field evaluation, while being the same for the other operations. It requires two Fourier transforms to go back and forth between space-time and the frequency domain, but those can be done by the Fast Fourier Transform (FFT) and should not be expensive.
The use of retarded potentials has been revived recently by Stupakov and Zhou [14], who find the longitudinal impedance for a line charge source in various idealized cases, assuming a parallel plate vacuum chamber. The impedance is in terms of double integrals, within a sum over image charges, which must be evaluated numerically. The local wake field is not obtained, only its integral over time. Our method generalizes directly to a higher dimensional charge distribution, and is more general even for a line charge in that it provides the local wake and allows the bunch profile to vary with s.
Our perturbative calculation of the Poynting vector for resistive wall heating follows an established idea [15]: replace H by its value H 0 for perfectly conducting boundaries, and replace E by its value from the resistive wall boundary condition, which is approximated by again using H 0 instead of H. As far as we know that scheme is always stated for single modes. Since we have a mode expansion in only one coordinate, we have to derive new formulas. On the horizontal walls we find interference between different vertical modes.
Section and subsection titles provide a synopsis of the paper. We include an appendix on the derivation of the resistive wall boundary condition, hoping to clarify certain issues that are not emphasized in textbook treatments [15][16][17][18][19]. In particular, we examine the basic assumption that variation of fields within the wall material is primarily in the direction normal to the wall.

II. WAVE EQUATION FOR THE SLOWLY VARYING AMPLITUDE IN ACCELERATOR COORDINATES
A. Fourier transforms in time and the y coordinate We work in standard accelerator coordinates (Frenet-Serret coordinates) defined in terms of a reference trajectory R 0 ðsÞ lying in a plane, and parametrized by its arc length s. Any spatial point in the laboratory system is represented as R ¼ R 0 ðsÞ þ xnðsÞ þ ye y where nðsÞ and e y are unit vectors normal to the unit tangent tðsÞ ¼ R 0 0 ðsÞ. To be definite we take the horizontal unit vector to be n ¼ e y × t. The vacuum chamber is to have a rectangular cross section with planar surfaces at y ¼ AEg, thus with full height h ¼ 2g. The vertical walls at are planar in straight sections and cylindrical with constant radius of curvature in bends. In this paper x þ and x − are constant. In a generalized treatment they can have a mild s-dependence. The chamber accommodates a beam centered at x ¼ y ¼ 0, following a sequence of straights and bends. Let Fðs; x; y; tÞ be any one of the six field components or a component of the charge/current density. Suppressing for the moment the variables ðx; yÞ we write F as the Fourier integral where βc is the velocity of the centroid (mean charge position) of the longitudinal charge distribution. By the Fourier inversion theorem which may also be written as an inversion with respect to z ¼ s − βct, This integral certainly converges in the situation we consider, with a single bunch making one pass through the system. At fixed s the field or source is negligible except for times close to the time of passage of the bunch, say when js − βctj is less than some multiple of the bunch length.
In definingF as the Fourier transform in time divided by e iks we hope to take out the fastest variation in s, leavingF as a slowly varying amplitude. By examining and solving numerically the differential equations for theF we shall show that this ambition can be realized for the parameters of our example. Then (2) will be a superposition of waves traveling in the positive s-direction, with slowly modulated amplitudes. It is also useful to view (2) as a description of the field in the beam frame at any fixed s, whether or not the amplitude is slowly varying. The where J and ρ are the current and charge densities of the beam. With these choices the Maxwell equations and boundary conditions are satisfied term-by-term in the sums over p. This follows from orthogonality,

B. Transformed Maxwell equations in Frenet-Serret coordinates
Next write the Maxwell equations using the standard expressions for divergence and curl in curvilinear coordinates. The metric tensor is diagonal with diagonal components ðg s ; g x ; g y Þ ¼ ðηðx; sÞ; 1; 1Þ; where κðsÞ and RðsÞ are the curvature and radius of curvature of the reference orbit at s. We substitute fields and sources in the form (5) and take the inverse Fourier transforms with respect to z ¼ s − βct and y to obtain the following system (in SI units).
These equations may be solved algebraically for all field components in terms ofÊ yp andĤ yp and their derivatives, yielding the resultŝ Here it is assumed that γ 2 p ≠ 0, a condition that is met in our calculations by a p-dependent low frequency cutoff or "shielding threshold." For a complete treatment allowing arbitrarily low frequencies one can give k a small positive imaginary part.
Equations similar to (17)- (20) are familiar in a scheme with separate Fourier transforms in s and t [10]. Fortunately, with only the one integral transform (2) we can still solve for all fields in terms ofÊ yp andĤ yp and their derivatives, but a new feature is that s-derivatives appear.
The fieldsÊ yp andĤ yp are obtained as solutions of two independent wave equations with sources. To derive the wave equations one may combine the transformed Maxwell equations as stated above, or proceed from the wave equations in Cartesian form and transform the differential operator to Frenet-Serret coordinates. The equation for The factor κ 0 ðsÞ in (22) is nonzero where the reference trajectory (which need not be an actual particle trajectory) changes from straight to curved or vice versa. If the change is abrupt at s ¼ s 0 then κ 0 ðsÞ contains δðs − s 0 Þ, and it is doubtful that the wave equation can be given a meaning in a neighborhood of that point. On the other hand, if we give κðsÞ a smooth transition over a distance comparable to a typical fringe field extent in a bending magnet, then each of the terms with κ 0 is small compared to the term immediately preceding it in (22). Accordingly we drop κ 0 terms but then allow κ to be a step function at bend-straight transitions, elsewhere in the equation. The solutionF p is required to be continuous at transitions.
Henceforth we drop the transverse currentsĴ xp ,Ĵ yp , but these could be restored in a more ambitious self-consistent scheme.

C. Slowly varying amplitude approximation and the simplified wave equation
We now proceed to the main approximation, which is to assume that the amplitudeF p in (5) is slowly varying as a function of s. We may state the criterion for slow variation in terms of a norm, for instance where dependence of f on variables other than the transverse coordinate x is suppressed. Then the requirement on F p ðk; s; xÞ is over the interval of integration ½s 0 ; s 1 . This does not make sense as k → 0, but there is a lower bound to the relevant k values, as will be shown presently. We shall actually monitor the condition (25) in our calculations, which apparently has not been done before in similar CSR studies. For convenience (25) is called the slowly varying amplitude (SVA) approximation or paraxial approximation. In our view the former name is more apt, since it reminds us of the only condition that need be enforced in the present framework. Now within a bend of constant bending radius R the simplified wave equation (22) takes the form The corresponding equation in a straight section is obtained in the limit R → ∞ as where γ is the Lorentz factor, 1 − β 2 ¼ 1=γ 2 . These equations are sometimes described as "parabolic," but that is a misnomer. They are of Schrödinger type owing to the factor i, with mathematical properties different from those of a proper parabolic equation. Next we choose a simple factored form for the charge density of the beam, good enough for the present limited study but capable of being generalized. With the corresponding current density it is ρðs; x; y; tÞ ¼ qλðs − βctÞδðxÞHðyÞ; Jðs; x; y; tÞ ¼ ðβcρ; 0; 0Þ; ð30Þ where q is the total charge. The continuity equation is satisfied. By (5), (6), and (7), the Fourier transform with respect to z and y iŝ ρ p ðk; xÞ ¼ qλðkÞH p δðxÞ; λðkÞ ¼ 1 2π Z e −ikz λðzÞdz; Thus the sources in (27) becomê S Ep ¼ qZ 0 α p cλðkÞH p δðxÞ; For an even distribution HðyÞ ¼ Hð−yÞ we have so that H p is zero for even p and alternates in sign for successive odd p. For H we try two examples at opposite extremes in their large-p behavior, a Gaussian of zero mean and r.m.s. deviation σ y ≪ h, and a square step with the same mean and deviation, thus for p odd For the longitudinal distribution we apply the result of a simulation for LCLS-II, to be described presently. A comparison can be made to the Gaussian and square step cases with formulaŝ

III. NUMERICAL SOLUTION OF THE SIMPLIFIED WAVE EQUATION
An elementary way to approach the solution of (26) or (28) is to discretize the right-hand side on a grid in x-space, representing the x-derivatives by finite differences. The discretization involves values of the solution at the boundaries,F p ðk; s; x AE Þ, which are to be fixed at values required by the boundary conditions at a perfect conductor. The equation is then regarded as a system of ordinary differential equations, with s as the independent variable, in the complex unknownsF p ðk; s; x i Þ; i ¼ 2; …; N − 1.
Here the x i are the interior points of the x-grid. The system is treated as an initial value problem, the initial value being the s-independent solution in an infinite straight wave guide.

A. A transformation to mollify the effective source
There is an impediment to discretization, however, due to the δðxÞ and δ 0 ðxÞ in the source terms (33). By a change of the dependent variable, this source can be replaced by a new effective source which behaves as θðxÞ near x ¼ 0, where θðxÞ is the Heaviside step function, A second transformation gives the continuous function xθðxÞ in the effective source, and successive transformations can make the source arbitrarily smooth. ForF p ¼Ê yp we write the expression in square brackets on the right-hand side of (26) or (28) so as to emphasize the x-dependence, suppressing variables ðk; pÞ. In terms of the differential operator the expression in square brackets is uðxÞ ¼Ê yp ðk; s; xÞ; where in the bend and in the straight section, with R ¼ ∞, Now define a new dependent variable u 1 ðxÞ by We see that ξ 00 1 ðxÞ ¼ c 1 δðxÞ cancels the original source c 1 δðxÞ and we gain a new effective source S 1 : Since ∂u=∂s ¼ ∂u 1 =∂s, the field u 1 satisfies the same partial differential equation as u, but with the new source which is much more suitable for discretization, being piecewise continuous with a jump of −c 1 að0Þ at x ¼ 0.
Numerical integration of the differential equations with this setup was found to be only partly successful, at least when done by the elementary finite difference method described below. An instability was encountered at large p in some cases, a behavior that could be traced to the jump in the source. To remove the jump we again change the dependent variable to u 2 defined by This yields a source behaving as xθðxÞ with It is clear that this process can be continued for additional smoothing. At the nth stage the source is It is determined from S n−1 through the transformation The divided difference in the first term in (49) is analytic at x ¼ 0, as a result of aðxÞ and bðxÞ being analytic at that point. A similar procedure works for the magnetic field, even though the source to start with is more singular. For F p ¼Ĥ yp the expression in square brackets on the right-hand side of (26) or (28) has the form uðxÞ ¼Ĥ yp ðk; s; xÞ; The first transformation to remove δ 0 also removes δ because of the special form of aðxÞ. Thus The δ drops out because að0Þ − 1=R ¼ 0 in both the bend and the straight. As in the discussion above, the second transformation will be and so on. Now let us summarize the net effect of two smoothing transformations, invoking the explicit forms of the coefficients a and b. The smoothed field ξ H Þ to make the full field ðÊ yp ;Ĥ yp Þ, and the effective source is denoted bỹ S p ¼ ðS Ep ;S Hp Þ. In the bend, The corresponding formulas for the straight section are obtained in the limit R → ∞, noting that bðxÞ tends to bð0Þ in the limit: Note thatÊ yp andĤ yp must be continuous at the transitions between bend and straight, while the corresponding u E , u H are not continuous. This must be kept in mind in designing the algorithm for s-integration.
In (18) and (20) we have the factorĴ sp − ∂ xĤyp , wherê J sp ¼ c 2 δðxÞ. Fortunately, the c 2 δðxÞ is cancelled by the term c 2 ∂ x θðxÞ in ∂ xĤyp as given by (55). Such a cancellation was noticed long ago in analytical models [10], but a good way to handle it in a numerical context was lacking before the present innovation.

B. Comparison to the beam field subtraction of Agoh and Yokoya
Agoh and Yokoya [3] also made a transformation to mollify the source, which is similar in spirit to ours but different in details. They work in ðx; yÞ-space, whereas we work just in x-space after Fourier analysis in y. They ⊥ where in their notation the so-called "beam field" is a particular solution of ð∂ 2 x þ ∂ 2 y ÞE ðbÞ ⊥ ¼ μ 0 λðkÞ½δ 0 ðxÞ; δ 0 ðyÞ. The effective source in the paraxial equation for E ðrÞ ⊥ , their equation (29), is then free of singularities. The analogous decomposition in our setup isÊ yp ¼ u 1 þ ξ 1 in (42), with ξ 1 being a particular solution of ∂ 2 x ξ 1 ¼ c 1 δðxÞ. Thus our lowest order smoothing is entirely analogous to the idea of Agoh and Yokoya, but there seems to be no analogy in higher orders.

C. Smoothing transformation for a beam of nonzero horizontal extent
A more realistic charge/current distribution will have a nonzero extent in x but could still be very narrow, for instance 20 μm in our example from LCLS-II. To handle that case numerically one can generalize the method described above. For an arbitrary source SðxÞ we wish to transform Φ ¼ Lu − S. The first transformation will be uðxÞ ¼ ξ 1 ðxÞ þ u 1 ðxÞ; Now only the integral of S appears in the effective source, a smoother and more extended function than S itself. Of course further transformations could provide additional smoothing, as we have seen.

D. Boundary conditions at the vertical walls
With perfect conductivity the boundary conditions at the vertical walls arê From (17) and (20) we see that these conditions are met if E yp satisfies a Dirichlet condition andĤ yp a Neumann condition, namelŷ The corresponding conditions on the smoothed fields u ¼ ðu E ; u H Þ follow from (54) and (55):

E. Finite difference scheme
We suppose that the field values are interpolated by 4th degree polynomials in x, and that derivatives are given by differentiating the interpolation. The 4th degree interpolation [20] of a function fðxÞ on a grid fx j g N j¼1 with uniform cell size Δx is with Lagrange polynomials The error ϵ is O½ðΔxÞ 5 , and is estimated in terms of the 5th derivative [20]. For evaluation at interior points of the grid respectively, for the necessary off-center interpolation. Differentiating (65) with respect to ξΔx gives the formulas for derivatives. Define Now we can write the discretized form of the wave equation (26) forÊ yp as follows, in terms of the smooth field u E : In view of (63) the boundary values that appear in these sums are where the inner and outer boundaries are at The equation forĤ yp in terms of the smooth field u H has the same form, with the appropriate definitions from (55), except that the boundary values are expressed in terms of interior values by discretizing the Neumann conditions (64): Solving for the boundary values we have These values allow the numerical derivatives, as in (69), to be expressed in terms of interior values of u H alone. In a straight section the discretized equation (68) reduces to.
with the definitions of (57) and the boundary conditions of (70). The equation for u H is the same with the definitions (58) and the boundary conditions of (72).
The following calculations will be for a single bend followed by a straight section. The generalization to treat an arbitrary sequence of bends and straights is obvious, and would not make a great complication in coding because one has only to specify the curvature κðsÞ to define the equations at any s.

F. Initial values for the evolution in s
The system of linear differential equations (68) is to be solved as an initial value problem. We take the initial value for s ¼ 0 at the beginning of the bend to be the steady-state field produced by the source in an infinitely long straight chamber. Thus the equation for an initial fieldF p is (28) Its general solution is a particular solution plus the general solution of the homogeneous equation,F in which A and B must be chosen to meet the boundary conditions. With the notation defined in (39) and (50) we haveŜ Ep ðxÞ ¼ c 1 δðxÞ,Ŝ Hp ðxÞ ¼ c 2 ðδðxÞ=R þ δ 0 ðxÞÞ. Evaluating the integral in (75) and applying the boundary conditions (62) we find Note that the definition θð0Þ ¼ 1=2 in (37) makeŝ For the numerical work it is essential that the definition of θð0Þ be the same in the initial condition as in the smoothing transformation.
The corresponding initial values of the other field components are derived from (17)-(20): where γ is the Lorentz factor. The mechanism for the expected small value ofÊ sp at large γ (in accord with the familiar disk-like picture of the field pattern) is the near cancellation of the terms fromÊ yp andĴ sp − ∂ xĤyp in (17). The cancellation becomes less precise during field evolution in the bend, butÊ sp is still a small difference of two large terms. A numerical difficulty arises in the application of (76) and (77) because of a close cancellation of large terms at large x ≈ x þ . The increasing part of the second term in (76) or (77), namely expðα p xÞ=2, cancels against a part of the first term. By some rearrangement we take out the cancelling terms and find the following formulas, suitable for numerical evaluation: Even after this step one must take care to avoid overflow or underflow in evaluation of the exponentials, by appropriate expansions.

G. Evolution in s
Suppressing irrelevant variables we write the system of differential equations for evolution of where u and f are vectors with N − 2 complex components, and f is linear in u. For the approximation at s ¼ s ðnÞ ¼ nΔs þ s ð0Þ we write u ðnÞ ≈ uðs ðnÞ Þ, where the integration step Δs is allowed to be different in bends from what it is in straight sections. We adopt the leapfrog integration rule, based on the central difference approximation to the derivative: To define u ð1Þ for the first step we use Euler's rule, As remarked above, the value of u at the end of a bend is not in general equal to the value of u at the beginning of a following straight, owing to a change in definition of u through source smoothing. Consequently, we use an Euler step to initialize a leapfrog integration in the straight, with the appropriate initial value defined by continuity of the physical (unsmoothed) field at the bend-straight transition.
Of course there are more powerful methods than the finite difference method for discretizing in x and the leapfrog method for s. We have chosen these simple schemes merely to make our strategies clear and to avoid complications in programming for this exploratory study. We have in fact compared results from a more sophisticated x-discretization using the discontinuous Galerkin method [7,9,21], as will be reported below. Future work should look for a method with a good compromise between speed and accuracy.

IV. POYNTING FLUX AT THE WALLS TO LOWEST ORDER
The Poynting vector E × H evaluated at a wall describes, through its outwardly directed normal component, the flow of energy into that wall, per unit area and per unit time. At a perfectly conducting wall E is normal to the wall while H is tangential, so the normal component of the Poynting vector vanishes. The resistive wall boundary condition (A13) implies a tangential component of E at the wall and a nonzero energy flow. We can calculate this flow to lowest order from a knowledge of H 0 , the magnetic field computed for perfectly conducting walls. We replace H by H 0 in both the second factor of the Poynting vector and in the boundary condition. In this approximation the Poynting vector S ¼ E × H at a point r ¼ ðs; x; yÞ on the wall is From here on we writeĤ forĤ 0 in accord with the notation of previous sections.
Since we are interested in the total energy loss we may integrate over t. Note that Hð−k;rÞ ¼Ĥðk;rÞ Ã ; ð87Þ Here the integrand has finite support in t because the fields follow the source, and are negligible for js − βctj greater than some length L, the maximum range of wake or predecessor fields. Now notice that since H satisfies the boundary condition for a perfect conductor, with zero normal component. Moreover, ð1 − iÞk 1=2 goes into its complex conjugate as k → −k, since k 1=2 → ijk 1=2 j as we have defined it in the complex plane in Appendix. Then in view of (87) the integral on k is twice the real part of the integral on positive k and We see that the time-integrated energy flux is solely along the normal direction and is positive toward the wall at all points (since n is directed inward toward the vacuum).
Next we wish to integrate (90) over one transverse dimension at the walls; namely, over y at x ¼ x AE for vertical walls and over x at y ¼ AEg for horizontal walls. By (5) and (6) At the vertical walls we can use the orthogonality of (7) to find the y-integral as At the horizontal walls the x-integral is dtSðs; x;AEg;tÞ with the same formula holding forĤ x . To find the total energy deposited in the walls the expressions (94) and (95) must be integrated over s using the numerical solutions for the tangential H fields.

V. TOTAL ENERGY RADIATED AND THE WAKEFIELD
Here we derive the formula for the total energy radiated, for comparison to the amount of energy absorbed in resistive walls. By conservation of energy this is just the It follows that the power radiated from all elements is For comparison to earlier work we are also interested in the longitudinal wakefield, Wðz;s;x;yÞ ¼ 2Re We evaluate this at x ¼ 0 and take its mean value with respect to the vertical charge distribution HðyÞ to obtain Wðz; sÞ ¼ 2gRe VI. NUMERICAL RESULTS

A. Parameters and bunch profile for LCLS-II
We present some numerical results, mostly with parameters anticipated for the final bend of the second bunch compressor (BC2) in the forthcoming Linac Coherent Light Source II (LCLS-II) at SLAC National Accelerator Laboratory. The beam is centered in the rectangular chamber of width w ¼ 2x þ and height h. The chamber dimensions, bending radius, and bend angle are w ¼ 5 cm; h¼ 2 cm; R¼ 12.9m; θ ¼ 42.5 mrad: The chamber material is copper; we take the conductivity to be σ ¼ 5.96 × 10 7 Ω −1 m −1 . The single-bunch charge, energy, and nominal r.m.s. beam dimensions are q ¼ 100 pC; E¼ 1.6 GeV; There is also a mode with q ¼ 300 pC and a longer bunch, which is of less interest for this study. The repetition rate can be as great as 1 mHz. We calculate the fields and energy loss for a single bunch. In view of the small value of σ x we represent the bunch charge density as in (30), with zero width in x. For the fields to be finite at the beam it is then necessary that there be a nonzero spread in the vertical density HðyÞ. We compare the Gaussian and square-step distributions with Fourier transforms (36). For the longitudinal distribution λðzÞ we apply the result of a realistic simulation, and compare the outcome to that for a Gaussian with the same σ z . A smoothed representation of the simulated distribution is shown in Fig. 1; it has zero mean. The simulation gave a histogram of 100 bins, which was smoothed by convolving with a quartic kernel having a half-width of 3 bins and zero slope at the ends. The real and imaginary parts of the Fourier transformλðkÞ of the distribution are shown in Fig. 2, along with the corresponding transform of a Gaussian with the same σ z . The smoothing has no discernible effect on the Fourier transform in the range of k plotted.
We set β ¼ 1, but our code allows β < 1 which is at least of interest for a general understanding of fields and code diagnosis, if not for immediate applications.

B. Initial conditions and convergence of the vertical mode expansion
We first show the behavior of the vertical fieldŝ E yp ðxÞ,Ĥ yp ðxÞ from which all other field components  (20). The integration on s begins at s ¼ 0, the beginning of the bend, with the steady state solutions for the straight pipe given in (76) and (77). In Fig. 3 we show the initialÊ yp ðxÞ for the first four modes (p ¼ 1, 3, 5, 7) and the sum over all modes, E y ðk; x; yÞ ¼ X pðoddÞ cos½α p ðy þ gÞÊ yp ðk; xÞ; ð107Þ evaluated at the upper boundary y ¼ g where the cosine is −1. The only k-dependence is through the dimensionless factorλðkÞ which is omitted in the plots. Convergence of the p-sum at x ¼ 0 happens only by virtue of the decay of H p , which is not appreciable until large p because of the small value of σ y . For the Gaussian we have convergence (as judged by graphical inspection) by p ¼ 139, whereas for the square step we go to p ¼ 1999. The limit is the same for the two cases, evidently because the field far from the beam at y ¼ g is not sensitive to the vertical charge density. The spikes at x ¼ 0 get narrower with increasing p, providing convergence at x ≠ 0 (but not uniform in x) without the help of the factor H p . Short of the limit there is a narrow spike in the sum, which alternates in direction (up or down) as each new mode is added.
evaluated at y ¼ 0 where the sine is equal to the alternating factor ð−1Þ ðp−1Þ=2 . SinceĤ yp has the same alternating factor arising from H p , the summand lacks the alternating sign that appeared in (108). Consequently the modes just add up coherently, giving the limit shown in Fig. 4 (right). Note that we have included the factor Z 0 because Z 0Ĥyp andÊ yp have the same units (volts) and can be compared in magnitude.

C. Characteristics and evolution of the various field components
Next we show the evolution of the initial fields in the bend, for a particular value of the wave number k within the important range of wave numbers. Figure 5 shows the real and imaginary parts of the sum to p ¼ 9 ofÊ yp ðk; s; xÞ at y ¼ g, for s at the end of the bend and kR ¼ 5 × 10 5 . The spike in the real part is inherited from the initial field, and as before it alternates in sign with the number of terms in the sum. As in the initial field, convergence is achieved only at very high p. The spike is absent in the imaginary part. Figure 6 shows the decomposition of ReÊ yp into the smooth numerically computed function Reu E and the given function ξ E , as in (54), for the case of Fig. 5. The benign behavior of Reu E illustrates the value of our change of dependent variable.
Fortunately, we can deal with convergence at high p because the high-p modes do not evolve appreciably in the bend and the following straight section. Starting at p ¼ 11 or so the change of a mode is so small as to be nearly invisible in a graph. Hence our procedure will be to represent the high-p part ofÊ y ,Ĥ y and their x-derivatives by their given initial values.
The number of oscillations increases with k. Increasing k by a factor of 4 we get the plot (of the real part) ofÊ yp in Fig. 7.
Next we show graphs for the other field components that enter the calculation, for the same parameters as used in Fig. 5. Figure 8 showsĤ y at y ¼ 0 at the end of the bend. The magnitude of the jump in the real part at x ¼ 0 will grow as the sum on p is extended, to approach the value seen in Fig. 4 (right). Figure 9 shows the longitudinal electric fieldÊ s evaluated at y ¼ 0. It shows no trace of a spike at x ¼ 0, even though it is constructed from (17) with contributions from E yp and ∂ xĤyp , both of which have big spikes at x ¼ 0.
The spikes cancel each other in a way that is certainly remarkable but mathematically obscure. One can understand the cancellation of high-p contributions, which is the same cancellation that occurs in the initial condition to give a residual of order 1=γ 2 seen in (78). The cancellation at small p is the puzzling point. In any event, just a few terms in the p-sum are needed to representÊ s and the longitudinal wakefield derived from it.
A colleague has asked whether this cancellation of peaks persists for a beam with β < 1. This is a worthwhile question for a future study with β < 1.
In Fig. 10 we have Z 0Ĥx evaluated at y ¼ g. In this component the spikes do not cancel, and the behavior of the spikes is like that ofÊ y ; a hint of this can be seen in the initial value (81). In computing the resistive wall loss on the horizontal walls the high-p sum must be included inĤ x to get convergence.
Finally in Fig. 11 we show Z 0Ĥs evaluated at y ¼ g. As in the case ofÊ s there is no spike at x ¼ 0, even though the separate terms in Eq. (19) do display spikes. Note thatÊ s and Z 0Ĥs are relatively small, in accord with their initial values, Oð1=γ 2 Þ and 0, respectively. The satisfaction of appropriate boundary conditions is apparent in the various graphs.
In Fig. 9 we see something like the whispering gallery picture, in which the field is concentrated near the outer wall and nearly absent at the inner wall. The concentration would be more pronounced if the bend were longer. Following theÊ s field into the following straight section we find that this pattern disappears, and the field spreads over the whole width of the chamber. Figure 12 and Fig. 13 give views at 1 m and 5 m along the straight. Notice a signal of increasing complexity in that the initial 90°phase shift between real and imaginary parts has disappeared by 5 m.

D. The p-dependent shielding cutoff at low frequencies
The present scheme works best at large k. As we see in (25), the slowly varying amplitude approximation stands to fail at sufficiently small k. Evidently a low frequency cutoff is required, but experience showed that this cannot be independent of p. We seek guidance on how to place the cutoff from the analytically soluble model of a complete circular torus with rectangular cross section [10]. The wave equation that defines that model is the Bessel equation with source, where definitions are the same as in (26) except that kR ¼ n is quantized to be an integer, so that solutions are periodic around the torus. The corresponding homogeneous equation has solutions J n ½γ p ðx þ RÞ, Y n ½γ p ðx þ RÞ and resonances occur near frequencies where J n ½γ p ðx þ þ RÞ ¼ 0 or J 0 n ½γ p ðx þ þ RÞ ¼ 0. Solutions at resonance are called whispering gallery modes and are concentrated near the outer wall and are very small at the inner wall. Bessel functions have oscillatory behavior, allowing zeros at the outer wall, only when their argument is greater than their order. The transition from exponential to oscillatory behavior, where argument equals order, coincides with vanishing of the coefficient ofF p in (109) at x ¼ x þ . This is analogous to the change in solutions of the harmonic equationü þ ω 2 u ¼ 0 when ω 2 changes sign. Thus a necessary condition for resonances is that Invoking the definitions and making an expansion for small x þ =R we cast this in the form This "shielding cutoff" for a toroidal or pillbox chamber was noted in [10] [Eq. 5.15 of SLAC-PUB-4562 or Eq. (90) of the published paper]. For p ¼ 1 and 2x þ ¼ h it agrees with the better known cutoff for the parallel plate model of the vacuum chamber, derived in [22]. A detailed analysis showed that there are no significant contributions to wakefields from frequencies below the resonance region [23]. In our single-pass system with moderate bend angle there are no sharp resonances, but there are broad peaks in the frequency spectrum of fields that apparently are vestiges of resonances, and also field patterns with some resemblance to whispering gallery modes. The resemblance to the full torus sharpens as the bend angle increases. Our example from LCLS-II has both small bend angle and large bending radius in comparison to most examples studied in the past, and is relatively remote in behavior from the full torus or "steady state." We should therefore check to see if the cutoff (111) looks reasonable by comparing field patterns below and above cutoff.
In Fig. 14 we show results for p ¼ 9, usually the highest necessary vertical mode, plotting the longitudinal electric field as a function of x as it looks at the end of the bend. Graph (b) is for kR at the cutoff given by (111). The pattern far above cutoff in graph (d) can be considered a whispering gallery mode. Such a pattern appears to be emerging in graphs (b) and (c), but is absent in graph (a) below cutoff. Moreover, graph (a) is smaller, somewhat noisy, and shows poorer convergence as the mesh in x is refined. Notice the 10-fold larger scale in graph (d), which is for k in the important spectral region for the wakefield. Similar results are found at smaller p.
We conclude that a cutoff slightly less than (111) will allow even incipient whispering gallery behavior, while excluding low frequency effects that will have negligible effect on wake fields and energy radiated. The excluded contributions are analogous to the small subresonant effects studied in [23]. To check the effect of frequencies somewhat below the cutoff (111), say down to the range of graph (a), we insert a reduction factor c r on the right-hand side of (111) and experiment with its value. We find that results for wakefields and energy radiated and absorbed are just the same for c r ¼ 0.5 and c r ¼ 1. Thus it could be said that the vanishing of low frequencies effects is abrupt, just as in the full torus and parallel plate models.
We incorporate the cutoff (111) and find that our code has excellent convergence properties in all the parameters controlling discretization. This was not always the case with a p-independent cutoff.

E. Energy radiated and deposited in resistive walls
Having indicated how the field components look in the frequency domain, we now turn to the integrals over all frequencies that give the energy radiated and absorbed in resistive walls. The blue curve in Fig. 15 shows the total energy radiated up to position s in the perfectly conducting model, as given by (102). The beginning of the bend is at s ¼ 0, where the fields are assumed to have the steady-state values for an infinite straight wave guide. The red curve Fig. 15 shows the energy deposited in resistive walls up to position s in the perturbative approximation. This is obtained from the sum of the s-integrals of (95) and (94), for horizontal and vertical walls respectively. The magenta and brown curves give the separate contributions of horizontal and vertical walls.
The slope of the blue curve is nearly zero at s ¼ 7.6 m, meaning that the decelerating field of CSR in the perfectly conducting model has nearly died out. This decay, a purely geometric effect, might be seen as a developing incoherence between different Fourier components of the longitudinal electric field at x ¼ 0. At large s any one Fourier component is spread over all x, does not decay with s, and acquires a random looking relation between its real and imaginary parts.
When the red curve crosses the blue curve at s ¼ s c ¼ 7.6 m all energy radiated at s < s c has been absorbed, to the accuracy of the perturbative model. In the present example, that happens to be almost all the energy that will be radiated. With a lower wall conductivity it would be somewhat less than that, since the absorbed energy is proportional to σ −1=2 .
Since a nonzero tangential H at the walls persists for s > s c the red curve continues to rise. This corresponds in exact physics to the fact that after all of the CSR is absorbed there is still a huge kinetic energy of the beam, some of which will be dissipated in the walls. Even though the perturbative model does not account self-consistently for decreasing energy of the beam, it is correct in predicting ongoing absorption in the walls. Asymptotically this should increase linearly with s, and indeed we find highly linear behavior of the red curve when the integration is extended by another 2m. The slope in that range is 5 μJ=m.
Since the total CSR energy deposited from one bunch is 28μJ, and there are about 10 6 bunches per second, the deposited power is about 28 W. There is 5 W=m additional dissipated power for s > s c .
It is interesting to look at the field patterns in the region where a major part of the radiation occurs, say for s < 2 m. In Fig. 16 we show contour plots in ðx; yÞ-space of H x ðs; x; y; tÞ (in arbitrary units) at the instant the bunch passes position s. We show this field component because its value on the horizontal walls is the largest contributor to wall heating. We plot in just the upper quarter of the chamber, g=2 < y < g, since the range of values is too great to make a good plot in the full cross section. Corresponding surface plots are presented in Fig. 17. The flux of energy into the upper horizontal surface, integrated over time, is shown as a function of x at various s in Fig. 18. This function is given by (95) without the integration on x. The contribution ofĤ s to these plots is negligible.

F. Longitudinal wakefield and its Fourier transform
Most earlier work on CSR has concentrated on fields at positions within the bunch. Usually only the longitudinal wakefield is computed, but on occasion transverse forces have been studied as well [13]. Following this tradition we plot in Fig. 19 the longitudinal wake field at the end of the bend, for the simulated bunch form and a Gaussian with the same σ z . The red curves represent the bunch form on an arbitrary scale. The head of the bunch is at positive z ¼ s − βct, and positive W corresponds to energy gain. We warn that this is just the electric field per unit charge at a fixed s as a function of z, not the integral up to s that is sometimes seen in the literature.
The corresponding Fourier transforms of WðzÞ are displayed in Fig. 20. The wakes at two values of s beyond the bend are shown in Fig. 21. With increasing s there is decreasing energy loss within the bunch, while field structures outside the bunch are built up.
The peak several bunch lengths in front of the bunch in Fig. 19 is a parameter-dependent feature, which is not seen in several published plots for which the bend angle was larger [4][5][6]. A far-forward peak has turned up, however, in other publications; see Fig. 3 and Fig. 4 in [3] and Fig. 17 in [24]. Furthermore, Stupakov and Emma did an analytic calculation of the wakefield in free space [25], which showed such a peak evolving with position in the bend.
To show how a similar evolution works in our example for LCLS-II, we take a Gaussian bunch rather than the simulated bunch, to get a cleaner plot which is easier to understand. We also extend the length of the bend by 50% to s b ¼ 0.825 m. In Fig. 22 we show the wake field at successive positions in the bend. At s ¼ 0.1375 m the wake (blue curve) looks like the derivative of the bunch form. Subsequently the forward peak in this curve moves farther forward, and then splits into two. The forward peak of the resulting pair moves ever farther from the bunch centroid, decreasing in height, and leaves the domain of the graph before s ¼ s b . Presumably, the missing far-forward peak in certain publications is just out of range of the graphs.

G. Code validation and timing
We have applied two codes developed independently by the two authors. The results presented above are from the first code, which is written in Fortran and implements the numerical methods and equations we have described. The second code is written in Matlab and is based on the same simplified wave equation and bunch description, but it makes a more sophisticated discretization in x, by the discontinuous Galerkin (DG) method [9,21]. The DG scheme allowed first order smoothing of the source, as in (42), whereas second order smoothing was required for the finite difference method. Results of the two codes for the wake of the Gaussian in Fig. 19 agreed perfectly. The DG method looks to be superior, but needs to be implemented in a faster programming language.
We have also made comparisons to results in the literature, for instance to the wakefield in Fig. 3 with a Gaussian bunch having σ z ¼ 0.5 mm. This example was for the steady state case, computed by thorough mode expansions and requiring attention to poles on the real axis in frequency. To approximate the steady state we take a large bend angle of π=2 and find the wake at the end of the bend shown in the blue curve of Fig. 23. The corresponding result of [4] is in red. Considering that the two calculations were done by vastly different methods, the fairly close agreement seems satisfactory. We regard our calculation as simpler; for one thing, it does not require principal-value integrations around poles.
For the finite difference code the parameters controlling discretization are Δs, Δx, p max , k max , Δk. By spot checks we have verified good convergence of various results in each of these parameters separately. The rate of convergence depends somewhat on the quantity computed. For the wakefield at the end of the bend, in the LCLS-II example, a well-converged run could have parameters around the following: where s b ¼ 0.55 m is the length of the bend. The integration step Δs has more to do with stability of the s-integration than with accuracy. If it is small enough to ensure stability then making it even smaller does not change results appreciably. Examining the right-hand side of (26) we see that an estimate for the Courant-Friedrichs-Lewy stability criterion is where α < 1 is to be determined empirically. We found that α ¼ 0.2 worked in a few cases. The value of Δs in (112) is for α ¼ 0.24 and k ¼ k min from the lowest shielding threshold. For simplicity our code takes all discretization parameters to be fixed, but it could be made more efficient by increasing Δs with k, while keeping Δx fixed.
With the control parameters of (112) the computation of the wake field in Fig. 19 (left) took 4.2 minutes on a PC (Intel i7-4790, 3.6 GHz) and included evaluation of the wakefield at 400 z-points at each of 400 different values of s ≤ s b . (It also included 400 evaluations of the energies radiated and deposited in resistive walls, not needed for beam dynamics.) These results were achieved with a suboptimal serial code. In parallel processing the s-integration could be done independently for each mode pair ðk; pÞ, with a trivial calculation of initial data.
To secure convergence of the wakefield within the bunch, the important region for beam dynamics, one can make a choice of discretization parameters more economical than (112), say by decreasing k max and increasing Δx, the latter allowing a bigger Δs. In any case the code timing seems very promising for an application of our algorithm as a field solver for a self-consistent macroparticle simulation, or even a Vlasov calculation.
At large s out to 8m we have to take a smaller Δk than in (112) in order to get a smooth curve of energy loss as in Fig. 15, say Δk ¼ k max =1000. Otherwise we get a curve with about 3% jitter at large s, meandering about the smooth curve.

H. Check of the slowly varying amplitude approximation
To verify the condition (25) for validity of the SVA approximation we take ∂Ê yp =∂s from the right-hand side of (68) or (73), and approximate the second derivative as a divided difference. Then we plot the ratio r E ðsÞ of the lefthand side of (25) to the right-hand side, and the similar ratio r H for the magnetic field. The ratios are largest at small k and large p, so we put p ¼ 9 and k at the shielding threshold (111) for that p. This "worst case" for r E is shown in Fig. 24 (left). In a more important range of smaller p and larger k the values shown in Fig. 24 (right) are typical. The large step in values occurs at the bend-tostraight transition. Overall it appears that the SVA approximation is very well justified, for bothÊ yp andĤ yp , for the mild bend of the present example. For a large bend angle and small bend radius the justification is not so clear. For the case of Fig. 23 with p ¼ 5 and k at the shielding threshold we find r E ¼ 0.13 near the beginning of the bend.

VII. ANTICIPATED ALGORITHM FOR AN EVOLVING SELF-CONSISTENT CHARGE/CURRENT DENSITY
Here we sketch a scheme to accommodate a general charge/current density evolving under its self-field, in place of the idealized model treated heretofore. Since the scheme has yet to be tested, changes may be needed after actual experience. Here we are concerned only with the effect of CSR and space charge on the particle distribution in phase space, not with resistive wall heating. We adopt a macroparticle description, while noting that a more ambitious Vlasov description could be pursued as well. We have in mind a single-pass system such as a chicane bunch compressor or a "dog leg" transport line, leaving aside for now a possible application to storage rings.
Every source and field component is regarded as a function of ðz; s; x; yÞ with z ¼ s − βct, where βc is the mean longitudinal particle velocity. As before we make Fourier developments in z and y, and all equations from (4) to (28) apply. Since we no longer invoke the simplified model of (30), the Fourier transforms of charge/current are now general, unfactored functions of s and x, and the current may have nonzero transverse components: where Sðz; x; yÞ is a normalized smoothing kernel [26], β i c is the velocity vector of the ith particle, and q is the total charge. The kernel will have small support in each of its arguments, extending to a small multiple of the average interparticle spacing in that argument.
Given the macroparticles at s, these densities will be computed on a grid in ðz; x; yÞ-space, uniform in each variable. Their Fourier amplitudes (30) can then be calculated by a two-dimensional FFT in z and y, thus for all values of k and p at once, for each x on the grid.
The densities will be updated only at intervals δs of s that are very much larger than the integration step Δs for solving the differential equations (68): δs ≫ Δs. Recall the example in (112), with 30000 steps in s through a bend. We hope that only a few hundred updates of the charge/current would be necessary in a bend.
Given the fieldsÊ yp ðk; s; xÞ,Ĥ yp ðk; s; xÞ, particle positions r i ðsÞ and velocities β i ðsÞc, the update of these quantities to their values at s þ δs proceeds as follows: (1) Compute the sums (115) and (116) on the grid in ðz; x; yÞ. (2) Compute the Fourier amplitudes (114) by an FFT for each x, and from those form the sources (23) for the differential equations. (3) Integrate the differential equations forÊ yp andĤ yp over an interval δs in steps of Δs, regarding the sources as independent of s on that interval. (4) Using the resultingÊ yp ðk; s þ δs; xÞ,Ĥ yp ðk; s þ δs; xÞ and Eqs. (17)- (21), form the Fourier amplitudes of the field components necessary to advance the particles, for instancê E sp ðk; s þ δs; xÞ. (5) Take the inverse FFT of the latter to get the Lorentz force in ðz; x; yÞ-space. (6) Apply the Lorentz force to advance the particle orbits by the equations of motion, with s as the time-like independent variable. With the resulting r i ðs þ δsÞ and β i ðs þ δsÞ, return to step (1) to continue the evolution in s.
Note that this describes a full field calculation, accounting for space charge effects along with CSR.
In many cases it may be important to apply the smoothing transformation as explained in Sec. III C. Then step (2) would be followed by the smoothing transformation, and in step (3) the differential equations would have the new effective source. For more information on using s as the independent variable for equations of motion, as well as information on distributions in the beam frame variable z, see Ref. [13].
In a typical realistic calculation, the initial particle distribution will come from an independent simulation, perhaps representing the particle distribution at the end of a linac. In such a case we propose to define corresponding initial field valuesÊ yp0 andĤ yp0 for our calculation as follows: define an s-independent charge/current equal to the charge/current for the given initial distribution, and takê E yp0 andĤ yp0 to be the steady state fields in an infinite straight wave guide for that constant charge/current. Those fields can be computed from (75).
On evidence from the present work we expect step (3), the field solution, to cause a minor part of the computational cost. We also conjecture that step (1) will dominate the total cost, for suitable values of the particle number N .

VIII. CONCLUSIONS AND OUTLOOK
We have described an effective scheme for fast numerical computation of CSR in a vacuum chamber of rectangular cross section. The central step is solving a simple system of linear ordinary differential equations which describe evolution in s of a slowly varying wave amplitude. The system is autonomous within a bend, Eq. (68), or within a straight section, Eq. (73). The source termS p is a continuous (or piecewise-continuous) function of x, obtained from the original line charge source by a change of dependent variable.
We have applied our method to find the pattern of energy deposited in resistive walls. To our knowledge, this effect of CSR was not previously studied. We have treated the Poynting flux to the wall just to lowest order in resistivity, which is good enough to establish that the effect is small. An interesting next step would be a direct solution of the Maxwell equations under the resistive wall boundary condition, a calculation required to find the rigorous resistive wall wake field. This has been done analytically for the toroidal model [10], and for special geometries with rectilinear beams [27,28], but never for CSR in a vacuum chamber with successive bends and straights.
In this study we have learned three basic techniques of our chosen computational scheme: effective source smoothing by a change of dependent variable, the proper treatment of high p, and the proper low frequency cutoff. It remains to refine the numerical integration algorithms in several directions, and to take advantage of parallel processing.
A promising project is to apply this method as a field solver for self-consistent macroparticle simulations with a very short bunch, for instance for fields in a chicane including space charge as well as CSR. This should include a test of kernel smoothing [26] as an elegant alternative to conventional particle-in-cell procedures for smoothing the charge/current density. Our method also offers opportunities to improve the study of CSR and beam dynamics in storage rings, especially to clarify the question of inter-bunch communication through long range wakefields from whispering gallery resonances [29].

ACKNOWLEDGMENTS
We thank James Ellison and Jack Bergstrom for a great deal of help and inspiration over many years. Paul Emma raised the question of resistive wall heating, and kindly provided a bunch simulation. Gennady Stupakov suggested the perturbative approach to energy flux, gave us good advice, and provided results of his code to compare with ours. The use of Frenet-Serret coordinates in this study came out of conversations with Rui Li. The

APPENDIX: RESISTIVE WALL BOUNDARY CONDITION
We recall the derivation of the resistive wall boundary condition, adapting it to our particular context. We find the wave equation for the magnetic field within the wall material, which is assumed to have magnetic permeability μ, electric permittivity ϵ, and conductivity σ, all independent of position and frequency. The basic assumption is that the variation of the field within the wall is by far the strongest in the direction normal to the wall. This picture can be checked a posteriori by first assuming it to be true, then deducing the consequent normal variation. This variation, characterized by a small penetration depth (skin depth) can be compared with estimates of variation in the tangential directions.
Invoking Ohm's Law J ¼ σE, we have the curl equations within the wall as Next we take the Fourier transform with respect to t to obtain The term −iωϵ from the displacement current is tiny in comparison to σ at the highest frequencies we consider, and will be dropped henceforth.
We define a positive depth coordinate ξ, the distance from the beginning of the wall to an interior point of the wall medium, and a unit vector n normal to the wall and directed from the wall toward the vacuum. At the horizontal walls ξ ¼ AEðy − gÞ, whereas at vertical walls ξ ¼ AEðx − x AE Þ. Then with the assumption of dominant normal variation the gradient is represented as We can then eliminateẼ in (A6) by taking the curl of the first equation and substituting the second: Since ∇ ·H ¼ ∂ðn ·HÞ=∂ξ ¼ 0, we have The general solution of this harmonic equation with complex frequency is H ¼ a þ expðξ=ΔÞ þ a − expð−ξ=ΔÞ; The a AE depend only on coordinates other than ξ. Since the solution must decay at large ξ > 0 we retain only the second term and choose the branch of the square root so that ReΔ > 0, namely as where the square root in (A10) is positive at positive real ω. We define this root in the complex ω-plane with a cut on the positive real axis. It then acquires a factor of i in analytic continuation to negative ω, so that Δ −1 has positive real part at negative as well as positive ω. The conventional skin depth d is defined by so that the field decays by a factor 1=e in a distance d. By (A9) we have ∂H=∂ξ ¼ −H=Δ which when substituted in the first equation of (A6) yields Taking the limit ξ → 0 in (A12) we have the resistive wall boundary condition, since there must be continuity with the fields in the vacuum. The Fourier transform (A5) with respect to time is related to the slowly varying amplitudeF by the phase factor expð−iksÞ=βc, which cancels out in (A12). Thus with ω ¼ βkc and the near-perfect approximation μ ¼ μ 0 , the boundary condition for the Fourier amplitudes used in this paper iŝ Eðk;s;x;yÞ ¼ ð1 − iÞ βZ 0 k 2σ 1=2 n ×Ĥðk;s;x;yÞ; ðA13Þ at every point ðs; x; yÞ on the boundary, with the unit normal n to the boundary directed toward the vacuum. The skin depth with μ ¼ μ 0 is d ¼ ð2=βZ 0 kσÞ 1=2 . To test the assumption of dominant normal variation, we assume that any tangential variation would not be faster within the wall than it is at the surface. We can then estimate the scale of transverse variations at the surface from the perfectly conducting model, and compare it to the skin depth. We first compare d to the scale of variation in the s-direction, which should be about 1=k. Since d decreases with increasing k, an upper bound on d=k −1 will be its value at k max , the largest relevant k. To decide on the latter we examine the Fourier spectrum of field components, for instance jĤ x ðk; s b ; 0; gÞj as a function of k as plotted in Fig. 25 (left), or the bunch spectrum in Fig. 2. The most important range of the spectrum is for kσ z < 3, but there are substantial contributions out to kσ z ¼ 8 or more. Taking k max σ z ¼ 8 we find For y-variation the corresponding ratio of interest is where we use the shielding threshold (111) for k min ðpÞ. Recall that p max ¼ 9 in our calculations. For corresponding estimates of variation in the x-direction, we refer to the numerical calculation ofĤ x ðk; s b ; xÞ as a function of x.
That has rapid oscillations when k is large, as shown in Fig. 25 (right) for kσ z ¼ 8. We define δxðkÞ as 1=4 of the period of the oscillations. A rough fit shows that δxðkÞ decreases more quickly than dðkÞ, at about the rate k −1.125 , which means that dðkÞ=δxðkÞ will have its maximum value at k max . Reading off δx from graphs for k ¼ k max ¼ 8=σ z we find dðkÞ δxðkÞ < 0.021: This decreases to 0.014 at kσ z ¼ 3. To summarize, it appears that the scale of tangential field variations is at most about 2% of the skin depth for the parameters of our example. This justifies the assumption of dominant normal variation, but not by the huge margin that might have been expected.