NEW METHOD OF CALCULATION OF THE WAKE DUE TO RADIATION AND SPACE CHARGE FORCES IN RELATIVISTIC BEAMS∗

Wepropose a new approach to the calculations of radiation and space charge longitudinal forces based on the use of the integrals for the retarded potentials. Our main result expresses the rate of change of particles’ energy through 2D (in a 2D model) or 3D integrals for a given orbit of the beam. It generalizes the 1D model and includes the transient effects at the entrance and the exit from the magnet. For a given beam line with known magnetic lattice, and a known distribution function of the beam, the calculation reduces to taking 2D or 3D integrals along the orbit.


INTRODUCTION
When the trajectory of a relativistic beam is bent by magnetic field, the beam radiates electromagnetic field and experiences a radiation reaction force. A popular 1D model for this force in the case of a circular motion, often called the coherent synchrotron radiation (CSR) wake, was developed in Refs. [1][2][3]. A generalization of this model for the case of a bending magnet of finite length has been done in Refs. [4,5] and is implemented in several computer codes 1 . These models are simple and easy to use but they miss an important part of the total force in relativistic beams moving in a curvilinear trajectory.
The attention to this force was attracted by M. Dohlus [6] in 2002, who pointed out that if the beam is compressed (either longitudinally or transversely) the energy of its Coulomb field changes and this should result in a change of the kinetic energy of the beam particles. A force that is responsible for this change can be called the compression force. Note that this force is different from the radiation reaction force, because the compression is a reversible process, and if the beam is decompressed, this force changes sign. It cannot be associated with what is conventionally called the space charge force because the latter typically scales as 1/γ 2 with γ the Lorentz factor. The compression-decompression effect occurs even in the limit γ → ∞ (hence, v = c), when the space charge force vanishes.
A part of this work has bean developed earlier by this author in collaboration with D. Ratner [7].
We use the CGS system of units throughout this paper. * Work supported by the Department of Energy, contract DE-AC03-76SF00515 † stupakov@slac.stanford.edu 1 Here, we only deal with the longitudinal part of the forces in the bunch that changes particles' energy in the beam; for the effect of the transverse force, see Ref. [8].

FORMULA FOR THE LONGITUDINAL WAKE
The beam is represented by its charge density ρ(r, t) that depends on time t and coordinate vector r, and its velocity (r, t), with the beam current density j given by the product j = ρ . Note that in this model assigning a particular value of at each point in the beam we neglect the spread in velocity due the the angular and energy spread in the beam -an approximation that is typically well satisfied for relativistic beams. For given functions ρ(r, t) and j(r, t) we then derive an equation for the electric field in the beam, E(r, t), and calculate the instantaneous energy change per unit time and per unit charge The result is expressed as an integral over the volume around the beam trajectory at preceding times t < t. It is remarkable that while the integrand in this integral has a singularity, it is integrable and is easily treated by standard numerical integration techniques. This is in contrast to the Green function approach in traditional simulations of 3D radiation forces that use a relativistic point charge fields with a strong nonintegral singularity at the location of the charge. We will loosely call P the longitudinal wake, although the classical wake fields are typically associated with the energy loss integrated over the beam path and the transverse cross section of the beam. Due to the lack of space we omit the derivation that expresses P in terms of the integrals over the whole space and give the final result: where β = /c. Note that due to the factor |r − r | −1 the integrand has a singularity at r → r, however, this singularity is integrable 2 . In many cases the last two integrals in Eq. (2) can be neglected. Indeed, the first integral involves the spacial derivative of ρ that can be estimated as |∂ r ρ| ∼ ρ/σ, where σ is the characteristic size of the beam. In the last two integrals we have the spacial and time derivatives of the 2 It is also integrable in a 2D model considered below, when the three dimensional integration velocity field in the beam, and these are estimated as |∂ r · β| ∼ 1/L and |∂ t ret β| ∼ c/L with L the external scale of the problem that is determined by the magnetic lattice (L can usually be associated with the radius of curvature of the beam orbit). If we assume that the size of the beam is much smaller than L, σ L, the last two lines in Eq.
(2) will be typically small. Neglecting them, we arrive to a much simpler expression for P, In next sections, we will show several illustrative examples of the wake calculation.

CSR WAKE ON A CIRCULAR ORBIT
The first example is the relativistic motion of a beam on a circular orbit of radius R specified by the following equations in the Cartesian coordinate system, where e 1 and e 2 are the unit vectors of an orthogonal Cartesian coordinate system in the plane of the orbit and s is the distance along the orbit. Assuming = c, the normalized velocity β is this model is equal to the unit tangential vector to the circle, β(s) = e 1 sin (s/R) + e 2 cos (s/R) (note that β does not depend on time). The horizontal coordinate x is measured in the orbit plane in the direction perpendicular to the orbit so that the vector r(x, s) is specified by r(x, s) = r 0 (s) + xn(s), where n(s) is the unit vector normal to the trajectory, n(s) = e 1 cos (s/R) − e 2 sin (s/R) . The beam distribution function is assumed a 2D Gaussian, with σ x and σ z the transverse and longitudinal rms sizes of the beam and Q the beam charge. For calculation of the wake we used Eq. (3) in which the integration goes in two dimensions with d 2 r = dx ds . We used the following set of parameters: R = 1 m, σ x = σ z = 1 mm. We first show in Fig. 1 the integrand in Eq. (3) for s = 0, t = 0 as a function of s and x . Note that transversely it is localized over the distance of several σ x , while longitudinally it extends over much larger length. This longitudinal extension is the formation length of the wake; it can be estimated from 1D CSR theory as (24R 2 σ z ) 1/3 ≈ 30 cm and shows a good agreement with the longitudinal extension of the integrand.
The plot of the wake function on the axis of the beam x = 0, at t = 0, is shown in Fig. 2. The dashed red line in this figure is calculated by 1D theory [2,3]. One can see a very good agreement of our calculation with 1D even though the bunch transverse and longitudinal sizes are equal.

TRANSIENT CSR WAKE OF A BEND MAGNET
In our second example, we consider the problem treated in Refs. [4,5]. In this problem we assume a bending magnet of finite length L with the beam trajectory inside the magnet, − 1 2 L < s < 1 2 L, given by the same Eq. (4). Outside of the magnet, the beam moves with = c along straight lines tangential to the circular orbit at the points of entrance and exit, respectively. The bunch distribution is given by Eq. (5) with s the length along the orbit and x the horizontal coordinate perpendicular to the orbit.
We have chosen the same set of parameters as in Ref. [5]: R = 1.5 m, σ z = 50 μm, Q = 1 nC and L = 25 cm. In the 1D model of [5] the parameter σ x = 0; in our calculations we have chosen σ x = σ z .
We first calculated the wake inside the bunch when it enters the bend from the straight line. The plot of this wake at various distances from the magnet entrance edge is shown in Fig. 3 by solid lines. For comparison, the dashed green lines show the result of 1D model computed in Ref. [5].
We have also calculated the wake in the bunch after it exits the bend magnet and continue motion also a straight

GAUSSIAN BEAMS AND BEAM OPTICS
In the two examples above the beam charge density and velocity were specified by equations which are very much evident due to simplicity of our model. In practical cases, where a beam can have an energy chirp and may travel through a complicated beam line, the calculation of the beam charge density and velocity is a more difficult problem. In this section, using the approach developed in Ref. [9], we will show how this problem is solved for the case of a beam with a Gaussian distribution function propagating through a system that can be described by linear optics. We limit our analysis to two dimensions -the horizontal and longitudinal ones.
To calculate ρ and we will use the formalism of the Vlasov equation. We use the notation x for the horizontal particle offset relative to the nominal orbit, θ = dx/ds is the angular slope of the orbit, η = ΔE/E is the relative energy deviation of the particle, z is the longitudinal coordinate of the particle in the bunch, and s is the path length along the nominal orbit. The beam distribution function F(x, θ, z, η, s) is a function of integrals of motion (see, e.g., [10]); it is normalized by ∫ Fdx dθ dz dη = 1. We will assume that F depends on the following three integrals of motions. The first one is the action connected to the betatron oscillation in the horizontal plane, where β f (s) is the beta function (not to confuse with the vector β = /c and its components β x and β s below), α(s) = − 1 2 dβ f /ds, D(s) is the dispersion function and D (s) = dD/ds. The second integral of motion is the energy deviation η, because we assume that the particle is not accelerated on the trajectory. Finally, the third integral, is obtained if one expresses the initial value of the coordinate z (say, at s = 0) though the values of z, η, x and θ on the orbit using the transport matrix R(s) (see details in [9]). Since this initial value is a constant, the resulting expression, which we denote by Z, is an integral of motion. Using the symplecticity of the matrix R, Z can be written as Z = z − R 56 (s)η + xR 26 (s)− θR 16 (s). By construction, R ij (s) with i j are equal to zero at s = 0, where the coordinate Z is equal to the longitudinal coordinate in the beam z.
For the equilibrium function F we choose a Gaussian distribution, where N b is the number of particles in the beam at the entrance (s = 0), is the horizontal emittance, σ z0 is the rms bunch length of the beam at s = 0 and σ η is the uncorrelated energy spread of the beam. The chirp parameter h in this equation accounts for the correlation between the position of the particle in the bunch and its energy. By integrating F we can find the two-dimensional charge density of the beam, ρ(x, z, s) = e ∫ Fdθ dη and its transverse velocity β x (x, z, s) = e ρ(x,z,s)

∫
Fθdθ dη, where we assume a relativistic beam, and use the approximation x ≈ cθ for the particle velocity. Because of the explicit dependence of the distribution function F versus the variables θ and η, this integration can be carried out analytically. Here we present the result of the integration for a special case 3 when R 16 (s) =R 26 (s) = 0, where and Σ 2 = D 2 σ 2 η σ 2 z0 + β f (hR 56 + 1) 2 σ 2 z0 + R 2 56 σ 2 η . (10) In Eqs. (8) we need to substitute z = s − ct; this makes ρ and β x functions of x, s and time t. After that, Eqs. (8) can be used to calculate the wake (2). Typically, we assume that β x 1 and treat β x as of the first order. Neglecting the second order terms we also use the approximation β z ≈ 1.