Algorithm for loading shot noise microbunching in multidimensional, free-electron laser simulation codes

We discuss the underlying reasoning behind and the details of the numerical algorithm used in the GINGER free-electron laser simulation code to load the initial shot noise microbunching on the electron beam. In particular, we point out that there are some additional subtleties which must be followed for multidimensional codes which are not necessary for one-dimensional formulations. Moreover, requiring that the higher harmonics of the microbunching also be properly initialized with the correct statistics leads to additional complexities. We present some numerical results including the predicted incoherent, spontaneous emission as tests of the shot noise algorithm’s correctness


I. INTRODUCTION AND MOTIVATION
In the past decade, there has been increasing interest worldwide in developing extremely short wavelength (e.g., l s # 1 nm), high instantaneous power (i.e., P $ 1 GW), free-electron laser (FEL) sources. Because of the lack of survivable, high reflectivity mirrors and available coherent "seed" sources at these wavelengths, the most popular approach to extreme ultraviolet (XUV) and x-ray wavelength FEL development is a single-pass amplifier configuration based upon the principle of self-amplified spontaneous emission (SASE). The initial radiation signal in a SASE FEL is incoherent spontaneous emission which originates with the temporally random microbunching present on the electron beam as it enters the magnetic undulator. That portion of the emission which falls within the spectral gain bandpass of the FEL and which also couples well to the fastest growing transverse mode then exponentially amplifies in the remaining portion of the undulator (or until saturation occurs).
A key issue in the numerical modeling of SASE FEL's is the proper simulation of the initial start-up signal and its underlying electron beam microbunching. The prodigal increase in computing power (accompanied by the nearly miraculous decrease in memory prices) has made it practical in some circumstances to consider modeling every single real electron in a few slippage length section of the electron beam. Nonetheless, it still remains most cost efficient for the great majority of investigations to use ϳ10 3 10 6 macroparticles in total, each of which represents hundreds to millions or more of real electrons. In this case, numerical codes must artificially construct macroparticle phase space distributions whose statistical properties mimic those of the actual electron beam as it enters the undulator.
This issue is not a new one for the FEL community. Probably the most widely known work on simulation of FEL start-up noise is that of Penman and McNeil [1] who discuss a particular algorithm used in their 1D code. Preceding this article, there had been work in the early 1980's by others such as Sprangle et al. [2] and Benson et al. [3,4] applicable to oscillator start-up, and independent, unpublished work in the mid-to-late 1980's by both the present author while at LLNL and by Bourianoff, Wong, and Rosenbluth at Austin Research Associates oriented toward multidimensional modeling of sideband instability start-up noise in extremely high power, single-pass FEL amplifiers. In part because multidimensional codes such as GINGER [5] and GENESIS [6] are in active use modeling ongoing SASE FEL experiments (e.g., LEUTL [7] and VISA [8]) and in part because the GINGER shot noise algorithm differs in some significant ways from that of Penman and Mc-Neil, the present author felt that it was timely to present its details.
There are two particular areas in which the GINGER shot noise algorithm diverges from the Penman-McNeil treatment. First, in a pure 1D formulation (i.e., one that ignores instantaneous energy spread and transverse emittance), the individual macroparticles (and the corresponding real electrons they represent) are essentially indistinguishable from each other, apart from their longitudinal coordinates. Consequently, the statistics of the microbunching are described by the entire number N of electrons in the some chosen temporal duration T [i.e., see Eq. (3) of [1] ]. In a multidimensional formulation (GINGER is a r 2 z axisymmetric code in space but effectively follows all three particle momentum components via p x , p y , g), the macroparticles can be arranged into a physically distinguishable collection of "beamlets" in five-space for which the relevant number N b is therefore much smaller than the total number of electrons N. This issue of distinguishability echoes a point made by Benson et al. who point out for a classical electron beam, the individual electrons comprise "a random but deterministic system." Consequently, as explained in detail in the following section, the GINGER macroparticles are subdivided into a discrete number of beamlets, each one of which has its own, statistically independent, longitudinal microbunching distribution spread among two or more individual macroparticle members.
The other area in which the current GINGER shot noise algorithm differs from previous treatments is that serious attention is paid to properly reproducing the statistics of microbunching at higher harmonics of the fundamental wavelength l s . For high gain FEL's, it has long been known that together with exponentially growing (with z) bunching at the fundamental resonant wavelength, a nonlinear effect leads to the microbunching at the higher harmonics also growing exponentially with a growth rate linearly proportional to harmonic number m. More recently, there has been interest in exploiting this phenomenon as a means to extract significant radiation power at wavelengths much shorter than l s and/or using multiple undulators in a harmonic "cascade" configuration to reach a very short output wavelength. To model such devices accurately, one must ensure that the code sets the initial microbunching (which effectively is a nuisance reducing the longitudinal coherence of the nonlinear harmonic bunching mechanism) at the appropriate level and with the correct statistics.
Mainly through the experience of trial and error, the author found that the original formulation in GINGER [9] suffered two basic problems with regard to higher harmonic microbunching. First, for any harmonic number m, one must ensure that in the absence of shot noise, the initial bunching through harmonic number m 1 1 must be identically zero because coupling in the exponential growth regime between radiation at the fundamental wavelength (m 1) and bunching at harmonic m M $ 1 leads to exponentially growing bunching at harmonics M 6 1. This implies that properly modeling the z behavior of the third harmonic requires that there be no unintentional bunching at the second and fourth harmonics. GINGER originally used sets of simple paired macroparticles (i.e., each particle at longitudinal phase u had an "identical twin" at phase u 1 p). Although this ensured no initial bunching at the fundamental wavelength (and all odd harmonics), the same is not necessarily true at even harmonics.
Furthermore, when either the instantaneous energy spread s g is significant and/or individual electron betatron motion (i.e., due to transverse emittance) introduces an important spread in p z , it is not sufficient to ensure that the global average of bunching at each harmonic m is initially zero. Rather, one must ensure locally in 5D phase space that the initial bunching is also zero. Otherwise, even in the absence of any shot noise, gain, or FEL radiation, one can find both locally and globally that the longitudinal microbunching can show oscillations with z arising from "interference" effects (due to unintended correlations) of macroparticles with different parallel velocities. This problem would arise even in 1D problems for which there was a significant spread in p k . A more elegant description of this difficulty may be found in Reiche et al. [10].
To prevent such problems, each beamlet in GINGER is now composed of 2M macroparticles (where M is the highest harmonic for which zero bunching is sought), each of which has the same initial five-space coordinate ͑x, y, p x , p y , g͒ but whose initial longitudinal coordinate u j ϵ 2pct j ͞l s sequentially differs by p͞M. As explained in the next section, when shot noise effects are desired, one must add to each macroparticle's "quiet start" u j a small du whose statistics are correct for each harmonic number 1 # m # M. Moreover, the du are completely uncorrelated, both between different beamlets and (when decomposed harmonically) between different harmonic numbers in a given beamlet. For electron beams whose distribution in five-space cannot be reduced (at least conceptually) to a product involving three or more d functions, the author believes that it is quite difficult to prevent distributions constructed from 1D shot noise algorithms similar to that of Penman and MacNeil from inadvertently developing unintended (and unphysical) correlations in phase space. A different solution to the interference effect mentioned above was implemented by Reiche et al. into the fully 3D GENESIS code. They also employ multimember beamlets but then divide the macroparticle phase space into a group of "stripes" in dg, each one of which (in the absence of shot noise) has its average microbunching at the fundamental wavelength set to zero. However, no particular effort was paid to the initialization of higher harmonics.
With this justification for the need for a shot noise algorithm which will behave well in a multidimensional FEL simulation code, we now give details of the present GINGER algorithm.

II. THE PRESENT GINGER CODE SHOT NOISE INITIALIZATION ALGORITHM
The structure of the GINGER code follows that developed by Quimby [11] in the early 1980's in his own sideband instability code. Namely, there is a temporal simulation window W enclosing the electron beam (or a portion thereof) and the FEL radiation field. This window is then separated into n s discrete longitudinal sections whose centers are evenly spaced in time with a separation Dt s W͞n s . GINGER applies the eikonal approximation for the radiation field evolution in time and space and also wiggler period averaging for the interaction between the electrons and the field. Consequently, in the frequency domain, the window W corresponds to an interval Dv 6pn s ͞W centered on the central angular frequency 2pc͞l s . For short wavelength FEL's where the gain bandpass Dl G ഠ rl o is normally quite small compared with the wavelength l o of peak gain, one generally picks l s l o and a window length cW which is much greater than the so-called "cooperation" length l c ഠ l o ͞4pr. Here r is the normal 070701-2 070701-2 FEL efficiency parameter [12]. One also usually chooses cW greater than the slippage length l o L w ͞l w , where L w and l w are the undulator's length and period, respectively. Simultaneously, one also picks n s sufficiently large that cDt s # l c ͞8. This is necessary to ensure that the effective spectral window of the simulation is much greater than the FEL's gain bandpass. When r # 10 22 (as is true for most sub-mm wavelength SASE devices), cDt s ¿ l o . Each radiation and electron beam slice should then be thought of as representing a statistical average of all of the individual optical periods and ponderomotive wells contained within the interval Dt s . This point is important vis-à-vis shot noise statistics because the n p macroparticles in each GINGER electron beam slice (usually 1024 -32 768 in number) represent all N electrons in the slice's temporal duration Dt s , not just one optical period l o ͞c (as adopted in [1]).
As mentioned above, in order to model FEL's for which transverse effects (e.g., emittance, betatron motion, diffraction) and longitudinal effects (e.g., instantaneous energy spread) are important, the n p macroparticles in each GINGER electron beam slice are subdivided into an ensemble of n b beamlets, each composed of n p ͞n b ϵ 2M $ 2 macroparticles. Each such beamlet represents N b N͞n b ͑I B Dt s ͞n b e͒ actual electrons, where I B is the total beam current and e is the electron charge. Each of the macroparticles corresponding to a given beamlet is initialized with the same five-space coordinate ͑x, y, p x , p y , g͒ but the different individual beamlets have unique (and physically distinguishable experimentally) coordinates. For a given four-space ͑x, y, p x , p y ͒ distribution and instantaneous energy spread distribution s g about a mean energyḡ, GINGER employs a quiet-start load using multiple base, bit-reversed numbers (see [13,14]). Because of the usual dominance of the fundamental radial mode in the exponential growth regime, one often can use as few as 64 beamlets and still obtain reasonably accurate results. However, if either the effective fractional energy spread (instantaneous and/or that due to transverse emittance) is comparable to the FEL efficiency parameter r or if the optical mode size is only a fraction of the focused electron beam size, it is safer to use a much larger number (e.g., $256) of beamlets.
Each macroparticle "l" has a longitudinal coordinate u l which represents its ponderomotive phase measured relative to a plane electromagnetic wave. This phase evolves following the standard FEL resonance relation where the rms vector potential a 2 w incorporates the fast wiggle motion, p 2 Ќ represents the remaining perpendicular momentum corresponding to the "slow" betatron motion (both normalized to m e c), and we have assumed g ¿ 1.
Note that in the absence of coupling to the radiation field (which causes dg͞dz to be a function of u), all individual macroparticles corresponding to a given beamlet will maintain an identical five-space coordinate and therefore the same value of du͞dz. Consequently, whatever the initial absolute value for the microbunching at harmonic m $ 1 is for a given beamlet of index k (e.g., j P j2M j1 exp͑imu j,k ͒j͞2M͒, in the absence of radiation it will remain constant with z. However, beamlets with different values of g and/or total perpendicular momenta will have different values of du͞dz. This leads to both the phase and amplitude of the total bunching (i.e., summed over the n p macroparticles of all the beamlets) fluctuating with z even in the absence of radiation unless the initial bunching value for each and every beamlet is identically zero. When shot noise is included, this fluctuation occurs physically and its statistics provide a good test of the simulation code (see Sec. III).
Temporarily postponing the consideration of shot noise microbunching, the 2M macroparticles in each beamlet k are loaded uniformly over a 2p interval in u, e.g., For M $ 2 neither the z evolution of the radiation power nor that of the longitudinal microbunching of individual beamlets is sensitive to the chosen initial offsets u 0,k . To deal with M 1 (and to improve postprocessor plotting aesthetics), however, GINGER uniformly increases u 0,k with index k over the interval ͓0, p͞M͔. Also note that with this load, for each individual beamlet there is no initial bunching through harmonic number M, whatever the choice of u 0,k . To model the effects of shot noise, it is necessary to add a small random du j,k to each of the 2M macroparticles in beamlet k. From one beamlet to the next, there should be no statistical correlation in du. Otherwise, the variations in du͞dz from beamlet to beamlet could lead to anomalous secular variations with z in the total bunching. Consequently, the random numbers used to generate du j,k should be completely independent of those used to generate du j,k 0 for k fi k 0 . Moreover, we request that, for a given beamlet, the averages of the complex bunching corresponding to different harmonic numbers also be statistically uncorrelated. To do this numerically within the code, we set Each a m,k and b m,k are independently picked from a random one-dimensional Gaussian distribution centered about zero whose rms width s m obeys The factor 2͞m 2 comes from the relation 070701-3 070701-3 in the limit mb m,k ø 1 which is obeyed for N b ¿ 1. The sum over sin 2 is identically M for m , M. For m M, the sum varies from 0 to 2M depending upon the value u 0,k ; for a uniform distribution of u 0,k over the interval ͓0, 2p͔ the sum has an average value of M. The equivalent relation is found for ͗sinmu͘ k except b m,k is replaced by a m,k and the function sin 2 by cos 2 . Consequently, the average expectation value for the length squared of the Since the phase and amplitude of this phasor are uncorrelated for different beamlets, the average expectation value for the rms bunching at harmonic m # M due to all n p macroparticles is as expected for N ¿ 1.

III. NUMERICAL EXAMPLES
As a partial test of the shot noise algorithm discussed above, we have used GINGER to study two particular problems. The first examines the z evolution of microbunching statistics at different harmonics in the absence of coupling to the emitted radiation field. The second example examines the output characteristics of spontaneous emission for low current beam with negligible FEL gain.

A. Example 1: Random microbunching in a drifting low current beam with large energy spread
In the absence of any significant radiation field and for curved undulator pole face focusing [which leads to a wiggle-period-averaged constant value of a 2 w ͑r͒ 1 p 2 Ќ for each electron], individual electrons will maintain both a constant g and parallel momentum p z . As mentioned in the previous section, all the electrons comprising each individual beamlet k will have an identical longitudinal phase drift du k dz . However, in the presence of shot noise, a large energy spread (e.g., 2N w s g ͞ḡ ¿ 1) among the different beamlets will lead to a large variation in their phase drift and thus a continually evolving change with z in the complex microbunching phasor sum b m ϵ n 21 p P l expimu l over all macroparticles. Nonetheless, there should be no To test the behavior in this regard of the GINGER shot noise algorithm, we simulated a 2.5-kA beam with s g ͞ḡ 0.02 traversing a (hypothetical) 256-m long undulator (l w 0.08 m, N w 3200, a w 1.24) to examine the z-dependent behavior of the macroparticle microbunching. The radiation field was artificially fixed at the subnanowatt level. In addition to the instantaneous energy spread, the beam was given a large transverse emittance (´N 3.16 3 10 23 m rad) which corresponded to an additional effective fractional energy spread of 0.010. With equal focusing in both transverse planes, b x b y 1.45 m. The 16 384 macroparticles were loaded with 2M 16 (i.e., 1024 beamlets) to represent the N ഠ 5.52 3 10 8 real electrons in a single optical wavelength of l s 10.6 mm. Figure 1 shows the fundamental bunching jb 1 j versus z computed 1-m intervals along the undulator. Statistically for N ¿ 1, one expects that for the fundamental and each harmonic, b rms ϵ p †jbj 2 ‡ ഠ N 21͞2 , b ϵ †jbj ‡ ഠ ͑4N͞p͒ 21͞2 , and that the standard deviation s jbj ϵ p †͑b 2b͒ 2 ‡ ഠ q 1 2 p 4 N 21͞2 . Table I shows these computed statistics for the 257-z bunching values for harmonics numbers 1 through 4, together with the theoretically expected value. One sees that the computed statistics lie quite close to the expected values. Correlation tests between jb m j and z and between the individual harmonics were negative, as expected.

B. Example 2: Incoherent spontaneous emission
In a situation where the total FEL gain is small, the radiation from an initially unbunched beam (save for shot noise microbunching) traversing a periodic undulator will 070701-4 070701-4 be dominated by incoherent synchrotron emission. This emission provides the radiation seed for high gain SASE devices. Consequently, it is important that FEL simulation codes properly compute this emission (or at least that portion which couples to the coherent FEL gain). As a test of the shot noise algorithm (and concurrently the GINGER radiation emission and r 2 z propagation algorithms which for brevity we do not discuss here), we set up a test case of a 40-period, linearly polarized undulator with l w 25 mm and a w 1 (i.e., K p 2 resonant at l s 100 nm with a 255-MeV (g 500), 1 A electron beam with´N ഠ 1 3 10 26 m rad and no instantaneous energy spread. The FEL parameter r # 4 3 10 24 was sufficiently low so that any FEL gain should be very small over the 1-m interaction length. The GINGER simulation employed 256 beam and radiation slices, each enclosing one optical wavelength of particles, and each with 2048 macroparticles with 2M 8 and 256 beamlets. Periodic boundary conditions were applied in time; since the total beam length exceeded six slippage lengths, the periodicity should have negligible effect upon the results. The results discussed below are composed of averages over eight separate runs with different random number seeds. Figure 2 shows the time-averaged simulation results for both the far field, on-axis radiation intensity dP͞dV and the total radiated power versus z; the figure also plots the theoretically expected values using formulas [5.49] and [5.50a] of Atwood [15]. For both quantities, the simulation results show the expected linear dependence with z. The on axis dP͞dV (which is dominated by spectral components within a few percent of l s ) shows quite good quantitative agreement between the theoretical value and the simulation result. However, the simulation's P͑z͒ is a factor of ϳ5X low, whose cause we tentatively attribute to the limited wavelength bandpass of the simulation (Dl 650 nm) which consequently neglects radiation emitted off-axis to the red side of l ഠ 150 nm. The simulation's fractional bandpass, however, does far exceed 1͞N w , which approximates the fractional spectral width corresponding to the central angular cone of the undulator radiation at exit. This may explain why the agreement in dP͞dV can be so good despite the discrepancy in the total emitted power. Figure 3 shows the output P͑l͒ in the central wavelength region. The output power integrated over a 64% bandpass of 0.19 W is quite close to the predicted, theoretical output power contained with the central cone [15], where ͓JJ͔ is the standard Bessel function difference term. Consequently, the good agreement in the central cone power and the on axis dP͞dV between theory and GINGER simulation results suggests that the total power discrepancy is due either to the simulation's limited spectral bandpass and/or maximum angular bandpass.

IV. DISCUSSION
The results of these examples, together with those of numerous additional simulations corresponding to other parameter regimes, lead to some measure of confidence that the multibeamlet, multiharmonic, shot noise algorithm employed in GINGER gives a reasonably accurate model for shot noise effects in FEL's, including those where energy spread and/or transverse effects play an important role. One advantage of this model is that it can straightforwardly model the shot noise corresponding to "peculiar" initial electron beam configurations: e.g., beams composed of two or more energy components (as can be generated by "beam compressor" elements), beams with "holes" or very unsymmetrical phase space distributions (as could be generated by cathode abnormalities and/or off-center, nonlinear focusing elements). A great deal of effort is now being applied toward "beginning-to-end" modeling of both actual devices such as VISA and LEUTL and proposed FEL's such as the LCLS, and it is undoubtedly safe to assert that the predicted six-space distribution of the electron beam at undulator entrance will not be a simple convolution of Gaussian distributions in energy, configuration, and momentum space. This belief, together with the growing interest in the behavior of the higher harmonics, makes it even more important that whatever the algorithm underlying the shot noise model, it be able to deal with a reality far more complex than simple one-dimensional, monoenergetic, laminar flow.