Investigation of terahertz radiation generation from laser-wakefield acceleration Open Access

THz radiation, with a frequency range of 0.1–10 THz, is a non-ionizing type of electromagnetic (EM) radiation, which holds great potential for many scientific and technological applications.^1–3 Its ability to penetrate through certain materials makes it ideal for imaging, sensing, and security applications.¹ THz radiation also interacts strongly with molecules, which makes it useful for spectroscopy and the analysis of various substances, including biological molecules.^4,5 In addition, the potential for high-speed data transmission using THz radiation may revolutionize communication and information technology.⁶ 

The generation of THz radiation is also a topic of significant recent interest. Different techniques have been developed to generate THz radiation; each of them has its own features.^2,3,7–9 Diodes,¹⁰ quantum cascade lasers (QCLs),¹¹ and high-speed transistors¹² have been used to produce continuous wave (CW) THz radiation.¹³ Furthermore, various techniques, including photo-mixing, frequency multiplication, parametric conversion, and backward wave oscillators, have been developed to generate tunable CW THz waves.^1–3 Another technique to generate powerful THz radiation involves laser-plasma interaction.^14–25 This method utilizes solid,^19,26 liquid,^21,27 or gas^{8,15,16,28,29} targets.

Among many laser-plasma-based THz sources, the laser-wakefield acceleration (LWFA) method, using a gas target, can produce high-power THz radiation.^{8,15,16,28,29} For instance, a magnetized plasma can stimulate the resonant excitation of an electromagnetic plasma wakefield.^13,30–33 This wakefield partially transitions to an electromagnetic state with a nonzero group velocity, leading to the emission of electromagnetic pulses. The frequency of these pulses closely aligns with the plasma frequency. A key facilitator of this process is the presence of a plasma magnetized transversely.^30,31 Another example of an LWFA-based THz source is the coherent transition radiation (CTR) emitted by relativistic electron bunches, produced in LWFA, exiting a plasma-vacuum boundary.^16,34,35 Because the electron bunches are ultrashort (<50 fs), they can radiate coherently in a wide bandwidth (1–10 THz), yielding high-intensity THz pulses.³⁴ Pukhov et al. proposed that electrostatic plasma wakefield can radiate at harmonics of plasma frequency in a plasma with a positive density gradient based on nonlinear currents.³⁶ On the other hand, Sheng et al. explored the generation of powerful coherent THz radiation from LWFA through a linear mode conversion when the laser pulse is obliquely incident to a density gradient of an inhomogeneous plasma.³⁰ Qu et al. suggested that THz radiation can be produced by electrons refluxing along the sheath of the wakefield cavity (or the bubble) formed behind the laser pulse. Such refluxing electrons, moving along this “bubble shell,” have similar dynamic trajectories and thus can emit coherent backward THz radiation.³⁷ 

Recently, high-power THz generation in LWFA at 1–10 THz with an energy of up to 4 mJ per pulse was demonstrated by using a 150-TW laser.¹⁵ This process has yielded a laser-to-THz conversion efficiency of $≈0.15%$⁠, an indication that LWFA can be a strong THz source. In this study, we present a particle-in-cell (PIC) simulation model, which makes it possible to capture and interpolate THz radiation and exploit both fixed and moving windows in the simulation. We also report on evaluating the spectrum of THz radiation emitted during the interaction, the shape of the THz profile, and the energy and radiation angle of the THz pulse extracted from PIC simulations.

Our study explores the radiation produced by LWFA in the THz range. We consider a high-intensity femtosecond laser pulse (I > 10¹⁸ W/cm²) propagating through a helium (He) gas jet [see Fig. 1(a)]. The gas, ionized by the leading edge of the pulse, becomes a plasma. Then the laser’s ponderomotive force pushes the electrons, creating a space charge behind the pulse. This results in the formation of a wakefield, which can be understood as an electron density wave co-propagating behind the laser pulse, as seen in the electron density map in Fig. 1(b).

To simulate this process, we use a three-dimensional (3D) gas density distribution. The longitudinal density profile is designed to mimic experimental interferometric data,³⁸ and the transverse one is chosen to reduce the computational resources. We assume that the laser pulse initially propagates through vacuum over a distance of 200 μm along the x-axis. The neutral gas density along the beam propagation direction starts from zero and gradually increases to n₀ = 2.5 × 10¹⁸ cm⁻³ (equivalent to a plasma electron density of n_e = 5 × 10¹⁸ cm⁻³) over a distance of 600 μm. Subsequently, there is a density plateau region extending over 1800 μm, followed by a reduction in density to zero over 600 μm. Finally, the laser propagates an additional 600 μm in the vacuum. In the transverse direction, we model a plateau of 150 μm at the center, which gradually decreases to zero over 50 μm. The thickness of the gas is chosen to be ten times the laser beam’s waist (w₀). Furthermore, we extend the transverse size up to r = 207 μm [see Fig. 1(c)] to ensure that the electric field reflected from the boundary is negligible compared to the main signal. Both longitudinal and transverse density profiles have exponential falloffs to prevent the occurrence of sharp spatial charge gradients.

During the simulation, we capture electromagnetic fields to analyze the generated THz radiation. However, due to memory limitations, it is not feasible to store electromagnetic data at every position and timestep. To overcome this issue, we leverage the azimuthal symmetry of the wakefield. As depicted in Fig. 1(d), we define two types of probes to record the electric and magnetic components during the interaction. One-dimensional (1D) probes are aligned parallel to the x direction (i.e., the laser propagation direction) with its length equal to the longitudinal length of the simulation box (L_x). The two-dimensional (2D) probes are placed perpendicular to the laser propagation direction with a length of L_r (vertical length of the simulation box), as shown in Fig. 1(d).

We used the PIC code SMILEI³⁹ in cylindrical coordinates with the field decomposed into a basis of two azimuthal harmonic modes. A Gaussian laser pulse propagates through He gas with a₀ = 1.9, τ = 27 fs, w₀ = 15 μm, and λ_l = 800 nm, where a₀ = (eE₀)/(m_eω_lc) represents the normalized vector potential, where e is the electron charge, E₀ is the electric field amplitude, m_e is the electron mass, ω_l is the laser angular frequency, and c is the speed of light. τ denotes the pulse duration (FWHM), w₀ represents the beam waist, and λ_l is the laser wavelength. The laser pulse is focused onto the edge of the plateau region of the medium filled with He. We used the Ammosov–Delone–Krainov (ADK) ionization model to calculate the ionization rate. The simulation box has perfectly matched layer (PML) boundary conditions, and all parameters are normalized based on the laser frequency. The spatial and temporal resolution is set to Δz = λ_l/32, Δr = λ_p/120 (λ_p plasma period), and Δt = 0.86 times the Courant–Friedrichs–Lewy (CFL) limit (Δt = 0.86 × Δt_CFL). We use both fixed-window and moving-window simulation boxes, with box sizes of 152 576 × 1664 and 24 576 × 1568 cells, respectively. Each cell contains eight macroparticles (He atoms).

To detect the radiation in the desired THz range (corresponding to wavelengths between 2997 and 29.97 μm, equivalent to 0.1–10 THz), the simulation box should be properly sized. For example, a 3-mm-wide box can capture radiation down to 0.1 THz, which roughly matches our interaction length. A window length of 100 and 600 μm covers frequencies above 3 and 0.5 THz, respectively. In order to reduce simulation costs in PIC simulations, we implemented a moving window size, having a length of 600 μm. We also use a fixed window as a reference and compare it to a moving window with a length of 600 μm. This comparison aims to determine the minimum window size required to capture the electromagnetic field in the THz region while achieving satisfactory agreement in terms of spectrum and total THz energy. Such a comparison is essential as it allows us to reduce computational costs without significant information losses. About the simulation cost, the fixed window ran on 512 core CPUs for 511 h, and the moving window ran on 256 core CPUs for 180 h. This means the simulation cost of the fixed window is about 4.5 times higher than that of the moving window.

Figure 2(a) illustrates a distribution of electric fields pointed in the x direction (E_x) after the laser pulse travels over 2700 μm within a fixed window. The inset window reveals a wakefield structure composed of five oscillations of E_x behind the laser pulse. It is worth noting that the wakefield structure dissipates after a few hundred micrometers of the laser propagation. As a result, the acceleration, the derivative of the electron velocity dv/dt, is zero or minimal in the dissipated plasma region, indicating that no or little radiation is produced after the wakefield dissipates based on the Lienard–Wiechert potentials.⁴⁰ 

We present a snapshot of the electric field (E_x) distribution within a moving window at three specific times [before reaching the plateau, during the plateau, and within the down-ramp density region, as shown in Figs. 2(b)–2(d)]. It is found that a window length of 600 μm is enough to capture the main fluctuations observed in the laser wakefield process. To record the radiation generated from LWFA, we use 1D and 2D probes in our simulation. These probes, which are parts of the SMILEI code, allow us to interpolate the electric and magnetic fields and the current density at the probe locations. As the electromagnetic fields are semi-symmetric⁴¹ about the x-axis, we place 1D probes on a plane of y–z = 0. We employ an electromagnetic field sampling rate of 0.21 fs corresponding to 4.7 × 10¹⁵ samples per second. As a result, our analysis remains valid within the frequency range of up to 2350 THz. This satisfies the Nyquist criterion as the sampling rate exceeds twice the highest frequency component present within the signal, ensuring the reliability of our analysis.

For the spectral analysis of the EM signals, we focus on the data from the fourth 1D probe, which lies in the vacuum region. To compare the fixed-window and moving-window simulations, we plot spatially resolved spectra of the magnetic fields, B_x and B_y, along the laser beam propagation direction (x-axis) in Fig. 3. In addition, the lineouts of the spectra are shown at x = 800, 1800, and 2400 μm for B_x and B_y in the THz range. To obtain the spectra, we treat each point along the probe line as a detector that captures the magnetic field over time. The signal is then Fourier-transformed to obtain its power spectrum.

As shown in Figs. 3(a) and 3(b), the spectral map of B_x is similar for both fixed and moving windows, and so is that of B_y [see Figs. 3(e)–3(g)]. However, we observe some frequency leakage at the end of the interaction in the moving windows. This is not ideal because the simulation does not account for all data points, particularly toward the end. Figure 3 also shows that THz radiation is emitted continuously along the LWFA propagation direction. The intensity of the THz signal increases when the density increases and continues to increase in the plateau region where a₀ increases owing to relativistic self-focusing. As the number of relativistic electrons increases due to the increased laser intensity, the THz signal also increases. When both the laser intensity and the electron density decrease, the THz signal also decreases. As indicated in the lineout plot, the intensity of B_y is greater than that of B_x. This implies that the longitudinal acceleration is more dominant over the transverse one as the magnitude of the magnetic field is related to dβ/dt and β of electrons in LWFA by Lienard–Wiechert potentials. The THz spectrum calculated from the PIC simulation is consistent with the experimental results reported in Ref. 15.

Using the 2D probes, we can also calculate the Poynting vector in the x-direction, for example, when the laser passes by two positions, x = 970 and 3600 μm, as shown in Figs. 4(a)–4(d). At x = 970 μm, shown in Fig. 4(a), the energy flux is negative, indicating that the THz radiation propagates in the backward direction. At x = 3600 μm, the flux is positive, as shown in Fig. 4(d), implying forward THz radiation. When we examine a lineout at y = 0 in Figs. 4(b)–4(e), we find that the amplitude of the energy flux for the second 2D probe (forward radiation) is higher than that of the first one (backward radiation), indicating that forward THz radiation is dominant. The backward radiation is mostly emitted from the density up-ramp area. The forward radiation originates mostly from the plateau area, yielding a conical radiation pattern. This is consistent with the radiation pattern created by relativistic electrons accelerated in the longitudinal direction.

The spectra of the backward and forward radiation can be retrieved from the time-varying S_x field at the crossing points of the third (P₃) and fourth (P₄) 1D probes with the 2D probes [through the crossing points shown in Figs. 4(a)–4(d)]. The Fourier transforms of the fields are shown in Figs. 4(c)–4(f). The backward signal [see Fig. 4(c)] has a noticeable peak at around 18 THz, while the forward signal is peaked around the zero (DC) frequency, as shown in Fig. 4(f). Compared to Fig. 3, Fig. 4(f) shows rather broad spectra. This is because Fig. 4(f) presents the THz spectra obtained from near the end of the simulation box where the laser quickly expands and defocuses after exiting the plasma. Because of this, low-frequency (<3 THz) components produced earlier significantly diverge and are not well detected on the finite forward 2D probe. Interestingly, we observe several peaks in the backward direction. This can possibly be associated with the electron movement around the sheath of the bubble.³⁷ Qu et al. modeled these electron trajectories as synchrotron radiation, of which the fundamental frequency is given by ω₀ = v/R, where v and R are the velocity of the electron and the radius of the bubble, respectively.³⁷ 

Compared to the CTR mechanism, THz generation in our case has two main differences. In the case of CTR, the main source of THz is the high-energy electrons accelerated by the LWFA mechanism. However, in our case, it is the low-energy electrons that radiate THz waves. In addition, THz waves are produced continuously along the plasma channel in our case, whereas the CTR model predicts THz emission mainly from the end of the plasma (plasma-vacuum boundary). We also note that the generated THz waves can propagate through and escape from the overdense plasma because of the following reasons. First, the plasma frequency is reduced to $ωp=ωp0/γ$ in the relativistic regime, where γ is the relativistic Lorentz factor. In our case, a₀ changes from 1.9 to 4, and the corresponding $γ=1+a02$ value ranges from 2.15 to 4. At n_e = 5 × 10¹⁸ cm⁻³, the relativistic plasma frequency is f_p = 10–13.6 THz. Therefore, THz waves at >10 THz can propagate in the plasma. In addition, the radius of the plasma is small (∼25 μm), and THz waves whose wavelength is much longer than ∼25 μm can penetrate through the overdense plasma and reach the vacuum.

Figure 5 illustrates the spectra of the Poynting vector components (x, y, z) spatially resolved along the laser beam propagation (x-axis). Here, we calculate the Poynting vector in time at the fourth probe and then apply a Fourier transformation to obtain the spectrum. Figure 5 shows the resulting Fourier transforms for both fixed and moving window scenarios, with the frequency range selected below 30 THz. The dashed line represents the longitudinal gas density, and the white line indicates the summation of energy flux corresponding to electromagnetic frequencies below 10 THz.

The spectrum shapes in the fixed windows [Figs. 5(a), 5(c), and 5(e)] are similar to those in the moving windows [Figs. 5(b), 5(d), and 5(f)]. However, we see some difference at around x = 1800 μm. This point corresponds to the restarting point of the simulation in the moving window scheme, where a slight change in the electromagnetic field information can be amplified during PIC.

Both schemes show that the energy flux of THz radiation propagating in the x direction, S_x, increases in response to changes in gas density and laser intensity. The signal S_x starts to increase when the laser’s focus coincides with the edge of the density plateau. S_x peaks at around x = 2200 μm and then decreases as the laser intensity diminishes. As the laser enters the down-ramp, the THz energy gradually decreases. The transverse energy fluxes, S_y and S_z, also similarly vary along the laser propagation direction, but with about ten times weaker magnitudes. This is because the Poynting vector of THz radiation is mostly directed in the forward direction (x-axis). Based on the simulation results, the angle of THz radiation can be estimated as $θ=tan−1((Sy2+Sz2)/Sx)$⁠, which in our case yields θ = 8.5°–10°. Alternatively, the radiation angle can be obtained from Fig. 5(a), where the onset of the THz signal is shifted by $≈1100μ$m from the gas jet edge. From the radial position of the fourth probe (r = 207 μm), we obtain θ ≈ 10°.

In our model, THz radiation arises from a multitude of electrons suddenly accelerated by the ponderomotive force and plasma wakefield. These electrons, initially formed by the laser field, are abruptly repelled by the laser’s ponderomotive force. When they return due to electrostatic force, they experience significant scattering at the rear of the plasma bucket where electron charge density peaks. This radiation can be coherent and scale with the square of the charge as the source dimensions (8–20 μm) are much smaller than the THz wavelength (1000–50 μm). Moreover, the source spans the entire plasma length, allowing for constructive interference of radiation along the propagation direction, thereby generating a high-energy coherent THz pulse in the far field. Importantly, our PIC simulation confirms that the total charge of low-energy (∼MeV) electrons is of a few nC, sufficient to produce multi-mJ THz in the forward direction. These results are consistent with recent experimental data.¹⁵ 

This research offers a comprehensive understanding of the generation of THz radiation in LWFA through the use of a PIC simulation model. The PIC simulation model developed in this study serves as a powerful tool for capturing and interpolating THz radiation, utilizing both fixed and moving windows. Notably, a moving window of 600 μm allows us to interpret the continuous generation of electromagnetic waves at THz frequencies. Our simulation results show how the amplitude of THz waves evolves with longitudinally varying gas densities and laser intensities. Moreover, this study also helps to better understand the energy transfer mechanism from a laser pulse to plasma and THz radiation during LWFA.

This work was supported by the Institute for Basic Science, Korea, under Grant No. IBS-R012-D1 and the Ultrashort Quantum Beam Facility (UQBF) Operation Program (Grant No. 140011) through APRI, GIST. M.R.-P. expresses his heartfelt gratitude to Dr. Chul Min Kim for his invaluable assistance and support in this work. Computational works for this research were performed on the IBS Supercomputer, Aleph, at the IBS Research Solution Center. In addition, the authors thank F. Pérez and F. Massimo from the SMILEI team for their helpful insights regarding the PIC simulation.

The authors have no conflicts to disclose.

Mohammad Rezaei-Pandari: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Mohammad Mirzaie: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Calin Ioan Hojbota: Formal analysis (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Ali Reza Niknam: Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Reza Massudi: Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Ki-Yong Kim: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Chang Hee Nam: Funding acquisition (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.