2.1 Theoretical model and simulation setup

When an intense laser pulse irradiates on a plasma, as shown in Figure 1(b), hot electrons whose beam length is close to the laser duration can be generated by laser ponderomotive force^[ Reference Wilks, Kruer, Tabak and Langdon ^{39 ]}. The electron beam transports in the plasma and passes through its transverse interfaces. As a result, transverse sheath fields ${E}_{\mathrm{s}}$ (more details in the Supplementary Material) can be generated and the maximum strength^[ Reference Daido, Nishiuchi and Pirozhkov ^{40 ]} on the plasma surfaces can be expressed as follows:

(1) $$\begin{align}{E}_{\mathrm{s}}=\sqrt{\frac{8\pi {n}_{\mathrm{e}}{T}_{\mathrm{e}}}{e}},\end{align}$$

where $e$ is the Euler’s number and ${n}_{\mathrm{e}}$ and ${T}_{\mathrm{e}}$ are the density and temperature of the hot electrons, respectively. Here, transverse sheath fields ${E}_{\mathrm{s}}$ higher than ${10}^{12}$ V/m can be induced near the plasma transverse interfaces^[ Reference Tarkeshian, Vay, Lehe, Schroeder, Esarey, Feurer and Leemans ^{41 ]}, as shown in Figure 1(c). Under the action of ${E}_{\mathrm{s}}$ , electrons whose kinetic energy is less than ${\varepsilon}_{\mathrm{e}}$ could be pulled back into the plasma and pass through the transverse interface on the other side^[ Reference Sentoku, Cowan, Kemp and Ruhl ^{42 ]}, where ${\varepsilon}_{\mathrm{e}}$ could be calculated by the following:

(2) $$\begin{align}{\varepsilon}_{\mathrm{e}}={\int}_0^{l_{\mathrm{e}}}\frac{E_{\mathrm{s}}}{\sin \phi } \mathrm{d}l.\end{align}$$

Here, $l$ is the transverse size of ${E}_{\mathrm{s}}$ , $\phi$ is the angle at which the electron enters the transverse sheath fields ${E}_{\mathrm{s}}$ and ${l}_{\mathrm{e}}$ is the transverse location where the electron could be pulled back into the plasma. An electron field derived from the simulation was used (more details in the Supplementary Material) for calculation, and then the relation between ${\varepsilon}_{\mathrm{e}}$ and ${l}_{\mathrm{e}}$ could be plotted in Figure 2 by employing ${E}_{\mathrm{s}}$ as shown in Figure 1(c). In the case that the electron perpendicularly penetrates ${E}_{\mathrm{s}}$ ( $\phi ={90}^{\circ }$ ), the threshold kinetic energy is about 6 MeV. Most of the electrons penetrate the transverse sheath fields with a much smaller $\phi$ , and the threshold kinetic energy could be much higher than 6 MeV. Some electrons, which are located in the forefront of the electron beam and experienced much weaker ${E}_{\mathrm{s}}$ , could not be pulled back to the plasma even with a much smaller kinetic energy^[ Reference Mora ^{43 ]}. According to the large difference in the dynamical behavior of the hot electrons, as shown in Figure 1(c) with red and blue colors, we could separate the hot electrons into two groups: the electrons in group A located in the front of the electron beam are marked with a red color and the electrons in group B located in the beam rear express are marked with a blue color.

Figure 2 In the case of the electron penetrating ${E}_{\mathrm{s}}$ with different $\phi$ (the angle at which the electron enters ${E}_{\mathrm{s}}$ ), the relation between the electron threshold kinetic energy ${\varepsilon}_{\mathrm{e}}$ and the transverse location ${l}_{\mathrm{e}}$ where the electron could be pulled back into the plasma.

Here, ${E}_{\mathrm{s}}$ , which propagates along the plasma surface^[ Reference Li, Yuan, Xu, Zheng, Sheng, Chen, Ma, Liang, Yu, Zhang, Liu, Wang, Wei, Zhao, Jin and Zhang ^{44 ]} at a velocity close to light speed, can guide the hot electrons in group B to oscillate transversely and propagate longitudinally, as shown in Figure 1(c). Every time the electron passes through the transverse interfaces of the plasma, THz radiation can be emitted. This behaves like a wiggler in traditional light sources^[ Reference Bilderback, Elleaume and Weckert ^{45 ]}. Therefore, the plasma with transverse sheath fields ${E}_{\mathrm{s}}$ can be considered as a wiggler for the electrons in group B. The reciprocating (wiggler) period of the electron can be expressed as follows:

(3) $$\begin{align}{\tau}_{\mathrm{r}}=\frac{4{l}_{\mathrm{s}}+2{l}_{\mathrm{t}}}{v_{\mathrm{t}}}.\end{align}$$

Here, ${l}_{\mathrm{s}}$ is the transverse distance between the plasma interface and the point where the electron turns around, ${l}_{\mathrm{t}}$ is the thickness of the plasma and ${v}_{\mathrm{t}}$ is the median transverse velocity of the electrons; ${l}_{\mathrm{s}}$ is determined by the transverse momentum of electrons and the strength of ${E}_{\mathrm{s}}$ . The radiation frequency $f$ is closely related to the wiggler frequency ${f}_{\mathrm{w}}=1/{\tau}_{\mathrm{r}}$ , and then we have the following:

(4) $$\begin{align}f=\frac{f_{\mathrm{w}}\left(2{\gamma}^2\right)}{1+{K}^2/2}=\frac{v_{\mathrm{t}}\left(2{\gamma}^2\right)}{\left(4{l}_{\mathrm{s}}+2{l}_{\mathrm{t}}\right)\left(1+{K}^2/2\right)},\end{align}$$

where $\gamma$ is the electron Lorentz factor and $K$ is the wiggle (undulator) parameter^[ Reference Jackson ^{46 ]}. Normally, ${l}_{\mathrm{s}}$ is several micro-meters^[ Reference Mora ^{43 ]}; while the changes in $\gamma$ , $K$ and ${v}_{\mathrm{t}}$ can be ignored for the same laser intensity and plasma density, then ${l}_{\mathrm{t}}$ is the only parameter that could be varied. Hence, one can regulate the frequency of the THz radiation by changing the plasma thickness ${l}_{\mathrm{t}}$ according to Equation (4). Since the length of the electron beam is close to the laser pulse duration, which is only 32 fs here, the radiation satisfies the time coherence condition at the THz frequencies^[ Reference Liao, Li, Liu, Scott, Neely, Zhang, Zhu, Zhang, Armstrong, Zemaityte, Bradford, Huggard, Rusby, McKenna, Brenner, Woolsey, Wang, Sheng and Zhang ^{32 ]}.

2.2 Electron kinetics

To generate electrons in group B, one should use a plasma with limited transverse size to ensure the occurrence of reciprocating motion. Simulations, with plasma of different lengths ${l}_{\mathrm{p}}$ (fixing the plasma thickness ${l}_{\mathrm{t}}$ to $30\;\mu \mathrm{m}$ ), were performed to see the effect of plasma length on the electron dynamics. From the simulations with ${l}_{\mathrm{p}}\ge 50\;\mu \mathrm{m}$ , it was found that more than $75\%$ of the laser energy was converted to hot electrons (electron kinetic energy ${e}_{\mathrm{k}}\ge 0.5$ MeV). Meanwhile, in the case of ${l}_{\mathrm{p}}=10\;\mu \mathrm{m}$ , there were few electrons participating in reciprocating motion and most of the electrons could be classified into group A.

When ${l}_{\mathrm{p}}$ was increased to $50\;\mu \mathrm{m}$ , one could easily distinguish the electrons in group B from those in group A, as plotted in Figures 3(a) and 3(b). Figure 3(a) shows a typical angular-spectra distribution of the electrons, with a large divergence angle, accelerated by the laser ponderomotive force^[ Reference Wilks, Kruer, Tabak and Langdon ^{39 ]}. This figure is very similar to the result of ${{l}_{\mathrm{p}}=10\;\mu \mathrm{m}}$ . Under the action of ${E}_{\mathrm{s}}$ , the reciprocating motion of the electrons is facilitated and obvious changes in the angular-spectra distribution can be seen in Figure 3(b). In this case, the electron beam had a divergence angle of approximately ${60}^{\circ }$ , which indicates that most of the electrons penetrate the transverse sheath fields with an angle of ${\phi \approx {60}^{\circ }}$ . Hence, the threshold kinetic energy of the electron ${\varepsilon}_{\mathrm{e}}$ is about 8 MeV, which is consistent with the theoretical value predicted by Equation (2), as shown in Figure 2. If ${l}_{\mathrm{p}}$ is further increased, such as ${l}_{\mathrm{p}}=200\;\mu \mathrm{m}$ , more electrons of relatively low energy in group A come into group B, as shown in Figure 3(c), and the number of electrons in group B increases with ${l}_{\mathrm{p}}$ . As a result, the electron number shown in Figure 3(d) is larger than that in Figure 3(b). It is noted that the electrons in group B were manipulated and collimated into a smaller divergence angle by ${E}_{\mathrm{s}}$ , and the electron threshold kinetic energy ${\varepsilon}_{\mathrm{e}}$ could be much higher, as predicted by Equation (2).

Figure 3 The angular-spectra distribution of the hot electrons. The electrons from a plasma length of ${l}_{\mathrm{p}}=50\;\mu \mathrm{m}$ (a) collected by a screen with a radius of $55\;\mu \mathrm{m}$ in the first 270 fs of the simulation could be classified to group A, and (b) the electrons behind group A can be assigned to group B. The electrons from a plasma length of ${l}_{\mathrm{p}}=200\;\mu \mathrm{m}$ (c) in group A collected by a screen whose radius is $205\;\mu \mathrm{m}$ in the first 770 fs and (d) the electrons in group B.

We used a semicircular receiving screen, whose radius was 250 μm with the center locating at the midpoint of the plasma right-hand interface, to collect the radiation field. Then, the angular-spectra distribution of the THz source could be obtained from the radiation field through Fourier transform^[ Reference Cai, Shou, Han, Huang, Wang, Song, Geng, Yu and Yan ^{47 ]}. In this work, the angular spectrum method^[ Reference Cai, Shou, Han, Huang, Wang, Song, Geng, Yu and Yan ^{47 –} Reference Ritter ^{49 ]} was employed to remove the near-field radiation, and then all the results of the THz source could be recognized as far-field radiation.

In the simulation with ${l}_{\mathrm{p}}=150\;\mu \mathrm{m}$ , the median longitudinal (transverse) velocity ${v}_{\mathrm{l}}$ ( ${v}_{\mathrm{t}}$ ) of the electrons in group B was 0.94c (0.85c). The transverse size of ${E}_{\mathrm{s}}$ was about $2\;\mu \mathrm{m}$ (at full width at half maximum (FWHM)), and therefore we could assume ${l}_{\mathrm{s}}\approx 2\;\mu \mathrm{m}$ . According to Equation (3), one can get ${\tau}_{\mathrm{r}}=270$ fs and ${\lambda}_{\mathrm{w}}\approx {\tau}_{\mathrm{r}}\times {v}_{\mathrm{l}}=79.4\;\mu \mathrm{m}$ . The wiggler period ${\tau}_{\mathrm{w}}={\tau}_{\mathrm{r}}$ corresponds to a wiggle frequency ${f}_{\mathrm{w}}=3.7$ THz, and the Lorentz factor $\gamma$ for the longitudinal velocity ${v}_{\mathrm{l}}$ was 2.93. By tracking the electromagnetic fields exerting on the electrons in group B, we got an average value of $K\approx 6.0$ , which decreased to 5.5 when ${l}_{\mathrm{s}}$ increased to $80\;\mu \mathrm{m}$ . According to Equation (4), the radiation frequency was ${f=3.34}$  THz, which is also the same as the center frequency of 3.35 THz from the simulation. Hence, the theoretical model can accurately predict the terahertz frequency obtained from the simulations.

2.3 Spectrum of terahertz pulses

Figures 4(a)–(c) show the angular-spectra distribution of the THz pulses from the plasma lengths ${l}_{\mathrm{p}}=50$ , 150 and $300\;\mu \mathrm{m}$ . In the condition of ${l}_{\mathrm{p}}$ being shorter than ${\lambda}_{\mathrm{w}}$ , the role of the electrons in group B could be ignored and the generation of the THz source is dominated by target rear CTR^[ Reference Liao, Li, Liu, Scott, Neely, Zhang, Zhu, Zhang, Armstrong, Zemaityte, Bradford, Huggard, Rusby, McKenna, Brenner, Woolsey, Wang, Sheng and Zhang ^{32 ,} Reference Liao, Liu, Scott, Zhang, Zhu, Zhang, Li, Armstrong, Zemaityte, Bradford, Rusby, Neely, Huggard, McKenna, Brenner, Woolsey, Wang, Sheng and Zhang ^{37 ,} Reference Liao, Li, Li, Su, Zheng, Liu, Wang, Hu, Yan, Dunn, Nilsen, Hunter, Liu, Wang, Chen, Ma, Lu, Jin, Kodama, Sheng and Zhang ^{50 ]} from the electrons in group A. Hence, the THz source had a large divergence angle (pointed at 88°) and a small energy conversion efficiency of 0.45%, as shown in Figure 4(a).

Figure 4 Simulation results from the plasma of different lengths ${l}_{\mathrm{p}}$ , while the thickness ${l}_{\mathrm{t}}$ was fixed to $30\;\mu \mathrm{m}$ . The angular-spectra distribution of the THz pulses from (a) ${l}_{\mathrm{p}}=50\;\mu \mathrm{m}$ , (b) ${l}_{\mathrm{p}}=150\;\mu \mathrm{m}$ and (c) ${l}_{\mathrm{p}}=300\;\mu \mathrm{m}$ . (d) The radiation field before filtering (blue line) and the field of the THz pulse (red line) collected at ${37}^{\circ }$ from the simulation of ${l}_{\mathrm{p}}=300\;\mu \mathrm{m}$ .

Every time the electron passes through the transverse interfaces of the plasma, THz radiation can be emitted and becomes stronger in the case of the electrons in group B oscillating for several periods. When ${l}_{\mathrm{p}}$ was enlarged to $150\;\mu \mathrm{m}$ , more than $1.3\%$ of the laser energy was converted to THz radiation and the characteristics of the narrow spectrum became obvious, as shown in Figure 4(b). The simulation results indicated that ${l}_{\mathrm{p}}\ge 150\;\mu \mathrm{m}$ was more conducive to obtain THz pulses of a narrow spectrum, as shown in Figures 4(b) and 4(c). As the electrons in group B are guided to oscillate along the plasma by ${E}_{\mathrm{s}}$ , the resulting THz pulse could be better collimated than the other intense laser–plasma-based THz sources^[ Reference Yi and Fülöp ^{26 ,} Reference Liao, Liu, Scott, Zhang, Zhu, Zhang, Li, Armstrong, Zemaityte, Bradford, Rusby, Neely, Huggard, McKenna, Brenner, Woolsey, Wang, Sheng and Zhang ^{37 ,} Reference Liao, Li, Li, Su, Zheng, Liu, Wang, Hu, Yan, Dunn, Nilsen, Hunter, Liu, Wang, Chen, Ma, Lu, Jin, Kodama, Sheng and Zhang ^{50 ]}. It was found that the THz pulse that pointed at $\theta ={37}^{\circ }$ could be collimated into a divergence angle of approximately ${20}^{\circ }$ , as shown in Figure 4(c). Since $\theta >1/\gamma$ , we call it a wiggler instead of an undulator. Then we collected the radiation field at ${37}^{\circ }$ from the simulation of ${l}_{\mathrm{p}}=300\;\mu \mathrm{m}$ , and filtered the THz signal (frequency range: 0.1–10 THz) from the radiation field, as shown in Figure 4(d). The THz field was higher than 80 GV/m, which was much greater than those from THz sources driven by similar laser parameters^[ Reference Yi and Fülöp ^{26 ,} Reference Liao, Liu, Scott, Zhang, Zhu, Zhang, Li, Armstrong, Zemaityte, Bradford, Rusby, Neely, Huggard, McKenna, Brenner, Woolsey, Wang, Sheng and Zhang ^{37 ]}, after the THz pulse propagated $250\;\mu \mathrm{m}$ off the plasma. Because the THz radiation comes from the wiggler process, the THz electric field waveform of multi-cycles^[ Reference Tian, Liu, Bai, Zhou, Sun, Liu, Zhao, Li and Xu ^{51 ]} was completely different from the half-cycle THz sources driven by intense lasers and plasma^[ Reference Yi and Fülöp ^{26 ,} Reference Liao, Liu, Scott, Zhang, Zhu, Zhang, Li, Armstrong, Zemaityte, Bradford, Rusby, Neely, Huggard, McKenna, Brenner, Woolsey, Wang, Sheng and Zhang ^{37 ,} Reference Liao, Li, Li, Su, Zheng, Liu, Wang, Hu, Yan, Dunn, Nilsen, Hunter, Liu, Wang, Chen, Ma, Lu, Jin, Kodama, Sheng and Zhang ^{50 ]}.

2.4 Tunable frequency

From Figures 3(a) and 3(c), one can see that the electrons in group A had a large divergence angle. The electron density ${n}_{\mathrm{e}}$ decreases with the increase of transmission distance, and ${E}_{\mathrm{s}}$ becomes weaker for the longer plasma according to Equation (1), in which case the electrons in group B are more difficult to pull back to the plasma. As a result of enhancing ${l}_{\mathrm{p}}$ , ${l}_{\mathrm{s}}$ becomes larger and more electrons no longer participate in the wiggler motion. Hence, the center frequency decreases and the growth of energy conversion efficiency gradually slows down with the enhancement of ${l}_{\mathrm{p}}$ . Simulations by varying ${l}_{\mathrm{p}}$ from 50 to $900\;\mu \mathrm{m}$ have demonstrated the effect on the center frequency, as shown in Figure 5(a). From the figure, one can see that the center frequency of the THz pulse was tunable from 3.5 to 2.6 THz by changing the plasma length ${l}_{\mathrm{p}}$ from 150 to $900\;\mu \mathrm{m}$ . On the other hand, the laser–THz energy conversion efficiency, which increased with ${l}_{\mathrm{p}}$ and became saturated after $600\;\mu \mathrm{m}$ , as shown in Figure 5(a), was 2.05%. The saturation length may be longer than $600\;\mu \mathrm{m}$ since a small portion of the THz source came out of the simulation window, in the case of longer ${l}_{\mathrm{p}}$ , due to the limitation of computational ability.

Figure 5 (a) The center frequency of the THz source and the laser–THz energy conversion efficiency from the simulations with plasmas of different lengths ${l}_{\mathrm{p}}$ from 50 to $900\;\mu \mathrm{m}$ , while ${l}_{\mathrm{t}}$ was fixed to $30\;\mu \mathrm{m}$ . (b) The center frequency of the THz source from the simulations (blue line) and the theoretical model of Equation (4) for the plasma thickness ${l}_{\mathrm{t}}$ changing from 20 to $80\;\mu \mathrm{m}$ (light black shadow).

Equation (4) indicates that we could regulate the center frequency of the THz pulse by changing the thickness of the plasma ${l}_{\mathrm{t}}$ . More simulations were performed by changing ${l}_{\mathrm{t}}$ from 20 to $80\;\mu \mathrm{m}$ (fixed ${l}_{\mathrm{p}}$ at $300\;\mu \mathrm{m}$ ). It was found that the laser–THz energy conversion efficiency was almost the same for different ${l}_{\mathrm{t}}$ , while the center frequency of the THz source decreased with the increase of ${l}_{\mathrm{t}}$ , as shown in Figure 5(b). For different ${l}_{\mathrm{t}}$ , ${v}_{\mathrm{t}}$ can be obtained from the corresponding simulations, and the simulation results also show ${l}_{\mathrm{s}}$ ranging from 0.5 to $5\;\mu \mathrm{m}$ . We assumed that $K$ decreased with the increase of ${l}_{\mathrm{s}}$ from 6.0 to 5.5 for the above cases, then we calculated the center frequencies of the THz pulse according to Equation (4) and plotted them in Figure 5(b). From the figure, one can see that the theoretical model developed here works very well in predicting the generation of center-frequency-tunable THz pulses with a laser-driven plasma wiggler.