Fine 3D Control of THz Emission in Air with Dual Femtosecond Laser Pre-pulses at Tunnelling Ionisation Regime

1. Introduction

High-intensity ultra-short laser pulses become inseparable from the future development of industrial applications in 3D nano-/micro-machining/printing due to Moore’s law-like increase of their average power over the last years from 2000 [1]. Deterministic energy deposition via highly nonlinear multi-photon and avalanche ionisation is key to 3D CNC machining metals and dielectrics from TW/cm² to PW/cm² intensity window [2]. The laser ablation thresholds in the case of ultra-short 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₂, the ablation threshold was 1.3 J/cm², i.e., 1.5 times smaller than for longer fs-pulses, and the most efficient rate of ablated volume per pulse was at 4 J/cm² fluence or 0.57 PW/cm² intensity for 7 fs pulses [3]. Formation of an over-critical (reflective) plasma mirror at the central part of the Gaussian intensity profile defined energy deposition efficiency and high precision ablation depth control with tens-of-nm per pulse. The ablated volume per energy saturates at the 30 fs pulse duration. Plasma mirror is proposed for compression of ultra-short laser pulses to reach the exawatt intensity scale [5]; 1 EW = ${10}^{18}$ W.

A 3D printing/polymerisation at ∼TW/cm² becomes ionisation controlled, and absorption defines polymerisation rather than the chemical composition of photo-initiators and two-photon absorbers. Deterministic energy deposition occurs at higher intensities corresponding to the tunnelling ionisation of molecules and atoms at $\sim 10$ PW/cm² [6]. In ambient air, focusing ultra-short pulses of a few tens of fs into spots of a few micrometres 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 sub-wavelength for emitted THz radiation.

THz spectral range $0.1-30$ THz (3 mm - 10 $\mu$m in wavelength) has fundamental importance for material science due to characteristic lattice vibrations of phonon modes. Among established semiconductor-based THz emitters/receivers working on photo-activated optically gated currents (Auston switch), a new emerging direction is using nano-thin topological insulators [8]. Another emerging direction is the use of dielectric breakdown as a tool to create short-lived ($\sim 1$ ps) and micro-localized (sub-wavelength at THz frequencies) currents for THz emission [9]. Coherent directional emitters at far-IR wavelengths can be realised using Wolf’s effect, where surface sub-wavelength structures of period $\Lambda \approx 10\mu$m for far-IR black body radiation ${\lambda }_{\mathrm{black}}$ are defined to extract emission of the surface phonon polariton (SPP) mode ${\lambda }_{\mathrm{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_{\mathrm{black}}\sin \theta =k_{\mathrm{SPP}}+2\pi /\Lambda$ [10]. In this way, acoustic/optical photon THz frequencies are extracted into free-space radiation, which is, additionally, coherent due to the sub-wavelength character of the grating.

The single-cycle high-intensity of electrical THz field of 1 MV/cm was demonstrated with fs tilted-front pulses using optical rectification in LiNbO₃ crystal [11]. Periodically polled lithium niobate was used to produce THz emission via walk-off of optical pump and THz wave [12]. High THz pulse energy of 10 $\mu$J at a peak field strength of 0.25 MV/cm, was produced by optical rectification when phonon-polariton phase velocity was matched with the group velocity of driving fs-laser pulse and reached high visible-to-THz photon conversion of 45% [13]. Tunable 10-70 THz emitters were realised via difference-frequency mixing of two parametrically amplified pulse trains from a single white-light seed in GaSe or AgGaS₂ crystals and reached 100 MV/cm peak field at ${10}^8$ W peak power (19 $\mu$J energy) using table-top 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. Spin-Hall currents in tri-layers of 5.8-nm-thick W/Co₄₀Fe₄₀B₂₀/Pt were used as THz broadband 1-30 THz emitters generating field strength of several 0.1 MV/cm and outperformed ZnTe $〈110〉$ emitters [15]. An electrical conductivity induced via the Zenner tunneling 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₂ by sub-4 fs pulse of 1.7 eV light, the currents were detected due to induced electron conductivity of 5 ${\Omega }^{-1}·$m⁻¹, which is 18 orders of magnitude higher than the static d.c. conductivity of amorphous silica [16]. Optical control of phonon modes and Zenner tunneling can be operated without the dielectric breakdown maintaining solid-state material. High-intensity ultra-short irradiation of gasses, liquids, and solids at higher intensities, when the dielectric breakdown occurs, can also be used for producing high-intensity THz fields [17]. By co-propagation of the first and second harmonics beams focused with $f=100$ mm lens into long filaments in air plasma, the AC-bias conditions, THz emitters were realised delivering peak intensities up to 10 kV/cm [18,19,20]. Polarisation of the emitted THz radiation can be set circular by introducing non-colinearity between the filaments [21]. Irradiation of He-jet in vacuum was shown to produce coherent THz emission via the transition radiation of ultra-short 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. The field of THz sources has become an active research frontier [23,24,25].

Here, we used a focused fs-laser pulse for air breakdown from sub-wavelength volumes, which can potentially lead to the engineering of high-intensity THz emitters due to the coherent addition of THz electrical fields and ultra-short duration of a few optical periods. Two pre-pulses used in this study revealed a position dependence on the THz emission and its polarisation, extending from the previous study where one pre-pulse was used [26]. Experimental conditions of focusing and THz detection were kept the same [26], and only two pre-pulses were introduced using a holographic beam shaping technique based on a spatial light modulator (LC- SLM) [27]. For consistency, the same pulse energies were used for pre- and main-pulses: 0.2 and 0.4 mJ (13 GW/pulse), respectively [26]. The range of THz emission was $0.1-2.5$ THz [26]. 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 three-pulse 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).

2. Experimental results

Two-pulse experiment with pre- and main pulses [26] revealed the required timing and spatial separation for the shocked-compressed air by the pre-pulse to be pushed into the position of the main pulse at $(x,y)=(0,0)$ for the maximum emission of THz radiation. Pre-pulse had to arrive 9.7 ns earlier to the focus of the main pulse from the position $\Delta y\sim 33\mu$m away for the maximum THz emission. Polarisation of THz radiation was found aligned along the direction of focal positions of pre-pulse and main pulse [26]. Here, we add a second pre-pulse of the same 0.2 mJ energy with the possibility of placing it freely in 3D space $(x,y,z)$. Figure 1 shows the scheme of the experiment and the diagnostics. Details of the experimental setup and the measurement methodology are described in 6. The effect of placement of the second pre-pulse on polarisation and intensity of THz radiation was systematically studied. Due to short 35 fs pulse ($\sim 11\mu$m in length) and z-position (along the beam propagation) change only tens-of-micrometers. Hence, the timing of both pre-pulses can be considered simultaneous, regardless of the actual z-coordinate of the pre-pulse 2 within an error of $\sim 0.1$ nm, which is by an order of magnitude smaller than the shockwave thickness.

2.1. Y-position of pre-pulse 2

Figure 2 shows the emission of THz radiation ($E_y$ component) for the situation when the pre-pulse 1 is focused at $\Delta y=+33\mu$m. At the same time, the second pre-pulse is focused at other $(0,\pm y)$ position with respect to the main-pulse focus at $(x,y)=(0,0)$; see Figure 1. Change of $E_y$ field orientation was observed over the measured range of THz frequencies as shown by the time domain transients in Figure 2(a). When two pre-pulses 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.

For two pre-pulse experiments, a dependence of THz emission on the offset $\Delta y_2$ of the second pre-pulse was determined. The delay of the first pre-pulse 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 peak-to-peak value (see inset in Figure 3). A pictorial diagram for each irradiation point is presented by the projection of THz E-field on $(E_x,E_y)$-plane. The projection shows the overall temporal evolution of THz electric field orientation during emission time.

One pre-pulse produced enhancement of THz emission by 13 times, i.e., the value $E_{\mathrm{peak}-\mathrm{to}-\mathrm{peak}}^2=53\mu$V² compared to the value $E_{\mathrm{peak}-\mathrm{to}-\mathrm{peak}}^2\approx 4\mu$V² without a pre-pulse. The $E_y$ component of the THz field becomes larger for offsets $|\Delta y_2|>33\mu$m compared to the enhancement produced by the first Conversely, when pre-pulse 2 was closer to the position $(0,0)$ where the main pulse is focused, the amplitude of THz field was reduced.

The second pre-pulse 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 pre-pulse 1 with respect to the main pulse focus. A different evolution in time for such 0.05 mJ second pre-pulse was observed. It showed a $\sim 6$ ps duration of larger negative values of $E_x$ component (see inset in Figure 3(c) 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 pre-pulse 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\left(t\right)\propto E_p^{1/4}\sqrt{t}$ according to the Sedov–von Neumann–Taylor model [28]. The major contribution to the THz enhancement was governed by the first pre-pulse with an energy of 0.2 mJ.

Interestingly, when pre-pulse 2 was on the same side as pre-pulse 1 at a distance from the main pulse position larger than 33 $\mu$m (pre-pulse-1), the $E^2$ enhancement was $\sim 15$ times (see Sec. Sec. 3 for discussion). In this case, the shock-affected volume from the pre-pulse 1 was crossed by the pre-pulse 2. When a shock wave propagates through a region of rarefied density ${\rho }_1$ formed by the other pre-pulse, the same cylindrical explosion model predicts scaling $R\propto {\rho }_1^{-1/4}$. When pre-pulse 2 arrived at $y=0$ to $-33\mu$m (the opposite side of the y-axis as compared to pre-pulse 1), it quenched the effect of THz enhancement. However, at even more negative positions of pre-pulse 2, the same enhancement was observed for arriving shocks to the position of the main pulse; see red-dashed lines drawn as eye guides in Figure 3(c).

2.2. X-position of pre-pulse 2

Next, the dependence of THz emission of the x-position of the second pre-pulse was determined, see Figure 4. The $(E_x,E_y)$ trajectories indicate the radial influence of the pre-pulse irradiation position with respect to the main pulse (the pre-pulse 1 was always at the $(0,33)\mu$m). The ellipticity of THz polarisation was larger for the pre-pulse 2 locations (larger $|\Delta x_2|$). This is explained by the optimal timing of 9.7 ns for the pre-pulse position at $\sim 33\mu$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\mu$m, the shock arriving along the x-axis does not fit the perfect timing for THz enhancement and alters the effect of the pre-pulse 1 arriving along a perpendicular direction. When $|x_2|>33\mu$m, pre-pulse 2 makes a small contribution to the enhancement determined by pre-pulse 1. The pre-pulses 1 and 2 arriving at perpendicular directions define the ellipticity of the THz radiation and its handedness (see Figure 4). The duration of the most intense THz radiation lasts approximately one cycle, corresponding to $2\pi$ rotation in the $(x,y)$ plane (see a side-view cross-section along the time axis in the inset of Figure 4).

The dependence of THz emission on the position of the pre-pulse along the Z axis is presented in B.

3. Theoretical description of light-air interaction and THz emission

3.1. Laser air ionization and electron heating

Here, a physical scenario of air breakdown by the pre-pulse 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 C.

Laser focusing in our experiment was carried out with an off-axis parabolic mirror with an effective numerical aperture $NA=0.125$ as determined from direct imaging of water-jet breakdown [29]. The focal spot area is estimated from the Airy disk radius $w_0=0.61\lambda /NA\approx 3.9\mu$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 120-150\mu$m [26]. It is close to the depth-of-focus $2z_R$, where the Rayleigh length $z_R=\pi w_0^2/\lambda \approx 60\mu$m. The focusing volume is considered as a cylinder of length $2z_R$, radius $w_0$ and volume $V_{\mathrm{las}}=5.7\times {10}^{-9}$ cm³.

The energy of the pre-pulse $E_1=0.2$ mJ and the main-pulse $E_m=0.4$ mJ correspond to the intensity $1.2\times {10}^{16}$ W/cm² and $2.4\times {10}^{16}$ W/cm², respectively. Under these conditions, the laser ponderomotive potential is much larger than the ionization 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 ionization rate $w_i\left(E\right)$ as a function of laser electric field $E\left(t\right)$ and atomic charge i is taken according to the Ammosov-Delone-Krainov (ADK) theory [30]. For the sake of simplicity, we assume a consecutive ionization from unperturbed atomic levels.

Electron ionization and heating are described by the kinetic equation for the distribution function of free electrons, $f_e={\sum }_j{f}_{e,j}$, where $f_{e,j}$ is the electron partial distribution function corresponding to the atom ionization from the level j:

$$\begin{array}{c}\hfill {\partial }_t{f}_{e,j}-\frac{e}{m_e}\mathbf{E}·{\partial }_{\mathbf{v}}{f}_{e,j}=(1-n_{e,j})w_j\left(E\right)\delta \left(\mathbf{v}\right),\end{array}$$

(1)

where $f_{e,j}$ is normalized by the density of neutral atoms, $n_a$, $E\left(t\right)$ is the laser electric field, $n_{e,j}=\int dv{f}_{e,j}$ is the probability of ionization from the level j, $\delta \left(\mathbf{v}\right)$ is the delta function and $\mathbf{v}$ is velocity. In the case of linear polarization, non-relativistic electric field, and neglecting electron collisions, the electrons are moving in the direction of the laser field, and the distribution function is one-dimensional.

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 ionization 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 dvv{f}_{e,j}$ and $n_{e,j}{v}_{e,j}^2=\int dv{v}^2{f}_{e,j}$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\partial }_t{n}_{e,j}& =(1-n_{e,j})w_j\left(E\right),{\partial }_t\left(n_{e,j}{v}_{e,j}\right)=-\frac{e}{m_e}{n}_{e,j}E,\hfill \\ & \frac{1}{2}{\partial }_t\left(n_{e,j}{m}_e{v}_{e,j}^2\right)=-e{n}_{e,j}{v}_{e,j}E.\hfill \end{array}\end{array}$$

(2)

These equations are solved numerically for each ionization level and each atom separately, and the average ionization and the average residual electron energy are evaluated as

$$\begin{array}{c}\hfill Z_e=\sum _j{n}_{e,j},{\epsilon }_e=Z_e^{-1}\sum _j\frac{1}{2}{n}_{e,j}{m}_e{v}_{e,j}^2.\end{array}$$

(3)

Two examples of numerical solutions are shown in Figure 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 [31].

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 pre-pulse ${\epsilon }_{\mathrm{ion}}=127.2$ eV, $Z_e=4.06$, and ${\epsilon }_e=35.7$ eV; for main pulse ${\epsilon }_{\mathrm{ion}}=244.8$ eV, $Z_e=4.85$, and ${\epsilon }_e=74.0$ eV. Multiplying the total energy loss per ion, ${\epsilon }_{\mathrm{tot}}={\epsilon }_{\mathrm{ion}}+Z_e{\epsilon }_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 pre-pulse propagating in the air at normal conditions, the electron density is $n_e=Z_e{n}_a=2.2\times {10}^{20}$ cm⁻³ and the energy loss is ${\epsilon }_{\mathrm{tot}}{n}_a{V}_{\mathrm{las}}=13.4\mu$J. For the main pulse, considering propagation along the shock front in air compressed by a factor $C_r=5.4$, the electron density is $n_e=2Z_e{C}_r{n}_{\mathrm{air}}=1.4\times {10}^{21}$ cm⁻³ and the laser energy loss is ${\epsilon }_{\mathrm{tot}}{C}_r{n}_a{V}_{\mathrm{las}}=160.9\mu$J.

3.1.1. Laser collisional absorption

In addition to the ionisation energy loss, the laser pulse transfers its energy to electrons due to the electron-ion 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$ [32]:

$$\begin{array}{c}\hfill {\nu }_e=2.9\times {10}^{-6}Z{n}_e{\Lambda }_e{T}_e^{-3/2}{\mathrm{s}}^{-1},{\Lambda }_e=24-\ln (n_e^{1/2}/T_e),\end{array}$$

(4)

where the $n_e=Z{n}_a$ is the electron density in cm⁻³ and electron temperature is in eV. For estimates of the collision frequency during the laser pulse, we take the electron ponderomotive energy instead of $T_e$; after the end of the laser pulse, the isotropisation is fast, and we take as $T_e$ the two-thirds of the residual electron energy.

For the pre-pulse, the electron collision frequency ${\nu }_e=0.95$ ps⁻¹, and the fraction of laser energy lost is 0.4%, which corresponds to the electron energy gain of $0.84\mu$J. For the main pulse, assuming that air is compressed 5.4 times, the electron collision frequency ${\nu }_e=2.5$ ps⁻¹, the fraction of laser energy lost is 7.2% and the electron energy gain is $28.7\mu$J.

Combining the ionisation energy, residual electron energy, and collisional heating for the pre-pulse, we estimate the total energy loss of $14.3\mu$J, electron density $2.2\times {10}^{20}$ cm⁻³ 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⁻¹, so isotropisation proceeds very rapidly after the end of the laser pulse.

For the main pulse, the total absorbed laser energy is $189.6\mu$J, the electron density is $1.4\times {10}^{21}$ cm⁻³, 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⁻¹.

3.2. Propagation of the heat wave and shock formation

After the laser pre-pulse, the thermal energy in plasma is $8\mu$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 [26], the shock propagates a distance of $33\mu$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)/(\gamma -1)=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²/s is the electron heat conductivity, and $c_{s0}=13$ km/s is the sound velocity. (See details of derivation in D). The shock formation time is 650 ps, and the initial shock position is $8.7\mu$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 Taylor-Sedov formula: $r\left(t\right)=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 [28]. A more detailed analysis [33] of the shock in air evaluates the front thickness as $2.5-3$ 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\mu$m. Consequently, we take the shock front thickness $d_{\mathrm{sh}}=0.2\mu$m as a reference value.

The authors of ref. [34] have reported on the studies of shock excitation in air under conditions similar to our experiment. They used a 0.5 PW/cm², 360 fs pulses at 1030 nm wavelength focused by $NA=0.4$ objective lens and showed initiation of air breakdown along $\sim 10\mu$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². Experimental results of shock wave propagation are well fitted by the equation of state of air $P=(\gamma -1){\mathcal{E}}_{in}\rho /{\rho }_0$, where ${\mathcal{E}}_{\mathrm{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 $\sim 30$ 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 cross-sections of tens-of-$\mu$m [34] and are promising for sub-wavelength sources of THz radiation. These conditions are similar to our study.

3.3. 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:

$$\begin{array}{c}\hfill {\mathbf{E}}_s={\eta }_e\mathbf{j}-\frac{1}{e{n}_e}\nabla \left(n_e{T}_e\right),\end{array}$$

(5)

where the first term on the right-hand 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_e{T}_e$.

The main laser pulse propagates parallel to the shock front (along the z-axis) and perpendicular to the direction of shock wave propagation (the x-axis). The density profile in the shock front $n_{e0}\left(x\right)$ is asymmetric: it has a fast rise at the front side and a slow decrease at the rear side.

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 non-parallel, 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.

3.3.1. Electrostatic diffusion current.

The electrostatic current is associated with the charge density perturbation described by the continuity equation $e{\partial }_t{n}_e=\nabla ·\mathbf{j}$. Inserting in this equation expression for the current from Eq. (5) and using the Poisson equation for the electric field, ${ϵ}_0\nabla ·\mathbf{E}=e(n_i-n_e)$, we obtain the following equation for the electron density:

$$\begin{array}{c}\hfill {\partial }_t{n}_e-D_e{\nabla }^2{n}_e=\frac{{\omega }_{\mathrm{pe}}^2}{{\nu }_e}(Z_e{n}_i-n_e),\end{array}$$

(6)

where $D_e=T_e/m_e{\nu }_e$ is the electron diffusion coefficient and ${\omega }_{\mathrm{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 photo-Dember effect in semiconductors [35,36]. 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 }_{\mathrm{pe}}$, which is about 3 fs. The value of the charge separation, $\delta n_e=n_e-Z_e{n}_i$, is estimated by comparing the second term on the left-hand side, $D_e{n}_{e0}/d_{\mathrm{sh}}^2$ with the term on the right-hand side. The value of charge separation is $\delta n_e/n_{e0}\simeq {\lambda }_{\mathrm{De}}^2/d_{\mathrm{sh}}^2$, where ${\lambda }_{\mathrm{De}}={({ϵ}_0{T}_e/e^2{n}_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 linearized and written for $\delta n_e=n_e-Z_e{n}_i$:

$$\begin{array}{c}\hfill {\partial }_t\delta n_e=-\frac{{\omega }_{\mathrm{pe}}^2}{{\nu }_e}\delta n_e+D_e{\nabla }^2{n}_{e0}.\end{array}$$

(7)

The second term on the right-hand side drives the charge separation because the diffusion coefficient increases with time during the laser heating. The first term on the right-hand side accounts for the charge relaxation due to collisions. This time is much shorter than the laser pulse duration. So, the time derivative on the left-hand side can be neglected; that is, the charge separation is almost instantaneously accommodated to the laser heating and $\delta n_e={\lambda }_{\mathrm{De}}^2{\nabla }^2{n}_{e0}$. The strength of the saturated electric field is $E_s\simeq T_e/e{d}_{\mathrm{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.

3.3.2. Rotational diffusion current.

The mechanism of generation of the rotational component of electric current is illustrated in Figure 6. It is called the Biermann battery effect [37,38] 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 time scale. The current is related to the magnetic field by Ampere’s equation, ${\mu }_0\mathbf{j}=\nabla \times {\mathbf{B}}_s$. The equation for the magnetic field is obtained by taking the curl of Ohm’s equation (5) and using Faraday’s equation for the curl of the electric field:

$$\begin{array}{c}\hfill {\partial }_t{B}_s-\frac{{\eta }_e}{{\mu }_0}{\nabla }_{\perp }^2{B}_s=\frac{1}{e{n}_e}\nabla T_e\times \nabla n_e,\end{array}$$

(8)

(We assume here a constant resistivity.) Since the gradients of the temperature and density on the right-hand side are in the $(x,y)$ plane, the diffusion proceeds in that plane, and only z-component of the magnetic field is generated.

The characteristic time of magnetic field generation, $t_{\mathrm{dif}}=w_0^2{\mu }_0/{\eta }_e$, is found by comparing the two terms in the left-hand side of this equation. The value of the magnetic field after the saturation is estimated by comparing the second term on the left-hand side with the term on the right-hand side: $e{B}_0=T_e{\mu }_0{w}_0/d_{\mathrm{sh}}{\eta }_e$, and the estimate of the current follows from the Ampere equation: $j_x\simeq T_e/e{d}_{\mathrm{sh}}{\eta }_e$.

Equation (8) is solved numerically by considering a Gaussian distribution for the temperature, $T_e\left(r\right)=T_{e0}\exp (-r^2/w_0^2)$, and an exponential density profile, $n_e\left(x\right)=2n_{e0}/[1+\exp (x/d_{\mathrm{sh}})]$ created at $t=0$. The coordinate origin corresponds to the axis of the main pulse. The magnetic field is approximated in the polar coordinates $(r,\phi )$ as $B_s=B_0\sin \phi b(t,r)$ where $b(r,t)$ describes the radial profile of magnetic field:

$$\begin{array}{c}\hfill t_{\mathrm{dif}}{\partial }_{t}b-w_0^2\left(b^{\prime \prime }+\frac{1}{r}{b}^{\prime }-\frac{b}{r^2}\right)=\frac{r}{w_0}{\mathrm{e}}^{-r^2/w_0^2}.\end{array}$$

(9)

The left panel in Figure 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_{\mathrm{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 Figure 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 E.

For the parameters of the main pulse and for the air compression 5.4, we find the diffusion time $t_{\mathrm{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². 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_{\mathrm{dif}}$, we obtain the energy of the magnetic field after the saturation: $W_B\simeq B_{\max }^2{V}_{\mathrm{las}}/2{\mu }_0$. For our parameters, the total magnetic energy is about $6.6\mu$J. This sets the upper limit of the emitted electromagnetic energy.

3.4. Electromagnetic emission

Electromagnetic emission can be deduced from the general formula for the vector potential induced by the time-dependent current,

$$\begin{array}{c}\hfill \mathbf{A}(\mathbf{r},t)=\frac{{\mu }_0}{4\pi }\int d{\mathbf{r}}^{\prime }\frac{\mathbf{j}({\mathbf{r}}^{\prime },t^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|},\end{array}$$

(10)

where $t^{\prime }=t-|\mathbf{r}-{\mathbf{r}}^{\prime }|/c$ is the delayed time. Considering the radiated field in the far field zone, $|\mathbf{r}-{\mathbf{r}}^{\prime }|\approx r-\mathbf{n}·{\mathbf{r}}^{\prime }$, where $\mathbf{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

$$\begin{array}{c}\hfill \mathbf{A}\left(\mathbf{r}\right)=\frac{{\mu }_0}{4\pi r}{\mathrm{e}}^{ikr}\int d{\mathbf{r}}^{\prime }{\mathbf{j}}_{\omega }\left({\mathbf{r}}^{\prime }\right){\mathrm{e}}^{-ik\mathbf{n}·{\mathbf{r}}^{\prime }},\end{array}$$

(11)

where ${\mathbf{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, $k{z}_0\ll 1$, we develop the exponent under the integral in Taylor series, ${\mathrm{e}}^{\xi }={\sum }_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{\mathbf{r}}^{\prime }\mathbf{j}\left({\mathbf{r}}^{\prime }\right)=0$ and also $\int d{\mathbf{r}}^{\prime }\mathbf{j}\left({\mathbf{r}}^{\prime }\right)\mathbf{n}·{\mathbf{r}}^{\prime }=0$. The first relation is evident since the current is a derivative of magnetic field: ${\mu }_0{j}_x={\partial }_y{B}_s$ and ${\mu }_0{j}_y=-{\partial }_x{B}_s$. The second relation is obtained by integrating by parts the left-hand side and knowing that the integral of the magnetic field over y is zero. Thus, a non-zero contribution comes from the quadrupole term $j=2$, which is written as:

$$\begin{array}{c}\hfill \mathbf{A}\left(\mathbf{r}\right)=\frac{k^2{n}_y}{4\pi r}{\mathrm{e}}^{ikr}\mathbf{n}\times \mathbf{z}\int d{\mathbf{r}}^{\prime }{y}^{\prime }{B}_{\omega }(x^{\prime },y^{\prime }),\end{array}$$

(12)

where $\mathbf{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 (A6). Taking then a curl of the vector potential, we have the magnetic field in the far field zone:

$$\begin{array}{c}\hfill \mathbf{B}(r,\mathbf{n})=\frac{i}{r}{\mathrm{e}}^{ikr}{b}_{\omega }{B}_{\max }{z}_R{n}_y\mathbf{n}\times (\mathbf{n}\times \mathbf{z}){\left(k{w}_0\right)}^3,\end{array}$$

(13)

where $b_{\omega }=t_{\mathrm{dif}}/(2+i\omega t_{\mathrm{dif}})$ is the Fourier component of the time-dependent factor in Eq. (A6). The last factor accounts for the suppression of the emission because of a small ratio of the emission zone to the emission wavelength, $k{w}_0\ll 1$.

Knowing the magnetic field in the far field, we calculate the electric field as $\mathbf{E}=-c\mathbf{n}\times \mathbf{B}$ and obtain the emission power in the solid angle $d\Omega$:

$$\frac{d{P}_{\mathrm{rad}}}{d\Omega }=\frac{r^2}{{\mu }_0}\mathbf{n}·(\mathbf{E}\times \mathbf{B}).$$

Integrating the emission power over time, we obtain the total emitted energy. Using the theorem of Parceval, $\int f^2\left(t\right)dt={\left(2\pi \right)}^{-1}\int {\left|f_{\omega }\right|}^{2}d\omega$ we replace the integral over time with the integral over frequency and obtain the spectrum of emission:

$$\begin{array}{c}\hfill \frac{d{W}_{\mathrm{rad}}}{d\Omega d\omega }=\frac{c}{\pi {\mu }_0}{B}_{\max }^2{z}_R^2{\left|b_{\omega }\right|}^2{\left(k{w}_0\right)}^6{n}_y^2(1-n_z^2).\end{array}$$

(14)

In the spherical coordinate system with the polar axis along the z direction and azimuthal angle $\phi$ in the $(x,y)$ plane with respect to the x axis, we have $n_z=\cos \theta$ and $n_y=\sin \theta \sin \phi$, so the angular dependence of the emission power is

$$d{P}_{\mathrm{rad}}/d\Omega \propto {\sin }^4\theta {\sin }^2\phi .$$

Since both electric and magnetic fields are proportional to $n_y$, the maximum of emission is in the y direction. The polarization is in the $\mathbf{n}\times \mathbf{z}$ direction; that is, the polarization 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{array}{c}\hfill S\left(\omega \right)=|b_{\omega }{|}^2{\left(k{w}_0\right)}^6=t_{\mathrm{dif}}^2{(\omega w_0/c)}^6/(4+{\left(\omega t_{\mathrm{dif}}\right)}^2).\end{array}$$

(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 ionization by laser and electron heating. There are several options to introduce the switch-on 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: $z_R/c\simeq 0.2$ ps. Since $c{t}_{\mathrm{dif}}\gg z_R$ the integral over the spectrum is dominated by the cut-off frequency ${\omega }_{\mathrm{cut}}\simeq c/z_R$, corresponding to the cut-off wavelength of $60\mu$m and frequency of 5 THz.

Integrating the emission power over the angles and spectrum, we obtain the total emitted energy:

$$\begin{array}{c}\hfill W_{\mathrm{rad}}=\frac{32c}{15{\mu }_0}{B}_{\max }^2{z}_R^2{\int }_0^{{\omega }_{\mathrm{cut}}}d\omega {\left|b_{\omega }\right|}^2{(\omega w_0/c)}^6\simeq \frac{32}{15\pi }{W}_B{\left(\frac{w_0}{z_R}\right)}^4.\end{array}$$

(16)

Comparing the radiated energy $W_{\mathrm{rad}}$ with the energy stored in the magnetic field in the focal zone, $W_B=B_{\max }^2{V}_{\mathrm{las}}/2{\mu }_0$, the factor ${(w_0/z_R)}^4={(\mathrm{NA}/0.61\pi )}^4=1.8\times {10}^{-5}$ defines the emission efficiency. For the magnetic field energy of $6.6\mu$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.

4. Discussion

The theoretical estimates for the frequency cut-off and the emitted energy are compatible with the observations. The analytical expressions for these parameters show the way of further improvements. The conversion efficiency of the magnetic field into radiation can be increased either by increasing the ratio $w_0/z_R$ with a tighter laser beam focusing or by creating a non-zero average magnetic field in the plasma. This will increase the conversion efficiency by a factor of ${(z_R/w_0)}^2\sim 100$. The structure of the magnetic field given by Eq. (A6) corresponds to the circular laser beam interacting with a planar shock. To produce a non-zero average magnetic field one may consider an elliptic laser beam with the ellipse axis at an angle (not equal to 0 or 90°) with respect to the x axis.

Another possibility of increasing emission consists of increasing the magnetic field energy. According to the analysis of Eq. (8), the generated magnetic field depends effectively on two plasma parameters: electron temperature $T_e$ and shock front width $d_{\mathrm{sh}}$. The magnetic field energy depends on electron temperature with the fifth power. So, even a modest increase in electron temperature by a factor of 1.5 will increase the emission energy by 10 times. This can be achieved by increasing the energy of the main laser pulse.

In conclusion, the theoretical analysis of THz emission in air with two laser pulses confirms the hypothesis proposed by Huang et al. [39] that the emission is produced from a plasma created by the main pulse in the shock front. The theoretical estimates of the plasma and shock parameters are in agreement with the measurements. However, the mechanism of electromagnetic emission is different from the diffusion current in semiconductors [36]. An electrostatic current related to the charge separation cannot be excited in plasma. Instead, a rotation current is created by non-parallel gradients of electron density and temperature. This current induces a magnetic field and produces a magneto-quadrupole electromagnetic emission. The emission wavelength is controlled by the length of the laser focal region and can be varied by changing the numerical aperture.

5. Conclusions and Outlook

By using two pre-pulses with the same timing at focal spots placed at different 3D locations with respect to the focus of the main pulse, it is possible to engineer air density gradients arriving at the main irradiation site. An addition of two perpendicular air density gradients ${\nabla }_y\rho$ and ${\nabla }_x\rho$ allows one to control the ellipticity and intensity of the emitted THz radiation. A single cycle circular polarisation THz radiation can be produced. By placing $d\sim 10\mu$m-diameter focal spots in 3D space with separations of $(3-5)d$ in space and $\sim 10$ ns in time, it is possible to synchronise the arrival of shock wave fronts travelling at high Mach $(3-10)\mu$m/ns speeds in the air for THz emission enhancement and polarisation control at the laser pulse intensities of a fully deterministic tunnelling ionisation ${10}^{16}$ W/cm².

The origin of THz emission is the current induced at the crossed plasma density and temperature gradients, called the Biermann battery effect [38]. This study clarified that the electrostatic current related to the charge separation cannot be excited in plasma. A coalescence of two shocks at the position of the main pulse produces stronger air compression than a single pulse, thus generating a stronger THz emission at a shorter wavelength. Since the electron current is proportional to the electron temperature and inversely proportional to the shock thickness, the way to increase the power of THz emission would be the increase of electron temperature and the strength of the shock.

After the collision of two shocks created by two pre-pulses, the interaction of the main laser pulse with the compressed air creates a plasma with asymmetric density distribution, which affects the intensity and polarisation of the THz emission in a two-fold way. First, the three-dimensional density asymmetry produces the anisotropy of pressure force and permits the current circulation on a picosecond time scale. Second, the shock density gradient is defined by the angle of collision of two shocks created by two pre-pulses. It is not necessarily aligned along the line connecting the pre-pulse and the main pulse, which affects the polarisation of the THz signal. The privileged direction of THz emission is in the plane perpendicular to the laser pulse propagation. This is the theoretical prediction that is verified in the experiment.

An example of a high-density micro-target for the generation of an extreme UV emission at 13.5 nm for the next generation of photolithography is realised by irradiation of 20-$\mu$m diameter Sn micro-spheres by a focused 10 $\mu$m wavelength CO₂ laser [40]. Similarly, micro-jets of Ga produced by laser pulse showed a strong enhancement of X-ray generation [41]. In water, micro-droplets and micro-sized gradients of compressed air serve the same purpose of delivering more material into laser pulse/beam, which results in plasma with higher density and temperature [42,43].

6. Materials and Methods

All the experiments were conducted in air under atmospheric pressure (1 bar) at room temperature (RT; 296 K). The experimental setup is described in detail elsewhere [29]. A pulsed femtosecond laser (35 fs, transform-limited, $\lambda =800$ nm, 1 kHz, Mantis, Legend Elite HE USP, Coherent, Inc.) was used and the output pulses were split into the pre-pulse 1 (linearly-polarized parallel to the y-axis, y-pol., 0.2 mJ), pre-pulse 2 (linearly-polarised parallel to the y-axis, y-pol., 0.2 mJ) and the main pulse (linearly-polarized parallel to the x-axis, x-pol., 0.4 mJ). THz emission is induced by tightly focussing the main pulse with an off-axis parabolic mirror (1-inch diameter, focal length $f=50.8$ mm, 47-098, Edmund Optics) into air with a time-delay, $\Delta t=9.7$ ns, between the main pulse and the pre-pulse arrivals. Under the main pulse focusing condition with effective numerical aperture $NA=0.125$, the laser focus area was estimated to be about $8\mu$m in diameter and $120\mu$m in length along the axial direction (depth-of-focus) [29].

Spatial offsets of the two pre-pulses along the x- and y-axes, $\Delta x$ and $\Delta y$, were made possible with a Liquid Crystal on Silicon-Spatial Light Modulator (LCOS-SLM, X15213, Hamamatsu Photonics). The homemade program used computer-generated holograms (CGHs) that form multiple focal points and display them on the LCOS-SLM. CGHs were optimized by the weighted Gerchberg-Saxton (WGS) algorithm [44] for the uniform intensity of the focal points. The pre-pulse 1 was set stationary at ($\Delta x_1$, $\Delta y_1$, $\Delta z_1$) = ($0\mu$m, $+33\mu$m, $0\mu$m) while the pre-pulse 2 was scanned in three dimensions along the x-, y-, or z-axis using 3D holographic focusing [45].

The detection of the THz wave emission is achieved by the standard electro-optic sampling method (see A), THz time-domain spectroscopy (THz-TDS), in the transmission direction along the z-axis with a $\langle 110\rangle$-oriented ZnTe crystal (1-mm thick, Nippon Mining & Metals Co., Ltd.) [9,29,46]. The measurements of the THz polarisation were carried out with a set of wire grids before the detection $\langle 110\rangle$-oriented ZnTe crystal following the standard procedure [47,48,49]. Both polarisations of the THz emission $(E_x,E_y)$ were determined from the event of laser irradiation by the main pulse via calculation of the measured fields $E_{\pm {45}^{\circ }}$ using $\pm {45}^{\circ }$ orientated wire-grid polarisers in respect to the x-polariser (high transmission) placed just before the $\pm {45}^{\circ }$ polariser. The effective repetition of the laser excitation is at 0.5 kHz to restore the same air conditions for subsequent laser shots. Two laser pulses (main) were required to determine $E_{-45}$ and $E_{+45}$ components separately. Transient of THz emission was established from separate laser shots for the -5 to 15 ps $(E_x,E_y)$-projections.