Influence of tunnel ionization to third-harmonic generation of infrared femtosecond laser pulses in air

Here we present an experimental as well as theoretical study of third-harmonic generation in tightly focused femtosecond filaments in air at the wavelength of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.5 \,\upmu \hbox {m}$$\end{document}1.5μm. At low intensities, longitudinal phase matching is dominating in the formation of 3rd harmonics, whereas at higher intensities locked X-waves are formed. We provide the arguments that the X-wave formation is governed mainly by the tunnel-like ionization dynamics rather than by the multiphoton one. Despite of this fact, the impact of the ionization-induced nonlinearity is lower than the one from bound–bound transitions at all intensities.

Here we present an experimental as well as theoretical study of third-harmonic generation in tightly focused femtosecond filaments in air at the wavelength of 1.5 µm . At low intensities, longitudinal phase matching is dominating in the formation of 3rd harmonics, whereas at higher intensities locked X-waves are formed. We provide the arguments that the X-wave formation is governed mainly by the tunnel-like ionization dynamics rather than by the multiphoton one. Despite of this fact, the impact of the ionization-induced nonlinearity is lower than the one from bound-bound transitions at all intensities.
Third harmonic generation (THG) and laser frequency tripling is one of the most important nonlinear optical phenomena. Its application areas range from development of coherent radiation sources in ultraviolet (UV) and vacuum UV (VUV) ranges 1,2 , nonlinear media and plasma characterization 3,4 , nonlinear spectroscopy and microscopy 5,6 to temporal laser pulse characterization 7,8 , white light filament formation and remote atmosphere sensing 9,10 .
THG in air is accompanied by plasma generation and filament formation. Kerr self-focusing leads to a beam collapse which is arrested by plasma defocusing and nonlinear losses 11 . An important part of this filament formation picture is a description of ionization process which is very often made in the terms of multiphoton dynamics. In multiphoton ionization regime, as contrasted to tunneling one, ionization takes place on slow time scale comparable to the pulse envelope and ionization rate grows as some power of intensity, whereas in the tunneling regime the dynamics inside the laser cycle plays the key role, and ionization is described as an exponential function of instantaneous intensity. The key quantity defining which of these two seemingly very different regimes plays a role, is so called Keldysh parameter γ = 2πτ I /T , describing the characteristic time at which the ionization occurs τ I = √ 2I P /E , in relation to a single optical cycle T (here the expression for τ I is written in atomic units, I P is the ionization potential and E is the electric field). The multiphoton range takes place for γ ≫ 1 , whereas for γ ≪ 1 the tunnel ionization is more appropriate description. For the filaments obtained from femtosecond Ti:Sapphire laser pulses with the central wavelength of about 800 nm, intensity clamping takes place at intensities of the order of 10 13 W/cm 2 , which corresponds to γ ≈ 2 , at the boundary of the multiphoton range. At longer pump wavelengths, one can approach this boundary ( γ = 1 ) or even cross it 12 . Such long wavelength filaments can be obtained using pump from high power optical parametric amplifiers, and are actively studied nowadays.
One of the major difficulties arising in the intermediate region γ ≈ 1 is that neither tunnel nor multiphoton formulas work well. Nevertheless, there are few ionization models which are able to fill the gap 13,14 . In particular, in the work of Yudin and Ivanov 14 an approach was developed which allowed to describe not only the scaling of ionization with the driving field intensity correctly in tunnel, multiphoton and transient regions, but also both the fast intra-cycle tunnel-like and slow intercycle multiphoton-like dynamics of ionization.
In the tunnel regime, ionization-based nonlinearity creates also harmonics of the fundamental, which are often called Brunel harmonics 15,16 . In the multiphoton regime, the Brunel harmonics are not generated, because in this case the ionization nonlinearity is "too slow" with respect to the cycle period of the third harmonic radiation 15 . Ionization-based radiation mechanism is related to transitions of electrons which were bounded Scientific Reports | (2020) 10:17437 | https://doi.org/10.1038/s41598-020-74263-x www.nature.com/scientificreports/ inside the atom before the pulse, but remain free after it. This is in contrast to transitions giving impact into χ (3) nonlinearity, corresponding to electrons which remain bound before, after, and during the pulse. To mention here is also the transitions resulting in so called high-harmonic generation (HHG), which start from bound state and end there, but the electron spends some time outside the atom. Third harmonic generated in filaments has impacts from the χ (3) -based nonlinearity as well as from ionization-based nonlinearity. Theoretical estimations 17 show that, whereas for many gases at 800 nm wavelength χ (3) -based mechanism should dominate, situation changes for longer wavelengths. In particular, in experimental work 16 arguments were provided, that at the mid-infrared wavelengths around 4 µm Brunel mechanism dominates for several harmonics, presumingly including the third one.
In general, there are number of works which study THG induced by infrared (IR) laser pulses 1,3-5,9,18-24 . One of the important pecularities by filament propagation is a generation of so-called X-waves. An X-wave is a multicolor conical wave with each cone travelling at the same group velocity. Such X-waves were first suggested to appear in the fundamental harmonic and are explained by spatio-temporal reshaping governed by dynamical nonlinear phase-matched four-wave mixing processes 25,26 . Later on, such characteristic X-shaped distributions were observed in the third harmonic [22][23][24] .
In this report, we present results of experimental and theoretical investigation of THG in air by IR femtosecond laser pulses at 1.5 µ m. Theoretically, we approached the process using two models: the multiphoton one 27 , which included only slow inter-cycle ionization and the nonadiabatic tunnel-based Yudin-Ivanov one 14 , demonstrating sub-cycle ionization features (see above). Comparison of the results of the two models with experiment shows that the tunnel-based model explains the observed experimental data much better. In particular, the tunnel-based model predicts higher clamping intensities, stronger influence of ionization on the dispersion of the fundamental and third harmonics and, as a result, a stronger divergence of X-waves. Moreover, in the Yudin-Ivanov model we compared the impacts of the bound-bound and χ (3) nonlinearities to the third harmonic.

Methods
Experimental set-up. For the experiments we have used a 1 kHz repetition rate femtosecond Ti:sapphire chirped pulse amplification laser system (Legend elite duo HE+, Coherent Inc.), delivering 35-40 fs (FWHM) light pulses centred at 790 nm with maximal pulse energy of 8 mJ. The laser output was used to pump an optical parametric amplifier (OPA, Light Conversion, Inc.), which produced the tunable infrared (IR) signal and idler wave pulses in the wavelength range between 1.2 and 1.6 µ m and 1.6 and 2.4 µ m, respectively. The duration and single-pulse energy of the IR signal wave pulses was about 40 fs and 1.3-1.75 mJ, respectively. The beam width at FWHM was about 4 mm. The IR signal wave pulses were focused in air by the lenses of various focal lengths, which resulted in a visible third harmonic generation.
Far-field patterns of third harmonic radiation were observed on a white screen placed at about 1 m beyond the focal plane of the focusing lens and registered by a CCD camera (Fig. 1). The spectra of this radiation have been recorded with a help of a fibre spectrometer sensitive in the visible and IR spectral ranges. Frequency-angular distributions of third harmonic pulses were recorded with an imaging spectrometer Andor SR-500i-D1 by placing its entrance slit in the focal plane of an objective lens of the short focal length (10-20 cm). Under such lens and slit configuration the spectrometer disperses the light energy across a wide spectral range and produces a twodimensional image at the output imaging plane (the vertical coordinate is proportional to the light direction). Assuming the axial symmetry of the emissions and by placing the CCD camera in the output imaging plane of the spectrometer, the frequency-angular distributions can be captured in a single shot.

Numerical modelling of THG in air.
Here, we describe two mathematical models of THG in air, which we call "MPI model" and "YI model". The MPI model is based on the slow varying envelope approximation (SVEA) for the field, and accounts only the multiphoton ionization of air, neglecting the tunnel ionization. In contrast, the YI model is free from SVEA, that is, formulated directly for the fast optical electric field, and uses the formula www.nature.com/scientificreports/ of Yudin and Ivanov (YI) 14 , which takes into account both slow multiphoton and fast tunnel components of the ionization dynamics. In the MPI model, plasma density depends on the light intensity and the Brunel harmonics are not generated. Third harmonic generation is predefined only by the third-order Kerr term. Because of SVEA, two separate governing equations for the fundamental harmonic (FH) and third harmonic (TH) are used. On the other hand, the YI model allows generation of the third harmonic also via plasma-based nonlinearity. In both models, so called slowly evolving wave 28 and paraxial approximations were used.
MPI model. In the MPI model, we solve the following governing equations for the fundamental harmonic (FH) and third harmonic (TH) waves: where E 1,3 (�, β x , β y , z) stands for the Fourier transform of the complex slow envelopes E 1,3 (t, x, y, z) of the FH (index 1) and TH (index 3) of the fast optical field E(t, x, y, z). Here z is the propagation distance, x and y are the transverse coordinates and t is time. � = ω − ω j0 , where ω j0 = 2πc/ j0 is the central cyclic frequency of the jth wave ( j = 1, 3 ) and ω is the cyclic frequency. j0 is the central wavelength and c is the speed of light. The first rhs term in Eqs. (1) and (2) describes free-space propagation in air. k 1z and k 3z are the z-projections of the wavevectors centered at the corresponding frequencies, e.g.
x + β 2 y )/(2k 1,30 ) − �/u 10 and u 10 is the group velocity of the fundamental wave in air. Here, k(ω) = ωn( )/c , = 2πc/ω is the wavelength, n is the wavelength-dependent refractive index determined by the Sellmeier equations for air 29 , and k 1,30 is the wavenumber at the central frequency ω 1,30 . The shift �/u 10 is introduced to consider both harmonics in the reference frame of the FH pulse.
The nonlinear terms are given by: The first term, P Kerr , includes self-and cross-phase modulations. It is the Fourier transform of the following expressions: where n 2 is the nonlinear refractive index of air, I j = cε 0 |E j | 2 /2 is the intensity ( j = 1, 3 ), ε 0 is the vacuum permittivity. P (a) and P (b) are the Fourier transforms of ε 0 n 2 E 3 E * 2 1 /2 and ε 0 n 2 E 3 1 /6 , respectively. The last term generates the third harmonic due to the χ (3) nonlinearity.
The plasma terms P (pl) 1,3 are related to the plasma density described by the equation: where ρ (n) is the normalized to the neutral density ρ 0 plasma density and W(t) is the generation rate. In MPI model, The plasma is generated both by FH ( j = 1 ) and TH ( j = 3 ). Two main constituents of air-oxygen (O) and nitrogen (N)-are taken into account. Coefficients Ŵ were calculated using the ADK-formula taken from 27 (in 27 , the misprints were corrected). The intensity dependence of the ionization rate was fitted by a line in the logarithmic scale and Ŵ was found, Fig. 2. The linear approximation was performed with the formula log 10 (W ADK ) = log 10 (Ŵ) + K log 10 (I) . K is the photon number needed for the multiphoton ionization. The coefficients are given in Table 1.
For the plasma term, the following formula is utilized: where 30 γ r = 1/200 fs −1 and U O and U N are the ionization energies of oxygen and nitrogen, respectively. These ionization potentials are also utilized when calculating W O,N (t).
The input FH beam. We assume no TH at the input and the FH is focused by a lens with the focus length f placed at z = 0 . Then, at z = 0 , the initial conditions read: for the MPI model, and E(t, x, y, z = 0) = E 1 (t, x, y, z = 0)e −iω 10 t for the YI model. Here, r 0 is the beam radius and τ is the pulse duration. E 0 is the peak amplitude and the corresponding input peak intensity is noted as I 0 .

Results
Equations (1,2) and (13) were simulated assuming the cylindrical symmetry of the problem, thus the fast Hankel transform for the space-domain was utilized 31 . For the integration, we used the Euler numerical scheme. For the initial conditions, FH at 10 = 1.5 µ m was assumed, with r 0 = 5 mm and τ = 30 fs. In the MPI model, the space window (0, 2] × r 0 was taken, with a grid consisting of 1500 points, as well as the temporal window [−7τ , 7τ ) , with a uniform grid containing 128 points. For the YI model, the same spatial grid, and a temporal uniform grid containing 512 points spanning the temporal window [−4τ , 4τ ) were used.  (Figs. 7, 8, 9). At the larger value of f and lower input energy, two maxima in the TH frequency-angular spectrum are observed, see Fig. 3a. These maxima mean the appearance of the one-color ring in the TH angular spectrum. The radius of the ring was found from the longitudinal phase-matching condition ( k 30 cos θ p − 3k 10 = 0 ) and is equal to θ p = 3.4 mrad 22,32 (see dashed lines in Fig. 3a). In this case, the plasma density is low and the intensity clamping is not yet achieved. From the curves in Fig. 3b we see that after passing the focus, the peak intensity of the TH wave (red dashed line) decreases faster than the intensity of the FH wave (blue dashed line) since the ring is formed in the TH spatial profile and consequently the central spot vanishes. At the focus, the TH beam obeys the central spot and its width is smaller than the width of the FH beam.
For the intensities shown in Fig. 3 both MPI and YI models give the same result. This is not surprising, since the peak intensity approaches 10 11 W/cm 2 in this case, that is, we are clearly in the multiphoton regime, where both MPI model and YI models being valid. The Keldysh parameter γ is equal to 1 at intensity around 3 × 10 13 W/cm 2 for the fundamental harmonic, so we expect that as we approach this intensity, MPI and YI models start to deviate, and the tunnel ionization which is not accounted in the MPI model, starts to play a role. www.nature.com/scientificreports/ When the input energy increases, the ring merges into a single central spot as seen in Fig. 4a-d, which is observed for both models, although for somewhat different input intensities. In the MPI model, the plasma density becomes saturated at high focal intensity, Fig. 4e, whereas the saturation in the YI model is not so prominent (compare black solid lines in Fig. 4e,f). As a result, the intensity clamping is observed in the MPI model. After passing the focus, the peak intensity of the TH beam falls with the same speed as the peak intensity of the FH beam, compare blue and red dashed lines in Fig. 4e,f. The experimentally recorded spectra of the TH show the same evolution from the ring-shaped spectrum to the central spot, as seen in Fig. 5. In Fig. 6, the experimental TH patterns are shown, which show also the merging. This merging of the ring into a central spot may be attributed to cross-and self-phase modulations of the TH wave. To check this, we repeated the simulation of YI model without the χ (3) term and the ring structure did not disappear (data not shown). We also checked the stability of the results with respect to the change of Ŵ and input intensity I 0 , e.g. we repeated the simulations first with 10 times larger Ŵ 1N and second with I 0 increased by 1 percent. In both cases, the results are undistinguishable from the initial data.
At the smaller focal length, f = 20 cm, larger peak intensity can be achieved. The TH spectra become noticeably blue-shifted, see Fig. 7a,c,e. In these figures, only the results for the YI model are shown and the comparison  www.nature.com/scientificreports/ with the MPI model will be provided below. The intensity clamping is observed, e.g. in Fig. 7b,d,f at the maximum FH intensity around 7 × 10 13 W/cm 2 for all three I 0 values. From the dashed lines in Fig. 7b,d,f we see that TH and FH intensities of the beams follow the same trend. We think that the signal at 400 nm (Fig. 7e) is a result of the cascaded process of mixing of TH at 500 nm and FH at 1.6 µm. This wavelength component in FH appears due to nonlinear broadening of the pump (see Fig. 8b,d).
In Fig. 8, one can see that the locked FH and TH X-waves travelling at the same group velocity are formed 22 . This result is obtained in both YI-and MPI-models, but the group velocities are different in both models (cf. Fig. 8a,b and Fig. 8c,d). In Fig. 8, we also depicted the X-wave dispersion curves by the dashed lines. The dispersion law is given by: . V is the group velocity of the X-wave. The group velocity V was varied by the step of 10 4 m/s to get the best coincidence with the numerically calculated contours. The resulting θ( ) following from Eq. (20) is depicted in Fig. 8 by dashed lines. Good coincidence with our experiment and also 22 is observed in the case of the YI model. In particular, we found the same group velocity V = 3.01 × 10 8 m/s as in 22 . In contrast, the group velocity obtained from the MPI model deviates from this value: V = 3.0017 × 10 8 m/s. This fact is an evidence of the influence of tunnel ionization on the X-wave formation. The experimental TH spectra are shown in Fig. 9. Comparison with the numerical results of the YI model in Fig. 7 shows a good agreement. To facilitate such comparison, the dispersion curve resulted from the YI model is shown on the experimental graph (Fig. 9c) by a dashed line. We note that in the theoretical spectra, the angle was calculated as θ = k ⊥ /k 30 , where k ⊥ = β 2 x + β 2 y is the transverse projection of the wavevector. The formula is only approximately correct. So, the experimental angle θ is slightly smaller than the theoretical one at the wavelengths shorter than 30 = 500 nm.
In order to examine the influence of tunnel ionization we provide an exemplary picture of the free electron density (curve 1) and ionization rate (curve 2) shown in Fig. 10a, at the center of the beam for the case of Fig. 8c,d. Here, the oscillations are prominent and the sub-cycle tunnel dynamics are dominant comparing to the multiphoton ones. The question may arise if the Brunel harmonics provided by plasma dominate also the third harmonic generated by χ (3) nonlinear term? To answer this question we switched off the χ (3) term in the YI model and repeated the simulation for the parameters of Fig. 7f. The results are presented in Fig. 10b. Comparison of these two figures shows that the TH intensity of the plasma radiation is by around two orders lower than the full TH intensity.

Discussion and conclusions
In conclusion, evolution of the third-harmonic spectrum from the ring-type to locked X-waves by increasing the focal pump intensity was analyzed both theoretically and experimentally. For the theoretical consideration, two models were used, the MPI model based on multiphoton ionization and YI model, working in multiphoton, tunnel as well as in intermediate ranges. We observed several distinct stages of evolution of TH beam. At low www.nature.com/scientificreports/ intensities, TH forms a ring, which radius is determined by the longitudinal phase-matching. This ring merges into a single central spot at higher intensities due to self-and cross-phase modulations. At even higher intensities, locked X-waves in the third-and fundamental harmonics are generated and co-propagate with the same group velocity. We found two different values of the group velocity of the X-waves: one from the MPI model and another one from the YI model. We observe that the latter value coincides with the experimental one and it is in agreement with the previous experimental work 22 , in contrast to the value given by the MPI model.