Methodical inaccuracy of the Z-scan method for few-cycle terahertz pulses

Modern sources of THz radiation generate high-intensity pulses allowing to observe nonlinear effects in this spectral range. To describe many nonlinear effects theoretically, it is necessary to know the nonlinear refractive index coefficient of optical materials. The work studies the applicability of the Z-scan method to determine the nonlinear refractive index coefficient in the THz frequency range for few-cycle pulses. We have discussed the correctness of the known Z-scan method for calculating the nonlinear refractive index coefficient for broadband THz radiation regarding number of cycles pulses have. We have demonstrated that the error in determining the nonlinear refractive index coefficient is always greater than 70% for true single-cycle pulses. With the increase in the number of oscillations to the measurement error shows strong dependence on the sample thickness and can vary from 2% to 90% regarding the parameters chosen. The fact that such radiation dispersion length is commensurate with the nonlinear length or even less than the latter results in the discrepancy mentioned. It is demonstrated that the decrease in the sample thickness leads to the reduction of the nonlinear refractive index coefficient determination error, and this error is <2% when the ratio between the sample thickness and the pulse longitudinal spatial size is ≤1. This can relate to the fact that the nonlinear effects in such a thin sample occur faster than the dispersion ones.

Currently, the study of the terahertz (THz) frequency range of optical radiation is one of the developing and promising scientific fields. The reason for active performing in the area lies in the broad applications of THz radiation in biology 1,2 , medicine 3,4 , security systems 5 , non-destructive evaluation 6,7 , wireless information transmission 8 . Until recently practical devices and systems were mainly developed using linear effects of THz wave optics 9 . At present, there are sources of pulsed THz radiation with energy in a single pulse up to 1 mJ have appeared 10,11 , i.e., the peak values of the electric field can reach values up to hundreds of MV/cm 12 . This allowed to move to profound research in the field of nonlinear THz optics [13][14][15] . High-energy sources were constructed using singe-and-double-color filamentation [16][17][18] , as well as using femtosecond pulses rectification in LiNbO 3 crystals 10,11,14,19 .
First experimental observations of such a classical nonlinear phenomenon as self-phase modulation in THz nonlinear optics 20 demonstrated significant changes in the refractive index of materials. The theoretical work 21 predicts that values of the fast-response nonlinear refractive index coefficient (n 2 ) of materials in the THz spectral range are several orders of magnitude higher than its values for ultrashort pulses of the visible and near-IR spectral ranges. Moreover, nonlinear effects are experimentally observed in semiconductors 22,23 in the THz frequency range, for which the theory was preliminarily built 24,25 using nonlinear THz optics. Regarding rotational states of gas molecules, it was shown, both experimentally and theoretically 26,27 , that the coherent control of THz radiation by populations of rotational sublevels and their relaxation is possible. The work 28 investigated nonlinear four-wave mixing in gases between intense ultrashort optical pulses and THz fields.
There are several techniques for measuring n 2 in IR and optical spectral ranges. One of the most common ones is the Z-scan method 29 . Used for monochromatic radiation, the technique can also be applied to femtosecond radiation, which radiation spectrum is quite wide 30 . Recently, the first works on preliminary experiments estimating n 2 coefficient using the Z-scan technique in THz spectral range clearly showed the presence of a significant nonlinear effect 31,32 . However, the application of the Z-scan method raises several questions since the spectrum of pulsed THz sources is even wider than femtosecond one 33 and its electric field can have few cycles only 34,35 . Thus, in the experimental works mentioned estimates of the nonlinear refractive index coefficient were made only,

Results
To analyze the applicability of the Z-scan method for broadband few-cycle THz pulses the results of the numerical simulation were compared with the analytical model of the method for monochromatic radiation (see Methods). It should be mentioned that the pulse duration does not appear in the equations, and therefore, the analytical curve does not depend on this pulse characteristic. The numerical simulation of the Z-scan method was conducted using the equations of intense light pulse propagation in a dielectric medium with normal group dispersion and nonresonant nonlinearity 36 (see Methods). Using this model the propagation of a spherical THz beam at two focal lengths was numerically simulated. For different positions of the nonlinear medium sample on the optical axis of propagation, the resulting THz pulses were detected through the aperture. Figure 1 shows the results of the numerical simulation of the typical Z-scan method curves (see Methods). They represent the dependence of the transmission through the closed aperture T on the sample position z for the following parameters: central wavelength of the THz radiation λ = 0.3 mm, period of oscillation T 0 = 1 ps, peak intensity in the caustic (a) I 0 = 3.1 ⋅ 10 8 W/cm 2 , and (b) I 0 = 8.3 ⋅ 10 8 W/cm 2 , durationτ 0 = 0.3, 1, 10T 0 , thickness of the ZnSe L = 0.3 mm. The value of the nonlinear refractive index coefficient n 2 = 4 ⋅ 10 −11 cm 2 /W 31 . As seen, the obtained curves in Fig. 1 qualitatively correspond to the Z-scan curves. For the central wavelength λ 0 = 0.3 mm the duration τ = 0.3 ps corresponds to the single oscillation of electromagnetic field 37 . In the simulation in Fig. 1(b) a different intensity value is used due to the fact that in the case of a true single-cycle pulse its waist diameter is 2 times less than the one of multi-cycle pulses and the focusing of such a pulse leads to its transformation into a 1.5-cycle pulse 37 . It happens since the central frequency shifts to higher values (see Methods) in the case of a true single-cycle pulse.
Noteworthy, the shorter the pulse duration is the less its peak-to-valley ratio corresponds to the one the analytical curve has, and therefore, the higher the inaccuracy of the obtained results of n 2 estimations is. It can be concluded that the Z-scan method introduces the most significant inaccuracy in the determination of n 2 when it comes to pulses of less than 2 oscillations. As mentioned above, for a true single-cycle pulse a change in its duration in the caustic and decrease in waist size are observed. Additionally, a shift of the central frequency to the higher values occures. All the facts mentioned result in the peak intensity in the caustic increase by 2.7 times (see Methods). Despite this fact, the peak-to-valley ratio of the latter is smaller than even for a two-cycle pulse, as seen in Fig. 1. The analytical expression gives the value of 4.0 ⋅ 10 −11 cm 2 /W for n 2 , whereas the numerical simulation curve leads to the value of 1.15 ⋅ 10 −11 cm 2 /W. Thus, the inaccuracy of the method, in this case, is more than 70%. Figure 2 illustrates the dependence of the magnitude of the error of n 2 calculation from the differential curve obtained by numerical simulation of the Z-scan method on the number of cycles of the THz pulse T 0 in the case of a fixed sample thickness L = 0.3 mm for different values of the wavelength λ 0 = 0.3 mm and 0.4 mm.
The results show that with a decrease in the number of pulse cycles when the pulse has less than 2 oscillations, the magnitude of n 2 measurement error crucially increases. The highest magnitude of the error corresponds to a true single-cycle pulse. As seen, the magnitude of the error practically does not change for different wavelengths. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
The reason for the discrepancy between the simulation results for pulses with less than 2 oscillations and the analytical calculations for monochromatic THz radiation may result from the fact that the dispersion effects on such a pulse are consistent with the nonlinearity 38 . For the few-cycle pulses, the dispersion length L disp commensurates the nonlinear length L nl or even less than the latter (see Methods) 38 . The calculation of the dispersion and the nonlinear lengths for few-cycle pulses was performed (see Methods). For example, in the case of the true single-cycle pulse with the central wavelength λ 0 = 0.3 mm they are L disp = 7.81 mm and L nl = 9.07 mm correspondingly. These values are larger than the sample thickness L = 0.3 mm. We assumed that an increase in the sample thickness would lead to a larger error in n 2 estimation while its decrease would lead to the smaller one. The latter can be related to thinning of the sample which causes the nonlinear properties of the medium to occur faster than the dispersive ones. To confirm this assumption, we performed a simulation for different wavelengths λ 0 = 0.2-0.6 mm and different sample thicknesses L = 0.06-4.5 mm. To colligate the results for all values of central wavelengths, the general parameter L/x was introduced. It is the ratio of the sample thickness L to the spatial size of the pulse and describes the number of pulses fitting inside the sample, was introduced. The dependence of the magnitude of n 2 measurement error on the ratio between the sample thickness and the spatial size of the pulse L/x is the same for all the wavelengths and is shown in Fig. 3.
It is seen that for broadband THz single-cycle pulsed radiation, the magnitude of n 2 measurement by the Z-scan method error increases with the L/x ratio. For the sample thickness corresponding to the L/x ratio less than 10, the rapid increase in n 2 measurement error is observed. For L/x = 10 n 2 measurement error is equal to more than 70%. Otherwise, when L/x > 10, the slope of the curve is smaller and approaches the values close to the saturation (90% approximately). Accordingly, the most accurate results are achieved in the case of the lowest value of the last relation. Since the size of the sample cannot be infinitely reduced, the best option is when the ratio ≤ L x / 1 and the Z-scan method has a negligibly small error of n 2 measurement. This may be explained by the fact that nonlinear effects occur faster than the dispersion one.  www.nature.com/scientificreports www.nature.com/scientificreports/

Methods
Analytical model of the Z-scan method for quasi-monochromatic radiation. The work 29 proposes to calculate the Z-scan curve differential (the transmission through the closed aperture T) using the following formula: 2 is the aperture linear transmittance; the transmitted power through the aperture gives T r a 0 0 0 0 The following values are used in these formulas: absorption coefficient α = 0.85 cm −1 , sample length L = 0.3 mm, central frequency of the radiation ν 0 = 1 THz (λ 0 = 0.3 mm), beam waist radius w 0 = 1.4 mm, aperture radius r a = 1.5 mm, radius of the THz beam w a = 12.5 mm, intensity of the THz beam in caustic I 0 = 3 ⋅ 10 8 W/cm 2 , nonlinear refractive index coefficient n 2 = 4 ⋅ 10 −11 cm 2 /W. This value of n 2 was obtained in experiment previously 31 .
Numerical simulation of the Z-scan method. In the paper 21 analytical model for calculating the non-resonant vibrational contribution to the nonlinear refractive index n 2 (Kerr coefficient) is developed. It is demonstrated that the value of this contribution in the THz spectral region can exceed the value of n 2 in the visible and near-IR spectral ranges by several orders of magnitude. Also, authors show that for the case of ultrashort optical pulses, including intense picosecond THz pulses, the dominant source of nonlinearity tends to be the low-inertia ones where nonlinear mechanisms are based on the nonlinear response of each molecule to the radiation field. Noteworthy, for the pulses in the THz range, one expects the dominant mechanism of nonlinearity to be associated with anharmonic vibrations of the lattice, unlike visible and IR, where nonlinearities has the electronic nature. These mechanisms are well-described with nonlinear refractive index coefficient n 2 as a material characteristic. These findings have now been confirmed experimentally 31,39,40 . Thus, the use of the Z-scan method is justified allowing to calculate the Kerr nonlinearity in the THz frequency range. Work 21 shows that n 2 dispersion of the vibrational nature of crystals can usually be neglected in the THz frequency range. n 2 dispersion can significantly increase in the region of two-photon resonance with the vibrational mode of the crystal or natural oscillation of the molecule in liquid. For example, in water, oscillations that determine its vibrational nonlinearity in the far IR range of the spectrum have frequency of 15.9 THz while central frequency of used THz radiation is 0.75 THz 39 . Accordingly, up to the frequency values of 8 THz a small dispersion of a vibrational nature in water is expected. Therefore, in our work, we assume that n 2 is considered as a constant value.
It is well known that the essence of the Z-scan method consists of induced narrowing and expansion of an intense spherical light beam when a nonlinear medium moves along the axis of its propagation and passes through the focus 29 . In this case, the nonlinear medium plays the role of a thin lens and its placement in or near the focus leads to a minimal change in the distribution of the beam field in the far field. The resulting characteristic Z-scan curve represents the peak and valley of the transmission of the nonlinear medium, from the magnitude of their difference it is possible to determine the value of the nonlinear refractive index coefficient. In this work, we performed a numerical simulation of the Z-scan method, an illustration of which is shown in Fig. 4(a). It illustrates the propagation of a spherical THz beam at two focal lengths. In this graphical representation, the minimum of the electric field corresponds to the color of blue, and the maximum of the electric field to red. For each position of the ZnSe crystal, the THz pulse was distributed through the air, the crystal itself and then further through the remaining distance in the air to be finally captured at the detector through the aperture. Figure 4(b-i) represents examples of single and multi-cycle pulses and a central wavelength of 300 μm. To represent the transverse size of pulse fractions of wavelengths calculated as r/λ 0 ratio are used as units. A pulse with a duration of 1 ps, in fact, has 1.5 cycles, whereas a true single-cycle pulse has a duration of 0.3 ps. It can be seen in Fig. 4(d) though, that the spectrum of a true single-cycle pulse shifts to the higher value area. Regarding the numerical simulation performed, this fact means that the central wavelength of the radiation becomes smaller and, therefore, the waist radius in the caustic also decreases. In turn, this leads to the peak intensity in the caustic raise by 1.5 times and should cause an increase in the peak-to-valley ratio of the Z-scan curve. A true single-cycle pulse corresponds to the ideal THz pulse 37 . However, generated THz pulses have the form of the first type in real experiments 41 . It should be noted that despite the focal length being much larger than the transverse size of the THz field (which is consistent with the case of the paraxial propagation), a very strong curvature of the wavefront takes place for a true single-cycle pulse. The latest analysis of nonlinear effects for nonparaxial beams of single-cycle THz pulses 42 shows that the paraxial approximation characterizes the dynamics of the transverse field component and strongly focused nonparaxial THz beams well. The appearance and nontrivial dynamics of the longitudinal field component are the difference in this case.
www.nature.com/scientificreports www.nature.com/scientificreports/ The numerical simulation of the Z-scan method was conducted for the following system parameters: sample under the investigation -ZnSe (its thickness L = 0.06-0.3 mm), focal length of the lens f = 50 cm, central wavelength of the THz pulse λ 0 = 0.2-0.6 mm, with period of oscilation T 0 = 0.75-2 ps, pulse duration δt = 0.3-10T 0 , the transverse width of the input beam D = 25 mm, size of the aperture A = 1.5 mm, the peak intensity in the caustic I 0 = 3.1 ⋅ 10 8 W/cm 2 . The dispersion parameters are specified in 43 and are approximated by the formula n 0 (ω) = N 0 + acω 2 , where N 0 = 3 is the refractive index at the "zero" frequency, and a = 6 ⋅ 10 −36 s 3 /m is the dispersion coefficient; n 2 value is taken from the paper 31 as the result of experimental data evaluation and is equal to n 2 = 4 ⋅ 10 −11 cm 2 /W. To simulate the propagation of THz pulse in the air and ZnSe we use the equations of an intense light pulse propagation in a waveguide dielectric medium with a normal group dispersion and nonresonant nonlinearity 36 : where the second item from the left part describes the dispersion of the linear polarization response of the electronic and vibrational nature, the third one describes the nonlinearity of the response of the electron nature, and the term in the right part describes the diffraction of the extremely short pulse. In this equation, z is the sample position measured with respect to the focal plane, Δ ⊥ -Transverse Laplacian, τ = t − n 0 z/c -time in the accompanying coordinate system, is the transverse radius of the input beam, R(x, y) = exp(−ik(x 2 + y 2 )/2f) is the transmission function of the spherical lens with the focal length f, k = n 0 ω/c is the wavenumber, x and y are the transversed coordinates.
During the data processing of an array of numbers characterizing the detected signal at each crystal position, the energy was integrated along the time axis for an aperture with the radius of r 0 = 0.75 mm and t = 2000 fs using the equation: www.nature.com/scientificreports www.nature.com/scientificreports/ It should be noted that few-cycle THz pulses are critically short regarding the number of oscillations. As known, the envelope approach was historically suggested as a method to analyze the evolution of pulses containing a large number of field cycles, therefore, to precisely describe the dynamics of few-cycle pulses additional equations are needed to be introduced. The latter leads to the mathematical model be more cumbersome. The field approach, in turn, allows to express all the necessary information about a pulse, including generation at triple and other high frequencies, in one equation only 36 . Therefore, the key problem of the theoretical study of their propagation laws in a nonlinear media is the need to improve existing 36 and elaborate new approaches to the analysis of the field dynamics and emission spectra associated with extremely short THz pulses features. the dispersion and the nonlinear lengths calculation. To perform the calculation of the dispersion and nonlinear lengths for the few-cycle pulses the following formulas have been used 37 : nl nl 2 is a modification of the refractive index at the central wavelength λ 0 due to dispersion, Δ = ⋅ n nI nl 1 2 2 is a nonlinearly induced change of the optical refractive index, ω = π T 0 2 0 is the central optical frequency, and N 0 is the refractive index at the "zero" frequency, a is the dispersion coefficient, which describe the dispersive properties of the medium: n 0 (ω) = N 0 + acω 2 .

Conclusion
The work has studied the Z-scan method applicability for the broadband pulsed THz radiation of various number of the field oscillations through the numerical simulation. Numerical Z-scan curves for broadband THz radiation of three and more cycles show good correspondence with the analytical curve for monochromatic THz radiation, while for true single-cycle pulses n 2 measurement error occurs constantly and is always high -70% approximately. Furthermore, regardless of the central wavelength of the radiation, an increase in the error of n 2 measurement with an increase in the sample thickness is observed for pulses with 1.5-3 cycles. n 2 measurement error can be as large as 90% in this case. This results from the fact that for such radiation the dispersion length is comparable or even less than the nonlinear length. Meanwhile, n 2 measurement error decreases along with the sample thickness. As a result, it has been shown that an essential parameter of the correctness of n 2 measurement for few-cycle pulses is the ratio L/x, where L is the sample thickness, and x is the spatial size of the pulse. It was found that the magnitude of the error increases together with the magnitude of the ratio. For L/x < 1 n 2 measurement error equals to less than 2%. Therefore, it is recommended to use samples which thickness does not exceed the longitudinal spatial size of the pulse to measure its nonlinear parameters using the Z-scan method. This is explained by the fact that the nonlinear effects occur faster than the dispersive ones in such a thin sample.
This study is useful for working in the field of nonlinear optics and using the Z-scan method to determine various nonlinear parameters of media in different spectral ranges.