A self-consistent numerical approach to track particles in FEL interaction with electromagnetic field modes

In this paper we present a novel approach to FEL simulations based on the decomposition of the electromagnetic field in a finite number of radiation modes. The evolution of each mode amplitude is simply determined by energy conservation. The code is developed as an expansion of the General Particle Tracer framework and adds important capabilities to the suite of well-established numerical simulations already available to the FEL community. The approach is not based on the period average approximation and can handle long-wavelength waveguide FELs as it is possible to include the dispersion effects of the boundaries. Futhermore, it correctly simulates lower charge systems where both transverse and longitudinal space charge forces still play a significant role in the dynamics. For free-space FEL interactions, a source dependent expansion approximation can be used to limit the number of transverse modes required to model the field profile and speed up the calculation of the system evolution. Three examples are studied in detail including a single pass FEL amplifier, the high efficiency TESSA266 scenario and a THz waveguide FEL operating in the zero-slippage regime.

In this paper we present a novel approach to FEL simulations based on the decomposition of the electromagnetic field in a finite number of radiation modes. The evolution of each mode amplitude is simply determined by energy conservation. The code is developed as an expansion of the General Particle Tracer framework and adds important capabilities to the suite of well-established numerical simulations already available to the FEL community. The approach is not based on the period average approximation and can handle long-wavelength waveguide FELs as it is possible to include the dispersion effects of the boundaries. Futhermore, it correctly simulates lower charge systems where both transverse and longitudinal space charge forces play a significant role in the dynamics. For free-space FEL interactions, a source dependent expansion approximation can be used to limit the number of transverse modes required to model the field profile and speed up the calculation of the system's evolution. Three examples are studied in detail including a single pass FEL amplifier, the high efficiency TESSA266 scenario, and a THz waveguide FEL operating in the zero-slippage regime.

I. INTRODUCTION
Numerical simulations have played a significant role in the development of X-ray Free Electron lasers [1]. As the theory underlying the FEL [2,3] only admits analytical solutions under strong approximations, accelerator physicists have over the years developed a well assorted suite of numerical approaches to better understand the details of the evolution of charged particles and electromagnetic fields in their interaction through magnetic undulators.
There are a large variety of FEL simulation codes and many good reviews on the subject have been given [4][5][6]. These range from fast one dimensional models (Perseo [7], Perave [8]) which help in quick design studies and can be used to explore time-dependent and non linear effects, to more complete 3D simulations (Ginger [9], Genesis 1.3 [10], Fast [11], Puffin [12], Minerva [13]) which include transverse effects and can simulate wakefields and complex beam distributions with correlations between the phase spaces. Each code has been (at least initially) developed to solve a particular FEL problem, but it has often been the case that, by comparing and understanding the various assumptions in each model, insights on the various physical processes taking place in an FEL system have been gained.
Here we discuss yet another instance of a three dimensional FEL simulation based on the decomposition of the electromagnetic field in a discrete set of transverse and frequency modes. In this respect it is more similar to the family of frequency-based codes like Puffin or Minerva. The code is built as an expansion of the widely available General Particle Tracer code for charged particle simulations [14]. In this sense, it can use a complete set of already built-in functions for beam transport and interface seamlessly with photoinjector [15] and CSR cal-culations [16]. This choice also brings several important advantages. The calculation does not resort to period averaging and a full (simulated or even measured) undulator field map can be used to move the particles. The effects of the interaction at the undulator entrance and exit can therefore be correctly evaluated. Furthermore, space charge effects are naturally incorporated, including the transverse space charge effects that at low beam energy play a significant role in the beam transport and evolution.
The code can be used to simulate both free-space and waveguide propagating electromagnetic fields and can take into account the dispersive properties of the medium. For free-space there is some freedom in choosing the basis for the field expansion, making it possible to take advantage of the Source Dependent Expansion [17,18] algorithm to reduce the number of modes needed to accurately describe the field and significantly speed up the calculation.
The paper is organized as follows. We first review the modal expansion and the equations implemented in the simulation [19]. We then make three different application examples. The first one is just a simple seeded FEL amplifier in vacuum (analyzed both in helical and planar geometry). The second one applies to the study of the system in the strong non linear regime and refers to the simulation of the TESSA266 experiment [20]. The final example is a waveguide THz FEL where the code is used to correctly simulate the zero-slippage amplification [21].

II. MODE EXPANSION
In order to self-consistently simulate the interaction between radiation and electrons, we begin with the Maxwell wave equation for the complex field amplitude wheren and x denote the polarization vector and transverse coordinates, respectively. Definingẑ as the direction of propagation, the polarization vector can be written in complex notation asn =x orn = (x ± iŷ)/ √ 2 for linearly and circularly polarized light. The polarization vector formalism is particularly convenient to unify the description of the planar and helical geometries. The time-averaged Poynting vector (representing the wave intensity) can be written in both cases as 0 c|E( x, z, t)| 2 /2.
If we write the scalar field amplitude in terms of its z-coordinate spatial Fourier transform the left hand side (LHS) of the equation can be rewritten as where we factored out the harmonic time-dependence and have neglected the second derivative of the slowly varying field amplitude, i.e. ∂ 2 E( x, k, t)/∂t 2 ω 2 E( x, k, t). The current density on the RHS can be written in complex notation using the particle positions and velocities where v j = √ 2K rms ce −ikuzj /γ jn represents the particle velocities in the undulator, K rms = eB rms /mck u is the root mean square (rms) undulator strength parameter, λ u = 2π/k u is the undulator period, e and m are the charge and mass of an electron, and γ j is the relativistic factor. Note that in most simulations a macroparticle model is used where one simulation particle represents multiple actual electrons in the beam. In this case, the sum in Eq. 4 will run over the macroparticle index.
Using nested Fourier transforms, we have The delta function allows easy integration over z . The time derivative is straightforward using chain rule with z j (t) after noticing that K and γ have a very slow dependence on z j ( dK dz k u and dγ dz k u ) and the transverse velocity is negligible.
Combining Eqns. 3 and 6, we can then rewrite Eq. (1) for the spatial frequency components of the field as where the source term is obtained by projecting the current density onton as Each spatial frequency component of the field can be further decomposed into an orthogonal mode basis labeled by index m and normalized such that is one of the complex mode solutions of the source-free wave equation (i.e. S = 0 in Eq. (7)) and A m is a normalization constant.
InsertingẼ( x, k, t) = m a m (t)Θ m ( x, k, t) into Eq. (7), we can multiply both sides of the equation by Θ * n and integrate over the transverse coordinates d x to obtain the mode amplitude excitation equation [19] where Θ m,j means evaluating the m-th mode at the jth particle position. As we sum over the particles, only the spatial frequencies that are nearly resonant with the particle speeds (β zj = β ph = ω/c(k + k u )) will contribute to a net energy exchange with the field so that the bracketed term can be approximated as 1. This mode excitation equation can also be independently derived from (and is fully consistent with) energy conservation. To see this, we write the energy of the system W using the spatial frequency Fourier transform of the electric field as After differentiating, we find The rate of change in the electromagnetic energy is the negative of the work done on the particles, where terms that do not satisfy the resonant condition average to zero in the particle sum. Equating the coefficients of a * m leads tȯ which matches our previous calculation in (9). In other words, the evolution of the amplitude of each electromagnetic mode in the system can be simply calculated by adding the energy changes induced by that mode on the particles.

A. GPT Numerical Implementation
In order to extend the capabilities of GPT to selfconsistently calculate the interaction with the radiation modes in the undulator, we based our development on the built-in function that computes the interaction with the modes of a gaussian optical resonator [22].
In the numerical model, the continuous integral of (2) is approximated using a discrete basis of spatial frequency modes (15) where the sum over index q includes both spatial frequencies and transverse modes. With respect to the previous section, u q and v q now represent the actual electric field amplitudes and have absorbed the user-defined mode separation interval ∆k and the 1/2π from the Fourier transform. Consequently, the source term in Eq. (7) also gains an additional factor of ∆k/2π.
In the input file, the user can specify the number of modes and the spatial frequency interval for the simulation. That choice of interval and associated spectral resolution should be taken judiciously to include the resonant frequency of the system and to correctly simulate the radiation bandwidth. Since the latter depends on various factors including the gain parameter, the length of the undulator, and the electron bunch length, it is always advisable to check the results for consistency and convergence as the number of modes and their separation is varied.
The choice of the spatial frequency interval defines the distance in the z-dimension L = 2π/∆k over which periodic boundary conditions are applied for the field. The frequencies ω q are determined from the longitudinal wavenumbers using the mode dispersion relation given by ω q = ck q in free space or ω 2 q = (k 2 mn + k 2 q )c 2 in a waveguide.
Writing the complex mode amplitude as Θ q = T q e iψq , we can then express the x and y component of the elec-tromagnetic field at time t at the particle locations as From these fields, the electromagnetic forces acting on the particles are computed at each time step. Particle velocities and positions are then used to self-consistently calculate the evolution of the amplitudes of each mode (u q and v q ) according to Eq. 9.
It is also possible to run the code in single frequency mode. In this case, the field is assumed to be perfectly periodic, with only one spatial frequency term in Eq. 15 and the time-averaged sum (now only running over the transverse modes) 0π ∆k q (u 2 q + v 2 q )A q corresponds to the total radiation field power.

B. Curved Parallel Plate Waveguide
The geometry of the interaction to be simulated determines the choice of the mode basis Θ m and the associated dispersion relation. An important application of our new code is the study of the evolution of an FEL system in a waveguide. The dispersive properties in the waveguide can not be easily modeled in conventional FEL codes which adopt a time-dependent (slice) model for the description of the radiation. For this case we can expand the field in the complete set of orthonormal modes for the particular waveguide cross-section under study. Here we focus on the TE modes of a curved parallel plate waveguide [21] where the fields can be written in terms of the longitudinal component of the magnetic field H z . The transverse wavenumber k mn for the modes in the waveguide is written as where b is the separation and R is the radius of curvature of the waveguide. The confocal case (i.e R = b) minimizes diffraction losses and is typically employed in practice [23], but in the numerical model these parameters can be chosen by the user separately. The dispersion relation is then expressed as Analytical expressions for the field E m,n (r ⊥ ) in the guide can be found in [23] and [24]. The longitudinal field Φ mn , corresponding to H z for TE modes and E z for TM modes can be written in terms of Hermite polynomials He m as where α mn (y) = 1 + 4 The transverse field components are then calculated as The effective mode area is hard-coded in the software.

C. Free space propagation Source Dependent Expansion
Another important case is where the waveguide boundaries are removed or very far away so that one can use free-space modes to describe the radiation field. Either Laguerre-Gaussian or Hermite-Gaussian modes can be used depending on the symmetry of the problem. Assuming azimuthal symmetry (i.e. r 2 = | x| 2 ), we start by writing the complex scalar field amplitude as a sum of different spatial frequency Laguerre-Gaussian modes, E( x, z, t) = n,m a n,m (t)Θ n.m (r, t)e iknz−iωnt (23) where we explicitly show that the sum index runs over the different spatial frequencies (n) and the transverse mode numbers (m). The modal basis for the field expansion can be written as Θ n,m (r, t) = 1 where L m is the Laguerre polynomial of order m, w n and α n indicate the waist size and the curvature of the phase fronts for the mode having spatial frequency k n , and ψ n (t) = arctan α n (t). In the case that no electron beam is present and the radiation is freely diffracting, w n (t) = w 0,n 1 + c 2 t 2 /z 2 r,n and α n (t) = ct/z r,n with the implicit frequency dependence in z r,n = k n w 2 0,n /2, the Rayleigh range of the nth-mode. The mode area normalization constants are The effectiveness of the Laguerre-Gaussian mode expansion depends critically on the choice of the waist size and location, and in the absence of any prior knowledge or extra information, the simulation should include a large number of transverse modes in order to accurately model the radiation field.
In many cases, as for example when the FEL is seeded with an external laser and the radiation transverse profile is mainly dominated by one or a few modes, it is a good approximation to truncate the sum to only include a small number of terms. To further minimize this number (and proportionally speed up the computational time), it is possible to take advantage of the source dependent expansion originally developed for the FEL framework by Sprangle et al. [17] where the waist size and location of the expansion are adjusted along the interaction.
Following the original work in [17] (recently revisited by Baxevanis et al. [18]), after plugging Eq. 23 into the inhomogeneous wave equation, we obtain a coupled system of differential equations for the mode amplitudes in terms of the projections of the source term onto the mode basis.
Using the definition of S from (8), it is possible to write the source projection moments F m,n in terms of sums over the particle (or macroparticle) coordinates. We can then solve for how w n and α n should vary in order to truncate the system at the desired order. For example, neglecting all m ≥ 1 we get ∂u n ∂t = G n (α n u n − v n ) + (u n B I,n + v n B R,n ) + F 0I,n ∂v n ∂t = G n (u n + α n v n ) + (v n B I,n − u n B R,n ) − F 0R,n ∂α n ∂t = 2(1 + α 2 n )c 2 ωw 2 n + 2B R,n − 2α n B I,n ∂w n ∂t = 2c 2 α n ω n w n − w n B I,n where G n = 2 1+α 2 n (B R,n − α n B I,n ). B n represents the correction to the mode waist and radius induced by the source and can be written as A closer inspection to Eq. 26 c and d indicates that c/|B n | is a distance which sets the scale for the variation of the mode radius. In multi-frequency simulations, the modes with small initial amplitudes cause the magnitude of B n to diverge. This is taken care of by setting a userdefined input parameter L thresh which limits the spot size variation along the interaction by setting B n = 0 whenever c/|B n | < L thresh .
The equations for radiation evolution are then selfconsistently solved with the GPT equations of motion for the macroparticles.
The general equations for complex mode evolution and B n with M spatial modes arė a n,m = [B I,n + α n G n + i(2m + 1)(G n − B R,n )] a n,m + imB n e 2iψn a n,m−1 + i(m + 1)B * n e −2iψn a n,m+1 − iF m,n a n,m≥M = 0 =⇒ B n = F M,n e −2iψn M a n,M −1 .
Higher order modes with small initial amplitudes are initially considered perturbations to the gaussian mode such that (27) still holds. Once the approximation |a 1 |/|a 0 | 1 breaks down (≈ .01), the correct definition of B n from (28) can be used without divergence or significant numerical noise. In practice, errors from the perturbative approximation are negligible since it is accurate far into the linear regime.

D. Quiet start
In multifrequency simulations where many longitudinal wavenumbers and corresponding frequencies are used to simulate the field along a finite length bunch, it is critical to pay attention to the details associated with loading the particle coordinates in the simulation. Because it is common to have a much smaller number of macroparticles than real number of electrons, the noise in the bunching source term can be unacceptably high, causing unphysical growth of the field along the undulator.
This problem is common and well discussed in the vast literature of simulations for FELs [25,26]. While there are a number of possible solutions, our situation is slightly complicated as we need to ensure that the intrinsic bunching is and remains very small for all of the discrete frequencies in the simulation. This first requires equally distributing particles in the z-coordinate over a length L = 2π/∆k. For example in Fig. 2 we show the input phase space when the simulation spans a bandwidth of 3 % around the central wavelength of 266nm. In this case, the beam longitudinal profile (a gaussian with rms bunch length 30 µm) is initialized by assigning a different charge weight to each macroparticle. When shot-noise effects are desired, each macroparticle's position is shifted by a small dz according to well described algorithms [27,28] to achieve the correct statistics.
In addition, it is important to make sure that the noise from other coordinates would not contribute to a growth of the bunching as the beam propagates in the absence of an interaction. This is taken care of by mirroring the energy, transverse coordinates, and momenta over a large number of 5D phase space bins. The number of bins (typically larger than 32) should be chosen such that bunching in the absense of an interation remains small for all the discrete frequencies included in the simulations.

III. EXAMPLES
We limit this discussion to three examples that highlight the main features of our approach, even though it is expected that the new code can be successfully applied to a variety of other situations. The first case considered is a classical single-pass FEL seeded amplifier which will enable a quantitative comparison with the semi-analytical M. Xie formulas [29] as well as with a traditional period-average code like Genesis for both planar and helical geometries. The second example is relevant to the TESSA266 experiment being planned at the LEA beamline at the APS linac in Argonne National Laboratory aiming at very high conversion efficiency at 266 nm [20]. This case serves to illustrate the capability of using a 3D magnetic field map for a fairly complicated segmented tapered undulator. The code compares well with a traditional FEL code like Genesis, even deep in the non-linear regime. The details of the beam transport (injection, entrance and exit sections and especially undulator break sections) can only be included in Genesis by using a linear beam transport approximation. GPT follows the evolution of the beam distribution along the beamline using field maps for all the magnetic elements (undulators, quadrupoles and phase shifter dipoles) and calculates energy exchange using the self-consistent interaction with the free-space modes. The results allow us to quantitatively include the effects of the entrance and exit sections (which add an effective 0.5 periods of interaction on each side of the undulator) and the trajectories after the prebuncher and in between the undulators. The final example is a waveguide THz FEL where GPT-FEL is used to correctly simulate the zero-slippage amplification. In this configuration, the strong dispersive properties of the guide affect the interaction which takes place in the zero-slippage regime. This scenario highlights a unique capability of our code which would be particularly challenging to simulate with traditional FEL codes.

A. FEL amplifier
The parameters for this example are reported in Table I and somewhat arbitrarily chosen to be similar to an un-tapered version of the TESSA266 experiment discussed below. The main differences are that a 200 period long undulator (with no break-section) is used for this example and the input seed power is lowered to 10 kW. An analytical model for the undulator magnetic field is used. The beam is transversely matched to the undulator natural focusing (equally distributed in the horizontal and vertical plane) so that its rms spot size remains nearly constant along the interaction. The main goal of this example is to benchmark GPTFEL against the fitting formulas for the 3D gain length of an untapered FEL amplifier and compare with a conventional FEL code like Genesis. We also used this example to evaluate the performance of the single mode SDE approximation versus a simulation with n = 11 azimuthally symmetric Laguerre Gaussian SDE modes to decompose the electromagnetic field. GPTFEL took 1.5 minute to simulate 76800 particles on an 8 processor for the single SDE mode and 5 minutes for 11 SDE modes. The time-independent, single frequency results for the planar and helical geometries are shown in Figure 3 and compared with Genesis 1.3. When using multiple spatial modes, the gain lengths in the planar and helical case are in good agreement (within 10 %) of the semianalytical and numerical model predictions. The radiation spot sizes defined by σ 2 r = 1 2 r 2 |E| 2 d 2 x |E| 2 d 2 x also closely follow the prediction. Note that while a single SDE mode is able to achieve qualitative results up to and near saturation, a larger number of spatial modes is required to correctly simulate the evolution of the radiation profile after saturation.
The multi-frequency simulation used an SDE gaussian mode for 31 spatial frequencies with a 6% bandwidth to simulate 128,000 particles in 23 minutes. The userdefined parameter L thresh limits the spot size variation along the undulator. Figure 4a shows a waterfall plot in the electron beam frame normalized at each z position to display the relative velocity of the radiation wavepacket, which is close to the beam velocity in the exponential regime and becomes superluminal in the non linear regime [30]. In Figure 4b, the spectrum just before saturation is shown as a function of L thresh normalized to the gain length. If an increased spectral resolution is required, computation time scales linearly with number of tracked modes.

B. TESSA266
In this next example we take advantage of the GPT functions to track the electron beam in the fairly complex transport line of the TESSA 266 experiment. The beamline includes a short, 8 period undulator followed by a 3 dipole chicane to convert the imprinted energy modulation into microbunching. Quadrupole doublets match the beam transversely into the focusing channel of each 0.96 meter, strongly tapered undulator section. A small dipole is placed between the second quadrupole doublet so that the three magnets can be used as a phase shifter between the undulator sections.
The GPT transport functions are used to set up the trajectory and the beam optics prior of turning on the seed and the FEL interaction module. Our timeindependent simulation of the TESSA266 beamline includes 21 higher order spatial modes to ensure an accurate modeling of the radiation profile. A 1 GW peak power input radiation pulse is focused at the entrance of the tapered undulator to a waist of 0.3 mm. The simulation is compared with Genesis results, but it should be noted that GPTFEL uses full 3D magnetic field maps for the undulators as well as for the dipoles and quadrupoles in the system. The magnetic field in the chicane dipoles is fine tuned to maximize the bunching and simultaneously optimize the injection phase of the beamlets relative to the radiation phase at the entrance of the tapered undu-lator. In Genesis, both the R56 and phase shifts are applied post-facto to the beam distribution at the entrance of the tapered undulator, explaining the large difference in the bunching factor evolution in Fig. 6a. In practice, the phase shifter between the tapered undulator sections had to be re-optimized to account for the additional slippage incurred by the beam when passing in the entrance and exit section of the wigglers. This is accomplished by horizontally shifting the quadrupoles in opposite directions to steer the beam and tuning the magnetic field amplitude of the dipole to recover a straight trajectory while maximizing the energy exchange in the second undulator.

C. Zero slippage THz FEL
A final example to showcase the capabilities of the new GPTFEL code is the simulation of a THz FEL operating in the zero-slippage regime [31]. The size of the waveg- GPTFEL correctly simulates the waveguide dispersive properties as shown in Fig. 8 by plotting the electric field at the entrance and exit of the 1 meter long waveguide system in the absence of strong interaction (i.e. for very low charge beams).
The parameters of this example are summarized in Table II. We have chosen a planar undulator geometry with equally distributed focusing in the horizontal and vertical plane. In this case, the largest coupling is obtained with the TE10 mode profile of the curved parallel plate waveguide. The beam is initialized at the entrance of the simulation with a large bunching factor (0.5) while we set the amplitude of the initial input seed to zero.
There are two main advantages of using the waveguide in this system. First, the waveguide maintains a constant radiation cross section along the interaction, avoiding diffraction effects. Second the waveguide's dispersive properties enable a zero slippage interaction. This large bandwidth interaction can drive the FEL with a much shorter beam because the slippage effects are effectively minimized and the radiation continues to interact and exchange energy with the particles even after a large number of periods. The simulation results are shown in Fig.  9 where the THz electric field waveform and the electron beam longitudinal phase space are shown to be temporally overlapping at the end of the undulator. Note that the system evolves in the non linear regime from the be- ginning as the electron beam enters the undulator with a very large bunching at the 1 THz resonant frequency induced by modulating the photocathode drive laser [32]. The undulator is linearly tapered starting from its half way point with a relative change in normalized vector potential K of 30 %/m to avoid saturation effects due to particles falling off the resonance curve. The efficiency of  conversion is above 10 % in this example.

IV. CONCLUSIONS AND OUTLOOK
A new approach for FEL simulations has been presented. The characteristic features are the decomposition of the field in a set of spatial and frequency modes and the integration with the GPT numerical integration engine which allows access and compatibility with a large number of beam transport designs and functions. There are a number of research opportunities which go beyond the scope of this paper but will be the subject of future studies, including a detailed study of the effects of the transverse space charge forces and an upgrade to include higher harmonic interactions. Parallelization of the code will allow much faster run times, increasing the number of macroparticles and modes that can be simulated. GPT-FEL is not expected to replace traditional approaches to FEL numerical simulations, but is intended to be a research tool to explore the interaction of relativistic electrons and electromagnetic waves in undulator systems in regimes where the approximations of standard FEL codes are questionable. The application of GPTFEL to dispersive systems allows for exploration of novel interaction regimes like the tapered waveguide THz FEL.