Enhanced coherent transition radiation from midinfrared-laser-driven microplasmas

We present a particle-in-cell (PIC) analysis of terahertz (THz) radiation by ultrafast plasma currents driven by relativistic-intensity laser pulses. We show that, while the I0 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\lambda }_{0}^{2}$$\end{document}λ02 product of the laser intensity I0 and the laser wavelength λ0 plays the key role in the energy scaling of strong-field laser-plasma THz generation, the THz output energy, WTHz, does not follow the I0 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\lambda }_{0}^{2}$$\end{document}λ02 scaling. Its behavior as a function of I0 and λ0 is instead much more complex. Our two- and three-dimensional PIC analysis shows that, for moderate, subrelativistic and weakly relativistic fields, WTHz(I0 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\lambda }_{0}^{2}$$\end{document}λ02) can be approximated as (I0λ02)α, with a suitable exponent α, as a clear signature of vacuum electron acceleration as a predominant physical mechanism whereby the energy of the laser driver is transferred to THz radiation. For strongly relativistic laser fields, on the other hand, WTHz(I0 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\lambda }_{0}^{2}$$\end{document}λ02) closely follows the scaling dictated by the relativistic electron laser ponderomotive potential \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr{F}}_{{\text{e}}}$$\end{document}Fe, converging to WTHz ∝ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{0}^{1/2}{\lambda }_{0}$$\end{document}I01/2λ0 for very high I0, thus indicating the decisive role of relativistic ponderomotive charge acceleration as a mechanism behind laser-to-THz energy conversion. Analysis of the electron distribution function shows that the temperature Te of hot laser-driven electrons bouncing back and forth between the plasma boundaries displays the same behavior as a function of I0 and λ0, altering its scaling from (I0λ02)α to that of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr{F}}_{{\text{e}}}$$\end{document}Fe, converging to WTHz ∝ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{0}^{1/2}{\lambda }_{0}$$\end{document}I01/2λ0 for very high I0. These findings provide a clear physical picture of THz generation in relativistic and subrelativistic laser plasmas, suggesting the THz yield WTHz resolved as a function of I0 and λ0 as a meaningful measurable that can serve as a probe for the temperature Te of hot electrons in a vast class of laser–plasma interactions. Specifically, the α exponent of the best (I0λ02)α fit of the THz yield suggests a meaningful probe that can help identify the dominant physical mechanisms whereby the energy of the laser field is converted to the energy of plasma electrons.

2 0 scaling. Its behavior as a function of I 0 and λ 0 is instead much more complex. Our two-and three-dimensional PIC analysis shows that, for moderate, subrelativistic and weakly relativistic fields, W THz (I 0 2 0 ) can be approximated as (I 0 λ 0 2 ) α , with a suitable exponent α, as a clear signature of vacuum electron acceleration as a predominant physical mechanism whereby the energy of the laser driver is transferred to THz radiation. For strongly relativistic laser fields, on the other hand, W THz (I 0 2 0 ) closely follows the scaling dictated by the relativistic electron laser ponderomotive potential F e , converging to W THz ∝ I 1/2 0 0 for very high I 0 , thus indicating the decisive role of relativistic ponderomotive charge acceleration as a mechanism behind laser-to-THz energy conversion. Analysis of the electron distribution function shows that the temperature T e of hot laser-driven electrons bouncing back and forth between the plasma boundaries displays the same behavior as a function of I 0 and λ 0 , altering its scaling from (I 0 λ 0 2 ) α to that of F e , converging to W THz ∝ I 1/2 0 0 for very high I 0 . These findings provide a clear physical picture of THz generation in relativistic and subrelativistic laser plasmas, suggesting the THz yield W THz resolved as a function of I 0 and λ 0 as a meaningful measurable that can serve as a probe for the temperature T e of hot electrons in a vast class of laser-plasma interactions. Specifically, the α exponent of the best (I 0 λ 0 2 ) α fit of the THz yield suggests a meaningful probe that can help identify the dominant physical mechanisms whereby the energy of the laser field is converted to the energy of plasma electrons.
Identifying the physical limits of this method of THz generation would provide a deeper understanding of the fundamental aspects of relativistic laser-plasma electrodynamics and would help answer an urgent call of optical science for the development of high-power THz sources. Specifically, the generic I 0 2 0 scaling of the kinetic energy www.nature.com/scientificreports/ that electrons tend to pick up from the driver field with intensity I 0 and wavelength λ 0 gives grounds to expect that longer-λ 0 driver pulses, such as those available as an idler-wave output of high-peak power optical parametric chirped-pulse amplifiers (OPCPAs) [37][38][39] , could help enhance THz generation in relativistic laser-plasma settings. However, relativistic-intensity laser fields can give rise to a complex interplay of concurrent highly nonlinear motions of charged particles in a plasma target. This complex dynamics imprints itself onto a highly sophisticated picture of secondary radiation [27][28][29][30][31][40][41][42][43][44] , giving rise to long-wavelength radiation with widely different spatial, spectral, temporal and polarization properties. The intensity of THz emission in this complex laser-plasma interaction scenario depends on a variety of factors whose relevance to THz generation is sometimes difficult to assess and whose significance is often difficult to quantify. The actual behavior of the THz yield as a function of I 0 and λ 0 in such settings can significantly differ from the naïve I 0 2 0 expectation [45][46][47][48] . Many of the pressing questions related to the brightness scaling of relativistic laser-plasma THz sources, including the scaling of the THz yield as a function of I 0 and λ 0 , remain to be answered.
Here, we seek to address these questions by resorting to a particle-in-cell (PIC) analysis of laser-plasma interaction dynamics driven by relativistic-intensity ultrashort laser pulses in a metal-film target. We show that, while the I 0 2 0 product of the laser intensity I 0 and the laser wavelength λ 0 plays the key role in the energy scaling of strong-field laser-plasma THz generation, the THz output energy, W THz , does not follow the I 0 2 0 scaling. Its behavior as a function of I 0 and λ 0 is instead much more complex. Our two-and three-dimensional PIC analysis shows that, for moderate, subrelativistic and weakly relativistic fields, W THz (I 0 2 0 ) can be approximated as (I 0 λ 0 2 ) α , with a suitable exponent α, as a clear signature of vacuum electron acceleration as a predominant physical mechanism whereby the energy of the laser driver is transferred to THz radiation. For strongly relativistic laser fields, on the other hand, W THz (I 0 2 0 ) closely follows the scaling dictated by the relativistic electron laser ponderomotive potential F e , converging to W THz ∝ I 1/2 0 0 for very high I 0 , thus indicating the decisive role of relativistic ponderomotive charge acceleration as a mechanism behind laser-to-THz energy conversion. Analysis of the electron distribution function shows that the temperature T e of hot laser-driven electrons bouncing back and forth between the plasma boundaries displays the same behavior as a function of I 0 and λ 0 , altering its scaling from (I 0 λ 0 2 ) α to that of F e , converging to W THz ∝ I 1/2 0 0 for very high I 0 . These findings provide a clear physical picture of THz generation in relativistic and subrelativistic laser plasmas, suggesting the THz yield W THz resolved as a function of I 0 and λ 0 as a meaningful measurable that can serve as a probe for the temperature T e of hot electrons in a vast class of laser-plasma interactions. Specifically, the α exponent of the best (I 0 λ 0 2 ) α fit of the THz yield suggests a meaningful probe that can help identify the dominant physical mechanisms whereby the energy of the laser field is converted to the energy of plasma electrons.

The physical model and numerical setting
We consider a high-intensity laser pulse with a central wavelength λ 0 that drives a plasma slab (Fig. 1a, b) with a width l x = 16λ 0 , a variable depth l y , and initial electron density n 0 = 4n c , where n c = m e ω 0 2 /(4πe 2 ) is the critical plasma density, ω 0 = 2πc/λ 0 is the central frequency of the laser driver, c is the speed of light in vacuum, m e and e are the electron mass and charge. The initial electron velocity distribution inside the plasma slab is Maxwellian with an electron temperature of 230 eV. Such parameters are representative of plasmas generated on thin solid targets by high-intensity ultrashort laser pulses 40,[49][50][51] .
We simulate the solutions to the Maxwell-Vlasov equations for this laser-plasma interaction setting 40,52 using a 2D Smilei PIC code 53 . The solutions to the Maxwell equations are computed in this code with the finitedifference time-domain method in a two-dimensional configuration space and three-dimensional momentum space (2D3V) with the electromagnetic field discretized on a Yee grid with space steps Δx = Δy = λ 0 /128 along the x and y coordinates chosen, respectively, within the plane of the target and along the normal to this plane. Each grid cell contains 49 electrons and 49 protons, totaling to over 10 7 particles of each sort within the entire simulation domain. Time derivatives in the Maxwell equations were computed on a time grid with a step Δt = T 0 /256, T 0 being the driver field cycle.
The driver field in our simulations is taken in the form of an ultrashort Gaussian pulse with a pulse width τ 0 , an electric field amplitude E 0 , central frequency ω 0 , and variable polarization. The laser beam makes an angle φ 0 = 45 0 with the normal to the plasma surface (Fig. 1a, b) and is focused into a beam waist with a diameter d 0 = 4λ 0 on the plasma surface. Providing a meaningful quantifier of the significance of relativistic effects in this laser-plasma interaction setting is the normalized vector potential a 0 = eE 0 /(m e ω 0 c), which can be viewed as a ratio of the velocity v 0 = eE 0 /(m e ω 0 ) of nonrelativistic laser-driven electron quiver motion to the speed of light. The driver field is injected into the simulation domain via the Silver-Müller-type boundary conditions. PIC simulations were performed using graphic-processor-unit laboratory clusters 54 , as well as a shared research facility of high-performance supercomputing resources at M.V. Lomonosov Moscow State University 55 .

THz supercontinua and high-harmonic radiation
As it drives the plasma target, the laser field gives rise to intense radial waves of secondary emission (Fig. 1a-j). The low-frequency part of this radiation field (E x , E y , B z ) is isolated in our calculations by applying a super-Gaussian filter F e (ω) = exp[− (ω/ω c ) q ] to the total, driver + secondary radiation field ( E x , E y , B z ) as shown in Fig. 1a-e. We set q = 12 and ω c = ω 0 /4 in calculations presented here to provide a steep passband edge at ν c = ω c / (2π) ≈ 20 THz. Fully capturing the field structure of low-ω secondary radiation in 2D3V simulations is the B z field component, whose xy-plane maps are shown in Fig. 1f-j.
As the high-intensity laser pulse reaches the plasma surface at around t ≈ 15T 0 (Fig. 1b), the laser field starts to drive plasma electrons, modulating the electron density n e (x, y, t) inside the plasma and giving rise to signature ripples on the front surface of the plasma target (Fig. 2a). The ponderomotive force induced by the laser pulse drives electron oscillations, making plasma electrons bounce back and forth between the plasma boundaries. www.nature.com/scientificreports/ In Fig. 2a, b, we follow several oscillatory electron y(t) traces that start at different (y, t) space-time points at the center of the plasma target, x = 0. Shown against the map of the electron density, n e (x = 0, y, t), these traces illustrate how the amplitude of electron oscillations tends to grow as the electrons pick more energy toward the trailing edge of the pulse, following a similar behavior of the amplitude of the ripples in n e (x = 0, y, t).
Each time the laser-driven electrons leave or re-enter the plasma target as a part of their oscillatory motion (e.g., at around t/T 0 ≈ 21 and 30 in Fig. 2b), they emit ultrashort bursts of secondary radiation ( Fig. 1h-j), giving rise to bright THz supercontinua along with an intense short-wavelength output, dominated by high-order harmonics (HH) of the driver. To calculate the spectrum of this radiation, we choose to work in polar coordinates r and θ, such that x = rcosθ and y = rsinθ, fix r at r = 30λ 0 , and integrate over the r = 30λ 0 , y > 0 semi-circumference, is the Fourier transform and Q(t, r, θ) = B z (t, r, θ) and E z (t, r, θ) for p-and s-polarized field, respectively. The spectrum S(ω) spans the entire THz, infrared, and visible regions, Figure 1. (a-j) Spatial distributions of (a-e) the overall field B z (x, y) and (f-j) the THz-filtered field B z (x, y) and (k-o) electron y-p y phase-space distributions f(y, p y ) for a plasma slab with l y = 4λ 0 driven by a p-polarized field with a 0 = 1 and τ 0 = 80 fs at t = 0 (a, f, k), 15T 0 (b, g, l), 30T 0 (c, h, m), 40T 0 (d, i, n), and 50T 0 (e, j, o). www.nature.com/scientificreports/ extending deep into the ultraviolet range. In the spatial maps of the THz-filtered field B z (x, y) ( Fig. 1f-j), this process shows up as the first wave of THz radiation, whose emission starts as soon as the driver pulse strikes the plasma surface (cf. Fig. 1b, c, g, h). Presented in Fig. 3a, b are the temporal traces of specular-reflected THz radiation (blue line) and the overall, THz + HH field (green line) whose xy-plane snapshots are shown in Fig. 3e, g. The troughs of the first-wave bursts of THz radiation (blue line) are highlighted in these traces with pink arrows. As can be seen from these traces and the spatial field structure, emission of the first THz burst is accompanied by HH generation (HHG). The HH field is emitted in the form of a train of attosecond pulses, tightly confined to the central, most intense part of the laser pulse. Polarization properties of the HH field are in full obedience to the selection rules for relativistic HHG 40,56 -a p-polarized driver generates p-polarized even-and odd-order HHs; an s-polarized driver gives rise to s-polarized odd HHs and p-polarized even harmonics, and a circularly polarized driver generates even and odd harmonics in both p and s polarization states.

Subcycle CTR waves
The spectral properties and the spatiotemporal structure of broadband radiation in Figs. 1, 2 and 3 are understood in terms of a physical scenario that combines CTR by electrons traversing the plasma boundary [23][24][25][26][27][28]30,31 with HHG by laser-driven relativistic electrons near the front plasma surface 32,[37][38][39] . In this study, we focus on the THz part of broadband radiation emitted as a part of this complex laser-plasma interaction scenario. Central to THz generation is the ponderomotive force induced by the laser driver, which accelerates plasma electrons, making them move toward the rear surface of the plasma slab. While the laser-driven oscillatory motion of relativistic electrons near the plasma surface gives rise to attosecond pulses of HH radiation, the ponderomotive force induced by the laser driver accelerates plasma electrons, making them move toward the rear surface of the plasma slab. This process is readily discernible in the electron y − p y phase-space distribution f(y, p y ) as shown in Fig. 1k-o for a 0 = 1. Specifically, at t ≈ 15T 0 (Fig. 1l), some of the electrons near the front surface of the plasma slab at y = 0 are seen to start gaining momentum from the laser driver (marked with an arrow in Fig. 1l), which has just reached the y = 0 plasma surface. Because the laser ponderomotive force pushes electrons toward the inside of the plasma slab, i.e., toward negative y, the momentum gained by the electrons in Fig. 1l is negative. Accordingly, at later t, the electrons that have picked up their momentum from the laser driver are found deeper within the plasma target (Fig. 1m), with their energy further increased via the action of the laser ponderomotive force. Temporal traces of specular-reflected THz radiation (blue line) and the overall, THz + HH field (green line), (c, d) spectra of specular-reflected THz radiation, and (e-g) the snapshots of the spatial distribution of the THz-filtered field B z (x, y) taken at t = 50T 0 (e, f) and 65T 0 (g) for a plasma slab with l y = 4λ 0 (a, c, e, f) and 6λ 0 (b, d, g) driven by a p-polarized laser pulse with a 0 = 1 and τ 0 = 80 fs in (a-e, g) a full PIC simulation and (f) PIC simulation where electron recirculation is suppressed by halting all the electrons reaching the rear plasma boundary. Shown with arrows in the temporal traces are the instantaneous THz + HH bursts of secondary radiation (pink arrows) and recurring CTR waves (blue arrows). www.nature.com/scientificreports/ Near the rear surface of the plasma slab, however, the electrostatic sheath field gradually decelerates and eventually reflects some of these electrons, preventing them from leaving the plasma slab. This effect is readily visible in Fig. 1m, where electrons with positive p y appear near the rear surface of the plasma slab at y = −4λ 0 and move toward its front side, in accordance with the positive sign of p y , building up their momentum and gaining their kinetic energy in the process of this motion (the arrow in Fig. 1m). By the time these electrons reach the front surface of the plasma target at y = 0, the energy of some of these electrons is high enough to overcome the binding forces that tend to keep the electrons within the plasma. When such electrons cross the plasma boundary, they emit an outgoing radial CTR wave, clearly visible in Fig. 1i, j. The electrons whose energy is not sufficient to make it to the other side of the plasma boundary are reflected, once again, to continue their circulation within the plasma (the arrow in Fig. 1n). Such recirculating electrons give rise to recurring subcycle bursts of THz radiation in the y > 0 semi-space (highlighted with blue arrows in Fig. 3a, b) each time they approach the y = 0 plasma boundary with a sufficiently high kinetic energy.
The spatiotemporal structure of the THz field (Figs. 1f-j, 3a, b, e, g) is fully consistent with this picture of THz generation. While the first burst of THz radiation is generated as soon as the driver pulse strikes the plasma surface (cf. Figure 1b, g), the highest intensity of THz radiation in the direction of specular reflection is achieved not within, but in the wake of the laser pulse, at the trough of the first radial subcycle CTR wave (Fig. 1i, j), emitted at t 1 ≈ 0.70 ps (blue arrow in Fig. 3a), as a completion of the first cycle of electron recirculation within the plasma target driven by a laser pulse with λ 0 = 3.9 μm. This first burst of CTR is followed by a second CTR wave, which is observed as a well-resolved subcycle waveform in the temporal trace of the THz field in Fig. 3a with a trough at t 2 ≈ 0.88 ps. As can be seen in the spatial maps of B z (x, y) in Figs. 1j and 3e, this trough is separated by Δr ≈ 17.5λ 0 from the trough of the first CTR wave. When the thickness of the plasma slab is increased by a factor of 1.5, both the spatial separation Δr between the troughs of the first and second CTR waves (Figs. 3e, 3g) and the delay time δt = t 2 − t 1 between these troughs in the temporal trace of the THz field (Fig. 3a, b) increase by a factor of ≈1.5.
Because THz CTR waves are emitted within time intervals when laser-driven electrons traverse plasma boundaries, the visibility of the spatial crest-and-trough structure of the THz output is sensitive to the sharpness of the plasma boundaries. To quantify this effect, Fig. 4a-d present the spatial maps of the THz radiation field calculated for a plasma slab with the initial electron density profile defined as n e (x, y, t = 0) = n 0 for − 4λ 0 ≤ y ≤ 0 and n e (x, y, t = 0) = n 0 exp(− y/L) for y ≥ 0 (Fig. 4a). As the extension of the plasma transition layer L (plasma gradient length scale) increases from L = 0 to L = λ 0 /20, the essential features in the crest-and-trough structure of the THz output show little to no variation. This suggests a rather comfortable margin of L within which simulations performed within the approximation of infinitely sharp plasma boundaries continue to provide an adequate physical picture of THz radiation. While for gas targets, steep plasma profiles with L ≤ λ 0 /20 are hard to achieve, in laser-plasma experiments with solid targets, plasma density profiles with L ≤ λ 0 /100 have been repeatedly demonstrated [57][58][59] and even plasma gradients as steep as L ≈ λ 0 /200 have been achieved 60 .
In the frequency domain, the recurring bursts of THz radiation manifest themselves as clearly resolved fringes in the spectrum of THz radiation (Fig. 3c, d). Unlike the field intensity of the THz laser-plasma output, which is found, in close agreement with earlier studies 31 , to be a very weak function of the plasma depth l y (Fig. 4e), the spacing between the fringes in THz spectra is defined by the time interval δt between the CTR bursts, or the separation Δr between the CTR waves, Δν = c/Δr, and is therefore highly sensitive to l y . Specifically, for a plasma www.nature.com/scientificreports/ target with l y = 4λ 0 , the CTR waves are separated by Δr ≈ 17.5λ 0 , giving rise to fringes with Δν ≈ 4.4 THz in the spectrum of THz radiation (Fig. 3c). As l y is increased up to 6λ 0 , the spacing between the fringes in the THz spectrum decreases by a factor of ≈1.5, becoming Δν ≈ 3.0 THz (Fig. 3d). As a useful verification of the role of recirculating electron currents for THz generation, we performed a series of PIC simulations in which electron recirculation was suppressed by artificially halting all the electrons reaching the rear plasma boundary. Any artificial THz emission that such a halting of electrons may cause does not disturb the THz field outside the plasma, because the plasma, whose plasma frequency is in the 100-THz range in this regime, totally screens such artificial radiation. Specifically, with all the electrons reaching y = −4λ 0 brought to a halt, the most intense component of THz radiation, as is readily seen from a comparison of Fig. 3e, f, is completely suppressed.

The I 0 λ 0 2 scalability of the THz output and the electron energy distribution function
We now focus on the behavior of the THz yield in the considered laser-plasma setting as a function of the driver wavelength λ 0 and field intensity I 0 . To this end, we set a detection point at ξ = 30λ 0 along the ξ-axis in the direction of specular reflection (Fig. 1a, b, f) and define the field intensity as I THz = [c/(2μ 0 )][B 0 (ξ = 30λ 0 )] 2 , where B 0 is the peak B z field at the trough of the first CTR wave (blue arrows in Fig. 3a, b) and μ 0 is the magnetic permeability. In Fig. 5a, the THz intensity I THz is plotted as a function of the driver field intensity I 0 for two values of the driver wavelength, representing two types of relativistic-intensity short-pulse laser sources-Ti: sapphire lasers, λ 0 = 0.8 μm 61,62 , and high-peak-power mid-infrared OPCPAs, λ 0 = 3.9 μm 63-65 . Simulations here are performed for a laser driver with a pulse width τ 0 = 80 fs and a beam-waist diameter d 0 = 4λ 0 . The laser slab in these simulations has an initial electron density n 0 = 4n c , a width l x = 16λ 0 , and a depth l y = 4λ 0 . The laser and plasma parameters τ 0 , d 0 , n 0 , l x , and l y are thus defined in such a way that they change with the driver wavelength λ 0 . Specifically, the pulse width τ 0 = 80 fs corresponds to ≈ 30 field cycles for a laser driver with λ 0 = 0.8 μm and ≈ 6 field cycles for a λ 0 = 3.9 μm driver. For both driver wavelengths, the growth of the THz intensity as a function of I 0 is seen to significantly slow down as the driver field intensity approaches the a 0 = 1 relativistic borderline (the vertical dashed lines in Fig. 5a). Simulations performed within the range of laser pulse widths τ 0 from 30 to 250 fs, initial electron densities n 0 from 4 to 80n c , and plasma slab depths l y from 0.8λ 0 to 15λ 0 show that, when varied within these ranges, neither of these parameters has any significant effect on the behavior of I THz as a function of I 0 and λ 0 (Fig. 4e). In agreement with earlier studies 66,67 , an increase in n 0 was found to give rise to a growth of computational noise. www.nature.com/scientificreports/ To gain insights into this behavior of I THz (I 0 ), we resort to the analysis of the energy distribution function of the recirculating electrons at the moment of time when these electrons reach the y = 0 plasma boundary. In Fig. 5b, we present such energy distribution functions for I 0 λ 0 2 ranging from 1.4 × 10 18 W μm 2 /cm 2 (a 0 ≈ 1) to 1.4 × 10 20 W μm 2 /cm 2 (a 0 ≈ 10). We then define the temperature of CTR-emitting electrons, T e , by fitting the high-energy tails of these functions with Maxwellian distribution functions (shown by the pink lines in Fig. 5b). Remarkably, in the limit of a 0 > 1, the temperature T e obtained via such a procedure tends to follow the same scaling as a function of I 0 λ 0 2 as the THz intensity I THz (cf. blue filled circles and red open circles in Fig. 5c). Moreover, both I THz (I 0 λ 0 2 ) and T e (I 0 λ 0 2 ) are seen to converge in the high-a 0 limit to the scaling of Φ(I 0 λ 0 2 ) = F e / (m e c 2 ) (dotted line in Fig. 5a, c), as dictated by the laser pondermotive potential F e = m e c 2 (γ − 1) of a relativistic electron with rest energy m e c 2 and relativistic Lorentz factor γ.
To relate these tendencies in the behavior of the THz intensity to the physical picture of THz generation, we search for the best (I 0 λ 0 2 ) α fit for T e (I 0 λ 0 2 ), allowing α to take different values below and above the a 0 = 1 relativistic borderline. As shown in earlier studies [66][67][68][69][70] , for high, yet nonrelativistic driver intensities, with a 0 ≤ 1 or a 0 ~ 1, the temperature of hot electrons, T e , tends to scale as (I 0 λ 0 2 ) α , with α depending on the plasma gradient scale length and the incidence angle of the laser beam, varying from its lower bound at 1/3 to α > 1 in the case of laser vacuum heating 68,69 . However, for driver intensities well above the relativistic borderline of a 0 = 1, the energy of laser radiation is absorbed predominantly via a ponderomotive acceleration of electrons, leading to a T e ∝ (I 0 λ 0 2 ) 1/2 scaling 69,70 . The α exponent is thus a meaningful probe that can help identify the physical scenario whereby THz radiation is generated as a part of ultrafast laser-plasma electrodynamics.
In Fig. 5c, we plot the intensity of THz radiation as a function of I 0 λ 0 2 along with two best (I 0 λ 0 2 ) α fits for the a 0 ~ 1 and a 0 > > 1 regimes. That the best (I 0 λ 0 2 ) α fit for I THz (I 0 λ 0 2 ) near the a 0 = 1 borderline (shown by dashed vertical lines in Fig. 5a, c) regime is achieved with α = 1.2 (green dash-dotted line) is indicative of the predominance of laser vacuum heating as a mechanism of laser electron acceleration in our laser-plasma setting. For a 0 > > 1, on the other hand, a search for the best (I 0 λ 0 2 ) α fit for I THz (I 0 λ 0 2 ) converges to α = 1/2 (purple dash-dotted line), as expected for CTR-emitting electrons generated mainly by the laser ponderomotive force [66][67][68][69][70] . Both the electron temperature T e and the THz yield are seen to monotonically increase with λ 0 (Figs. 5a, 5c). At I 0 ≈ 8 × 10 17 W/cm 2 , the THz field in the trough of the first CTR wave induced by a λ 0 = 3.9 μm driver is seen to be an order of magnitude higher than the respective THz field from laser plasmas driven by λ 0 = 0.8 μm laser pulses (Fig. 5a). As the THz yield continues to grow with I 0 in the a 0 > 1 regime, THz fields with amplitudes above 0.2 MV/m can be generated from laser plasmas driven by mid-infrared pulses with a sufficiently high field intensity (I 0 ≥ 7 × 10 18 W/cm 2 for λ 0 = 3.9 μm in Fig. 5a). While the THz yield and the electron temperature T e in subrelativistic, a 0 < 1 plasmas with L ≤ λ 0 /20 follows the (I 0 λ 0 2 ) α scaling with α > 1 (Fig. 5a, c), as an indication of laser vacuum heating 68,69 , for plasmas with significantly longer L, with L being as long as a few λ 0 , plasma electrons gain much of their energy via resonant absorption 40 , leading to T e ∝ (I 0 λ 0 2 ) α with α ≈ 1/3 40,66 . For strongly relativistic plasmas, on the other hand, with a 0 > > 1, both the THz yield and the electron temperature T e are found to show little to no deviations from the (I 0 λ 0 2 ) 1/2 scaling regardless of L. While a p-polarized laser field can drive plasma electrons via both the ponderomotive force and vacuum heating, an s-polarized driver does not provide an effective field component that would induce vacuum-heating acceleration. The intensity of THz radiation from laser plasmas driven by a moderate-a 0 s-polarized laser pulse is therefore much lower (purple triangles in Fig. 5a) than the intensity of THz radiation generated by a p-polarized driver with the same a 0 (blue boxes in Fig. 5a). Yet, in the high-a 0 regime, where a ponderomotive acceleration dominates over vacuum heating, I THz (I 0 λ 0 2 ) is seen to converge to the scaling of Φ(I 0 λ 0 2 ) = F e /(m e c 2 ) (blue and purple dotted lines in Fig. 5a) regardless of the polarization of the driver.

Three-dimensional analysis
The purpose of the analysis presented in this section is to show that the key tendencies in the behavior of the THz laser-plasma output as a function of the parameters of the laser driver, such as I 0 and λ 0 , are in no way limited to the specific 2D oblique-incidence laser-plasma interaction geometry examined in the previous sections, but should be viewed instead as signatures of laser-driven plasma electrodynamics representative of a vast, rather general class of laser-plasma interaction settings. To reach this goal, we resort to a 3D PIC modeling of laser-plasma interactions driven by a relativistic-intensity laser driver propagating along the normal to the plasma surface. Due to its cylindrical symmetry, the field of such a driver can be decomposed into azimuthal modes 53 . With the y-axis chosen along the field propagation direction (Fig. 6a, b), as in the previous sections, a laser-driver field polarized along the z-axis is recognized as a pure azimuthal mode of order m = 1. Such a field is conveniently described in cylindrical coordinates r, θ, y (Fig. 6a) and is fully defined by m = 1 azimuthal-mode components E 1y (y, r), E 1r (y, r), E 1θ (y, r), B 1y (y, r), B 1r (y, r), and B 1θ (y, r). As such a laser field drives ultrafast electron oscillations within a plasma target, it gives rise to a pure m = 0 azimuthal-mode of THz radiation 53 with cylindrical field components E 0y (y, r), E 0r (y, r), E 0θ (y, r), B 0y (y, r), B 0r (y, r), and B 0θ (y, r). The azimuthal-mode decomposition of the fields and electron currents for laser-plasma interactions of such a symmetry thus helps to substantially reduce the computation complexity of the problem while preserving its 3D nature 53,71 . The plasma target in these simulations is taken in the form of a l y x l r = 2.5λ 0 × 25λ 0 cylinder with n e = 5n c (Fig. 6b), supporting the overall cylindrical symmetry of the problem.
In Fig. 6c-e, we present typical ry-maps of the E 0y (y, r), E 0r (y, r), and B 0θ (y, r) field components in the THz laser-plasma output. Clearly seen in these maps are two forward-propagating THz CTR waves, emitted by the rear plasma boundary, and one backward THz CTR wave, radiation by the front plasma boundary. Also readily discernible are spherical THz waves emitted by plasma edges. Radiation of THz waves of this class has been earlier identified as antenna-like emission by shielding electron currents flowing along the plasma surfaces as a part of laser-driven ultrafast charge-carrier dynamics inside the plasma 31  www.nature.com/scientificreports/ such dynamics is shown by red and blue coding in Fig. 6c-e. Similar to electric currents in an impulsively driven antenna 72 , such electron currents provide a source of outgoing spherical waves of electromagnetic radiation (clearly seen in Fig. 6c-e, as well as in Fig. 7a-j) as they are brought to a halt near the edges of the plasma target 31 . Shown in Fig. 7a-j are the ry-maps of the B 0θ (y, r) field component in THz radiation from the plasmas driven by a laser field with the laser vector potential ranging from a 0 = 1 in Fig. 7a to = 10 in Fig. 7j. As a 0 increases, antenna-like THz emission of laser-driven plasmas is seen to become more and more isolated in space from THz CTR waves, allowing an accurate separation of radiation energies, W 1 and W 2 , in the antenna-emission and CTR-wave components of the THz laser-plasma output. In Fig. 7k, we present these energies along with the total energy of the THz output, W THz = W 1 + W 2 , as functions of a 0 . The THz output energy is seen to grow with a 0 , reaching the level above 10 mJ for a 0 ~ 10.
The key tendencies in the behavior of W THz as a function of a 0 , are fully consistent, as a comparison of Figs. 5a, b, 7k, l shows, with the behavior of the THz radiation output in the 2D laser-plasma interaction geometry. Similar to the 2D laser-plasma setting, THz radiation in 3D simulations is emitted in both forward direction and in reflection. However, because the THz energy cannot be calculated in a 2D setting, where one dimension is missing, 2D scaling laws are formulated for the THz intensity rather than the THz energy. In the 3D setting, on the other hand, where the full THz radiation energy W THz can be calculated, analysis of the THz energy as a function of I 0 and λ 0 becomes possible.
Similar to the THz output in 2D simulations (Fig. 5a, b), the THz energy W THz in Fig. 7k is seen to rapidly grow with a 0 in the range of moderate a 0 , indicating vacuum heating as a predominant physical mechanism whereby the energy of the laser driver is transferred to THz radiation via electron acceleration. As a 0 grows above the a 0 > 1 level, W THz slows its growth, converging to a much more gently sloping Φ(I 0 λ 0 2 ) asymptotics (shown by the dashed line in Fig. 7k), thus showing once again a close similarity with 2D simulations (cf. Figures 5a, b  and 7k). As the laser-driver field approaches the a 0 ~ 10 level, W THz (a 0 ) follows the ∝ I 1/2 0 0 scaling with a very high accuracy, indicating, in full agreement with the 2D analysis, the decisive role of relativistic ponderomotive charge acceleration as a physical mechanism behind laser-to-THz energy conversion.

Conclusion
To summarize, we have shown that, while the I 0 2 0 product of the laser intensity I 0 and the laser wavelength λ 0 plays the key role in the energy scaling of strong-field laser-plasma THz generation, the THz output energy, W THz , does not follow the I 0 2 0 scaling. Its behavior as a function of I 0 and λ 0 is instead much more complex. Our two-and three-dimensional PIC analysis shows that, for moderate, subrelativistic and weakly relativistic fields, W THz (I 0 2 0 ) can be approximated as (I 0 λ 0 2 ) α , with a suitable exponent α, as a clear signature of vacuum electron acceleration as a predominant physical mechanism whereby the energy of the laser driver is transferred to THz radiation. For strongly relativistic laser fields, on the other hand, W THz (I 0 2 0 ) closely follows the scaling dictated by the relativistic electron laser ponderomotive potential F e , converging to W THz ∝ I 1/2 0 0 for very high I 0 , thus indicating the decisive role of relativistic ponderomotive charge acceleration as a mechanism behind laser-to-THz energy conversion. Analysis of the electron distribution function shows that the temperature of hot laser-driven electrons bouncing back and forth between the plasma boundaries displays the same behavior as a function of I 0 and λ 0 , altering its scaling from (I 0 λ 0 2 ) α to that of F e , converging to W THz ∝ I 1/2 0 0 for very high I 0 . Specifically, as can be seen from the comparison of THz yields attainable from laser-plasma sources driven with 0.8-and (c-e) Spatial distributions of the E 0y (y, r) (c), E 0r (y, r) (d), and B 0θ (y, r) (e) components of the THz field for t = 90T 0 and a 0 = 1. The dashed line is the FWHM radius of the driver laser beam propagating in space with no plasmas. Also shown is laser-driven electric charge separation in shielding plasma currents (red and blue coding is for the positive and negative charges, respectively), giving rise to antenna-like emission of THz radiation from plasma edges. www.nature.com/scientificreports/ 3.9-μm laser pulses (green and blue boxes in Fig. 5a along with the respective Φ(I 0 λ 0 2 ) asymptotes), the use of a 3.9-μm, sub-100-fs output of the latest-generation high-power OPCPAs 37-39,63-65 can significantly enhance THz generation from relativistic laser-plasma settings relative to plasmas driven by standard, 0.8-μm Ti: sapphire laser pulses. Relativistic field intensities have been already achieved for such sources 63,64 . The work to even higher I 0 is in progress. These findings provide a clear physical picture of THz generation in relativistic laser plasmas, suggesting the THz yield W THz resolved as a function of I 0 and λ 0 as a meaningful measurable that can serve as a probe for the temperature T e of hot electrons in a vast class of laser-plasma interactions.

Data availability
All data generated or analyzed during this study are included in this published article.