Microbunching instability in single-pass systems using a direct two-dimensional Vlasov solver

We apply a recently developed Vlasov solver to the study of the microbunching instability generated by shot noise in the beam delivery systems of x-ray free electron lasers (FELs). We discuss two lattices presently under consideration for the FEL Fermi project at Elettra and show that at least one of the two lattices appears capable of delivering a beam with the desired quality in the longitudinal phase space.


I. INTRODUCTION
The microbunching instability in linacs can limit the performance of single-pass x-ray FELs by significantly degrading the beam quality [1]. The instability stems from small irregularities in the longitudinal charge density, which can be amplified by self-field induced energy variations when the beam travels through dispersive regions. Accurate and efficient modeling of this instability is important for linac design but is challenging because of the high resolution needed to capture small fluctuations in phase space. Macroparticle simulations are valuable but using a number of macroparticles significantly smaller than the bunch population introduces spurious noise that may overshadow the small fluctuations responsible for the genuinely physical instability [2].
Immune to this problem direct methods [3][4][5][6][7][8][9] to solve the Vlasov equation offer an interesting alternative. In [10] we proposed a 2D Vlasov solver tailored to the specific features of beams (like a strong energy/position correlation) in the delivery systems for x-ray FELs. Direct methods are generally more computationally intensive than macroparticle simulations but the strong dependence of the computational load on the dimensionality of the problem makes them particularly attractive for 2D simulations.
Being 2D, the solver we developed would seem to apply only to beams with vanishing transverse emittance. At first this would not appear to be a significant limitation as the microbunching instability is essentially a longitudinal phase-space phenomenon, where longitudinal charge density fluctuations give rise to an electric self-field affecting the beam energy through the longitudinal component. However, the longitudinal slippage, which is central to the development of the instability, can be substantially affected by a finite horizontal emittance. A realistic representation of the beam dynamics should then include some account of the horizontal degree of freedom. In [10] we proposed a model for the smearing effect of a transverse emittance on the microbunching instability consisting of a low-pass filter applied to the kernel for the evaluation of the collective force. Effectively, this enables the modeling of beams with finite transverse emittance by a 2D solver. The generally good agreement found with the 4D linear theory, in the regime where linear theory applies, has encouraged us to pursue application of this solver to the study of systems of practical interest.
The core of the paper (Sec. IV) is a discussion of two linac designs proposed for the Fermi@Elettra FEL [11] and their response to the microbunching instability seeded by shot noise-the random charge density fluctuations caused by the granularity of the elementary charge and the most fundamental source of the microbunching instability. Our main result is an estimate of the beam (uncorrelated) energy spread at the exit of the linac showing that at least one of the two lattices considered meets the required maximum 150 keV specification.
As for the rest of the paper, in Sec. II we summarize the main aspects of the method introduced in [10] and report on some features added to the solver, like the option to enforce periodic-boundary conditions. In the same section we also report testing of the 2D Vlasov solver by comparison against the 4D linear theory. In Sec. III we describe the model adopted to represent shot noise and, finally, in the Appendix we show how the existing 4D linear theory for the microbunching gain function was extended to include acceleration.

II. OUTLINE OF THE 2D SOLVER AND VALIDATION AGAINST LINEAR THEORY
The method proposed in [10] is based upon the replacement of the coordinate pair z; E with the pair (ẑ z,Ê E ÿ ), where the correlation function, or energy ''chirp,'' ẑ; s has been subtracted. Here E and z denote the particle energy and the longitudinal position relative to the reference orbit (z > 0 for particles in the bunch head).
The method in [10] was developed with the possibility of representing the entire longitudinal phase space of a beam in mind. This may not be needed, though, if one is exclusively interested in the study of the microbunching instability. We can substantially reduce the computation time by choosing to follow only the dynamics of particles belonging to a window of length L b 2l b , smaller than the bunch length. In practice, this can be accomplished by imposing the periodic-boundary condition inẑ. The presence of the boundaries will not alter the description of the beam dynamics significantly if L b is sufficiently larger than the characteristic length scale of the microbunching instability.
In the calculations presented in this paper, we always enforce periodic-boundary conditions. The initial beam density fẑ;Ê, is represented on a Cartesian grid with cell sizes ẑ and Ê spanning the region ÿl 0 b ; l 0 b ÿÊ 0 max ;Ê 0 max , whereÊ 0 max is chosen so that fẑ;Ê is negligible for jÊj >Ê 0 max . The normalization is such that N b fẑ;ÊẑÊ represents the number of electrons contained in the cell ẑÊ centered at ẑ;Ê, and N b is the total number of electrons contained in the region spanned by the grid.
Propagation of the density function is done by interleaving kicks under the collective force and advancements carried out under the action of the external forces. At each step, the values of the density function off-grid points are determined by local interpolation using cubic polynomials. During propagation the support of the beam is adapted to follow compression inẑ and dilation inÊ, i.e. at the current position s along the lattice the domain spanned by the grid is max , where C Cs 0 ; s is the compression factor from the start s 0 to s.

A. Mapping under collective force
The mapping for collective kicks separated by s is given by Ê 0 Ê Fẑ; ss; (1) where the collective force in frequency space is defined in terms of the impedance (per unit length)Ẑk as andk 2 ÿ1 R l b ÿl b dze ÿikz z is the Fourier integral of the normalized charge density z R dEfz; E, with normalization R l b ÿl b dzz 1, and k n 2n=L b is the wave number.
The smearing effect of a finite horizontal emittance is modeled by the insertion of a low-pass filter e ÿk 2 2" x H q , yielding a cutoff Here H x D 2 2 x DD 0 x D 0 2 is the dispersion invariant and x , x , x the familiar Twiss functions. The only sources of collective effects considered for the present study are coherent synchrotron radiation (CSR), treated in the free-space model of Refs. [10,12] (1D line charge in uniform circular motion with no account of transients), and space charge. The model adopted for the latter is given by [10,13,14] where K 1 x is the modified Bessel function and Z 0 120 , the vacuum impedance. This formula applies to a bunch with transversally uniform density and circular cross section of radius r b in free space and yields the electric field on the beam axis. To model a transversally Gaussian beam with rms sizes x and y , we set r b 1:7 x y =2, as suggested in [13]. For a discussion of the range of validity of this impedance, see [15].

B. Mapping under external force
We assume ultrarelativistic motion and neglect any slippage along the longitudinal coordinate not due to dispersive effects.
Propagation under external fields in a bending magnet from location s to location s 0 > s is provided by the mapping ẑ 0 ẑ Cs; s 0 Ê E r s dR 56 s; s 0 ; Ê 0 ÊCs; s 0 ; where Cs; s 0 1=1 hsdR 56 s; s 0 is the compression that the beam undergoes from s to s 0 , dR 56 s; s 0 R 56 s 0 ÿ R 56 s the increment of the R 56 entry of the transfer matrix describing slippage by off-momentum particles. A linear energy/position correlation is assumed, z; s zhsE r s E r s, where E r s is the design beam energy. In a bend, the linear chirp evolves according to hs 0 hs=1 hsdR 56 s; s 0 , while in an rf structure, where the energy of the reference particle varies according to 1 s hsE r s evolves as In the rf structures as well as in all other lattice elements, with the exception of bends, the dynamics in the coordi-nates ẑ;Ê in the absence of collective effects is simply the identity (ẑ 0 ẑ,Ê 0 Ê).

C. Validation against linear theory
We investigated the accuracy of our model for emittance-induced smearing of microbunching (the lowpass filter mentioned in Sec. II A) by comparing the smallamplitude gain function calculated from the numerical solution of the Vlasov equation and the 4D linear theory (see Appendix A). The gain function between locations s i and s f along the linac is defined as the ratio A f =A i =Cs i ; s f , where A i is the amplitude of a small sinusoidal perturbation at s i evolving into a perturbation of amplitude A f at s f , and Cs i ; s f is the compression factor. In [10] we reported good agreement for a section of linac encompassing a single bunch compressor. We have later verified that the agreement remains very satisfactory through longer sections spanning two bunch compressors. This is illustrated in Figs. 1-3, referring the Fermi lattices that will be discussed in detail in Sec. IV. A disagreement between numerical and linear theory becomes noticeable (but is still acceptable) only when collective effects outside the bunch compressors are turned off, as in Fig. 3. Incidentally, notice the relative unimportance of CSR in comparison to space charge ( Fig. 3 vs Fig. 1) and the considerable lower gains offered by the one-BC lattice compared to the two-BC lattice ( Fig. 2 vs Fig. 1). Also, from these pictures one can appreciate the beneficial effect of a larger beam energy spread.

III. MODEL OF SHOT NOISE
Shot noise, caused by the granularity of the elementary electron charge, is an unavoidable and the most fundamental source of undesired fluctuations. Other sources (e.g. noise in the photo-gun laser) may be significant but will not be considered in this study.
We model shot noise in the electron beam by applying a random perturbation to a smooth density function f 0 ij f 0 z i ; E j defined on the grid nodes i; j at the start of the simulation. In the calculations for this paper f 0 z; E is assumed to be uniform in z and Gaussian in E.
Let N ij be the number of electrons occupying the area zE of phase space centered on the grid node z i ; E j . We regard N ij as a stochastic process [16] obeying the Poisson statistics with average hN ij i N b f 0 ij zE and variance hN ij ÿ hN ij i 2 i 1=2 hN ij i 1=2 . If N ij is sufficiently large we can write where ij is an uncorrelated, bivariate normal stochastic process with zero average and variance equal to unity h 2 ij i ij . Dividing both sides of the above equation by N b zE yields the following prescription for the perturbed beam density: We end this section with two observations. An obvious limitation in the above model is represented by the use of a finite-density grid. While the power spectrum of shot noise is uniform, a finite grid cell size introduces a natural highfrequency cutoff. Care has to be taken to make sure that the grid density can support the physically meaningful part of the spectrum of the perturbation. Fortunately, very highfrequency components of the noise spectrum are not expected to undergo appreciable magnification because of the presence of mixing mechanisms smearing the microbunching (finite energy spread and transverse emittance). Therefore, excluding the high-frequency part of the noise spectrum should have limited consequences on the accuracy of the model.
A second observation concerns the potential problems posed to the Vlasov solver by propagation of a nonsmooth density. A key element of a solver is an interpolation scheme to reconstruct the value of the density function off-grid points after each time step. The interpolation algorithm (in our case using cubic polynomials [10]) presupposes a certain smoothness of the underlying function in order for the interpolation errors to remain bounded. Clearly the random perturbation we impose to the initial density violates this smoothness assumption and one may fear that it could result into an unacceptable error. A possible solution to this difficulty is to have the noisy density function first defined on a coarser grid and then extended to a denser grid by interpolation before starting propagation by the Vlasov solver. In our code we implemented this option but have found that, when applied, it did not affect the results significantly.

IV. APPLICATION TO FERMI@ELETTRA
We illustrate the application of our Vlasov solver to the study of two lattice designs under consideration for the Fermi project. Fermi is a tunable, soft x-ray, seeded FEL 4th generation light source presently under design at the Elettra laboratory in Trieste with a target radiation wavelength in the 10 -100 nm range [11].
Starting from a photoinjector, a 1.2 GeV linac will feed a beam with assumed E 150 keV maximum uncorrelated energy spread and transverse (normalized) rms emittance below 1:5 m to two distinct undulator lines. A 0.8 kA peak current will be needed to meet the desired brightness specifications.
Other important beam requirements (concerning, for example, the linearity of the energy chirp) are not affected by (or have any consequence on) microbunching and are not of concern in our study.
The desired peak current will be achieved by using bunch compressors to enhance the beam density by about a factor 10. At this time two alternate lattices are being considered employing one and two bunch compressors. In either case it is anticipated that a laser heater will be installed to increase the beam energy spread as a way to control the microbunching instability [13,17]. In our simulations we have yet to implement a model for the beamlaser interaction in the heater and the energy distribution of the initial beam is assumed to be Gaussian.
In the simulations, the beam is followed from the start of the main linac at about 96 MeV beam energy. The only collective effects considered are those relevant for microbunching, i.e. space charge and coherent synchrotron radiation. We propagate the beam density function through the linac with the exact account of the linear optics for determination of the transfer matrix and local beam transverse sizes (needed for evaluating the space-charge force). We assume that through the linac the beam maintains a Gaussian transverse density with constant emittances. In the simulations discussed here " x " y 1 m (normalized rms emittances).
The effect of rf wakefields is neglected. This does not represent a serious omission as the spectrum of the impedance for the rf structures is not expected to overlap significantly with the regions of the spectrum relevant for the microbunching instability. However, rf wakes affect the evolution of the beam energy chirp [18]. To compensate for the absence of the rf wakes, we adjusted the rf cavities in the lattice so as to give the beam the z=E correlation required to achieve the desired compression.

A. The two-BC Lattice
A two bunch-compressor lattice represents the current baseline for the Fermi project [19]. The section of the linac considered in our simulations is approximately 160 m long including the spreader (however, in the calculation discussed here the bends in the spreader were replaced by drifts). The two bunch compressors (which we will refer to as BC1 and BC2) are placed at locations along the linac where the beam energy is about 220 and 600 MeV, respectively. They yield an overall 10.5 compression factor partitioned about equally between the two (3.5 and 3, respectively).
Snapshots of the longitudinal phase space (uncorrelated energy deviation vs longitudinal position) at selected locations along the lattice are shown in Figs. 4 and 5 together with the normalized longitudinal density. In this example the initial uncorrelated energy spread is E0 13 keV and the final peak current I f 1 kA, 25% larger than the Fermi specifications. The three sets of pictures in Fig. 4 are taken at the start of the linac, exit of BC1, and at the entrance of BC2. The random noise present in the initial distribution can be seen to have evolved into a 2% charge density fluctuation by the exit of BC1. At this location the energy modulation is still small compared to the beam energy spread but it increases noticeably by the time the beam enters BC2 (top-right picture). A cursory inspection of this picture shows a dominant modulation somewhat below ' 10 m, consistent with a maximum gain at wavelength ' 24 m=3:5 7:1 m indicated by linear theory (see the gain function in Fig. 6). Notice that between the two bunch compressors the charge density profile remains unchanged as the beam is ''frozen'' outside the dispersive regions of the lattice.
The left phase-space picture of Fig. 5 at the exit of BC2 shows clear evidence of instability saturation (the folding of the beam density in phase space). The instability saturation leaves behind fairly large fluctuations (about 5%) in the charge density, with a frequency spectrum roughly delimited from above by the wavelength of the maximum gain predicted by linear theory ( ' 8 m). These relatively large fluctuations remain unchanged to the end of the linac, and result into additional space-charge induced energy modulation. At the exit of BC2 where the energy E r ' 600 MeV and the effective transverse radius is r b ' 200 m, the space-charge impedance (per unit length) at ' 8 m is aboutẐ ' 90 =m. A relative modulation in the charge density of A 0:05 results into an energy change of the order E keV=sm ' 10 ÿ3 AIAmpereẐ =m ' 4:5 keV=m, which projected over the remaining 60 m of linac adds up to a 270 keV energy modulation. The amplitude of the actual energy modulation gained after BC2 as determined from the simulations turns out to be a bit smaller, as the doubling of the beam energy by the end of the linac tames the effect of space charge (at the end of the linac where E r 1:2 GeV, E keV=s m ' 1:4 keV=m), but can be clearly seen in the phase space in the top-right picture of Fig. 5. At the end of the linac the uncorrelated rms energy spread averaged over z as determined from the beam 2D distribution equals about 190 keV, or 50% larger than the 136 keV 10:5 13 keV beam energy spread induced by compression alone (i.e. in the absence of collective effects). About half of this incremental energy spread is accumulated in the linac section past BC2.
Such a large energy spread at extraction would not be acceptable for Fermi. One might hope to reduce the effect of the instability by further increasing the energy spread E0 at the start of the linac with a modified tuning of the laser heater. However, because a lower bound to the energy spread at extraction is the zero-current limit E C E0 , it is clear that a compression factor C ' 10:5 leaves limited room for maneuver. Indeed, the systematic study of E vs E0 reported in Fig. 7 shows that for I f 1 kA, E 190 keV is close to the minimum attainable energy spread. Figure 7 reports the average energy spread at extraction from 10 different realizations of random perturbation to the initial beam density modeling shot noise. The error bars delimit the range of the results. At the design peak current, I f 0:8 kA, the minimum attainable energy spread is lower but still above the 150 keV target value.

B. The one-BC lattice
In an attempt to reduce microbunching, a second lattice design has been proposed [19,20] consisting of a single bunch compressor providing at once the desired factor 10 compression. The version discussed in this section is about 200 m long and includes the small chicane for laser heating (which was excluded in the lattice considered in IVA). The beam interacts with the laser in the middle of the laser-  heater chicane but in our simulations the beam is assumed to have from the start the rms energy spread resulting from the interaction with the laser. Also, in contrast to the simulations in IVA, the dipoles in the spreader are treated as actual bending elements.
The single bunch compressor of this lattice is comparable with BC1 of the two-BC lattice of Sec. IVA and provides only a slightly larger R 56 . The larger compression is therefore mostly due to a larger energy chirp. The location of the bunch compressor within the linac and the beam energy are also comparable to those of BC1.
Full compression of the beam early on in the linac causes an unfavorable enhancement of the self-fields, as they scale with the beam peak current. However, this is abundantly compensated by a larger uncorrelated energy spread experienced by the beam over a longer section of the linac.
Evidence of the better standing of the one-BC lattice is already apparent from the linear-theory gain curves of Fig. 2. The linear gain is an order of magnitude smaller than in the two-BC lattice for comparable initial energy spread (Fig. 2 vs Fig. 1).
In Fig. 2 the additional gain experienced by microbunching past the bunch compressor (the darker vs the lighter dots) is due to the dispersive region of the spreader. We should mention that care was taken to minimize the variation of R 56 through the spreader by suitable setting of the spreader dipoles. In an earlier version of this lattice, before such an optimization was made, the gain function through the end of the machine turned out to be considerably larger than in Fig. 2. (A detailed discussion of the spreader design and its impact on microbunching will be reported elsewhere.) Phase-space snapshots are shown in Fig. 8 for the case of a beam with I f 1 keV peak current and E0 7 keV initial energy spread (about half the energy spread used for the calculations of Figs. 4 and 5). The final energy spread is seen to remain below E 120 keV and the charge den- sity fluctuations at a modest 1% level. Also, observe the absence of saturation (the energy modulation in phase space remains ''upright'', top-right picture).
The generally more attractive behavior shown by the one-BC lattice is confirmed by the study of the rms uncorrelated energy spread at extraction vs the initial beam rms energy spread (Fig. 9) showing that energy spreads smaller that 110 keV would be within reach at I f 0:8 kA. The figure indicates an optimum tuning of the laser heating in the neighborhood of 7-8 keV.

V. CONCLUSION
This paper contains the first attempt to apply direct Vlasov solver methods to the study of the microbunching instability in single-pass systems.
We have discussed two lattices proposed for the Fermi project at Electra and limited our study to consideration of the microbunching instability stemming from shot noise. This is the most fundamental but not the only source of undesired charge density fluctuations, and therefore our results should be interpreted as providing, within the limitations of our model, a lower bound to quality degradation experienced by the beam. In the current baseline lattice (with two bunch compressors), we found that the microbunching instability from shot noise alone would be sufficient to cause an energy spread at extraction larger than the desired 150 keV. In contrast, the one bunch-compressor lattice we have considered would meet the required specifications and leave a comfortable buffer provided that the laser heater be tuned to generate an initial energy spread of 7-8 keV. Additional studies are needed to verify that the desired beam quality would be maintained in the transverse phase space as well.
Because of a number of approximations involved in the model (including a simplified treatment of the longitudinal space charge), some caution should be exercised at this time in assessing our results.
A validation of our solver against macroparticle simulations is in our plans. We have yet to carry out a detailed comparison. However, preliminary contacts with results from IMPACT simulations using 1B macroparticles for the Fermi two-BC lattice are quite encouraging [21]. We should point out that a meaningful comparison with simulations employing a limited number of macroparticles could be made as well provided that in our Vlasov solver the amplitude of the initial random perturbation to the beam density be adjusted to reflect the shot noise of the macroparticle distribution.
While we do not expect the level of accuracy of a 2D Vlasov solver to be the same as that of multibillion macroparticle simulations in a full 6D phase space, we should emphasize that the value of our model is in its simplicity and speed of execution [22], which should make it a useful tool for optimization and comparative studies.

ACKNOWLEDGMENTS
We would like to thank W. Fawley, R. Warnock, and, in particular, A. Zholents for many fruitful discussions. This work was supported by Department of Energy Contract No. DE-AC02-05CH11231.

APPENDIX A: LINEAR THEORY IN THE PRESENCE OF ACCELERATION
The bunching function for a coasting beam with density function fz x; p x ; z; in a 4D phase space is defined as the Fourier integral bk; s Z e ÿikz fz; sz: In linear approximation, a sinusoidal charge perturbation of amplitude bk 0 ; s 0 at s s 0 with wave number k 0 will evolve into a sinusoidal modulation of amplitude bk; s and wave number ks Cs 0 ; sk 0 at s > s 0 , where Cs 0 ; s is the compression experienced by the beam from s 0 to s. A Volterra-type integral equation for bk; s was derived in [12] (see also [23] whose notation we follow more closely) in the absence of acceleration. In this Appendix we show how with few modifications the same equation can be extended to include the more general case with acceleration. Consider the unperturbed dynamics (no collective effects) in terms of the horizontal coordinate x, the longitudinal position z, and the energy deviation E ÿ E r s=E 0 scaled with respect to the beam energy E 0 E r s 0 at s s 0 : x 00 rsx 0 k x sx Rs where rs E ÿ1 r dE r =ds is the (relative) energy gain due to acceleration and qs E cav E 0 L cav ! rf c cos s , with both rs and qs vanishing everywhere outside the rf structures.
The motion in these coordinates is non-Hamiltonian. However, it is well known that we can recover canonical equations of motion if rs 1 by introducing the new coordinatex xE 0 =E r ÿ1=2 . Upon inserting this expression into (A2) and (A3) and neglecting slow varying terms, we find Denote withR s!s 0 the transfer matrix yielding the solutions of (A4)-(A6) in terms of the coordinatesẑ x;p x dx=ds;ẑ z;.
The transfer matrix has the general form Written in terms of the new coordinatesẑ and transfer matrixR, the integral equation for the bunch function is formally the same as in [12,23] In the above expressions we used the shorthand notation R ij s R s 0 !s ij , and Cs Cs 0 ; s. The beam peak current Is and wave number ks at s scale according to Cs: ks Csk 0 and Is CsI 0 , where k 0 and I 0 are the values at the start of the linac. The relativistic factor 0 s 0 and Twiss functions x0 x s 0 , x0 x s 0 are also understood to be the values at s s 0 .
Finally, the inhomogeneous term b 0 ks; s on the righthand side of (A8) has the expression b 0 ks; s b 0 k 0 ; s 0 exp ÿ k 2 sR 2 56 s 2 0 2 exp ÿ ks 2 H s" x0 2s : (A13) The dispersion invariant H can be written as H s s 0 x0R51 s ÿ x0R52 s 2 R 2 52 s x0 :