**1. Introduction**

Free-electron laser (FEL) oscillators (FELO) have been part of the overall research activity since the beginnings of the field. An FEL oscillator consists of an undulator placed within an optical resonator, also known as an optical cavity, that typically consists of two spherical mirrors separated by a distance *L*cav. Note that the undulator is not necessarily centered within the optical resonator. Furthermore, an electron beam is directed into the cavity through the undulator, and out of the cavity, as is schematically shown in Figure 1, which is taken from Reference [1]. The optical pulse generated by the electrons streaming through the undulator is extracted from the cavity either by a partially reflective (i.e., transmissive) mirror or through a hole in the mirror. Observe that multiple optical pulses may be contained within the resonator depending upon the repetition rate of the electron bunches and the roundtrip time for the optical pulses to circulate through the resonator. We also note that, when the growth of the optical field exponentiates in a single pass through the undulator, i.e., in case of a high-gain FEL [2], a large fraction of the light can be extracted from the resonator.

The first FEL oscillator driven by a relativistic electron beam was realized in 1977 [3], shortly after the coherent amplification of radiation of the same wavelength was experimentally demonstrated [4]. Coherent undulator radiation was first demonstrated using the so-called Ubitron, which was configured both as amplifier and oscillator, and which is nowadays considered to be the first (non-relativistic) FEL [5]. In the early days of FEL development, the technology available to generate high-quality, high-brightness electron beams, that are necessary for efficient FEL operation, was relatively primitive by today's

**Citation:** van der Slot, P.J.M.; Freund, H.P. Three-Dimensional, Time-Dependent Analysis of Highand Low-*Q* Free-Electron Laser Oscillators. *Appl. Sci.* **2021**, *11*, 4978. https://doi.org/ 10.3390/app11114978

Academic Editor: Giuseppe Dattoli

Received: 6 May 2021 Accepted: 24 May 2021 Published: 28 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

standards. At that time, the limited available gain, which in part may be due to space restrictions, favoured oscillator configurations to bring the laser into saturation and provide fully coherent output, as mirrors were readily available all the way from the visible to well into the microwave spectral range. Although we limit ourselves here to infrared wavelengths or shorter, exciting FEL research also takes place at longer wavelengths.

**Figure 1.** Schematic representation of a typical layout for an FEL oscillator. The electron beam consisting of electron pulses is guided through the undulator and around the mirrors. The optical pulse that is amplified within the undulator is recirculating within the resonator and extracted through one of the mirrors. Reproduced with permission from [1]. Copyright Springer, 2018.

The unique properties of FEL oscillators, in particular to fill the spectral gap for coherent light sources in the visible to microwave spectral range, resulted quickly in the establishment of several facilities for applied research. The characteristics of user facilities using FEL oscillators depends, in part, on the accelerator technology used, which includes electrostatic acceleration [6–9], energy recovery linear accelerators (linacs) [10–15], nonrecovery linacs [16–30], microtrons [31,32] and synchrotrons [33–37]. Research performed at these user facilities includes, amongs<sup>t</sup> others, biological and material science [38]. A more complete overview of the various FEL oscillators, together with the main system parameters, can be found in Ref. [39].

As mentioned, the optical properties of FEL oscillators are in part determined by the electron accelerator used. For example, electrostatic accelerators can only accelerate electrons up to a kinetic energy of a few MeV and, consequently, the spectral range is limited to microwave and THz frequencies [39]. The much higher electron energies of storage rings allow the FEL oscillators to operate in the visible down to the deep UV wavelengths [39], while RF linacs can, in principle, accelerate electrons to still higher energies, from a few MeV up to tens of GeV, allowing FELs to cover a large fraction of the electromagnetic spectrum, from microwaves all the way down to hard X-rays [39].

Furthermore, electrostatic accelerators are capable of long to continuous wave electron pulses, that can result in single longitudinal mode oscillation, with a relative linewidth of about 10−<sup>8</sup> over a time of 5 μs, which was set by the voltage droop present in the accelerator [6] at the time of measurement. On the other hand, the short duration of, and spread in, electron energy for the electron pulses produced by RF linacs typically leads to a relative linewidth of the order of 10−<sup>4</sup> to 10−<sup>2</sup> [13,29,40–53], unless additional spectral filtering is used. For example, using the Vernier effect of a double cavity resonator to reduce the number of intra-cavity longitudinal modes and subsequent external filtering, an RF-linac-based FEL oscillator was able to provide a single longitudinal mode [54]. FEL oscillators based on storage rings provide a lower relative linewidth within the range 10−<sup>5</sup> to 5 × 10−<sup>4</sup> [34,55–63], due to the higher electron beam quality of storage rings.

Finally, the average and peak power are, to a large extent, determined by the accelerator technology. Superconducting energy recovery linacs are the primary drivers for high average power FEL oscillators [64–70]; however, other types of FEL oscillators have used various techniques to increase the output pulse energy including, but not limited to, optical klystrons [34,36,48,57,58,60,63,71–76], cavity dumping [77–81], pulse stacking in an external cavity [82–84], electron beam energy ramping [85,86] or dynamic cavity desynchronization [87–92]. Tapering of the undulator [93–95] has also been used to improve the

peak power in the pulse [96–100], however, tapering of the undulator in FEL oscillators leads to complicated dynamics and a strong dependence on system parameters [101–104].

The wavelength tuning of FEL oscillators is typically obtained by selecting a set of electron beam energies and varying the undulator gap. Operation at other wavelength ranges is possible through harmonic generation [105–107] and such harmonics were quickly observed [40]. Lasing on the third harmonic was first demonstrated using the Standford MARK-III FEL oscillator [108], quickly followed by others [43,48,58,61,97,109–121]. The simultaneous generation of multiple wavelengths may be beneficial for certain applications. The fixed relation between the fundamental and its harmonics can be overcome using different undulators in a single optical cavity [19,122–125], two undulators driven by a single electron beam at a single beam energy, but each having its own optical resonator [126] or having different beam energies [127,128] to allow operation at two independently chosen wavelengths.

FEL oscillators share many characteristics with amplifier systems, however, there are also many differences. For example, slippage between the optical and electron pulse, in combination with cavity desynchronization, will effect the FEL oscillator dynamics. Here, a cavity is considered synchronized when it is tuned to a length where the roundtrip time for the optical pulse matches an integer multiple of the time separating two subsequent electron pulses. Maxima in the gain and extraction efficiency are found for different cavity desynchronizations. In particular, dynamic cavity desynchronization [87] can be used to first set a cavity length corresponding to maximum gain for a quick initial growth of the optical pulse energy and then switched to a cavity length for maximum extraction efficiency that is obtained at small cavity desynchronization. The latter is typically also associated with the generation of short optical pulses [19,88,92,129–134]. For short pulse FEL oscillators, i.e., where the root-mean-square (rms) width of the electron bunch is smaller than the total slippage distance, the oscillator is operating in the superradiant regime [132,134–142]. An analytic theory for this regime predicts a scaling of the optical pulse energy *EL* with the bunch charge *q* as *EL* ∝ *q*3/2 and a scaling of the temporal width *τL* of the optical pulse as *τL* ∝ *q*<sup>−</sup>1/2 [138], which has been experimentally confirmed [139], although higher extraction efficiencies have also been observed for a perfectly synchronized cavity length, where the analytic theory breaks down [92,132]. Furthermore, short pulse FEL oscillators operating in the large-slippage regime also show a stable oscillation in the macropulse power, known as limit-cycle oscillations, for certain cavity lengths detuned from perfect synchronization, which has been observed both experimentally [102,137,143–146] and in numerical simulations [136,138,147–149].

Techniques well known from laser physics have also been applied to FEL oscillators, either experimentally or in a conceptual study, such as injection seeding [150–155] and mode locking [133,151,152,156]. On the other hand, the unique nature of the gain process in FELs allows for the micro-structuring of the gain medium [157–159] to enhance the gain, which is not readily available for ordinary lasers. One important application of this technique is to transfer coherence from a low-frequency optical beam to a higher frequency optical beam. This process is known as high-gain harmonic generation (HGHG) [160]. This technique has been successfully used to generate fully coherent light at soft X-ray frequencies [161]. Alternative methods currently investigated, in particular to reach even shorter wavelengths, and rely on Bragg reflections from atomic layers using appropriated crystals [162].

The drive to higher gain per unit length and to shorter wavelengths has pushed accelerator development; in particular, research into photo-cathode based electron sources [163,164] that typically produce higher quality electron beams has contributed to today's state-ofthe-art accelerator technology. The availability of electron beams with improved beam quality has driven FEL development over the years, especially the capability to generate shorter wavelengths. Due to a lack of suitable mirrors, initial focus was on FELs relying on self-amplified spontaneous emission (SASE), with a proof-of-principle experiment at 530 nm [165], followed by demonstrations at shorter wavelengths [166]. To overcome poor

temporal coherence, HGHG is used to provide completely coherent output pulses down to 4-nm wavelength [167], while for even smaller wavelengths, FEL oscillators using the above-mentioned Bragg mirrors are considered [162].

Along with increasingly sophisticated FEL oscillator designs and experiments, analytical theories [101,138,168–171] have been developed for an initial assessment of oscillator performance. An initial theoretical description of the FEL gain was quantum mechanical in nature [172,173], however, it was quickly demonstrated that the gain could be described using classical theory [174]. Furthermore, the more complete simulation codes have also become adept at treating increasingly complicated designs and include both steady-state and time-dependent simulations in one- and three-dimensions.

In the remainder of this paper, we will focus on the simulation of a few different FEL oscillators. To model these oscillators, we have selected the FEL gain code MIN-ERVA, which models the light amplification within the undulator line, and the optical wavefront propagation code OPC, which models the light propagation in the remainder of the resonator. The mathematical formulations of these codes are presented in Section 2. Subsequently, we present some considerations on resonators relevant for FEL oscillators in Section 3. We continue with a discussion on a low-gain/high-*Q* oscillators, taking the Infrared Demo and Infrared Upgrade experiments performed at the Thomas Jefferson National National Accelerator Facility (JLab) as examples. Both are discussed in Section 4. We compare this with a high-gain/low-*Q* oscillator, also known as a regenerative amplifier (RAFEL), also operating in the infrared, and investigate the performance of such a system in the soft X-ray spectral range. These systems are described in Section 5. We conclude with a summary and discussion in Section 6.

#### **2. Numerical Formulation**

Simulations of FEL oscillators have appeared in the literature by including optical propagation algorithms to existing FEL simulation codes, as well as self-contained FEL oscillator codes [103,149,175–182]. In this paper, we describe the use of an existing FEL simulation code to treat the interaction within the undulators and to link that to a code specifically designed to propagate the optical field through various resonator configurations. Such optics codes are sufficiently general to be able to treat a variety of resonator designs, and many may have been originally created to deal with other types of lasers. In this simulation environment, the FEL code hands off the optical field at the resonator exit to the optics code, which then propagates the field around the resonator and back to the undulator entrance, after which it is handed off to the FEL code for another pass through the undulator.

Among the earliest applications of this method was to simulations of the IR-Demo [64,65], as described in Section 4.1, and the 10-kW Upgrade experiments at JLab [183], which employed the MEDUSA FEL code [184] and OPC [185] and successfully reproduced many of the essential features of the experiment. OPC has also been interfaced with the GENESIS FEL code [186] and the PUFFIN FEL code [187] to study VUV FEL oscillators [188,189]. In this paper, we discuss the linkage of the MINERVA simulation code with OPC.

As will be described below, while MINERVA represents the optical field as a superposition of Gaussian modes with two independent polarization directions, OPC propagates the field on a grid. This necessitates the mapping of the Gaussian modes in MINERVA at the undulator exit to the grid supported by OPC, and the decomposition of the optical field in OPC at the undulator entrance back into Gaussian modes for both polarization components.

#### *2.1. The MINERVA Simulation Code*

MINERVA is a time-dependent, three-dimensional simulation code that models the interaction between electrons and a co-propagating optical field through an undulator line, which may include strong-focusing quadrupoles and/or dipole chicanes [1,190–192]. The optical fields are described using a superposition of Gauss-Hermite modes, using two independent polarization directions taken to be in the *x*- and *y*-directions with the *z*-direction taken along the axis of the undulator line. The Gauss–Hermite modes constitute

a complete set, and together with the two independent polarizations, are able to describe the generation of arbitrary polarizations that might arise for any given undulator geometry. The vector potential representing the optical field in terms of the Gauss–Hermite modes is given by

$$\begin{split} \delta \mathbf{A}(\mathbf{x},t) &= \mathfrak{e}\_{\mathbf{x}} \sum\_{\substack{l,n=0\\h\neq 1}}^{\infty} \mathfrak{e}\_{l,n,h}^{(\mathbf{x})}(\mathbf{x},t) \Big( \delta A\_{l,n,h}^{(1,\mathbf{x})}(\mathbf{z},t) \sin(\mathfrak{q}\_{h}^{(\mathbf{x})}(\mathbf{x},t)) + \delta A\_{l,n,h}^{(2,\mathbf{x})}(\mathbf{z},t) \cos(\mathfrak{q}\_{h}^{(\mathbf{x})}(\mathbf{x},t)) \Big) + \\ & \quad \mathfrak{e}\_{\mathbf{y}} \sum\_{\substack{l,n=0\\h\neq 1}}^{\infty} \mathfrak{e}\_{l,n,h}^{(\mathbf{y})}(\mathbf{x},t) \Big( \delta A\_{l,n,h}^{(1,\mathbf{y})}(\mathbf{z},t) \sin(\mathfrak{q}\_{h}^{(\mathbf{y})}(\mathbf{x},t)) + \delta A\_{l,n,h}^{(2,\mathbf{y})}(\mathbf{z},t) \cos(\mathfrak{q}\_{h}^{(\mathbf{y})}(\mathbf{x},t)) \Big), \end{split} \tag{1}$$

where (*l*, *n*) are the transverse mode numbers, *h* is the harmonic number and *<sup>δ</sup>Al*,*n*,*<sup>h</sup>* are the mode amplitudes. Furthermore, **e**ˆ*j* is a unit vector in the *j*-direction, *j* = *x*, *y* and

$$w\_{l, \mathbf{u}, \mathbf{h}}^{(j)}(\mathbf{x}, t) = \frac{w\_{0, \mathbf{h}}^{(j)}}{w\_{\mathbf{h}}^{(j)}(z - z\_{0}, t)} e^{-r^{2}/(w\_{\mathbf{h}}^{(j)}(z - z\_{0}, t))^{2}} H\_{l}(\zeta\_{\mathbf{x}}^{(j)}(\mathbf{x}, z, t)) H\_{\mathbf{n}}(\zeta\_{\mathbf{y}}^{(j)}(y, z, t)) \tag{2}$$

are the transverse eigenfunctions for polarization direction *j r* = *x*<sup>2</sup> + *y*2, *Hl* and *Hn* are the Hermite polynomials of order *l* and *n*, respectively, *ζ*(*j*) *x* (*<sup>x</sup>*, *z*, *t*) = √<sup>2</sup>*x*/*w*(*j*) *h* (*z* − *z*0, *t*), *ζ*(*j*) *y* (*y*, *z*, *t*) = √<sup>2</sup>*y*/*w*(*j*) *h* (*z* − *z*0, *t*) and *w*(*j*) 0,*h* and *w*(*j*)(*z* − *z*0, *t*) are the beam radii in the waist at *z* = *z*0 and at location *z* along the optical axis. Note, *e*(*j*) *l*,*n*,*h* and *w*(*j*) are dependent on time through the source dependent expansion (SDE) [193], as different parts of the optical pulse will experience different gain. Finally, the optical phases *ϕ*(*j*) *h* (**<sup>x</sup>**, *t*) are given by

$$\varphi\_{\rm li}^{(j)}\left(\mathbf{x},t\right) = h(k\_0 z - \omega\_0 t) + a\_{\rm li}^{(j)}\left(z,t\right) \left(\frac{r}{w\_{\rm li}^{(j)}\left(z - z\_{0\prime} t\right)}\right)^2,\tag{3}$$

where *k*0 = *ω*0/*c* is the vacuum wavenumber belonging to the carrier frequency *ω*0, *c* is the speed of light in vacuo and *α*(*j*) *h* (*z*, *t*) is the curvature of the wavefront for harmonic *h*. In this formulation, we assume that the mode amplitudes, curvature of the phase fronts and radii of the modes are slowly varying functions of *z* and *t*. The total optical power is given by

$$P(z,t) = \sum\_{\substack{l,n=0\\h=1}}^{\infty} P\_{l,n,h}(z,t) = \frac{m\_c^2 c^5}{8\epsilon^2} k\_0^2 \sum\_{\substack{l,n=0\\h=1}}^{\infty} 2^{l+n-1} l! n! \sum\_{j=x,y} (w\_{0,h}^{(j)})^2 \left( (\delta a\_{l,n,h}^{(1,j)}(z,t))^2 + (\delta a\_{l,n,h}^{(2,j)}(z,t))^2 \right), \tag{4}$$

where *<sup>δ</sup>a*(*<sup>i</sup>*,*j*) *<sup>l</sup>*,*n*,*<sup>h</sup>*(*<sup>z</sup>*, *t*) = *eδA*(*<sup>i</sup>*,*j*) *<sup>l</sup>*,*n*,*<sup>h</sup>*(*<sup>z</sup>*, *<sup>t</sup>*)/*mec*<sup>2</sup> are the components of the normalized mode amplitudes, *i* = 1, 2, *j* = *x*, *y* and *me* is the electron rest mass.

The optical fields are driven by accelerated electrons, as they co-propagate through the magnetic fields that are placed along the transport line. Of main interest here is the amplification of the optical field by the electrons within the undulators that can be placed in arbitrary configuration along the electron beam transport line. Within the slowly varying phase and amplitude approximation, the evolution of the normalized mode amplitudes *<sup>δ</sup>a*(*<sup>i</sup>*,*j*) *<sup>l</sup>*,*n*,*<sup>h</sup>*(*<sup>z</sup>*, *t*) is given by

$$\frac{d}{dz}\begin{pmatrix}\delta a\_{l,n,h}^{(1,j)}(z,t)\\\delta a\_{l,n,h}^{(2,j)}(z,t)\end{pmatrix} + K\_{l,n,h}^{(j)}(z)\begin{pmatrix}\delta a\_{l,n,h}^{(2,j)}(z,t)\\-\delta a\_{l,n,h}^{(1,j)}(z,t)\end{pmatrix} = \begin{pmatrix}S\_{l,n,h}^{(1,j)}(z,t)\\S\_{l,n,h}^{(2,j)}(z,t)\end{pmatrix},\tag{5}$$

where *d*/*dz* is the convective derivative and *K*(*j*) *<sup>l</sup>*,*n*,*<sup>h</sup>*(*<sup>z</sup>*, *t*) is given by

$$K\_{l,n,h}^{(j)}(z,t) = (1+l+n)\left[\frac{a\_h^{(j)}(z,t)}{w\_h^{(j)}(z,t)}\frac{dw\_h^{(j)}(z,t)}{dz} - \frac{1}{2}\frac{da\_h^{(j)}(z,t)}{dz} - \frac{1+(a\_h^{(j)}(z,t))^2}{k\_0(w\_h^{(j)}(z,t))^2}\right].\tag{6}$$

The source terms *<sup>S</sup>*(*<sup>i</sup>*,*j*) *<sup>l</sup>*,*n*,*<sup>h</sup>*(*<sup>z</sup>*, *t*) in Equation (5) are given by

$$
\begin{pmatrix} S\_{l,n,h}^{(1,j)}(z,t) \\ S\_{l,n,h}^{(2,j)}(z,t) \end{pmatrix} = \frac{1}{\pi} \frac{\omega\_b^2}{k\_0 c^2} \frac{1}{2^{l+n-1}} \frac{1}{l! n!} \frac{1}{(w\_{0,h}^{(j)})^2} \left\langle e\_{l,n,h}^{(j)}(\mathbf{x}) \frac{v\_j(z,t)}{|v\_z(z,t)|} \left(\frac{-\cos(\varrho\_h^{(j)}(z,t))}{\sin(\varrho\_h^{(j)}(z,t))}\right) \right\rangle,\tag{7}
$$

where *ωb* is the nominal plasma frequency, *vj* is the component of the electron velocity in direction *j*, *vz* is the longitudinal velocity component and (...) is an average over the electron distribution given by

$$\begin{aligned} \langle\langle\dots\rangle\rangle &= \\ \int\_0^{2\pi} \frac{d\psi\_0}{2\pi} \int\_1^\infty \frac{d\gamma\_0}{\sqrt{\pi/2}\Lambda\gamma} \frac{e^{-(\gamma\_0-\gamma\_{\rm avg})^2/2\Lambda\gamma^2}}{1+\text{erf}(\gamma\_{\rm avg}/\sqrt{2}\Lambda\gamma)} \iint \frac{d\mathbf{x}\_0 d\mathbf{y}\_0}{2\pi\sigma\_r^2} \iint \frac{d\mathbf{p}\_0 d\mathbf{p}\_{\mathbf{y}0}}{2\pi\sigma\_p^2} e^{-r^2/2\sigma\_r^2} e^{-\frac{\mathbf{p}\_{\perp}^2}{\mathbf{p}\_{\perp}\mathbf{p}}/2\sigma\_p^2} \langle\dots\rangle. \end{aligned} \tag{8}$$

In Equation (8), *γ*avg and Δ*γ* are the relativistic factors corresponding to the initial average electron energy and energy spread, and *σr* and *<sup>σ</sup>p* describe the width of the initial distribution function in transverse and momentum phase space, respectively.

The particle distribution described by Equation (8) is supplemented in the simulation code MINERVA, with methods to import 6-dimensional particle distributions which may be obtained from various sources, such as electron beam transport models, and this allows full start-to-end simulation capabilities.

To minimize the number of modes required to accurately describe the amplified optical field within the undulators, the so-called source dependent expansion [193] is used. Within this approximation, the spot size and curvature of the eigenmodes for each of the polarization directions *j* = *x*, *y* are allowed to evolve according to

$$\frac{dw\_h^{(j)}(z,t)}{dz} = \frac{2a\_h^{(j)}(z,t)}{hk\_0 w\_h^{(j)}(z,t)} - w\_h^{(j)}(z,t) \mathcal{Y}\_h^{(j)}(z,t) \tag{9}$$

and

$$\frac{1}{2}\frac{da\_h^{(j)}(z,t)}{dz} = \frac{1 + (a\_h^{(j)}(z,t))^2}{hk\_0(w\_h^{(j)}(z,t))^2} - X\_h^{(j)}(z,t) - a\_h^{(j)}(z,t)Y\_h^{(j)}(z,t),\tag{10}$$

where

$$\begin{aligned} X\_h^{(j)}(z,t) &= \\ \frac{2}{(\delta a\_{0,0}^{(j)}(z,t))^2} \Big[ \left( S\_{2,0,h}^{(1,j)}(z,t) + S\_{0,2,h}^{(1,j)}(z,t) \right) \delta a\_{0,0,h}^{(1,j)}(z,t) - \left( S\_{2,0,h}^{(2,j)}(z,t) + S\_{0,2,h}^{(2,j)}(z,t) \right) \delta a\_{0,0,h}^{(2,j)}(z,t) \Big] \end{aligned} \tag{11}$$

and

$$\begin{split} Y\_{h}^{(j)}(z,t) &= \\ & - \frac{2}{(\delta a\_{0,0}^{(j)})^2} \Big[ \Big( \mathcal{S}\_{2,0,h}^{(1,j)}(z,t) + \mathcal{S}\_{0,2,h}^{(1,j)}(z,t) \Big) \delta a\_{0,0,h}^{(1,j)}(z,t) + \Big( \mathcal{S}\_{2,0,h}^{(2,j)}(z,t) + \mathcal{S}\_{0,2,h}^{(2,j)}(z,t) \Big) \delta a\_{0,0,h}^{(2,j)}(z,t) \Big]. \end{split} \tag{12}$$

Furthermore, in Equations (11) and (12), (*a*(*j*) 0,0,*<sup>h</sup>*(*<sup>z</sup>*, *t*))<sup>2</sup> = (*a*(1,*j*) 0,0,*<sup>h</sup>*(*<sup>z</sup>*, *t*))<sup>2</sup> + (*a*(2,*j*) 0,0,*<sup>h</sup>*(*<sup>z</sup>*, *t*))2. Note that SDE (Equations (9)–(12)) recovers vacuum diffraction when no electron beam is present [*X*(*j*) *h*(*z*, *t*) = *Y*(*j*) *h*(*z*, *t*) ≡ 0].

The source terms in Equation (5) depend on the evolution of the electron coordinates and velocities that are given by the Newton–Lorentz equations

$$\frac{d\mathbf{x}}{dz} = \frac{v\_x}{v\_z} \, \tag{13}$$

$$\frac{dy}{dz} = \frac{v\_y}{v\_z} \,\tag{14}$$

$$\frac{d\mathbf{p}}{dz} = -\frac{\varepsilon}{v\_z} \left[ \delta \mathbf{E} + \frac{\mathbf{v}}{c} \times (\mathbf{B} + \delta \mathbf{B}) \right],\tag{15}$$

and the evolution of the ponderomotive phase *ψ* given by

$$\frac{d\psi}{dz} = k + k\_u - \frac{\omega}{v\_z}.\tag{16}$$

Using the Coulomb gauge, the optical fields *δ***E** = −1*c ∂δ***A***∂t* and *δ***B** = ∇ × *δ***A** are derived from the optical vector potential given by Equation (1), while **B** is the static magnetic field produced by the magnetic elements along the electron transport line. The model presented here ignores space-charge effects, however, this can be easily incorporated [1].

The dynamical equations for the particles and fields are integrated simultaneously using a 4th order Runge–Kutta algorithm. Hence, the number of equations in the simulation is *<sup>N</sup>*equations = *<sup>N</sup>*slices[6*N*particles + <sup>4</sup>(*<sup>N</sup>*modes + *<sup>N</sup>*harmonics)], where *N*slices is the number of slices in the simulation, and *<sup>N</sup>*particles is the number of particles in each slice, *N*modes is the total number of modes in the fundamental and all the harmonics, and *N*harmonics is the number of harmonics. The Runge-Kutta algorithm allows the step size to change so optimized step sizes can be used in each magnetic element or in the drift spaces while this imposes no limitation on the placement of components along the electron beam path.

The formulation self-consistently tracks the generation of the optical field with arbitrary polarizations depending on the undulator configuration. The polarization state of the output light, therefore, can be determined by calculation of the Stokes parameters [194]. Measurements of the polarization in FELs have been characterized by the Stokes parameters [195,196].

Magnetic field elements, such as undulators, quadrupoles and dipoles, can be placed in arbitrary sequences to specify a variety of different transport lines, and the gap lengths between different undulators can be varied as well. Field configurations can be set up for single or multiple undulator segments and can contain quadrupoles placed between the segments or be superimposed on the lattice to create a FODO lattice. Magnetic elements can also be dipole chicanes to model optical klystron, high-gain harmonic generation (HGHG) or phase shifters. See Section 2.2 for analytical models for the various static magnetic field elements.

#### *2.2. Analytical Models for Static Magnetic Fields*

The various types of undulators are modeled using three-dimensional analytical representations. Two models are available for the planar undulator, representing a flat-pole face and a parabolic-pole face with weak two-plane focusing, respectively. The flat-poleface undulator with the magnetic field oscillating in the *y*-plane is described by

$$\mathbf{B}\_{\mathsf{u}}(\mathbf{x}) = B\_{\mathsf{u}}(z) \left[ \left( \sin(k\_{\mathsf{u}}z) - \frac{\cos(k\_{\mathsf{u}}z)}{k\_{\mathsf{u}}B\_{\mathsf{u}}(z)} \frac{dB\_{\mathsf{u}}}{dz} \right) \mathbf{\hat{e}}\_{\mathsf{y}} \cosh(k\_{\mathsf{u}}y) + \mathbf{\hat{e}}\_{z} \sinh(k\_{\mathsf{u}}y) \cos(k\_{\mathsf{u}}z) \right], \tag{17}$$

where *Bu* is the amplitude of the undulator field, *ku* = 2*π*/*λ<sup>u</sup>*, *λu* is the undulator period and *Bu*(*z*) describes the entrance and exits taper at the ends of an undulator segment. The field described in Equation (17) is both curl- and divergence-free when the amplitude is constant. The analytic entrance and exit taper function *Bu*(*z*) is given by

$$B\_{\rm u}(z) = \begin{cases} B\_{\rm u0} \sin^2 \left( \frac{k\_{\rm tr} z}{4 \text{N}\_{\rm tr}} \right) & 0 \le z \le N\_{\rm tr} \lambda\_{\rm u} \\ B\_{\rm u0} & N\_{\rm tr} \lambda\_{\rm u} \le z \le L\_{\rm tr} \\ B\_{\rm u0} \cos^2 \left( \frac{k\_{\rm u} (L\_{\rm tr} - z)}{4 \text{N}\_{\rm tr}} \right) & L\_{\rm tr} \le z \le L\_{\rm u} \end{cases} \tag{18}$$

where *Bu*0 is the undulator amplitude in the homogeneous part of the undulator, *Lw* is the length of the undulator segment, *Ntr* is the number of undulator periods in the transition region and *Ltr* = *Lw* − *Ntrλu* is the location of the start of the exit taper. The field in the taper regions has zero divergence, and the *z*-component of the curl also vanishes. The transverse components of the curl do not vanish, but are of the order of 1/(*kuBu*(*z*))*dBu*/*dz*, which is usually small.

A parabolic-pole-face undulator with the magnetic field orientation in the *y*-plane and weak focusing in the *x*- and *y*-planes is described by the analytical model

$$\mathbf{B}\_{\rm u}(\mathbf{x}) = B\_{\rm u}(z) \left[ \left( \cos(k\_{\rm u} z) - \frac{\sin(k\_{\rm u} z)}{k\_{\rm u} B\_{\rm u}(z)} \frac{dB\_{\rm u}}{dz} \right) \mathbf{\dot{e}}\_{\perp}(\mathbf{x}, y) - \sqrt{2} \mathbf{\dot{e}}\_{z} \cosh\left(\frac{k\_{\rm u} x}{\sqrt{2}}\right) \sinh\left(\frac{k\_{\rm u} y}{\sqrt{2}}\right) \sin(k\_{\rm u} z) \right], \tag{19}$$

where

$$\mathfrak{e}\_{\perp}(x,y) = \mathfrak{e}\_{x}\sinh\left(\frac{k\_{\rm u}x}{\sqrt{2}}\right)\sinh\left(\frac{k\_{\rm u}y}{\sqrt{2}}\right) + \mathfrak{e}\_{y}\cosh\left(\frac{k\_{\rm u}x}{\sqrt{2}}\right)\cosh\left(\frac{k\_{\rm u}y}{\sqrt{2}}\right) \tag{20}$$

and the taper function *Bu*(*z*) is again given by Equation (18). As in the case of the flat-poleface model, Equation (19) is divergence free and the *z*-component of the curl also vanishes. Both, the flat-pole-face and parabolic-pole-face ideally produce linearly polarized light.

Circular polarized light is produced by helical undulators that are described by

$$\begin{split} \mathbf{B}\_{\rm u}(\mathbf{x}) &= 2B\_{\rm u}(z) \left[ \left( \cos(\chi) - \frac{\sin(\chi)}{k\_{\rm u}B\_{\rm u}(z)} \frac{dB\_{\rm u}}{dz} \right) I\_{1}^{'}(k\_{\rm u}r) \mathbf{\hat{e}}\_{r} - \\ & \left( \sin(\chi) - \frac{\cos(\chi)}{k\_{\rm u}B\_{\rm u}(z)} \frac{dB\_{\rm u}}{dz} \right) \frac{I\_{1}(k\_{\rm u}r)}{k\_{\rm u}r} \mathbf{\hat{e}}\_{\theta} + I\_{1}(k\_{\rm u}r) \sin(\chi) \mathbf{\hat{e}}\_{z} \right] \end{split} \tag{21}$$

in cylindrical coordinates *r*, *θ*, *z*, where *χ* = *kuz* − *θ*, *I*1 is the regular Bessel function of the first kind, the prime () indicates the derivative of the function with respect to its argumen<sup>t</sup> and *Bu*(*z*) is again given by Equation (18).

Finally, a magnetic field with varying degrees of ellipticity is produced by an APPLE-II undulator, which can be approximated by a superposition of two flat-pole-face undulators that are oriented perpendicularly to each other, and one of which can be displaced with respect to the other along the axis of symmetry. As such, the field is represented by

$$\mathbf{B}\_{\rm ll}(\mathbf{x}) = B\_{\rm ll}(z) \left[ \left( \sin(k\_{\rm u}z + \phi) - \frac{\cos(k\_{\rm u}z + \phi)}{k\_{\rm u}B\_{\rm u}(z)} \frac{dB\_{\rm u}}{dz} \right) \cosh(k\_{\rm u}x) \mathbf{\dot{e}}\_{x} + \right. \tag{22}$$

$$\left( \sin(k\_{\rm u}z) - \frac{\cos(k\_{\rm u}z)}{k\_{\rm u}B\_{\rm u}(z)} \frac{dB\_{\rm u}}{dz} \right) \cosh(k\_{\rm u}y) \mathbf{\dot{e}}\_{y} + \tag{22}$$

$$\left( \sinh(k\_{\rm u}x)\cos(k\_{\rm u}z + \phi) + \sinh(k\_{\rm u}y)\cos(k\_{\rm u}z) \right) \mathbf{\dot{e}}\_{z} \Big]\_{}$$

where *φ*, 0 ≤ *φ* ≤ *π*, represents the phase shift between the two linear undulators and as before, *Bu*(*z*) is given by Equation (18). This model is valid near the axis the axis of symmetry. The ellipticity, *ue* of the APPLE-II undulator is given by

$$\mu\_{\varepsilon} = \frac{1 - |\cos(\phi)|}{1 + |\cos(\phi)|}. \tag{23}$$

We note that the APPLE-II undulator configured as linear undulator has *ue* = 0, while *ue* = 1 for a helical configuration.

The remaining static magnetic field models are used to describe quadrupole magnets, using

$$\mathbf{B}\_{\rm Q}(z) = B\_{\rm Q0}(z) \left( y \mathbf{\hat{e}}\_x + x \mathbf{\hat{e}}\_y \right), \tag{24}$$

where *BQ*0(*z*) is the constant field gradient over some range *z*1 ≤ *z* ≤ *z*2, and dipole magnets using a field model that is described by a constant field oriented perpendicularly to the axis of symmetry over some range *z*3 ≤ *z* ≤ *z*4. Both, the field models for quadrupole and dipole magnets, have hard-edge field transitions and are curl- and divergence-free over the range where these are defined.

#### *2.3. Performance of the MINERVA FEL Code*

To illustrate the performance of the MINERVA FEL code, we briefly compare the predictions of this code with experimental data from two SASE FEL experiments. We consider the "Sorgente Pulsata ed Amplificata di Radiazione Coerente" (SPARC) experiment, a SASE FEL located at ENEA Frascati [197], and the Linac Coherent Light Source (LCLS), which is a SASE FEL at the Stanford Linear Accelerator Center [198].

#### 2.3.1. The SPARC SASE FEL

The experimental parameters of SPARC [197] are as follows. The electron beam energy was 151.9 MeV, with a bunch charge of 450 pC, and a bunch width of 12.67 ps. The peak current was approximately 53 A for a parabolic temporal bunch profile. The *x*and *y*-emittances were 2.5 mm-mrad and 2.9 mm-mrad respectively, and the rms energy spread was 0.02 percent. There were six undulators, each of which were 77 periods in length, with a period of 2.8 cm and an amplitude of 7.88 kG. Each undulator was modeled using Equation (17) with one period for the entrance up-taper and another for the exit down-taper. Furthermore, in the simulation, eight undulators were used to show saturation of the system. The gap between the undulators was 0.4 m in length and the quadrupoles (0.053 m in length with a field gradient of 0.9 kG/cm), forming a strong focusing lattice were located 0.105 m downstream from the exit of the previous undulator. Note that the quadrupole orientations were fixed and did not alternate. The electron beam was matched into the undulator/focusing lattice. The resonance occurred at a wavelength of 491.5 nm. In the experiment, the pulse energies were measured in the gaps between the undulator segments.

A comparison of the pulse energy as found in the simulation (blue line) and from the experiment (red markers) is shown in Figure 2, where the simulation result is averaged over 20 simulation runs with different noise seeds. This yields convergence to better than 5 percent. Energy conservation in the simulation is maintained to within better than one part in 104. Each marker represents a single measurement that is repeated several times, while the error bar indicates the error in the measurement. The experimental data are courtesy of L. Giannessi. Agreement between the simulation and the measured performance is excellent.

A comparison between the evolution of the relative linewidth as determined from simulation (blue line) and by measurement (red markers, data courtesy of L. Giannessi) is shown in Figure 3. Agreement between the simulation and the measured linewidth is within about 35 percent after 15 m. As shown, the predicted linewidths are in substantial agreemen<sup>t</sup> with the measurements.

**Figure 2.** Simulated (blue line) and measured (red crosses) pulse energy as a function of the propagation distance *z* for the SPARC experiment. The different points correspond to individual measurements and the error bar indicates the error in each of the measurements. The simulation result is an average over 20 noise realizations.

**Figure 3.** Simulated (blue line) and measured (red data points) relative linewidth as a function of the propagation distance *z* for the SPARC experiment. The error bar indicates the error in the measurement. The simulation result is, on average, over 20 noise realizations.

#### 2.3.2. The LCLS SASE FEL

To illustrate the performance of MINERVA at much shorter wavelengths, we also briefly compare with experimental data from LCLS [198], which is a SASE FEL user facility that became operational in 2009, operating at a 1.5 Å wavelength.

To operate at 1.5 Å, LCLS employs a 13.64 GeV/250 pC electron beam with a flat-top temporal pulse shape of 83 fs duration. The normalized emittance (*x* and *y*) is 0.4 mm-mrad and the rms energy spread is 0.01 percent. The undulator line consisted of 33 segments with a period of 3.0 cm and a length of 113 periods. In the simulation, each segmen<sup>t</sup> is modeled using Equation (17), with one period each in entry and exit tapers. A mild down-taper in field amplitude of −0.0016 kG/segment starting with the first segmen<sup>t</sup> (with an amplitude of 12.4947 kG and *K*rms = 2.4748) and continuing from segmen<sup>t</sup> to segmen<sup>t</sup> was used. This is the so-called gain taper. The electron beam was matched into a FODO lattice consisting of 32 quadrupoles, each with a field gradient of 4.054 kG/cm and a length of 7.4 cm. Each quadrupole was placed a distance of 3.96 cm downstream from the end of the preceding undulator segment, and the gap lengths between the undulators followed a repetitive sequence of short (0.48 m)-short (0.48 m)-long (0.908 m).

The LCLS produces pulses of about 1.89 mJ at the end of the undulator line [198], and saturation is found after about 65–75 m along the undulator line. A comparison between the measured pulse energies (red circles), as obtained by giving the electrons

a kick to disrupt the FEL process, and the simulation (blue) is shown in Figure 4. The experimental data are courtesy of P. Emma and H.-D. Nuhn at SLAC, and the simulation results represent an average over an ensemble of 25 runs performed with different noise seeds. As shown in Figure 4, the simulations are in good agreemen<sup>t</sup> with the measurements in the exponential growth region with close agreemen<sup>t</sup> for the gain length. The simulation exhibits saturation at the same distance as the experiment in the range of 65–75 m at a pulse energy of 1.5 mJ. After saturation, in view of the gain taper, the pulse energy grows more slowly to about 2.02 mJ at the end of the undulator line, which is approximately 8 percent higher than the observed pulse energy.

**Figure 4.** Simulated (blue line) and measured (red circles) optical pulse energy as a function of the propagation distance *z* for the LCLS experiment. The simulation result is on average over 25 noise realizations.

The agreemen<sup>t</sup> between simulation and experiment for the pulse energy is poorer during the early stages of the interaction. This may be due to a variety of reasons. On the experimental side, as the pulse energy grows by 5 to 6 orders of magnitude from the initial shot noise to saturation, it is difficult to calibrate the detectors for the low pulse energies at the early stages of the interaction. Moreover, while kicking the electrons provides a fast measurement method, it is accompanied by a larger background signal and the possibility of restarting the FEL process downstream the undulator line when the kick is performed at the beginning of the undulator line [198]. On the simulation side, there may be some inaccuracies in the shot noise algorithm or the finite number of optical modes used that underestimates the initial noise level.

#### *2.4. The Optical Propagation Code OPC*

When modeling FEL oscillators, the electron–light interaction within the undulator segments needs to be modeled, as well as the propagation of the light through the remainder of the oscillator. The light propagators internal to FEL gain codes, including that in MINERVA, are typically not very efficient when considering free-space propagation over large distances and interaction with various optical elements.

The wavefront propagator [185,199] is based upon the scalar paraxial Helmholtz equation [200]

$$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} - 2ik\_0 \frac{\partial u}{\partial z} = 0,\tag{25}$$

where *u* is the complex scalar wave amplitude that describes the transverse profile of a time-harmonic optical beam, *k*0 = 2*π*/*λ*0, and *λ*0 is the free-space wavelength. By applying a two-dimensional spatial Fourier transform, Equation (25) transforms into

$$(k\_x^2 + k\_y^2)\ddot{\imath} - 2ik\_0 \frac{\partial \dddot{\imath}}{\partial z} = 0,\tag{26}$$

where

$$\overline{u}(k\_x, k\_y, z) \equiv \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} dx dy \, u(x, y, z) e^{-i(k\_x x + k\_y y)}\tag{27}$$

is the two-dimensional Fourier transform of *<sup>u</sup>*(*<sup>x</sup>*, *y*, *<sup>z</sup>*). Given a known profile *<sup>u</sup>*0(*<sup>x</sup>*, *y*) = *<sup>u</sup>*(*<sup>x</sup>*, *y*, 0) at *z* = 0, e.g., at the exit of the undulator line, the analytical solution in reciprocal space to Equation (26) is given by

$$\vec{u}(k\_x, k\_y, z) = \vec{u}\_0(k\_x, k\_y) e^{iz\lambda \rho \left(k\_x^2 + k\_y^2\right)/\pi}.\tag{28}$$

The solution in normal space at location *z* is found by the inverse spatial Fourier transform of Equation (28),

$$u(x,y,z) \equiv \frac{1}{4\pi^2} \int \int \int dk\_x dk\_y \,\vec{u}(k\_{x\prime}k\_{y\prime}z)e^{i(k\_x x + k\_y y)} \,\mathrm{.}\tag{29}$$

This wavefront propagation is known as spectral propagation and is very suitable for propagation over short and large distances. However, when propagating over long distances, care has to be taken that the complete optical field remains contained within the transverse domain of the two-dimensional spatial Fourier transform (Equation (27)), e.g., as a result of a fast Fourier transform implementation, to avoid artificial reflections from the edges of this domain.

By substituting Equation (28) into Equation (29), using the initial profile *<sup>u</sup>*0(*<sup>x</sup>*, *y*) and Equation (27), changing the order of integration and performing the integration over *kx* and *ky*, we obtain the Fresnel diffraction integral [201]

$$u(x,y,z) = -\frac{i}{\lambda\_0 z} \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} d\xi d\eta \, u\_0(\xi, \eta) e^{i\pi[(x-\xi)^2 + (y-\eta)^2]/\lambda\_0 z}.\tag{30}$$

Note that the Fresnel diffraction integral can be derived using various methods, including the free-space Green's function [202], or by applying the paraxial approximation to the spherical wavelets in Huygens' integral [200]. It can also be shown that propagation through a cascaded set of paraxial optical components described by an overall ray-optical ABCD matrix can be described by a single Huygens' integral, given by [200]

$$u(x,y,z) = e^{-ik\_0z} \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} d\xi d\eta \, \mathcal{K}(x,y,\xi,\eta) u\_0(\xi,\eta),\tag{31}$$

where Huygens' kernel *K* is given by

$$K(\mathbf{x}, \mathbf{y}, \xi, \eta) \equiv \frac{i}{\lambda\_0 \sqrt{B\_\mathbf{x} b\_\mathbf{y}}} e^{-i \frac{\pi}{B\_\mathbf{x} \lambda\_0} (A\_\mathbf{x} \xi^2 - 2\xi \mathbf{x} + D\_\mathbf{x} \mathbf{x}^2)} e^{-i \frac{\pi}{B\_\mathbf{y} \lambda\_0} (A\_\mathbf{y} \eta^2 - 2\eta \eta + D\_\mathbf{y} \eta^2)},\tag{32}$$

and *Ai*, *Bi*, *Ci* and *Di* are the components of the system ABCD matrix for *i* = *x*, *y*. The ABCD matrix can describe real or complex orthogonal paraxial optical systems that may contain astigmatism [200]. Note, in the presence of apertures, the optical field must be propagated from aperture to aperture, since apertures cannot be included in the ABCD matrix.

An efficient calculation of the integrals appearing in the spectral propagation method and Fresnel's diffraction integrals relies on fast Fourier transforms. By applying a transformation to Equation (31), it can be put into a form that also allows the use of fast Fourier transforms to efficiently evaluate the integrals. We define *Mx* = *<sup>a</sup>*2/*<sup>a</sup>*1, where *a*1 and *a*2 are measures for the size of the optical field, such that the optical field becomes negligible for |*x*| ≥ *<sup>a</sup>*1,2 at the input and output plane, respectively. Similarly, *My* = *<sup>a</sup>*3/*<sup>a</sup>*4 where *a*3 and *a*4 are measures for the size of the optical field, such that the optical field becomes negligible for |*y*| ≥ *<sup>a</sup>*3,4 in the input and output plane, respectively. Using the transformation *x* = *a*1*x*, *ξ* = *<sup>a</sup>*2*ξ*, *y* = *<sup>a</sup>*3*y*, *η* = *<sup>a</sup>*4*η*,

$$w\_0(\boldsymbol{\xi}^{\prime}, \boldsymbol{\eta}^{\prime}) = \sqrt{a\_1 a\_3} \mu\_0(\boldsymbol{\xi}^{\prime}, \boldsymbol{\eta}) e^{-i \pi \frac{(A\_{\overline{x}} - M\_{\overline{x}}) \hat{\boldsymbol{\xi}}^2}{B\_{\overline{x}} \lambda\_0}} e^{-i \pi \frac{(A\_{\overline{y}} - M\_{\overline{y}}) \boldsymbol{\eta}^2}{B\_{\overline{y}} \lambda\_0}},\tag{33}$$

and

$$w(\mathbf{x}', \mathbf{y}') = \sqrt{a\_2 a\_4} \mu(\mathbf{x}, \mathbf{y}) e^{i\pi \frac{(D\_\mathbf{x} - M\_\mathbf{x}^{-1})\mathbf{x}^2}{B\_\mathbf{x}\lambda\_0}} e^{i\pi \frac{(D\_\mathbf{y} - M\_\mathbf{y}^{-1})\mathbf{y}^2}{B\_\mathbf{y}\lambda\_0}}.\tag{34}$$

Equation (31) can be transformed into

$$w(\mathbf{x}', \mathbf{y}') = i \sqrt{\mathcal{N}\_{\mathbf{c}, \mathbf{x}} \mathcal{N}\_{\mathbf{c}, \mathbf{y}}} \int \prod\_{-1}^{1} d\xi' d\eta' \, \mathcal{K}(\mathbf{x}', \mathbf{y}', \xi', \eta') \text{vol}(\xi', \eta') \tag{35}$$

with the kernel *K* given by

$$K(\mathbf{x'}, \mathbf{y'}, \mathbf{\tilde{y'}}, \boldsymbol{\eta'}) = e^{i\pi N\_{\mathbf{c}, \mathbf{x}} (\mathbf{x'} - \mathbf{y'})^2} e^{i\pi N\_{\mathbf{c}, \mathbf{y}} (\mathbf{y'} - \boldsymbol{\eta'})^2} \,. \tag{36}$$

In Equations (35), *Nc*,*<sup>x</sup>* and *Nc*,*<sup>y</sup>* are equivalent collimated Fresnel numbers [200] defined as

$$N\_{\mathfrak{c},\mathfrak{x}} = \frac{M\_{\mathfrak{x}}a\_1^2}{B\_{\mathfrak{x}}\lambda\_0}, \; N\_{\mathfrak{c},\mathfrak{y}} = \frac{M\_{\mathfrak{y}}a\_3^2}{B\_{\mathfrak{y}}\lambda\_0}.\tag{37}$$

We shall refer to Equation (35) as the modified Fresnel diffraction integral, and this integral can also be efficiently evaluated using fast Fourier transforms, with the added benefit that independent magnification factors in the *x*- and *y*-directions can be applied in going from the input to the output plane. This means that the output plane can grow or shrink in size with the optical beam when it expands or contracts in propagating from the input to the output plane that are defined by the system ABCD matrix used in the propagation.

#### 2.4.1. Optical Elements

To allow the modeling of FEL oscillators, the optical field needs to interact with various types of optical elements, such as lenses, mirrors or diaphragms. Diaphragms are implemented as hard edge apertures that reduce the optical field to zero outside the aperture. Various apodization functions are available that can be applied to the edge of the aperture. Lenses and mirrors are implemented as elements that apply a phase shift to the optical field. This is done by multiplying the optical field by *<sup>e</sup>*<sup>−</sup>*iq*(*<sup>x</sup>*,*y*), where *q*(*<sup>x</sup>*, *y*) is the local phase shift in the transverse plane. For example, a thin lens is modeled as *q*(*<sup>x</sup>*, *y*) = *<sup>k</sup>*0(*x*<sup>2</sup> + *y*<sup>2</sup>)/2 *f* with *f* as the focal strength of the lens. Mirrors are modeled as thin lenses using *f* = *R*/2, where *R* is the radius of curvature of the mirror. More complicated optical elements as thick lenses or a combination of lenses can be implemented by determining the overall ABCD matrix and using the modified Fresnel diffraction integral (Equation 35). OPC also allows for more complicated phase masks that can be created using Zernike polynomials [194]. These can be used to implement not only various types of mirror and lens aberrations, but also mirror distortion resulting from thermal loading [203,204] (see Section 2.4.2 for more detail). In the latter case, OPC can set a scaling factor for the aberrations depending on the (average) optical power loading of the mirror.

Typically, it is assumed that the FEL gain bandwidth is sufficiently small that dispersive effects within the optical elements can be neglected. However, this is not the case when crystal Bragg mirrors are used, e.g., to model oscillator configurations at X-ray wavelengths. For such mirrors, the angle of reflection and mirror loss strongly depend on the type of crystal used, its orientation and the X-ray photon energy [205]. To properly model the reflection of Bragg mirrors, the incident optical field needs to be Fourier transformed into the frequency domain to handle the different photon energies in the optical field, and a

two-dimensional spatial Fourier transform of the incident field is needed to deal with the various angles of incidence on the mirror.

Several optical elements can be combined to form a more complex optical component, e.g., by combining a mirror with a diaphragm element, extraction of radiation from a resonator through a hole in one of the mirrors can be modeled. Another example is to use an external finite element program to simulate mirror surface distortion due to thermal loading of the mirror. This surface distortion can then be converted to a phase mask, the amplitude of which can be dynamically scaled at each roundtrip to model mirror distortion during the start-up of the oscillator.

Finally, the collection of wavefront propagators together with the optical components and interfaces with several FEL gain codes is known as the optical propagation code (OPC). As OPC is specifically designed to work several FEL gain codes to model FEL oscillators, some of the optical components allow for the forking of the optical propagation path. For example, at the mirror used for coupling the light out of the resonator, the tracking of the light within the resonator can temporarily be suspended to propagate the light outside the resonator to some diagnostic point where several optical properties of the light can be obtained. Afterwards, the propagation of the light within the resonator can be continued. This allows simultaneous monitoring of the light properties at some external diagnostic point, as well as the build-up of the light within the resonator.

#### 2.4.2. Modeling Mirrors

The default optical components used by OPC are ideal thin lenses and spherical mirrors. However, by using amplitude and phase masks, additional optical components can be added, in particular more realistic lenses and mirrors. Here, we limit ourselves to phase masks that are generated using the circle polynomials *R*|*m*<sup>|</sup> *n* of Zernike [194]. These polynomials are used to generate a phase difference *dθ* defined on a transverse plane that is applied to the optical field as

$$d\theta = D\_{mm} R\_n^{|m|}(\rho) \times \begin{cases} \cos \left( m\phi \right) & m \le 0 \\ \sin \left( m\phi \right) & m < 0 \end{cases} \tag{38}$$

where *ρ* is a scaled radial distance, *ρ* = *x*<sup>2</sup> + *<sup>y</sup>*2/*ρc* with *ρc* some characteristic length scale, *φ* is the angle tan <sup>−</sup><sup>1</sup>(*y*/*x*) and *Dmn* is the amplitude of the polynomial. The indices *m*, *n* describe the type of aberration, for example, *m* = 0 and *n* = 4 correspond to spherical aberration and *m* = 1 and *n* = 3 to coma [194]. Therefore, adding a phase mask and applying it at the location of a lens of mirror adds the associated aberrations to create a more realistic model of the component.

It is also possible to use these Zernike polynomials to generate new optical components. For example, adding the polynomials *m*, *n* = 2, 2 with *m*, *n* = 0, 2 with appropriate amplitudes generates a cylindrical lens of a certain focal strength. This can be used to correctly model the performance of a spherical mirror with a non-zero angle of incidence, e.g., when used in a ring resonator, to model for the different foci in tangential and sagittal direction.

Another application is to use the Zernike polynomials to model mirror surface distortions, in particular, due to thermal loading of the mirror. Here, we assume that round mirror is uniformly cooled at its circumference and that the heat loading is axi-symmetric. A finite element program can be used to calculate the steady state deviation *<sup>δ</sup>z*(*r*/*<sup>r</sup>*0), where *r*0 is some characteristic length, that results from the heat generation in the mirror due to some absorbed power *P*abs. A fit using *R*<sup>0</sup> *n*(*r*/*<sup>r</sup>*0) Zernike polynomials with *n* taken even is then used to determine the amplitudes to be used in generating the phase mask. To convert the mirror distortion *δz* into a localized change in phase for the optical field, we set

$$d\theta(r/r\_0) = -\frac{4\pi\delta z(r/r\_0)}{\lambda},\tag{39}$$

where the minus sign is required to comply with the phase advance used in OPC and the phase change is proportional to twice the mirror displacement. It is reasonable to assume that a stable transverse field distribution is established after starting from noise when the intracavity optical starts to heat the mirror. Therefore, the phase distortion due to mirror heating can be obtained by scaling Equation (39) with a factor *Pi*(*t*)*τ*/(*<sup>P</sup>*abs*T*rep) with *P*abs the absorbed power used to calculate *δ<sup>z</sup>*, *Pi* the instantaneous absorbed power, *τ* the duration of the optical pulse and *<sup>T</sup>*rep the repetition time for the optical pulses.

#### **3. Some Considerations for Optical Resonators Used in Free-Electron Lasers**

Most optical resonators used in laser oscillators consist of two spherical mirrors separated by a distance *L*cav, as schematically shown in Figure 1. When the radius of curvature *Rc* of the mirrors is equal, the resonator is called concentric when *L*cav = 2*Rc* and confocal when *L*cav = *Rc*. Other configurations with an unequal radii of curvature of different distances between the mirrors are known as generalized resonators.

Most of the low-gain FEL oscillators employ stable concentric resonators to ensure that the optical pulse train does not "walk" out of the cavity. Consider a generalized resonator with a cavity length *L*cav, as shown in Figure 5, which shows the radii of curvature, *Rc*,*<sup>i</sup>*(*<sup>i</sup>* = 1, <sup>2</sup>), for the two mirrors, as well as the distances *zi* from the mode waist to each of the mirrors. The transverse optical modes of the cold resonator, i.e., without a gain medium (electron beam) present, can be described as either Gauss–Hermite *cf.* Equations (1)–(3) or Gauss–Laguerre modes.

**Figure 5.** Schematic representation of a generalized two-mirror resonator, defining some of its parameters.

Defining

$$g\_i = 1 - \frac{L\_{\text{cav}}}{R\_{c,i}} \tag{40}$$

then the Rayleigh range *zR* for the Gaussian eigenmodes in the cavity is given by

$$z\_R = L\_{\rm cav} \sqrt{\frac{g\_1 g\_2 (1 - g\_1 g\_2)}{g\_1 + g\_2 - 2g\_1 g\_2}} \tag{41}$$

Since the Rayleigh range must be positive, we must have 0 ≤ *g*1*g*2 ≤ 1, which is the condition for resonator stability [200]. This means that Gaussian beams can only be eigenmodes of stable resonators. The stability condition implies that

$$L\_{\rm cav} < R\_1 + R\_2. \tag{42}$$

We also note that the waist of the Gaussian mode is at a distance

$$z\_1 = L\_{\rm cav} \frac{g\_2(1 - g\_1)}{g\_1 + g\_2 - 2g\_1 g\_2} \tag{43}$$

from mirror 1 and at a distance *z*2 = *L*cav − *z*1 from mirror 2. Note that the sign convention used for the ray-optic ABCD formalism [200] applies here as well, meaning that *Rc* > 0 for a concave mirror and *Rc* < 0 for a convex mirror and that the distances *zi* have a sign.

The Gauss–Hermite transverse modes as described by Equations (1)–(3) are eigenmodes of an empty generalized resonator with a particular Rayleigh range and location of the waist defined by the resonator. Whenever the resonator contains internal apertures clipping the optical field, partial waveguiding, or if the light is extracted via a hole, a single transverse Gaussian mode can no longer be an eigenmode of the resonator. The Fox–Li method [206] can be used to find the eigenmodes of such resonators [207].

The resonant frequencies, or longitudinal modes, of the resonator are those frequencies for which the total roundtrip phase is an integer multiple of 2*π*. As a result, the resonant frequencies differ for the different transverse optical modes [200]. The frequency spacing between two subsequent resonances of the same transverse mode is known as the freespectral range Δ*ν*FSR. FEL resonators typically have a Δ*ν*FSR that is much smaller than the gain bandwidth, due to the large mirror separation used. Consequently, most FEL oscillators are operating on multiple longitudinal modes.

As shown in Figure 1, the resonator may contain more than one optical pulse if the roundtrip time in the resonator is larger than the time between subsequent electron bunches. The resonator design must ensure synchronism between the optical pulses and the electron bunches and optimal performance is typically found near the synchronous or zero-detuning cavity length, where the incoming electron bunch coincides with the returning optical pulse train. This zero-detuning cavity length, *L*0, is given by

$$\frac{\mathcal{M}}{f\_{\text{rep}}} = \frac{2L\_0 - L\_u}{c} + \frac{L\_u}{v\_{\mathcal{S}}} \, \, \, \tag{44}$$

where *M* is the number of simultaneous electron bunches in the optical cavity, *f*rep is the repetition rate of the electron bunches, *Lu* is the undulator length and *vg* is the group velocity of the light within the undulator. Solving for *L*0 gives

$$L\_0 = \frac{cM}{2f\_{\rm rep}} - \frac{L\_u}{2} \frac{c}{v\_\%} \left(1 - \frac{v\_\%}{c}\right),\tag{45}$$

which reduces to the well-known expression

$$L\_0 = \frac{c\mathcal{M}}{2f\_{\rm rep}},\tag{46}$$

for low gain oscillators where *vg* ≈ *c*. As will be discussed in the Section 5 on highgain/low-*Q* oscillators, exponential gain in the undulator results in a reduction in the group velocity, which has the effect of shortening the zero-detuning cavity length compared to that for low-gain/high-*Q* oscillators (Equation (46)).

In general, the dynamics in an oscillator are set by the interplay between the instantaneous gain and loss, where in a free-electron laser, the gain is provided by the electrons streaming through the undulator and loss is due to out-coupling of the radiation and other loss mechanisms that may be present within the resonator. The total loss of the resonator determines the quality factor *Q* of the resonator [200]. To obtain laser oscillation, the roundtrip small-signal gain has to be higher than the loss per roundtrip. The so-called threshold gain is the roundtrip small-signal gain needed to just balance the loss per roundtrip. When the light builds up inside the resonator, the gain is reduced and a stationary state is obtained when the saturated gain equals the loss per roundtrip. Note that

this is the same condition as for threshold, i.e., the saturated gain is equal to the threshold gain. If we denote the output power on the *n*th pass as *Pn*, then the power after the (*n*+1)th pass is *Pn*+<sup>1</sup> = (1 − *L*)(1 + *<sup>G</sup>*)*Pn*. Equilibrium is characterized by *Pn*+<sup>1</sup> = *Pn*; at which point, *G* = *L*/(1 − *<sup>L</sup>*).

In FELs, where the bunch charge/peak current is relatively small and the gain cannot reach the exponential regime in a single pass through the undulator, the resonator *Q* must be high enough that the oscillator can lase. In such cases, most of the energy extracted from the electrons remains circulating in the resonator and a small fraction of the power is outcoupled. This implies that the loss, which is often dominated by the fraction of radiation that is coupled out of the resonator, but also includes losses at the mirrors or within the resonator, can impose difficulties if mirror absorption is large enough to result in excessive heating or degradation of the reflectivity. In such cases, the mirrors must be cooled or otherwise protected. For example, storage ring FELs have a low small-signal gain and thus require high-quality mirrors that are susceptible to the higher-harmonics generated. Short-wavelength X-ray FEL oscillators also have a low small-signal gain and the available X-ray mirror technology is pushed to its limits to produce mirrors with sufficiently low loss. High average power operation may require the cryogenic cooling of the mirrors to prevent mirror distortion due to thermal loading.

The efficiency of a low-gain/high-*Q* oscillator is predicted to be *η* = 1/2.4 *Nu* [1], where *Nu* denotes the number of uniform periods in the undulator; hence, high efficiency requires relatively short undulators. Because of this, oscillator design is a balance between having a sufficiently long undulator for the gain to exceed the losses, set by the quality *Q* of the resonator, but not so long that the efficiency is negatively impacted.

#### **4. Low-Gain/High-***Q* **Oscillators**

In this section, we describe simulations of two low-gain/high-*Q* infrared FEL oscillator experiments conducted at JLab: the IR-Demo [50,64,65] and the 10-kW Upgrade [11]. Both of these experiments were conceived as demonstrations for high average power infrared FELs based on energy recovery linacs.

#### *4.1. The IR-Demo Experiment at JLab*

The IR Demo experiment [50,64,65] was based on a superconducting energy recovery linac at JLab that produced 0.4 ( ±0.1) ps rms electron bunches, with energies of about 38 MeV bunch charges of 60 pC (60 A peak current) at a repetition rate of 18.7 MHz corresponding to an average current of 1.2 mA. The transverse emittance of the bunches was 7.5(±1.5) mm-mrad (rms) and the rms energy spread was about 0.25 percent. The experiment uses a 42-period flat-pole-face planar undulator with a period of 2.7 cm and an on-axis amplitude of 5.56 kG, yielding an undulator strength parameter of *K*rms = 0.99. Therefore, the simulation model uses Equation (17) to model the undulator with one period entry and exit tapers. The resonant wavelength was 4.8 μm and a concentric resonator was used with a optical cavity length of about 8 m [208]. The radii of curvature of the mirrors were 4.045 m with a cold cavity Rayleigh range of 40 cm. Transmissive out-coupling through the downstream mirror was used and had a transmittance of approximately 10 percent. The outer radius of the mirrors was 2.54 cm. OPC uses the modified Fresnel propagator (Equation (35)) to handle the divergence and convergence of the optical beam inside the resonator using a fixed number of grid points.

Simulations of this experiment were conducted with MEDUSA/OPC. MEDUSA [184] and MINERVA both employ Gaussian representations for the optical field and integrate the three-dimensional Lorentz force equations for the electron trajectories without performing an average over the wiggle-motion. However, MEDUSA is the more primitive code and MINERVA contains many additions and incorporates superior algorithms not present in MEDUSA. Nevertheless, in cases where the two codes have been compared, their predictions for FEL performance are in agreemen<sup>t</sup> to within about 10 percent.

The evolution of the output energy found in simulation versus pass is shown in Figure 6 for a cavity length near the peak of the detuning curve of 8.01049 m. Observe that saturation is found after about 70 passes at a pulse energy of 0.17 mJ, corresponding to an average power of about 318 W, at a repetition rate of 18.7 MHz. This represents an average efficiency of about 0.70 percent, which is close to the theoretically predicted efficiency of about 1.0 percent for an undulator with 40 uniform periods.

**Figure 6.** Simulated pulse energy as a function of roundtrip number *n* at a near optimal cavity detuning of *L*cav = 8.01049 m for the IR Demo experiment. Remaining parameters as described in the text.

The average power found in simulation near the peak in the detuning curve is consistent with the observation of 311 W during CW operation of the IR Demo [50,64,65]. This is also shown in Figure 7, which shows the detuning curves found in simulation (blue) and during pulsed operation of the IR Demo (green). The average power during CW operation is indicated by the dashed line. Observe that the zero-detuning point in the simulation is shifted by about 5 μm from that observed in the experiment. This may be due to an uncertainty in the cavity length measurement by up to 10 μm. The full width of the detuning curve is found to be within about 30–35 μm.

**Figure 7.** Tuning curve for the IR Demo experiment including measurements from CW and pulsed runs (green) and from simulation (blue). Other parameters as described in the text.

Higher bunch charges were used in the pulsed runs than in the CW runs in the IR Demo. Lower bunch charges were used in the CW runs in order to reduce distortion due to mirror heating. Indeed, mirror distortion made it difficult to obtain a detuning curve during CW operation. Mirror heating was not a serious problem in the pulsed runs. When the gain to loss ratio is small, as in the CW runs and in the simulation, it is expected that the tuning curve for a low-gain/high-Q oscillator will display a peak at the zero detuning point, followed by a rapid decline as the cavity length decreases, and this is what is found in simulation. However, when the gain-to-loss ratio increases, then the detuning curve will exhibit a shoulder, as indicated in the detuning curves for the pulsed runs. This shoulder will be much more prominent in the RAFEL simulations discussed below.

#### *4.2. The 10-kW Upgrade Experiment at JLab*

This experiment represented an upgrade to the original IR-Demo experiment [65] and numerous elements represent upgrades. For example, the accelerating modules were upgraded to achieve higher energy, bunch charge, and average current. Effort was also made to reduce/mitigate the beam breakup instability in view of the higher average current. The original undulator was replaced to achieve high gain and resonance at a shorter wavelength. In order to handle the higher power, the resonator was lengthened to reduce the mirror loading and cryogenic, edge-cooling was used for the mirrors. As a result, the kinetic energy was increased to 115 MeV, while the energy spread of 0.3 percent remained approximately the same. The charge of the electron bunch was almost doubled to 115 pC with a pulse length of 390 fs. The normalized emittance of 9 mm-mrad in the wiggle plane and 7 mm-mrad in the plane orthogonal to the wiggle plane was similar as in the IR demo experiment, and the maximum repetition rate of 74.85 MHz for the electron bunches remained unchanged. Note, that although the IR-Demo experiment was capable of running at 74.85 MHz, the results reported in Section 4.1 corresponded to operation at 17.85 MHz. To operate at somewhat shorter wavelengths, the planar undulator, which is modeled in the simulation using Equation (19) with one period up- and down-taper, had a longer period of 5.5 cm, a total of 30 periods, and a peak on-axis magnetic field of 3.75 kG. The electron beam was focused into the undulator with the focus at the center of the undulator. The concentric resonator was also updated [209] and had a length of about 32 m with a cold-cavity Rayleigh length of 0.75 m. The total loss in the resonator was 21 percent with about 18 percent out-coupled per pass from the downstream mirror. For these settings, the wavelength was 1.6 μm.

In simulating this experiment, the number of particles in MINERVA was 5832 per slice, while the separation between slices was 5.4 fs. The number of optical modes was dynamically adjusted each roundtrip to accommodate the evolution of the optical field inside the resonator. OPC uses, again, the modified Fresnel propagator (Equation (35)) to handle the divergence and convergence of the optical beam inside the resonator.

The length of the optical cavity must be selected so that the returning optical pulse is in synchronism with the electron bunches. The roundtrip time for the optical pulses in the cavity is *τr* = 2*L*cav/*c* and the separation between electron bunches is *<sup>τ</sup>*sep = 1/ *f*rep, where *L*cav is the cavity length and *f*rep is the electron bunch repetition rate. Perfect synchronism (referred to as zero-detuning) is obtained when *τr* = *<sup>M</sup>τ*sep, which leads to Equation (46). Here, *M* is the number of optical pulses in the cavity. In this case, there were 16 optical pulses in the cavity and the zero-detuning length is *L*0 = 32.041946079 m. The cavity detuning curve obtained from simulations is shown in Figure 8 as a function of the difference between the cavity length *L*cav and the zero-detuning length. With a maximum pulse energy of 0.194 mJ and a repetition rate of 74.85 MHz, we find that the maximum output power of 14.52 kW occurs for a positive detuning of 2 μm and is close to the measured value of 14.3 ± 0.72 kW [69]. As a result, the predicted extraction efficiency is about 1.4 percent, which is close to the theoretical value of *η* = 1.7 percent. We remark that the previous simulation of this experiment with MEDUSA/OPC [183] yielded an average output power of 12.3 kW, and the present formulation is in better agreemen<sup>t</sup> with the experiment than in the earlier simulation. As in the previous simulation [183], the roughly triangular shape of the detuning curve is also in agreemen<sup>t</sup> with the experimental observation.

**Figure 8.** Tuning curve for the IR upgrade experiment. Simulation parameters as described in the text.

The temporal profiles of the optical pulse at the undulator entrance and exit as well as that of the electron bunch current are shown in Figure 9 for the zero-detuning cavity length after pass 100, which corresponds to a stable, saturated steady-state. Observe that the electron bunch is centered in the time window, which has a duration of 1.4 ps. This is at zero-detuning, as indicated by the fact that the incoming optical pulse at the undulator entrance is in close synchronism with the electron bunch. It is also evident that the center of the optical pulse advances by about 0.16 ps, as it propagates through the undulator, and this is in good agreemen<sup>t</sup> with the theoretical slippage estimate of *Nuλ*/*c* for a low-gain FEL, where *Nu* is the number of periods in the undulator. Finally, it should be remarked that this is in the steady-state regime where the losses in the resonator and the out-coupling are compensated by the gain in the undulator.

**Figure 9.** Temporal profiles of the power in the optical pulse after 100 passes at the undulator entrance (green) and exit (blue), as well as the current in the electron bunch (right axis, red) for the IR upgrade experiment. The cavity length is *L*cav = 32.041946079 m and other parameters, as described in the text.

#### **5. High-Gain/Low-***Q* **Oscillators**

An alternate approach to oscillator design is to use a high-gain undulator where the radiation grows exponentially on a single pass through the undulator [2]. Because the gain is high, the permissible resonator loss can be relatively high; hence, a large fraction of the power can be coupled out of the resonator. This type of oscillator has been referred to as a regenerative amplifier FEL (RAFEL), and the concept has been experimentally demonstrated at the Los Alamos National Laboratory [210]. Hence, a RAFEL may also be thought of as a low-*Q* oscillator and has advantages, both for (1) high power designs, since the mirror loading can be kept below mirror deformation or damage thresholds, and (2) for VUV and X-ray oscillators.

RAFELs differ from low-gain oscillators in a number of ways. One difference is that, as shown in Madey's theorem [211], a low-gain oscillator exhibits no gain directly on the resonance. In contrast, the growth rate in the exponential gain regime has a peak onresonance, and this is reflected in the wavelengths excited in a RAFEL. A second difference is the overall efficiency. The saturation efficiency, *η*, of a low-gain oscillator is *η* = 1/2.4*Nu*, where *Nu* is the number of periods in the undulator. Since the radiation exponentiates in each pass through the undulator in the RAFEL, the efficiency is given by that found in the high-gain Compton regime, where *η* = *ρ* where *ρ* is the Pierce parameter. A third difference is in the linewidth, which scales inversely with *Nu* in a low-gain oscillator, but which is given by the linewidth of the exponential interaction in the high-gain Compton regime. A fourth difference is in the longitudinal and transverse mode structure, which is determined largely by the resonator properties in a low-gain oscillator. In a RAFEL, by contrast, the exponential gain leads to saturation in a very small number of passes through the resonator, and the mode structure is largely governed by the interaction in the undulator. A fifth difference is in the effect of slippage. Slippage in a low-gain oscillator scales with *Nu*. However, the high-gain in a RAFEL results in a reduction in the group velocity, such that slippage scales with *Nu*/3. However, one point of similarity that the RAFEL shares with low-gain oscillators, is the presence of limit-cycle oscillations.

In this section, we discuss simulations of an infrared RAFEL with the intention of illustrating many of the general properties of a RAFEL and how it compares both to lowgain/high-*Q* oscillators and SASE FELs [212] and simulations of an X-ray RAFEL concept, making use of hole out-coupling [213]. We note that the infrared RAFEL simulations were performed with MEDUSA/OPC, while the X-ray RAFEL simulations were performed with MINERVA/OPC.

#### *5.1. An Infrared RAFEL*

The electron beam, undulator, and resonator parameters are summarized in Table 1. Observe that the temporal profile of the electron bunch is parabolic with a full width of 1.2 ps, and that the undulator is a two-plane-focusing (i.e., parabolic pole face) design. Consequently, Equation (17) is used to model the undulator with the first and last periods tapered up and down to model the injection and ejection of the beam. Since there is exponential growth and, hence, the optical guiding of the radiation, we use a matched beam in the undulator. The resonator is concentric with the power coupled out through a 5.0 mm hole in the downstream mirror, which provides for a typical average out-coupling of about 97 percent. Given the repetition rate, *f*rep, the nominal zero-detuning cavity length, *L*0, is 6.85239904 m when *M* = 4 and assuming *vg* = *c* (see Equation (46)). The temporal window is an important numerical consideration, and must be chosen to be large enough to accommodate the maximum cavity detuning length that is consistent with pass-to-pass amplification so that the optical pulse remains within the time window for all the usable choices of cavity length. In practice, for this example, we choose a temporal window of 4.0 ps and include 182 temporal slices, which corresponds to the inclusion of one temporal slice every three wavelengths. OPC uses the modified Fresnel propagator (Equation (35)) to handle the divergence and convergence of the optical beam inside the resonator, using a grid with a fixed number of grid points.


**Table 1.** Electron beam, undulator and resonator parameters for the IR RAFEL.

We first consider the single-pass gain because this will affect the performance of the RAFEL. Since the undulator is long enough to achieve exponential growth, we show the gain length, *LG*, found in simulation versus the undulator field strength in Figure 10.The optimal (i.e., minimal) gain length of 0.176 m occurs for an on-axis undulator field strength of about 6.7 kG, which corresponds to an rms undulator strength parameter of *K*rms = 1.06. This is in good agreemen<sup>t</sup> with the prediction based on the parameterization of the interaction by Ming Xie [214].

**Figure 10.** Gain length *LG* versus on-axis undulator field amplitude *Bu* for the IR RAFEL. Other parameters as in Table 1.

It should be noted that the well-known resonance condition *λ* = *<sup>λ</sup>u*(<sup>1</sup> + *<sup>K</sup>*2rms)/2*γ*<sup>2</sup> predicts an undulator field of about 6.81 kG (*K*rms = 1.08) at a 2.2 μm wavelength. This shift in the resonance is due to three-dimensional effects and is in disagreement with the shift in the resonance associated with low-gain oscillators. The gain length has implications over the permissible range of *K*rms for which the RAFEL will operate (i.e., over which there is pass-to-pass amplification). Since the RAFEL will saturate when the gain balances the loss, and the loss for the resonator is about 97 percent, this implies that the RAFEL will operate as long as the single-pass gain exceeds about 3200–3300 percent. In order to identify this range more closely, we perform multi-pass simulations and take the average gain over the first 10 passes. We take an average because there are fluctuations in the gain on a pass-to-pass basis (i.e., limit-cycle like oscillations), which will be discussed in more detail below. The average gain is shown as a function of the on-axis undulator field under the assumption of a cavity length of 8.65238 in Figure 11. This represents a cavity detuning with respect to the zero-detuning length of Δ*L*cav = <sup>−</sup>8*λ*. It is clear from Figure 11 that the gain is relatively constant over the range of about 6.65–6.85 kG (*K*rms = 1.054–1.085) and falls off rapidly as the field diverges outside this range, which is consistent with the behavior of the gain length shown in Figure 10. The cutoff for a gain of about 3300 percent occurs for field levels of about 6.518 kG (*K*rms = 1.03) at the low end and 6.878 kG (*K*rms = 1.09) at the high end, and we do not expect the RAFEL to function outside of this range of undulator fields.

**Figure 11.** Average gain produced by the IR RAFEL over the first 10 passes versus the undulator field *Bu* for a cavity length of *L*cav = 8.62538144 m (Δ*L*cav = <sup>−</sup>8*λ*. Other parameters as in Table 1.

In order to demonstrate how the saturation efficiency of a RAFEL differs from that of a low-gain oscillator, it is instructive to compare the RAFEL with an equivalent SASE FEL. The performance of the RAFEL is shown in Figure 12, where we plot the average pulse energy (blue circles) in the steady-state as a function of the on-axis undulator field for the same cavity detuning (Δ*L*cav = −8*λ*), as used for Figure 11. The error bars characterize the limit-cycle oscillations. The RAFEL reaches its peak pulse energy for an undulator field of 6.678 kG (*K*rms = 1.058) and falls to zero outside the range predicted in Figure 10. We also plot the equivalent SASE saturated pulse energy (red triangles). In order to deal with the statistical fluctuations inherent in the SASE output, a large number of runs were made with different shot noise distributions, and found the average (red triangles) and standard deviations (error bars). Note that the SASE results represent the pulse energies over whatever length of undulator is required to reach saturation. Observe that (1) the RAFEL configuration saturates with a higher pulse energy than the SASE configuration, (2) the fluctuations in the RAFEL in the steady-state regime are comparable in magnitude to the statistical fluctuations found in SASE, however, the fluctuations of the RAFEL are deterministic in nature, and (3) the FWHM of the tuning range in *K*rms is comparable for both the RAFEL and SASE configurations.

The RAFEL saturates with about a 0.28 mJ pulse energy, which corresponds to an extraction efficiency of about 0.64 percent. This compares well with the empirical formula [214] that predicts a saturation efficiency of about 0.76 percent. In contrast, the saturation efficiency of a low-gain oscillator is predicted to be *η* = 1/(2.4*Nu*) = 0.43%.

**Figure 12.** Output pulse energy for the IR RAFEL (blue circles) and an equivalent SASE FEL (red triangles) as a function of the undulator field strength *Bu*. The cavity detuning for the IR RAFEL is Δ*L*cav = <sup>−</sup>8*λ*, while the other simulation parameters are given in Table 1. The error bars indicate energy fluctuations due to limit-cycle-like oscillations for the RAFEL and pulse-to-pulse rms energy fluctuations in case of SASE.

The spectral linewidth of the RAFEL also differs from that of a low-gain oscillator. The full width of the spectrum for a typical low-gain oscillator is given by Δ*ω*/*ω* = 1/*Nu* = 0.01 for the example under consideration. This can be translated into a tuning range over the undulator field, as follows

$$
\left|\frac{\Delta B\_{\rm u}}{B\_{\rm u}}\right| = \frac{1 + K\_{\rm rms}^2}{K\_{\rm rms}^2} \left|\frac{\Delta \omega}{\omega}\right|.\tag{47}
$$

This implies a full width tuning range of Δ*Bu* = 0.063 kG (Δ*K*rms = 0.001), which is narrower than what we find in the simulation. The relative SASE linewidth is given by (Δ*ω*/*ω*)rms = *ρ* [215] where *ρ* denotes the Pierce parameter. Here, *ρ* = 0.0097 and (Δ*ω*/*ω*)rms = 0.0097. Converting this to a tuning range in the undulator field and going from the rms width to a FWHM tuning range, we obtain (Δ*Bu*/*Bu*)FWHM = 0.022, which compares well with the simulation results that give (Δ*Bu*/*Bu*)FWHM = 0.019. Hence, the RAFEL behaves more like a SASE FEL than a typical low-gain oscillator in regards to the spectral linewidth.

Another way in which the RAFEL differs from a low-gain oscillator is in the cavity detuning. In a low-gain FEL oscillator, the zero-detuning length is obtained by assuming that the group velocity *vg* equals the speed of light in vacuo *c* throughout the resonator. However, *vg* is reduced in a RAFEL by the interaction in the undulator, and results in smaller synchronous cavity length. This makes the cavity detuning dependent on the gain of the FEL. The change in group velocity also affects slippage. In a low-gain oscillator, the group velocity reduction is small and the slippage is one wavelength per undulator period; hence, the slippage distance is *<sup>l</sup>*slip = *Nuλ*. However, the slippage per undulator period is reduced in a high-gain RAFEL, or in any FEL where there is exponential growth because they have medium decreases, in both the phase and group velocities. The reduced phase velocity results in the optical guiding of the radiation, while the reduced group velocity results in less slippage. It has been shown that *<sup>l</sup>*slip = *Nuλ*/3 at the resonant wavelength [216]. For the example under consideration, this yields *<sup>l</sup>*slip = 72 μm, which is much less than the slippage length of 220 μm if the RAFEL behaved as a low-gain oscillator.

In order to estimate the effect of this on the detuning length, we note that the zerodetuning length is found by equating the roundtrip time of the radiation through the cavity with the spacing between electron bunches (1/*f*rep), cf. Equation (45). As a result, in the

high gain regime where *vg* = *c*/(<sup>1</sup> + 1/3*γ*<sup>2</sup>), the difference in synchronous cavity length Δ*L*0 for a high-gain RAFEL and a low-gain FEL oscillator is given by Equation (45) as

$$
\Delta L\_0 = \frac{L\_\text{ll}}{6\gamma^2} (1 + K\_{\text{rms}}^2) = -\frac{N\_\text{ll}\lambda}{3}.\tag{48}
$$

This is comparable to what is found in the simulation.

The detuning curve found in the simulation is shown in Figure 13, where we plot the output pulse energy versus cavity detuning for an undulator field of 6.658 kG (*K*rms = 1.055). Here, we define the cavity detuning relative to the nominal zero detuning length (Equation (46) with *M* = 1), so that Δ*L*cav = *L*cav − *L*0. As shown in Figure 13, we find a full width detuning range of about 50–110 μm and an FWHM detuning range of about 40 μm, which are in reasonable agreemen<sup>t</sup> with the estimate based on the one-dimensional analysis of slippage.

**Figure 13.** Cavity detuning curve of the IR RAFEL for an on-axis undulator field of *Bu* = 6.658 kG (*K*rms = 1.055), while the other parameters are given in Table 1.

The temporal evolution of the pulse energy is shown for an undulator field of 6.658 kG (*K*rms = 1.055) in Figure 14, where we plot the pulse energy versus pass number through the undulator for the choice of several cavity detunings that samples the complete detuning curve. It is clear that significant fluctuations are found over a large range of detunings and that both the magnitude and period of the fluctuations decrease as the magnitude of the detuning increases, although the magnitude of the fluctuations decreases as well near the zero-detuning length.

**Figure 14.** Temporal evolution of the pulse energy for cavity detunings of Δ*L*cav = 0, <sup>−</sup>18*λ*, <sup>−</sup>33*λ*, −43*λ* and <sup>−</sup>50*λ*. The on-axis undulator field is *Bu* = 6.658 kG (*K*rms = 1.055), while the other parameters are given in Table 1.

The fluctuations seen in simulation can be rapid and irregular. There are two possible explanations for this. One is that due to the high gain and high out-coupling, small changes in the mode structure from pass to pass can result in relatively large changes in the gain and, hence, the pulse energy. These "small" changes can include variations in the transverse mode structure (both in terms of the modal decomposition and spot size), and the temporal pulse shape. The second explanation, related to the first, is that since we have employed hole out-coupling, these relatively small changes in the transverse mode structure at the mirror can give rise to large differences in the out-coupling of the optical mode. It is not surprising, therefore, that the magnitude of the fluctuations varies depending on the cavity detuning. In Figure 15, we show the variation in the rms magnitude of the fluctuations in the out-coupled pulse energy as a function of the cavity detuning. It is clear from Figure 15 that the fluctuation level is relatively constant, at about the 0.03 mJ level over most of the detuning range, but with rapid declines at the end of the detuning range. Furthermore, the oscillation period is of the order of a few passes through the resonator.

**Figure 15.** Standard deviation of the optical pulse energy as a function of cavity detuning Δ*L*cav/*λ* for a undulator field strength of *Bu* = 6.658 kG ( *K*rms = 1.055), while the other parameters are given in Table 1.

Fluctuations/oscillations have been observed in low-gain oscillators and are referred to as limit-cycle oscillations. The observation of limit-cycle behavior in the low-gain FELIX FEL oscillator corresponds to an oscillation period of [143]

$$
\Delta \tau = -\tau\_{\text{slip}} \frac{L\_{\text{cav}}}{\Delta L\_{\text{cav}}},
\tag{49}
$$

where *<sup>τ</sup>*slip = *<sup>l</sup>*slip/*<sup>c</sup>* is the slippage time. For the case of FELIX, *L*cav = 6 m, *λ* = 40 μm, and *Nu* = 38, and the cavity detuning ranges over about 160 μm. As a result, *<sup>τ</sup>*slip = 5.1 ps and *τr* (=2*L*cav/*c*) = 40 ns is the nominal roundtrip time; hence, this implies that the limit cycle oscillation occurs over a period of about 3 μs or 75 passes for a cavity detuning of −100 μm. Incontrast,ifweapplytheslippagetimeforthehighgainRAFELunderconsideration

 *NuλNuλ*

$$
\Delta \tau = -\frac{N\_\text{u} \lambda}{3 \Delta L\_{\text{cav}}} \frac{L\_{\text{cav}}}{c} = -\frac{\tau\_r}{2} \frac{N\_\text{u} \lambda}{3 \Delta L\_{\text{cav}}}.\tag{50}
$$

As such, we expect the oscillation period to occur on the scale of a small number of passes for the indicated cavity detuning range. This is indeed what is observed in Figure 14. For example, the oscillations occur approximately every 2–4 passes for Δ*L*cav/*λ* = −8, which is consistent with Equation (50). However, there is not a grea<sup>t</sup> deal of variation with detuning possible when the oscillations occur on such a fast time scale, and we must take Equations (49) and (50) as approximate measures of the oscillation period. Still, the observed oscillation period is well described by the formula for the oscillation period for limit-cycle oscillations found in low-gain oscillators when the appropriate slippage is

taken into account. For a high-gain RAFEL, the much lower slippage results in very short oscillation periods.

The limit-cycle-like oscillations in the RAFEL are correlated with the fluctuations/ oscillations in the transverse mode structure. The transverse mode structure in a low-gain oscillator is largely (but not completely) determined by the mode structure in the cold cavity, since the optical guiding of the radiation in the undulator is weak. This is not the case in a RAFEL where the mode is guided through the undulator. As a result, the mode structure that forms as the RAFEL saturates differs substantially from the cold cavity modes, and our choice of a Rayleigh range of 0.5 m serves mainly to determine the radii of curvature of the mirrors. Since the radiation is guided in the undulator, the spot size at the undulator exit may be smaller than it would be in the cold cavity, which means that the Rayleigh range of the radiation as it exits the undulator is smaller than it would be in the cold cavity. This implies, in turn, that the optical mode will expand more rapidly as it propagates to the downstream mirror. Alternatively, decomposing the smaller spot size at the undulator exit in cold cavity modes, necessarily leads to higher order transverse modes in the optical field. After propagating to the outcoupler, the superposition of these modes determine the fraction of the optical field coupled out through the hole, and, similarly, after propagation to the undulator entrance, the superposition sets the field profile at undulator entrance. Small variations in the exponential growth rate, e.g., due to changing coupling of the electrons to the optical field at the undulator entrance, lead to relatively larger effects on the optical guiding of the radiation. This in turn changes the spot size at the undulator exit and, hence, the energy coupled out of the resonator and the spot size at the undulator entrance.

The oscillation in mode size has also been observed in simulation of low-gain oscillators [183] where the magnitude of the oscillation depends on the amount of optical guiding, which may vary between the FEL gain codes used [204]. Furthermore, it is found that mirror aberrations, e.g., spherical aberration, also effect the variation in optical mode size observed from pass to pass [204].

This is illustrated in Figure 16, where we plot the pass-to-pass variation in the width of the optical mode on the downstream and upstream mirrors for *Bu* = 6.585 kG ( *K*rms = 1.043) and Δ*L*cav = <sup>−</sup>8*λ*. It is clear that both the spot size and the fluctuations of the spot size on the upstream mirror are greater than those on the downstream mirror, due to the optical properties of the resonator. At saturation, the location of the smallest optical beam size moves over the axis of the undulator and this changes the optical magnification. Consequently, the size of the optical field at both mirrors as well as at the entrance of the undulator changes from pass to pass, as can be observed in Figure 16.

**Figure 16.** Variation in the width (rms diameter) of the optical beam incident on the upstream mirror (red) and downstream mirror (blue) as a function of the number of passes for a cavity detuning of Δ*L*cav = −8*λ* and an undulator field strength of *Bu* = 6.658 kG ( *K*rms = 1.055). The other parameters are given in Table 1.

The transverse mode structure is not only a result of optical guiding; it is also affected by the hole out-coupling. Consider the case of *Bu* = 6.658 kG ( *K*rms = 1.055) and Δ*L*cav = −18 *λ*. The cross-section of the field delivered to the undulator entrance on pass 60 is shown in Figure 17, where we plot the normalized power in the *x*-direction (i.e., the wiggle plane). Observe that the bulk of the power is at the edge of the optical field, but there is a spike at the center that seeds the subsequent pass through the undulator. Despite the multiple peaks in the cross section at the undulator entrance, the strength of the interaction in the undulator yields a near-Gaussian mode peaked on-axis at the undulator exit, as shown in Figure 18a. What has happened is that the interaction with the electron beam, which has a diameter of about 0.784 mm, essentially amplifies and guides the central peak shown in Figure 17, while the power in the wings falls outside the electron beam and is not amplified. This near-Gaussian mode then propagates to the downstream mirror, during which it expands by about a factor of three, as shown in Figure 18b, where the FWHM is about 3.9 mm in width. The FWHM of the modal superposition at the undulator exit is about 1.2 mm.

**Figure 17.** Cross-section of the optical field at the undulator entrance on pass 60 for a cavity detuning Δ*L*cav = −18*λ* and an undulator field strength of *Bu* = 6.658 kG ( *K*rms = 1.055). The other parameters are given in Table 1.

Ignoring the hole in the out-coupling mirror, the waist (0.59 mm) of the fundamental cold cavity mode is designed to be about √2 times the matched electron beam radius in the undulator. The FWHM of the fundamental cold cavity mode at the undulator exit and downstream mirror are 1.84 and 4.82 mm, respectively. We thus observe that the optical mode in the RAFEL expands faster from the undulator exit to the downstream mirror than the fundamental cold cavity mode (factor 3.25 and 2.62 respectively). That the RAFEL mode size is still smaller at the downstream mirror is the result of a balance between the faster expansion and smaller spot size of the RAFEL optical mode at the undulator exit compared to the cold cavity mode. The smaller spot size at the undulator exit is due to gain guiding as described above. Both the smaller spot size at the undulator exit and faster expansion of the RAFEL optical beam again indicate that at the undulator exit, the optical field consists of fundamental and higher order cold-cavity modes. Note, a fundamental Gaussian beam having a waist at the undulator exit with the same size as the RAFEL optical beam would have a Rayleigh range of 1.49 m. Finally, the cross section of the mode incident on the upstream mirror is shown in Figure 18c. After reflection from the upstream mirror and propagation through the resonator to the undulator entrance, this field results in a modal pattern similar to that shown in Figure 17.

(**c**) Incident on upstream mirror.

**Figure 18.** Cross-section of the optical field at various location on pass 60 for a cavity detuning Δ*L*cav = −18*λ* and an undulator field strength of *Bu* = 6.658 kG (*K*rms = 1.055). The other parameters are given in Table 1.

Now consider the formation of temporal coherence. Since the RAFEL starts from shot noise on the beam, the initial growth of the mode starts from spiky noise and, just as in a SASE FEL, develops coherence as it propagates through the undulator. However, the undulator is not long enough to reach saturation on a single pass, so that the evolution to temporal coherence is expected to develop over multiple passes. In Figure 19, we plot the temporal pulse shapes found on the first pass (a) at z = 0.5 m (blue) and 1.0 m (red), and (b) at z = 2.0 m (blue) and 2.4 m (red), the latter being at the undulator exit, for an on-axis undulator field amplitude *Bu* = 6.658 kG (*K*rms = 1.055). Figure 19 shows the full simulation time window, and it should be noted that the electron beam is centered in the time window with a full (parabolic) width of 1.2 ps. The pulse shows many spikes at 0.5 m, but that it has coalesced into about 5 spikes after 1.0 m. The development of temporal coherence continues, until at the undulator exit at 2.4 m, only two spikes remain.

The multi-pass development of temporal coherence after the first pass depends strongly on the cavity detuning. The temporal pulse shapes on the 60th pass at the undulator exit for *Bu* = 6.658 kG (*K*rms = 1.055) and for cavity detunings of Δ*L*cav = 0, <sup>−</sup>8*λ*, −18*λ* and −43*λ* are shown in Figure 20. As demonstrated previously, the synchronized cavity length for a RAFEL is shorter than the synchronized cavity length of a low-gain FEL oscillator. Therefore, as we have used the speed of light in vacuo, instead of the actual group velocity to define the cavity detuning, we note that Δ*L*cav = 0 actually corresponds to a cavity length that is larger than the synchronized length for the high-gain RAFEL. Consequently, the returning optical pulse will lag behind the center of the electron bunch. The pulse will be amplified as it propagates through the undulator, but as shown in Figure 19,

the two spikes formed over the first pass remain. This behavior is also found for small detunings, as shown in Figure 20b for Δ*L*cav = <sup>−</sup>8*λ*.

(**b**) Locations *z* = 2.0 (blue) and 2.4 m (red).

**Figure 20.** Temporal pulse shape at the exit of the undulator after 60 passes, and for various cavity detunings. The undulator field strength is *Bu* = 6.658 kG (*K*rms = 1.055). The other parameters are given in Table 1.

If the detuning is closer to the center of the detuning range, then the synchronism between the returning optical pulse and the electrons is a better match, and the multiple spikes are "washed out". This is shown in Figure 20c, where Δ*L*cav = −18 *λ* and we see that a broader pulse has formed. As the cavity length is decreased further, the optical pulse arrives increasingly near the head of the electron bunch. In this case, there is not enough gain to wash out the multi-spike character of the signal. This is shown in Figure 20d for Δ*L*cav = −43 *λ*, where the total power (or pulse energy is much reduced).

#### *5.2. An X-ray RAFEL*

The impetus driving research into X-ray free-electron laser oscillators (XFELOs) and RAFELs is the character of SASE emission. Since there are no seed lasers available at X-ray wavelengths, fourth generation FEL light sources [198,217–221] at these wavelengths are based on SASE over a single pass through a long undulator line. While pulse energies of the order of 2 mJ have been achieved at Ångstrom to sub-Ångstrom wavelengths, with the potential to reach multi-TW peak powers [222–227], SASE exhibits shot-to-shot fluctuations in the output spectra and power of about 10–20 percent. For many applications, these fluctuations are undesirable, and efforts are underway to find alternatives.

The utility of an XFELO has been under study for a decade [162,228–233], making use of resonators based upon Bragg scattering from atomic layers within diamond crystals [234–238]. The development of these crystals is a major breakthrough in the path toward an XFELO. Estimates indicate that using a superconducting RF linac producing 8 GeV electrons at a 1 MHz repetition rate is capable of producing 10<sup>10</sup> photons per pulse at a 0.86 Å wavelength, with a FWHM bandwidth of about 2.1 × 10−7. This design is consistent with the LCLS-II High Energy Upgrade [239]. As a consequence, an XFELO on a facility, such as the LCLS-II and LCLS-II-HE, is expected to result in a decrease in SASE fluctuations in the power and spectrum, and to narrow the spectral linewidth.

As with the majority of FELOs to date, the aforementioned XFELOs use low gain/high-*Q* resonators with transmissive out-coupling through thin diamond crystals [229]. Potential difficulties with low-gain/high-*Q* resonators derive from sensitivities to electron beam properties, mirror loading and alignments. In addition, transmissive out-coupling with high intra-cavity power can result in mirror damage. While experiments show that diamond crystals can sustain relatively high thermal and radiation loads [238], transmissive out-coupling cannot be easily achieved at the photon energies of interest here. Because of this, X-ray RAFELs [213,240,241] using a variety of out-coupling schemes are receiving a grea<sup>t</sup> deal of attention, and we consider a RAFEL design using a pinhole in one of the mirrors [210,213,240] here.

Optics propagation codes such as OPC must treat reflections from the diamond crystal Bragg mirrors, where the mirror losses and angles of reflection depend on the crystal orientation/geometry and the X-ray photon energy. X-ray Bragg mirrors typically have a very narrow reflection bandwidth and a narrow angle of acceptance [205]. For the X-ray RAFEL, and for computational efficiency, a temporal Fourier transform is applied at the beginning of the optical path, when the optical field is passed from MINERVA to OPC and the propagation is performed in the wavelength domain, i.e., each wavelength is independently propagated through the resonator using the modified Fresnel propagator (Equation (35)). The inverse Fourier transform is calculated at the end of the optical path, before the field is handed back to MINERVA. As the optical field inside the cavity is typically not collimated, a spatial Fourier transform in the transverse coordinates is calculated for each of the wavelengths when a Bragg mirror is encountered. Each combination of transverse and longitudinal wavenumber corresponds to a certain photon energy and angle of incidence on the Bragg mirror, and these parameters are used to calculate the complex reflection and transmission coefficients of the Bragg mirror [205]. After applying the appropriate parameter to the optical field, depending on whether it is reflected or transmitted, the inverse spatial Fourier transform is calculated and the field is propagated to the next optical element along the path until the end is reached.

Here, we consider a six-crystal, tunable, compact cavity [236], as illustrated in Figure 21, which shows both a top and a side view. The crystals are arranged in a non-coplanar (3-D) scattering geometry. There are two backscattering units comprising three crystals ( *C*1, *C*2, and *C*3) on one side of the undulator and three crystals ( *C*4, *C*5, and *C*6) on the other side. Collimating and focusing elements are shown as CRL1,2, which could be grazing incidence mirrors, but are represented in Figure 21 by another possible alternative—compound refractive lense [242,243]. In each backscattering unit, three successive Bragg reflections take place from three individual crystals to reverse the direction of the beam from the undulator. Assuming that all the crystals and Bragg reflections are the same, the Bragg angles can be chosen within the range 30º < *θ* < 90º; however, Bragg angles close to *θ* = 45º should be avoided to ensure high reflectivity for both linear polarization components, as the reflection plane orientations for each crystal change. The cavity allows for tuning the photon energy in a large spectral range by synchronously changing all Bragg angles. In addition, to ensure constant time of flight, the distance *L* (which brackets the undulator), and the distance between crystals, as characterized by *H*, have to be changed with *θ*. The lateral size *G* is kept constant as the resonator is tuned.

**Figure 21.** Schematic view of a 6-mirror ring resonator for an X-ray FEL oscillator using plane Bragg mirrors *C*1... *C*6 and compound refractive lenses, CRL1 and CRL2 for focusing. After [236].

Because the *C*1*C*6 and *C*3*C*4 lines are fixed, intracavity radiation can be out-coupled simultaneously for several users at different places in the cavity, although we only consider out-coupling through *C*6 at the present time. Out-coupling through crystals *C*3 and *C*6 are most favourable, since the direction of the out-coupled beams do not change with photon energy, but out-coupling for more users through crystals *C*1 and *C*4 is also possible. Such multi-user capability is in stark contrast with present SASE beamlines, which support one user at a time. We consider that the electron beam propagates from left to right through the undulator and the out-coupling is accomplished through a pinhole in the first downstream mirror ( *C*6).

Consider the LCLS-II beamline [239] with the HXR undulator corresponding to an electron energy of 4.0 GeV, a bunch charge in the range of 10–30 pC with an rms bunch duration (length) at the undulator of 2–173 fs (0.6–52 μm) and a repetition rate of 1 MHz. The peak current at the undulator is 1000 A, with a normalized emittance of 0.2–0.7 mmmrad, and an rms energy spread of about 125–1500 keV. The HXR undulator [239] is a plane-polarized, hybrid permanent magne<sup>t</sup> undulator with a variable gap, a period of 2.6 cm, and a peak field of 10 kG. Each HXR undulator has 130 periods. In the simulation, each undulator is modeled using Equation (17) and we consider that the first and last period

describe an entry/exit taper. There is a total of 32 segments that can be installed. The break sections between the undulators are 1.0 m in length and contain steering, focusing and diagnostic elements, although we only consider the focusing quadrupoles in the simulation, which we position in the center of the breaks. The quadrupoles are assumed to be 7.4 cm in length with a field gradient of 1.71 kG/cm.

A fundamental resonance at 3.05 keV (=4.07 Å) implies an undulator field of 5.61 kG. We assume that the electron beam has a normalized emittance of 0.45 mm-mrad and a relative energy spread of 1.25 × <sup>10</sup>−4, corresponding to the nominal design specification for LCLS-II. This yields a Pierce parameter of *ρ* = 5.4 × 10−4. In order to match this beam into the undulator/FODO line, the initial beam size in the *x*(*y*)-direction is 37.87 (31.99) μm with Twiss *αx* = 1.205 (*<sup>α</sup>y* = −0.8656). Note that this yields Twiss *βx* = 24.9 m and *βy* = 17.80 m.

The resonator dimensions were fixed by means of estimates of the gain using the Ming Xie parameterization [214], and MINERVA simulations indicated that about 40–60 m of undulator would be required to operate as a RAFEL. As such, the distance, L, between the two mirrors framing the undulator was chosen to be 130 m, which is also the distance separating the two mirrors on the back side of the resonator (elements *C*3 and *C*4). In studying the cavity tuning via time-dependent simulations, these two distances are allowed to vary while holding fixed the configurations of the backscattering units. The compound refractive lenses, which are modeled as thin lenses by OPC, are placed symmetrically around the undulator and are designed to place the optical focus at the center of the undulator in vacuo. In this study, the focal length is approximately 94.5 m.

In order to out-couple the X-rays through a transmissive mirror at the wavelength of interest, the diamond crystal would need to be impractically thin (about 5 μm); hence, we consider out-coupling through a hole in the first downstream mirror ( *C*6). We consider all the mirrors to be 100 μm thick. Due to the high computational requirements of time-dependent simulations, we begin with an optimization of the RAFEL with respect to the hole radius and the undulator length using steady-state (i.e., time-independent) simulations.

The choice of hole radius is important, because if the hole is too small then the bulk of the power remains within the resonator while if the hole is too large then the losses become too grea<sup>t</sup> and the RAFEL cannot lase. The results for the optimization of the hole radius indicate that the optimum hole radius is 135 μm, which allows for 90 percent out-coupling, where we fixed the undulator line to consist of 11 HXR undulator segments. This is shown in Figure 22, where we plot the output power as a function of pass number for the optimum hole radius and the variation in the saturated power with the hole radius (inset).

**Figure 22.** Output power versus number of passes through the resonator for a whole radius of 135 μm. The inset shows the output power as a function of hole radius. Other simulation parameters are described in the text.

A local optimization on the undulator length for a hole radius of 135 μm is shown in Figure 23, where we plot the peak recirculating power (left axis) and the average output power (right axis). The error bars in the figure indicate the level of pass-to-pass fluctuations in the power which is generally smaller than the level of shot-to-shot fluctuations in SASE. Note that while this represents steady-state simulations, the average power is calculated under the assumption of an electron bunch with a flat-top temporal profile having a duration of 24 fs which yields a duty factor of 2.4 × 10−8. Each point in the figure refers to a given number of HXR undulators ranging from 9–13 segments. It is evident from the figure that the optimum length is 47.18 m, corresponding to 11 segments.

**Figure 23.** Recirculating output power versus length of the undulator *Lu* for a hole radius of 135 μm. Other simulation parameters as described in the text.

Time-dependent simulations were performed under the assumption of electron bunches with a flat-top temporal profile having a full width duration of 24 fs and a peak current of 1000 A. This corresponds to a bunch charge of 24 pC. The detuning curve defining what cavity lengths are synchronized with the repetition rate of the electrons is shown in Figure 24. For simplicity, we take the synchronous cavity length *L*0 as *L*0 = *c*/ *f*rep (*cf.* Equation (46)), where a single optical pulse is assumed to be inside the cavity and the effect of gain on the group velocity is ignored. Here *L*0 = 299.7924580 m. The range of cavity lengths with overlap between the optical pulse and an electron pulse with a length of about 7.2 μm is given by *L*0 − 7.2 μm< *L*cav < *L*0 + 7.2 μm, where *L*cav is taken here to be the total roundtrip length of the cavity.

**Figure 24.** Output energy as a function of the cavity detuning Δ*L*cav = *L*cav − *L*0 for a hole radius of 135 μm. The simulation parameters are as described in the text.

The evolution of the output energy at the fundamental and the 3rd harmonic, and the spectral linewidth of the fundamental, vs. pass is shown in Figure 25 for a detuning of 5 μm, which is close to the peak in the detuning curve (Figure 24). While it is not evident in the figure, the rms fluctuation in the energy from pass to pass is of the order of about 3 percent, which is lower than the shot-to-shot fluctuations expected from SASE. At least as important as the output power is that the linewidth contracts substantially during the exponential growth phase and remains constant through saturation. Starting with a linewidth of about 3.7 × 10−<sup>4</sup> after the first pass corresponding to SASE, the linewidth contracts to about 6.0 × 10−<sup>5</sup> at saturation. The SASE linewidth after the first pass through the undulator is slightly smaller than the predicted linewidth based on 1-D theory [240], which is approximately 4.5 × 10−4. Hence, the RAFEL is expected to have both high average power and a stable narrow linewidth, which is to a large extent determined by the spectral filtering of the Bragg mirrors.

**Figure 25.** Evolution of the fundamental pulse energy (blue, left) and the 3rd harmonic (green, left), as well as the relative linewidth (red, right) when the cavity detuning is Δ*L*cav = 5 μm and for an electron beam energy spread Δ*Eb*/*Eb* = 1.25 × 10−4. Other parameters as described in the text.

The 3rd harmonic grows parasitically from high powers/pulse energies at the fundamental in a single pass through the undulator [184] and has been shown to reach output intensities of 0.1 percent that of the fundamental in a variety of FEL configurations, and this is what we find in the RAFEL simulations. As shown in Figure 25, the 3rd harmonic intensity remains small until the fundamental pulse energy reaches about 1 μJ after which it grows rapidly and saturates after about 12 passes. This is close to the point at which the fundamental saturates as well. The saturated pulse energies at the 3rd harmonic reach about 0.067 μJ. Given a repetition rate of 1 MHz, this corresponds to a long-term average power of 67 mW.

The reduction in the linewidth after saturation shown in Figure 25 indicates that a substantial level of longitudinal coherence has been achieved in the saturated regime. The RAFEL starts from shot noise on the beam during the first pass through the undulator, and longitudinal coherence develops over the subsequent passes. Hence, we expect that the temporal profile of the optical field will exhibit the typical spiky structure associated with SASE at the undulator exit after the first pass and this is indeed depicted in Figure 26a which shows the temporal profile at the undulator exit after the first pass. The number of spikes expected, *<sup>N</sup>*spikes, is given approximately by *<sup>N</sup>*spikes = *lb*/(2*πlc*), where *lb* is the rms bunch length and *lc* is the coherence length. For the present case, *lb* = 7.2 μm and *lc* = 60 nm; hence, we expect that *<sup>N</sup>*spikes = 19. We observe about 14 spikes in Figure 26a which is in reasonable agreemen<sup>t</sup> with the expectation. Note that the time axis encompasses the time window used in the simulation. As indicated in Figure 25, the linewidth after the first pass is of the order of 4.3 × 10−<sup>4</sup> which is relatively broad and corresponds to the interaction due to SASE. This is reflected in the output spectrum from the undulator after the first pass, which is shown in Figure 26b.

**Figure 26.** Temporal pulse shape and spectrum of the optical pulse at the exit of the undulator for the first pass. Other parameters as described in the text.

The spectral narrowing that is associated with the development of longitudinal coherence as the interaction approaches saturation results in a smoothing of the temporal profile. This is illustrated in Figure 27a, where we plot the temporal profiles of the optical field at the undulator exit corresponding to passes 12–16, which are after saturation, has been achieved (left axis). As shown in the figure, the temporal pulse shapes from pass-to-pass are relatively stable and exhibit a smooth plateau with a width of about 23–24 fs, which corresponds to, and overlaps, the flat-top profile of the electron bunches, which is shown on the right axis. Significantly, the smoothness of the profiles corresponds with the narrow linewidth and contrasts sharply with the SASE output after the first pass through the undulator (see Figure 26b). Both the pass-to-pass stability and smoothness of the output pulses contrast markedly with the large shot-to-shot fluctuations and the spikiness expected from the output pulses in pure SASE.

(**a**) Temporal profile for the optical pulse after different number of passes and current pulse (red).

(**b**) Spectrum at the exit of the undulator after 16 passes

**Figure 27.** Temporal pulse shape for the electron bunch and optical pulse after different number of passes and the spectrum at pass number 16 of the optical pulse at the exit of the undulator. Other parameters as described in the text.

> The narrow relative linewidth in this regime of about 7.3 × 10−<sup>5</sup> at the undulator exit, as shown in Figure 27b after pass 16, as well as the smooth temporal profiles shown in Figure 27a, are associated with longitudinal coherence after saturation is achieved.

#### **6. Summary and Conclusions**

The development of both low-gain/high-*Q* and high-gain/low-*Q* oscillators is an active and ongoing area of research into FEL sources over virtually the entire electromagnetic spectrum [30,86,91,92,241,244–248]. Both types of FEL oscillators have their own merits. For example, high-gain/low-*Q* are typically better suited for higher peak power operation, due to a reduced mirror loading [203] and higher extraction efficiency, as found in Section 5. On the other hand, properties like gain bandwidth, and stability of the output, are closer to SASE FELs than low-gain/high-*Q* FEL oscillators. However, we have not ye<sup>t</sup> explored the full range of feedback, and increasing the feedback may shift the performance in this respect closer to that of low-gain/high-*Q* oscillators.

As we have described in this paper, grea<sup>t</sup> progress has been made in our ability to simulate both the interaction through the undulator and in multi-pass configurations through the undulator and the resonator. This was demonstrated in Section 4 by the excellent agreemen<sup>t</sup> found between the simulations and the IR Demo and 10-kW Upgrade experiments at JLab, after which the simulations of both an infrared and an X-ray RAFEL were described in Section 5.

The combination of a three-dimensional, time-dependent FEL simulation code with an optics propagation code provides enormous flexibility in being able to study a wide range of undulator and optical resonator configurations. In particular, the use of a general FEL code such as MINERVA allows the simulation of variable polarization (Section 2.1) governed by the choice of undulator geometries (Section 2.2) and ultra-short electron bunches [249] down to the neighborhood of the cooperation length in the undulators, which defines the limit of applicability of the slowly varying amplitude and phase approximation. The use of OPC permits the simulation of arbitrary optical resonators, including thermal heating and distortion as well as Bragg reflections.

In summary, the numerical algorithms that have been implemented for the simulation of FEL oscillators are now sufficiently mature to serve as reliable design tools for virtually any conceivable undulator and oscillator configuration.

**Author Contributions:** Both authors contributed equally to conceptualization, methodology, software, validation, analysis, writing—original draft preparation, writing—review and editing, and visualization. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work described was supported over the course of time by multiple contracts with the U.S. Office of Naval Research, the Department of Defense Joint Technology Office, and the Department of Energy.

**Data Availability Statement:** Data is available from the authors on request.

**Acknowledgments:** The authors would like to thank numerous staff members at JLab for many helpful discussions about the IR Demo and 10-kW Upgrade experiments. We would also like to thank Heinz-Dieter Nuhn and Paul Emma for providing assistance with understanding the data from the LCLS and thank Luca Giannessi for assistance concerning the SPARC data. The work described was supported over the course of time by multiple contracts with the U.S. Office of Naval Research, the Department of Defense High Energy Laser Joint Technology Office, and the Department of Energy.

**Conflicts of Interest:** The authors declare no conflict of interest. Furthermore, the funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
