HTML

Highintensity ultrashort laser pulses become inseparable from the future development of industrial applications in 3D nano/micromachining/printing due to Moore’s lawlike increase of their average power over the last years from 2000^{1}. Deterministic energy deposition via highly nonlinear multiphoton and avalanche ionisation is key to 3D CNC machining metals and dielectrics from TW/cm^{2} to PW/cm^{2} intensity window^{2}. The laser ablation thresholds in the case of ultrashort 10–50 fs pulses attracted interest due to better control of laser machining depth and removed volume as observed in dielectrics^{3}, however absent in metals^{4}. For example, in fused silica SiO_{2}, the ablation threshold was 1.3 J/cm^{2}, i.e., 1.5 times smaller than for longer fspulses, and the most efficient rate of ablated volume per pulse was at 4 J/cm^{2} fluence or 0.57 PW/cm^{2} intensity for 7 fs pulses^{3}. Formation of an overcritical (reflective) plasma mirror at the central part of the Gaussian intensity profile defined energy deposition efficiency and high precision ablation depth control with tensofnm per pulse. The ablated volume per energy saturates at the 30 fs pulse duration. Plasma mirror is proposed for compression of ultrashort laser pulses to reach the exawatt intensity scale^{5}; 1 EW = 10^{18} W.
A 3D printing/polymerisation at ~ TW/cm^{2} becomes ionisation controlled, and absorption defines polymerisation rather than the chemical composition of photoinitiators and twophoton absorbers. Deterministic energy deposition occurs at higher intensities corresponding to the tunnelling ionisation of molecules and atoms at ~ 10 PW/cm^{2}^{6}. In ambient air, focusing ultrashort pulses of a few tens of fs into spots of a few micrometers in diameter is possible by reflective optical elements, which deliver intensities sufficient for the tunnelling ionisation^{7}. Such conditions were explored in this study for THz radiation using instantaneous air breakdown along pulse propagation throughout the focus, which was subwavelength for emitted THz radiation.
THz spectral range 0.1–30 THz (3 mm – 10 μm in wavelength) has fundamental importance for material science due to characteristic lattice vibrations of phonon modes. Among established semiconductorbased THz emitters/receivers working on photoactivated optically gated currents (Auston switch), a new emerging direction is using nanothin topological insulators^{8}. Another emerging direction is the use of dielectric breakdown as a tool to create shortlived (~ 1 ps) and microlocalised (subwavelength at THz frequencies) currents for THz emission^{9}. Coherent directional emitters at farIR wavelengths can be realised using Wolf's effect, where surface subwavelength structures of period $ \Lambda\approx 10\; {\text μm} $ for farIR black body radiation $ \lambda_{\rm{black}} $ are defined to extract emission of the surface phonon polariton (SPP) mode $ \lambda_{\rm{SPP}} $ at an angle $ \theta $ from the normal to the crystal surface according to the phase matching condition for wavevectors $ k=2\pi/\lambda $ along the surface: $ k_{\rm{black}}\sin\theta = k_{\rm{SPP}} + 2\pi/\Lambda $^{10}. In this way, acoustic/optical photon THz frequencies are extracted into freespace radiation, which is, additionally, coherent due to the subwavelength character of the grating.
The singlecycle highintensity electrical THz field of 1 MV/cm was demonstrated with fs tiltedfront pulses using optical rectification in LiNbO_{3} crystal^{11}. Periodically poled lithium niobate was used to produce THz emission via walkoff of the optical pump and THz wave^{12}. High THz pulse energy of 10 μJ at a peak field strength of 0.25 MV/cm, was produced by optical rectification when phononpolariton phase velocity was matched with the group velocity of driving fslaser pulse and reached high visibletoTHz photon conversion of 45%^{13}. Tunable 10–70 THz emitters were realised via differencefrequency mixing of two parametrically amplified pulse trains from a single whitelight seed in GaSe or AgGaS_{2} crystals and reached 100 MV/cm peak field at 10^{8} W peak power (19 μJ energy) using tabletop laser system^{14}. Such high powers are usually produced by free electron lasers (FELs). The THz radiation is related to fast electrical current as per the Fourier transform, i.e., 1 ps current generates 1 THz radiation. SpinHall currents in trilayers of 5.8nmthick W/Co_{40}Fe_{40}B_{20}/Pt were used as THz broadband 1–30 THz emitters generating field strength of several 0.1 MV/cm and outperformed ZnTe $ \left<110\right> $ emitters^{15}. An electrical conductivity induced via the Zenner tunnelling can be harnessed to induce electrical transients in semiconductors and dielectrics by strong optical fields. When an optical field of 170 MV/cm was applied to the Auston switch with a 50 nm gap on SiO_{2} by sub4 fs pulse of 1.7 eV light, the currents were detected due to induced electron conductivity of 5 Ω^{–1}·m^{–1}, which is 18 orders of magnitude higher than the static d.c. conductivity of amorphous silica^{16}. Optical control of phonon modes and Zenner tunnelling can be operated without the dielectric breakdown maintaining solidstate material. Highintensity ultrashort irradiation of gasses, liquids, and solids at higher intensities, when the dielectric breakdown occurs, can also be used for producing highintensity THz fields^{17}. By copropagation of the first and second harmonics beams focused with f = 100 mm lens into long filaments in air plasma, the ACbias conditions, THz emitters were realised delivering peak intensities up to 10 kV/cm^{18−20}. Polarisation of the emitted THz radiation can be set circular by introducing noncollinearity between the filaments^{21}. Irradiation of Hejet in vacuum was shown to produce coherent THz emission via the transition radiation of ultrashort electron bunches produced by the laser wakefield, which, in turn, was produced by 8 TW peak power pulses of 800 nm/50 fs^{22}. The short length of the electron bunch was responsible for the coherent emission. It was shown recently that the use of prepulses can provide control and stability of the wakefield formation^{23}. It can also be useful for THz generation using transition radiation. The field of THz sources has become an active research frontier^{24−26}.
Here, we used a focused fslaser pulse for air breakdown from subwavelength volumes, which can potentially lead to the engineering of highintensity THz emitters due to the coherent addition of THz electrical fields and ultrashort duration of a few optical periods. Two prepulses used in this study revealed a position dependence on the THz emission and its polarisation, extending from the previous study where one prepulse was used^{27}. Experimental conditions of focusing and THz detection were kept the same^{27}, and only two prepulses were introduced using a holographic beam shaping technique based on a spatial light modulator (LCSLM)^{28}. For consistency, the same pulse energies were used for pre and mainpulses: 0.2 and 0.4 mJ (13 GW/pulse), respectively^{27}. The range of THz emission was 0.1–2.5 THz^{27}. The experimental results were interpreted with a theoretical model accounting for the air ionisation and free electron heating, generation and propagation of a shockwave, and spontaneous generation of the magnetic field due to coupled thermal and density gradients at the shock front. Interestingly, two or threepulse irradiation with air breakdown and shock formation is approaching the frequencies <0.1 THz of 5G network 24.25–52.6 GHz (Technical Specification Group Radio Access Network TS 38.104).

Twopulse experiment with pre and main pulses^{27} revealed the required timing and spatial separation for the shockedcompressed air by the prepulse to be pushed into the position of the main pulse at $ (x,y) = (0,0) $ for the maximum emission of THz radiation. Prepulse had to arrive 9.7 ns earlier to the focus of the main pulse from the position $ \Delta y\sim 33\; {\text{μm}} $ away for the maximum THz emission. Polarisation of THz radiation was found aligned along the direction of focal positions of prepulse and main pulse^{27}. Here, we add a second prepulse of the same 0.2 mJ energy, which arrives at the same time as the first prepulse, with the possibility of placing it freely in 3D space $ (x,y,{\textit z}) $. Fig. 1 shows the scheme of the experiment and the diagnostics. Details of the experimental setup and the measurement methodology are described in Sec. “Materials and methods”. The effect of placement of the second prepulse on polarisation and intensity of THz radiation was systematically studied. Due to short 35 fs pulse (~120 μm in length) and zposition (along the beam propagation) change only tensofmicrometers. Hence, the timing of both prepulses can be considered simultaneous, regardless of the actual zcoordinate of the prepulse 2 within an error of ~0.1 nm, which is by an order of magnitude smaller than the shockwave thickness.
Fig. 1 Schematic diagram of the experimental setup with temporal ($ \Delta t $ = 9.7 ns) and spatial offsets of the the prepulse 1 ($ \Delta x_1 $, $ \Delta y_1 $, $ \Delta {\textit z}_1 $), and the prepulse 2 ($ \Delta x_2 $, $ \Delta y_2 $, $ \Delta {\textit z}_2 $) irradiation positions. THz electrical fields $ E_x, E_y $ are calculated from $ E_{+45} $ and $ E_{45} $ measurements, where $ \pm 45^\circ $ is the angle of the polariser wiregrids set out for the polarisation measurement: $ E_x=E_{+45}+E_{45} $ and $ E_y=E_{+45}E_{45} $. The prepulse position was also tuned to place focus at the upstream (positive $ \Delta {\textit z} > 0 $) and downstream $ \Delta {\textit z} <0 $ locations; see bottomright inset. The standard polarised THzTDS detection system was based on a pair of wiregrid polarisers in front of ZnTe crystal, a quarterwave plate (QWP), a Wollaston prism (WP) and a balanced photodetector (BPD). Lowerinset shows 3D positions of the focal region of both prepulses with shockwave fronts (geometry is based on optical shadowgraphy imaging).
$ Y $position of prepulse 2 Fig. 2 shows the emission of THz radiation ($ E_y $ component) for the situation when the prepulse 1 is focused at $ \Delta y = +33\; {\text {μm}} $. At the same time, the second prepulse is focused at other $ (0,\pm y) $ position with respect to the mainpulse focus at $ (x,y)=(0,0) $; see Fig. 1. Change of $ E_y $ field orientation was observed over the measured range of THz frequencies as shown by the time domain transients in Fig. 2a. When two prepulses were simultaneously irradiated at opposite (with respect to $ y=0 $) positions, the cancellation of THz $ E_y $ components occurred. This confirms that the air density gradient induced by the shock front is more important than the air mass density, which is increased by two colliding shocks.
Fig. 2 The THzTDS waveforms of a $ y $component $ E_y $ of the single prepulse when the delay time between the prepulse and main pulse irradiation is at 9.7 ns. The spatial offset for the prepulse irradiation is ($ \Delta x $, $ \Delta y $) = ($ 0\; {\text μ} $m, $ +33\; {\text μ}$m) for the blue solid line or ($ 0\; {\text μ} $m, $ 33\; {\text μ}$m) for the red dotted line. b The $ x $component $ E_x $ of the dual prepulses when ($ \Delta x_1 $, $ \Delta y_1 $, $ \Delta {\textit z}_1 $) = ($ 0\; {\text μ} $m, $ +33\; {\text μ}$m, $ 0\;{\text μ} $m) and ($ \Delta x_2 $, $ \Delta y_2 $, $ \Delta {\textit z}_2 $) = ($ 0\; {\text μ} $m, $ 33\; {\text μ} $m, $ 0\; {\text μ} $m).
For two prepulse experiments, a dependence of THz emission on the offset $ \Delta y_2 $ of the second prepulse was determined. The delay of the first prepulse of 9.7 ns was optimised, so the THz emission was the most intense. The THz intensity is calculated from the $ (E_x,E_y) $ components as peaktopeak value (see inset in Fig. 3). A pictorial diagram for each irradiation point is presented by the projection of THz Efield on $ (E_x,E_y) $plane. The projection shows the overall temporal evolution of THz electric field orientation during emission time.
Fig. 3 a THz polarisation evolution in time shown as a $ (E_x,E_y) $ map. b Geometry of main pulse (m) and prepulses 1,2 (p1,p2). Prepulse 1 and the main pulse were fixed, and prepulse 2 was scanned along the $ y $axis. c The THzpolarisation $ E_{{\rm{peak}}{\rm{to}}{\rm{peak}}}^2 $ intensity as a function of the vertical offset for the second prepulse, $ \Delta y_2 $. $ E_{{\rm{peak}}{\rm{to}}{\rm{peak}}}^2 $ is as shown in the inset and is equal to $ \Delta E_x^2 $ + $ \Delta E_y^2 $. The delay time between the prepulse and main pulse irradiation is 9.7 ns. The spatial offset for the first prepulse irradiation is fixed at ($ \Delta x_1 $, $ \Delta y_1 $) = ($ 0\; {\text μ} $m, $ +33\; {\text μ} $m) while the position of the second prepulse irradiation, $ \Delta y_2 $ and $ \Delta x_2 $ was varied from 0 to $ 30\; {\text μ} $m. The blue THzpol. plot indicates the results from two prepulse irradiation with both prepulses at 0.2 mJ. The black THzpol. plot indicates the results from two prepulse irradiation with the prepulses at 0.2 mJ and the second prepulse at 0.05 mJ. The red THzpol. plot indicates single prepulse irradiation condition with the prepulse at 0.2 mJ and ($ \Delta x $, $ \Delta y $) = ($ 0\; {\text μ} $m, $ +33\; {\text μ} $m). The main pulse energy is 0.4 mJ. The red dotted line indicates the THz emission from a single main pulse irradiation in air.
One prepulse produced enhancement of THz emission by 13 times, i.e., the value $ E_{{\rm{peak}}{\rm{to}}{\rm{peak}}}^2 = 53\; {\text μ} $V^{2} compared to the value $ E_{{\rm{peak}}{\rm{to}}{\rm{peak}}}^2 \approx 4\; {\text μ} $V^{2} without a prepulse. The $ E_y $ component of the THz field becomes larger for offsets $ \Delta y_2>33\; {\text μ} $m compared to the enhancement produced by the first. Conversely, when prepulse 2 was closer to the position $ (0,0) $ where the main pulse is focused, the amplitude of THz field was reduced.
The second prepulse of four times smaller energy of 0.05 mJ, had no significant effect on the THz emission, regardless of the position at $ y<0 $, opposite side from prepulse 1 with respect to the main pulse focus. A different evolution in time for such 0.05 mJ second prepulse was observed. It showed a $ \sim 6 $ ps duration of larger negative values of $ E_x $ component (see inset in Fig. 3c and colour asymmetric $ (E_x,E_y) $ markers). This is probably due to a later arrival and shorter contribution of a shock wave triggered by a four times lower prepulse energy. Indeed, the radius $ R $ of a cylindrical shock wave created by the laser energy deposition on the axis $ E_p $ increases with time $ t $ as $ R(t)\propto E_p^{1/4}\sqrt{t} $ according to the Sedov–von Neumann–Taylor model^{29}. The major contribution to the THz enhancement was governed by the first prepulse with an energy of 0.2 mJ.
Interestingly, when prepulse 2 was on the same side as prepulse 1 at a distance from the main pulse position larger than 33 μm (prepulse1), the $ E^2 $ enhancement was $ \sim 15 $ times (see Sec. “Theoretical description of lightair interaction and THz emission” for discussion). In this case, the shockaffected volume from the prepulse 1 was crossed by the prepulse 2. When a shock wave propagates through a region of rarefied density $ \rho_1 $ formed by the other prepulse, the same cylindrical explosion model predicts scaling $ R\propto\rho_1^{1/4} $. When prepulse 2 arrived at $ y = 0 $ to $ 33\; {\text μ} $m (the opposite side of the $ y $axis as compared to prepulse 1), it quenched the effect of THz enhancement. However, at even more negative positions of prepulse 2, the same enhancement was observed for arriving shocks to the position of the main pulse; see red dotdashed line drawn as eye guides in Fig. 3c.
$ X $position of prepulse 2 Next, the dependence of THz emission of the $ x $position of the second prepulse was determined, see Fig. 4. The $ (E_x,E_y) $ trajectories indicate the radial influence of the prepulse irradiation position with respect to the main pulse (the prepulse 1 was always at the $ (0,33)\; {\text μ} $m). The ellipticity of THz polarisation was larger for the prepulse 2 locations (larger $ \Delta x_2 $). This is explained by the optimal timing of 9.7 ns for the prepulse position at $ \sim 33\; {\text μ} $m from the main pulse position. This is valid along any direction perpendicular to beam propagation in the $ (x,y) $ plane. When $ x_2<33\; {\text μ} $m, the shock arriving along the $ x $axis does not fit the perfect timing for THz enhancement and alters the effect of the prepulse 1 arriving along a perpendicular direction. When $ x_2>33\; {\text μ} $m, prepulse 2 makes a small contribution to the enhancement determined by prepulse 1. The prepulses 1 and 2 arriving at perpendicular directions define the ellipticity of the THz radiation and its handedness (see Fig. 4). The duration of the most intense THz radiation lasts approximately one cycle, corresponding to $ 2{\text π} $ rotation in the $ (x,y) $ plane (see a sideview crosssection along the time axis in the inset of Fig. 4).
Fig. 4 The THzpolarisation $ E_{{\rm{peak}}{\rm{to}}{\rm{peak}}}^2 $ [μV^{2}] intensity as a function of the horizontal offset for the second prepulse, $ \Delta x_2 $. The delay time between the prepulse and the main pulse is 9.7 ns. The spatial offset for the first prepulse irradiation is ($ \Delta x_1 $, $ \Delta y_1 $) = ($ 0\; {\text μ} $m, $ +33\; {\text μ} $m), while the second prepulse is focused at different $ \Delta x_2 $ position with $ \Delta y_2 $ = $ 0\; {\text μ} $m (see topright inset for geometry schematics). The RGB projections ($ E_x,E_y $) of THz plots are from two prepulse irradiation with both prepulses at 0.2 mJ. The THzpol. plot, when the prepulse is at 0.2 mJ and ($ \Delta x $, $ \Delta y $) = ($ 0\; {\text μ} $m, $ +33\; {\text μ} $m), indicates single prepulse irradiation condition with the main pulse energy is at 0.4 mJ. The red dotted line indicates the THz emission from single main pulse irradiation in air at $ E^2\approx 4\; {\text μ} $V^{2}. The boxedarea (yellow) is $ 6\times 6\; {\text μ} $V for the $ (E_x,E_y) $projection maps located at the crossmarkers, at which the total intensity $ E_{{\rm{peak}}{\rm{to}}{\rm{peak}}}^2 = E_x^2 +E_y^2 $ is plotted after power calibration. The 3D projection insets illustrate RGB time evolution over a span −5 ps to 15 ps and reveal the sense of polarisation rotation (looking into the beam or time axis); it is plotted for the closest projection RGB markers, respectively.
The dependence of THz emission on the position of the prepulse along the z axis is presented in Appendix “Zposition of prepulse 2”.

Laser air ionisation and electron heating Here, a physical scenario of air breakdown by the prepulse and main pulse and the dynamics of the plasma formation are analysed. The laser pulse parameters, the air ionisation potentials and the main experimental results are summarised in Appendix “Input parameters and experimental results for theoretical estimations”.
Laser focusing in our experiment was carried out with an offaxis parabolic mirror with an effective numerical aperture $ NA = 0.125 $ as determined from direct imaging of waterjet breakdown^{30}. The focal spot area is estimated from the Airy disk radius $ w_0=0.61\lambda/NA\approx 3.9\; {\text μ} $m, which is close to optical observation by shadowgraphy at low pulse energies. The axial extension of the air breakdown zone in the experiment is $ \sim 120150\; {\text μ} $m ^{27}. It is close to the depthoffocus $ 2{\textit z}_R $, where the Rayleigh length $ {\textit z}_R = \pi w_0^2/\lambda \approx 60\; {\text μ} $m. The focusing volume is considered as a cylinder of length $ 2{\textit z}_R $, radius $ w_0 $ and volume $ V_{\rm{las}} =5.7\times 10^{9} $ cm^{3}.
The energy of the prepulse $ E_1=0.2 $ mJ and the mainpulse $ E_m=0.4 $ mJ correspond to the intensity $ 1.2\times 10^{16} $ W/cm^{2} and $ 2.4\times 10^{16} $ W/cm^{2}, respectively. Under these conditions, the laser ponderomotive potential is much larger than the ionisation potential of atoms, the air is ionised in the tunnelling regime, and a plasma is created in the focusing volume.
The plasma creation is estimated for the air consisting of 20% of oxygen and 80% of nitrogen molecules. The ionisation rate $ w_i(E) $ as a function of laser electric field $ E(t) $ and atomic charge $ i $ is taken according to the AmmosovDeloneKrainov (ADK) theory^{31}. For the sake of simplicity, we assume a consecutive ionisation from unperturbed atomic levels.
Electron ionisation and heating are described by the kinetic equation for the distribution function of free electrons, $ f_e=\sum\limits_j f_{e,j} $, where $ f_{e,j} $ is the electron partial distribution function corresponding to the atom ionisation from the $ j $level:
$$ \partial_t f_{e,j} \frac{e}{m_e}{\bf{E}}\cdot \partial_{\bf{v}}f_{e,j}= (1n_{e,j}) w_j(E) \,\delta({\bf{v}}) $$ (1) where $ f_{e,j} $ is normalised by the density of neutral atoms, $ n_a $, $ E(t) $ is the laser electric field, $ n_{e,j}= \int dv\,f_{e,j} $ is the probability of ionisation from the level $ j $, $ \delta({\bf{v}}) $ is the delta function and $ {\bf{v}} $ is velocity. In the case of linear polarisation, nonrelativistic electric field, and neglecting electron collisions, the electrons are moving in the direction of the laser field, and the distribution function is onedimensional.
This equation describes the production of free electrons and their energy increase. The energy is taken from the laser field; therefore, the laser absorption and atom ionisation proceed simultaneously. Integrating kinetic equation, one obtains equations for the three moments of electron distribution function, $ n_{e,j} $, $ n_{e,j} v_{e,j}=\int dv\,v\,f_{e,j} $ and $ n_{e,j} v_{e,j}^2 =\int dv\,v^2 f_{e,j} $:
$$ \begin{aligned} {\partial_t n_{e,j}}&{=(1n_{e,j} ) w_j (E)} \\{ \partial_t (n_{e,j} v_{e,j} ) }&{=\frac e{m_e} n_{e,j} E} \\ {\frac12 \partial_t (n_{e,j} m_e v_{e,j}^2 ) }&{=en_{e,j} v_{e,j} E } \end{aligned} $$ (2) These equations are solved numerically using Wolfram Mathematica package for each ionisation level and each atom separately, and the average ionisation and the average residual electron energy are evaluated as
$$ Z_e= \sum\limits_j n_{e,j}, \qquad \varepsilon_e= Z_e^{1}\sum\limits_j \frac12 n_{e,j} m_e v_{e,j}^2 $$ (3) Two examples of numerical solutions are shown in Fig. 5 for $ j=2 $ (left) and $ j=5 $ for the main pulse and the nitrogen atom. The electron residual energies are on the order of ionisation potential and significantly smaller than the ponderomotive energy. This is in agreement with the theoretical analysis^{32}.
Fig. 5 Ionisation probability (red) and the electron kinetic energy for the parameters of the main pulse interacting with the nitrogen ion and liberation an electron from the level $ j=2 $ (left) and 5 (right). Pulse duration is 35 fs at FWHM; it is shown with the dashed lines.
The total ionisation energy per atom (assuming 80% nitrogen + 20% oxygen), the number of free electrons per atom, and the residual average electron energy are: for the prepulse $ \varepsilon_{\rm{ion}}=127.2 $ eV, $ Z_e=4.06 $, and $ \varepsilon_e=35.7 $ eV; for main pulse $ \varepsilon_{\rm{ion}}=244.8 $ eV, $ Z_e =4.85 $, and $ \varepsilon_e= 74.0 $ eV. Multiplying the total energy loss per ion, $ \varepsilon_{\rm{tot}}= \varepsilon_{\rm{ion}}+ Z_e\varepsilon_e $, by the ion density and the focal volume, we obtain the free electron density and the laser energy loss for ionisation and electron heating. For the prepulse propagating in the air at normal conditions, the electron density is $ n_e=Z_e n_a= 2.2\times 10^{20} $ cm^{−3} and the energy loss is $ \varepsilon_{\rm{tot}} n_a V_{\rm{las}}= 13.4\,{\text μ} $J. The main pulse is propagating along the shock front in compressed air. The compression factor $ C_r = 5.4 $ was measured in Ref. 27. Therefore, the electron density in the focal zone is $ n_e = 2Z_e C_r n_{\rm{air}} = 1.4 \times 10^{21} $ cm^{−3}, and the laser energy loss is $ \varepsilon_{\rm{tot}} C_r n_a V_{\rm{las}} = 160.9\,{\text μ} $J.
Laser collisional absorption In addition to the ionisation energy loss, the laser pulse transfers its energy to electrons due to the electronion collisions. As it is a weaker effect, we account for it as an addition to ionisation using standard expressions for the electron collision frequency $ \nu_e $ and Coulomb logarithm $ \Lambda_e $^{33}:
$$ \begin{align} {\nu_e} & {= 2.9\times 10^{6} Z n_e \Lambda_e T_e^{3/2}\,{\rm{s}}^{1} } \\ {\qquad \Lambda_e} &{ =24\ln(n_e^{1/2}/T_e)} \end{align} $$ (4) where the $ n_e=Zn_{a} $ is the electron density in cm^{3} and electron temperature is in eV. For estimates of the collision frequency during the laser pulse, we account for the electron quiver velocity being larger than the electron thermal velocity. According to Ref. 34, under such conditions, the electronion collision frequency depends on the electron ponderomotive energy instead of the electron temperature. Conversely, after the end of the laser pulse, the electron energy is strongly reduced, and the collision frequency increases. The electron distribution is isotropized on a time scale of less than 1 ps. Since that time is shorter than the time of THz emission, we assume that electron energy is equally distributed in three directions and take $ T_e $ as twothirds of the residual electron energy.
For the prepulse, the electron collision frequency $ \nu_e =0.95 $ ps^{−1}, and the fraction of laser energy lost is 0.4%, which corresponds to the electron energy gain of $ 0.84\,{\text μ} $J. For the main pulse, assuming that air is compressed 5.4 times, the electron collision frequency $ \nu_e=2.5 $ ps^{−1}, the fraction of laser energy lost is 7.2% and the electron energy gain is $ 28.7\,{\text μ} $J.
Combining the ionisation energy, residual electron energy, and collisional heating for the prepulse, we estimate the total energy loss of $ 14.3\,{\text μ} $J, electron density $ 2.2\times 10^{20} $ cm^{−3} and electron temperature of 26.6 eV. Collisional heating contributes about 10% to the electron temperature at the end of the laser pulse. The electron collision frequency is 72.6 ps^{−1}, so isotropisation proceeds very rapidly after the end of the laser pulse.
For the main pulse, the total absorbed laser energy is $ 189.6\,{\text μ} $J, the electron density is $ 1.4\times 10^{21} $ cm^{−3}, and the electron temperature is 64.1 eV. Collisional heating contributes about 23% to the electron temperature at the end of the laser pulse. The electron collision frequency is 147.4 ps^{−1}.
Propagation of the heat wave and shock formation After the laser prepulse, the thermal energy in plasma is $ 8\,{\text μ} $J. It is larger but comparable to the shock energy measured in the experiment. The plasma pressure is 9.3 kbar. It leads to plasma expansion and the formation of a cylindrical shock on a time scale of a few hundred ps. According to the measurements^{27}, the shock propagates a distance of $ 33\,{\text μ} $m in 10 ns, and the air compression in the shock front is $ \rho/\rho_0\simeq 5.4 $, which is slightly less than the maximum compression $ (\gamma+1)/(\gamma1)=6 $ for the air adiabat $ \gamma=1.4 $.
The parameter that defines the shock formation is $ \kappa_{e0}/c_{s0}w_0 $ where $ \kappa_{e0}= 4500 $ cm^{2}/s is the electron heat conductivity, and $ c_{s0}= 13 $ km/s is the sound velocity. (See details of derivation in Appendix. "Propagation of the heat wave and shock formation". The shock formation time is 650 ps, and the initial shock position is $ 8.7\,{\text μ} $m. The initial shock velocity is $ c_{s0} (\kappa_{e0}/c_{s0}w_0)^{1/5}=7 $ km/s, which corresponds to the Mach number $ M\simeq 20 $. These values are the initial conditions for the cylindrical shock that can be described by the TaylorSedov formula: $ r(t) = r_s (t/t_s)^{1/2} $. The Mach number decreases as a square root of time.
The thickness of the shock front is needed for further analysis of mechanisms of electromagnetic emission. A general theory of shock waves predicts the front width on the order of a few molecule mean free paths^{29}. A more detailed analysis^{35} of the shock in air evaluates the front thickness as $ 2.53 $ mean free paths for the shock Mach number $ M\gtrsim 3 $. The mean free path of a molecule in air at normal conditions is $ 0.066\,{\text μ} $m. Consequently, we take the shock front thickness $ d_{\rm{sh}} = 0.2\,{\text μ} $m as a reference value.
The authors of Ref. 36 have reported on the studies of shock excitation in air under conditions similar to our experiment. They used a 0.5 PW/cm^{2}, 360 fs pulses at 1030 nm wavelength focused by $ NA = 0.4 $ objective lens and showed initiation of air breakdown along $ \sim 10\; {\text μ} $m long focus. A shock wave formation started 0.1 ns after the shot. The experimentally determined threshold intensity of the air breakdown was 80 TW/cm^{2}. Experimental results of shock wave propagation are well fitted by the equation of state of air $ P=(\gamma1){\cal{E}}_{in} \rho/\rho_0 $, where $ {\cal{E}}_{\rm{in}} $ is the internal energy per unit volume, which is driving explosion. Initial pressures reaching kbar level were reduced to tens of bar after propagation of a few micrometers; the Mach number $ \sim30 $ and velocity of shock wave were determined from a side view imaging. This paper also shows that a tighter focusing can localise the shocked region to small volumes within crosssections of tensofμm^{36} and are promising for subwavelength sources of THz radiation. These conditions are similar to our study.
Diffusion currents Since the electromagnetic emission is produced on a time scale larger than the electron collision time, the electron inertia can be neglected, and Ohm's law describes the current excitation after the end of the laser pulse:
$$ {\bf{E}}_s = \eta_e {\bf{j}}  \frac{1}{e n_e} \nabla (n_e T_e) $$ (5) where the first term on the righthand side describes the electric resistivity due to electron collisions, $ \eta_e=m_e \nu_e/e^2 n_e $, and the second term represents the source related to the gradient of the electron pressure $ n_eT_e $.
The main laser pulse propagates parallel to the shock front, along the zaxis, and it is perpendicular to the direction of shock wave propagation, which is along the $ x $axis, see Fig. 6. The density profile in the shock front $ n_{e0}(x) $ is asymmetric: it has a fast rise at the front side and a slow decrease at the rear side.
Fig. 6 Twopulse enhancement of THz generation in air by the spontaneous magnetic field $ B_{\textit z} $ due to coupled nonparallel gradients of $ \nabla T_e $ and $ \nabla n_e $.
The electron temperature is also inhomogeneous; it is larger on the laser axis and decreases with the radius $ r =\sqrt{x^2 +y^2} $ in the plane $ (x,y) $ perpendicular to the laser propagation direction. Since the electron density and temperature gradients are nonparallel, the electric current has two components: an electrostatic component in the direction of the pressure gradient and a rotational component, which generates a magnetic field perpendicular to the density and temperature gradient. These currents are considered in the next two sections.
Electrostatic diffusion current. The electrostatic current is associated with the charge density perturbation described by the continuity equation $ e\partial_t n_e=\nabla\cdot {\bf{j}} $. Inserting in this equation expression for the current from Eq. 5 and using the Poisson equation for the electric field, $ \epsilon_0 \nabla\cdot {\bf{E}}=e(n_i n_e) $, we obtain the following equation for the electron density:
$$ \partial_t n_e D_e \nabla^2 n_e=\frac{\omega_{\rm{pe}}^2}{ \nu_e} (Z_e n_in_e)$$ (6) where $ D_e =T_e/m_e \nu_e $ is the electron diffusion coefficient and $ \omega_{\rm{pe}} $ is the electron plasma frequency. This equation describes the generation of an electrostatic field due to the diffusion of electrons heated by the laser pulse, which is similar to the photoDember effect in semiconductors^{37, 38}. Before the laser pulse arrival, the air is neutral. The laser pulse creates free electrons and heats them. Under the pressure gradient, the electrons start moving, while the ions are left behind because of their large mass. Therefore, the electron mobility creates an electric current, which induces the charge separation and electric field. Eventually, the electron diffusion stops when the pressure gradient is compensated by the electrostatic field.
We neglected in Eq. 6 excitation of plasma oscillations assuming that the time of electron heating equal to the laser pulse duration (about 30 fs) is much longer than the plasma wave period, $ 2\pi/\omega_{\rm{pe}} $, which is about 3 fs. The value of the charge separation, $ \delta n_e=n_e Z_en_i $, is estimated by comparing the second term on the lefthand side, $ D_e n_{e0}/d_{\rm{sh}}^2 $ with the term on the righthand side. The value of charge separation is $ \delta n_e/n_{e0} \simeq \lambda_{\rm{De}}^2/d_{\rm{sh}}^2 $, where $ \lambda_{\rm{De}}= (\epsilon_0 T_e/e^2n_e)^{1/2} $ is the electron Debye length. Since the Debye length in the experiment is on the order of 1 nm, about two orders of magnitude smaller than the shock front thickness, the charge separation of very small, about four orders of magnitude smaller than the plasma density. Therefore, Eq. 6 can be linearised and written for $ \delta n_e=n_e Z_en_i $:
$$ \partial_t\delta n_e =\frac{\omega_{\rm{pe}}^2}{ \nu_e}\, \delta n_e+ D_e\nabla^2 n_{e0} $$ (7) The second term on the righthand side drives the charge separation because the diffusion coefficient increases with time during the laser heating. The first term on the righthand side accounts for the charge relaxation due to collisions. This time of the charge relaxation can be estimated by comparing the term on the lefthand side of Eq. 7 with the first term on the righthand side. This time is much shorter than 1 ps. So, the time derivative on the lefthand side can be neglected; that is, the charge separation is almost instantaneously accommodated to the laser heating and $ \delta n_e = \lambda_{\rm{De}}^2 \nabla^2n_{e0} $. The strength of the saturated electric field is $ E_s \simeq T_e/e d_{\rm{sh}} \simeq 3 $ MV/cm is quite large, but no current is created.
This analysis confirms that electron diffusion cannot produce electromagnetic emission in the experiment.
Rotational diffusion current. The mechanism of generation of the rotational component of electric current is illustrated in Fig. 6. It is called the Biermann battery effect^{39, 40} in plasma physics. The crossed gradients of temperature and density induce a rotational electric current, which is not related to the charge separation and evolves on a longer timescale. The current is related to the magnetic field by Ampere's equation, $ {\text μ}_0{\bf{j}}=\nabla\times{\bf{B}}_s $. The equation for the magnetic field is obtained by taking the curl of Ohm’s Eq. 5 and using Faraday's equation for the curl of the electric field:
$$ \partial_t B_s \frac{\eta_e}{{\text μ}_0}\nabla_\perp^2 B_s =\frac1{e n_e} \nabla T_e \times \nabla n_e $$ (8) (We assume here a constant resistivity.) Since the gradients of the temperature and density on the righthand side are in the $ (x,y) $ plane, the diffusion proceeds in that plane, and only zcomponent of the magnetic field is generated.
The characteristic time of magnetic field generation, $ t_{\rm{dif}}=w_0^2{\text μ}_0/\eta_e $, is found by comparing the two terms in the lefthand side of this equation. The value of the magnetic field after the saturation is estimated by comparing the second term on the lefthand side with the term on the righthand side: $ eB_0=T_e {\text μ}_0 w_0/d_{\rm{sh}}\eta_e $, and the estimate of the current follows from the Ampere equation: $ j_x \simeq T_e/e d_{\rm{sh}}\eta_e $.
The magnetic field is approximated in the polar coordinates $ (r,\varphi) $ as $ B_s= B_0 \sin\varphi\, b(t,r) $ where $ b(r,t) $ describes the radial profile of magnetic field:
$$ t_{\rm{dif}}\, \partial_t b w_0^2 \left(b^{\prime\prime}+ \frac1r\, b'  \frac b{r^2} \right) = \frac{r}{w_0} {\rm{e}}^{r^2/w_0^2} $$ (9) This equation is obtained from Eq. 8 by multiplying it by $ \sin \varphi $ and integrating over $ \varphi $ from 0 to $ 2{\text π} $, assuming that $ w_0 \gg d_{\rm{sh}} $. Eq. 9 is solved numerically using Wolfram Mathematica package with the initial condition $ B_s(t=0)=0 $ by considering a Gaussian distribution for the temperature, $ T_e(r)=T_{e0}\exp(r^2/w_0^2) $, and an exponential density profile, $ n_e(x) =2n_{e0}/[1+\exp(x/d_{\rm{sh}})] $. The coordinate origin corresponds to the axis of the main pulse. The left panel in Fig. 7 shows the dependence of the dimensionless function $ b(r,t) $ on the radius and time. The magnetic field increases linearly with time for $ t\lesssim t_{\rm{dif}} $, and then it is saturated at the level $ b_{\max}\simeq 0.14 $. The radial profile has a maximum at $ r\simeq w_0 $. The spatial distribution of the magnetic field in the plane $ (x,y) $ perpendicular to the laser propagation direction is shown in the right panel in Fig. 7. It has a dipole structure with a positive lobe for $ y>0 $ and a negative lobe for $ y<0 $. So, the average magnetic field is zero. The lines of constant magnetic field correspond to the current streamlines. The current is directed along the $ x $ direction near the axis, and the current loops are closed. An approximate analytical expression for the magnetic field is given in Appendix. “Approximate solution for the magnetic field generated at the shock front”.
Fig. 7 a Time and radial dependence of the dimensionless function $ b(r,t) $ obtained from a numerical solution of Eq. 9. b Distribution of the magnetic field in the plane $ (x,y) $ obtained from the solution of Eq. 8 at time $ t=t_{\rm{dif}} $. Arrows show the direction of the electric current. The origin $ x=0 $ and $ y=0 $ correspond to the laser beam axis and the maximum of the plasma density gradient.
For the parameters of the main pulse and for the air compression $ C_r=5.4 $, we find the diffusion time $ t_{\rm{dif}} = 4.7 $ ps, the maximum value of the magnetic field $ B_{\max}=b_{\max}B_0=54 $ T and the maximum current density $ j_{\max}=1.2 $ GA/cm^{2}. The lifetime of that current and magnetic field is limited by heat diffusion. It is evaluated in the preceding section to be on the order of a few hundred ps. Performing integration over the focal volume for $ t\gg t_{\rm{dif}} $, we obtain the energy of the magnetic field after the saturation: $ W_B\simeq B_{\max}^2 V_{\rm{las}}/2{\text μ}_0 $. For our parameters, the total magnetic energy is about $ 6.6\,{\text μ} $J. This sets the upper limit of the emitted electromagnetic energy.
Electromagnetic emission Electromagnetic emission can be deduced from the general formula for the vector potential induced by the timedependent current,
$$ {\bf{A}}({\bf{r}},t) = \frac{{\text μ}_0}{4\pi} \int d{\bf{r}}' \,\frac{{\bf{j}}({\bf{r}}',t')}{{\bf{r}}{\bf{r}}'} $$ (10) where $ t'= t{\bf{r}}{\bf{r}}'/c $ is the delayed time. Considering the radiated field in the far field zone, $ {\bf{r}}{\bf{r}}'\approx r  {\bf{n}}\cdot{\bf{r}}' $, where $ {\bf{n}} $ is the unit vector in the emission direction, and performing the Fourier transform in time, we have then the following expression for the vector potential
$$ {\bf{A}}({\bf{r}}) = \frac{{\text μ}_0}{4\pi r}{\rm{e}}^{ikr} \int d{\bf{r}}' \,{\bf{j}}_\omega ({\bf{r}}')\,{\rm{e}}^{ik{\bf{n}}\cdot{\bf{r}}'} $$ (11) where $ {\bf{j}}_\omega $ is the Fourier component of the current and $ k=\omega/c $ is the wave number of emission. Since the emission zone is smaller than the wavelength of emission, $ kz_0 \ll 1 $, we develop the exponent under the integral in Taylor series, $ {\rm{e}}^\xi=\sum\nolimits_j \xi^j/j! $, and the terms in this series correspond to multipole expansion.
The terms $ j=0 $ and $ j=1 $ do not contribute to the radiation because $ \int d{\bf{r}}'\, {\bf{j}}({\bf{r}}')=0 $ and also $ \int d{\bf{r}}'\, {\bf{j}}({\bf{r}}') \,{\bf{n}}\cdot{\bf{r}}'=0 $. The first relation is evident since the current is a derivative of magnetic field: $ {\text μ}_0 j_x=\partial_y B_s $ and $ {\text μ}_0 j_y= \partial_x B_s $. The second relation is obtained by integrating by parts the lefthand side and knowing that the integral of the magnetic field over $ y $ is zero. Thus, a nonzero contribution comes from the quadrupole term $ j=2 $, which is written as:
$$ {\bf{A}}({\bf{r}}) = \frac{k^2 n_y}{4\pi r}{\rm{e}}^{ikr} {\bf{n}}\times {\bf{z}} \int d{\bf{r}}' y' B_\omega(x',y') $$ (12) where $ {\bf{z}} $ is the unit vector in the laser propagation direction, and $ B_\omega $ is the Fourier component of magnetic field $ B_s $. The integral over the magnetic field at the source can be performed using the analytic approximation for the magnetic field Eq. S.6. Taking then a curl of the vector potential, we have the magnetic field in the far field zone:
$$ {\bf{B}}(r,{\bf{n}}) =\frac{i}{r} {\rm{e}}^{ikr} b_\omega B_{\max} z_R n_y {\bf{n}}\times({\bf{n}}\times {\bf{z}})\, (kw_0)^3 $$ (13) where $ b_\omega= t_{\rm{dif}}/(2+i\omega t_{\rm{dif}}) $ is the Fourier component of the timedependent factor in Eq. S.6. The last factor accounts for the suppression of the emission because of a small ratio of the emission zone to the emission wavelength, $ kw_0 \ll 1 $.
Knowing the magnetic field in the far field, we calculate the electric field as $ {\bf{E}}=c {\bf{n}}\times{\bf{B}} $ and obtain the emission power in the solid angle $ d\Omega $:
$$ \frac{dP_{\rm{rad}}}{d\Omega} = \frac{r^2}{{\text μ}_0}{\bf{n}}\cdot({\bf{E}}\times{\bf{B}}) $$ Integrating the emission power over time, we obtain the total emitted energy. Using the theorem of Parceval, $ \int f^2(t)\,dt= (2{\text π})^{1} \int f_\omega^2d\omega $ we replace the integral over time with the integral over frequency and obtain the spectrum of emission:
$$ \frac{dW_{\rm{rad}}}{d\Omega d\omega} =\frac{c}{\pi{\text μ}_0} B_{\max}^2 z_R^2 b_\omega^2 (kw_0)^6 n_y^2(1n_z^2) $$ (14) In the spherical coordinate system with the polar axis along the z direction and azimuthal angle $ \varphi $ in the $ (x,y) $ plane with respect to the $ x $ axis, we have $ n_{\textit z}=\cos\theta $ and $ n_y=\sin\theta \sin\varphi $, so the angular dependence of the emission power is
$$ dP_{\rm{rad}}/d\Omega \propto \sin^4\theta \sin^2\varphi $$ Since both electric and magnetic fields are proportional to $ n_y $, the maximum of emission is in the $ y $ direction. The polarisation is in the $ {\bf{n}}\times{\bf{z}} $ direction; that is, the polarisation is perpendicular to the observation direction in the plane defined by the observation direction and the $ x $ axis. The function under the integral represents the emission spectrum:
$$ \begin{align} {S(\omega)} & {= b_\omega^2 (kw_0)^6 }\\ & {=t_{\rm{dif}}^2 (\omega w_0/c)^6/(4+(\omega t_{\rm{dif}})^2)} \end{align} $$ (15) Divergence of the spectrum at high frequencies is due to the discontinuity in the time dependence of the magnetic field $ B_s $ at time $ t=0 $ corresponding to the shock ionisation by laser and electron heating. There are several options to introduce the switchon time of the magnetic field generation: the laser pulse duration, the electron collision time, and the time of laser pulse propagation along the focal volume. Among the three options listed above, the longest one is the latter: $ {\textit z}_R/c\simeq 0.2 $ ps. Since $ ct_{\rm{dif}}\gg {\textit z}_R $ the integral over the spectrum is dominated by the cutoff frequency $ \omega_{\rm{cut}}\simeq c/{\textit z}_R $, corresponding to the cutoff wavelength of $ 60\,{\text μ} $m and frequency of 5 THz.
Integrating the emission power over the angles and spectrum, we obtain the total emitted energy:
$$ \begin{align} {W_{\rm{rad}}} & {=\frac{32c}{15{\text μ}_0} B_{\max}^2 {\textit z}_R^2 \int_0^{\omega_{\rm{cut}}} d\omega b_\omega^2 (\omega w_0/c)^6} \\ & {\simeq \frac{32}{15{\text π}} W_B \left(\frac{w_0}{{\textit z}_R}\right)^4 } \end{align} $$ (16) Comparing the radiated energy $ W_{\rm{rad}} $ with the energy stored in the magnetic field in the focal zone, $ W_B=B_{\max}^2 V_{\rm{las}}/2{\text μ}_0 $, the factor $ (w_0/{\textit z}_R)^4= ({\rm{NA}}/0.61{\text π})^4=1.8\times 10^{5} $ defines the emission efficiency. For the magnetic field energy of $ 6.6\,{\text μ} $J, the radiated energy is 81 pJ, and the radiation power is about 400 W. That estimation is about four times larger than the measured emitted energy of 22 pJ, which can be considered as a reasonable agreement.