1. INTRODUCTION

In recent decades, terahertz (THz) science and technology has become a cutting-edge research topic attracting intense interest from all over the world [1]. One of the frequently used methods to generate THz waves in a laboratory environment is by means of femtosecond laser filamentation in air [2–25]. Hence, significant efforts have been made to investigate the effect of various experimental parameters on THz wave generation during filamentation, including the pumping laser pulse wavelength [2,3], duration [4], polarization [5–8], and energy [4,5,9–15]; the species and pressure [10] of the gas interacting with the pumping laser; the relative phase [7–9] between the fundamental laser pulse and its second harmonic in case of two-color laser pumping; and also the frequency doubling beta-barium borate (BBO) crystal’s rotation angle and BBO-to-focus distance [5–7,10–13], and even the length of the laser filament [7,13,16], etc.

However, some of the reported characteristics of THz radiation from femtosecond laser filaments remain controversial. For example, as for the polarization of the THz emission from a single-color laser filament, radial and elliptical polarizations have been suggested by Refs. [17–19], respectively. In the two-color case, there are also different points of view that the THz wave could be considered as linear [20] or elliptically polarized [7,21]. Contradictory reports also lie on the frequency-dependent angular distribution of the emitted THz wave during two-color filamentation. Reference [22] demonstrated a highly directional THz beam with a divergence angle and a central dip, which is wider at higher THz bands. In contrast, an opposite trend that the dip tends to decrease with the increasing THz frequency is shown in Ref. [23].

The diverse experimental observations have made the understanding of the underlying physical mechanism of THz wave generation remain fuzzy. Two major physical models have been often referred to account for the THz emission during filamentation. A microscopic polarization scheme attributes the THz radiation to the free electron drifting current driven by the femtosecond laser field [11], in which process both the plasma generation and the subsequent electron velocity matter [24]. An alternative four-wave mixing (4WM) model [14] based on the perturbation theory has been suggested in order to resolve the difficulty of the photocurrent (PC) model in explaining the unusual polarization properties of the THz emission from the filament [5,7,20,25]. Nevertheless, the intrinsic connection between the microscopic PC scheme and the macroscopic 4WM model remains mysterious. Recently, Andreeva et al. [26] have proved that both 4WM (in neutrals) and PC (in plasma) models contribute to the THz radiation during two-color laser filamentation, which has offered a valuable new clue to solve the longstanding problem of the underlying mechanism of THz wave generation induced by laser-plasma interaction.

In this work, it has been observed that the plasma density ($N_e$) varied with the BBO crystal’s rotation angle $\theta$, which is an aspect to which little attention has been paid so far. By studying the correlated relationship between $N_e$ and THz wave generation during filamentation, the reactions of plasma on THz spectral amplitude and width have been demonstrated. The one-dimensional negative dielectric (1DND) optical waveguide model [27] has been proposed to explain the novel experimental observations of THz-plasma interaction. At last, a three-step model, i.e., 4WM+PC+1DND, has been built to reveal the nature of THz wave generation. According to our results, the combination of 4WM and PC models, which also takes into account of 1DND model as the post-process, could interpret the characteristics of the THz pulses generated in both orthogonal polarization directions.

2. EXPERIMENTAL RESULTS

A. Two-Dimensional Detection of the THz Pulse Amplitude and Theoretical Fittings

The experimental setup, schematically shown in Fig. 1, adapted the typical THz wave generation and detection scheme by using filamentation. In brief, a 1-kHz, 800-nm, 50-fs full width at half-maximum (FWHM) Ti:sapphire laser pulse with 1 mJ/pulse was focused ($f=175\mathrm{mm}$) and passed a 100-μm-thick type I phase-matched BBO crystal, creating a two-color plasma filament ($\sim 3\mathrm{mm}$ in length) in air. The emitted THz radiation was collected by two identical off-axis parabolic mirrors onto a 1.5-mm-thick ZnTe crystal, performing standard electric-optic sampling (EOS) measurements [20]. The detected THz pulse has the characteristics of a single cycle and peak frequency at $\sim 0.7\mathrm{THz}$. $X–Y–Z$ is our laboratory coordinate system. The $X$ and $Y$ axes are parallel and perpendicular to the pump beam polarization, respectively, and the $Z$ axis corresponds to the laser propagation direction. Two orthogonal components of THz transient, namely, $E_{\mathrm{THz},X}$ and $E_{\mathrm{THz},Y}$, were detected by setting a wire-grid polarizer between the two parabolic mirrors to project $E_{\text{THz}}$ onto the $X$ or $Y$ axis. During measurements, the BBO crystal rotation angle $\theta$ was varied, which was defined as the angle between the $X$ axis and the extraordinary refractive index axis ($e$ axis) of the BBO crystal.

 

Fig. 1. Schematic of the experimental setup and the coordinate system defined in this work. Inset: the recorded (a) $E_{\mathrm{THz},X}$, (b) $E_{\mathrm{THz},Y}$, and the corresponding (c) $E_{\mathrm{THz},o}$, (d) $E_{\mathrm{THz},e}$ as a function of $\theta$.

Download Full Size | PDF

The detected $E_{\mathrm{THz},X}$ and $E_{\mathrm{THz},Y}$ with respect to $\theta$ are shown in Figs. 1(a) and 1(b), respectively. Weak THz emission by the optical rectification in BBO crystal has been carefully subtracted. By further projecting $E_{\mathrm{THz},X}$ and $E_{\mathrm{THz},Y}$ onto the $o$, $e$ axes, $E_{\mathrm{THz},o}$ and $E_{\mathrm{THz},e}$ can be retrieved as $E_{\mathrm{THz},o}=-E_{\mathrm{THz},X}\sin \theta +E_{\mathrm{THz},Y}\cos \theta$ and $E_{\mathrm{THz},e}=E_{\mathrm{THz},X}\cos \theta +E_{\mathrm{THz},Y}\sin \theta$, which are shown in Figs. 1(c) and 1(d), respectively. It can be clearly seen that $E_{\mathrm{THz},o}$ is much larger than $E_{\mathrm{THz},e}$ for almost all $\theta$ values, consistent with Ref. [20], which suggests that the polarization of the THz emission from two-color filament is mainly parallel to the ordinary refractive index axis ($o$ axis) of the BBO crystal. For the sake of simplifying analysis, the peak-to-peak amplitudes of $E_{\mathrm{THz},o}$ and $E_{\mathrm{THz},e}$, namely, $A_{\mathrm{THz},o}$ and $A_{\mathrm{THz},e}$, are displayed in Figs. 2(a) and 2(b) as black circles.

 

Fig. 2. From top to bottom: simulation outcomes (lines) of 4WM, PC, 4WM + PC, and 4WM + PC + 1DND, together with experimental results (black circles) for (a) $A_{\mathrm{THz},o}$ and (b) $A_{\mathrm{THz},e}$.

Download Full Size | PDF

Then, the well-known 4WM and PC models have been applied to estimate $A_{\text{THz}}(\theta )$. The 4WM model describes the THz pulse generation by mixing with a fundamental laser and its second harmonic ($\omega$ and $2\omega$) under the effect of ${\chi }^{(3)}$, namely $E_{\text{THz}}\propto {\chi }^{(3)}{E}_{\omega }{E}_{\omega }{E}_{2\omega }$ [14]. The PC model considers THz wave emission from free electron oscillation in the two-color asymmetric laser field [10,11]. Here, tunneling ionization dominates the free electron generation process rather than multiphoton ionization [10,11,25], since the Keldysh parameter [25] of the fundamental light at 800 nm is less than 1, and the main contribution of the 400 nm wave is to make a mixed asymmetric field (Appendix A.2.2). Compared with thorough numerical simulations based on a complex model, such as unidirectional pulse propagation equation (UPPE), which would lead to more precise quantitative results, simplified models, such as 4WM and PC models, are also highly valued because they could not only efficiently identify the major underlying physical mechanism involved during the THz wave generation by the filament, but also quickly provide physical explanations about the experimental results without going through time-consuming rigorous numerical simulations.

The computed results given by 4WM and PC models are shown in Fig. 2 as red and blue lines, respectively (details in Appendices A.1–A.4). Immediately, one can see that both models failed to fit the experimental results well, mainly for $A_{\mathrm{THz},e}$ in Fig. 2(b). The experimental $A_{\mathrm{THz},e}$ has a weak fluctuation around zero, while the 4WM model predicts all-positive THz amplitude in the direction of $e$ axis (red line). Besides, the simulation outcomes of the PC model (blue line) have obvious discrepancies with experimental $A_{\mathrm{THz},e}$ in the vicinity of $\theta =90°$ and 270°.

Although the above inconsistency has been noticed by different groups [5,7,25], its physical origin still remains secret. Recently, it has been revealed that during two-color filamentation in air, both neutrals and plasma contribute to THz radiation with an energy ratio of $15:85$, which could be explained by the 4WM and PC models, respectively [26]. In this work, since 4WM and PC are both local physical models, which could not give rise to precise quantitative results because of the unknown proportionality constant, the above reported ratio has been adopted as a fitting parameter in order to well reproduce the experimental results. Therefore, following Ref. [26], 4WM and PC simulation results (red and blue lines in Fig. 2) have been added together with amplitude ratio of $0.42:1$ (square root of $15:85$), and then normalized. Better fittings have been achieved as shown as green lines in Fig. 2, further proving the validity of this two-process theory. However, considerable differences between the experimental and simulated $A_{\mathrm{THz},e}$ in the vicinity of 90° and 270° still exist. Although a better fitting could be achieved with the ratio of $4\mathrm{WM}:\mathrm{PC}=2:1$ (see Fig. 13 in Appendix A.5), this fitting parameter has been discarded since it is in contradiction to the fact that the contribution by 4WM should be much weaker than that from the plasma (PC), especially in the low THz band [26]. Inspired by Ref. [26], one will be curious if other physical process(es) could be involved in THz wave generation during filamentation, accounting for the remaining discrepancies. For example, Ref. [28] has pointed out that, in view of practical applications, the current physical models need to include the macroscopic effects such as beam propagation and dispersion in the plasma.

B. Detection of the Plasma Density and THz Beam Confinement Inside the Filament

Note that the variation of plasma density ($N_e$) as a function of $\theta$ is an aspect to which little attention has been paid so far. Therefore, while recording THz waveforms by EOS setup, measurements on the nitrogen ($N_2$) fluorescence signal, which was emitted from the plasma column and is a good indication of $N_e$ [29], have also been carried out with $\theta$ varied from 0° to 180°. In detail, the longitudinal distribution of $N_2$ fluorescence signal was characterized by side imaging of the filament with a $16\times$ microscopy objective and a charge-coupled device (CCD) camera. A blue band-pass filter with high transmissivity at 337 nm was positioned in front of the objective to eliminate the scattered fundamental wavelength.

The maximum values in each CCD image have been displayed in Fig. 3(a) as black open circles, along $\theta$ ranging from 0° to 180°. One can see that the evolution of $N_e$ has a central peak located at $\theta =\sim 90°$. Off center, the $N_e$ value sharply decreases and then increases to sub-maximums at $\theta =0°$ and 180°. Furthermore, $N_e(\theta )$ has been fitted based on the well-known static tunneling ionization model (see Appendix A.2), as shown as the red line in Fig. 3(a). The simulated $N_e$ values are between $0.87\times {10}^{17}$ and $3.25\times {10}^{17}{\mathrm{cm}}^{-3}$, with two maximums at $\theta =90°$ and 270° and two sub-maximums at $\theta =0°$ (360°) and 180°, well reproducing the characteristics of the experimental $N_e(\theta )$.

 

Fig. 3. (a) $N_e$ as a function of $\theta$. Black open circles: experiment; red line: simulation. (b) KE measured THz beam diameter $d$ inside the plasma column and its frequency dependence at three $\theta$ values.

Download Full Size | PDF

After carefully reading Figs. 2(b) and 3(a), it can be seen that the peaks of $N_e$ are recorded at the same $\theta$ value with the locations of the maximum discrepancies between experimental and simulated $A_{\mathrm{THz},e}$. This observation guides us to pay attention to the correlation between $A_{\text{THz}}(\theta )$ and $N_e(\theta )$. A question is naturally raised: could it be possible that the interaction between the THz wave ($A_{\text{THz}}$) and plasma ($N_e$) results in the inconsistences as shown in Fig. 2(b)?

In order to clarify this issue, knife-edge (KE) measurements [30] on the THz beam diameters $d$ inside the plasma column have been performed around $\theta =90°$. The maximum $d$ at each $\theta$ is shown in Fig. 3(b) as blue bars, whose values are much smaller than the THz wavelength ($\sim 430\mathrm{\mu m}$, at peak frequency of 0.7 THz), and closely follow the evolution trend of $N_e$ within the same range of $\theta$ in Fig. 3(a). Moreover, the frequency dependence of $d$ has also been studied. Three representative results recorded at $\theta =90°$, 120°, and 150°, as displayed in Fig. 3(b), have indicated that $d$ decreases with the increasing THz frequency. Although this feature of THz emission has been reported in the far field after the end of the filament [31,32], here we focus on the similar phenomenon inside the plasma filament. The above experimental observations can be explained by the spatial confinement effect of THz wave inside the plasma filament [27], which actually provides the experimental bases for THz wave interaction with the air plasma. Next, the imaged $N_e$ distributions have been studied in detail. Two CCD images of the filaments have been shown in Fig. 4(a) as insets, whose maximums are pointed by arrows ${\mathrm{b}}_1$ and ${\mathrm{b}}_2$ in Fig. 3(a). Along vertical dashed lines, the radial $N_e$ distributions are extracted and shown in Fig. 4(a) as black open circles, together with Gaussian fittings as red lines. From ${\mathrm{b}}_1$ to ${\mathrm{b}}_2$, one can see the peak $N_e$ values increased by 3.7 times (from $0.87\times {10}^{17}$ to $3.25\times {10}^{17}{\mathrm{cm}}^{-3}$), while the FWHM of $N_e(r)$ changed slightly (from 51 to 45 μm). Therefore, in the following, the distributions of $N_e(r)$ have been considered as Gaussian profiles with FWHM of 50 μm, whose peak values come from the simulated results in Fig. 3(a) (red line).

 

Fig. 4. (a) Side CCD images of the plasma filaments [insets ${\mathrm{b}}_1$ and ${\mathrm{b}}_2$, with maximums pointed by arrows in Fig. 3(a)], and the corresponding radial $N_e$ distributions (experiment denoted by black open circles versus Gaussian fitting denoted by red lines). (b) $\text{Re}[{ϵ}_r]$ distribution with the definition of $2a$.

Download Full Size | PDF

C. THz Wave-Plasma Interaction and THz Spectrum Modulation

Then, the radial distribution of the real part of relative dielectric constant ${ϵ}_r$ inside the plasma filament can be calculated with $N_e(r)$ by

Next, following Ref. [27], the radial THz mode field intensity distributions ${|E_r|}^2$ given by the 1DND model have been calculated for 0.35, 0.7, and 1.4 THz in two cases of ${\mathrm{b}}_1$ and ${\mathrm{b}}_2$ (see Appendix A.6). The results are shown in Figs. 5(a)–5(c), respectively, where THz energy was strongly constrained within the plasma region (mostly at $r=a$) [27,34]. For quantification, $\eta$ has been defined as the confined THz amplitude ratio, calculated by the square root of the ratio of THz power within $r<a$ (the confinement region) to that inside the whole cross-section of the plasma filament (see Appendix A.7).

 

Fig. 5. Normalized ${|E_r|}^2$ distributions given by 1DND model in case of $N_e={\mathrm{b}}_1$ and ${\mathrm{b}}_2$ for (a) 0.35, (b) 0.7, and (c) 1.4 THz.

Download Full Size | PDF

The variations of $\eta$ as a function of $\theta$ at 0.35, 0.7, and 1.4 THz have been calculated and shown in Fig. 6(a). Clearly, $\eta$ decreases when $N_e$ is larger (from ${\mathrm{b}}_1$ to ${\mathrm{b}}_2$). On the other hand, $\eta$ at higher THz frequency decreases even more (e.g., ${\eta }_{1.4\mathrm{THz}}$, as blue dotted line).

 

Fig. 6. (a) Variations of the confined THz amplitude ratio $\eta$ along $\theta$ at 0.35, 0.7, and 1.4 THz, respectively. (b) FWHM of THz spectra as a function of $\theta$ (inset) and two representative THz spectra with FWHM pointed by green arrows.

Download Full Size | PDF

In view of the above $\eta$ variations [Fig. 6(a)], it can be concluded that when $N_e$ increases, the decreasing of $\eta$ would very likely result in the decreasing of THz amplitude (due to diffraction) and the narrowing of THz spectral width, mainly at a higher frequency band. This is actually what can be seen in Fig. 6(b), where the FWHMs of THz spectra have been extracted as a function of $\theta$ (0°–90°) and shown as inset. One can see that the FWHM increased from $\theta =0°$ to $\theta =\sim 45°$ (${\mathrm{b}}_1$) and then decreased till $\theta =90°$ (${\mathrm{b}}_2$). Within the same range of $\theta$, the opposite tendency of $N_e$ variation can be found in Fig. 3(a). Moreover, two THz spectra at ${\text{b}}_1$ and ${\text{b}}_2$ (pointed by green arrows) are shown in Fig. 6(b) as a black solid line and red dashed line, respectively. One can see the width of the former is obviously larger than that of the latter, especially towards high THz frequency.

As discussed before, the reaction of plasma density ($N_e$) on the generated THz spectral amplitude and width via 1DND model has been demonstrated. The additional THz amplitude loss with the increasing of $N_e$, i.e., the THz propagation effects inside the filament, cannot be neglected. However, apparently, neither the 4WM nor PC model has covered this issue, which might be the reason for their failing to perfectly interpret the experimental results in Fig. 2(b). Thus, finally, $\eta (\theta )$ at 0.7 THz based on the 1DND model has been used to modify the $4\mathrm{WM}+\mathrm{PC}$ simulations by multiplying the red dashed line in Fig. 6(a) with the green lines in Fig. 2 (please also see Appendix A.7). The calculation results are normalized and shown in Fig. 2 as pink lines. It can be found that not only the distinguished disagreements between experimental and simulated $A_{\mathrm{THz},e}$ at $\theta =90°$ and 270° have been significantly reduced, but also better fittings of $A_{\mathrm{THz},o}$ (especially around $\theta =90°$ and 270°) have been achieved.

3. CONCLUSION

In summary, the current work has been devoted to shedding light on the underlying mechanism of THz wave generation during filamentation by bridging the well-known 4WM and PC models with the THz–plasma interaction. That is, when considering the influence of plasma density $N_e$ on THz spectral amplitude and width simulated by the 1DND model, as well as 4WM and PC models, the generated orthogonal polarization properties of THz wave during two-color filamentation can be well interpreted.

Compared with 3D numerical simulations [35,36] or a single physical model, the current three-step model could be a better balance between realism and simplicity. It has divided the whole dynamic process of THz wave generation into three separate effects, and described each of them with different (but more proper) physical models: (i) 4WM for Kerr nonlinearity of neutral air molecules from the self-focusing laser pulse, (ii) PC for free electron oscillating as the plasma builds up, and (iii) 1DND for macroscopic propagation effect (spatial confinement) of the generated THz wave inside the plasma column. The unique advantage of this multi-process model is that one can get a clear picture of each stage’s performance, which will significantly promote understanding of the nature of THz emission from laser-induced plasma.

It is also worth mentioning that the influence of THz–plasma interaction on the emitted THz wave features still needs to be systemically investigated in the future, e.g., under different focal conditions. These studies will potentially help to resolve the aforementioned experimental contradictions lying in the references in Section 1 (Introduction).

APPENDIX A

1. RECONSTRUCTION OF TWO-COLOR LASER FIELD AT THE FOCUS

When using the static tunneling ionization model and PC model, the knowledge about the temporal evolution of the two-color laser field at the focus is essential in order to describe the THz wave generation through laser–plasma interaction. Due to the birefringence of BBO crystal and dispersion in air, the laser field distortion during the propagation prior to focus has to be taken into account. The two-color laser field at the focus could be split as the following:

(A1)$$\mathit{E}(t)=E_o(t)\hat{\mathit{o}}+E_e(t)\hat{\mathit{e}}=E_{\omega ,o}(t)\hat{\mathit{o}}+E_{\omega ,e}(t)\hat{\mathit{e}}+E_{2\omega ,e}(t)\hat{\mathit{e}},$$

where the subscripts $\omega$ and $2\omega$ denote the fundamental laser (FL) and its second harmonic (SH), respectively, and the subscripts $o$ and $e$ indicate the polarizations parallel to the ordinary and extraordinary refractive index axes of the BBO crystal, respectively (Fig. 7).

 

Fig. 7. Coordinate system defined in this appendix.

Download Full Size | PDF

In the coordinate moving at the group velocity of $E_{\omega ,o}$, the expressions of $E_{\omega ,o}$, $E_{\omega ,e}$, and $E_{2\omega ,e}$ are

(A2)$$\begin{gathered} E_{\omega ,o}(t)=E_{\omega }\exp (-\frac{t^2}{a^2})·\sin \theta ·\cos [\omega t-(n_{\omega ,o}-n_{g,\omega ,o})\frac{\omega }{c}{L}_1-(n_{\omega ,\text{air}}-n_{g,\omega ,\text{air}})\frac{\omega }{c}{L}_2], \\ {E}_{\omega ,e}(t)=E_{\omega }\exp [-\frac{{(t-{\tau }_{\omega })}^2}{a^2}]·\cos \theta ·\cos [\omega (t-{\tau }_{\omega })-(n_{\omega ,e}-n_{g,\omega ,e})\frac{\omega }{c}{L}_1-(n_{\omega ,\text{air}}-n_{g,\omega ,\text{air}})\frac{\omega }{c}{L}_2], \\ {E}_{2\omega ,e}(t)=\beta E_{\omega }^2\exp [-\frac{2{(t-{\tau }_{2\omega })}^2}{a^2}]{\sin }^2\theta ·\cos [2\omega (t-{\tau }_{2\omega })-(n_{2\omega ,e}-n_{g,2\omega ,e})\frac{2\omega }{c}{L}_1-(n_{2\omega ,\text{air}}-n_{g,2\omega ,\text{air}})\frac{2\omega }{c}{L}_2+\phi ], \end{gathered}$$

where $E_{\omega }$ is the peak amplitude of the fundamental laser field, $a$ is half the pulse duration of the amplitude envelop at $1/e$ level connected to the FWHM of the pulse intensity via $a=\mathrm{FWHM}/\sqrt{2\ln 2}$, e.g., $a=42.47\mathrm{fs}$ when $\mathrm{FWHM}=50\mathrm{fs}$. $\beta$ is a proportional factor related to SH conversion efficiency, and $\phi$ is the relative phase between $\omega$ and $2\omega$. $L_1$ and $L_2$ denote the thickness of BBO crystal (100 μm) and the distance from the back surface of BBO to the focus ($\sim 5\mathrm{cm}$). ${\tau }_{\omega }$ and ${\tau }_{2\omega }$ account for the group delays of $E_{\omega ,e}$ and $E_{2\omega ,e}$ in BBO crystal ($L_1$) and air ($L_2$) compared with $E_{\omega ,o}$. ${\tau }_{\omega }$ is given by

(A3)$${\tau }_{\omega }=\frac{L_1}{c}(n_{g,\omega ,e}-n_{g,\omega ,o})+\frac{L_2}{c}(n_{g,\omega ,\text{air}}-n_{g,\omega ,\text{air}})=\frac{L_1}{c}(n_{g,\omega ,e}-n_{g,\omega ,o}),$$

where the group refractive indices of $E_{\omega ,e}$ and $E_{\omega ,o}$ are $n_{g,\omega ,e}$ and $n_{g,\omega ,o}$ in BBO crystal and both $n_{g,\omega ,\text{air}}$ in air. ${\tau }_{2\omega }$ is computed as

(A4)$${\tau }_{2\omega }=\frac{L_1}{c}(n_{g,2\omega ,e}-n_{g,\omega ,o})+\frac{L_2}{c}(n_{g,2\omega ,\text{air}}-n_{g,\omega ,\text{air}}),$$

where the group refractive indices of $E_{2\omega ,e}$ and $E_{\omega ,o}$ are $n_{g,2\omega ,e}$ and $n_{g,\omega ,o}$ in BBO crystal, $n_{g,2\omega ,\text{air}}$ and $n_{g,\omega ,\text{air}}$ in air. In order to derive the refractive indices and group indices in Eqs. (A2)–(A4), the following Sellmeier equations have been implemented in BBO crystal [37],

(A5)$$\{\begin{array}{c}{n}_o^2=2.7405+\frac{0.0184}{{\lambda }^2-0.0179}-0.0155{\lambda }^2\\ n_e^2=2.3730+\frac{0.0128}{{\lambda }^2-0.0156}-0.0044{\lambda }^2\end{array},$$

and in air [38],

(A6)$$(n_{\text{air}}-1)\times {10}^8=8342.1+\frac{2406030{\lambda }^2}{130{\lambda }^2-1}+\frac{15996{\lambda }^2}{38.9{\lambda }^2-1},$$

where $\lambda$ is in μm. Furthermore, for a BBO crystal cut at an angle $\alpha$ (= 29.8° in our experiment), the extraordinary refractive index of a normally incident light is obtained by

(A7)$$\frac{1}{n_e^2(\alpha )}=\frac{{\cos }^2\alpha }{n_o^2}+\frac{{\sin }^2\alpha }{n_e^2}.$$

By substituting all the experimental parameters into Eqs. (A1)–(A4), the temporal laser field distribution at the focus could be retrieved.

2. SIMULATION OF PLASMA DENSITY (${\mathsf{N}}_{\mathsf{e}}$) BASED ON STATIC TUNNELING IONIZATION MODEL

The plasma (free electron) generation dynamic plays a crucial role when dealing with THz wave generation by femtosecond laser filament in air [11,12]. However, the variation of $N_e$ is an aspect to which little attention has been paid so far. For instance, $N_e$ is sensitive to BBO crystal rotation angle. Nevertheless, it is barely systematically investigated in previous reports. In our experiment, tunneling ionization dominates the free electron generation process [11,39]. For a 400-nm laser whose ponderomotive (quiver) energy is low, making the value of Keldysh parameter ${\mathrm{\gamma }}_{400\mathrm{nm}}>1$, the multiphoton ionization model might be more appropriate. However, the main contribution of 400 nm is to make a mixed two-color asymmetric laser field (please see Figs. 10 and 11), which can strip off bound electrons and produce a non-vanishing transverse plasma current. This photocurrent surge, occurring on the timescale of the photo-ionization, can produce an electromagnetic pulse at THz frequencies. The ionization rate $R(t)$ could be commonly described by the static tunneling ionization model following Ref. [25],

(A8)$$R(t)=\frac{a}{A(t)}\exp [-\frac{b}{A(t)}],$$

where $A(t)=|E(t)|/E_a$ is the electric field in atomic units, and

(A9)$$\{\begin{array}{c}a=4{\omega }_a{r}_H^{5/2}\\ b=(2/3)r_H^{3/2}\\ E_a={\kappa }^3{m}^2{e}^5/{\hslash }^4\approx 5.14\times {10}^{11}\mathrm{V}/\mathrm{m}\\ {\omega }_a={\kappa }^{2}m{e}^4/{\hslash }^3\approx 4.13\times {10}^{16}/\mathrm{s}\end{array}.$$

${\omega }_a$ corresponds to the atomic frequency unit; $r_H=U_{\text{ion}}/U_H$ indicates the ionization potential of the gas molecules under consideration relative to that of hydrogen ($U_H=13.6\mathrm{eV}$) and $\kappa ={(4\pi {ϵ}_0)}^{-1}$. In the simulations, the air is considered as being composed of 78% nitrogen and 22% oxygen with $U_{\text{ion},N_2}=15.6\mathrm{eV}$ and $U_{\text{ion},{\mathit{O}}_2}=12.1\mathrm{eV}$. Hence, the electron density $N_e(t)$ is given by

(A10)$$\{\begin{array}{c}d{N}_e(t)=d{N}_e{(t)}_{{\mathrm{N}}_2}+d{N}_e{(t)}_{{\mathrm{O}}_2}\\ d{N}_e{(t)}_{N_2}=R{(t)}_{N_2}[N_{0N_2}-N_e{(t)}_{N_2}]dt\\ d{N}_e{(t)}_{{\mathit{O}}_2}=R{(t)}_{O_2}[N_{0{\mathit{O}}_2}-N_e{(t)}_{{\mathit{O}}_2}]dt\end{array}.$$

2.1. Confirming the Validity of Simulation Parameter Values

As discussed before, a two-color laser field and the free electron generation mechanism have been built up. Experimental plasma density distribution with respect to $\theta$, namely, $N_e(\theta )$, shown in Fig. 3(a) as black open circles, can be reproduced by time integral of $N_e(t)$ in Eq. (A10) at a varied $\theta$. However, three parameters related to the above-mentioned processes, namely, $E_{\omega }$, $\beta$, and $\phi$ need to be determined in order to understand the plasma density variation. $\beta$ has been approximated to be 0.34 according to the measured SH pulse energy after BBO crystal. $\phi$ could be set to be 180°, which accounts for the phase difference between the second harmonic pulse and the fundamental pulse when the phase-matching condition is satisfied inside the BBO crystal [40]. Then $E_{\omega }$ could be determined by fitting the calculation results given by Eq. (A10) with the experimental results shown in Fig. 3(a) (also in Fig. 8(b) as blue open squares).

 

Fig. 8. (a) Simulated $N_e$ distribution along $\theta$, considering $E_{\omega }$ varying from $1\times {10}^{10}$ to $10\times {10}^{10}\mathrm{V}/\mathrm{m}$. (b) Four representative $N_e(\theta )$ curves for $E_{\omega }=1\times {10}^{10}\mathrm{V}/\mathrm{m}$ (bold black line), $E_{\omega }=2.75\times {10}^{10}\mathrm{V}/\mathrm{m}$ (thin red line), $E_{\omega }=5\times {10}^{10}\mathrm{V}/\mathrm{m}$ (gray dashed line), and $E_{\omega }=7\times {10}^{10}\mathrm{V}/\mathrm{m}$ (green dotted line).

Download Full Size | PDF

The obtained calculation results are shown in Fig. 8(a), in which one can see that a variation of $E_{\omega }$ from $1\times {10}^{10}$ to $10\times {10}^{10}\mathrm{V}/\mathrm{m}$ will not change the overall trend of the evolution of $N_e(\theta )$, but only vary the relative ratio between the maximum and minimum of $N_e$. Four representative $N_e(\theta )$ curves have been indicated in Fig. 8(b) for comparison. The best fitting of experimental $N_e(\theta )$ (blue open squares) could be found when $E_{\omega }=2.75\times {10}^{10}\mathrm{V}/\mathrm{m}$, corresponding to a laser intensity of $1\times {10}^{14}\mathrm{W}/{\mathrm{cm}}^2$ at the focus [10], calculated by $I={|E|}^2/(2\eta )$, where $\eta =377\mathrm{\Omega }$ is the optical impedance in air. The simulation outcome at $E_{\omega }=2.75\times {10}^{10}\mathrm{V}/\mathrm{m}$ is depicted in Fig. 3(a) as the red line.

Moreover, measurements on the power of the fundamental light (800 nm) and its second harmonics (400 nm) have also been carried out as a function of $\theta$. The results are shown in Fig. 9 as green open circles and blue open squares, respectively. In the same figure, black and red lines are the corresponding simulated results given by Eq. (A2) using the same $E_{\omega }$, $\beta$, and $\phi$ as above. The fittings are acceptable, further clarifying the validity of the formulas and parameters set for the calculations in Appendices A.1 and A.2.

 

Fig. 9. Experimental and simulated power of 800 nm and 400 nm after the BBO crystal.

Download Full Size | PDF

2.2. Simulated Laser Field, ${\mathsf{N}}_{\mathsf{e}}$ and Their Asymmetries

Here is an example of the symmetry-broken laser field (obtained at $\mathit{\theta }=55°$) by mixing 800-nm and 400-nm lasers, as shown in Fig. 10(a), which is given by Eq. (A2), i.e., $|E_{\text{total}}(t)|=|E_{\omega ,o}(t)\hat{\mathit{o}}+E_{\omega ,e}(t)\hat{\mathit{e}}+E_{2\omega ,e}(t)\hat{\mathit{e}}|$. One can see it is asymmetric across individual optical cycles. Similar laser field distributions can also be obtained when $\theta =125°$, 235°, and 305°.

 

Fig. 10. (a) Example of asymmetric laser field $|E_{\text{total}}(t)|$ obtained at $\theta =55°$. (b) The corresponding ionization result of $\text{d}{N}_e(t)/\text{d}t$.

Download Full Size | PDF

Then, $|E_{\text{total}}(t)|$ is substituted into Eq. (A8) to calculate the corresponding $\mathrm{d}{N}_e(t)/\mathrm{d}t$ given by Eq. (A10). The result is shown in Fig. 10(b). One can clearly see the asymmetry of temporal ionization induced by the asymmetric $|E_{\text{total}}(t)|$. Due to this asymmetric ionization, THz yield in the $o$ direction, i.e., $A_{\mathrm{THz},o}(\theta )$, reached maximums at $\theta =55°$, 125°, 235°, and 305° [Fig. 2(a)].

By contrast, when $\theta =0°$ (360°), 90°, 180°, and 270°, both $|E_{\text{total}}(t)|$ and $\mathrm{d}{N}_e(t)/\mathrm{d}t$ are temporally symmetric as shown in Figs. 11(a) and 11(b), respectively. Thus $A_{\mathrm{THz},o}(\theta )$ reached minimums at these $\theta$ values [Fig. 2(a)].

 

Fig. 11. (a) Example of symmetric laser field $|E_{\text{total}}(t)|$ obtained at $\theta =0°$. (b) The corresponding ionization result of $\text{d}{N}_e(t)/\text{d}t$.

Download Full Size | PDF

One may be curious about the difference between $A_{\mathrm{THz},o}(\theta )$ and $A_{\mathrm{THz},e}(\theta )$ distributions in Figs. 2(a) and 2(b), especially the number of maximums (peaks and troughs). Reminding ourselves that, based on the PC model (in the following Appendix A.3), the THz emission is induced by not only the asymmetry of ionization [$|E_{\text{total}}(t)|$ or $\mathrm{d}{N}_e(t)/\mathrm{d}t$], but also the asymmetry of the electron-driving laser field [$E_o(t)$ or $E_e(t)$]. Here, $E_o(t)=E_{\omega ,o}(t)$ ($o$-component of fundamental light) is always temporally symmetric regardless of $\theta$. Thus $A_{\mathrm{THz},o}(\theta$) depends only on the asymmetry of $\mathrm{d}{N}_e(t)/\mathrm{d}t$ [or $|E_{\text{total}}(t)|$].

However, in the $e$ direction, $E_e(t)=E_{\omega ,e}(t)+E_{2\omega ,e}(t)$ (the combination of $e$-component of fundamental light and its second harmonic), whose asymmetry is also dependent on $\theta$. Thus $A_{\mathrm{THz},e}(\theta )$ is more complex compared with $A_{\mathrm{THz},o}(\theta )$, e.g., with more maximums at about $\theta =45°$, 90°, 135°, 225°, 270°, and 315° as shown in Fig. 2(b).

3. SIMULATION OF THz WAVE GENERATION BASED ON PHOTOCURRENT (PC) MODEL

With the plasma density $N_e(t)$ obtained by Eqs. (A1)–(A10), the electron drifting current $J(t)$ can be computed as $J(t)=-{\int }_{-\infty }^{\infty }e{N}_e(t_0)v(t;t_0)\mathrm{d}{t}_0$, where $v(t;t_0)$ denotes the electron velocity under the laser field given by $v(t;t_0)=-\frac{e}{m}{\int }_{t_0}^{t}E(t^{\prime })\mathrm{d}{t}^{\prime }$. The electromagnetic radiation $E_{\mathrm{EM}}(t)$ is determined by $E_{\mathrm{EM}}(t)\propto \partial J(t)/\partial t$, and the radiation spectrum is given by the Fourier transformation of $E_{\mathrm{EM}}(t)$. The corresponding THz waveform $E_{\text{THz}}(t)$ emitted from the two-color laser field is retrieved from the radiation spectrum by low-frequency filtering ($<2\mathrm{THz}$ in our experimental condition). It is noteworthy that since $E(t)$ can be divided into a pair of orthogonal components, i.e., $E_o(t)$ and $E_e(t)$ [Eq. (A1)], we actually calculated $E_{\mathrm{THz},o}(t)$ and $E_{\mathrm{THz},e}(t)$, whose polarizations are parallel to the $o$ axis and $e$ axis of the BBO crystal, by substituting $E(t^{\prime })$ in $v(t;t_0)=-\frac{e}{m}{\int }_{t_0}^{t}E(t^{\prime })\mathrm{d}{t}^{\prime }$ with $E_o(t^{\prime })$ and $E_e(t^{\prime })$.

4. SIMULATION OF THz WAVE GENERATION BASED ON FOUR-WAVE MIXING MODEL

For the sake of convenience, the BBO crystal coordinate system, i.e., $o–e–Z$, has been used here. Thus, the electric fields of the THz wave, $2\omega$ and $\omega$, are denoted in the new coordinate system by $E_{\text{THz}}$, $E_{2\omega }$, and $E_{\omega }$, respectively. These coordinates are sketched in Fig. 7, where $E_{2\omega ,e}\propto {(E_{\omega }\sin \theta )}^2$, $E_{\omega ,e}\propto E_{\omega }\cos \theta$, and $E_{\omega ,o}\propto E_{\omega }\sin \theta$. The 4WM model gives

(A11)$$E_{\mathrm{THz},i}\propto \sum _{jkl}{\chi }_{ijkl}^{(3)}{E}_{2\omega ,j}{E}_{\omega ,k}{E}_{\omega ,l}(i,j,k,l=o,e).$$

Recalling that the electric field of the input pump beam $\omega$ will have two complex components, $E_{\omega ,o}$ and $E_{\omega ,e}$, with relative phase shift induced by the birefringence of the crystal. This phase term has been dropped since it is a constant in principle.

Besides, in Eq. (A11), the relative phase between the two colors has been set always as a constant of 0, for the reason as follows. The variation of the relative phase between the two colors along the propagation distance is mainly associated with the chromatic dispersion. By taking into account the refractive indexes of two colors, the corresponding coherence length (dephasing length) is about 25 mm [23], which is longer than the filament length ($\sim 3\mathrm{mm}$) in this work. In addition, according to Aközbek et al.’s work [41], the actual coherence length would be much longer than the one estimated in the linear case. Therefore, we consider it is reasonable to ignore the phase delay variation between two colors in our estimation.

On the other hand, the polarization of the $2\omega$ wave is only along the $e$ axis of the BBO crystal, i.e., $E_{2\omega ,e}$. Therefore, consider the simplified form of Eq. (A11),

(A12)$$E_{\mathrm{THz},i}\propto \sum _{kl}{\chi }_{iekl}^{(3)}{E}_{2\omega ,e}{E}_{\omega ,k}{E}_{\omega ,l}(i,k,l=o,e).$$

Thus, $E_{\mathrm{THz},i}$ has three components related with ${\chi }_{ieeo}+{\chi }_{ieoe}$, ${\chi }_{ieee}$, and ${\chi }_{ieoo}$ as follows:

(A13)$$\begin{gathered} E_{\mathrm{THz},i}\propto ({\chi }_{ieeo}^{(3)}+{\chi }_{ieoe}^{(3)})E_{2\omega ,e}{E}_{\omega ,e}{E}_{\omega ,o} \\ +{\chi }_{ieee}^{(3)}{E}_{2\omega ,e}{E}_{\omega ,e}{E}_{\omega ,e} \\ +{\chi }_{ieoo}^{(3)}{E}_{2\omega ,e}{E}_{\omega ,o}{E}_{\omega ,o} \\ =({\chi }_{ieeo}^{(3)}+{\chi }_{ieoe}^{(3)})E_{\omega }^4{\sin }^3\theta \cos \theta \\ +{\chi }_{ieee}^{(3)}{E}_{\omega }^4{\sin }^2\theta {\cos }^2\theta \\ +{\chi }_{ieoo}^{(3)}{E}_{\omega }^4{\sin }^4\theta (i=o,e). \end{gathered}$$

The calculation results of the above three components of $E_{\mathrm{THz},i}$ (without ${\chi }^{(3)}$) are shown in Fig. 12.

 

Fig. 12. Calculation results of $E_{2\omega ,e}{E}_{\omega ,k}{E}_{\omega ,l}/E_{\omega }^4$.

Download Full Size | PDF

Compared with the variation tendency of the experimental $A_{\mathrm{THz},o}$ and $A_{\mathrm{THz},e}$ (black solid circles in Fig. 2), we finally get

(A14)$$\begin{gathered} E_{\mathrm{THz},o}\propto ({\chi }_{oeeo}^{(3)}+{\chi }_{oeoe}^{(3)})E_{2\omega ,e}{E}_{\omega ,e}{E}_{\omega ,o} \\ =({\chi }_{oeeo}^{(3)}+{\chi }_{oeoe}^{(3)})E_{\omega }^4{\sin }^3\theta \cos \theta \\ \Rightarrow A_{\mathrm{THz},o}(\theta )\propto ({\chi }_{oeeo}^{(3)}+{\chi }_{oeoe}^{(3)}){\sin }^3\theta \cos \theta , \end{gathered}$$

(A15)$$\begin{gathered} E_{\mathrm{THz},e}\propto {\chi }_{eeee}^{(3)}{E}_{2\omega ,e}{E}_{\omega ,e}{E}_{\omega ,e}+{\chi }_{eeoo}^{(3)}{E}_{2\omega ,e}{E}_{\omega ,o}{E}_{\omega ,o} \\ \approx {\chi }_{eeee}^{(3)}{E}_{2\omega ,e}{E}_{\omega ,e}{E}_{\omega ,e}({\chi }_{eeee}^{(3)}\gg {\chi }_{eeoo}^{(3)}[9]) \\ ={\chi }_{eeee}^{(3)}{E}_{\omega }^4{\sin }^2\theta {\cos }^2\theta \\ \Rightarrow A_{\mathrm{THz},e}(\theta )\propto {\chi }_{eeee}^{(3)}{\sin }^2\theta {\cos }^2\theta . \end{gathered}$$

Noting that $({\chi }_{oeeo}+{\chi }_{oeoe})/{\chi }_{eeee}$ is about $5:1$ [20], the calculated results of Eqs. (A14) and (A15) have been shown in Fig. 2 as red lines.

5. COMBINED SIMULATION VIA PC AND 4WM MODELS

As for the combination of 4WM and PC models (without 1DND model), except for the amplitude ratio of $0.42:1$ suggested by Ref. [26], other ratios have also been tried to reproduce the experimental results. The best fitting can be seen in Fig. 13 as the red lines when the amplitude ratio of 4WM/PC equals $2:1$. However, this disagrees with the fact that the contribution by 4WM is much weaker than that from the plasma (PC), especially in the low THz band [26]. Even if this was allowed, one can see the fitting of $A_{\mathrm{THz},o}$ in Fig. 13(a) has little improvement (red solid line versus green dashed line). Therefore, another physical model has been additionally considered, i.e., the 1DND model, as described in Figs. 5 and 6(a).

 

Fig. 13. Results of simulations via the combination of 4WM and PC models with different amplitude ratios of $0:1$, $1:1$, and $2:1$, for (a) $A_{\mathrm{THz},o}$ and (b) $A_{\mathrm{THz},e}$, along with the experimental results (black solid circles).

Download Full Size | PDF

6. ONE-DIMENSIONAL NEGATIVE DIELECTRIC WAVEGUIDE MODEL

The $\text{Re}[{ϵ}_r]$ distribution in Fig. 4(b) (black line) has been simplified as an ND homogeneous cylindrical structure with diameter of $2a$, ${ϵ}_{r1}=\text{Re}{[{ϵ}_r]}_{\min }$ for $r<a$ (core), and ${ϵ}_{r2}=1$ for $r>a$ (cladding). The electromagnetic field in a 1DND pin is assumed to be propagating in the $z$ direction in the form of $\exp [i(\omega t-\beta z)]$. In this situation, only the zeroth-order TM mode has been considered, since the TE mode is not a propagation mode in a 1D optical waveguide [42]. By analytically solving Maxwell’s equations in cylindrical coordinates, the TM-mode field in the ND core ($r<a$) with ${ϵ}_{r1}<0$ is given by

(A16)$$E_{z1}=A{I}_0({\gamma }_{1}r),E_{r1}=\frac{i\beta }{{\gamma }_1}A{I}_1({\gamma }_{1}r),H_{\varphi 1}=\frac{i\omega {ϵ}_{r1}{ϵ}_0}{{\gamma }_1}A{I}_1({\gamma }_{1}r),$$

and in the cladding ($r>a$) with ${ϵ}_{r2}=1$, the field is

(A17)$$\begin{gathered} E_{z2}=B{K}_0({\gamma }_{2}r), \\ {E}_{r2}=-\frac{i\beta }{{\gamma }_2}B{K}_1({\gamma }_{2}r), \\ {H}_{\varphi 2}=-\frac{i\omega {ϵ}_{r2}{ϵ}_0}{{\gamma }_2}B{K}_1({\gamma }_{2}r), \end{gathered}$$

where $A$ and $B$ are any constants, $r$ is the radius direction, ${ϵ}_0$ is permittivity in vacuum, $I_v$ and $K_v$ are the $v$th-order modified Bessel functions, and ${\gamma }_j$ ($j=1$, 2) is defined as ${\gamma }_j={[{\beta }^2-{ϵ}_{rj}{(\omega /c)}^2]}^{1/2}$. From the boundary conditions of $E_{z1}=E_{z2}$ and $H_{\varphi 1}=H_{\varphi 2}$ at $r=a$, the characteristic equation has been achieved as

(A18)$$\frac{{\gamma }_1{I}_0({\gamma }_{1}a)K_1({\gamma }_{2}a)}{{\gamma }_2{I}_1({\gamma }_{1}a)K_0({\gamma }_{2}a)}=-\frac{{ϵ}_{r1}}{{ϵ}_{r2}}.$$

Based on Eqs. (A16) and (A17), THz mode field intensity ${|E_r|}^2$ inside the plasma filament can be calculated.

It is also worth mentioning that, in Figs. 5 and 6, the FWHM of radial distribution of $N_e(r)$ has been unified as 50 μm for simplification. In this case, $\mathrm{\Delta }\eta ={\eta }_{b2}/{\eta }_{b1}=\sim 0.66$ at 0.7 THz, as shown as the red dashed line in Fig. 6(a). If precise values of 51 and 45 μm were used for FWHM at $N_e={\mathrm{b}}_1$ and ${\mathrm{b}}_2$, respectively, the corresponding $\mathrm{\Delta }{\eta }_{0.7\mathrm{THz}}$ would slightly change to $\sim 0.67$. This will not considerably alter the simulation results of 4WM + PC + 1DND models (pink lines in Fig. 2).

Compared with our 1DND model, numerical simulations (e.g., UPPE [35,36]) are more precise in quantity. However, a physical model has essential difference from numerical simulations. For instance, without fully running the complex numerical simulation, the 1DND model is able to tell immediately that inside the filament, the THz pulse beam diameter approximately coincides with the position where the plasma frequency is equal to the THz frequency [27]. This foresight with great guiding significance obviously cannot be comfortably forecasted by numerical simulations without large-scale computational resources and mass mathematical program operations.

7. MATHEMATICAL BASIS OF THE THREE-STEP MODEL

Here is the detailed calculation process. First, the generated THz amplitude as a function of BBO rotation angle $\theta$ has been calculated via the 4WM and PC models, respectively. The results are written as $A_{\mathrm{THz}}^{4\mathrm{WM}}(\theta )$ and $A_{\mathrm{THz}}^{\mathrm{PC}}(\theta )$. Then, based on Ref. [26], the combination of 4WM and PC models could be given by

(A19)$$A_{\mathrm{THz}}^{4\mathrm{WM}+\mathrm{PC}}(\theta )=0.42·\frac{A_{\mathrm{THz}}^{4\mathrm{WM}}(\theta )}{A_{\mathrm{THz}}^{4\mathrm{WM}}{(\theta )}_{\max }}+1·\frac{A_{\mathrm{THz}}^{\mathrm{PC}}(\theta )}{A_{\mathrm{THz}}^{\mathrm{PC}}{(\theta )}_{\max }}.$$

Next, the spatial confinement of the THz wave inside the plasma filament occurs. This is induced by the self-interaction between the generated THz wave [$A_{\mathrm{THz}}^{4\mathrm{WM}+\mathrm{PC}}(\theta )$] and the free electrons. The 1DND model has been proposed to reveal the underlying mechanism. Via 1DND model, THz-mode field intensity $I_{\mathrm{THz}}^{1\mathrm{DND}}$ has been found to be mainly constrained within the cross section of the filament (Fig. 5). A coefficient $\eta$ has been defined as the ratio of the confined THz amplitude (within $r<a$, the confinement region) to the total THz amplitude within the whole cross section of the plasma filament,

(A20)$${\eta }_{\mathrm{THz}}(\theta )=\sqrt{\frac{{\int }_0^{r<a}{I}_{\mathrm{THz}}^{1\mathrm{DND}}(r,\theta )·2\pi r·\mathrm{d}r}{{\int }_0^{\infty }{I}_{\mathrm{THz}}^{1\mathrm{DND}}(r,\theta )·2\pi r·\mathrm{d}r}}.$$

At last, the coupled result of 4WM + PC, i.e., $A_{\mathrm{THz}}^{4\mathrm{WM}+\mathrm{PC}}(\theta )$ in Eq. (A19) needs to be modified by ${\eta }_{\mathrm{THz}}(\theta )$ in Eq. (A20) with the following formula

(A21)$$A_{\mathrm{THz}}^{4\mathrm{WM}+\mathrm{PC}+1\mathrm{DND}}(\theta )=A_{\mathrm{THz}}^{4\mathrm{WM}+\mathrm{PC}}(\theta )·{\eta }_{\mathrm{THz}}(\theta ),$$

in order to remove the lossy THz amplitude within $r>a$ (the diffraction region).

8. JUSTIFICATION OF IGNORING THE PROPAGATION EFFECTS

Similar to Kim et al.’s [11] and Nguyen et al.’s [43] works, nonlinear propagation effects and self-actions along the $z$ axis, including self-phase modulation (SPM), have been ignored in our analysis. This simplification would be reasonable based on the following considerations.

In our experiment, the filament length is relatively short ($\sim 3\mathrm{mm}$), as shown in the inset of Fig. 4(a). SPM may not play a significant role because the broadening of the pump pulse spectrum induced by SPM is determined mainly by the propagation distance. Quantitatively, the nonlinear phase is given by

(A22)$$\mathrm{\Delta }\phi =k_0·n_{2}I·z.$$

By substituting typical values ($k_0=2\pi /\lambda$, $\lambda =800\mathrm{nm}$, $n_2=2.8\times {10}^{-19}{\mathrm{cm}}^2/\mathrm{W}$, $I{=5\times 10}^{13}\mathrm{W}/{\mathrm{cm}}^2$, $z=3\mathrm{mm}$) into Eq. (A22), $\mathrm{\Delta }\phi \sim 0.33$. This deduction has been confirmed by the previous work of Andreeva et al. [26]. Under similar experimental conditions as our work, based on the simulation of UPPE, Andreeva et al. have demonstrated that the THz signal induced by the Kerr nonlinearity could be quantitatively reproduced by the simple 4WM model. It implies that the spectral broadening of the pump pulse ($\omega$) and the variation of the relative phase between $\omega$ and $2\omega$, resulted from the SPM (or other nonlinear propagation effects) may not be significant. Hence, nonlinear propagation effects [44] and self-actions [45] along the $z$ axis, including SPM and self-steepening, have been ignored in this work.

Funding

National Key R&D Program of China (2017YFC0821300); National Natural Science Foundation of China (NSFC) (11704252, 11574160, 61722111); National Program on Key Basic Research Project of China (973 Program) (2014CB339802, 2014CB339806); Major National Development Project of Scientific Instrument and Equipment (2017YFF0106300, 2016YFF0100503); Tianjin Research Program of Application Foundation and Advanced Technology (15JCZDJC31700); Key Scientific and Technological Project of Science and Technology Commission of Shanghai Municipality (STCSM) (15DZ0500102); Shanghai Leading Talent (2016-019); Young Yangtze River Scholar (Q2016212); Tianjin Special Program for Talent Development.

REFERENCES

1. M. Tonouchi, “Cutting-edge terahertz technology,” Nat. Photonics 1, 97–105 (2007). [CrossRef]  

2. M. Clerici, M. Peccianti, B. E. Schmidt, L. Caspani, M. Shalaby, M. Giguère, A. Lotti, A. Couairon, F. Légaré, T. Ozaki, D. Faccio, and R. Morandotti, “Wavelength scaling of terahertz generation by gas ionization,” Phys. Rev. Lett. 110, 253901 (2013). [CrossRef]  

3. S. Wu, J. Liu, S. Wang, and Y. Zeng, “Influence of optical wavelength on terahertz radiation from laser-induced air plasma,” Chin. Opt. Lett. 11, 101402 (2013).

4. T. Bartel, P. Gaal, K. Reimann, M. Woerner, and T. Elsaesser, “Generation of single-cycle THz transients with high electric-field amplitudes,” Opt. Lett. 30, 2805–2807 (2005). [CrossRef]  

5. A. Houard, Y. Liu, B. Prade, and A. Mysyrowicz, “Polarization analysis of terahertz radiation generated by four-wave mixing in air,” Opt. Lett. 33, 1195–1197 (2008). [CrossRef]  

6. Y. Minami, T. Kurihara, K. Yamaguchi, and T. Suemoto, “High-power THz wave generation in plasma induced by polarization adjusted two-color laser pulses,” Appl. Phys. Lett. 102, 041105 (2013). [CrossRef]  

7. D. Dietze, J. Darmo, S. Roither, A. Pugzlys, J. N. Heyman, and K. Unterrainer, “Polarization of terahertz radiation from laser generated plasma filaments,” J. Opt. Soc. Am. B 26, 2016–2027 (2009). [CrossRef]  

8. J. Dai, N. Karpowicz, and X. C. Zhang, “Coherent polarization control of terahertz waves generated from two-color laser-induced gas plasma,” Phys. Rev. Lett. 103, 023001 (2009). [CrossRef]  

9. X. Xie, J. Dai, and X. C. Zhang, “Coherent control of THz wave generation in ambient air,” Phys. Rev. Lett. 96, 075005 (2006). [CrossRef]  

10. K. Y. Kim, A. J. Taylor, J. H. Glownia, and G. Rodriguez, “Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions,” Nat. Photonics 2, 605–609 (2008). [CrossRef]  

11. K. Y. Kim, J. H. Glownia, A. J. Taylor, and G. Rodriguez, “Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields,” Opt. Express 15, 4577–4584 (2007). [CrossRef]  

12. M. Kress, T. Löffler, S. Eden, M. Thomson, and H. G. Roskos, “Terahertz-pulse generation by photoionization of air with laser pulses composed of both fundamental and second-harmonic waves,” Opt. Lett. 29, 1120–1122 (2004). [CrossRef]  

13. Y. Liu, A. Houard, M. Durand, B. Prade, and A. Mysyrowicz, “Maker fringes in the terahertz radiation produced by a 2-color laser field in air,” Opt. Express 17, 11480–11485 (2009). [CrossRef]  

14. D. J. Cook and R. M. Hochstrasser, “Intense terahertz pulses by four-wave rectification in air,” Opt. Lett. 25, 1210–1212 (2000). [CrossRef]  

15. F. Blanchard, G. Sharma, X. Ropagnol, L. Razzari, R. Morandotti, and T. Ozaki, “Improved terahertz two-color plasma sources pumped by high intensity laser beam,” Opt. Express 17, 6044–6052 (2009). [CrossRef]  

16. Y. S. You, T. I. Oh, and K. Y. Kim, “Mechanism of elliptically polarized terahertz generation in two-color laser filamentation,” Opt. Lett. 38, 1034–1036 (2013). [CrossRef]  

17. C. D’Amico, A. Houard, M. Franco, B. Prade, A. Mysyrowicz, A. Couairon, and V. T. Tikhonchuk, “Conical forward THz emission from femtosecond-laser-beam filamentation in air,” Phys. Rev. Lett. 98, 235002 (2007). [CrossRef]  

18. Y. Zhang, Y. Chen, C. Marceau, W. Liu, Z.-D. Sun, S. Xu, F. Théberge, M. Châteauneuf, J. Dubois, and S. L. Chin, “Non-radially polarized THz pulse emitted from femtosecond laser filament in air,” Opt. Express 16, 15483–15488 (2008). [CrossRef]  

19. Y. Chen, C. Marceau, W. Liu, Z.-D. Sun, Y. Zhang, F. Théberge, M. Châteauneuf, J. Dubois, and S. L. Chin, “Elliptically polarized terahertz emission in the forward direction of a femtosecond laser filament in air,” Appl. Phys. Lett. 93, 231116 (2008). [CrossRef]  

20. Y. Zhang, Y. Chen, S. Xu, H. Lian, M. Wang, W. Liu, S. L. Chin, and G. Mu, “Portraying polarization state of terahertz pulse generated by a two-color laser field in air,” Opt. Lett. 34, 2841–2843 (2009). [CrossRef]  

21. Y. Chen, C. Marceau, S. Génier, F. Théberge, M. Châteauneuf, J. Dubois, and S. L. Chin, “Elliptically polarized terahertz emission through four-wave mixing in a two-color filament in air,” Opt. Commun. 282, 4283–4287 (2009). [CrossRef]  

22. H. Zhong, N. Karpowicz, and X. C. Zhang, “Terahertz emission profile from laser-induced air plasma,” Appl. Phys. Lett. 88, 261103 (2006). [CrossRef]  

23. Y. S. You, T. I. Oh, and K. Y. Kim, “Off-axis phase-matched terahertz emission from two-color laser-induced plasma filaments,” Phys. Rev. Lett. 109, 183902 (2012). [CrossRef]  

24. I. Babushkin, S. Skupin, A. Husakou, C. Koehler, E. Cabrera-Granado, L. Berge, and J. Herrmann, “Tailoring terahertz radiation by controlling tunnel photoionization events in gases,” New J. Phys. 13, 123029 (2011). [CrossRef]  

25. H. G. Roskos, M. D. Thomson, M. Kreß, and A. T. Löffler, “Broadband THz emission from gas plasmas induced by femtosecond optical pulses: from fundamentals to applications,” Laser Photon. Rev. 1, 349–368 (2007). [CrossRef]  

26. V. A. Andreeva, O. G. Kosareva, N. A. Panov, D. E. Shipilo, P. M. Solyankin, M. N. Esaulkov, P. González de Alaiza Martínez, A. P. Shkurinov, V. A. Makarov, L. Bergé, and S. L. Chin, “Ultrabroad terahertz spectrum generation from an air-based filament plasma,” Phys. Rev. Lett. 116, 063902 (2016). [CrossRef]  

27. J. Zhao, W. Chu, Z. Wang, Y. Peng, C. Gong, L. Lin, Y. Zhu, W. Liu, Y. Cheng, S. Zhuang, and Z. Xu, “Strong spatial confinement of terahertz wave inside femtosecond laser filament,” ACS Photon. 3, 2338–2343 (2016). [CrossRef]  

28. K. Y. Kim, “Generation of coherent terahertz radiation in ultrafast laser-gas interactions,” Phys. Plasmas 16, 056706 (2009). [CrossRef]  

29. S. Xu, Y. Zhang, W. Liu, and S. L. Chin, “Experimental confirmation of high-stability of fluorescence in a femtosecond laser filament in air,” Opt. Commun. 282, 4800–4804 (2009). [CrossRef]  

30. J. Zhao, W. Chu, L. Guo, Z. Wang, J. Yang, W. Liu, Y. Cheng, and Z. Xu, “Terahertz imaging with sub-wavelength resolution by femtosecond laser filament in air,” Sci. Rep. 4, 3880 (2014). [CrossRef]  

31. A. Gorodetsky, A. D. Koulouklidis, M. Massaouti, and S. Tzortzakis, “Physics of the conical broadband terahertz emission from two-color laser-induced plasma filaments,” Phys. Rev. A 89, 033838 (2014). [CrossRef]  

32. P. Klarskov, A. C. Strikwerda, K. Iwaszczuk, and P. U. Jepsen, “Experimental three-dimensional beam profiling and modeling of a terahertz beam generated from a two-color air plasma,” New J. Phys. 15, 075012 (2013). [CrossRef]  

33. A. Houard, Y. Liu, B. Prade, V. T. Tikhonchuk, and A. Mysyrowicz, “Strong enhancement of terahertz radiation from laser filaments in air by a static electric field,” Phys. Rev. Lett. 100, 255006 (2008). [CrossRef]  

34. J. Zhao, Y. Zhang, Z. Wang, W. Chu, B. Zeng, W. Liu, Y. Cheng, and Z. Xu, “Propagation of terahertz wave inside femtosecond laser filament in air,” Laser Phys. Lett. 11, 095302 (2014). [CrossRef]  

35. L. Bergé, S. Skupin, C. Köhler, I. Babushkin, and J. Herrmann, “3D numerical simulations of THz generation by two-color laser filaments,” Phys. Rev. Lett. 110, 073901 (2013). [CrossRef]  

36. P. G. de Alaiza Martínez, X. Davoine, A. Debayle, L. Gremillet, and L. Bergé, “Terahertz radiation driven by two-color laser pulses at near-relativistic intensities: competition between photoionization and wakefield effects,” Sci. Rep. 6, 26743 (2016). [CrossRef]  

37. Handbook of Optics, 2nd ed. (McGraw-Hill, 1994), Vol. 2.

38. W. J. Smith, Modern Optical Engineering, 3rd ed. (McGraw-Hill, 2000).

39. T. I. Oh, Y. S. You, and K. Y. Kim, “Two-dimensional plasma current and optimized terahertz generation in two-color photoionization,” Opt. Express 20, 19778–19786 (2012). [CrossRef]  

40. R. W. Boyd, Nonlinear Optics (Academic, 2003).

41. N. Aközbek, A. Iwasaki, A. Becker, M. Scalora, S. L. Chin, and C. M. Bowden, “Third-harmonic generation and self-channeling in air using high-power femtosecond laser pulses,” Phys. Rev. Lett. 89, 143901 (2002). [CrossRef]  

42. J. Takahara, “Negative dielectric optical waveguides for nano-optical guiding,” in Plasmonic Nanoguides and Circuits, S. I. Bozhevolnyi, ed. (Pan Stanford Publishing, 2008), pp. 46–47.

43. A. Nguyen, P. G. de Alaiza Martínez, I. Thiele, S. Skupin, and L. Berge, “THz field engineering in two-color femtosecond filaments using chirped and delayed laser pulses,” New J. Phys. (accepted), doi: 10.1088/1367-2630/aaa470. [CrossRef]  

44. A. Nguyen, P. G. de Alaiza Martinez, J. Déchard, I. Thiele, I. Babushkin, S. Skupin, and L. Bergé, “Spectral dynamics of THz pulses generated by two-color laser filaments in air: the role of Kerr nonlinearities and pump wavelength,” Opt. Express 25, 4720–4740 (2017). [CrossRef]  

45. E. Cabrera-Granado, Y. Chen, I. Babushkin, L. Bergé, and S. Skupin, “Spectral self-action of THz emission from ionizing two-color laser pulses in gases,” New J. Phys. 17, 023060 (2015). [CrossRef]