A self-consistent numerical approach to track particles in FEL interaction with electromagnetic field modes
AA self-consistent numerical approach to track particlesin FEL interaction with electromagnetic field modes
A. Fisher and P. Musumeci
University of California at Los Angeles, Los Angeles, CA, 90066
S.B. Van der Geer
Pulsar Physics, Eindhoven, The Netherlands (Dated: August 6, 2020)In this paper we present a novel approach to FEL simulations based on the decomposition of theelectromagnetic field in a finite number of radiation modes. The evolution of each mode amplitudeis simply determined by energy conservation. The code is developed as an expansion of the GeneralParticle Tracer framework and adds important capabilities to the suite of well-established numericalsimulations already available to the FEL community. The approach is not based on the periodaverage approximation and can handle long-wavelength waveguide FELs as it is possible to includethe dispersion effects of the boundaries. Futhermore, it correctly simulates lower charge systemswhere 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 limitthe number of transverse modes required to model the field profile and speed up the calculation ofthe 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-slippageregime.
I. INTRODUCTION
Numerical simulations have played a significant role inthe development of X-ray Free Electron lasers [1]. As thetheory underlying the FEL [2, 3] only admits analyticalsolutions under strong approximations, accelerator physi-cists have over the years developed a well assorted suiteof numerical approaches to better understand the detailsof the evolution of charged particles and electromagneticfields in their interaction through magnetic undulators.There are a large variety of FEL simulation codes andmany good reviews on the subject have been given [4–6].These range from fast one dimensional models (Perseo[7], Perave [8]) which help in quick design studies and canbe 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 includetransverse effects and can simulate wakefields and com-plex beam distributions with correlations between thephase spaces. Each code has been (at least initially) de-veloped to solve a particular FEL problem, but it hasoften been the case that, by comparing and understand-ing the various assumptions in each model, insights onthe various physical processes taking place in an FELsystem have been gained.Here we discuss yet another instance of a three dimen-sional FEL simulation based on the decomposition of theelectromagnetic field in a discrete set of transverse andfrequency modes. In this respect it is more similar tothe family of frequency-based codes like Puffin or Min-erva. The code is built as an expansion of the widelyavailable General Particle Tracer code for charged parti-cle simulations [14]. In this sense, it can use a completeset of already built-in functions for beam transport andinterface seamlessly with photoinjector [15] and CSR cal- culations [16]. This choice also brings several importantadvantages. The calculation does not resort to periodaveraging and a full (simulated or even measured) undu-lator field map can be used to move the particles. Theeffects of the interaction at the undulator entrance andexit can therefore be correctly evaluated. Furthermore,space charge effects are naturally incorporated, includ-ing the transverse space charge effects that at low beamenergy play a significant role in the beam transport andevolution.The code can be used to simulate both free-spaceand waveguide propagating electromagnetic fields andcan take into account the dispersive properties of themedium. For free-space there is some freedom in choos-ing the basis for the field expansion, making it possibleto take advantage of the Source Dependent Expansion[17, 18] algorithm to reduce the number of modes neededto accurately describe the field and significantly speed upthe calculation.The paper is organized as follows. We first review themodal expansion and the equations implemented in thesimulation [19]. We then make three different applica-tion examples. The first one is just a simple seeded FELamplifier in vacuum (analyzed both in helical and planargeometry). The second one applies to the study of thesystem in the strong non linear regime and refers to thesimulation of the TESSA266 experiment [20]. The finalexample is a waveguide THz FEL where the code is usedto correctly simulate the zero-slippage amplification [21].
II. MODE EXPANSION
In order to self-consistently simulate the interactionbetween radiation and electrons, we begin with the a r X i v : . [ phy s i c s . acc - ph ] A ug Maxwell wave equation for the complex field amplitude (cid:18) ∇ ⊥ + ∂ ∂z − c ∂ ∂t (cid:19) E ( (cid:126)x, z, t ) = µ ∂(cid:126) J ( (cid:126)x, z, t ) · ˆ n ∗ ∂t (1)where ˆ n and (cid:126) x denote the polarization vector and trans-verse coordinates, respectively. Defining ˆ z as the direc-tion of propagation, the polarization vector can be writ-ten in complex notation as ˆ n = ˆ x or ˆ n = (ˆ x ± i ˆ y ) / √ (cid:15) c | E ( (cid:126)x, z, t ) | / E ( (cid:126)x, z, t ) = 12 π (cid:90) ∞−∞ ˆ E ( (cid:126)x, k, t ) e ikz − iωt dk, (2)the left hand side (LHS) of the equation can be rewrittenas LHS = 12 π (cid:90) ∞−∞ (cid:18) ∇ ⊥ − k + ω /c + 2 iωc ∂∂t (cid:19) ˆ E ( (cid:126)x, k, t ) e ikz − iωt dk (3)where we factored out the harmonic time-dependence andhave neglected the second derivative of the slowly varyingfield amplitude, i.e. ∂ E ( (cid:126)x, k, t ) /∂t (cid:28) ω E ( (cid:126)x, k, t ).The current density on the RHS can be written in com-plex notation using the particle positions and velocities (cid:126) J ( z, t ) = (cid:88) j q j (cid:126) v j δ ( (cid:126)x − (cid:126)x j ( t )) δ ( z − z j ( t )) , (4)where (cid:126) v j = √ K rms ce − ik u z j /γ j ˆ n represents the particlevelocities in the undulator, K rms = eB rms /mck u is theroot mean square (rms) undulator strength parameter, λ u = 2 π/k u is the undulator period, e and m are thecharge and mass of an electron, and γ j is the relativisticfactor. Note that in most simulations a macroparticlemodel is used where one simulation particle representsmultiple actual electrons in the beam. In this case, thesum in Eq. 4 will run over the macroparticle index.Using nested Fourier transforms, we have RHS = µ π ∂∂t (cid:20)(cid:90) (cid:90) ∞−∞ (cid:126) J ( (cid:126)x, z (cid:48) , t ) e − ikz (cid:48) dz (cid:48) e ikz dk (cid:21) · ˆ n ∗ . (5)The delta function allows easy integration over z (cid:48) . Thetime derivative is straightforward using chain rule with z j ( t ) after noticing that K and γ have a very slow depen-dence on z j ( dKdz (cid:28) k u and dγdz (cid:28) k u ) and the transversevelocity is negligible. RHS = − iµ π (cid:90) ∞−∞ (cid:88) j q j cβ z,j ( k u + k )( (cid:126) v j · ˆ n ∗ ) × δ ( (cid:126)x − (cid:126)x j ) e − ikz j + ikz dk (6) Combining Eqns. 3 and 6, we can then rewrite Eq. (1)for the spatial frequency components of the field as (cid:18) ∇ ⊥ − k + ω c + 2 iωc ∂∂t (cid:19) ˆ E ( (cid:126)x, k, t ) = S ( (cid:126)x, k, t ) (7)where the source term is obtained by projecting the cur-rent density onto ˆn as S ( (cid:126)x, k, t ) = (cid:88) j − iµ cq j β z,j ( k u + k )( (cid:126) v j · ˆ n ∗ ) δ ( (cid:126)x − (cid:126)x j ) e − ikz j + iωt . (8)Each spatial frequency component of the field canbe further decomposed into an orthogonal mode ba-sis labeled by index m and normalized such that (cid:82)(cid:82) Θ ∗ m Θ n d(cid:126)x = δ mn A m where Θ m ( (cid:126)x, k, t ) is one of thecomplex mode solutions of the source-free wave equation(i.e. S = 0 in Eq. (7)) and A m is a normalization con-stant.Inserting ˜ E ( (cid:126)x, k, t ) = (cid:80) m a m ( t )Θ m ( (cid:126)x, k, t ) into Eq.(7), we can multiply both sides of the equation by Θ ∗ n and integrate over the transverse coordinates (cid:82)(cid:82) d(cid:126)x toobtain the mode amplitude excitation equation [19]˙ a m = − (cid:88) j q j (cid:15) A m (cid:20) cβ z,j ( k u + k ) ω (cid:21) ( (cid:126) v j · ˆ n ∗ )Θ ∗ m,j e − ikz j + iωt (9)where Θ m,j means evaluating the m-th mode at the jthparticle position. As we sum over the particles, only thespatial frequencies that are nearly resonant with the par-ticle speeds ( β z j = β ph = ω/c ( k + k u )) will contribute toa net energy exchange with the field so that the bracketedterm can be approximated as 1.This mode excitation equation can also be indepen-dently derived from (and is fully consistent with) energyconservation. To see this, we write the energy of the sys-tem W using the spatial frequency Fourier transform ofthe electric field as W = 12 (cid:15) π (cid:90) (cid:88) m | a m | A m dk. (10)After differentiating, we find dWdt = (cid:90) (cid:88) m a ∗ m π (cid:20) ˙ a m (cid:15) A m (cid:21) dk + c.c. (11)The rate of change in the electromagnetic energy is thenegative of the work done on the particles, (cid:88) j (cid:126) F j · (cid:126) v j = − (cid:88) j q j (cid:60) ( E ( (cid:126)x, z, t )ˆ n ) · (cid:60) ( (cid:126) v j )= (cid:88) j q j E ∗ ( (cid:126)x, z, t )ˆ n ∗ · (cid:126) v j ) + c.c. (12)= (cid:90) (cid:88) m a ∗ m π (cid:88) j q j ∗ m,j e − ikz j + iωt (ˆ n ∗ · (cid:126) v j ) dk + c.c. (13)where terms that do not satisfy the resonant conditionaverage to zero in the particle sum. Equating the coeffi-cients of a ∗ m leads to˙ a m = − (cid:88) j q j (cid:15) A m ( (cid:126) v j · ˆ n ∗ )Θ ∗ m,j e − ikz j + iωt (14)which matches our previous calculation in (9). In otherwords, the evolution of the amplitude of each electro-magnetic mode in the system can be simply calculatedby adding the energy changes induced by that mode onthe particles. A. GPT Numerical Implementation
In order to extend the capabilities of GPT to self-consistently calculate the interaction with the radiationmodes in the undulator, we based our development onthe built-in function that computes the interaction withthe modes of a gaussian optical resonator [22].In the numerical model, the continuous integral of (2)is approximated using a discrete basis of spatial frequencymodes (cid:126) E ( (cid:126)x, z, t ) = (cid:88) q ( u q + iv q )Θ q ( (cid:126)x, k, t ) e ik q z − iω q t ˆ n (15)where the sum over index q includes both spatial frequen-cies and transverse modes. With respect to the previoussection, u q and v q now represent the actual electric fieldamplitudes and have absorbed the user-defined mode sep-aration interval ∆ k and the 1 / π from the Fourier trans-form. Consequently, the source term in Eq. (7) also gainsan additional factor of ∆ k/ π .In the input file, the user can specify the number ofmodes and the spatial frequency interval for the simu-lation. That choice of interval and associated spectralresolution should be taken judiciously to include the res-onant frequency of the system and to correctly simulatethe radiation bandwidth. Since the latter depends onvarious factors including the gain parameter, the lengthof the undulator, and the electron bunch length, it is al-ways advisable to check the results for consistency andconvergence as the number of modes and their separationis varied.The choice of the spatial frequency interval definesthe distance in the z-dimension L = 2 π/ ∆ k over whichperiodic boundary conditions are applied for the field.The frequencies ω q are determined from the longitudinalwavenumbers using the mode dispersion relation givenby ω q = ck q in free space or ω q = ( k mn + k q ) c in awaveguide.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 E x ( (cid:126)x j , z j , t ) = (cid:88) q T q ( u q cos φ q − v q sin φ q ) | ˆ n · (cid:126) x | E y ( (cid:126)x j , z j , t ) = − (cid:88) q T q ( u q sin φ q + v q cos φ q ) | ˆ n · (cid:126) y | (cid:126) B = 1 ω q ˆ k q × (cid:126) E (16)where φ q = k q z j − ω q t + ψ q .From these fields, the electromagnetic forces acting onthe particles are computed at each time step. Particlevelocities and positions are then used to self-consistentlycalculate 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 frequencymode. In this case, the field is assumed to be perfectlyperiodic, with only one spatial frequency term in Eq. 15and the time-averaged sum (now only running over thetransverse modes) (cid:15) π ∆ k (cid:80) q ( u q + v q ) A q corresponds to thetotal radiation field power. B. Curved Parallel Plate Waveguide
The geometry of the interaction to be simulated deter-mines the choice of the mode basis Θ m and the associateddispersion relation. An important application of our newcode is the study of the evolution of an FEL system ina waveguide. The dispersive properties in the waveguidecan not be easily modeled in conventional FEL codeswhich adopt a time-dependent (slice) model for the de-scription of the radiation. For this case we can expandthe field in the complete set of orthonormal modes for theparticular waveguide cross-section under study. Here wefocus on the TE modes of a curved parallel plate waveg-uide [21] where the fields can be written in terms of thelongitudinal component of the magnetic field H z . TE11 y ) x E y xy , , () () E x xy , , () () + dd TE10 E FIG. 1: TE10 and TE11 y-component of the electricfield for a curved parallel plate waveguide. The TE10mode is the one that has the largest FEL coupling toextract energy from a relativistic electron beam.The transverse wavenumber k mn for the modes in thewaveguide is written as k mn = 1 b (cid:18) nπ + (2 m + 1) tan − b √ Rb − b (cid:19) (17)where b is the separation and R is the radius of curvatureof the waveguide. The confocal case (i.e R = b ) minimizesdiffraction losses and is typically employed in practice[23], but in the numerical model these parameters can bechosen by the user separately.The dispersion relation is then expressed as k z ( ω ) = (cid:114) ω c − k mn . (18)Analytical expressions for the field E m,n ( r ⊥ ) in theguide can be found in [23] and [24]. The longitudinalfield Φ mn , corresponding to H z for TE modes and E z forTM modes can be written in terms of Hermite polyno-mials He m asΦ mn = e − β mn x α mn ( y ) (cid:112) α mn ( y ) He m (cid:32) β mn x (cid:112) α mn ( y ) (cid:33) e ± ik z z (cid:20) cossin (cid:21) (cid:20) k mn y + 2 β mn yx k mn α mn ( y ) − ( m + 12 ) arctan 2 β mn yk mn (cid:21) (19)where α mn ( y ) = 1 + 4 β mn y k mn β mn = (cid:115) k mn √ Rb − b . (20)The transverse field components are then calculated as E ( x,y ) = − ik mn (cid:18) k z ∂E z ∂ ( x, y ) ± ωµ ∂H z ∂ ( y, x ) (cid:19) . (21)The effective mode area A mn = (cid:82) |E mn ( r ⊥ ) | dr ⊥ | E peak | (22)is hard-coded in the software. C. Free space propagation Source DependentExpansion
Another important case is where the waveguide bound-aries are removed or very far away so that one can usefree-space modes to describe the radiation field. EitherLaguerre-Gaussian or Hermite-Gaussian modes can beused depending on the symmetry of the problem. As-suming azimuthal symmetry (i.e. r = | (cid:126)x | ), we start by writing the complex scalar field amplitude as a sum ofdifferent spatial frequency Laguerre-Gaussian modes, E ( (cid:126)x, z, t ) = (cid:88) n,m a n,m ( t )Θ n.m ( r, t ) e ik n z − iω n t (23)where we explicitly show that the sum index runs over thedifferent spatial frequencies (n) and the transverse modenumbers (m). The modal basis for the field expansioncan be written asΘ n,m ( r, t ) = 1 (cid:112) α n ( t ) L m (cid:18) r w n ( t ) (cid:19) e − r /w n ( t ) × e iα n ( t ) r /w n ( t ) − i (2 m +1) ψ n ( t ) (24)where L m is the Laguerre polynomial of order m , w n and α n indicate the waist size and the curvature of thephase fronts for the mode having spatial frequency k n ,and ψ n ( t ) = arctan α n ( t ). In the case that no electronbeam is present and the radiation is freely diffracting, w n ( t ) = w ,n (cid:113) c t /z r,n and α n ( t ) = ct/z r,n withthe implicit frequency dependence in z r,n = k n w ,n / A n,m = πw ,n / . (25)The effectiveness of the Laguerre-Gaussian mode ex-pansion depends critically on the choice of the waist sizeand location, and in the absence of any prior knowledgeor extra information, the simulation should include alarge number of transverse modes in order to accuratelymodel the radiation field.In many cases, as for example when the FEL is seededwith an external laser and the radiation transverse pro-file is mainly dominated by one or a few modes, it is agood approximation to truncate the sum to only include asmall number of terms. To further minimize this number(and proportionally speed up the computational time), itis possible to take advantage of the source dependent ex-pansion originally developed for the FEL framework bySprangle et al. [17] where the waist size and location ofthe expansion are adjusted along the interaction.Following the original work in [17] (recently revisitedby Baxevanis et al. [18]), after plugging Eq. 23 intothe inhomogeneous wave equation, we obtain a coupledsystem of differential equations for the mode amplitudesin terms of the projections of the source term onto themode basis. F m,n = c ω n πw ,n (cid:90) S ( r )Θ ∗ m ( r ) d(cid:126)x Using the definition of S from (8), it is possible to writethe source projection moments F m,n in terms of sumsover the particle (or macroparticle) coordinates.We can then solve for how w n and α n should vary inorder to truncate the system at the desired order. Forexample, neglecting all m ≥ ∂u n ∂t = G n ( α n u n − v n ) + ( u n B I,n + v n B R,n )+ F I,n ∂v n ∂t = G n ( u n + α n v n ) + ( v n B I,n − u n B R,n ) − F R,n ∂α n ∂t = 2(1 + α n ) c ωw n + 2 B R,n − α n B I,n ∂w n ∂t = 2 c α n ω n w n − w n B I,n (26)where G n = α n ( B R,n − α n B I,n ). B n represents thecorrection to the mode waist and radius induced by thesource and can be written as B n = F n e − iψ n /a n . (27)A closer inspection to Eq. 26 c and d indicates that c/ | B n | is a distance which sets the scale for the variationof the mode radius. In multi-frequency simulations, themodes with small initial amplitudes cause the magnitudeof B n to diverge. This is taken care of by setting a user-defined input parameter L thresh which limits the spotsize variation along the interaction by setting B n = 0whenever c/ | B n | < L thresh .The equations for radiation evolution are then self-consistently solved with the GPT equations of motionfor the macroparticles.The general equations for complex mode evolution and B n with M spatial modes are˙ a n,m = [ B I,n + α n G n + i (2 m + 1)( G n − B R,n )] a n,m + imB n e iψ n a n,m − + i ( m + 1) B ∗ n e − iψ n a n,m +1 − iF m,n a n,m ≥ M = 0 = ⇒ B n = F M,n e − iψ n M a n,M − . (28)Higher order modes with small initial amplitudes are ini-tially considered perturbations to the gaussian mode suchthat (27) still holds. Once the approximation | a | / | a | (cid:28) ≈ . B n from(28) can be used without divergence or significant nu-merical noise. In practice, errors from the perturbativeapproximation are negligible since it is accurate far intothe linear regime. D. Quiet start
In multifrequency simulations where many longitudi-nal wavenumbers and corresponding frequencies are usedto simulate the field along a finite length bunch, it iscritical to pay attention to the details associated withloading the particle coordinates in the simulation. Be-cause it is common to have a much smaller number of macroparticles than real number of electrons, the noisein the bunching source term can be unacceptably high,causing unphysical growth of the field along the undula-tor.This problem is common and well discussed in thevast literature of simulations for FELs [25, 26]. Whilethere are a number of possible solutions, our situationis slightly complicated as we need to ensure that the in-trinsic bunching is and remains very small for all of thediscrete frequencies in the simulation. This first requiresequally distributing particles in the z-coordinate over alength L = 2 π/ ∆ k . For example in Fig. 2 we show theinput phase space when the simulation spans a band-width of 3 % around the central wavelength of 266nm.In this case, the beam longitudinal profile (a gaussianwith rms bunch length 30 µ m) is initialized by assigninga different charge weight to each macroparticle. Whenshot-noise effects are desired, each macroparticle’s posi-tion is shifted by a small dz according to well describedalgorithms [27, 28] to achieve the correct statistics.In addition, it is important to make sure that the noisefrom other coordinates would not contribute to a growthof the bunching as the beam propagates in the absenceof an interaction. This is taken care of by mirroring theenergy, transverse coordinates, and momenta over a largenumber of 5D phase space bins. The number of bins (typ-ically larger than 32) should be chosen such that bunch-ing in the absense of an interation remains small for allthe discrete frequencies included in the simulations. III. EXAMPLES
We limit this discussion to three examples that high-light the main features of our approach, even though itis expected that the new code can be successfully ap-plied to a variety of other situations. The first caseconsidered is a classical single-pass FEL seeded ampli-fier which will enable a quantitative comparison with thesemi-analytical M. Xie formulas [29] as well as with a tra-ditional period-average code like Genesis for both planarand helical geometries. The second example is relevantto the TESSA266 experiment being planned at the LEAbeamline at the APS linac in Argonne National Labo-ratory aiming at very high conversion efficiency at 266nm [20]. This case serves to illustrate the capability ofusing a 3D magnetic field map for a fairly complicatedsegmented tapered undulator. The code compares wellwith a traditional FEL code like Genesis, even deep inthe non-linear regime. The details of the beam transport(injection, entrance and exit sections and especially un-dulator break sections) can only be included in Genesisby using a linear beam transport approximation. GPTfollows the evolution of the beam distribution along thebeamline using field maps for all the magnetic elements(undulators, quadrupoles and phase shifter dipoles) andcalculates energy exchange using the self-consistent inter-action with the free-space modes. The results allow us to -400 -200 0 200 400 z(nm) G a mm a FIG. 2: left) Longitudinal phase space distribution with quiet loading for time-independent (i.e. single frequency)simulation. right) Longitudinal phase space distribution for multifrequency simulation. Particles are color coded bytheir charge weight. The projection onto the z-axis shows the Gaussian current profile.quantitatively include the effects of the entrance and exitsections (which add an effective 0.5 periods of interactionon each side of the undulator) and the trajectories afterthe prebuncher and in between the undulators.The final example is a waveguide THz FEL whereGPT-FEL is used to correctly simulate the zero-slippageamplification. In this configuration, the strong disper-sive properties of the guide affect the interaction whichtakes place in the zero-slippage regime. This scenariohighlights a unique capability of our code which wouldbe particularly challenging to simulate with traditionalFEL codes.
A. FEL amplifier
The parameters for this example are reported in Ta-ble I and somewhat arbitrarily chosen to be similar toan un-tapered version of the TESSA266 experiment dis-cussed below. The main differences are that a 200 periodlong undulator (with no break-section) is used for thisexample and the input seed power is lowered to 10 kW.An analytical model for the undulator magnetic field isused. The beam is transversely matched to the undula-tor natural focusing (equally distributed in the horizon-tal and vertical plane) so that its rms spot size remainsnearly constant along the interaction. The main goal ofthis example is to benchmark GPTFEL against the fit-ting formulas for the 3D gain length of an untapered FELamplifier and compare with a conventional FEL code likeGenesis. We also used this example to evaluate the per-formance of the single mode SDE approximation versus asimulation with n = 11 azimuthally symmetric Laguerre Gaussian SDE modes to decompose the electromagneticfield. GPTFEL took 1.5 minute to simulate 76800 par-ticles on an 8 processor for the single SDE mode and 5minutes for 11 SDE modes.TABLE I: Parameters for the 266 nm FEL amplifiersimulation. Electron BeamEnergy 375.5 MeVEnergy Spread 0.1 %RMS Bunch length 20 µ m (cid:15) n,x , (cid:15) n,y · mrad I peak σ x , σ y µ mRadiation λ
266 nmInput Power 10 kWRayleigh Length 1.41 mWaist location 0 mUndulator K rms λ u The time-independent, single frequency results for theplanar and helical geometries are shown in Figure 3 andcompared with Genesis 1.3. When using multiple spa-tial modes, the gain lengths in the planar and helicalcase are in good agreement (within 10 %) of the semi-analytical and numerical model predictions. The radia-tion spot sizes defined by σ r = (cid:82) r | E | d x (cid:82) | E | d x also closelyfollow the prediction. Note that while a single SDE modeis able to achieve qualitative results up to and near sat-uration, a larger number of spatial modes is required toFIG. 3: A comparison of GPTFEL running with SDE versus Genesis 1.3. a) The predicted gain length for theplanar amplifier is 0.287 m. Simulating with SDE and a single spatial mode overshoots by 16%. Running with 11SDE spatial modes reduces the error to 5.9%. b) The predicted gain length for the helical amplifier is 0.224 m.Simulating 1 and 11 SDE modes leads to errors of 15% and 8.2%, respectively.FIG. 4: GPTFEL results for 31 spatial frequencies, each with a single gaussian transverse mode. a) Waterfall plot ofnormalized power. b) Spectrum at P=0.1 GW for different thresholds on SDE interaction. ∆ is the ratio of L thresh to the theoretical gain length. Numerical errors occur when ∆ (cid:47) L thresh should be an order ofmagnitude larger than the theoretical gain length for convergent results.correctly simulate the evolution of the radiation profileafter saturation.The multi-frequency simulation used an SDE gaussianmode for 31 spatial frequencies with a 6% bandwidthto simulate 128,000 particles in 23 minutes. The user-defined parameter L thresh limits the spot size variationalong the undulator. Figure 4a shows a waterfall plotin the electron beam frame normalized at each z po-sition to display the relative velocity of the radiationwavepacket, which is close to the beam velocity in theexponential regime and becomes superluminal in the nonlinear regime [30]. In Figure 4b, the spectrum just beforesaturation is shown as a function of L thresh normalized to the gain length. If an increased spectral resolution isrequired, computation time scales linearly with numberof tracked modes. B. TESSA266
In this next example we take advantage of the GPTfunctions to track the electron beam in the fairly com-plex transport line of the TESSA 266 experiment. Thebeamline includes a short, 8 period undulator followed bya 3 dipole chicane to convert the imprinted energy mod-ulation into microbunching. Quadrupole doublets matchFIG. 5: Energy exchange and spotsizes in the first two tapered undulators of the TESSA beamline.FIG. 6: Bunching and Phase Space from the TESSA beamline.the beam transversely into the focusing channel of each0 .
96 meter, strongly tapered undulator section. A smalldipole is placed between the second quadrupole doubletso that the three magnets can be used as a phase shifterbetween the undulator sections.The GPT transport functions are used to set upthe trajectory and the beam optics prior of turning onthe seed and the FEL interaction module. Our time-independent simulation of the TESSA266 beamline in-cludes 21 higher order spatial modes to ensure an accu-rate modeling of the radiation profile. A 1 GW peakpower input radiation pulse is focused at the entrance ofthe tapered undulator to a waist of 0.3 mm. The simu-lation is compared with Genesis results, but it should benoted that GPTFEL uses full 3D magnetic field maps forthe undulators as well as for the dipoles and quadrupolesin the system. The magnetic field in the chicane dipoles isfine tuned to maximize the bunching and simultaneouslyoptimize the injection phase of the beamlets relative tothe radiation phase at the entrance of the tapered undu- lator. In Genesis, both the R56 and phase shifts are ap-plied post-facto to the beam distribution at the entranceof the tapered undulator, explaining the large differencein the bunching factor evolution in Fig. 6a. In practice,the phase shifter between the tapered undulator sectionshad to be re-optimized to account for the additional slip-page incurred by the beam when passing in the entranceand exit section of the wigglers. This is accomplishedby horizontally shifting the quadrupoles in opposite di-rections to steer the beam and tuning the magnetic fieldamplitude of the dipole to recover a straight trajectorywhile maximizing the energy exchange in the second un-dulator.
C. Zero slippage THz FEL
A final example to showcase the capabilities of the newGPTFEL code is the simulation of a THz FEL operatingin the zero-slippage regime [31]. The size of the waveg-FIG. 7: TESSA Beamlineuide is chosen in order to match the group velocity of theradiation with the electron beam longitudinal velocity in-side the undulator. This increases the bandwidth of theresonant interaction and extracts a significant amount ofenergy from very short electron beams.GPTFEL correctly simulates the waveguide dispersiveproperties as shown in Fig. 8 by plotting the electric fieldat the entrance and exit of the 1 meter long waveguidesystem in the absence of strong interaction (i.e. for verylow charge beams).The parameters of this example are summarized in Ta-ble II. We have chosen a planar undulator geometry withequally distributed focusing in the horizontal and verti-cal plane. In this case, the largest coupling is obtainedwith the TE10 mode profile of the curved parallel platewaveguide. The beam is initialized at the entrance of thesimulation with a large bunching factor (0.5) while weset the amplitude of the initial input seed to zero.There are two main advantages of using the waveguidein this system. First, the waveguide maintains a con-stant radiation cross section along the interaction, avoid-ing diffraction effects. Second the waveguide’s dispersiveproperties enable a zero slippage interaction. This largebandwidth interaction can drive the FEL with a muchshorter beam because the slippage effects are effectivelyminimized and the radiation continues to interact and ex-change energy with the particles even after a large num-ber of periods. The simulation results are shown in Fig.9 where the THz electric field waveform and the electronbeam longitudinal phase space are shown to be tempo-rally overlapping at the end of the undulator. Note thatthe system evolves in the non linear regime from the be- -0.01 -0.005 0 0.005 0.01-1.5-1-0.500.511.5 T H z e l e c t r i c fi e l d ( M V / m ) z (m) initialfinal FIG. 8: GPT time-dependent simulation temporal fieldprofile at the entrance and at the exit of the 1 m longwaveguide. The shift in the peak corresponds to thegroup velocity difference from the speed of light whichis matched to the electron beam longitudinal velocity inthe undulator in the zero-slippage regime. Helicalgeometry. Radiation spectrum and temporal profile ofthe pulse along the undulator. Final longitudinal phasespace.ginning as the electron beam enters the undulator witha very large bunching at the 1 THz resonant frequencyinduced by modulating the photocathode drive laser [32].The undulator is linearly tapered starting from its halfway point with a relative change in normalized vectorpotential K of 30 %/m to avoid saturation effects due toparticles falling off the resonance curve. The efficiency of0 Energy vs. z -200-1000100200
THz field
THz spectrum
Beam long phase space E n e r g y ( m J ) z(m)Frequency (THz) z(m)z(m) E l e c t r i c F i e l d ( M V / m ) G a mm a FIG. 9: a) THz pulse energy along the undulator. b) THz waveform at the undulator exit. c) THz spectrum. d)Longitudinal phase space of electron beam.TABLE II: Parameters for high efficiency THzamplifier.
Electron BeamEnergy 10.2 MeVEnergy Spread 1.25 %Bunch length 2000 µ m I peak
60 A (cid:15) n,x , (cid:15) n,y · mrad σ x , σ y µ mUndulator and waveguide K rms λ u b R conversion is above 10 % in this example. IV. CONCLUSIONS AND OUTLOOK
A new approach for FEL simulations has been pre-sented. The characteristic features are the decomposi-tion of the field in a set of spatial and frequency modesand the integration with the GPT numerical integrationengine which allows access and compatibility with a largenumber of beam transport designs and functions. Thereare a number of research opportunities which go beyondthe scope of this paper but will be the subject of futurestudies, including a detailed study of the effects of thetransverse space charge forces and an upgrade to includehigher harmonic interactions. Parallelization of the codewill allow much faster run times, increasing the number ofmacroparticles and modes that can be simulated. GPT-FEL is not expected to replace traditional approaches toFEL numerical simulations, but is intended to be a re-search tool to explore the interaction of relativistic elec-trons and electromagnetic waves in undulator systems in1regimes where the approximations of standard FEL codesare questionable. The application of GPTFEL to disper-sive systems allows for exploration of novel interactionregimes like the tapered waveguide THz FEL.
ACKNOWLEDGMENTS
The authors would like to thank Luca Giannessi andAvraham Gover for their helpful discussions. This workwas supported by DOE grant No. de-sc0009914. [1] B. W. McNeil and N. R. Thompson, Nature photonics ,814 (2010).[2] Z. Huang and K.-J. Kim, Physical Review Special Topics-Accelerators and Beams , 034801 (2007).[3] C. Pellegrini, A. Marinelli, and S. Reiche, Reviews ofModern Physics , 015006 (2016).[4] S. Biedron, Y. Chae, R. J. Dejus, B. Faatz, H. Fre-und, S. Milton, H.-D. Nuhn, and S. Reiche, Nuclear In-struments and Methods in Physics Research Section A:Accelerators, Spectrometers, Detectors and AssociatedEquipment , 110 (2000).[5] S. Reiche, Proceedings of FEL2010, Malm¨o, Sweden (JA-CoW, 2010), Vol. MOOCI1 , 165 (2010).[6] L. Giannessi, Physical Review Special Topics-Accelerators and Beams , 114802 (2003).[7] L. Giannessi, in Proceedings of the free-electron laser con-ference (2006).[8] C. Emma, N. Sudar, P. Musumeci, A. Urbanowicz, andC. Pellegrini, Physical Review Accelerators and Beams , 110701 (2017).[9] W. M. Fawley, A user manual for GINGER and itspost-processor XPLOTGIN , Tech. Rep. (Lawrence Berke-ley National Lab.(LBNL), Berkeley, CA (United States),2002).[10] S. Reiche, Nuclear Instruments and Methods in PhysicsResearch Section A: Accelerators, Spectrometers, Detec-tors and Associated Equipment , 243 (1999).[11] E. Saldin, E. Schneidmiller, and M. Yurkov, NuclearInstruments and Methods in Physics Research Section A:Accelerators, Spectrometers, Detectors and AssociatedEquipment , 233 (1999).[12] L. Campbell and B. McNeil, Physics of Plasmas ,093119 (2012).[13] H. Freund, P. van der Slot, et al. , in Proc. of 36th Int.Free-Electron Laser Conf.(Basel,) (2014) p. 408.[14] M. De Loos and S. Van Der Geer, in (1996) p. 1241.[15] I. V. Bazarov and C. K. Sinclair, Physical Review SpecialTopics-Accelerators and Beams , 034202 (2005).[16] A. Brynes, in (JACOW Pub-lishing, Geneva, Switzerland, 2019) pp. 578–583.[17] P. Sprangle, A. Ting, and C. Tang, Physical Review A , 2773 (1987). [18] P. Baxevanis, R. D. Ruth, and Z. Huang, Physical Re-view Special Topics-Accelerators and Beams , 010705(2013).[19] A. Gover, R. Ianconescu, A. Friedman, C. Emma, N. Su-dar, P. Musumeci, and C. Pellegrini, Reviews of ModernPhysics , 035003 (2019).[20] Y. Park, R. Agustsson, T. Campese, D. Dang, P. Den-ham, I. Gadjev, C. Hall, A. Murokh, P. Musumeci, N. Su-dar, et al. , in (JACOW Pub-lishing, Geneva, Switzerland, 2019) pp. 730–733.[21] E. Curry, S. Fabbri, P. Musumeci, and A. Gover, NewJournal of Physics , 113045 (2016).[22] M. De Loos, C. van der Geer, S. Van der Geer, A. van derMeer, D. Oepts, and R. W¨unsch, Nuclear Instrumentsand Methods in Physics Research Section A: Acceler-ators, Spectrometers, Detectors and Associated Equip-ment , 97 (2003).[23] T. Nakahara and N. Kurauchi, IEEE Transactions onMicrowave Theory and Techniques , 66 (1967).[24] E. J. C. Snively, Electron-THz Wave Interactions in aGuided Inverse Free Electron Laser (University of Cali-fornia, Los Angeles, 2018).[25] H. Freund, L. Giannessi, and W. Miner Jr, Journal ofApplied Physics , 123114 (2008).[26] B. McNeil, M. Poole, and G. Robb, Physical ReviewSpecial Topics-Accelerators and Beams , 070701 (2003).[27] W. M. Fawley, Physical Review Special Topics-Accelerators and Beams , 070701 (2002).[28] C. Penman and B. McNeil, Optics communications ,82 (1992).[29] M. Xie, Nuclear Instruments and Methods in Physics Re-search Section A: Accelerators, Spectrometers, Detectorsand Associated Equipment , 59 (2000).[30] X. Yang, N. Mirian, and L. Giannessi, Physical ReviewAccelerators and Beams , 010703 (2020).[31] E. Snively, J. Xiong, P. Musumeci, and A. Gover, Opticsexpress , 20221 (2019).[32] P. Musumeci, R. Li, and A. Marinelli, Physical reviewletters106