Dynamics of Longitudinal Phase-Space Modulations in an rf Compressor for Electron Beams

Free-electron lasers operating in the UVor x-ray radiation spectrum require peak beam currents that are generally higher than those obtainable by present electron sources, thus making bunch compression necessary. Compression, however, may heighten the effects of collective forces and degrade the beam quality. In this paper we provide a framework for investigating some of these effects in rf compressors by focusing on the longitudinal dynamics of small-amplitude density perturbations, which have the potential to cause the disruptive appearance of the so-called microbunching instability. We develop a linear theory valid for low-to-moderate compression factors under the assumption of a 1D impedance model of longitudinal space charge and provide validation against macroparticle simulations.

up to a factor 3 or 4). A second approximation results from adopting a 1D model for the space charge impedance describing collective effects, which at smaller wavelenghts and low energy is known to be problematic[? ? ]. Finally a coasting beam approximation is assumed, which limits the applicability to particle dynamics in the longitudinal core of the physical electron bunches.
To validate the analytical model and provide a check on the simplifying approximations, we also present a comparison against macroparticle simulations. We find that a reasonable agreement with the theory can be obtained over a wide spectrum of perturbation wavelengths provided that certain empirical modifications in the parameters for the longitudinal space charge impedance be introduced. We ascribe the need for this adjustment to the limitations of the adopted 1D model for space charge.
The paper is organized as follows. After deriving in Sec. II the transfer map for single-particle motion through an rf structure in the linear approximation, in Sec. ?? we present the small-amplitude perturbation theory. The main result is the derivation of an integral equation obeyed by the beam bunching function expressing the evolving amplitude of an initial sinusoidal density modulation through the compressor. We follow with a discussion of the numerical model used for the macroparticle simulations and a comparison simulations against theory (Sec. ??).

II. SINGLE-PARTICLE DYNAMICS THROUGH AN RF COMPRESSOR
We assume that compression takes place in a travelling wave rf structure. The longitudinal motion of an electron in such a structure is described by the Hamiltonian H = m 2 c 4 + p 2 s c 2 − eφ (s, t) where −e is the electron charge, s the longitudinal coordinate, and φ = (E 0 /k rf ) cos(k rf s − ω rf t + ψ 0 ) the electric potential, yielding the longitudinal electric field E s = −∂φ/∂s = E 0 sin(k rf s − ω rf t + ψ 0 ). For simplicity, in the following we assume a phase velocity ω rf /k rf = c for the travelling wave. The canonical equations are readily written: ds/dt = p s c 2 / m 2 c 4 + p 2 s c 2 and dp s /dt = −eE 0 sin ψ, where we have introduced the phase ψ = k rf s − ω rf t + ψ 0 .
We then change the dynamical coordinates from s and p s to t and the relativistic factor γ = m 2 c 4 + p 2 s c 2 /mc 2 , while turning s into the independent variable. The resulting equations for t(s) and γ(s) read where we have introduced the parameter α = eE 0 /(mc 2 k rf ). The orbit for the reference particle is a particular solution of (1) and (2), which we denote as (t r , γ r ). The orbit of any other particle can be described in terms of the deviation variables ∆t = t − t r and ∆γ = γ − γ r , where (t, γ) are also solutions of (1) and (2). It is convenient to introduce the space separation ∆z = z − z r between an electron and the reference particle, in place of ∆t. We have ∆z = −cβ(s)∆t, where cβ(s) is the beam velocity. The negative sign results from the convention that for a particle in the head of the bunch ∆z > 0. In the following for the purpose of determining ∆z from ∆t, we will assume that the beam is sufficiently relativistic so that β 1, yielding ∆z −c∆t.
The equations of motion for the deviation variables ∆z and ∆γ read The first-order solution of the above linear system can be expressed in terms of the transfer matrix M : with initial condition M (s 0 ) = 1 and matrix A defined by Incidentally, the form of the matrix (6) allows us to infer immediately that the transformation M is symplectic.  (1) and (2). The left picture shows the evolution of the relativistic γ-factor, the right picture the relative longitudinal distance between the reference electron and an imaginary particle travelling at the speed of light. Examples of solutions of (1), (2) and (5), for the reference orbit and entries of the transfer matrix M are shown in Fig.'s ?? and ?? respectively. Here and in the examples to follow we make reference to a 3 m long, S-band, SLAC-type, travelling wave (TW) section currently in use for the SPARC experiments [? ]. The rf structure frequency is 2856 MHz (corresponding to λ rf = 10.5 cm or k rf =59.8 m −1 ). We assume a E 0 = 25 MeV/m peak accelerating gradient (yielding α=0.82). The rf structure is preceded by a 0.6 m long drift and the electron beam injected with 5.6 MeV kinetic energy (γ r0 12). The reference orbit is best represented in terms of the variable ζ r , defined as ζ r = s − ct r , (Fig. ??, picture to the right). This quantity expresses the relative longitudinal distance between the reference electron and an imaginary particle travelling at the speed of light. The initial conditions were specified so that as the reference particle enters the structure (s = 0.6 m) we have ζ r = 0, i.e. the rf phase is ψ = ψ 0 .

III. LINEAR THEORY
In the following we will refer to the dynamical variables in vector notation as x = (∆z, ∆γ). Moreover, to simplify notation, we will use z to denote ∆z and denote ∆γ with p, i.e. x = (z, p).
We assume that the beam distribution f (x 0 ; s 0 ) = f 0 (x 0 ; s 0 ) + f 1 (x 0 ; s 0 ) at the entrance of the rf compressor s = s 0 consists of a zero-order smooth density, uniform in z and Gaussian in p with a chirp h 0 and a first-order perturbation f 1 (x 0 ; s 0 ). The normalization of the distribution function is chosen so that n 0 dz ∞ −∞ dp 0 f (x; s) yields the number of particles in the interval dz, where n 0 is the beam line density. Let us consider the evolution of the unperturbed beam density first. The beam density function f 0 (x s ; s) at s is related to the beam density f 0 (x 0 ; s 0 ) at s 0 by f (x s , s) = f (M −1 x s ; s 0 ). Here and in the following we use M , M (s), or M (s 0 → s) interchangeably to denote the transfer matrix solution of (5) from s = s 0 to s and use the notation M (s → s) to represent the matrix advancing the linear solutions of (3) and (4) The normalized charge density evolves from ρ(z 0 ; The last equality in the above equation identifies the compression factor C = C(s). In writing (??) we made use of the symplectic property of the matrix M . The effect of collective forces described by an impedance Z(k; s) (to be specified later) is to change the particle energy according to where I 0 = ecn 0 is the electron beam current, I A = ce/r e 17kA the Alfvén current, and Z 0 the vacuum impedance, withρ with the collective force F depending on the FTρ 1 of the first-order density perturbation. Starting from (??) our goal is to derive an equation forρ 1 . To this end it is convenient to think of both sides of (??) as functions of the dynamical variables at current position s: where The next step is to make a more definite assumption about the form of the initial perturbation. We assume an expression of the form i.e. consisting of a sinusoidal perturbation to the charge density (while the p density is the same as in the unperturbed distribution). As usual in this kind of calculations the physically meaningful component is the real part of the complex quantity (??).
We are now ready to integrate both sides of Eq. (??) over the phase space after multiplying by e −ikz s . By definition, the FT of the LHS of Eq. (??) yieldsρ 1 (k; s). The FT of the first term on the RHS, which we denote with I 1 , requires more work. First, we carry out the transformation of variables x s → x 0 and exploit symplecticity d 2 x s = d 2 x 0 to write where z s (x 0 ) = [M x 0 ] 1 = M 11 z 0 + M 12 p 0 . We then we insert expression (??) into (??) and introduce the change of variable t = p 0 − h 0 z 0 , yielding with δ(·) denoting the Dirac function.
As for the second term on the RHS of Eq. (??) we have: The second equality above follows from a change of variables, the third from an integration by parts and z s ( where we have made use of f 0 (x s ; s ) = f 0 (x 0 ; s 0 ). Observe that in the argument of the exponential function of the above expression we have Using the expression (??) for the collective force we find Again, we introduce the change of variable p 0 → p 0 = t + hz 0 and observe that The integral on z 0 in (??) yields a Dirac function, allowing for a straightforward integration on k . After integration on t, we obtain Finally, upon combining the various expressions evaluated so far we arrive at the following integral equation for the Fourier components of the first-order density perturbation: with kernelK (s , s) = 4πi where I(s ) = I 0 C(s ) is the beam current at s . We look for solutions of (??) in the space of generalized functions (distributions). Such solutions will have the form with the ordinary functions b obeying the equation with kernel K(s , s) obtained from (??) with the substitution k → k 0 C(s ). The integral equation (??) is our main result. Somewhat unexpectedly, but not too surprisingly, (??) turns out to be formally identical to the equation describing bunching in magnetic compressors[? ? ]. We define the linear gain as the ratio of the amplitude of the perturbation at the exit s = s f to that at the entrance (s = s 0 ): To summarize the meaning of the calculation carried out in this Section: in linear approximation a sinusoidal perturbation with wavenumber k 0 to a beam charge density at the entrance of an rf compressor will maintain its sinusoidal form while the wavenumber evolves according to k = C(s)k 0 . The quantity b(k; s) or 'bunching function', representing the amplitude of the sinusoidal perturbation relative to the local beam density ρ(s) = C(s)ρ(s 0 ), is determined by solving the integral equation ( where Z 0 120π is the vacuum impedance, and I 1 and K 1 are the modified Bessel functions of first and second kind. Expression (??) is obtained[? ] from transverse averaging of the longitudinal component of the electric filed of an infinitely long beam with circular cross-section of radius r b (and uniform transverse density) perturbed by a small longitudinal modulation of wavenumber k. In principle, the radial dependence of the longitudinal electric field could be accounted for in the present framework but at the cost of increasing the dimensionality[? ] of the integral equation (??). Notice that in (??) the dependence on the longitudinal coordinate s is through the parameter r b = r b (s), assumed to be known.
In Fig.'s ?? and ?? we show a numerical example for a model of beam line consisting of a 1 m drift followed by a 3 m rf compressor. The beam injected at 5.6 MeV energy with I 0 = 100 A peak current and vanishing uncorrelated energy spread. We assume a beam transversely uniform with circular cross section with radius r b = 0.5 mm remaining constant through the beam line. The compression factor at the exit of the rf structure is 1.87. The gain curve as a function of the wavelength of the initial perturbation is reported in Fig. ??, and shows a maximum value of about 0.7 in the λ 300µm region. The gain is <1 in the whole range of modulation wavelengths considered. Gains larger than unity could be observed but only in the presence of substantially higher current or compression factor.
Because of space charge an initial modulation in the charge density will induce an energy modulation along the beam. We are interested in determining to first order the amplitude of this latter modulation. To this end we introduce a new set of dynamical variables to remove the (in general nonvanishing) correlation energy-positioñ Following the assumptions made in the previous section, we consider a beam distribution function f = f 0 + f 1 in the form of the sum of a uniform zero-order density f 0 and a first-order perturbation f 1 . The amplitude of the energy modulation is obtained by taking the FT the energy centroid along the beam p = f (z,p; s)pdp/ f (z,p; s)dp, a function ofz: where the last equation follows from first-order approximation and The normalization factor f 0 (z,p; s)dp is just the compression factor C(s). To calculate f 1 (z,p; s)pdp we make use of the first-order perturbation equation (??) for f 1 , yielding After some some algebra, here omitted for brevity, we find FIG. 5: Amplitude of the energy modulation amplitude along the rf compressor resulting from an initial density modulation with 10% amplitude relative to an initial current I 0 = 100A. Other beam and rf compressor parameters as in Fig. ?? and ??.
with the bunching factor b(k 0 C(s ); s ) determined by solving Eq. (??). In Fig. ?? we show two examples of energy modulations resulting from sinusoidal density perturbations with initial relative amplitude A = 0.1 and wavelengths λ = 75 and 150 µm.

IV. VALIDATION AGAINST MACROPARTICLE SIMULATIONS
As a way to validate the model presented in the previous section, we carried out macroparticle simulations using the code TSTEP[? ], a derivative of PARMELA [? ]. We considered the evolution of a 1 nC idealized flat-top beam spanning a 10 • rf phase at 2856 MHz (corresponding to about 3 mm and I 0 = 100 A peak current) and a range of initial charge-density sinusoidal perturbations with wavelength λ between 50 and 300µm. The upper limit of this range is determined by the need to consider wavelengths small enough compared to the bunch length so that the coasting beam approximation assumed in the analytical model may apply. The smallest modulation wavelength is limited by the statistical noise associated with use of a relatively small number of macroparticles. In the simulations we used up to 4.5×10 6 macro-particles, resulting in an acceptable compromise between numerical accuracy and computation time. The amplitude of the initial sinusoidal perturbation was set to A = 10%.
We considered the 3 m long rf structure mentioned at the end of Sec. II preceded by a 0.6 m long drift. The presence of the drift is a realistic feature of any physical set-up as the rf structure for compression compressor requires a certain separation from the exit of the gun. However, in these simulations, the exact value of the drift length was chosen to correspond (in the range of wavelength we considered and adopted initial beam conditions) to roughly a quarter wavelength of longitudinal plasma oscillation as this choice tends to maximize the amplitude of the modulation amplitude at the exit of the compressor (having started with pure density perturbations, i.e. no initial energy modulation).
The linear gain is calculated as the ratio between the amplitudes of the charge density perturbation at the exit of the compressor and entrance of the leading 0.6 m drift. Care was taken to limit the analysis of the numerical data to the core of the bunch to minimize edge effects. The space-charge forces were calculated by solving the Poisson equation on a grid with a 5 mm longitudinal span and number of mesh cells varying between 1200 and 2400 along the longitudinal coordinate.
We start the simulations with a beam with transverse uniform density and circular cross-section of initial radius r b = 2σ x = 2σ y = 2 mm and a convergent envelope dr b /ds < 0. It turns out that the plasma oscillation wavelength has a fairly strong dependence on the beam transverse radius r b and is therefore affected by the exact value of the initial beam convergence. This is exemplified in Fig. ?? where we show how different choices of the initial transverse conditions for the beam envelope affect the subsequent evolution of the beam radius (picture to the left) and the location of the first minimum for the linear gain for an initial perturbation of wavelength λ = 100µm (figure to the right). The minimum in the linear gain corresponds to a quarter of plasma oscillation wavelength from the entrance of the drift, which is where the initial density modulation is converted into energy modulation. See evidence of this in Fig.??. For the remaining simulations presented in this paper we adopted the initial conditions for the beam envelope corresponding to the solid line in the left picture of Fig. ??. For simplicity, in the simulations we did not include solenoidal focusing along the rf structure, which would be required for emittance compensation [? ], and let the beam expand freely.
Incidentally, notice that in the figures in this section starting from Fig. ?? the variable z in the abscissa is the same as the arclength coordinate s introduced in Sec. II.
Compression is controlled by moving the linac rf phase away from the crest toward the zero crossing of the rf field. In this section we refer to the rf phase as defined as Ψ = −90 • − ψ, where ψ is the phase defined in Sec. II).
To highlight the importance of a full account of space-charge effects we show a comparison between results obtained with two different settings where collective forces are turned on and off in the rf structure (while they are included through the leading drift in both cases), see Fig. ??. The initial rf phase is Ψ 0 = −82 • corresponding to a compression factor ∼2). It is seen that neglecting space charge in the rf structure (data points with round markers) leads to a gross overestimate of the gain.
As the beam size variation affects in a sensitive way the results (see Fig. ??), it is important that when we make comparison with the theory we account for the evolution of the beam radius along s. To this end a high-order polynomial interpolation for the rms transverse size as a function of s was carried out from data points extracted from the simulations. From analysis of the simulation data it was found that the best agreement between the analytical model of Sec. ?? and simulations is obtained when in the expression for the impedance (??) we use the relationship r b = aσ x with the factor a defined as a = 1.95−0.001×λ[µm] instead of a = 2, as it would be expected for transversely uniform beams. The good agreement between simulations (solid line) and theory (dots), including the empirically adjusted factor a is shown in Fig. ??. The good agreement also extends to the determination of the compression factor, Fig. ??, (with the compression factor from the simulation data calculated as the ratio σ z0 /σ z between the initial σ z0 and final σ z rms bunch lengths).
Further simulation-vs.-theory comparisons are reported in Fig.'s ?? and ??. In particular in Fig. ?? we show the linear gain at the exit of the compressor over a range of perturbation wavelengths for a uncompressed (Φ 0 = 0 • ) and compressed C=2 (Φ 0 = −82 • ) beam. Finally in Fig. ?? we report the evolution of the amplitude of the energy modulation induced by an initial density modulation with 74 µm wavelength. The energy modulation amplitude was retrieved from the simulation data by first removing the correlation phase-energy in a window selected around the beam core in the longitudinal phase space and then carrying out a sinusoidal fit. An example of longitudinal beam phase space exhibiting the energy modulation is shown in Fig. ??. Notice that as it evolves the beam develops a finite energy spread due to the radial dependence of the space-charge fields, which is missed by the 1D model of impedance adopted in our linear model. Nonetheless, the linear model, including the empirical adjustment of the a-factor mentioned earlier, appears to quite well the amplitude of the energy modulation.

V. CONCLUSIONS
In this paper we have derived a linear theory for the amplitude gain of charge density modulations of a beam passing through an rf compression system and have discussed a comparison against results from macro-particle simulations carried out with the code TSTEP. A satisfactory agreement with the linear theory is found when a parameter in the adopted 1D model for space charge effects is adjusted to account for an empirically determined dependence on the perturbation wavelength.
The numerical solutions of the linearized equations and the TSTEP simulations indicate that for parameters of interest in typical FEL applications, the amplification of small initial density perturbations through an rf compressor tends to be quite modest (relative to the peak current) if not outright smaller than unity (i.e. implying damping of the initial perturbation) even in the absence of any uncorrelated energy-spread induced mixing. This result is not completely unexpected. We know, for example, that in a low-energy beam drifting in free space, and hence without compression, the amplitude of longitudinal plasma oscillations remains constant. It turns out that for moderate compression, velocity bunching changes the dynamics of longitudinal plasma oscillations from the case of free-space drifting only mildly -the main difference being in the adiabatic lengthening of the plasma oscillation wavelength and reduced strength of the space charge forces as the beam undergoes acceleration during compression. These results do not imply that the dynamics of small-amplitude density perturbations in the rf compressor should be neglected altogether as these perturbations can seed an instability downstream if further compression by magnetic chicanes is applied. Indeed, a scenario in which rf compression is supplemented by magnetic compression represents the most likely mode of operation envisioned for the accelerator drivers of FEL-based 4th generation light sources [? ]. We expect that the theory elaborated in this paper will represent a useful tool in the evaluation of the compression schemes for these FEL sources. Finally, we should mention the possible relevance of our finding in connection with recent investigations of noise suppression in FEL injectors, see Gover et al.[? ] and Nauser et al. [? ].