**1. Introduction**

This article describes an important mechanism associated with the physics of the free electron laser (*FEL*) regarding the nonlinear harmonic generation. The topic will be discussed by entering deeply into the mathematical aspects of what is currently referred as the high gain *FEL* equation.

A *FEL* is a laser device employing a beam of relativistic electrons injected into an undulator magnet, where it undergoes transverse oscillations and emits bremsstrahlung radiation, inversely proportional to the electrons energy square [1]. The *FEL* devices can be operated in the oscillator configuration by storing the emitted radiation in an optical cavity (whenever appropriate mirrors to confine the "light" are available). The radiation moves back and forth inside the resonator and interacts, at each entrance inside the undulator, with a freshly injected e-beam, thereby becoming amplified. The intracavity intensity grows, round trip after round trip; it is amplified by the process of stimulated bremsstrahlung emission and exhibits the same pattern of conventional laser oscillators; it undergoes exponential growth, and eventually saturation occurs, when the gain is reduced by the nonlinear contributions to the level of the cavity losses.

In the region of the spectrum where efficient mirrors are not available, cavity-less operation is an obliged step. In this regime, the gain of the system should be sufficiently large to bring the device to saturation in one undulator passage. These *FEL* devices, nowadays devoted to the production of intense X-ray beams, require long undulators (hundreds of meters) and high energy (up to tens of GeV)/high brightness electron beams (namely, beams with high intensity peak current, low energy, angular and transverse spatial dispersions). The "single pass" *FEL* devices may be operated in the amplifier configuration, when coherent input seeds are available; if not, the system works in the so-called self-amplified spontaneous emission (*SASE*) regime. The electrons emit ordinary synchrotron radiation at the entrance of the undulator, which reinteracts with the electrons, becomes amplified and eventually reaches saturation.

The paradigmatic steps leading the *FEL* to the generation of coherent laser-like radiation are the beam energy modulation, the bunching, the exponential growth and the saturation.

These phases are common to all other free electron devices (such as gyrotrons and coherent auto-resonance masers), and most of the relevant theoretical and mathematical

**Citation:** Dattoli, G.; Di Palma, E.; Licciardi, S.; Sabia, E. Free Electron Laser High Gain Equation and Harmonic Generation. *Appl. Sci.* **2021**, *11*, 85. https://dx.doi.org/10.3390/ app11010085

Received: 21 October 2020 Accepted: 17 December 2020 Published: 24 December 2020

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

**Copyright:** © 2020 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/).

descriptions can be framed within the same context [2]. Many *FELs* are presently operating and have operated in the past [3]. They are covering the electromagnetic spectrum from *THz* region to X-rays. The last generation of *FEL* sources (sometimes referred as fourth generation synchrotron radiation sources) has opened up unprecedent scenarios in different research fields because of the radiation characteristics, in terms of brightness (ten orders of magnitude above the ordinary sources) and the shortness of the pulses [4–7].

The growth of the *FEL* field is ruled, in the small signal regime (namely, the dynamical conditions, when the field intensity is not large enough to induce nonlinear effects), by an equation known as the *FEL* integral equation; it is derived from the linearization of the pendulum equation, reported below [8–16].

$$\begin{aligned} \frac{d^2}{d\tau^2}\zeta &= \mid a \mid \cos\left(\zeta + \phi\right) \,, & \qquad \frac{d}{d\tau}a &= -j\left(e^{-i\frac{\pi}{b}}\right) ,\\ \pi &= \frac{z}{L'} & \qquad j &= 2\pi g\_0 \,, & \qquad a = \mid a \mid e^{i\phi} \,, & \qquad \nu &= \frac{d}{d\tau}\zeta \end{aligned} \tag{1}$$

where *z* is the longitudinal coordinate; *L* is the undulator length; *g*0 is the small signal gain coefficient; and *a* and *j* are Colson's dimensionless amplitude and current respectively, *ζ*, *ν* the *FEL* longitudinal phase space variables and the brackets, denote the average of the phase space distribution. The use of dimensionless variables has several advantages. They constitute a set of comprehensive quantities such as *a* and *g*0, allowing one to merge all the significant *FEL* parameters into a single variable, providing quantities of interest (such as the gain and saturation intensity; see below) crucial in the descriptions and designs of *FEL* devices. The relevant small signal limit is achieved through the following approximations [17,18].

(a) From the first of Equation (1), we obtain the lowest order (in the field strength) expansion of *ζ*, namely

$$
\zeta \simeq \zeta\_0 + \nu\_0 \tau + \delta \zeta, \qquad \qquad \delta \zeta = \int\_0^\tau (\tau - \tau') \text{Re} \Big( a(\tau') e^{i(\zeta\_0 + \nu\_0 \tau')} \Big) d\tau' \tag{2}
$$

where *ν*0 = 2*πN ω*0−*ω ω*0 is the detuning parameter.

(b) Inserting *δζ* in the second of Equation (1) and averaging on the electron-field phase *ζ*0, we ge<sup>t</sup>

$$\begin{split} \frac{d}{d\tau}a &= -2\pi g\_0 b\_1 e^{-i\nu\_0 \tau} + i\pi g\_0 \int\_0^\tau (\tau - \tau') a(\tau') e^{-i\nu\_0(\tau - \tau')} d\tau' + \\ &+ i\pi g\_0 b\_2 \int\_0^\tau (\tau - \tau') a^\*(\tau') e^{-i\nu\_0(\tau + \tau')} d\tau' \end{split} \tag{3}$$

where the averages have been taken by considering the e-beam mono-energetic and without spatial and angular dispersion, thereby getting

$$\begin{aligned} \langle \kappa(\zeta\_0) \rangle\_{\zeta\_0} &= \frac{1}{2\pi} \int\_0^{2\pi} f(\zeta\_0) \kappa(\zeta\_0) d\zeta\_0, \\\ f(\zeta\_0) &= \sum\_{n=-\infty}^{\infty} b\_n e^{in\zeta\_0}, \qquad \qquad b\_n = \langle e^{-in\zeta\_0} \rangle\_{\zeta\_0}. \end{aligned} \tag{4}$$

The physical meaning of Equation (3) is transparent; the three terms on the right-hand side account for three different regimes:

(i) Absence of an initial bunching (*f*(*ζ*0) constant); in this case Equation (3) reduces to

$$\frac{d}{d\tau}a = i\pi \text{g}\_0 \int\_0^\tau (\tau - \tau') a(\tau') e^{-i\imath\_0(\tau-\tau')} d\tau'.\tag{5}$$

The underlying physics is that energy modulation and consequent bunching are due to the input coherent seed *a*0.

(ii) Nonconstant *f*(*ζ*0) and emergence of non zero *bn* coefficients; the field may grow in a seedless mode, induced by the initial bunching coefficient, as illustrated below.

In order to understand the role of the bunching coefficients during the early stages of the *FEL* interaction, we must treat Equation (5) in terms of a naïve expansion in the small signal gain coefficient, limiting ourselves to the second-order expansion in the small and neglecting *b*2; then we find

$$\begin{aligned} a(\tau) &\simeq a\_0 + \mathcal{g}\_0 a\_1 + \mathcal{g}\_0^2 a\_2, & a\_0(\tau) = 0, \\ \frac{d}{d\tau} a\_1 &= -2\pi b\_1 e^{-i\nu\_0 \tau}, & \frac{d}{d\tau} a\_2 = i\pi \int\_0^\tau (\tau - \tau') a\_1(\tau') e^{-i\nu\_0(\tau - \tau')} d\tau' \end{aligned} \tag{6}$$

and the relevant solution is

$$a(\tau) \simeq -2\pi b\_1 \lg \left( \frac{\sin(\frac{\nu\_0 \tau}{2})}{\frac{\nu\_0}{2}} e^{-i\nu\_0 \frac{\tau}{2}} - \frac{1}{2} \pi g\_0 \frac{(i\nu\_0^2 \tau^2 + 4\nu\_0 \tau - 6i)e^{-i\nu\_0 \tau} + 6i + 2\nu\_0 \tau}{\nu\_0^4} \right). \tag{7}$$

According to the previous equation, the field grows initially because of being triggered by the bunching coefficient. It provides a kind of coherent spontaneous emission, which in the second phase is responsible for the onset of the exponential growth regime (see Figure 1), where we have reported the combined effects of bunching and initial seed.

**Figure 1.** Square modulus of the dimensionless amplitude vs. the detuning parameter for *g*0 = 0.5 and different values of *τ*: (0.5, 0.6, 0.7, 0.8, 0.9, <sup>1</sup>). The intensity grows with incresing dimensionless time.

The previous results have been obtained using a perturbative expansion, which is not strictly necessary, since the *FEL* integral equation can be reduced to the ordinary differential equation reported below.

$$\left[\mathcal{D}\_{\tau}^{3} + 2i\upsilon\_{0}\mathcal{D}\_{\tau}^{2} - \upsilon\_{0}^{2}\mathcal{D}\_{\tau}\right]a(\tau) = i\pi\mathfrak{g}\_{0}\Big(a(\tau) + b\_{2}\ a^{\*}(\tau)e^{-2i\imath\_{0}\tau}\Big), \qquad \mathcal{D}\_{\tau} = \frac{d}{d\tau},\tag{8}$$

obtained after noting that integrals of type

$$i\pi g\_0 \int\_0^\tau (\tau - \tau') a(\tau') e^{-i\eta \cdot (\tau - \tau')} d\tau' = i\pi g\_0 e^{-i\eta \cdot \tau} \int\_0^\tau \left( \int\_0^{\tau'} a(\tau') e^{i\eta \cdot \tau''} d\tau'' \right) d\tau' \tag{9}$$

appearing in Equation (6) are double integrals which can be eliminated by keeping two successive derivatives [17,18]. It is evident that, for negligible *b*2 (we have checked that the second-order bunching coefficient does not produce any appreciable contribution and can be safely neglected), Equation (8) reduces to a naïve third-order ordinary differential equation, leading to the exponential growth also referred as *FEL* instability, the characteristic of which occurs in any Free Electron device [19–25]. (This type of instability is common to any free electron device (including gyrotrons and cyclotron auto-resonance masers (*CARM*)). This was widely established before the studies developed in the previously quoted references. Some interesting papers within this respect are listed in the bibliography.)

It is worth noting that, even though not explicitly contributing to Equation (8), the bunching coefficient *b*1 appears in the initial conditions provided by

$$a(0) = a\_{0\prime} \qquad \bigcirc\_{} \quad \bigcirc\_{} a \mid\_{\left[\tau=0\right)} = -2 \text{ }\pi \text{ }\gspace{0.1cm} \text{D}\_{\tau} \quad \bigcirc\_{} \quad \bigcirc\_{} \mid\_{\left[\tau=0\right)} = 2 \text{ }\pi \text{ }i \text{ }\nu\_{0} \text{ }g\_{0} \text{ }b\_{1} . \tag{10}$$

We have so far fixed the mathematical formalism to treat the problems associated with the effect of the bunching on the high gain (and not only) *FEL* evolution. The paper outline is reported below.

In Section 2 we discuss solution methods to deal with either the integral and thirdorder differential equation.

In Section 3 we deal with the extension of the integral equation to the harmonic generation and include effects associated with nonlinear regime.

Section 4 is finally devoted to application of the formalism for the design of *FEL*, exploiting segmented undulator devices.

#### **2. Algorithmic and Analytical Solutions of the** *FEL* **Integral Equation**

The equations we have dealt with in the previous section describe the 1 *D* small signal *FEL* dynamics, with the assumption of an ideal, sufficiently long beam such that short pulse effects can be neglected. The last assumption allows one to ignore the possible detrimental effects to the slippage, which will be considered later in this paper.

Within the present context, they are sufficiently general to develop useful considerations on the role of a pre-bunched e-beam on the *FEL* dynamics. We have also noted that Equation (3) can be reduced to a straightforward cubic equation, that analytical solutions are possible and that interesting information is obtained by the use of a perturbative expansion. The analysis has been so far limited (see Equation (7)) to a second-order expansion in terms of the small signal gain coefficient. The perturbative solution of Equation (5) can be obtained by including the higher order terms in *g*0, which yields the following recursion:

$$\frac{d}{d\tau}a\_n = i\tau \int\_0^\tau (\tau - \tau')a\_{n-1}(\tau')e^{-i\imath\psi(\tau-\tau')}d\tau', \qquad \qquad a\_n(0) = \delta\_{n,0}.\tag{11}$$

We will go back to perturbative treatments in the forthcoming sections of this article; here we consider the use of the other means based allowing either algorithmic and analytical solutions.

Regarding the first, we consider a method outlined in a paper by F. Ciocci et al. [26]; the procedure will be referred to as the Ciocci algorithm (*CA*). The technique is used via the following steps:

(i) The solution of the cubic equation can be written as

$$a(\tau) = e^{-i\nu \otimes \tau} \sum\_{j=1}^{3} \kappa\_j e^{-i\delta \nu\_j \tau} \tag{12}$$

where *δνj* are the roots of the third-degree algebraic equation

$$
\delta\nu^2(\nu\_0 + \delta\nu) = \pi \mathfrak{g}\_0;\tag{13}
$$

(ii) The amplitudes *κj* are fixed by the initial conditions, as shown below.

$$\begin{aligned} a(0) &= a\_0 = \sum\_{j=1}^3 \kappa\_j, \\ \frac{d}{d\tau} a\Big|\_{\tau=0} &= -i\nu\_0 a\_0 - i\sum\_{j=1}^3 \kappa\_j \delta \nu\_j = -2\pi g\_0 b\_1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

(iii) Equation (12) is cast in the more convenient form

$$a(\pi) = e^{-i\mathbf{r}\cdot\mathbf{r}} \sum\_{m=0}^{\infty} \frac{(-i)^m}{m!} \alpha\_m \pi^m, \qquad \qquad \alpha\_m = \sum\_{j=1}^3 \kappa\_j \left(\delta \nu\_j\right)^m. \tag{15}$$

(iv) The coefficient *αm* is obtained from Equation (13). It can be written as

$$
\pi \mathbf{v}\_0 \left(\delta v\_j\right)^2 + \left(\delta v\_j\right)^3 = \pi \mathbf{g}\_0 \tag{16}
$$

which after multiplying by *κj* and summing on the index *j*, yields the recursion

$$
\mathfrak{a}\_3 = \pi \mathfrak{g}\_0 \mathfrak{a}\_0 - \nu\_0 \mathfrak{a}\_2. \tag{17}
$$

and which for the higher order terms (*m* > 2) is generalized as

$$
\mathfrak{a}\_{m} = \pi \mathfrak{g}\_{0} \mathfrak{a}\_{m-3} - \mathfrak{v}\_{0} \mathfrak{a}\_{m-1} \tag{18}
$$

The solution of the cubic equation via *CA* is easily implemented numerically; it is extremely fast and yields the results contained in Figure 2 reporting the cases of seeded operation and bunching.

The intensity evolution vs. *z* for a seedless evolution (*α*(0) = 0, *b*1 = 0) and for the amplifier configuration is reported in Figure 2a,b. The *ν*0, *τ* surfaces yield an idea of the small signal intensity evolution, and further comments will be provided below.

Analytic solutions can be obtained too; they require some algebra, and in particular the use of Cardano's method to ge<sup>t</sup> the solution of Equation (13). Putting everything together we find:

(a) The Fang–Torre formula (this Formula appeared in an unpublished manuscript by H. Fang, and was later derived with minor refinements by A. Torre, and reported in [27]) valid for *a*0 = 0, *b*1 = 0:

$$\begin{split} a(\tau, v\_{0}) &= \frac{a\_{0}}{3(\nu\_{0} + p + q)} e^{-\frac{2\tau}{3}\eta\_{0}\tau} \left\{ (-\nu\_{0} + p + q) e^{-\frac{i}{\hbar}(p+q)\tau} + 2(2\nu\_{0} + p + q) e^{\frac{i}{\hbar}(p+q)\tau} \right. \\ & \cdot \left[ \cosh\left(\frac{\sqrt{3}}{6}(p - q)\tau\right) + i \frac{\sqrt{3}\nu\_{0}}{p - q} \sinh\left(\frac{\sqrt{3}}{6}(p - q)\tau\right) \right] \right\}, \\ p &= \left[ \frac{1}{2}(r + \sqrt{d}) \right]^{\frac{1}{2}}, \qquad q = \left[ \frac{1}{2}(r - \sqrt{d}) \right]^{\frac{1}{2}}, \\ r &= 27\pi\varepsilon g\_{0} - 2v\_{0}^{3}, \qquad d = 27\pi\varepsilon g\_{0} (27\pi\varepsilon g\_{0} - 4v\_{0}^{3}). \end{split} \tag{19}$$

The square modulus of *<sup>a</sup>*(*<sup>τ</sup>*, *<sup>ν</sup>*0) yields the intensity growth as a function of the dimensionless time and of the detuning parameter. Regarding Equation (19) we find

$$\begin{split} |a(\mathbf{r},\boldsymbol{\upnu}\_{0})|^{2} &= \frac{|\,^{\text{d}}\boldsymbol{\upnu}\_{0}|^{2}}{9(\boldsymbol{\upnu}\_{0}+\boldsymbol{s}\_{+})^{2}} \Bigg{{}^{2} \Big{(}4(2\boldsymbol{\upnu}\_{0}+\boldsymbol{s}\_{+})^{2} \Bigg{[}\Big{(}\boldsymbol{\upnu}\_{0}\boldsymbol{s}\_{-}\big{)}^{2} + \frac{3\boldsymbol{\upnu}\_{0}^{2}}{s\_{-}^{2}}\sinh\Big{(}\frac{\sqrt{3}}{6}\boldsymbol{s}\_{-}\big{)}^{2} \Bigg{]} + \\ &+ (-\boldsymbol{\upnu}\_{0}+\boldsymbol{s}\_{+})^{2} + 4(-\boldsymbol{\upnu}\_{0}+\boldsymbol{s}\_{+})(2\boldsymbol{\upnu}\_{0}+\boldsymbol{s}\_{+}) \Bigg{[}\cos\left(\frac{\boldsymbol{s}\_{+}}{2}\boldsymbol{\uppi}\right)\cosh\Big{(}\frac{\sqrt{3}}{6}\boldsymbol{s}\_{-}\big{)} + \\ & \qquad - \frac{\sqrt{3}\boldsymbol{\upnu}\_{0}}{s\_{-}}\sin\Big{(}\frac{\boldsymbol{s}\_{+}}{2}\boldsymbol{\uppi}\big{)}\sinh\Big{(}\frac{\sqrt{3}}{6}\boldsymbol{s}\_{-}\big{)} \Bigg{]} \Bigg{]}, \\ \boldsymbol{\upnu}\_{0} &= \boldsymbol{u} + \boldsymbol{e} \end{split} \tag{20}$$

*s*± = *p* ± *q*.

We have very loosely described in introductory remarks the *SASE* regime or seedless operation. In this case the field grows from a kind of prebunching characterizng the beam itself [28,29].

(b) The field growing from a bunching coefficient associated with the electron distribution and the solution of the evolution problem reads:

$$\begin{split} a(\tau,\upsilon\_{0}) &= \frac{2\tau\xi\_{0}b\_{1}}{\Delta}e^{-i\frac{2}{3}\upsilon\_{0}\tau} \left\{ \mu e^{-\frac{i}{3}(p+q)\tau} + \\ & - \left(\chi\_{-}\cosh\left[\frac{\sqrt{3}}{6}(p-q)\tau\right] + \chi\_{+}\sinh\left[\frac{\sqrt{3}}{6}(p-q)\tau\right]\right) e^{\frac{i}{6}(p+q)\tau} \right\}, \\ \Delta &= \frac{\sqrt{3}}{9}(p^{3}-q^{3}), \qquad \mu = i\frac{\sqrt{3}}{9}[\upsilon\_{0}-(p+q)](p-q), \\ \chi &= \frac{1}{6}\left[\upsilon\_{0}+\frac{1}{2}(p+q)-i\frac{\sqrt{3}}{2}(p-q)\right]\left[p+q+i\frac{\sqrt{3}}{3}(p-q)\right], \qquad \chi\_{\pm} = \chi\pm\chi^{\*}. \end{split} \tag{21}$$

Before proceeding further, we write Equation (21) in a form more appropriate for high gain/*SASE FEL* operation. To that end, we remind the reader that the small signal gain coefficient and Pierce parameter *ρ* (see references [8–16]) are linked by the identity

$$
\rho = \frac{\sqrt[3]{\pi \mathcal{G}\_0}}{4\pi N},
\tag{22}
$$

where *N* is the number of undulator periods. The intensity field growth is ruled by the so-called gain length

$$L\_{\mathcal{S}} = \frac{\lambda\_{\mathcal{U}}}{4\pi\sqrt{3}\rho'} \tag{23}$$

which is the pivotal quantity controlling the small signal high gain dynamics; its inverse defines the growth per unit length. With these remarks in mind, we note that:

(i) Quantities such as *ν*0 and *τ* should be replaced by

$$\nu\_0 \to \nu\_0 = \frac{\nu\_0}{4\pi N \sqrt{3}\rho} = \frac{1}{2\sqrt{3}\rho} \frac{\omega\_0 - \omega}{\omega\_0}, \qquad \qquad \tau \to \tilde{z} = \frac{z}{L\_\circ} \tag{24}$$

and it is worth noting that the number of undulator periods does not appear anymore in the definition of the new dimensionless quantities (see below for further comments);

(ii) Quantities involving the product of *p*, *q* and *τ*, written in the new variables, read

$$p\tau = \sqrt{3} \left( \sqrt[3]{1 + \frac{\sqrt{1 - 4\vec{v\_0^3}}}{1 - 2\vec{v\_0^3}}} \right) \overline{z}, \qquad q\tau = \sqrt{3} \left( \sqrt[3]{1 - \frac{\sqrt{1 - 4\vec{v\_0^3}}}{1 - 2\vec{v\_0^3}}} \right) \overline{z}. \tag{25}$$

It is finally important to stress that the Colson's dimensionless amplitude in terms of dimensional quantities reads

$$a = 2\sqrt{2}\pi\sqrt{\frac{I}{I\_s}}\tag{26}$$

where *I* is the field intensity (power/surface) and *Is* is the saturation intensity [30–33], a quantity which, in the theory of lasers, is exploited as a reference quantity to fix the onset of saturation. Even though its use is more appropriate for *FEL* operating in the oscillatory configuration, it can be exploited within the context of cavity-less operation too. Its dependence on the *FEL* parameters is specified by (practical units)

$$\left[I\_s \left[\frac{MW}{cm^2}\right] = 6.9312 \cdot 10^2 \left(\frac{\gamma}{N}\right)^4 \left(\lambda\_u [cm] K f\_b(\xi)\right)\right]^{-2}, \qquad \zeta = \frac{1}{4} \left(\frac{K^2}{\left(1 + \frac{K^2}{2}\right)}\right). \tag{27}$$

A more appropriate expression for the high gain/*SASE* regime can be obtained by expressing the second of Equation (1) in the new dimensionless coordinates, thereby finding

$$\begin{split} \frac{d}{d\tilde{z}}\tilde{a} &= -\frac{2}{\sqrt{3}} \langle e^{-i\tilde{\xi}} \rangle, \\ |\det| &= 2\sqrt{2}\pi\sqrt{\frac{I}{I\_s}} \qquad \qquad \qquad I\_s \left[\frac{MW}{cm^2}\right] = 6.9312 \cdot 10^2 \frac{(4\pi\gamma\rho)^4}{\left[\lambda\_u \left[cm\right] K f\_b(\xi)\right]^2}. \end{split} \tag{28}$$

For actual values of the *FEL* parameters entering in the definition of ˜*Is*, we obtain for it reference numbers not dissimilar from the *FEL* saturated power, as discussed later in this paper.

The use of the previous variables allows one to cast the intensity evolution in the dimensional form (valid for *ν*0 = 0) [17,18]

$$\begin{aligned} A(\tilde{z}) &= \frac{1}{9} \left[ 3 + 2 \cosh\left(\frac{z}{L\_{\mathcal{S}}}\right) + 4 \cos\left(\frac{\sqrt{3}}{2} \frac{z}{L\_{\mathcal{S}}}\right) \cosh\left(\frac{z}{2L\_{\mathcal{S}}}\right) \right], \\ I(\tilde{z}) &= I\_0 A(\tilde{z}), \quad \qquad \qquad I\_0 = \frac{I\_s}{8\pi^2} \parallel a\_0 \parallel^2. \end{aligned} \tag{29}$$

which has been largely exploited in the past [17,18] to describe the small signal *FEL* evolution; early derivations of similar equations trace back to the early eighties of the 20th century [34].

In Figure 3 we have plotted Equation (29) and the mere exponential behavior, which allows the conclusion that the region before the onset of the exponential growth, usually called lethargic region, lasts about two gain lengths, namely,

$$\mathbf{10^{-2}}\begin{bmatrix}\hline\\\hline\\\hline\\\hline\\\hline\\\mathbf{10^{-3}}\\\hline\\\mathbf{10^{-3}}\\\hline\\\hline\\\hline\\\hline\\\hline\\\hline\\\hline\\\hline\\\hline\\\hline\\\end{bmatrix}\mathbf{10^{-3}}\mathbf{10^{-3}}$$

$$Z\_L \simeq 2.1 \ L\_{\S}.\tag{30}$$

**Figure 3.** Comparison between Equation (29) and the pure exponential behavior for *Lg* = 0.4.

In any *FEL* device the small signal gain is a crucial parameter; it is naively defined as the relative intensity variation:

$$G = \frac{|\|\vec{a}(\vec{z}, \vec{v}\_0)\|^2 - |\|\vec{a}\_0\|^2|}{\|\vec{a}\_0\|^2}. \tag{31}$$

In the case of low gain, its maximum value is proportional to the small signal gain coefficient through the well known identity [17,18]

$$G^\* = 0.85 \,\text{g}\_{0^\circ} \tag{32}$$

Such a quantity measures the amplification of the input seed at the end of the undulator. We can, however, extend the concept by considering the same gain measurement at different points inside the undulator, as shown in Figure 4, after a number of periods

$$N\_z = \frac{z}{\lambda\_u}.\tag{33}$$

We can therefore write *G*∗ in terms of *ρ*, *z* as (see Equation (32))

$$G^\* = \frac{0.85}{\pi} \left(\frac{z}{\sqrt{3}}\right)^3. \tag{34}$$

The low gain regime relies on the assumption that the laser amplitude remains constant during the transit inside the undulator. The approximation breaks when higher order terms in *g*0 are to be taken into account to specify the maximum gain at a given point in *z*. These corrections, derived in [35], can be written as

$$\begin{split} G^\* &= \frac{|\,\tilde{a}(\tilde{z}, \tilde{v}\_0^\*)\,|^2 - |\,\tilde{a}\_0\,|^2}{|\,\tilde{a}\_0\,|^2} \simeq \\ &\simeq \frac{0.85}{\pi} \left(\frac{\tilde{z}}{\sqrt{3}}\right)^3 + \frac{0.19}{\pi^2} \left(\frac{\tilde{z}}{\sqrt{3}}\right)^6 + \frac{4.23 \cdot 10^{-3}}{\pi^3} \left(\frac{\tilde{z}}{\sqrt{3}}\right)^9 + o(\rho^{12}) \end{split} \tag{35}$$

whose meaning is better explained in Figure 4, where we have also reported the departure of the gain line shape from the anti-symmetric curve, characterizing the early *FEL* experiments [36,37].

**Figure 4.** Power density (a.u.) vs. *z* and gain function "measured" at different points inside the undulator. Above the last point, where the dotted curve starts, the small signal gain is not properly defined because the small signal approximation does not strictly apply. The curves have been derived for the parameters *ρ*0 = 6.8 · 10−<sup>4</sup> and *λu* = 5 cm. The small signal gain coefficient goes from 5.35 · 10−<sup>3</sup> to 5.35.

In this section we have provided a fairly general analysis describing the solution of the *FEL* high gain small signal equation.

In the forthcoming section we discuss how the technique we have discussed so far can be extended to the treatment of *FEL* induced higher order harmonics.

#### **3. High Gain FEL Equations and Harmonic Generation**

In this section we discuss the mechanism of harmonic generation in *FEL* devices, which is a process "naturally" associated with the evolution of the fundemental harmonic. The importance of these further contributions in developing tools extending the capabilities of *FELs* was soon recognized in the community, but the practical implementations took some time.

The pendulum equations, including the higher order harmonics, can be written as (differently from the index *n* in the perturbative terms discussed in Equation (11); in this case *n* is just the order of the harmonic) [38]

$$\begin{aligned} \frac{d^2}{d\tau^2} \zeta &= \sum\_{n=0}^{\infty} |\!/a\_n \mid \cos(\psi\_n)\,, & \qquad \psi\_n = n\zeta + \phi\_{n\prime} \\ \frac{d}{d\tau} a\_{\hbar} &= -j\_{\hbar} \langle e^{-in\zeta} \rangle\_{\prime} \qquad & j\_{\hbar} = 2\pi \text{g}\_{0,n} \end{aligned} \tag{36}$$

where *an* =| *an* | *eiφ<sup>n</sup>* is the dimensionless amplitude of the *n*-th order harmonics and *g*0,*n* is the *n*-th harmonic small signal gain parameter to be specified later [39].

$$\begin{split} \frac{d}{d\tau}a\_n &= -2\pi g\_0 b\_n e^{-i\nu\_n \tau} + i\pi \eta g\_n \int\_0^\tau (\tau - \tau') a\_n(\tau') e^{-i\nu\_n(\tau - \tau')} d\tau', \\ \nu\_n &= 2\pi N \frac{n\omega\_0 - \omega}{n\omega\_0}. \end{split} \tag{37}$$

It is evident that the same considerations developed for the fundamental harmonic (*n* = 1), hold in the present case too, the only significant difference being that the *n*-the harmonic Pierce is now specified by

$$\rho\_n^\* = \sqrt[3]{n}\,\rho\_{n\prime} \qquad \rho\_n = \rho \left(\frac{f\_{b,n}}{f\_b}\right)^{\frac{2}{3}}, \quad f\_{b,n} = f\_{\frac{n-1}{2}}(n\_b^{\gamma}) - f\_{\frac{n+1}{2}}(n\_b^{\gamma}), \qquad f\_{b,1} = f\_b. \tag{38}$$

The gain length determining the harmonic linear growth is

$$L^\*\_{\mathbb{g},n} = \frac{\lambda\_{\mathbb{f}}}{4\pi\sqrt{3}\rho^\*\_{\mathbb{n}}} = \frac{1}{\sqrt[3]{n}} L\_{\mathbb{g},n}.\tag{39}$$

where *Lg*,*<sup>n</sup>* is defined in terms of *ρn* as

$$L\_{\mathbb{S}^M} = \frac{\lambda\_u}{4\pi\sqrt{3}\rho\_n} \tag{40}$$

The role of *<sup>L</sup>*<sup>∗</sup>*g*,*<sup>n</sup>* in the harmonic generation is clarified below; it is just a consequence of the assumption (summarized in Equation (36), where the electron field phase *ψn* is written as *nζ*) that the process is dominated by the bunching induced by the fundamental harmonic.

It is worth noting that the correction term √3 *n* in the definition of the harmonic Pierce parameter was initially suggested in references [17,18,39] as a consequence of an accurate analysis of the numerical data from the code PROMETEO [40].

In Figure 5 we have reported the harmonic intensity growth vs. the longitudinal coordinate; the relevant behavior is the same as for *n* = 1, reproduced by equations of the type (29) with *<sup>L</sup>*<sup>∗</sup>*g*,*<sup>n</sup>* in place of *Lg*.

**Figure 5.** Harmonic evolution and nonlinear harmonic generation which is triggered after the dashed line, where the bunching effects become significant.

However, this is only a part of the story. The harmonic intensities (3, 5) follow a growth rate characterized by two different regimes. After the lethargic region, the growth rates is initially specified by the inverse of the harmonic gain length *L*∗ *g*,*<sup>n</sup>* and then by an abrupt increase of the growth ruled by *nL*−<sup>1</sup> *g* . This behavior is the signature of a complex mechanism underlying the harmonic generation emission.

We have so far considered the harmonic linear dynamics, in which within certain limits (to be discussed in the final section), the evolutions of the associated intensities are independent. It should be noted that harmonics too can be triggered by a bunching coefficient, as implicitly contained in Equation (37). As it is well known, the bunching grows with the intensity itself; the dynamics of each harmonic are the result of intrigued feedback between induced bunching and higher order-bunching (namely, at higher values of *n*). We crudely describe the process by noting that each harmonic, while the intensity grows, induces further bunching, determining the onset of the sub-harmonics *m* of the *n*-th harmonic—namely, the first harmonic indices higher order harmonics and each of them determines a further set of harmonics. In terms of bunching coefficients, we have: The fundamental harmonics (*n* = 1) contributes to the increase of the bunching of order 1 and of other higher order bunching coefficients. The same holds for the higher order harmonics; the third is, e.g., responsible for the bunching of its own harmonics (3, 6, 9, ... with respect to the fundamental).

The bunching process is not "free," and indeed an energy spread, associated with gain degradation and saturation mechanisms, is produced. This phenomenology is not considered within the present context.

The inclusion of nonlinear terms in our analysis occurs by keeping intensity dependent terms [17,18,41] which allow the following redefinition of the *FEL* integral equation:

$$\begin{split} \frac{d}{d\tau}a &= -2\pi g\_0 e^{-i v\_0 \tau} \sum\_{m=-\infty}^{\infty} (-i)^m b\_{1-m} \frac{J\_m(|\,\,\Pi\,\vert)}{|\,\,\Pi\,\vert^m} \left( \int\_0^\tau (\tau-\tau') a(\tau') e^{i v\_0 \tau'} d\tau' \right)^m, \\ \int\_0^\tau (\tau-\tau') a(\tau') e^{i v\_0 \tau'} d\tau' &= |\,\,\Pi\,\vert \, e^{i\psi} \end{split} \tag{41}$$

where

$$J\_m(\lfloor \ell \Pi \rfloor) = \sum\_{r=0}^{\infty} \frac{(-1)^r}{r!(m+r)!} \left(\frac{\lfloor \ell \Pi \rfloor}{2}\right)^{m+2r}.\tag{42}$$

If we consider a perturbative solution of the above equation in terms of the fundamental intensity, we can write the first three terms as

$$\begin{split} \frac{d}{d\tau}a\_{l} &= -2\pi g\_{0}b\_{1}e^{-i\nu\_{0}\tau} + i\pi g\_{0} \int\_{0}^{\tau} (\tau - \tau')a\_{l}(\tau')e^{-i\nu\_{0}(\tau-\tau')}d\tau', \\ \frac{d}{d\tau}a\_{nl,2} &= -2\pi g\_{0}e^{-i\nu\_{0}\tau} \frac{b\_{-1}}{8} \left( \int\_{0}^{\tau} (\tau - \tau')a\_{l}(\tau')e^{i\nu\_{0}\tau'}d\tau' \right)^{2}, \\ \frac{d}{d\tau}a\_{nl,3} &= -2\pi g\_{0}e^{-i\nu\_{0}\tau} \frac{b\_{-2}}{2^{3}3!} \left( \int\_{0}^{\tau} (\tau - \tau')a\_{0}(\tau')e^{i\nu\_{0}\tau'}d\tau' \right)^{3} \end{split} \tag{43}$$

which have been derived using the limit lim |Π|→0 *Jm*(| Π |) | Π |*m* = 1 *m*!2*<sup>m</sup>* . The subscripts *l* and *nl* stand for linear and nonnlinear, respectively. The *nl*-contributions stand for the onset of the nonnlinear harmonic generation.

Considering a linearly polarized undulator not providing coupling to even on axis harmonic generation, we obtain from the third Equation (43)

$$\begin{split} a\_{nl,3} &= -2\pi g\_0 \frac{b\_{-2}}{2^3 3!} \int\_0^\tau e^{-i\mathbf{v}\_0 \cdot \mathbf{r}'} h\_3(\mathbf{v}\_0, \mathbf{r}') d\tau' \, , \\ h\_3(\mathbf{v}\_0, \mathbf{r}) &= \left( \int\_0^\tau (\tau - \tau') a\_0(\tau') e^{i\mathbf{v}\_0 \cdot \mathbf{r}'} d\tau' \right)^3 \, . \end{split} \tag{44}$$

Keeping the fast growing root only and assuming *ν*0 = 0, we find a further contribution to the harmonic intensity depending on the power of the fundamental and characterized (as already anticipated) by the gain length

$$L\_{\rm g}^{\left(n\text{lkg}\right)} = \frac{L\_{\rm g}}{n} \tag{45}$$

where *nlhg* stands for nonnlinear-harmonic generation.

The conclusion of this discussion is that there are two distinct contributions to the higher order harmonic emission. The first is a kind of lasing on harmonic; this concept was originally introduced in [42] by Colson et al., where the possibility of exploiting this mechanism to sustain higher order harmonics lasing in the oscillator configuration has been discussed. The second provides the well-known nonlinear harmonic generation mechanism; it goes beyond the linear analysis; its characterizing features and scaling relations are further discussed in the final section of this paper.

We have commented on the importance of the harmonic generation mechanism, which has opened the possibility of extending the relevant performances in terms of tunability range and not only.

We want to mention an important application of the harmonic generation in *FEL* which has determined a significant improvement of the relevant performances. It has already been underscored that running the *FEL* in the amplifier configuration is feasible if a coherent seed at certain wavelengths is available. This happens with *FEL* harmonics too, which can be triggered using the seeds from harmonic generation in gas (see references [43,44] for an adequate description). In these devices, the beam of laser (a Ti:sapphire laser, for example), is focused into a xenon gas cell, where high harmonic generation occurs. The output seed beam is then spatially and temporally overlapped to the electron beam moving in an undulator with the period chosen to match the seed frequency. The *FEL* harmonics then grow according to the mechanisms we have just outlined. The successful implementation of this concept has allowed the possibility of extending the *FEL* tunability of more than an order of magnitude (for more details, see [43] and references therein). A further idea, concerning the use of segmented undulators, after the beam self induced prebunching, is discussed in the forthcoming section.

#### **4. Inhomogeneous Broadening Partial Amplitudes and High Gain** *FEL* **Equation**

The high gain integral equation becomes slightly more complicated if the effect of non ideal e-beam qualities is included. The averages in Equation (2) need to be extended to the beam distribution (energy, spatial, angular, etc.). With these premises and limiting ourselves to the effect of the energy spread only, we find [45]

$$\frac{d}{d\tau}a = i\pi g\_0 \int\_0^\tau a(\tau - \tau')\tau' e^{-i\nu\_0 \tau' - \frac{1}{2}(\pi \mu\_t \tau')^2} d\tau', \qquad \mu\_\varepsilon = 4\text{N}\sigma\_\varepsilon \tag{46}$$

obtained after taking the average of a Gaussian relative energy distribution with r.m.s. *σε*.

The solution of Equation (46) cannot be obtained in analytical terms; however, the method of partial amplitude expansion pioneered by A. Segreto in references [46,47] yields a fairly useful analytical mean allowing the solution of, e.g., Equation (46) in the form (we use the *ν*˜0, *z*˜ variables, more appropriate for the forthcoming discussion)


$$\begin{split} \bar{a}(\boldsymbol{z},\boldsymbol{v}\_{0}) &= \bar{a}\_{0} \left( 1 + \sum\_{k=1}^{\infty} i^{k} \left( \frac{\boldsymbol{\bar{z}}}{\sqrt{3}} \right)^{3k} g\_{k}(\boldsymbol{z},\boldsymbol{v}\_{0},\boldsymbol{\bar{\mu}}\_{\ell}) \right), \\ \bar{g}\_{k}(\boldsymbol{z},\boldsymbol{v}\_{0},\boldsymbol{\bar{\mu}}\_{\ell}) &= \frac{M\_{k}\gamma\_{k}}{\Sigma\_{k}} \exp \left[ -\frac{(\bar{\mu}\_{\ell}\boldsymbol{\bar{z}})^{2} (\beta\_{k} + i\underline{\boldsymbol{v}}\_{k}(\boldsymbol{v}\_{0}\boldsymbol{\bar{z}})) + 12(i\lambda\_{k}(\bar{\boldsymbol{v}}\_{0}\boldsymbol{\bar{z}}) + (\bar{\boldsymbol{v}}\_{0}\boldsymbol{\bar{z}})^{2}) \right], \\ \boldsymbol{M}\_{k} &= \frac{1}{(3k)!}, \qquad \gamma\_{k} = (3k+1)\sqrt{\frac{3k+2}{2k(k+1)}}, \qquad \Sigma\_{k} = \sqrt{\gamma\_{k}^{2} + \frac{1}{12}(\bar{\mu}\_{\ell}\boldsymbol{\bar{z}})^{2}}, \\ \boldsymbol{\beta}\_{k} &= 2\frac{4k+1}{k+1}, \qquad \lambda\_{k} = 2(3k+1)\frac{3k+2}{k+1}, \qquad \tilde{\zeta}\_{k} = 2\frac{(2k+1)(k-1)}{(k+1)(3k+1)}. \end{split} \tag{47}$$

The previous equation can be profitably used to ge<sup>t</sup> information of a practical nature. The low gain expansion (*k* = 1) yields for the intensity growth

$$|\vec{a}(\vec{z}, \vec{v}\_0)|^2 \simeq |\vec{a}\_0|^2 \left[ 1 + 2 \frac{M\_1}{3\sqrt{3}\Sigma\_1} \vec{z}^3 \cdot \exp\left( -\frac{\beta\_1 \vec{\mu}\_\varepsilon^2 + \vec{v}\_0^2}{24\Sigma\_1^2} \vec{z}^2 \right) \sin\left( \frac{\vec{v}\_0 \left( 12\lambda\_1 + \xi\_1 (\vec{\mu}\_\varepsilon \vec{z})^2 \right)\_{\vec{z}}}{24\Sigma\_1^2} \vec{z} \right) \right] \tag{48}$$

which provides an idea of the interplay between inhomogeneous broadening and the other mechanisms contributing to the laser evolution. In Figure 6 we have reported the analogues of Figure 2b for different values of the inhomogeneous broadening parameter. The plot yields quite a good description of how the intensity growth is diluted by the presence of a non zero energy spread.

Other effects due to angular and spatial distribution of the electron beam can be introduced too. From the conceptual point of view, they do not imply any problem, but for the introduction of a further convolution term in Equation (46). The induced detrimental effects are expressed through appropriate *μ*-parameters, analogous to those associated with the energy spread. They are listed in [18].

It is important to stress that, from the practical point of view, the effect of non perfect electron beam qualities produces an increase of the gain length, which can be parameterized as [17,18]

$$L\_{\mathcal{K}}(\vec{\mu}\_{\varepsilon}) = \left(1 + 0.185 \frac{\sqrt{3}}{2} \vec{\mu}\_{\varepsilon}^{2}\right) L\_{\mathcal{K}} \tag{49}$$

The interest in the method of partial amplitudes stems also for the fact that it can be generalized to more complicated situations, involving, e.g., the pulse propagation contributions.

**Figure 6.** Intensity vs. *z*˜, *ν*˜ for inhomogeneous parameters *μ*˜*ε* = 1, *g*0 = 100, *a*0 = 1.

In the following we make the assumption that the electron bunch is sampled into a series of slices with a longitudinal size of the order of a coherence length [48–50]:

$$
\lambda\_{\mathfrak{c}} = \frac{\lambda}{4\pi\sqrt{3}\rho}.\tag{50}
$$

During the interaction the radiation is mismatched with respect to the electron bunch reference by a quantity depending on the slippage length. The computation of this aspect of the problem requires a certain amount of computation which, for the homogeneous case (no energy spread) has been accomplished in [51] and for the present case is summarized below.

According to the aforementioned reference and to the analytical steps outlined in [46,47], the problem is solved, in relatively easy terms, by replacing the detuning parameter with the operator

$$
\vec{v}\_0 \to \vec{v}\_0 = \vec{v}\_0 + i\frac{\lambda\_c}{L\_\mathcal{g}}\partial\_\zeta \tag{51}
$$

where *ζ* is the bunch coordinate. The replacement (51) in the partial amplitudes accounts for the action of the *FEL* dynamics on the optical packet, expressed by the simple Gaussian (to avoid any misunderstanding, we emphasize that *ζ* now indicates the optical bunch coordinate and should not be confused with the pendulum coordinate in Equation (1))

$$a\_0(\zeta) = \frac{1}{\sqrt[4]{2\pi\sigma\_{\zeta}^2}} e^{-\left(\frac{\zeta}{2\sigma\_{\zeta}^2}\right)^2}.\tag{52}$$

The action of the partial amplitudes on the initial packet is specified by

$$\begin{split} g\_{k}(\mathbb{Z}, \widehat{v}\_{0\prime}, \widehat{\mu}\_{\varepsilon}) a\_{0}(\zeta) &= g\_{k}(\mathbb{Z}, \mathbb{Z}, \widehat{v}\_{0\prime}, \widehat{\mu}\_{\varepsilon}) e^{\left(\operatorname{\mathbb{S}}\_{k} \partial\_{\zeta} + D\_{k} \partial\_{\zeta}^{2}\right)} a\_{0}(\zeta), \\ S\_{k} &= \frac{\lambda\_{\varepsilon}}{2 \Sigma\_{k}^{2} L\_{\mathbb{S}}} \left\{ \left[ \frac{(\pi \mu\_{\varepsilon} \varepsilon)^{2} \xi\_{k}}{12} + \lambda\_{k} \widetilde{\varepsilon} \right] - i \widetilde{v}\_{0} \right\}, \qquad \qquad D\_{k} = \frac{1}{2 \Sigma\_{k}^{2}} \left( \frac{\widetilde{\varepsilon} \lambda\_{\varepsilon}}{L\_{\mathbb{S}}} \right)^{2}. \end{split} \tag{53}$$

The action of the exponential operator containing the shift (first-order derivative) and diffusive operators (second-order derivative) on the initial packet is provided by

$$\mathcal{L}\left(\mathcal{S}^{\mathcal{S}\_{\mathbb{K}}+D\_{k}\partial\_{\zeta}^{2}}\_{\zeta}\right)a\_{0}(\zeta) = \frac{1}{\sqrt[4]{2\pi}}\sqrt{\frac{\sigma\_{\zeta}}{\sigma\_{\zeta}^{2}+D\_{k}}}\exp\left[-\frac{1}{4}\frac{\left(\zeta+\zeta\_{k}-i\frac{\tilde{v}\_{0}\lambda\_{c}\mathfrak{T}}{2\Sigma\_{k}L\_{\chi}}\right)^{2}}{\sigma\_{\zeta}^{2}+D\_{k}}\right].\tag{54}$$

A result which allows a few speculations:

(a) The wave packet is shifted back by a quantity

$$\zeta\_k = \frac{1}{24\Sigma\_k} \left[ (\pi\mu\_\ell)^2 \zeta\_n + 12\lambda\_k \right] \frac{\lambda\_\ell}{L\_\S} \,. \tag{55}$$

This means that the optical packet is moving back with respect to the bunch frame propagating at velocity *c*, the associated group velocity is

$$
\sigma\_{\mathcal{S},k} = \left(1 - \frac{1}{c} \frac{d}{dt} \zeta\_k\right) \mathbf{c}\,. \tag{56}
$$

The radiation moves slower than *c*, by the effect of the interaction itself. The gain dilution due to the energy spread counteracts the velocity slow down with respect to the case of negligible energy spread.

The complex refraction index derived from Equation (56) can be written as

$$m\_{\mathbb{S},k} \simeq 1 + \frac{1}{c} \frac{d}{dt} \zeta\_k \,. \tag{57}$$

(b) The optical packed undergoes a longitudinal diffusion specified by

$$
\sigma\_k = \sqrt{D\_k} = \frac{1}{\sqrt{2\Sigma\_k}} \left(\frac{\overline{z}\lambda\_\varepsilon}{L\_\S}\right) \tag{58}
$$

which ensures that in one gain length the additional packet width is essentially a coherence length. This effect has been studied in the theory of mode-locked *FEL* operation within the context of the oscillator theory; the relevant discussion can be found in [2].

In the forthcoming section we apply the result obtained so far to the slice evolution.

#### **5. Slice Evolution**

The concept of slice phase space is a by-product of the *SASE FEL* physics. It is indeed associated with the fact that, in these devices, the combination of mechanisms such as gain, slippage and finite coherence length, determines a local interaction, because the radiation experiences only a portion of the beam, having the longitudinal extension of a coherence length (see Figures 7 and 8). The interaction is therefore sensitive to the slice-brightness, which is characterized by the relevant six dimensional phase space distribution [52–58].

The analysis associated with the transverse phase space and with the relevant electronbeam transport properties has already been accomplished in a number of papers, which have clarified a significant deal of the physical and practical issues regarding the evolution and interplay between slices and *FEL* power output and performances.

In this section we clarify the effect of longitudinal phase space, by the use of the formalism we outlined in the previous section. We sample the optical bunch as indicated in Figure 7 and characterize each slice with a progressive number [59,60].

**Figure 7.** Slice number sampling and shape, using the e-bunch distribution as the reference frame.

**Figure 8.** Bunched envelope, associated slices and superimposed radiation.

Each slice is assumed to be Gaussian with an amount of charge *Qs* and an associated current

$$I\_s = \frac{Q\_s}{\sqrt{2\pi}\tau\_c}, \qquad \qquad \tau\_c \simeq \frac{L\_c}{c} \equiv coherence - time. \tag{59}$$

A natural assumption is that the amount of charge follows the Gaussian shape of Figure 8, namely,

$$Q\_{\mathfrak{s}} = Q\_{s^\*} \mathfrak{e}^{-(s-s^\*)^2 \delta^2}, \qquad \qquad \delta = \frac{2\pi \mathcal{L}\_{\mathfrak{c}}}{\sigma\_{\mathfrak{s}}} \tag{60}$$

with *δ*−<sup>1</sup> being the number of slices inside the bunch and *s*<sup>∗</sup> representing the slice with largest current. With these assumptions we can write the slice Pierce parameter as

$$
\rho\_s = \rho\_{s^\*} e^{-\frac{(s-s^\*)^2 \delta^2}{3}}.\tag{61}
$$

As we already stressed, each slice carries its own energy spread, and therefore we ge<sup>t</sup> a corresponding inhomogeneous parameter specified by

$$
\tilde{\mu}\_{\mathfrak{s}} = \tilde{\mu}\_{\mathfrak{s}^\*} e^{-\frac{(s-s^\*)^2 \delta^2}{3}}.\tag{62}
$$

The energy spread corresponding to each slice may be completely random and that with largest current may carry the worst spread.

By taking into account these effects, we obtain the results reported in Figures 9 and 10. The first refers to the case in which the slice with maximum current coincides with that at the center of the distribution, and the energy spread of each slice is chosen to be completely random. The second accounts for a configuration in which the maximum current is associated with a slice belonging to the rear part of the bunch.

Further comments including the elements for a more accurate computation of the slice distribution, are discussed in the forthcoming section dedicated to concluding remarks.

**Figure 9.** Slice evolution inside the undulator at different *Z*; the energy spread has been randomly chosen and the current follows the bunch Gaussian distribution with *s*∗ = 0, *λu* = 0.02, *ρ*0 = 0.001, *με* = 0.0004 and *ν* = 0. The initial distribution *z* = 0 of the slices is the same as in Figure 8. The vertical axis is in log-scale; the insets are plotted in linear scale.

**Figure 10.** Same as Figure 9—all slices with the same energy spread and the maximum current associated with the rear part of the bunch (*s*∗ = 1). Slice evolution vs *z* and intensity evolution vs. *z*. The slices are characterized by the same amount of charge and different values of the energy spread with *λu* = 0.015, *ρ*0 = 0.0011, *σε* = 0.05, *<sup>σ</sup>d* = 0.002, *δ* = 0.01 and *I*0 = 10−1. The intensity growth does not refer to a single slice but to the relevant average. The slices have a flat distribution in *z*, and each is characterized by an initial seed of 0.1 W. The vertical axis is in log-scale; the insets are plotted in linear scale.

#### **6. Final Remarks**

This paper has treated different topics in the small signal evolution of high gain *SASE* devices. We have just touched on the nonlinear effects leading to saturation. It has been underscored that the relevant pattern is accompanied by mechanisms of higher order bunching leading to the process of higher order harmonic emission.

From the phenomenological point of view, saturation can be included by modifying the small signal power evolution in Equation (29) as

$$I\_{\mathbb{S}}(\tilde{z}) = I\_0 \frac{A(\tilde{z})}{1 + \frac{I\_0}{I\_F}(A(\tilde{z}) - 1)}, \qquad I\_{\mathbb{F}} = \sqrt{2}\rho P\_{\mathbb{E}} \tag{63}$$

where *PE* is the electron beam power intensity and *I*0 is the input seed intensity. *I*0 is usually specified by the rule of thumb that a value usually 10<sup>8</sup> below the saturated power density is chosen. The intensity *IF* yields the saturated power. The term at the denominator accounts for the saturation mechanisms, using a kind of logistic model, as displayed in Figure 11. The last figure needs its own mention: it shows the power growth in the deep saturated regime. Equation (63) reproduces well the evolution up the maximum and then remains fixed without following the oscillations (see Figure 12) due to the post saturation dynamics, characterized by e-beam rotation in phase-space, induced by an exchange of power between laser field and electrons.

A more complete view is reported in Figure 13. The physical content can be noted as reported below.

The figure exhibits the first and third harmonics, along with the associated wave packet evolution (upper and lower panels). We have assumed the growth from a coherent seed, and no spiking appears before the saturation. With the onset of the power oscillations, the optical packets start to be characterized by the appearance of side bands, which degrade the laser pulse itself in correspondence with the minima of the oscillations [61–63]. The effects of beam qualities can be also included by replacing *<sup>A</sup>*(*z*) with approximants discussed in the previous section.

**Figure 11.** Same parameters as in Figure 10 with the inclusion of the saturation. The slices have a flat distribution in *z*, and each is characterized by an initial seed of 0.1 W. The vertical axis is in log-scale; the insets are plotted in linear scale.

**Figure 12.** Power growth and oscillation after saturation (no specific parameters of the simulation have been inserted because we are just interested to the oscillating behavior, characterizing any *SASE FEL* operating with constant parameters).

**Figure 13.** First and third harmonic power growth vs. *z*, along with the power oscillations and associated changes in the laser packet distribution.

The pivotal element of the discussion of this paper has been centered around the *FEL* high gain equation, which, including the necessary modifications, is able to describe different physical situations regarding the *FEL* phenomenology. The *FEL* high gain equation has been derived from the *FEL* pendulum equation after a suitable linearization. A different treatment can be, however, exploited: it employs a Hamiltonian picture and the associated Liouville equation. The by-product of this treatment (see, e.g., [9,10,17,18,49,50,52–60,64]) yields results not dissimilar from those discussed in the previous sections, and the bunching

coefficient of the *n*-th harmonic induced by the fundamental behaves like | *a* |*n*. Such an effect is of fundamental importance in the handling of the so-called segmented magne<sup>t</sup> configurations [50,65]. *FEL* adopting this type of solution exploits the growth of the bunching inside a first section, which is interrupted and connected to a second undulator, tuned at one of its harmonics. The result is the growth of a coherent signal emerging from the beam bunching acquired in the first section. Methods based on the joint use of the solution of the high gain *FEL* equation, on the study of the relevant Liouville equation and on the use of massive simulation codes, allow the possibility of deriving useful scaling relations concerning the onset and saturation of the *FEL* signal from a suitably pre-bunched beam [65].

$$\begin{aligned} \Pi\_{\mathfrak{n}}(z) &= \Pi\_{\mathbb{O},\mathfrak{n}} \frac{e^{\frac{z}{L\_{\mathfrak{g}}^{(n)}}}}{1 + \frac{\Pi\_{\mathbb{O},\mathfrak{n}}}{\Pi\_{\mathbb{F},\mathfrak{n}}} \left(e^{\frac{z}{L\_{\mathfrak{g}}^{(n)}}} - 1\right)} \mathsf{I} & \qquad \begin{aligned} \mathsf{L}\_{\mathfrak{g}}^{(n)} &= \frac{L\_{\mathfrak{g}}}{n}, \\ \Pi\_{\mathfrak{g}}^{(n)} &= \frac{L\_{\mathfrak{g}}}{n}, \end{aligned} \tag{64}$$
 
$$\Pi\_{\mathbb{O},\mathfrak{n}} = c\_{\mathfrak{n}} \left(\frac{P\_{\mathbb{O}}}{\mathfrak{g}\rho P\_{E}}\right)^{n} \Pi\_{\mathbb{F},\mathfrak{n}\prime} \qquad \qquad \Pi\_{\mathbb{F},\mathfrak{n}} = \frac{1}{\sqrt{n}} \left(\frac{f\_{\mathbb{O},\mathfrak{n}}}{n!f\_{\mathbb{O},\mathbb{I}}}\right)^{2} P\_{\mathbb{F}} \end{aligned} \tag{64}$$

where the coefficients *cn* are just numerical characteristics for each harmonic. Equation (64) yields the growth of the nonlinear part of the *n*-th harmonic. The relevant gain length is decreased by a factor corresponding to the order of the harmonics; Π*F*,*<sup>n</sup>* is the harmonic saturated power. The number of photons/seconds at the *n*-th harmonic can accordingly be calculated as

$$\dot{N}\_n = \frac{\Pi\_{F,n}}{n\hbar\omega} = E\_n \dot{N}\_{1\prime} \qquad \qquad \qquad E\_n = \frac{1}{n\sqrt{n}} \left(\frac{f\_{b,n}}{nf\_{b,1}}\right)^2 \tag{65}$$

where *En* represents a kind of efficiency displaying how the fundamental harmonic contributes to the higher order harmonics.

We have so far given a general discussion of analytical and numerical methods to treat the *FEL* intensity evolution in the high gain/*SASE* regime. We have just touched on the problems associated with short pulse effects, which will be treated elsewhere in a dedicated monography.

**Author Contributions:** Conceptualization, G.D., E.D.P., E.S.; methodology, G.D., E.D.P., E.S. software, E.D.P., E.S.; validation, G.D., E.D.P., S.L., E.S.; formal analysis, G.D., E.D.P., S.L., E.S.; data curation, G.D., E.D.P., S.L., E.S. writing–original draft preparation, G.D., E.D.P., E.S.; writing–review and editing, S.L.; visualization, G.D., E.D.P., S.L., E.S. supervision, G.D., E.D.P., S.L., E.S.; project administration, E.D.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The work of S. Licciardi was supported by an Enea Research Center individual fellowship.

**Conflicts of Interest:** The authors declare no conflict of interest.
