High resolution simulation of beam dynamics in electron linacs for x-ray free electron lasers

In this paper we report on large-scale high resolution simulations of beam dynamics in electron linacs for the next-generation x-ray free electron lasers (FELs). We describe key features of a parallel macroparticle simulation code including three-dimensional (3D) space-charge effects, short-range structure wakefields, coherent synchrotron radiation (CSR) wakefields, and treatment of radio-frequency (rf) accelerating cavities using maps obtained from axial field profiles. We present a study of the microbunching instability causing severe electron beam fragmentation in the longitudinal phase space which is a critical issue for future FELs. Using parameters for a proposed FEL linac at Lawrence Berkeley National Laboratory (LBNL), we show that a large number of macroparticles (beyond 100 million) is generally needed to control the numerical macroparticle shot noise and avoid overestimating the microbunching instability. We explore the effect of the longitudinal grid on simulation results. We also study the effect of initial uncorrelated energy spread on the final uncorrelated energy spread of the beam for the FEL linac.


I. INTRODUCTION
Generation of high-power coherent radiation in x-ray free electron lasers (FELs) requires high peak current, low emittance, and low energy spread multi-GeV electron beams.Typically, low peak current, low emittance, and low energy spread beams are produced in the electron gun while a high peak current is obtained through various stages of bunch compression during acceleration in the linac to a final energy.Unfortunately, collective effects driven by space charge, linac structure wakefields, and coherent synchrotron radiation (CSR) tend to cause an unacceptable degradation of the beam emittance and energy spread thus limiting the maximum achievable peak current.Previous studies have shown [1][2][3][4][5] that the accumulated effect of the longitudinal space-charge force along the linac driven by small density fluctuations can cause large energy modulations along the bunch.After passing through a bunch compressor, the energy modulation is further amplified in the presence of CSR and is converted into larger density modulations.This phenomenon, known as ''microbunching instability,'' results in the appearance of large current fluctuations and significant growth of the uncorrelated and correlated energy spread, where the latter is seen as electron beam fragmentation in the longitudinal phase space.Often, the gain of the instability is large enough to amplify to unacceptable level the inevitable density fluctuations present in the electron beam due to shot noise [6].
The microbunching instability has been studied using both a direct Vlasov solver [7,8] and macroparticle tracking [9].Direct Vlasov solvers do not suffer from the spurious noise characteristic of macroparticle sampling and may provide an accurate description of the microbunching instability both in the linear and nonlinear regimes.However, in the implementation presently available this method is limited to a simplified two-dimensional phase-space computational model.
Macroparticle codes provide well established and successful methods for studying beam dynamics but they tend to overestimate the effect of instabilities that are sensitive to small fluctuations of the beam density, like those induced by shot noise.These will be artificially magnified by a factor ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N=N mp q in a simulation employing N mp macroparticles to represent a beam of N electrons.Therefore, an obvious way to improve the accuracy of the simulations is to increase N mp .Fortunately, progress in code development and the availability of high-power computational resources are now making simulations with billions of macroparticles (a number comparable to a real electron bunch population) quite practical.We should add that important aspects of the microbunching instability, like the determination of the gain function, can still be addressed in a meaningful way using a limited number of macroparticles, as demonstrated in [9].Moreover, it is possible that quietstart techniques to reduce the unphysical component of shot noise, which have proved quite successful for the studies of beam dynamics in FELs [10], could also be applied to the beam dynamics in a linac.However, our preliminary attempts to experiment with these techniques did not yield satisfactory results and eventually motivated us to pursue the more straightforward strategy of increasing the number of macroparticles used in the simulations.Another advantage of using a large number of macroparticles (close to real number of electrons in a beam) is that PHYSICAL REVIEW SPECIAL TOPICS -ACCELERATORS AND BEAMS 12, 100702 (2009) 1098-4402=09=12(10)=100702 (11) 100702-1 Ó 2009 The American Physical Society the macroparticle distribution at the end of the linac can be directly used for next stage FEL simulation without worrying about proper modeling of initial shot noise of the beam.One purpose of the present paper is to report on the latest developments of the IMPACT code [11], which has enabled such simulations and is now being used in design studies of a proposed x-rays FEL linac at the Lawrence Berkeley National Laboratory [12].In this paper we also present a discussion for the need of large-scale high resolution simulations of high brightness beams, and show specific applications to beam dynamics studies for the next-generation x-ray FEL linac.The organization of this paper is as follows: the computational model used in the beam dynamics simulation is presented in Sec.II; the choice of numerical parameters is discussed in Sec.III; large-scale simulation of the proposed LBNL FEL linac is presented in Sec.IV; and the summary is given in Sec.V. Finally, a macroparticle up-sampling scheme that reduces the macroparticle shot noise from an initial distribution with a smaller number of macroparticles while maintaining the global properties of the original distribution is described in the Appendix.

II. COMPUTATIONAL MODEL
The beam dynamics simulation of the FEL linac was performed using a parallel macroparticle tracking code, IMPACT [11].The IMPACT code is an object-based parallel particle-in-cell code to simulate high intensity, high brightness beam transport mainly in rf linear accelerators.It uses a split-operator method to separate the particle motion due to the given external fields from the particle motion due to the collective self-consistent space-charge fields and wakefields [13].

A. Modeling of rf structures
Beam dynamics inside an rf cavity, including longitudinal accelerating fields and transverse focusing fields, is modeled using linear transfer maps and multiple reference particles to accommodate nonlinear rf effects along the bunch discussed at the end of this section.The vector potentials associated with an axis-symmetric rf cavity, specified through the order needed to compute linear transverse dynamics, are where "ðzÞ denotes the spatial part of the electric field at r ¼ 0, and where a superscript prime denotes d=dz, ! is the angular frequency of the rf field, c is the speed of light in vacuum, and is the initial phase.Using the longitudinal position z as the independent variable, the dynamics can be described in terms of canonical coordinates and momenta ðx; p x ; y; p y ; t; p t Þ and the Hamiltonian is given by Hðx; p x ; y; p y ; t; p t ; zÞ where t denotes the particle's arrival time at the location z, its canonically conjugate momentum p t is the negative of the particle's total energy, and q is the charge of the particle.One can next define a set of canonical variables that are dimensionless deviations from a reference trajectory [x ¼ p x ¼ y ¼ p y ¼ 0, t 0 ðzÞ, p t0 ðzÞ] and one can obtain the corresponding Hamiltonian for these new variables.Expanding the new Hamiltonian around the reference trajectory, the lowest order terms are quadratic, and they govern the linear dynamics.The linear transfer matrix corresponding to the quadratic Hamiltonian can furthermore be factored as [14] " where l ¼ c=!, ! is the scaling angular frequency, m is the rest mass of a particle, and superscripts i and f denote initial and final values, respectively.In the preceding formulas, Á is given by and z0 is velocity of the reference particle on axis, and s is the synchronous phase of the reference particle given by s ¼ !t 0 ðzÞ þ .Here, the time of flight t 0 and the relativistic factor 0 of the reference particle are given by The quantities a j , b j , c j , and d j (j ¼ x, y, ) satisfy the following equations: with initial values a j ¼ d j ¼ 1, b j ¼ c j ¼ 0, and j (j ¼ x, y, ) are given by As a test example, we transport a test particle through an rf cavity using the above linear transfer map and compare it with direct numerical integration of the Hamiltonian equations of motion.The computed trajectories from both methods are shown in Fig. 1.The transfer map method agrees well with direct numerical integration in this example.Using above transfer maps for a single reference particle, we slice the beam longitudinally into multiple slices so that each slice has an individual reference particle through the rf cavity.The particle coordinates are converted into the multiple reference particle system coordinates in front of the rf cavity, tracking through the rf cavity, and then converted back to the original coordinates.Using multiple reference particles helps capture the effect of nonlinear rf acceleration.

B. Treatment of space-charge effects
The 3D space-charge forces are calculated at each time step from the solution of the Poisson equation in the beam frame using a convolution of the electron charge density with the Green function for the scalar potential with open boundary conditions.This convolution is calculated numerically on a 3D grid using an integrated Green function method with an fast-Fourier-transform (FFT)-based calculation of the cyclic summation on a doubled computational domain [15].The space-charge fields are Lorentz transformed back to the laboratory frame and used along with the corresponding azimuthal magnetic self-field to advance the particle momentum.
As a test of the space-charge model, we computed the energy modulation of an initial 120 MeV round uniform electron beam with r b ¼ 200 m radius and I ¼ 120 A uniform current perturbed by a A ¼ 5% modulation, and zero initial emittance propagating through three-meter drift space.Figure 2 shows the amplitude of energy modulation as a function of distance in comparison with an often-used analytical model of longitudinal space-charge impedance [6], where Z 0 ' 120 Ohms is the vacuum impedance, K 1 the modified Bessel function of the second kind, the relativistic factor, and k ¼ 2= the perturbation wave number.This analytical model presupposes that the longitudinal component of the electric field across the beam is uniform and equal to the value on the beam axis.This is a good approximation if the wavelength of the perturbation as measured in the beam comoving frame is large compared to the beam transverse radius (i.e.kr b = ( 1).The relative energy modulation induced over the distance s by the model impedance (19) starting from a relative current sinusoidal modulation A is obtained from the formula ÁE=E ¼ 4A jZðkÞj Z 0 ðsIÞ=ðI A Þ, where I A ' 17 kA is the Alfve ´n current.However, as shown in the figure, at smaller wavelengths this analytical model (dashed lines) tends to overestimate the energy modulation when this is averaged over the beam transverse density [16].Indeed, a better approximation of the numerical result (see solid lines in the figure) is given by the average yielding where I 0 and I 1 are the modified Bessel functions of the first kind.Note that Eq. ( 19) corresponds to the on-axis impedance ZðkÞ ¼ Zðk; r ¼ 0Þ.

C. Computation of CSR wakefields
The CSR is modeled as a one-dimensional (1D) longitudinal wakefield as described in Refs.[17][18][19], where the change of energy due to CSR is given for steady state radiation by The above integration is separated into two integrals: one is from À1 to s À h, and the other is from s À h to s [20].Here, h is the mesh size for numerical integration.To calculate the first integral numerically, we carry out explicit low-pass filtering of the charge density and utilize a low order trapezoidal rule numerical integration scheme to s À h.The numerical filter used here is given by This local filter is computationally inexpensive to apply and also conserves total charge.Figure 3 shows the transfer function of the filter in the frequency domain.The transfer function of the filter approaches zero for frequency close to the Nyquist frequency.In the high statistic simulation such as using billion macroparticles, this filter can be turned off without affecting the final simulation results.For the second integral within the last mesh spacing, h, of the singularity, we replace ðsÞ with a second-order polynomial fit to the smoothed density, and then evaluate the integral over ½s À h; s analytically.Figure 4 shows the results of testing the CSR routine by simulation of the passage through a single 50 cm bend magnet of a distribution composed of two million particles all having equal energy, but with a sinusoidal charge density modulation in the z direction.The modulation wavelength is 20 microns with 20% amplitude.The beam has 0.8 nC charge with 230 MeV energy.The bunch length is 0.2 mm and the bending radius is 7.11 m.The solid blue line is the theoretical prediction, seen to be in very good agreement with the numerical simulation result, shown in red.Here, the theoretical prediction is based on direct analytical integration of Eq. ( 22) with a sinusoidal current modulation on top of a flattop distribution.In separate tests we verified that the local energy spread in the final distribution is due to a nonzero transverse emittance (0.015 mm mrad) of the initial distribution.What appears to be a solid horizontal line at Á ¼ 0 is the initial distribution footprint in the longitudinal phase space.

D. Computation of structure wakefields
Besides the CSR wakefield, the short-range longitudinal wakefield (monopole) and transverse wakefield (dipole) are also included in the computational model of the IMPACT code.The effective forces from wakefields are given by where W T is the transverse wake function in the time domain, s max is the bunch head location, W L is the longitudinal wake function in the time domain, and is the line density function of the beam.For the effective force in the transverse y direction, x is replaced by y in Eq. ( 24) and the corresponding transverse wake function is used.The transverse and longitudinal wake functions are provided externally according to different accelerating structures.To compute the effective forces from wakefields more efficiently, we extend the above integrals into a full domain convolution as where s min is the bunch tail location and ðsÞ ¼ xðsÞðsÞ for transverse wake ðsÞ for longitudinal wake: (28) The above convolution has the same form as the spacecharge field computing in an open boundary condition and can be calculated using the same FFT-based method on a doubled computational domain as the space-charge force calculation [15].This reduces the computational cost from the OðN 2 Þ by direct summation of the above wake function integral to the O½N logðNÞ by using the FFT-based method.

III. CHOICE OF NUMERICAL PARAMETERS
When carrying out macroparticle simulations of the beam dynamics in an FEL linac, the choice of numerical parameters, e.g., number of macroparticles, can significantly affect the final outcome because of the high sensitivity of the microbunching instability to small density fluctuations.In this study we modeled two cases: one with 2 keV uncorrelated energy spread at the beginning of the linac and one with 5 keV.The 2 keV case is expected to exhibit more microbunching instability [6].In what follows we show a few simulation results of the electron beam dynamics in the proposed FEL linac at the Lawrence Berkeley National Laboratory using an initial square root parabolic current distribution.The only parameter that we varied for each case was the number of macroparticles.Specific information regarding linac configuration and electron beam parameters can be found in [12] and will be discussed in Sec.IV.We characterize the microbunching instability using two parameters: the rms slice energy spread and the rms fluctuation of the energy centroid of each slice.The slices were chosen to be comparable to the size of the grid used in the Poisson solver, i.e., of the order of a fraction of a m.The rms energy spread for each ith slice is defined as where hÁi i denotes averaging over the electrons in the ith slice.The rms fluctuation of the slice energy centroid is defined as 2 flct ¼ hðhEi i À " EÞ 2 i, where " E is the average beam energy obtained after subtraction of long-scale smooth energy variations caused by the nonlinearity of the rf waveform and rf structure wakefields.Similarly to i , we exclude tails while calculating flct .The fluctuation of the slice energy centroid at short wavelength may interfere with the possible application of seeding techniques for x-ray production and induce an undesired increase in the bandwidth of the resulting x-ray pulses [21].Figure 5 shows the slice (i.e.uncorrelated) energy spread averaged over all slices in the bunch core at the end of the linac as a function of the number of macroparticles.The error bars in the figure give a measure of the variation of i from slice to slice, i.e., rms value of the slice energy spread variation.For the case with 5 keV initial energy spread (red data points), we observe that the uncorrelated energy spread at the end of the linac begins to saturate at 100 million macroparticles, while it continues to decrease in the 2 keV case (green data points).
Figure 6 shows the final rms energy fluctuation flct as a function of the number of macroparticles employed in the simulation.For the initial 5 keV uncorrelated energy spread case, the amplitude of the energy fluctuations ap-proaches a saturation value of about 120 keV when the number macroparticles exceeds one billion.For the initial 2 keV uncorrelated energy spread, the amplitude of the energy fluctuations continues to decrease up to the largest number of macroparticle used, five billion, which is the same as the physical number of electrons contained in the 0.8 nC electron beam.
Figures 5 and 6 suggest that an adequate number of macroparticles is needed to obtain a generally fair assessment of the energy spread at the end of the linac in the presence of microbunching instability starting from the shot noise.The exact number of macroparticles to be used may vary from case to case and will depend on the gain of the microbunching instability for the specific linac lattice under consideration.
The numerical grid used in the calculation of collective effects can provide artificial smoothing of the numerical electron density distribution compared with the real distribution.As the number of macroparticles approaches the physical number of electrons in a bunch, the finite number of numerical grid points being used may cause some unphysical smoothing of the real shot noise inside electron beam.For example, using 2048 numerical grid points in our application supports a minimum wavelength of 3 m.Any noise with a wavelength below that wavelength will be suppressed by the grid.Fortunately, the gain of the microbunching instability falls off quickly at short wavelengths due to the smoothing effect of finite energy spread.Figure 7 shows the gain function of the microbunching instability using the linear theory [4] together with the numerical gain computed using macroparticle simulation [22].It is seen that the gain of the microbunching instability has a peak around 200 m and decreases quickly towards short wavelength.Further increasing the number  The gain curve for the microbunching instability through a portion of beam line considered in our studies for an FEL linac at LBNL (see text).Linear theory (solid line) is compared to results from macroparticle simulations (dots).The reported value of is for the modulation wavelength before compression. of grid points beyond 2048 will not change the initial shot noise level in the current profile around low frequency where the gain curve is peaked, and hence the final energy spread driven by the microbunching instability.Incidentally, the above figure also shows a fairly good agreement between the numerically calculated gain and that from the linear theory.The slight discrepancies between the two functions could be due to the approximation involved in the space-charge impedance model used in the linear theory.To test the effects of the numerical grid, we also carried out self-consistent one billion macroparticle simulations of electron transport through the linac using two sets of numerical grid points.Figure 8 shows the final uncorrelated energy spread with 2048 and 4096 grid points.It is seen that both sets of grid points give nearly the same final energy spread.The above numerical example used an initial 5 keV uncorrelated energy spread.For a lower initial uncorrelated energy spread, the peak of the microbunching instability gain could move towards shorter wavelength, i.e., higher frequency.In this case, a larger number of numerical grid points would be needed to resolve the shorter wavelength in order to capture the physics of the microbunching instability.

IV. HIGH RESOLUTION SIMULATION OF A PROPOSED FEL LINAC AT LBNL
A design effort is now underway for a proposed soft xray FEL array at LBNL.We have carried out selfconsistent macroparticle tracking of the driver linac starting from the end of the injector to the end of the electron beam switch yard, right before the FEL, using the IMPACT code with one billion macroparticles and an initial square root parabolic current distribution.Figure 9 shows a schematic plot of the linac [12].It consists of a laser heater, two superconducting rf accelerating linac sections, a longitudinal phase-space linearizer, a bunch compressor, and an electron beam switch yard (spreader).While the real machine includes a laser heater, the interaction of the laser light with electrons was not simulated in this study.The accelerating structure used in this study is based on a 1.3 GHz nine-cell superconducting rf cavity since this linac is designed for 1 MHz electron beam repetition rate.The energy of the electron beam at the entrance to linac is about 40 MeV.The initial peak current is about 70 A with a total charge of 0.8 nC.The initial transverse rms beam size is 0.42 mm with a normalized rms emittance of 0.75 mm mrad.The electron beam is accelerated to 260 MeV at the end of the first linac (linac1) with an energy gain of 13:5 MeV=m within the accelerator cavities.It is followed by the harmonic linearizer (HL) which is the second superconducting linac producing a maximum accelerating gradient of 5 MeV=m at 3.9 GHz.The HL decelerates the electron beam slightly while (together with L1) linearizing the energy chirp of the bunch in front of the bunch compressor (BC).The BC compresses the electron bunches and increases the electron peak current up to about 1 kA.A single chicane bunch compressor is used in this design.Previous studies of the FERMI@ELETTRA FEL linac suggested that a single chicane could give a smaller final uncorrelated energy spread than the double chicane design [8,23,24].The bunch compressor consists of four rectangular bending magnets of length 0.25 meters.The bending angle for each magnet is 0.12 degrees.The separation between the first magnet and the second magnet is 4.5 meters.The separation between the second magnet and the third magnet is 1.14 meters.After the bunch compressor, the electron beam is accelerated from about 250 MeV to a final energy of 2.4 GeV before entering a beam switch yard, i.e., a spreader.The spreader includes ten branch beam line lattices.Each lattice has two parts, the beam take-off section and the FEL fan distribution section.Each part is built as a triple bend achromat.Here, a kicker, a septum and off-set quadrupoles are used together in place of one bending magnet in the beam take-off section.These achromats are also made isochronous for the nominal case in order to minimize the effects of the microbunching instability [25].A more detailed description of the spreader design is given in [26].The uncorrelated energy spread in an electron beam can smear out coherent growth of the microbunching instability.A larger initial energy spread leads to less growth of the microbunching instability.Figure 10 shows the final uncorrelated rms energy spread Figure 11 shows the longitudinal phase-space distribution at the end of the linac from a simulation using one billion macroparticles (left) and five billion macroparticles (right) with 5 keV uncorrelated energy spread, and 0.75 mm mrad transverse emittance.It is seen that the full width of the final relative energy spread is on the order 10 À4 for a major part of beam (the rms relative energy spread is 4 Â 10 À5 ).Such a small energy spread will be especially important for seeded FEL utilizing a high gain harmonic generation technique [21].There is a small energy modulation in the longitudinal phase space with a period of around 15 m.After accounting for the compression factor, this modulation period is in a good agreement with the gain function of the microbunching instability from the linear theory.This modulation comes from the initial shot noise in a 0.8 nC beam sampled with one billion macroparticles, and is further amplified through the linac due to the microbunching instability.With an increasing number of macroparticles, the initial shot noise is reduced.The final phase-space modulation as seen in the right plot of Fig. 11 using five billion macroparticles becomes smaller than that in the case of one billion macroparticles.
Figure 12 shows the current profile at the end of the linac.The peak of the current profile reaches about 1.2 kA with a reasonably flat distribution along a major part of the beam.
The simulations shown above were carried out on a Cray-XT4 parallel computer, Franklin, at the National Energy Research Scientific Computing Center.Here, we have used 512 processors (cores) for one billion macroparticle simulations.Each processor has a theoretical peak performance of 9:2 GFlop= sec and 2 GB of memory.The computing node is connected to a dedicated SeaStar2 router through HYPERTRANSPORT with a 3D torus topology to ensure high performance, low-latency communication for message passing interface jobs [27].The computing time for such a simulation is about two hours.The total memory usage of the simulation is about 112 GB.The rms  information of the beam is calculated internally during the process of simulation.The slice information such as slice emittance or energy spread and current profile are calculated internally at a given location using the macroparticle distribution.SNAPSHOT is used to output a fraction of randomly sampled macroparticles at a given location.The phase-space plots are made after projecting onto a two-dimensional grid and plotted using software such as MATLAB.

V. SUMMARY AND DISCUSSIONS
In this study we have shown that large-scale, high resolution simulation is needed to accurately model electron beam transport in the presence of the microbunching instability when designing a linac for future FEL applications.We have presented the latest developments of the IMPACT code that have made large-scale macroparticle simulation possible.Numerical parameters such as number of macroparticles and number of computational grid points need to be chosen carefully to simulate initial shot noise in the electron bunch and to capture the detailed physics of the microbunching instability.As an application, we performed simulations with one billion macroparticles of electron beam transport through a designed soft x-ray FEL linac at LBNL.The effects of initial uncorrelated energy spread were also studied.Using 5 keV initial energy spread resulted in the lowest level of energy spread at the end of the linac.
The microbunching instability is primarily driven by longitudinal space charge that is responsible for the growth of the uncorrelated energy spread at the end of the linac [6].In this study, we focused on the longitudinal space-charge effect and turned off the transverse space-charge effects inside the code.Including the transverse space-charge effects led to rms mismatch near the entrance of the linac and more than 30% emittance growth.Optimization of the linac lattice design including 3D space-charge effects is still in progress and will be reported elsewhere.At present, billion macroparticle beam dynamics simulations of the proposed FEL linac at LBNL have already shown that a reasonable beam quality can be achieved at the end of the linac with about 100 keV energy spread at 2.4 GeV energy, and 1.2 kA peak current.
In this study, we also carried out billion particle simulations starting from distributions with a smaller number of particles obtained from different codes and using the upsampling scheme described in the Appendix.The upsampling scheme is based on local particle phase-space repopulation.It produces reasonable smoothness to the original macroparticle distribution while preserving the global current profile and energy-phase correlations.There are a number of external parameters provided by the code users to control the extent of smoothness and the global structure of the up-sampled distribution.Currently, the choice of those parameters is based on a trial-and-error method.This up-sampling scheme is still under investigation for further improvement.

ACKNOWLEDGMENTS
We would like to thank W.

APPENDIX: MACROPARTICLE UP-SAMPLING OF INITIAL DISTRIBUTION
The FEL linac simulation starts with an initial electron beam coming out of a photoinjector.Often, this initial macroparticle distribution is produced with a low number of macroparticles, i.e., as small as a few hundred thousand or a few million.For an electron beam with about 1 nC charge (i.e.6.25 billion electrons), this small number of macroparticles will cause significantly larger numerical shot noise than the real shot noise among electrons.In order to reduce the numerical noise in the initial macroparticle distribution, a much larger number of macroparticles, e.g., one billion macroparticles, is used in the IMPACT code to repopulate the original particle distribution (this process is also called up-sampling).In this study, we used a local smoothing scheme instead of a global smoothing scheme in which the original distribution is fitted using some known functions and a large number of macroparticles are sampled from the new distribution function.Using a local smoothing scheme tries to preserve the structure in the original distribution while providing sufficient smoothing.However, the choice of the smoothing parameters (box size) described in the following is still based on a trial-and-error method.In the local smoothing scheme, a six-dimensional box centered at the original macroparticle phase-space location is used to up-sample the initial distribution.Inside each box, N new =N old macroparticles are generated in the sixdimensional phase space.Here, N new is the total number of macroparticles for the up-sampled distribution and N old is the total number of macroparticles in the original distribution.A uniform distribution is assumed inside the fourdimensional transverse phase space of the box.The new sampled particle coordinates in the x À p x phase space will be where " x old and " p old x are the original particle coordinates, " x new and " p new x are the new up-sampled particle coordinates, s x and s px are box sizes in the x and p x dimensions, and R x i and R px i are uniformly distributed random numbers between À1 and 1 with i ¼ 1; Á Á Á ; N new =N old .The same equation as above holds for the y-py phase space with the label x replaced by y.
The choice of the box size determines the smoothness and the structure of the new up-sampled distribution.If the box sizes are equal to zero, the new particles will be at the same phase-space location as the original particles, and there is no smoothing at all in the new up-sampled distribution and the full structure of the original distribution is maintained.On the other hand, if the box size is too large, the new distribution will have different phase-space structures compared with the original distribution due to the assumption of uniform distribution inside the box.We normally choose the box size so that the rms information such as emittances calculated from the up-sampled distribution is within a few percent of that from the original distribution.
In the longitudinal phase plane of the box, a linear current density distribution is assumed within the box based on the current profile from the original particle distribution.The correlated normalized energy deviation with respect to phase is obtained from the original distribution using a cubic spline.Here, the correlated energy deviation denotes the functional dependency of the normalized energy deviation with respect to the bunch length (phase) in the longitudinal phase space.To find the correlated energy deviation function " p corr t ðÞ, the original particle distribution is divided into a number of slices along the phase coordinate.The mean energy deviation is calculated for each slice.The functional dependence of the energy deviation with respect to phase is obtained using the cubic spline of the mean energy deviation value of each slice.The uncorrelated energy denotes the difference between the energy deviation of an individual particle and the correlated energy deviation at a given bunch phase location.The phase coordinate of each new particle, i , is generated from a uniform distribution centered around the original particle phase coordinate in the same way as the transverse coordinates.The normalized energy deviation " p new t of each new particle inside the box is given by where " p corr t ð i Þ is the correlated energy deviation at phase location i , p t is the uncorrelated energy deviation computed from the original particle centered at each box, and E i is a new uncorrelated energy deviation generated from the sampling of a Gaussian distribution with an rms amplitude provided by users.Often E is used to study the effects of initial uncorrelated energy spread on the microbunching instability and to define the optimal acceptable initial energy spread in the design.This term is set to zero if no artificial heating is included in the simulation.The box size of the up-sampling is controlled by the code user and should be chosen to produce smoothness of the resulting particle distribution representative of the real number of electrons while maintaining similar global properties (e.g.current profile) of the original distribution.If the box size is too large, the up-sampled distribution could be different from the original distribution due to the uniform or linear density assumption inside the box.On the other hand, if the box size is too small, the up-sampled macroparticles will be localized and will not provide sufficient smoothness to the original distribution.Figure 13 shows the initial current profile fluctuation (I=I 0 ) of a Gaussian distribution with two million particle direct sampling, from one billion particle up-sampling with 20 m longitudinal box size, with 300 m longitudinal box size, with 500 m longitudinal box size, and from one billion particle direct sampling.The direct two million particle sampling has a much larger current fluctuation as expected.This fluctuation is brought down with larger up-sampling box size.However, as the box size becomes too large (e.g.500 m), the upsampled current profile starts becoming noticeably different from the original distribution.From our experience, it is reasonable to choose the longitudinal box size to be large enough to cover the wavelength corresponding to the peak of the microbunching instability gain curve.This reduces the level of microbunching instability from the initial shot noise of the particle distribution.

FIG. 1 .
FIG. 1. (Color) Horizontal and longitudinal trajectories through an rf cavity from the linear transfer map and from direct numerical integration.
FIG. 3. (Color)The transfer function of the filter in frequency domain.

FIG. 4 .
FIG. 4. (Color) Longitudinal phase-space distribution after a bending magnet for an initial sinusoidal density modulation.Numerical result (red), analytical prediction (blue), and the initial distribution (green) are shown.

FIG. 6 .
FIG. 6. (Color) The rms energy fluctuation at the exit of the linac as a function of the number of macroparticles used in the simulations.Two examples of beams with 5 and 2 keV initial uncorrelated energy spread are compared.
FIG. 7. (Color)The gain curve for the microbunching instability through a portion of beam line considered in our studies for an FEL linac at LBNL (see text).Linear theory (solid line) is compared to results from macroparticle simulations (dots).The reported value of is for the modulation wavelength before compression.

FIG. 5 .
FIG. 5. (Color) Slice (or uncorrelated) energy spread at the exit of the linac as a function of the number of macroparticles used in the simulations.The reported slice energy spread is averaged over all the slices in the beam core.Two examples with 5 and 2 keV initial uncorrelated energy spread are compared.

FIG. 11
FIG. 11. (Color) Longitudinal phase space at the end of the linac using one billion (left) and five billion macroparticles (right) with an initial 5 keV rms uncorrelated energy spread.
Fawley and G. Penn for useful discussions.This research was supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.This research used resources of the National Energy Research Scientific Computing Center.

FIG. 12 .
FIG. 12. (Color) Current profile at the end of the linac with initial 5 keV rms uncorrelated energy spread.
FIG. 13. (Color) Initial current profile fluctuation of a Gaussian distribution from two million particle direct sampling, from one billion particle up-sampling with 20 m longitudinal box size, with 300 m longitudinal box size, with 500 m longitudinal box size, and from one billion particle direct sampling.
FIG.2.(Color) Amplitude of energy modulations as a function of distance induced by space charge starting from initial 5% current sinusoidal modulations with wavelengths ¼ 15 m, 30m, and 50 m for a beam drifting in free space with initial uniform transverse density, circular cross section, and vanishing emittance.The results from IMPACT simulations (dots) are compared to the longitudinal space-charge impedance models of Eq. (19), dashed lines, and Eq.(21), solid lines.