The demand for electromagnetic waves in the terahertz (THz) range is significant, spanning various applications, such as molecular spectroscopy, compact electron acceleration, medical imaging, and security checking.1^–4 Over the past two decades, remarkable advancements have been achieved in the development of diverse THz radiation sources. Notably, electron accelerator-based THz sources have garnered substantial global interest due to their potential to generate THz radiation with narrow bandwidth, high energy, exceptional brightness, and impressive repetition rates ranging from kilohertz to megahertz.5^–7 Despite the existence of various techniques for generating narrowband THz radiation,8^–10 the utilization of electron bunch trains holds immense potential for producing highly intense and narrowband THz radiation based on Smith–Purcell radiation,11 undulator radiation,12 etc. The crucial factor in achieving coherent THz radiation lies in the precise manipulation of the beam phase-space distribution, enabling the formation of microbunches with picosecond (ps) or even sub-picosecond spacing.13^–15 This approach offers promising prospects for advancing the field of THz generation.

The generation of the desired bunch train has been a subject of intense investigation, for example, directly modulating the drive laser,13^,16^–19 transforming transverse modulation to a longitudinal distribution,20^,21 using self-modulation instability in a linear plasma wake,22 converting energy modulation induced by wakefields in dielectric-lined or corrugated waveguides to density bunching,23^–27 and transforming laser-induced energy modulation to a density distribution.11^,28^–30 For the bunch train, the bunching factor $b$, defined as $b(k)={\int }_{-\infty }^{\infty }I(z)e^{-ikz}\mathrm{d}z/{\int }_{-\infty }^{\infty }I(z)\mathrm{d}z$ with the beam current distribution $I(z)$ and the modulation wave number $k$, is used to reflect the coherence of the longitudinal density distribution of electron beams. The coherent THz radiation energy is proportional to the square of the bunching factor $b$. Therefore, the larger $b$, the stronger the radiation energy. However, in the above methods, even under conditions of theoretically ideal simulations, the bunching factor $b$ is limited to $\sim 0.4$, which restricts the acquisition of higher radiation energy.

An alternative method of inducing energy and density modulation in an electron beam is to utilize the plasma medium. Recently, the use of uniform plasma-based modulators has been proposed to generate ultrahigh peak-current ($\sim 10\mathrm{kA}$ level) electron bunch trains with tunable ps spacing and high bunching factors.14 In the scheme, a high-charge (several nC) and high-current electron beam traverses a uniform plasma section and initiates a periodic “sawtooth” nonlinear plasma wake when the beam density is higher than the plasma density. This interaction leads to energy gain or loss for electrons at different positions, creating a sawtooth-shaped energy modulation along the bunch. Subsequently, a magnetic compressor is employed to efficiently convert this energy modulation into beam density modulation, resulting in the formation of microbunches with an impressive bunching factor of $\sim 0.8$.14 Although the scheme offers clear advantages, to excite a high-quality nonlinear plasma wakefield, the electron beam should be focused to a very small transverse size [$\sim 5\mu \mathrm{m}$ in root mean square (RMS)] and its current distribution is preferred to be a flattop distribution, which is hard to realize in practical experiments.

Here, we propose a novel two-bunch scheme for the generation of high-power coherent THz radiation. In this scheme, an additional short-pulse drive beam is added before the main beam as the excitation source of nonlinear plasma wake. This innovative approach can significantly loosen the requirements on main beam parameters (transverse size, longitudinal current distribution, etc.) and improve the bunching factor to $\sim 0.94$, thereby optimizing the energy conversion efficiency of THz generation. As passive modulation, compared with dielectric-lined or corrugated waveguide-based energy modulators, using plasma wakefield to implement energy modulation has unique advantages in tuning modulation frequency, avoiding radio-frequency (RF) breakdown, providing an ideal sawtooth energy modulation, etc. The nonlinear plasma wake has very good wakefield properties within the ion channel, including a constant longitudinal field in the transverse dimension and a uniform focusing gradient along the longitudinal dimension. These properties effectively preserve the energy spread and emittance of the beam slice.

As shown in Fig. 1(a), two electron beams enter a uniform plasma sequentially. The first beam, referred to as the “drive beam,” is characterized by relatively low or moderate charge but high current ($\sim \mathrm{kA}$) for its short bunch length (<ps). The drive beam boosts a brief duration and exceptionally high current, and even if its transverse size may be relatively large, its density significantly surpasses that of the surrounding plasma, resulting in an ideal sawtooth-shaped wakefield. The second beam, referred to as the “main beam,” features a relatively long bunch length ($\sim \mathrm{ps}$) and low current. Following the drive beam in plasma that excites a strong nonlinear wakefield, the main beam would not only experience the linear wakefield produced by itself but also gain modulation from the nonlinear wakefield generated by the former. Given that the current/density of the drive beam is much higher than that of the main beam, the nonlinear wakefield excited by the drive beam dominates the energy modulation process while the linear wakefield self-excited by the main beam plays a little role; therefore, there are no stringent demands on the transverse size or current profile of the main beam. Then, a small magnetic chicane located next to the plasma section is used to detour the electron bunch, and the sawtooth energy modulation of the main beam is converted to a density distribution with a very high bunching factor. Since the two beams have slight differences in energy, their different motion paths and deflection angles in the chicane could be applied to block the drive beam and have the main beam enter the subsequent section for radiation by introducing a metal (such as tungsten) slit as an energy filter.

To illustrate the key physics in the above scheme, a detailed theoretical analysis is carried out. For a uniform plasma section with density $n_p$ and a relativistic electron beam with density $n_b$, length $L_b$, and RMS spot size ${\sigma }_r$ (transverse Gaussian profile), nonlinear plasma wakefields are excited and plasma bubbles are formed when $n_b>n_p$.31 The shape of the excited bubble is represented by the trajectory of the innermost particle. The longitudinal wakefield $E_z$ is found to be proportional to the product of the local radius of the ion channel $r_b$ and the slope $d{r}_b/d\xi$,32$$E_z(\xi )\approx \frac{1}{2}{r}_b\frac{d{r}_b}{d\xi },$$(1)where $\xi =ct-z$ with the light speed in vacuum $c$ and the electron beam moving direction $+z$.

The normalized maximum radius of the ion channel is $k_p{r}_m\approx 2.58\sqrt{\mathrm{\Lambda }}$,14 in which $k_p=\sqrt{n_p{e}^2/m{\epsilon }_0{c}^2}$ is the plasma wavenumber, $\mathrm{\Lambda }={\int }_0^{\infty }{k}_{p}r(\frac{n_b}{n_p})\mathrm{d}{k}_{p}r=(\frac{n_b}{n_p})k_p^2{\sigma }_r^2=2I_b/I_A$ is the normalized charge per unit length, $e$ is the electron charge, $m$ is the rest mass of the electron, ${\mathit{\epsilon }}_0$ is the vacuum dielectric constant, $I_b=Q_c/L_b$ is the peak current of the beam, and $I_A\approx 17\mathrm{kA}$ is the Alfven current.

In this case, the trajectory of the inner particles of the ion channel can be simplified as32$$(1+\frac{1}{4}{r}_b^2)\frac{d^2{r}_b}{\mathrm{d}{\xi }^2}+\frac{1}{2}{r}_b{\left(\frac{d{r}_b}{d\xi }\right)}^2+\frac{1}{2}{r}_b=\frac{\lambda (\xi )}{r_b},$$(2)where the drive term $\lambda (\xi )$ in the beam region is $\mathrm{\Lambda }$ and outside the beam region is 0. Normalized units are applied here with the lengths normalized to skin depth $k_p^{-1}$, densities to the plasma density $n_p$, charges to the electron charge $e$, and fields to the cold nonrelativistic wave-breaking limit $E_p=m{k}_p{c}^2/\mathrm{e}$.

As the beam enters into the plasma section, a series of bubbles are excited along the propagating direction. In each bubble, we have $d{r}_b/d\xi =0$ at the maximum radius position $(\xi ,r_b)=({\xi }_{m,i},r_m)$, in which $i$ is the bubble number, and $m$ indicates the maximum radius position. Applying Taylor expansion at $\xi ={\xi }_{m,i}$, we have $$r_b(\xi )\approx r_m+\frac{1}{2}\frac{d^2{r}_b}{\mathrm{d}{\xi }^2}{\xi }_{m,i}(\xi -{\xi }_{m,i}{)}^2,$$(3)and thus obtain $$\frac{d{r}_b}{d\xi }\approx \frac{d^2{r}_b}{\mathrm{d}{\xi }^2}{\xi }_{m,i}(\xi -{\xi }_{m,i}),$$(4)and $$E_z(\xi )\approx \frac{1}{2}{r}_b\frac{d{r}_b}{d\xi }\approx \frac{1}{2}{r}_m\frac{d^2{r}_b}{\mathrm{d}{\xi }^2}{\xi }_{m,i}(\xi -{\xi }_{m,i}).$$(5)

From Eq. (2), we have $$\frac{d^2{r}_b}{\mathrm{d}{\xi }^2}({\xi }_{m,i})=-\frac{\frac{1}{2}-\lambda ({\xi }_{m,i})/r_m^2}{1+\frac{1}{4}{r}_m^2}{r}_m,$$(6)and thus the slope of $E_z$ at $\xi ={\xi }_{m,i}$ is $$\frac{d{E}_z}{d\xi }({\xi }_{m,i})=-\frac{1}{2}{r}_m^2\frac{\frac{1}{2}-\lambda ({\xi }_{m,i})/r_m^2}{1+\frac{1}{4}{r}_m^2}.$$(7)

Since the main beam is applied to generate THz radiation, the analysis focuses on its energy modulation. When the two beams enter the uniform plasma section, the drive beam moving ahead exists mainly at the first bubble, while the main beam exists at the second, third, fourth, and the following bubbles. Therefore, we focus on the field distributions of the second and the following bubbles. Since the beam current of the drive beam is much higher than that of the main beam, the radius of the first bubble $r_{m,1}$ is the largest, and $r_{m,i}(i\ge 2)$ is relatively small, namely, $r_{m,1}>r_{m,i}(i\ge 2)$. Assuming $r_{m,2}\approx r_{m,3}\approx r_{m,4}=\alpha ·r_{m,1}$ ($0<\alpha <1$), the value of $\alpha$ could be obtained through numerical simulations. Since the main beam for generating radiation exists in the second and the following bubbles, we focus on the field distributions in these bubbles. Taking the second bubble as an example, Eq. (7) is still applicable; that is, the slope of the longitudinal electric field is $$\frac{d{E}_z}{d\xi }({\xi }_{m,2})=-\frac{1}{2}{r}_{m,2}^2\frac{\frac{1}{2}-\lambda ({\xi }_{m,2})/r_{m,2}^2}{1+\frac{1}{4}{r}_{m,2}^2},$$(8)where $\lambda ({\xi }_{m,2})=2I_{\text{witness}}/I_{\mathrm{A}}$, and $r_{m,2}\approx \alpha ·r_{m,1}=\alpha ·2.58\sqrt{\frac{2I_{\text{drive}}}{I_{\mathrm{A}}}}$.

To further demonstrate the practical validity of the above scheme, we show here a typical example through self-consistent three-dimensional (3D) simulations utilizing a set of thoroughly benchmarked codes. The propagation of the two beams within the uniform plasma is simulated using the particle-in-cell code QuickPIC.33^–35 Subsequently, the macroparticles of the beam are imported into the general particle tracer (GPT) code,36 enabling comprehensive tracking of the beams’ six-dimensional (6D) phase-space dynamics in the chicane.

In our simulations, the drive beam is characterized by a beam current of 2 kA, beam charge of 1 nC, and flattop current profile of 0.5 ps ($L_b=0.15\mathrm{mm}$). The main beam with the energy of 135 MeV has a beam current of $I_{\text{witness}}=0.3\mathrm{kA}$, beam charge 1.5 nC, and 5 ps flattop current profile ($L_b=1.5\mathrm{mm}$). The simulation parameters are shown in Table 1 in detail. Two such electron beams with the above characteristics and a suitable time delay can be generated using typical photoinjectors together with the following injector linac sections. 18^,37^,38

The plasma density distribution after the two beams propagating into the plasma is shown in Fig. 1(b). Since the current of the drive beam is much higher than that of the main beam, the radius of the first bubble $r_{m,1}$ is the largest, and $r_{m,N}$$(N\ge 2)$ are relatively small. The excited longitudinal wakefield $E_z$ in $x\text{-}\xi$ plane and the on-axis lineout of $E_z$ clearly show the feature of a sawtooth waveform in Fig. 1(c). With the drive beam current of $I_{\text{drive}}=2\mathrm{kA}$, we can have $r_{m,1}=2.58\sqrt{\frac{2I_{\text{drive}}}{I_{\mathrm{A}}}}\approx 1.24$. In the second bubble, $\lambda ({\xi }_{m,2})=\frac{2I_{\text{witness}}}{I_{\mathrm{A}}}\approx 0.035$, $\alpha \approx 0.8$ (that is $r_{m,2}\approx 1$), thus the slope of $E_z$ (normalized to $m_e{\omega }_p^2/e$) is $-0.19$. In the simulation, as in Fig. 1(c), $r_{m,1}\approx 1.23$ and the slope of $E_z$ is $-0.20$, which are both in good agreement with the theoretical results. A lead collimator is placed between the plasma and the chicane, leading to a prebunching effect. The electrons in the region of positive $E_z$ slope have a relatively large transverse momentum $p_x$, which can be removed during further transport with the prebunching effect. Such an effect will further reduce the DC component after the chicane and increase the bunching factor. By assuming that electrons with $|p_x|>1\mathrm{mc}$ can be removed, the charge-loss ratio is relatively low ($\sim 30\%$).

We then use the particle tracking code GPT to further track the beam 6D phase-space dynamics in the chicane. In the simulation, the $R_{56}$ represents the change in the average beam position along the beamline per unit change in the particle’s momentum,39 which is applied to quantify the dispersion effect in the chicane. Dependence of longitudinal position on momentum is mainly caused by transverse dispersion ($\eta$) in bending magnets, which results in different path lengths for particles with different momenta. The simulated $R_{56}$ value of “full compression” is $-0.96\mathrm{mm}$; the beam longitudinal phase space is shown in Fig. 2(a). The beam energy modulation is converted into density modulation, resulting in the microbunches with an ideal bunching factor. Each microbunch stands nearly upright, which leads to a high peak current. The corresponding current profile shown in Fig. 2(b) reveals that the bunch length of each microbunch reduces to $\sim 50\mathrm{fs}$ (RMS) and the peak current reaches $\sim 4.5\mathrm{kA}$. Thus, the corresponding bunching factor $b$ at the fundamental and harmonic frequencies could be calculated and are shown in Fig. 2(c). The $b$ value reaches 0.94 at the fundamental frequency and over 0.5 at high harmonic frequencies, which is larger than those of previous methods with no stringent demands.

Sending the bunch train through a helical undulator, the radiation at the THz frequency region could be emitted. With undulator period ${\lambda }_u=20\mathrm{cm}$ and the peak magnetic field $B_u=1.1\mathrm{T}$, the undulator strength $K=0.934{\lambda }_u{B}_u=20.5$, and the radiation can be resonant at 1.13 THz (corresponding radiation wavelength ${\lambda }_r=264\mu \mathrm{m}$ with the plasma density $n_p=1.6\times {10}^{16}{\mathrm{cm}}^{-3}$).

As a longitudinally dispersive element, the undulator causes overcompression or debunching of a fully compressed electron beam and leads to a slight reduction in the bunching factor and thus the output THz radiation. To obtain a high average bunching factor within the undulator, the beam is undercompressed in the chicane to obtain a bunching factor of 0.85 at the entrance of the undulator. As the beam propagates in the undulator, it will be first further compressed to the full-compression and then overcompressed case. The simulation of the bunch trains passing through the undulator is carried out using the well-benchmarked code GENESIS,40^,41 and the space charge effect is also included. As shown in Fig. 3, simulation results show that with a 5 m-long undulator, the emitted THz pulse energy can reach as high as 2.4 mJ. The radiation is also highly directional, which is of great advantage in advanced scientific explorations. When the beam charge increases, the energy increases and reaches a high level accordingly, since the bunching factor is above $\sim 0.9$ when loosening the requirements on the two bunches.

Lowering the energy of the electron beam is beneficial for improving the energy conversion efficiency. Therefore, we further explored the feasibility of the scheme at lower electron energies. Reducing the energy of the dual bunch in the previous section to 50 MeV with other parameters unchanged, the simulated results of the beam are shown in Fig. 4. With the same plasma length (3.7 mm), the absolute value of energy modulation received by the electron beam is the same, but the relative energy spread is obviously increased [see Fig. 4(a); compared with that in Fig. 2(a)] due to the decrease of the initial energy. In Fig. 2(a), the longitudinal phase space is “rigid,” while in Fig. 4(a) there is an obvious curvature, which results from the nonlinear process of magnetic compression. When the relative energy spread increases, the nonlinearity of the beam longitudinal phase space becomes apparent. In this case, the bunching factor is 0.85 at the fundamental frequency.

In order to reduce the influence of relative energy spread, we decrease the length of the uniform plasma. Figure 5(a) shows the beam phase space with a plasma length decreased to 2.7 mm. Compared with Fig. 4(a), the relative energy spread decreases, and the longitudinal phase space is more rigid. Figures 5(b) and 5(c) show the current distribution and bunching factor $b$ under different plasma lengths, with the red line indicating that the highest bunching factor $b$ reaches 0.9. Through simulations, radiation energy in the same order can still be obtained. Therefore, by reducing the length of the plasma, the energy spread of the beam can be reduced and the bunching factor $b$ could be improved through the magnetic compression process. With reduced plasma length and energy spread, the radiation energy could be further increased by 10%.

In conclusion, we propose and study an approach leveraging uniform plasma to sequentially modulate two bunches, aimed at enhancing THz radiation generation. Through theoretical analysis and comprehensive simulations, it is revealed that, by adding an additional short-pulse drive beam before the main beam as the excitation source of nonlinear plasma wake, the bunching factor of the main beam can be further increased to 0.94 even though with a relatively low charge, low current, and large-size electron beam. Consequently, significant strides have been made in boosting THz radiation energy while mitigating constraints, such as beam charge and transverse beam size. The findings of this research not only contribute to the scientific understanding of THz radiation generation but also pave the way for potential breakthroughs in the practical utilization of this valuable radiation.

Biographies of the authors are not available.