Spectral dynamics of THz pulses generated by two-color laser filaments in air: The role of Kerr nonlinearities and pump wavelength

We theoretically and numerically study the influence of both instantaneous and Raman-delayed Kerr nonlinearities as well as a long-wavelength pump in the terahertz (THz) emissions produced by two-color femtosecond filaments in air. Although the Raman-delayed nonlinearity induced by air molecules weakens THz generation, four-wave mixing is found to impact the THz spectra accumulated upon propagation via self-, cross-phase modulations and self-steepening. Besides, using the local current theory, we show that the scaling of laser-to-THz conversion efficiency with the fundamental laser wavelength strongly depends on the relative phase between the two colors, the pulse duration and shape, rendering a universal scaling law impossible. Scaling laws in powers of the pump wavelength may only provide a rough estimate of the increase in the THz yield. We confront these results with comprehensive numerical simulations of strongly focused pulses and of filaments propagating over meter-range distances.


Introduction
Laser filaments produced by ultrashort light pulses proceed from the dynamic balance between Kerr self-focusing and plasma generation [1,2]. The interplay of these nonlinear effects contributes to broaden the pulse spectrum, promote self-compression [3] and self-guide intense optical wave packets over remote distances for, e.g., sensing applications [4]. Filamentation of pulses with different frequencies has recently been proposed as an innovative way to downconvert optical radiation into the THz range [5] and create broadband THz sources remotely [6]. Terahertz pulse generation by two-color (fundamental and its second harmonic) laser filaments mainly results from photocurrents [7,8], through which photoionization in tunneling regime drives the conversion process [9]. However, because THz radiation is also emitted by optical rectification through four-wave mixing [10], Kerr self-focusing can contribute to the overall THz yield [11]. Kerr-driven THz sources are expected to emit on-axis, mostly during the early self-focusing stage preceding ionization of air molecules [12], and their characteristic spectrum exhibits a parabolic distribution being maximum at non-zero frequency [10]. At clamping intensities > 50 TW/cm 2 from which Kerr self-focusing is stopped by plasma defocusing, THz emission is dominated by photocurrents and the pulse spectrum peaks at smaller THz frequencies [13].
Several important issues still need to be addressed in this physics, such as the role of the Raman-delayed part of the Kerr nonlinearity, which arises due to the excitation of rotational and vibrational transitions of the molecular constituents of air [14]. Another point is the impact of longer pump wavelengths (λ ), in particular in the near-infrared domain ∼ 1.6 µm that may be preferred for, e.g., ocular safety reasons. Antecedent studies [13,15] reported impressive growths of THz energy yields scaling as λ α with α ≃ 4 − 5. More recent ones [16,17], however, showed that, although single-color laser pulses produce THz energy yields increasing like λ 4 , no similar conclusion could be inferred for two-color pulses.
The present paper aims at clarifying the role of the pump wavelength and of the Kerr nonlinearities (instantaneous and delayed) in filament-driven THz pulse generation. We display numerical evidence that four-wave mixing impacts the THz generation process over long propagation distances, even in intensity regimes where the photocurrent mechanism is the dominating THz emitter. We also demonstrate that the Raman-delayed Kerr nonlinearity does not contribute as a THz source. By means of the local current (LC) model [18], we moreover explain the variations in the THz field strength with respect to the fundamental wavelength of the optical radiation. In this respect we emphasize the role of the electron current component associated with the high-frequency laser pulse, whose fundamental at longer wavelength affects more significantly the THz spectrum. Importantly, the scaling of the THz field strength with the fundamental pump wavelength is shown to vary with the relative phase between the fundamental and second harmonic, the pulse envelope and duration, so that no universal λ -dependent scaling is achievable with a two-color pulse. Despite this, the conversion efficiencies reported from the LC model show that the THz yield is roughly scaling as λ α with α > 4 for small relative phases between the two colors and α > 2 on the average. For focused pulses [15], these results are confirmed by direct simulations employing a unidirectional pulse propagator [19,20]. Smaller gain factors are achieved by meter-range two-color filaments due to the generation of weaker plasma densities.
The paper is organized as follows: Section 2 proposes a one-dimensional (1D) approach combining known laser-driven THz sources. It recalls that, in the range of intensities reached by two-color filaments in air, photoionization and to a lesser extent the Kerr nonlinearity are the principal players in THz generation. Section 3 discusses analytical estimates of the laser-to-THz conversion efficiency when taking a Raman-delayed nonlinearity into account and when increasing the fundamental laser wavelength. Section 4 verifies our analytical statements through three-dimensional (3D) comprehensive numerical simulations for both focused and filamentary pulses.

Transverse versus longitudinal THz fields -A 1D Approach
Before proceeding with full 3D simulations, we find it instructive to examine the dynamics of THz fields produced in-situ, i.e., inside the plasma channel created by two-color filaments in air. For simplicity, we use a reduced model discriminating THz transverse (x) from longitudinal (z) currents for a laser pulse polarized along the x-axis. Discarding ion motions, the current density induced by the photoionized electrons is J ≃ −eN e v e , where N e and v e are the free electron density and velocity. Following Sprangle et al. [22], the electron current obeys the following equation set in the non-relativistic interaction regime: where ν c is the electron collision rate equal to 10 ps −1 , B is the magnetic field associated with the electric field E that includes both the laser field E L e x and secondary (THz) radiation such as Ẽ =Ẽ x e x +Ẽ z e z . Assuming singly-ionized gases at moderate intensities I 0 ≤ 10 15 W/cm 2 , the growth of electron density is governed over femtosecond time scales by where N a is the initial gas density and W (E) is a field-dependent ionization rate, for instance the quasi-static tunnel (QST) rate [21] W with U i being the ionization energy of the gas, ν a = 4.13 × 10 16 Hz, E a = 5.14 × 10 11 V/m and U H = 13.6 eV is the ionization potential of hydrogen. For air composed of 80% of N 2 and 20% of O 2 , we shall start by considering ionization of oxygen molecules only, as their ionization potential (12.1 eV) is lower than that of nitrogen (15.6 eV). Photoionization of nitrogen molecules will be addressed later when considering, e.g., a field-dependent version of the cycle-averaged ionization rate derived by Perelomov, Popov and Terent'ev (PPT) [24].
Equations (1) and (2), combined with Maxwell-Ampère equation ∇ × B = c −2 ∂ t E + µ 0 J, readily provide the propagation equation for E: where ω pe = e 2 N e /ε 0 m e is the electron plasma frequency (ε 0 = 1/µ 0 c 2 ). We here neglect loss currents due to photoionization, which are small for our pump pulse configurations and in gas-based plasmas.
For technical convenience, Eq. (5) is reduced to a one-dimensional, z-propagating model. Discarding temporarily the linear (chromatic) and nonlinear polarizations of air molecules, we omit the diffraction operators (∂ x = ∂ y = 0), yielding As recently derived in [26], further approximations can be applied such as assuming a laser field propagating with the sole variable (z − ct). For the intensity range < 10 15 W/cm −2 , we can furthermore neglect the longitudinal current component compared to its transverse counterpart that obeys like the laser pulse. By splitting J x = J L +J x on the expansion E x = E L +Ẽ x , Eq. (6) reads as The equation for the longitudinal radiated fieldẼ z is easily expressed as Equations (8) and (9) both involve the current density computed on the laser field, J L = ε 0 (∂ t + ν c ) −1 ω 2 pe E L . The transverse radiated fieldẼ x mainly proceeds from this current density [16,18]. With two colors, the product in ∂ t J L [Eq. (8)] between the steplike increase of N e (t) and the fast oscillations of E L (t) acts as an efficient converter to low frequencies. By comparison, Eq. (9) describes longitudinal plasma oscillations that develop over longer time scales after the laser field has interacted with the gas, as previously established in [22,23].
Including the optical polarization of a noble gas is straightforward by replacing into the Maxwell-Ampère equation the electric field E by the displacement vector D = ε 0 E + P L + P NL , where P L and P NL refer to linear and nonlinear polarization, respectively. In scalar description, this amounts to adding into the right-hand side of Eq. (6) the term −ε −1 0 ∂ 2 t (P L + P NL ). The Fourier transform of P L involves the first-order frequency-dependent susceptibility χ (1) (ω) entering the optical linear index n(ω) = [1 + χ (1) (ω)] 1/2 . P NL involves the third-order susceptibility χ (3) responsible for four-wave mixing and Kerr self-focusing (in multidimensional media).
For molecular gases, the Kerr response admits a fraction x k (0 ≤ x k ≤ 1) of delayed contribution due to Raman scattering by rotational molecular transitions. Assuming that the laser is not resonant with the transition frequencies [14], stimulated Raman scattering usually affects the total time-dependent refraction index of the medium. Its corresponding polarization component describes a delayed-Kerr nonlinearity. Written with the real electric field E x [2], the overall nonlinear polarization can be expressed as where Here, symbol * means convolution product in time, Θ(t) is the usual Heaviside step function, τ 1 and τ 2 represent the rotational Raman time and dipole dephasing time, respectively. In air those quantities take the values τ 1 ≈ 62 fs and τ 2 ≈ 77 fs for the fraction x k = 0.5 [14,25]. We shall assume that these values still hold for pump wavelengths up to 2 µm.
To start with, we integrate the (1+1)-dimensional Eqs. (6) and (7) using the initial condition where I 0 is the mean pump intensity, r is the relative intensity ratio of the second harmonic, ϕ is the relative phase between the fundamental pulse with carrier frequency ω 0 and its second harmonic. For I 0 = 100 TW/cm 2 and r = 0.15, Figs. 1(a,b) show the transverse and longitudinal secondary fields filtered in a 80-THz frequency window for a two-color pulse propagating in air at ambient pressure. The full width at half maximum (FWHM) duration of the pump pulse centered at 800-nm wavelength is τ p = 40 fs with half-duration for the 2ω 0 component. Simulations are performed using the quasi-static tunneling ionization, accounting or not for the Kerr nonlinearity, which we here consider instantaneous (x k = 0) with a nonlinear index equal to 10 −19 cm 2 /W. For simplicity we first ignore linear dispersion whose action will be examined in Figs. 2(e,f). At such intensities O 2 molecules undergo most of the ionization events for a neutral density N a = 5.4 × 10 18 cm −3 . With two-color pulses, the transverse field at the distance z = 1 cm is found to be orders of magnitude larger than the longitudinal one [ Fig. 1(a)]. Plasma wakefield effects characterize the longitudinal field [ Fig. 1(b)], with a long plasma wave formed behind the pulse head and oscillating at plasma frequency, as shown in the spectrum (see inset). Note that the plasma frequency varies due to small changes in N e caused by Kerr-induced self-steepening. At larger propagation distances, the transverse THz field increases, whereas the longitudinal field decreases even more (not shown). The Kerr response, although of minor role in the conversion process, increases the peak of the transverse THz field to some extent. Hence, at air-based filament intensities ∼ 100 TW/cm 2 , only the transverse secondary fields generated through photocurrents and four-wave mixing appear to be relevant in a two-color pulse configuration (see also [26]).
Figures 2(a-d) detail the influence of the Kerr (blue dash-dotted curves), plasma (red dashed curves) and combined Kerr-plasma effects (black solid curves) when solving Eq. (6) for a fundamental pulse at 800 nm coupled with its second harmonic interacting with ambient air over short propagation distance (z = 1 cm). The left-hand side (LHS) column depicts the transverse THz field distributions computed from an inverse Fourier transform of the  6) and (7), and filtered in a 80-THz window for a 800+400-nm pulse with I 0 = 100 TW/cm 2 , ϕ = π/2 and r = 15%. The red dashed curve displays the THz field computed with the plasma response only (P); the solid curve corresponds to both Kerr (instantaneous) and plasma responses (K+P). Inset in (b) details the THz spectrum of the longitudinal field. whole field within the frequency window ν < 80 THz. The right-hand side (RHS) column shows corresponding spectra. At low intensities (25 TW/cm 2 ), which are characteristic of a self-focusing regime prior to ionization, the THz field exhibits a temporal profile shaped by the four-wave mixing contribution, here optimized with an initial zero relative phase between the two colors [10]. However, the field maximum reaches much higher values in the presence of photoionization [ Fig. 2(a)]. The same conclusion holds in the reciprocal configuration, when the laser pulse evolves in self-channeling regime favoring THz emissions by photocurrents with ϕ = π/2 [7,28]. The THz field profile follows that created by a plasma source alone, but it appears noticeably amplified by the Kerr term [ Fig. 2(c)], which also agrees with Ref. [12]. Along propagation this optical nonlinearity contributes through self-, cross-phase modulations and self-steepening to broaden the pump and low-order (second as well as third) harmonics.
Keeping a relatively small amplitude as pure THz emitter [see Fig. 2(d)], the Kerr response mainly acts by changing the pulse spectrum. More phase-matched components increase the contribution of four-wave mixing to the THz spectrum and can constructively interfere with those driven by plasma generation. In Figs. 2(b,d), we refind the typical Kerr spectral signature in ν 2 [11,13] that amplifies larger frequencies in the THz spectrum. We can also invoke the self-action of a THz seed originating from one potential source (e.g., Kerr or plasma) and augmenting the overall THz spectrum [27]. Note the important increase in the radiated field emerging at low laser intensity in Fig. 2(b). We observed that variations in the phase angle remain small and cannot solely switch on plasma emitters. So we suspect here that non-trivial couplings between the backward and forward THz components may enhance the THz yield in the low-frequency limit, which is beyond the scope of the present paper.
So far, linear dispersion has been discarded over propagation ranges ≤ 1 cm along which its action is usually expected to be small in gases [2]. However, over comparable ranges and depending on the pump wavelength, it may already significantly impact the relative phase ϕ that conditions the efficiency of the THz emitters. To illustrate this point, Figs. 2(e,f) display the evolution of the phase angle ϕ numerically extracted from the forward component of the electric field along the z axis. The forward electric field is computed from the 1D version of the UPPE model [see Eq. (33)] for our two pulse configurations that now undergo air dispersion as modeled in Ref. [29]. At low intensity [25 TW/cm 2 - Fig. 2(e)], for which the plasma response is small, linear dispersion drags the relative phase out of its initial value upon short distances ∼ 1.5 cm by a phase shift comparable to that driven by the Kerr nonlinearity. A similar conclusion applies to the plasma regime [100 TW/cm 2 - Fig. 2(f)]. Despite the smallness of its coefficients in air [29], linear dispersion induces a long phase mismatch between the 800-nm and 400-nm pulse components, which, combined with the nonlinearities, is able to drive a phase shift close to π/2 over 1.5 cm that cannot be neglected. This constraint is relaxed to some extent for longer pump wavelengths.
In summary, the above results show a net influence in the THz yield when accounting for the Kerr nonlinearity along cm-propagation ranges. The Kerr response directly alters the pump spectrum as well as variations in the phase angle between the ω − 2ω pulse components. Together with linear dispersion, this impacts the THz conversion efficiency, which is mainly driven by the photocurrent mechanism at clamping intensity.

Impact of delayed Kerr nonlinearities and longer pump wavelength
Below we address the influence of the Raman-delayed Kerr nonlinearity on the laser-to-THz conversion efficiency and we quantify the increase in THz generation when the fundamental wavelength belongs to the mid-infrared range. Since from Fig. 1 we expect no significant action from the longitudinal field E z , we only focus on the transverse field E x , whose subscript x is henceforth omitted.
Using the input two-color pulse (13), we can evaluate the low-frequency part of the overall Kerr response Eq. (10), when both envelope functions E ω 0 ,2ω 0 (t) take the value unity for the sake of simplicity. Cook and Hochstrasser's result [10] is easily recovered for the instantaneous part of the Kerr polarization yielding the direct-current (dc) contribution THz emission due to an instantaneous Kerr response is maximum for the phase offset ϕ = 0 [π].
Adding the Raman contribution now leads us to evaluate the integral After several trigonometric simplifications, we can extract from the previous expression a low-frequency (dc) contribution reading as where The values of T 1 and T 2 provide an optimal phase for THz generation when ϕ = tan −1 (T 2 /T 1 ). For pump wavelengths λ 0 comprised between 0.8 and 3 µm, we find that |T 1 | ≤ 0.021 and |T 2 | ∼ 0.02|T 1 |. It is thus reasonable to neglect T 2 , so that the low-frequency part of the total Kerr response involving the Raman-delayed component simplifies into By analogy with Cook and Hochstrasser's analysis, this formula indicates that the optimal phase difference for THz generation by the Kerr response remains ϕ = 0 [π]. Maximum THz generation is provided when there is no Raman nonlinearity (x k = 0). P Raman decreases the THz yield by a correction (T 1 ) of the percent order, which is confirmed by Figs. 3(a,b) for different values of the fraction x k . THz spectra and fields are computed at z = 0 from Eq. (8) including the nonlinear polarization. They indeed decrease in amplitude when x k is augmented. This property reflects the fact that the nonlinear integrand in Eq. (15) acts over relatively long relaxation times τ 1 − τ 2 ∼ 60 − 80 fs, along which the high-frequency laser oscillations cancel each other over the integration in time. The resulting integral is slowly varying and P Raman thus barely contributes to the THz spectrum.
This result confirms the behavior expected from envelope-like unidirectional models [2] for which the Raman nonlinearity is assumed insensitive to the pump harmonics and oscillates likely to the optical field. It justifies that we can employ Eq. (10) within a field description and still applies to alternative formulations of the rotational Raman scattering [30]. Besides the net decrease by the fraction x k , the impact of the overall Kerr source is also expected to decrease from the loss of self-and cross-phase modulations that affect the pump spectrum through the Kerr response. This aspect is detailed in Figs. 3(c,d), where we have compared the THz spectra obtained from inserting into the source terms ∂ t J and/or ∂ 2 t P NL some temporal profiles further computed in Section 4 in focused geometry [ Fig. 3(c)] and filamentation regime [ Fig. 3(d)] at the distance of maximum THz generation. Although this procedure loses the memory of the THz yield accumulated along previous distances, including that prior to ionization, it clearly indicates that in plasma regime the Kerr source remains minor compared to that associated to photocurrents. The action of the Raman nonlinearity mostly manifests by changing the pulse spectrum, which conditions the photocurrents.

Increase in the pump wavelength
References [13,15] reported an impressive increase in the THz conversion efficiency when doubling the fundamental wavelength for equal FWHM durations. For single-color pulses of same energy, we can easily expect that doubling λ 0 for THz generation triggered in tunneling regime keeps the final electron density N e unchanged, but it doubles the free electron velocity through its evident dependency on 1/ω 0 . Thereby the electron current density is doubled. With two colors involving a dominant fundamental pulse at ω 0 , the same scaling holds. However, N e noticeably increases, e.g., by a factor ∼ 2 at 100 TW/cm 2 , for two superimposed colors (ϕ = 0). With a π/2 relative phase, the pump field exhibits a temporal asymmetry around the electric field maxima and the current density J(t) develops a low-frequency component due to the stepwise increase of the electron density. This component is then the major THz source [28].
Ignoring again the influence of the envelopes (E ω 0 ,2ω 0 = 1), the ionization instants of a twocolor pulse can be approximated by [18] at leading order in a 2ω 0 /a ω 0 ≪ 1. We assume equal density jumps occurring over a large number of ionization events N ≫ 1. In the THz frequency range ω ≪ ω 0 , it is then straightforward to evaluate Equations (28) and (29) show that ∂ t J A and ∂ t J B dominate for ϕ = 0 and ϕ = π/2, respectively. The ionization steps δ N [Eq. (27)] increase, in the limit W τ n ≪ 1, linearly with the ionization duration τ n ≈ τ 1 . For a fixed pulse duration, the number of optical cycles is halved when one doubles the pump wavelength and since τ n ∝ ω −1 0 , one has τ 2λ 0 /τ λ 0 ∼ 2. The number of ionization events decreases accordingly, i.e., N 2λ 0 /N λ 0 = 1/2.

For a better understanding of ∂ t J A , it may be instructive to rewrite Eq. (24) in the form
where r f denotes the free electron position (r f = 0 at −∞): In Figs. 4(a,b), stricto-sensu, THz emission should not be zero due to the pulse envelope [18]. Considering the influence of bounded envelopes, Fig. 5 shows the ratio between the THz field induced by a two-color pump-pulse with fundamental wavelength at 1600 nm and one with fundamental wavelength at 800 nm. This ratio is evaluated through an ordinary least squares method applied to the THz field profiles calculated numerically from the LC model for ionizing beam intensities. This method mostly captures the ratio between the THz field maxima along the time axis. For comparison, the theoretical ratio |Ẽ 2λ 0 /Ẽ λ 0 | inferred from the inverse Fourier transform of Eq. (30) is plotted as a black solid line. It provides a gain factor that varies between 4 (ϕ = 0) and 2 (ϕ = π/2) with π-periodicity. For Gaussian pulses [E ω 0 ,2ω 0 (t) = exp (−2 2β −1 ln 2t 2β /τ 2β ω 0 ,2ω 0 ) with β = 1], this behavior is verified by the LC results (red curves) of Fig. 5(a), despite minor variations caused by envelope effects. The ratio |Ẽ 2λ 0 /Ẽ λ 0 | remains less sensitive to the pulse duration and the fundamental wavelength than to the relative phase ϕ. Maximum gain factor in the THz field amplitudes is obtained for ϕ ≈ 0, which underlines the role of the current density J A directly connected to the laser field. Opting next for 4th-order super-Gaussian profiles (β = 4), |Ẽ 2λ 0 /Ẽ λ 0 | again varies with the relative phase, but it strongly evolves and even exceeds the value 5 for ϕ ≤ π/10 when decreasing the FWHM duration [ Fig. 5(b)]. We attribute these changes to the steepness of the envelope, which makes the first ionization events not exactly located at the same times for 800-nm and 1600-nm pump pulses. For our laser parameters, the relative intensity ratio r appears to have a limited impact on the gain performances. Figures 5(c,d) illustrate the increase in the THz energy yield for 60-fs FWHM pump duration and intensities ∼ 200 TW/cm 2 rather reached in focusing geometries [15]. Solid lines refer to the computed energy values, while the dashed lines are fitting curves in λ α . We can observe that exponents α > 4 fit for small relative phases. A π/2 phase angle, however, renders the J B contribution dominant and thus decreases this exponent. So, even though λ -dependent scalings reported in [15] are possible, they are not generic as the gain factors are highly sensitive to the relative phase between the two colors, the shape of the pulse envelopes and their durations. Note from Figs. 5(c,d) that the THz pulse energy is much larger with a π/2 phase angle than with a null phase. This means that in the situation where J A is dominant (φ → 0), the overall THz spectrum is much smaller than when J B prevails for different phase values.

Comparison with unidirectional pulse propagation simulations
The previous properties are now checked by direct 3D numerical computations. Our reference model is the unidirectional pulse propagation equation (UPPE) [19,20] that governs the whereÊ(k x , k y , z, ω) is the Fourier transform of the laser electric field with respect to x, y, and t. The first term on the right-hand side of Eq. (33) describes linear dispersion and diffraction of the pulse. The termF NL =P NL + iĴ/ω + iĴ loss /ω contains the third-order nonlinear polarization with Kerr index n 2 = 3χ (3) /4n 2 0 cε 0 [n 0 = n(ω 0 )], the electron current J and a loss term J loss due to ionization [2,3]. Compared with Ref. [20], the denominator of the nonlinear term reduces to 2k(ω), assuming ω 2F NL relevant only for k(ω) = n(ω)ω/c ≫ k 2 x + k 2 y .
The coming analysis aims at first validating our theoretical expectations using the simple QST model for a single species (O 2 ) and classical values for the Kerr indices, in both focused and parallel propagation geometries. Next, more elaborated ionization models and recently measured Kerr coefficients in the mid-infrared will be employed to countercheck our findings.

Validation of theoretical issues
Equation (33) is here solved for experimental configurations close to those examined in [15], i.e., for focused pulses with f -numbers > 10 ( f -number refers to the ratio between the focal length and the FWHM input beam diameter). The initial relative phase ϕ is set equal to zero and the fundamental pump wavelengths of the two-color pulses are 800 nm and 1600 nm.
In a first set of simulations we choose f /# = 42 for the beam width w 0 = 500 µm and focal length f = 2.5 cm. Our two-color pulses have 200 µJ in energy and the FWHM duration of the pump pulse is 60 fs. About 7 % of the laser energy is contained in the second harmonic. Dispersion curve in air for the refractive optical index n(ω) is again taken from Ref. [29]. Kerr indices are chosen as n 2 ≃ 1.2 × 10 −19 cm 2 /W following Refs. [32,33]. At atmospheric pressure, N a = 5.4 × 10 18 cm −3 for O 2 molecules and the critical power for self-focusing, defined by P cr ≃ λ 2 0 /2πn 0 n 2 , is P cr = 8.5 GW at 800 nm and 35.1 GW at 1600 nm. Although strongly focused, our ultrashort pulses promote single-ionization events for peak intensities < 300 TW/cm 2 near focus due to local defocusing by the generated plasma. Therefore, we can here use the QST ionization rate (4) applied to oxygen molecules only.
Simulations have been performed with a time window of 1.22 ps, a temporal step ∆t = 75 attoseconds and transverse resolution of ∆x = ∆y ≈ 3 µm. Figure 6 shows the peak electron density reached near focus, the variations of the relative phase between the fundamental and second harmonic along z [ Fig. 6(a)], and the THz energy contained in our numerical box (3 × 3 mm 2 ) [ Fig. 6(b)]. THz radiation is collected within a 80-THz-large frequency window. Cyan/magenta curves ignore the delayed Raman nonlinearity; blue/red curves include it for comparison (x k = 0.5). There is a limited influence of the Kerr response in tight focusing regime, but its related THz conversion efficiency is clearly diminished by the delayed nonlinearity for the reasons given in Section 3. Concerning the wavelength dependency, the maximum intensity achieved near focus decreases with λ 0 , as the beam waist w f ≈ w 0 f /z 0 becomes proportional to the pump wavelength when the Rayleigh length z 0 = πw 2 0 /λ 0 is much larger than the focal distance f [35]. Consequently, since the QST rate (4) does not depend on the pump wavelength, the peak electron density decreases in turn [ Fig. 6(a)]. The relative phase ϕ covers the full interval [0, 2π] over the 4-cm-long propagation range. It experiences a Gouy phase shift up to π near focus, supplemented by another π phase shift induced by linear dispersion for the 800-nm pump [see Fig. 2(f)]. In Fig. 6(b) we observe a net increase of the maximum THz energy produced at z ≃ f when λ 0 is augmented.
The experimental data points of Ref. [15] are recalled by solid red circles in Fig. 6(c), which we compare with the THz pulse energy at focus. Despite differences between the original experiment and our laser parameters, the THz yield only evaluated from two pump wavelengths follows a comparable growth. Accounting for Raman scattering helps to reach a better agreement with the experimental data. For our THz window of 80 THz, a rough fitting curve indicates a growth rate in λ α with α ≈ 3.5, i.e., 2 < α ≤ 4 in agreement with Fig. 5(c,d), keeping in mind the variations in the relative phase ϕ shown in Fig. 6(b). Shortening this window to 20 THz as measured in [15] does not noticeably change this scaling, as the THz spectra emitted around focus are self-contained in the frequency domain ν < 20 THz [ Fig.   Fig. 6. Focused two-color Gaussian pulses in ratio r = 7.4%. (a) Maximum electron density (left-hand side axis, solid curves) and relative phase between the two pulse components (right-hand side axis, dashed curves), (b) THz energy yield (ν ≤ 80 THz) for the pump wavelengths 800 nm (blue/cyan curves) and 1600 nm (red/magenta curves), with and without the Raman nonlinearity. (c) THz energy vs pump wavelength. Cyan diamonds: no Raman; blue squares: With Raman. Green crosses × report THz gain factors in filamentation regime with no Raman nonlinearity as promoted in Fig. 7; Green symbols + report gain factors in filamentation regime with Raman nonlinearity as given by Fig. 9. The solid red circles recall the experimental data points of Ref. [15]. (d) Normalized THz spectra computed near focus. 6(d)]. The spectral distributions remarkably support the comparison with the measurements of [15], in particular by achieving a maximum at ν = 5 THz. The influence of the plasma volume is here limited: Visual inspection of the numerical data indeed revealed plasma channels being of comparable dimensions when using a pump wavelength of either 800 nm or 1600 nm, i.e., the plasma volumes only vary within a factor 0.8 -1.3 from density levels > 10 15−17 cm −3 .
We now employ Eq. (33) to describe THz generation by two-color filaments operating with two different pump wavelengths over long distances. The UPPE is integrated for two-color Gaussian pulses with input power P in = 34 GW, beam waist w 0 = 400 µm and FWHM durations τ ω 0 = 40 fs (τ 2ω 0 = τ ω 0 /2, r = 3.4%) in collimated propagation. Kerr indices are unchanged and we keep the quasi-static tunneling rate (4). For reasons of computational cost, the numerical resolution has been decreased to ∆t = 99 attoseconds and ∆x = ∆y ≈ 9 µm, which was checked to introduce no significant variations in the THz spectra and fields. The selected THz window is still ν ≤ 80 THz. Colored curves refer to wave propagation without (cyan/magenta curves) and with (blue/red curves) Raman-delayed Kerr nonlinearity. Bottom: Corresponding THz pulse energy yield accumulated along z inside a 80-THz-large frequency window. (c,d) On-axis THz fields and intensity spectra near the distance of maximum THz energy yield. Figures 7(a,b) illustrate the peak intensity and maximum electron density reached along meter-range distances in two-color filamentation regime. Light (cyan and magenta) curves show a propagation for which no delayed Kerr nonlinearity is accounted for (x k = 0). Dark (blue and red) curves include the Raman-delayed nonlinearity in ratio x k = 0. 5 [25, 37, 38]. From Fig. 7(a) it is clear that the Raman term weakens the contribution of the instantaneous Kerr response, which results in (i) a longer self-focusing distance, (ii) a decreased clamping intensity, and thereby (iii) an enhanced self-guiding range. Consistently, the peak plasma density decreases in turn and extends over longer distances. As expected, the bottom row of Figs. 7(a,b) displays a net decrease of the THz energy emitted along propagation. With a pump wavelength of 1600 nm, at equal energy content, the input power becomes closer to a single critical power in air and self-channeling with no delayed-Kerr response favors a more extended filamentation range. Peak densities decrease from 10 17 cm −3 to 10 16 cm −3 , which should weaken the high gain factors achieved in focused geometry. Indeed, the THz yield is only stronger (without Raman), locally by a factor ∼ 1.52 at maximum emission, as indicated by the green crosses of Fig. 6(c). Note that the THz energy with a 800-nm pump wave can escape early our numerical box (2.4 × 2.4 mm 2 ), so the previous evaluation is performed at the maximum of THz energy. Performances in the THz gain factor with an increased pump wavelength in collimated propagation thus appear weaker than in focused geometry. When adding the Raman-delayed nonlinearity, the available power contributing to the instantaneous Kerr response is subcritical (∼ 0.7), which prevents the two-color pulse from self-focusing and exceeding the ionization threshold. As a result, no plasma generation takes place and only a residual THz emission occurs due to four-wave mixing. Figure 7(c,d) depicts THz fields propagated over several tens of cm and their spectra for 800-nm and 1600-nm pump pulses. The selected distances correspond to the range of maximum THz energy shown in Figs. 7(a,d), bottom. One clearly sees that the Raman-delayed response contributes locally to diminish the THz conversion efficiency, because the pump dynamics is changed and in particular the peak intensity is reduced. In addition, a longer pump wavelength promotes the formation of a supercontinuum linking the tail of the fundamental pulse spectrum to the THz spectrum [ Fig. 7(d)], which can justify a stronger influence of the current component J A . THz fields with ∼ 0.1 GV/m amplitudes are achieved at both pump wavelengths with a clear amplification at z = 60 cm from the 1600-nm pump.

Generalization for different medium parameters
We now test our previous findings for more complex ionization rates applying to the two major air species and consider different bound electron responses.
We first choose the same f -number ∼ 14 as in Clerici et al.'s experiments [15]. For numerical reasons we limit the initial beam width to w 0 = 150 µm for a focal length f = 2.5 mm. Our two-color pulses have 400 µJ in energy with 5.2 % being injected into the second harmonic. These simulations include both Kerr and Raman nonlinearities and our selected THz window is again 80 THz. The simulations use a temporal step ∆t = 75 attoseconds and a transverse resolution of ∆x = ∆y = 0.88 µm. For completeness we also included ionization of nitrogen molecules using either a QST rate (U i = 15.6 eV, Z * = 1) or a field-dependent PPT rate for two species, adopting Talebpour et al.'s charge numbers Z * O 2 = 0.53, Z * N 2 = 0.9 [36]. When using a QST rate for both O 2 and N 2 molecules, the peak intensity and electron density reach 650 TW/cm 2 and 2.7 × 10 19 cm −3 (complete single ionization) near focus with a 800-nm pump, respectively [see Fig. 8(a)]. With the PPT rate, the effective charge numbers being less than unity promote weaker ionization rates [34], which increases the maximum pulse intensity at focus. Again a complete single ionization of both molecular species is achieved at 800 nm.  Fig. 8(b)]. Differences due to the choice of the ionization model are limited. A λ -dependent scaling of the THz yield appears closer to λ 2 than λ 4 . Figure 8(c) displays the spectra at focus. It is interesting to observe that the computed THz spectra now shifts their maximum to ν ≃ 30 THz, and not to ν = 5 THz as reported in Ref. [15]. This discrepancy may be attributed to the fact that the ABCD technique THz energy yield for ν ≤ 80 THz as a function of the pump wavelength for a focused beam and two species ionized with the QST rate (blue curves) and with an instantaneous PPT rate (green curves). Scaling curves in λ α shown as dashed curves are evaluated through least-square fitting (see legend). The solid red circles recall the experimental data points of Ref. [15]. (c) THz spectra computed at focus. Inset shows maximum THz fields. used in [15] can barely measure frequencies above ∼ 20 THz (see also [12]). Inset details the growth of the maximum THz electric field, which increases quasi-linearly with the pump wavelength. Comparing Figs. 6(c) and 8(b) demonstrates how sensitive the gain curves and THz spectra can be when one varies the laser parameters in focusing geometries.
Finally, to evaluate the influence of the nonlinearity coefficients, we present in Fig. 9 the peak intensity and plasma density, THz energy, spectra and fields of the same femtosecond pulses as in Fig. 7 subject to stronger self-focusing with higher Raman-delayed responses, n 2 = 3.79 × 10 −19 cm 2 /W, x k = 0.79 as recently measured in [39, 40] for 800-nm pump pulses. With 1600-nm pump pulses, following Ref.
[41], we selected the Kerr index values n 2 = 3.72 × 10 −19 cm 2 /W, x k = 0.78. For completeness, we employed an instantaneous PPT rate with effective charge numbers Z O 2 = 0.53 and Z N 2 = 0.9 [36]. Unlike in Fig. 7, the 1600-nm pump is here able to trigger a self-focusing sequence, starting with an input power ratio over critical equal to 3.1 that corresponds to an effective power ratio of about 2.5. Figure 9(b) shows the THz energy evolving with the propagation distance. Over the plasma zone (0.1 ≤ z ≤ 0.7 m), the laser-to-THz conversion efficiency with a 1600-nm pump wave compared to a 800-nm pump is smaller than in the previous filamentary configuration [see green symbols + in Fig. 6(c)], which can be attributed to the larger drop in the peak plasma densities reached in Fig. 9(a). The maximum THz yield achieved with the 1600-nm pump is smaller than that reached in Fig. 7. Figures 9(c,d) detail the THz fields and spectra at the distances of maximum THz energy emission. Computed for low frequencies ≤ 80 THz, on-axis THz fields created with the 800-nm pump prevail. On the whole, the main trend reported about Fig. 7 is retrieved: Compared to the THz gain factors achieved in focused geometry, doubling the pump wavelength for two-color filaments propagating over meter-range distances in air may not significantly increase the THz conversion efficiency.

Conclusion
In summary, we have theoretically studied the influence of long pump wavelengths belonging to the range 0.8-2 µm in THz emissions caused by two-color laser pulses through air photoionization. We also cleared up the action of Raman-delayed and instantaneous Kerr nonlinearities of air molecules on the laser-to-THz conversion efficiency. Optical nonlinearities contribute to THz generation, in particular prior to ionization, but rotational Raman scattering leads to weaken the THz energy yield. At clamping intensity, direct laser-to-THz conversion via fourwave-mixing is weak compared to the photocurrent mechanism. However, Kerr-induced prop-agation effects such as cross-phase modulation and self-steepening have significant impact on THz generation. Furthermore, increasing the pump wavelength can dramatically enhance the THz energy yield. We showed that the THz energy gain factor cannot be quantitatively formulated with a simple power law in λ α due to the influence of the relative phase between the two colors and their pulse envelopes. However, powers of growth rates between 2 and 5 can be justified from the local current model, mainly depending on the relative phase between the two colors. Scalings in ∼ λ 2−3.5 have been extracted in focused propagation geometries through 3D comprehensive UPPE simulations, which faithfully reproduce experimental measurements of THz pulse energies. Similar gain factors can, however, barely be reached in filamentation geometry that clamps the intensity at smaller values and features much lower peak plasma densities at longer wavelengths. These results should help anticipate the THz gain factors achieved with mid-infrared laser systems used in future experiments.