Widely tunable electron bunch trains for the generation of high-power narrowband 1–10 THz radiation

Terahertz radiation enables the coherent excitation of many fundamental modes, such as molecular rotation, lattice vibration and spin precession. To excite the transient state of matter far out of equilibrium, high-power and tunable narrowband terahertz radiation sources have been in great demand for years. However, the terahertz source available at present cannot meet these tunability and power demands, leaving a large scientific gap that has yet to be fully explored. Here we convert the energy modulation induced by the nonlinear longitudinal space charge force to a density modulation and experimentally demonstrate the generation of electron bunch trains with modulation frequencies that are adjustable between 1 and 10 THz. The electron bunch trains can be directly used to produce tunable high-power narrowband terahertz radiation that fully covers the long-standing ‘terahertz gap’, which will open up many more possibilities in terahertz science. A tunable terahertz radiation pulse is demonstrated based on a linear accelerator. The emission frequency of this terahertz radiation is tunable from 1 to 10 THz by changing the bunching frequency of a 34 MeV electron beam. The pulse energy is at the submillijoule level.

Terahertz radiation enables the coherent excitation of many fundamental modes, such as molecular rotation, lattice vibration and spin precession. To excite the transient state of matter far out of equilibrium, high-power and tunable narrowband terahertz radiation sources have been in great demand for years. However, the terahertz source available at present cannot meet these tunability and power demands, leaving a large scientific gap that has yet to be fully explored. Here we convert the energy modulation induced by the nonlinear longitudinal space charge force to a density modulation and experimentally demonstrate the generation of electron bunch trains with modulation frequencies that are adjustable between 1 and 10 THz. The electron bunch trains can be directly used to produce tunable high-power narrowband terahertz radiation that fully covers the long-standing 'terahertz gap', which will open up many more possibilities in terahertz science.
The search for high-intensity and tunable narrowband terahertz sources is motivated by broad scientific and industrial applications, such as the transient state control of materials 1-8 , imaging for security or medical applications 9 and time-domain spectroscopy 10 . State-of-the-art laser-based terahertz sources have limitations in the frequency tunability [11][12][13][14][15][16] or pulse energy 17 . Although submillijoule narrowband terahertz generation has been reported at fixed frequencies 16 , the pulse energy of tunable terahertz sources remains at the microjoule level 17 . Electron-accelerator-based terahertz sources stand out due to their exceptionally high peak or average power and wide tunability 18,19 . However, the peak power rapidly decreases as the radiation frequency increases due to the deteriorating longitudinal coherence [20][21][22] . Another approach commonly utilized in the present terahertz free-electron lasers (FELs) is the use of an optical resonator to accumulate the radiation energy extracted in the low-gain regime, although this does require a high repetition rate and multiple separate beamlines to cover the entire 'terahertz gap' [23][24][25][26] , which leads to a sophisticated design and a high maintenance cost.
Recently, terahertz electron bunch trains that consist of multiple equally spaced microbunches with picosecond and subpicosecond intervals have been widely studied to generate high-peak-power and narrowband terahertz radiation. Typical methods used to generate terahertz electron bunch trains include maintaining the initial density modulation imposed on the beam at the photocathode 27,28 , exchanging the modulated transverse profile of the beam with the longitudinal profile 29,30 , and transforming the energy modulation produced by the laser-beam interaction [31][32][33] , self-excited wakefield [34][35][36] or space charge force 37-39 into a density modulation. However, all these methods suffer from space charge effects, a nonlinear correlated energy spread or a high demand on the experimental setup, making obtaining of a high bunching factor and coverage of frequencies higher than 3 THz difficult. So far, the experimental generation of widely tunable terahertz electron bunch trains has remained a great challenge.
In this Article we propose a simple but effective mechanism to produce widely tunable electron bunch trains by transforming the sawtooth-shaped energy modulation that is produced by the nonlinear majority of the electrons reside between the spikes and form a large current pedestal, restricting the bunching factor. Moreover, when a dispersive device is applied to tune the spike intervals, the current spikes tend to be smeared out due to the relatively large energy spread, resulting in a limited frequency-tuning range.
To increase the bunching factor and broaden the frequency-tuning range, we propose to exploit the NLSC force of the modulated beam that has experienced one half-period of NLSCO to impose quasilinear negative energy chirps (the beam tail has a higher energy than the head) between the NLSCO-induced current spikes. The protruding density spikes provide a net space charge force exerted on the beam itself (see Fig. 1a). To ensure a relatively constant NLSC force during beam transport, we accelerate the beam to tens of megaelectronvolts to freeze the longitudinal distribution.
In contrast to classical sinusoidal energy modulation, which involves only the fundamental harmonic of the external field, the NLSC force of the current spikes will impose a sawtooth-shaped energy modulation on the beam due to the coherent interference of all harmonic orders (see equation (2)). Here, negative energy chirps will appear between the current spikes with positive chirps at the spikes. After passing through a magnetic chicane, the majority of the electrons (between the NLSCO-induced current spikes), which previously contributed to degradation of the bunching factor, will pile up at the same longitudinal position and form microbunches with an ultrashort bunch length due to the linearity of the sawtooth-shaped energy modulation.
Since the bunch train is obtained by converting the energy modulation, imposing an additional correlated energy chirp with a linac (linear particle accelerator) before the chicane will stretch or compress the longitudinal space charge (NLSC) force into a sharp density modulation, and we experimentally demonstrate that the frequency-tuning range can cover the entire terahertz gap (that is, 1-10 THz). The electron beam with a sharp density modulation can then be used to produce submillijoule-level terahertz sources that are tunable from 1 to 10 THz.
An electron beam with a strong initial longitudinal density modulation will experience NLSC oscillation (NLSCO). Here, we suppose that the beam has an initial longitudinal density modulation n(z,0) = n 0 (1 + 2b 0 cos(kz)), where b 0 = | | | ∫n(z,0)e ikz dz ∫n(z,0)dz | | | is the bunching factor at the wavenumber k, n 0 is the electron density and z is the coordinate along the bunch. The longitudinal density distribution (n(z,0)) and the longitudinal space charge force (F z (z,0)), after time t with phase advance ω p t = π, can be written as 38 where m e is the electron rest mass, ω p = √ e 2 n 0 ε 0 m e γ 3 is the plasma frequency, in which e is the electron charge, ε 0 is the permittivity of free space and γ is the Lorentz factor, and J m (x) is the mth-order Bessel function of the first kind. Figure 1a illustrates the longitudinal density distribution and the longitudinal component of the space charge force after the beam has undergone one half-period of NLSCO. Periodic current spikes appear in the beam longitudinal distribution. However, the projected density profile measured by the deflecting system and the magnetic spectrometer with the chicane off (c) and on (d). e, Longitudinal density distribution and projected density profile of the bunch train measured by the deflecting system with the magnetic spectrometer off. The linac phase was set on-crest and the measured modulation frequency was 2.8 THz.

Results
Our experiments were carried out at the Accelerator Laboratory of Tsinghua University. The experimental setup is shown in Fig. 1b. By controlling the solenoid strength after the gun, the electron beam experienced one half-period of NLSCO when entering the linac and formed periodic current spikes. The beam was then accelerated by the linac to freeze the longitudinal density distribution. Under the strong focusing of solenoid 2, the longitudinal space charge force induced strong energy modulation on the electron beam during a short drift after the linac. After passing through a magnetic chicane, the electrons re-grouped to form microbunches. A deflecting cavity was used to measure the longitudinal density distribution. Combined with an energy spectrometer downstream, the longitudinal phase space was measured. Figure 1c shows the measured energy modulation with a shallow density modulation when the chicane was off. After turning on the magnetic chicane to provide a longitudinal dispersion R 56 of −4 cm, the energy modulation was transformed into a density modulation, as shown in Fig. 1d.
To obtain higher temporal resolution, the spectrometer was closed off and the beam was monitored via screen 2. The measured longitudinal distribution of the terahertz electron bunch train is shown in Fig. 1e. The Fourier transform of the longitudinal density distribution gives the bunching factor of the bunch train. To increase the bunching factor, the modulation depth of the drive laser and the focusing strength of the solenoids were optimized. Figure 2a shows the simulated and measured optimized bunching factor data for different laser pulse widths. A shorter laser pulse width increases the proportion of electrons between the NLSCO-induced current spikes, and thus increases the peak current of the newly formed microbunches. However, further decreasing the pulse width of the laser will make the spike shape deviate from a Gaussian distribution and will affect the distribution of the longitudinal space charge field, which lowers the bunching and is observed in the simulation. However, for our current laser system, the shortest laser pulse width on the photocathode is approximately 300 ± 17 fs (FWHM (full-width at half-maximum)), which has not reached the optimized value for the highest bunching factor.
The linac phase was set off-crest to impose a positive or negative energy correlation along the beam. By fine-tuning the chicane strength, the bunching factor was optimized while the beam length and the microbunch intervals were stretched or compressed. The modulation frequency is tuned from 1 to 7 THz when the linac phase is varied from 20° to −35° off the crest. The measured bunching factors are shown in Fig. 2b. When the linac phase is further decreased to compress the bunch train to frequencies higher than 7 THz, the nonlinearity of the energy correlation along the entire beam becomes important, limiting the final optimized bunching factor. This nonlinearity can be linearized via an additional harmonic cavity or self-wake in a wakefield structure. Here, we only increased the chicane strength at the −35° linac phase to tune the bunch trains to 10 THz. We utilized the code ASTRA 40 to simulate the complete beam dynamics from the photocathode to the end of the beamline with space charge effects considered. The bunching factor from the simulated longitudinal distribution after the chicane and the simulated transverse distribution after the deflector are shown in Fig. 2b. The bunching factor measured by the deflector is lower due to the limited time resolution. We found in the experiments that when the beam length is compressed to tune the bunching to higher frequencies, the time resolution worsens, resulting in larger discrepancies between the measured and simulated bunching factors at frequencies higher than 4 THz. In our case, the bunching factor at 7 THz is measured to be approximately 0.15. When we further increase the frequency by compressing the beam to make it shorter, bunching cannot be clearly observed by deflector measurements. To demonstrate the microbunch formation at higher frequencies, we measured the coherent transition radiation from the microbunches, which will be shown later. Figure 3a-c illustrate the beam longitudinal distributions when the linac phase is set at 5°, −12° and −25°, respectively. The corresponding modulation frequencies are 1.9 ± 0.15 THz, 4.1 ± 0.3 THz and 6.6 ± 0.5 THz, measured by the deflecting system. The corresponding measured root-mean-squared lengths of the microbunches are approximately 50 ± 4 fs, 29 ± 2 fs and 30 ± 2 fs. Considering the time resolution, the estimated root-mean-squared length of the microbunches approaches 20 fs.
The coherent transition radiation emitted by the electron bunch trains when passing through a thin aluminium film placed after the chicane was detected using a Golay cell detector. A Michelson interferometer was installed before the Golay cell to obtain the interference signal. Figure 4a illustrates the typical results of autocorrelation measurements, measured by the deflector, for bunch trains with a bunching frequency of 2.8 ± 0.2 THz. With the direct-current offset subtracted, the Fourier transform of the autocorrelation signal gives the spectrum of the terahertz radiation, as shown in Fig. 4a. The central frequency is located at 2.84 THz, coinciding with the measurements from the deflector, and the bandwidth (FWHM) is approximately 9.3%. The beam length varies from ~16 to ~1.6 ps when the bunching frequency is tuned from 1 to 10 THz. The terahertz signal contains the contributions from both the microbunches and the current pedestals, which will be simultaneously detected by the Golay cell. Given that the bunch length is approximately a few picoseconds, the current pedestals generate radiation at frequencies mainly below 0.5 THz. To measure the narrowband terahertz radiation with only the frequency components from the microbunches, we inserted four narrow-bandpass filters with central frequencies of 1, 2.8, 5.1 and 9.8 THz before the Golay cell. The measured coherent transition radiation energy shown in Fig. 4b demonstrates the tunability of the terahertz electron bunch train from 1 to 10 THz. Table 1 summarizes the experimental parameters of the bunch trains.

Discussion
The generated terahertz electron bunch train that is tunable from 1 to 10 THz can produce high-intensity narrowband terahertz radiation when passing through an undulator. The terahertz radiation energy was simulated using the FEL code Genesis 41 . Considering a helical undulator with a period λ u = 10 cm, period number N w = 30 and undulator parameter K u varying from 0.6 to 4.2, the radiation energy produced by the terahertz electron bunch trains with a 1 nC beam charge is shown in Fig.  4c. The electron bunch trains can produce 1-10 THz tunable coherent undulator radiation with a pulse energy of hundreds of microjoules. The radiation bandwidth varies from 4 to 8% due to the different beam energy chirps.
In summary, we have demonstrated the experimental generation of terahertz electron bunch trains that are tunable from 1 to 10 THz with a measured bunching factor of up to 0.35 by exploiting the NLSC force. The NLSC force excited by the bunching structures is naturally phase-locked with the beam, which avoids the precise timing control necessary in laser-beam interaction schemes. The experimental setup is compact with no requirements for additional customized designs or complicated structures. This mechanism can also be directly applied to currently operating high-repetition-frequency accelerators without modifications to produce narrowband terahertz radiation with a high average power. Our scheme for high-intensity narrowband terahertz radiation generation that is tunable from 1 to 10 THz marks a great advance in the development of narrowband terahertz sources and provides a powerful tool for broad scientific and industrial applications.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41566-022-01131-7.

Experimental beamline
A 266-nm-wavelength ultraviolet laser from a 30 TW laser system drives the photoinjector to produce a premodulated electron beam. A 1.6-cell S-band radiofrequency gun accelerates the electron beam to approximately 5 MeV before injection into a 39-cell S-band linac with a peak gradient of 25 MV m −1 . Two solenoids are installed in the beamline to focus the beam, with one at the gun exit and the other along the linac. The electron beam is accelerated to 34 MeV after exiting the linac. Two quadrupole doublets are installed after the linac to control the beam size when the beam passes through a 1.4-m-long magnetic chicane. The R 56 of the chicane can be tuned from 0 to −15 cm. A 7-cell S-band deflecting cavity combined with four quadrupoles is installed after the chicane to measure the beam longitudinal distribution. The beam energy can be measured using a magnetic spectrometer installed at the end of the beamline.

Photocathode drive laser generation and measurements
The scheme for generating drive laser pulse trains is shown in Extended Data Fig. 1. An 800 nm wavelength laser is frequency doubled using a 0.2-mm-thick β-barium borate (β-BBO) crystal, followed by further mixing of the second harmonic with the residual fundamental in a second 0.2-mm-thick β-BBO crystal to produce the frequency-tripled laser. The pulse width of the generated 266 nm laser can be changed by adjusting the pulse width of the 800 nm laser or the time delay between the fundamental and the second harmonic. Note that the pulse shape of the 266 nm laser can vary for different time delays between the fundamental and the second harmonic, which will affect the distribution of the longitudinal space charge field of the emitted electron beam. The 266 nm laser then passes through four α-barium borate (α-BBO) crystals to create a train of 16 equally spaced laser pulses with periods of 0.44 ps. After the α-BBO crystals, there is still an 18-m-long distance to the photocathode, where the laser is transported in air and passes through three 3-mm-thick fused silica windows and one 3-mm-thick fused silica lens. The UV longitudinal profile is measured via difference-frequency mixing of the 266 nm and 800 nm lasers using a 0.1-mm-thick β-BBO crystal. To correctly consider the dispersion brought in during the transportation to the photocathode, we inserted three 3-mm-thick fused silica windows, one 3-mm-thick fused silica lens and two 4-mm-thick fused silica windows on the path of the 266 nm laser towards pulse width measurements. The group delay dispersion from two 4-mm-thick fused silica windows and a 2-m-long distance from the exit of the α-BBO crystals to the laser pulse width measurement point is approximately equivalent to the group delay dispersion from the 18-m-long distance transportation in air. By adjusting the pulse width of the 800 nm laser and the time delay between the fundamental and the second harmonic, we found that the shortest 266 nm laser pulse width is approximately 300 ± 17 fs (FWHM).

Electron beam longitudinal diagnostics
The longitudinal profile of the electron beam is measured using a deflecting system after the chicane. The resolution of the deflecting system can be written as Δ = , where σ y,0 is the vertical beam size when the deflecting cavity is closed off, W is the beam energy, ω is the angular frequency of the deflecting cavity, e is the electron charge, V def is the deflecting voltage and R 34 is the element of the transport matrix that represents the relation between the vertical divergence and the vertical displacement. The focusing strengths of the four quadrupoles after the deflecting cavity are optimized to minimize σ y,0 /R 34 . According to the simulations performed using the code ELEGANT 42 , when the beam is monitored via screen 2 with the spectrometer turned off, σ y,0 /R 34 can be optimized to 3 × 10 −5 when the beam has no energy chirp. For our 34 MeV beam energy, 2,856 MHz radiofrequency cavity and 3 MV deflecting voltage, the resolution of our deflecting system is approximately 20 fs. However, when the beam is monitored via screen 3 with the spectrometer turned on, the best resolutions considering both the time and energy resolutions are 42 fs and 25 keV according to the measurements. When calibrating the length and energy for each pixel, we take the average and standard deviation of ten shots for each measurement.

Terahertz interferometry
The transition radiation of the electron bunch trains is reflected out of the vacuum cavity through a TPX (poly(4-methylpent-1-ene)) window. An off-axis parabolic mirror (gold plated, 6 inch focal length) is installed outside the vacuum cavity to collect the terahertz radiation and redirect it to a Michelson interferometer. A beamsplitter separates the radiation into two pulses, one reflected by a fixed mirror and the other reflected by a movable mirror. The two pulses recombine after being reflected and are focused onto a Golay cell detector, which can record the pulse energy. By adjusting the position of the movable mirror to change the optical path difference, the autocorrelation of the radiation pulse is obtained. With the direct-current offset subtracted, the Fourier transform of the autocorrelation data gives the power spectrum of the terahertz radiation.

Simulation methods
The particle-tracking code ASTRA is used to simulate the beam dynamics from the photocathode to the end of the beamline with the space charge effects considered. The time-dependent FEL code Genesis is utilized to simulate the undulator radiation from the electron bunch trains. The Genesis simulations consider the evolution of the beam phase space in the radiator and the space charge effects. The 6D beam distribution from the ASTRA simulations is directly imported into Genesis and is used to load the phase space.

Data availability
The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.