Bunched beam superradiance is the coherent spontaneous radiation emission of a bunch of electrons, taking place when the bunch duration σ t is shorter than the optical period ( 2 π / ω ) of the radiation [1]. This process is analogous to the superradiance of an ensemble of dipole-excited molecules proposed first by Dicke [2]. When this condition is satisfied, the radiation emission of the bunch is proportional to the number of electrons in the bunch, squared (N²). This is a substantial enhancement of many orders of magnitude in comparison to spontaneous emission by a long duration beam that is proportional to the number of electrons in the bunch (N).

Here we extend earlier theory of bunched-beam superradiance [1, 3], and apply it to evaluate and optimize the spectrum and energy of a compact THz FEL shown schematically in Figure 1. This experimental setup is based on the design of the Israeli hybrid 6 MeV photocathode RF gun [4]. A picosecond electron bunch emitted from the cathode is accelerated and chirped within the short (0.6 m) hybrid accelerator, and then compressed through free drift compression, velocity bunching (ballistic compression) [5–7] along a 3.4 m beamline to σ t ≈ 70 f s , as predicted by GPT simulations including space charge effects [8, 9] (see the lower panel of Figure 1). The superradiant THz pulse is generated in an over-moded rectangular waveguide placed within a short planar magnetic undulator (0.8 m).

We present here the formulation for calculating the spectral energy and total emission energy of superradiant waveguide FEL. Using the formulation of [1], the spectral energy per emitted radiation mode q in the case of a perfectly bunched electron beam (zero length bunch—farther down we generalize to a finite length bunch) is: d W q d ω S P − S R = N 2 e 2 Z q 16   π a ¯ w β z γ 2 L w 2 A e f f , q sinc 2 θ L w 2 (1) and its integrated energy is: W q = N 2 e 2 Z q 16   π a ¯ w β z γ 2 L w 2 A e f f , q ω 0 , q N w (2)

These equations are valid for any complete orthogonal mode expansion in free space or in a waveguide. In the case of conventional free space FEL Z q = μ 0 / ε 0 , in a waveguide it is the waveguide mode impedance. The effective area of the mode is A e f f , q = ∬ E q ⊥ x , y 2 d x d y E q x x 0 , y 0 2 , where (x₀, y₀) are the transverse coordinates of the beam propagation axis. In these equations N is the number of particles, c is the speed of light, e is the electron charge, γ is the beam energy Lorentz factor, the rms average of the undulator parameter a_w = e B w m c k w is a ¯ w (m is the electron mass, B_w is magnetic field amplitude). β_z is the longitudinal velocity of the beam inside the undulator, L_w is the undulator length, N_w is the number of undulator periods. The FEL detuning parameter is defined by: θ = ω v z − k z , q − k w (3)

Where k w = 2 π λ w is the periodic undulator wavenumber, k z , q = ω c 2 − k ⊥ , q 2 is the longitudinal wavenumber of mode q, where k ⊥ , q is the transverse wave number of the mode.

The spectral energy of the radiation mode (Eq. (1)) depends on frequency through the tunning parameter θ q ω . The center frequency of the spectral radiation (1) corresponds to θ q ω = 0 . In a waveguide, this equation has two solutions [10]. We refer here to the higher frequency (+) solution only: ω 0 , q = γ z 2 β z c k w 1 ± β z 2 − k ⊥ , q γ z k w 2   (4)

The formulation so far is valid for any waveguide mode expansion. We consider an FEL configuration based on a rectangular waveguide and planar undulator polarized in the y direction (see Figure 2). For this configuration we find it most productive to use a waveguide expansion set of LSM (Longitudinal Section Magnetic) modes [11]. This is an alternative to the conventional TE TM mode expansion. The LSE^x modes (characterized by E_x = 0) and LSM^x modes (characterized by H_x = 0), are a complete set of orthogonal modes, equivalent to the {TE, TM} mode expansion. In a planar undulator configuration, shown in Figure 2, the wiggling of the electron beam, and thus the excitation current, are in the x dimension. Therefore, in the {LSE^x, LSM^x} expansion, we can eliminate the LSE^x modes, since they cannot be excited by the beam current, and thus we manage to describe the radiation field in terms of half the number of modes of the degenerate {TE, TM} mode expansion. We therefore calculate only the excitation of the LSM^x modes given by [11]: E z = ⁡ sin k x x sin k y y E x = k y 2 + k z 2 i k x k z cos k x x sin k y y E y = − k y i k z sin k x x cos k y y H z = − ω ϵ k y β k x cos k x x cos k y y H y = − i ω ϵ k x cos k x x sin k y y (5) k z , m n = ω c 2 − k ⊥ , m n 2 k ⊥ , m n = k y n 2 + k x m 2 k xm = m π a   m = 0 , 1 , 2 . . . k yn = n π b   n = 0 , 1 , 2 . . . Z LSM = Z 0 k y 2 + k z 2 k k z (6) ω is the frequency of the radiation mode, ε is the permittivity inside the waveguide. Z₀ = μ 0 / ε 0 is the impedance of free space k = ω / c , k_z,mn (ω) is the longitudinal wavenumber of the mode (also the dispersion relation of the mode). For an axi-symmetric electron beam propagating along the waveguide axis, we can refer to the symmetry of the modes and eliminate all the modes that are null on axis. Therefore, only the modes m = 0,2,4… n = 1,3,5…that have finite amplitude E x a 2 , b 2 ≠ 0 can be excited.

Substituting the fields of the excitable modes (5) in the excitation Equation 1, we calculated the spectral energy of the emitted radiation for the first lower 9 modes of the waveguide. For the designed example parameters, listed in Table 1, the radiation energies of these modes are listed in the third row of Table 2.

In practice the beam transverse dimension and its longitudinal bunching are not ideal. We considered here also the case of a beam of finite transverse and longitudinal dimensions assuming gaussian distribution. The effects of the transverse distribution of the beam can be taken into account by convolving Eq. (1) with the normalized transverse distribution of the beam profile, which gets modified because of the (x₀, y₀) dependances of the effective mode area A e f f , q . This was calculated in Supplementary Appendix B for the LSM^x modes. In Figure 3 we show the spectrum of the first 9 modes of the waveguide computed from Eq. (1) with the parameters of Table 1 including the effect of transverse beam dimensions for σ⊥ = 0.5 [mm]. The bunch is considered to be ideally bunched in the longitudinal dimension (σ_t = 0). In the fourth row of Table 2 we show the reduced emission energy of the excited modes assuming a transverse beam profile of standard deviation σ ⊥ = 0.5 m m . The emission energy is reduced relative to the ideally narrow beam σ ⊥ = 0 given in the third line.

The total emitted superradiant energy and spectral energy of the modes would be reduced also by the finite longitudinal dimension of the beam. This is a result of destructive interference between the radiation wave packets emitted by distributed electrons when the short bunch condition is not satisfied. The spectral energy and total superradiant radiant energy are reduced by a bunching factor relative to the ideal zero-length bunch (eq. 1,2). For a bunch longitudinal profile represented by a Gaussian distribution of standard deviation σ t emitting radiation at frequency f, the bunching factor is given by [1]: M b f | 2 = e − 2 π f σ t 2 (7) and the mode energy is: W m n σ t = W m m σ t = 0 M b 2 . (8)

The bunching factor Eq. 7 is plotted in Figure 4 as a function of the emission frequency for various values of σ t . Clearly, the drop of the bunching parameter at high frequencies determines the upper frequency limit of the superradiant source. The figure indicates that for a bunch of standard deviation 70 fs, the superradiant Terahertz FEL is limited to operation below f = 3 THz.

We compare the above analytical calculation results for the superradiant energy emission with numerical computations based on the UCLA GPTFEL code [12, 13] for a rectangular waveguide model. We performed the computations for the parameters of the ORGAD FEL configuration (Table 1). The comparison is made for the fundamental mode of the rectangular waveguide TE01 which is identical with the LSM01 mode used in the analytical calculation. Figure 5 displays the numerically simulated spectral energy of the mode TE01. The frequency spectrum in this model is presented in turns of longitudinal modes [12], and the beam at the entrance to the undulator is modeled in terms of a gaussian distribution. The spectral energy distribution computed by the numerical code is in good agreement with the corresponding curve of mode LSM01.

In Table 3 we compare the analytically calculated emission energy of the fundamental mode energy to the numerical simulation results for different values of the standard deviation of the beam bunch distribution σ_t. The calculated bunching factor (Eq. 7) and reduced emission energy (Eq. (8)) are listed in the third and fourth columns of Table 3. The numerically computed results, listed in the fifth column, are about 50% higher. They also indicate diminishing of the emission energy of the fundamental mode (centered at f = 3.4 THz) for σ t ≈ 70   f S .

Since the high frequency operation of the superradiant terahertz source is limited by the size (duration) of the electron beam bunch within the length of the undulator, it is important to control the free drift self-compression of the beam from the gun to the undulator and attempt to minimize the bunch duration nearly at the center of the undulator. In Figure 6 we show the evolution in phase space of a Gaussian beam starting from the buncher (modulator) section of the gun, right after full acceleration, and then going through free drift transport up to the undulator. In our linear beam transport model, the phase-space area of the beam is conserved under linear transport transformations.

We model the initial distribution in γ-t phase-space of the accelerated beam (before bunching) in terms of a Gaussian function: f i γ , t = 1 2 π σ t 0 σ γ 0 e − t − t 0 2 2 σ t 0 2 + − γ − γ 0 2 2 σ γ 0 2 = 1 2 π σ t 0 σ γ 0 e − ϕ γ , t (9) σ t 0 is the beam standard deviation duration, at the entrance of the gun. σ γ 0 is the intrinsic (uncorrelated) energy spread (standard deviation) in the gun. As shown in Figure 6 the beam is chirped in the modulation section of the gun at a rate S = Δ γ / Δ t . Following this step, the beam is compressed through free space drift and passage through the undulator. The compression is characterized by the longitudinal dispersion factor (compaction parameter) R₅₆. After the modulation transformation (chirping in the gun): γ − γ 0 →   Δ γ + S Δ t (10) ϕ γ , t = t − t 0 2 2 σ t 0 2 + Δ γ + S Δ t 2 2 σ γ 0 2 (11)

We Consider beam transport from the accelerator to the center of the undulator. In a linear phase-space dynamics model the beam phase-space distribution stretches in the time dimension according to t − t 0 → Δ t + R 56 c   Δ γ γ 0 (12)

The energy dispersion compaction parameter R₅₆ is the result of free drift along a distance L from the gun to the undulator entrance and a subsequent drift within periodic magnetic field along half the length of the undulator: R 56 L a L w 2 = 1 γ 0 2 L a ∫ 0 L w 2 a ⊥ 2 z ′ d z ′ a ⊥ z ′ = e / m c ∫ 0 z B ⊥ z ′ d z ′ (13) where B_┴ is the undulator magnetic field. Substituting Eq. (12) in Eq. (11) we obtain: ϕ L , γ , t = Δ t + R 56 c Δ γ γ 0 2 2 σ t 0 2 + Δ γ + S Δ t + R 56 c Δ γ γ 0 2 2 σ γ 0 2 (14)

We receive an oblique ellipse contour: ϕ L , γ , t = 1 2 σ t 0 2 + S 2 2 σ γ 0 2 Δ t 2 + 1 2 σ γ 0 2 + R 56 2 2 c 2 γ 0 2 σ t 0 2 + S 2 R 56 2 2 c 2 γ 0 2 σ γ 0 2 − R 56 S c γ 0 σ γ 0 2 Δ γ 2 + R 56 c γ 0 σ t 0 2 − S σ γ 0 2 + R 56 S 2 c γ 0 σ γ 0 2 Δ t Δ γ (15)

Requiring attainment of a minimal duration bunch at the center of the undulator, we observe that this happens when the ellipse, keeping constant area, becomes erect. This corresponds to a requirement that the mixed term in Eq. 15 nulls. This corresponds to a condition for R₅₆: R 56 = c γ 0 σ t 0 2 S S 2 σ t 0 2 + σ γ 0 2 (16) or alternatively, in order for the beam waist to fall inside the undulator one needs to satisfy a condition on the modulation (chirp rate) coefficient: S = c γ 0 1 + 1 − 4 R 56 2 σ γ 0 2 c 2 γ 0 2 σ t 0 2 2 R 56 ≈ c γ 0 R 56 (17)

At the waist, the ellipse is erect: ϕ w L w , γ , t = Δ t 2 2 σ t w 2 + Δ γ 2 2 σ γ w 2 (18) where σ t w  is the bunch duration at waist and σ γ w is the energy spread at waist.

Substituting (17) or (16) in (15), the radii of the erect ellipse are found to be: 1 σ t w 2 = 1 σ t 0 2 + S 2 σ γ 0 2 (19) 1 σ γ w 2 = 1 σ γ 0 2 + R 56 2 c 2 γ 0 2 σ t 0 2 + S 2 R 56 2 c 2 γ 0 2 σ γ 0 2 − 2 R 56 S c γ 0 σ γ 0 2 (20)

Since we expect compression - σ t w ≪ σ t 0 , we can approximate Eq. (19) a to: σ t w ≈ σ γ 0 S (21)

We conclude that the bunch duration at the waist is limited by the intrinsic energy spread (defined as the uncorrelated beam energy spread at its waist before modulation).

Note that in the framework of the linear model, the phase space area is conserved throughout the beam transport from the gun to the wiggler, σ t 0 σ γ 0 = σ t w σ γ w , as can be confirmed by substitution of equation 20, 19, or 21 and 17.

We follow the phase space evolution of the beam for the example of the ORGAD accelerator using the linear transformation phase-space model. Bear in mind that this model is limited, does not include space-charge effects and nonlinear phase space evolution. Also, it is assumed that the amplitude of the superradiant radiation field is low (below saturation), and the beam dynamics is independent of the emitted radiation, namely, the emitted radiation does not act back on the particles. Therefore the linear model can serve only as a preliminary guide for the design of the superradiant FEL, and further on it should be checked and compared to the results of detailed GPT transport simulations that include consideration of space-charge effects. Such simulations are shown in Supplementary Appendix A for the ORGAD parameters where the beam transport was adjusted to produce a ribbon beam in the undulator in order to reduce space-charge effects.

We calculate the compaction factor from the end of the gun to the center of the undulator. This is given as a sum of the free drift compaction factor from the gun to the undulator and the compaction factor of half the length of the undulator R₅₆ ^(a)+R^(b) ₅₆ = R₅₆. As derived in Supplementary Appendix C, for γ₀ = 13.1 and L_e = 3.24[m] L e is the drift length from the gun (z = 0) to the entrance of the undulator), R 56 a = L e γ 0 2 = 1.89 c m , R₅₆ ^(b) = 0.33[cm] for half of the undulator length z_hw = L w 2 = 40 c m , and a_w = 0.915. Finally, the compaction factor from the end of the gun to the center of the undulator is R₅₆ ^(a)+R^(b) ₅₆ = R₅₆ = 2.22[cm].

For these parameters we calculate the modulation coefficient desired for attaining a waist at the center of the undulator. The starting parameters for the linear phase-space evolution are taken from the GPT simulation of the gun section (see Supplementary Appendix A): the intrinsic energy spread is σ_γ0 = 0.01 (σ_E0 = 5 keV), and the beam bunch size at the gun is σ_t0 = 1[ps]. Using Eq. (17), this results in S = 1.77*10¹¹[1/sec], which is fairly consistent with the computed chirp rate shown in the last panel of Supplementary Appendix Figure SA1. According to eq. 21 and eq. (20), these parameters should enable attaining a beam waist of σ t w ≈ 60 [fs] at the center of the undulator length with σ γ w = 0.12 .

This result of an ideal linear transport model is an underestimate relative to the result of a full GPT simulations with the parameters of the ORGAD hybrid photocathode gun. The GPT simulation along the entire beam transport results in a beam bunch waist σ_tw = 70[fs] at the center of the undulator (see lower part of Figure 1 and Supplementary Appendix A). This is slightly bigger than the estimate of the linear model.

In order to compare the beam size evolution in the linear phase space model to the numerical GPT simulation we trace back the phase-space ellipse evolution from the waist location at the center of the undulator backward.

We start from Eq. 14 and assume that the phase-space ellipse is erect at the waist point in the middle of the undulator z = L 0 , evolves with z (within and before the undulator) according to: ϕ γ , t , z = t − t 0 + R 56 L 0 − R 56 z c   γ − γ 0 γ 0 2 2 σ t w 2 + γ − γ 0 2 2 σ γ w 2 = 1 (22)

The compaction factor dependence on z starting from the exit of the gun z = 0 is (Supplementary Appendix C): R 56 z = 1 γ 0 2 z a w 2 2 z − L e η z − L e (23) where η is the step function, and L 0 = L e + L w / 2 .

From here we derive the size of the beam (the projection of the ellipse on the time axis t [7]: σ t z = σ t w 1 + R 56 L 0 − R 56 z c γ 0 2 σ γ w 2 σ t w 2 1 2 (24)

This curve is shown in Figure 7 for the earlier computed parameters σ t w ≈ 60   fs , σ γ w = 0.12 .

In Figure 7 we also show the corresponding results of GPT simulations for the same parameters overlayed over the same coordinate axis. The simulated curve deviates approximately by 20% at the start of the drift section (exit of the gun). The discrepancy is attributed to space charge expansion of the bunch along the drift section and in the undulator.

Bunched beam superradiance is an attractive concept for attaining intense THz radiation in a compact FEL scheme [1, 3, 14–16]. Here we presented theory and design consideration of a superradiant waveguide FEL based on a compact (60 cm long) hybrid photocathode RF gun. For the configuration of the ORGAD accelerator, an 80 cm long undulator and a modest beam bunch charge of 20 pC, we predict for an ideally short bunch, emission of about W = 60 nJ radiation at the fundamental transverse radiation mode LSM01 at frequency 3.4 THz and nearly 1 μJ in 9 transverse modes in the frequency range 2.5–3.5 THz. To attain these energies, one is required to attain compression of the beam to short duration - σ t < 1 / 2 π f . For a nonoptimized design based on the ORGAD accelerator parameters and GPT simulations, including space-charge effects, a bunch duration of σ t = 70   f t S was calculated, which corresponds to reduction of the radiation emission by a bunching factor of 0.1 to 0.3 for the different modes. In the presented design the bunch duration was limited by space charge effects in the undulator section. It may be possible to mitigate these space-charge effects by proper shaping of the beam transverse and longitudinal profiles at the gun and at the entrance to the undulator [14] and thus restore the emission energies predicted for an ideal beam.

It is instructive to compare the parameters of the superradiant THz FEL design to the alternative scheme of enhancing spontaneous undulator radiation by SASE [17]. The superradiance enhancement factor relative to spontaneous emission is N = Q/e, the number of electrons in the beam. This factor is 10⁸ for Q = 20 pC. Larger enhancement factors are attainable by the SASE scheme, but this requires much longer undulator, higher beam energy and higher beam charge [17, 18]. In PITZ the measured terahertz energy was two orders of magnitude (tens of μ J ) higher in the same spectral range. However, this required threefold larger beam energy, fourfold longer undulator and two orders of magnitude higher charge. The bandwidth of an isolated single mode of the superradiant FEL design example is 2.5%, comparable to the SASE spectrum, but if one considers the multimode spectrum shown in Figure 3, it lies in the wide range 2.8–3.5 THz. For many applications of diagnostics and radiation effect research, particularly in university laboratories, only moderate radiation energy is required and accelerator dimensions and costs matter. For such applications, superradiant FEL may have an advantage.

We presented an analytical model evaluation for the evaluation of superradiannt energy and spectral distribution [1, 3], and compared the results to numerical computations using UCLA GPTFEL code. We found fair agreement between the radiation spectrum and the integrated radiation energy per mode computed using the code and using the analytical expression.

We presented a linear model for tracing the beam compression dynamics in energy-time phase-space evolution from the electron gun through free drift up to and within the undulator. This model is a useful tool for preliminary design of a superradiant FEL, however, it has limited validity because it does not take into consideration space-charge effects and nonlinear distortion of phase-space trajectories transformation. The model helped to recognize that the attainment of short bunches in the undulator and high frequency operation is limited also by the intrinsic (uncorrelated) energy spread of the gun. In realistic design, the predictions of the linear phase-space evolution model must be checked by numerical simulations that take into account apace charge effects. In the present paper we backed up the analytical model design example by GPT numerical computation that include space-charge effects. For the presented design example space charge effect deviations of the beam trajectories were moderate (20%).

Within the validity limits of the linear model, we found that the minimal temporal waist size of the beam σ t w in a given beam transport configuration is given by σ t w ≈ σ γ 0 S - proportional to the intrinsic energy spread of the gun σ γ 0 . Thus, to operate the superradiant FEL at higher frequency it is necessary to minimize the intrinsic energy spread of the gun. For the example based on the ORGAD accelerator, where space charge effects are marginal, we found that a minimal bunch temporal size of σ t w = 70   f s at the center of the undulator is attainable and is consistent with an intrinsic energy spread in the gun of σ_E0 = 5 keV. This is quite consistent with the GPT simulation results of the Gun (Supplementary Appendix Figure SA1). Deviation may be expected also because of the neglect of nonlinear effects in the phase-space evolution along the transport line that are not taken into account in the linear phase space model. These values of the intrinsic energy spread in the gun are similar to measurements of uncorrelated energy spread of photocathode RF-gun injectors in other laboratories [19]. Measurements at SwissFEL 110 m distance from the rf gun showed 15 keV for 200pC, and 6.5 keV for 10pC, while the simulation results predict well below 1 keV [20]. The discrepancy between the measurement and simulations is related to IBS (Intrabeam scattering) and MBI (microbunching instability) between the gun and the measurement point. Measurements at European XFEL, 40 m from the gun area, for beam charge of 250 pC show slice energy spread of 6 keV [21]. More recent publications reported slice energy spread of 2 KeV for 250pC, 20 m from the rf gun at the PhotoInjector Test facility at DESY Zeuthen (PITZ) [22].

The bunch duration in the ORGAD Accelerator design is also bound by the energy spread of the photocathode gun. GPT simulations in the gun section indicate intrinsic energy spread in the gun of ∼5 keV for a beam of charge of 20 pC that should be attributed to space charge effects in the early acceleration stages in the gun. We conclude that for the design example parameters, both the intrinsic energy of the gun (due to space charge effects) and space charge effects in the undulator limit the attainment of beam compression below σ t w ≈ 70   f t S and consequently the bunching factor (Eq. (7)) would diminish the superradiant radiation at frequencies beyond 3 THz. Attainment of higher frequency superradiance in future superradiant FEL designs requires technological advance in reduction of the intrinsic energy spread of the gun or enhancement of the modulation (chirp) factor S, (that would require also a shorter drift length to the undulator) and mitigation of space charge effects in the undulator region, possibly by optimizing the dimensions and distribution shape of the beam [14].