## Abstract

We theoretically investigate the propagation of broadband single-cycle terahertz (THz) pulses through a medium with a nonlinear optical response. Our model takes into account non-paraxial effects, self-focusing and diffraction, as well as dispersion, in both the linear and nonlinear optical regimes. We investigate the contribution of non-instantaneous Kerr-type nonlinearity to the overall instantaneous and delayed Kerr effect at the THz frequencies. We show how increasing the nonlinear relaxation time and its dispersion modifies the THz pulse after the propagation through a transparent medium. We also discuss the effect of linear dispersion on self-action during the pulse propagation.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Thanks to the recent developments in the terahertz (THz) coherent radiation sources, THz pulses with higher intensities are becoming routinely accessible [1–3]. These sources have enabled a range of nonlinear THz time-domain spectroscopy (THz-TDS) techniques where nonlinear light-matter interactions can be controlled by temporally trapping the generated field within a single-cycle pulse, and then by coherently releasing the high peak field into the matter. The broadband coherent spectral content of these pulses is desirable in many imaging and spectroscopic applications [4,5].

One of the earliest studies on the far-infrared nonlinear interactions was the observation of the ionization of high-lying Rydberg states in *Sodium* (Na) Rydberg atoms using half-cycle THz pulses [6]. Also, self-phase modulation (SPM) due to the nonlinear response of free electrons to the intense single-cycle THz pulse was reported at the THz frequencies [7]. In [7] the authors performed a nonlinear THz-TDS experiment in transmission configuration and reported on the observation of SPM and saturable absorption in an n-doped bulk semiconductor GaAs, resulting in a 200 fs group delay and 32$\%$ increase in transmission, respectively. Other THz-related works probing the nonlinear optical response have focused on THz-induced carrier multiplication via impact ionization [8–10] and THz saturable absorption and higher-harmonic generation by hot electrons [11–16]. For instance, the authors in [11] exposed an n-type bulk GaN semiconductor to a high-field THz radiation and found coherent emission centered at 2 THz with picosecond decoherence time. They also confirm that this nonlinear response grows superlinearly with the field peak value and saturates as the electric field goes beyond a certain value. They associated this emission with the impurity transitions in the semiconductor material.

Another physical process that can contribute to the third-order nonlinear optical interactions at THz frequencies is molecular vibrations. Since the time scale of the THz pulse is on the order of the relaxation constants of molecular vibrations in most solids, a strong nonlinear polarization can be produced as a result of coherent oscillation of phonons excited by the THz pulse. Previous theoretical and experimental studies focusing on the THz frequency range have shown that, by operating close to the vibrational lattice resonance, one can create a non-instantaneous nonlinear contribution to the Kerr effect in addition to the instantaneous electronic response, which could result in an enhancement of the overall Kerr effect by several orders of magnitude [17–19]. However, it is not clear what the share of the delayed response in the strength of the overall Kerr effect is. The focus of the present study is on the characteristic time scales for different contributions to the overall Kerr effect during the nonlinear optical interactions at the THz frequencies. We do so by adopting a model taking into account both contributions of the Kerr effect corresponding to an instantaneous and retarded third-order nonlinear responses. This model is similar to the conventional model for the Raman effect in silica fibers [20]. Since lattice-driven nonlinearities are slower than the electronic responses, they are modelled in a similar manner to the stimulated Raman scattering (SRS) in silica fiber where there is a long interaction length with optical pulses resulting in important changes in the optical spectrum observed at the fiber output.

In order to accurately model the electric field propagation in the THz range, one needs to use the proper approximations while solving Maxwell’s equations. In the many-cycle regime of optical pulses where there are several carrier oscillations within the envelope of the pulse, the slowly varying envelope approximation (SVEA) is valid [21]. The SVEA neglects the second order spatial derivative of the field in the wave equation and leaves only the first derivative with respect to the spatial coordinate along the propagation direction. For both ultra-short optical pulses and broadband THz transients, where the pulse envelope contains only a few optical cycles, the SVEA is no longer valid and fails to describe the propagation of such pulses [21,22].

In this study, we use a unidirectional pulse propagation equation (UPPE) [23,24] which gives a seamless transition from Maxwell’s equations to an optical pulse propagation model without making the SVEA. The UPPE is a non-paraxial version of the Forward Maxwell Equation (FME) which is itself a paraxial equation. However, both UPPE and FME allow for diffraction in the transverse plane as well as for modelling space-time focusing [25]. This comprehensive model helps with building a rich environment in which one can see and account for the impact of different linear and nonlinear effects on the electric field, as well as on the spatial profile of the THz field during its propagation. This knowledge can be useful in the interpretation of THz-TDS measurements, and it can help with estimating the material parameters when comparing the simulation and experimental results. Both linear and nonlinear propagation effects of a single-cycle THz pulse are accounted for in our model. To quantify the effect of dispersion on the induced nonlinearity, we compare two scenarios: an optical pulse propagating inside (i) a material with strong normal dispersion and (ii) a material with anomalous dispersion.

In section 2, we present the detail of the theoretical propagation model as well as the nonlinear response model. In section 3, we show the modelling results and discussion. First, we look at the electric field time profile of the beams and compare the instantaneous and non-instantaneous contributions to the Kerr effect. Secondly, we study the self-focusing of the beam that arises as a result of an interplay between linear dispersion and nonlinear effects in the THz regime, and, finally, we discuss frequency-resolved diffraction of the beam. The conclusion is presented in section 4.

## 2. Theoretical model

The derivation of the UPPE equations has been done in [23,26] using two different approaches based on the temporal (time-propagated UPPE) and spatial evolution. The spatial coordinate in the latter is usually chosen as $z$ ($z$-propagated UPPE). The initial condition for a time-propagated UPPE is in general the electric and magnetic fields’ amplitudes of a pulse at an initial moment as functions of spatial coordinate. On the other hand, the initial condition for the $z$-propagated UPPE is the electric and magnetic fields’ amplitudes at a fixed location as functions of time. This requires a complete knowledge of the total field at all the past and future times, which makes it hard to implement without simplifying assumptions. Although we use the UPPE solver that applies the actual time-propagated approach based on a direct integration of Maxwell’s equations, here we present an intuitive (with approximations) derivation of the other similar approach ($z$-propagated UPPE) to be able to introduce the UPPE on the terms that we are interested in. This procedure has been given in detail in [23,24].

In order to derive the $z$-propagated UPPE, we start from the most general form of the Maxwell’s equations. Combining the Maxwell-Ampère and Maxwell-Gauss differential equations, one can obtain:

#### 2.1 Non-paraxial regime

In general, the condition $|k|\;>>\;|k_\perp |$ used in the FME given by Eq. (5) represents a fair approximation, especially for the visible and near-infrared ranges. However, in THz frequency range, where the propagating pulse contains a wide-band spectrum, meaning that the spectral pulse width and the central wavenumber are comparable, this condition is often not valid. Hence, in order to get a more accurate picture of the THz beam focusing, it is desirable to go back to the non-paraxial regime where the extent of the angular spectrum of the field is non-negligible compared to the wavenumber of the beam along the propagation axis [28,29]. With this in mind, we take Eq. (2) and rewrite it in a more general form:

Following the procedure similar to that used for deriving Eqs. (3)–(5), the UPPE with no paraxial approximations can be obtained:#### 2.2 Nonlinear response modelling

Let us consider the propagation of a broadband THz pulse in a nonlinear optical material. The induced polarization can be represented as the sum of the linear and nonlinear contributions. The nonlinear polarization can be related to the electric field $E(t)$ and to the time-dependent susceptibility as [20]

Here $\chi _0^{(3)}$ is a scaling factor, and $g(t)$ is the response function given by The first term in Eq. (9) describes the instantaneous Kerr effect, while the second term describes the delayed Raman effect, modeled by a single Lorentzian function, centered at the phonon resonance frequency Here $\omega _R$ is the angular frequency of the phonon resonance and $\delta _R$ is the resonance linewidth. The parameter $\alpha$ ($0\leq \alpha \leq 1$) in Eq. (9) determines the ratio of the Kerr susceptibility to the overall susceptibility including both Kerr and Raman effects. In other words, By setting $\alpha$ to 1, one can eliminate any delayed response and observe the instantaneous Kerr response alone, while $\alpha =0$ gives the delayed Kerr response, only. Using Eqs. (8) and (9), we can study the impact of the parameter $\alpha$, describing the role of vibrational modes in the overall nonlinear response, on the induced polarization and the transmitted THz pulse. Also, in order to obtain the Raman response function ($g_R(\omega )$), we take the reported values for the linear refractive indices (see Fig. 1) and absorption coefficient of the materials at the THz frequencies, and we calculate the $\chi ^{(3)}$ dispersion using the two-level atom approximation [21]. In such a way, we obtain the material-specific response time ($1/\delta _R$) and phonon resonance frequency ($\omega _R$) from the $\chi ^{(3)}$ dispersion curve, and we plug them into Eq. (10).## 3. Results and discussion

We use the gUPPEcore package [30] to solve UPPE describing the propagation of a THz pulse within a dielectric material. In order to be consistent with most of the practical THz pulses, we set the spectral range of our input pulse between 0 and 3 THz. For the simulations, we have included both the linear and nonlinear instantaneous and Lorentz-dispersive effects. Also, we work in a non-paraxial regime, and there is no SVEA. We let the beam propagate through two different materials having linear optical properties of two broadly used semiconductors: silicon (Si) and zinc selenide (ZnSe), and we compare the results. The reason behind this material selection is the totally different linear dispersions of the two materials in the frequency range of interest. Group velocity dispersion (GVD) of both materials is given in Fig. 1. Si shows a weak anomalous GVD in the operation frequency range (0-3.5 THz) [31], while ZnSe exhibits a normal GVD in the same frequency range, together with a strong phonon resonance at $\sim 6.6$ THz [32]. In the following, we systematically study the influence of different mechanisms on the temporal, spatial and spectral profile of the propagating THz pulse.

#### 3.1 THz electric field

We studied the transition of light-matter interaction from the linear ($n_2=0$) to an arbitrarily high nonlinear regime ($n_2 = 1\times 10^{-15}\,\textrm {m}^2/\textrm {W}$) for different strengths of the Raman effect, described by the parameter $\alpha$. For the nonlinear study, we first set $\alpha$ to 1 to consider the instantaneous Kerr contribution only. After that, we set $\alpha$ to 0.9 to add a small contribution from the Raman (non-instantaneous) effect. The intensity of the input pulse was set to $I=1\times 10^{13}\,\textrm {W}/\textrm {m}^2$, while the thicknesses of both samples were set equally to 2 cm. Also, the input pulse duration was set to $\tau =0.5~\textrm {ps}$, and the initial beam waist was $w_0=1.5$ cm. The temporal profiles of the input THz electric field after its propagation through the two different samples are shown in Fig. 2. For Si (Fig. 2(a)), the pulse affected by the instantaneous Kerr effect experiences around 3.4 ps time delay due to the refractive index change of at least $\Delta n=n_2I=0.01$. Moreover, the Raman contribution introduces an additional time shift to the propagating pulse affected by the Kerr nonlinearity. This time shift becomes slightly larger when the response time of the phonons changes from 10 to 100 fs (solid violet line in Fig. 2(a)). The phonon resonance frequency of Si is set to 3 THz. On the other hand, the phonon resonance of ZnSe is calculated to be at 6.6 THz with a linewidth of 1 THz, which gives a response time of $1~\textrm {ps}$. As shown in Fig. 2(b), pulses affected by Kerr and Kerr+Raman effects experience 240 fs and 500 fs time shifts, respectively. Since there is a strong dispersion associated with ZnSe, the time delays in this case are smaller compared to the time delays of the pulses in Si.

#### 3.2 Self-focusing

Self-focusing, which is the spatial analog of self-phase modulation (SPM), depends on the linear dispersion during the propagation of an optical pulse in a material. It is, therefore, necessary to consider the combined effects of SPM and GVD on the pulse evolution. It can be shown that a medium with a normal GVD introduces a time broadening to the propagating pulse, while the pulse broadening is less pronounced in a medium with an anomalous GVD, whereby optical soliton propagation is possible [33]. This could be understood by looking at the effect of dispersion on the newly generated frequencies due to the SPM, in the two different dispersion regimes. In the normal dispersion regime, the lower frequencies generated at the leading edge tend to travel faster than the higher frequency components at the trailing edge of the pulse, hence, SPM leads to enhanced rate of broadening of the pulse. However, in the anomalous dispersion regime, trailing edge of the pulse travels faster than the leading edge which results in the much smaller rate of the pulse broadening than expected in the absence of SPM. We have studied the effect of GVD on the self-focusing of THz pulses in the two materials with different GVDs. Figures 3 and 4 show the intensity distribution of the THz field in Si and ZnSe, respectively. The subplots show the pulses at the beginning, at the middle, and at the end of the sample (left to right) and under different polarization models of linear, Kerr and Kerr+Raman (top to bottom). It is shown in Fig. 3 that, during the nonlinear propagation (second and third rows), the initial spatial Gaussian profile of the pulse does not remain Gaussian in Si and breaks up into a strong self-focused trailing peak and a weak leading pulse. However, this is not the case with ZnSe material (See Fig. 4). In ZnSe, due to the 200 times stronger GVD parameter, the bending of spatio-temporal profile does not occur. In other words, normal dispersion compensates for the nonlinear refraction, inhibiting the self-focusing effect. The described interplay of SPM and GVD can be quantified by taking the ratio of the dispersion length to the nonlinear length: $L_{\textrm {D}}/L_{\textrm {NL}}$ [33]. This ratio governs the relative importance of GVD and SPM effects during the pulse propagation. If this ratio is much greater than one, SPM dominates, otherwise, GVD dominates. Here, we calculate $L_{\textrm {D}}/L_{\textrm {NL}}$ for both Si and ZnSe to be able to see the relative strength of GVD over SPM in the two samples. The dispersion length ($L_{\textrm {D}}$) can be defined as the length at which the effective pulse width increases by a certain factor, and is related to the GVD parameter as

Here, $\tau =0.5$ ps is the initial time duration of the pulse, and $\beta _2$ is the GVD parameter (in $\textrm {s}^2/\textrm {m}$). The nonlinear length ($L_{\textrm {NL}}$) is the effective propagation length at which the acquired nonlinear phase shift of the pulse is one radian. It depends on the peak power P of the pulse as where $\gamma$ is the nonlinear coefficient, defined as [33] Here $n_2$ is the nonlinear refractive index, c is the speed of light in vacuum and $A_{\textrm {eff}}$ is the effective mode area, which in our case is $A_{\textrm {eff}}=7.07\times 10^{-4}\,\,\textrm {m}^2$. The peak power of the pulse is $P=7.07\times 10^9$ W. We calculated $\gamma$ to be $3\times 10^{-8}\,\,\textrm {(Wm)}^{-1}$ at the central frequency 1 THz. Therefore, for Si we get the ratio of the dispersion to the nonlinear lengths to be $L^{\textrm {Si}}_{\textrm {D}}/L^{\textrm {Si}}_{\textrm {NL}}\approx 520$. The same calculation for ZnSe yields $L^{\textrm {ZnSe}}_{\textrm {D}}/L^{\textrm {ZnSe}}_{\textrm {NL}}\approx 1$. This explains why self-focusing does not occur in ZnSe (Fig. 4): due to the larger contribution from GVD.#### 3.3 Diffraction and spectral broadening

Due to the super-long wavelength regime of THz radiation, it is more diffractive than optical or near-infrared radiation. We take a closer look at the diffraction of the beam and see how linear and nonlinear behaviour of the material can modify its diffraction pattern. Figure 5 shows the diffraction of the spectral components of the THz beam in Si in arbitrary logarithmic units.

Higher frequency components naturally tend to localize closer to the beam axis than the lower frequency components (see the outward curvature of the spectral distribution). Spectral broadening due to SPM is also shown in the same figure. One can see that, except for the linear case (Figs. 5(a)–5(c)), there is a strong spectral broadening occurring as a result of the nonzero third-order susceptibility. Moreover, by comparing the second and third rows of Fig. 5, we note that although Kerr strength in the third row is 10$\%$ less than the second row, Raman has substituted the missing Kerr effect to reproduce the same amount of spectral broadening. Thus, Raman effect contributes equally to the frequency generation as the Kerr effect. This spectral broadening is more pronounced for quasi-longitudinal rays compared to the off-axis rays located far from the beam axis.

Figure 6 shows the spectral power of the transverse components of the THz beam in ZnSe in arbitrary logarithmic units. In this case, linear absorption strongly suppresses the higher frequency portion of the spectrum. As mentioned in the previous part, the efficiency of SPM is low due to the strong normal GVD of the ZnSe in this range. Spectral broadening still occurs, but it is less effective than that acquired by the beam in Si. For each plot in the same row, it is seen that the spectrum shrinks as the beam moves forward. Also, unlike the linear case, Kerr and Kerr+Raman propagations demonstrate some spectral broadening which is considerably less efficient than the broadening of the spectrum in Si. Moreover, spectral broadening is the same for the Kerr and Kerr+Raman effects, which implies the equal contribution of Raman effect and spontaneous Kerr effect to the spectral broadening.

## 4. Conclusion

Nonlinear THz spectroscopy and coherent control of the vibrational states of an arrangement of atoms or molecules in different crystalline structures is an evolving field waiting for more efficient and more accurate modelling techniques to better understand light-matter interaction at the THz frequencies. Also, finding some promising materials that show highly nonlinear behaviour in the THz frequency range seems to be an emerging quest. We presented a model of THz pulse propagation with the emphasis on some of the unique properties of low frequency THz beams, such as non-paraxiality, as well as the single-cycle regime of the pulse. We verified that there is a strong role played by the phonon contribution to the overall Kerr-type THz nonlinear optical interactions. This model helps with the extraction of $n_2$ values during the interpretation of nonlinear THz-TDS measurements. It also helps with finding highly nonlinear solids that show high $n_2$ levels due to the proximity of their phonon responses at the THz frequencies.

## Funding

Natural Sciences and Engineering Research Council of Canada; Canada Research Chairs.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **J. Hebling, K. Yeh, M. C. Hoffmann, B. Bartal, and K. A. Nelson, “Generation of high-power terahertz pulses by tilted-pulse-front excitation and their application possibilities,” J. Opt. Soc. Am. B **25**(7), B6–B19 (2008). [CrossRef]

**2. **C. Vicario, B. Monoszlai, and C. P. Hauri, “GV/m Single-Cycle Terahertz Fields from a Laser-Driven Large-Size Partitioned Organic Crystal,” Phys. Rev. Lett. **112**(21), 213901 (2014). [CrossRef]

**3. **A. Gopal, S. Herzer, A. Schmidt, P. Singh, A. Reinhard, W. Ziegler, D. Brömmel, A. Karmakar, P. Gibbon, U. Dillner, T. May, H-G. Meyer, and G. G. Paulus, “Observation of Gigawatt-Class THz Pulses from a Compact Laser-Driven Particle Accelerator,” Phys. Rev. Lett. **111**(7), 074802 (2013). [CrossRef]

**4. **K. P. Cheung and D. H. Auston, “Excitation of Coherent Phonon Polaritons with Femtosecond Optical Pulses,” Phys. Rev. Lett. **55**(20), 2152–2155 (1985). [CrossRef]

**5. **K. P. Cheung and D. H. Auston, “A novel technique for measuring far-infrared absorption and dispersion,” Infrared Phys. **26**(1), 23–27 (1986). [CrossRef]

**6. **R. R. Jones, D. You, and P. H. Bucksbaum, “Ionization of rydberg atoms by subpicosecond half-cycle electromagnetic pulses,” Phys. Rev. Lett. **70**(9), 1236–1239 (1993). [CrossRef]

**7. **D. Turchinovich, J. M. Hvam, and M. C. Hoffmann, “Self-phase modulation of a single-cycle terahertz pulse by nonlinear free-carrier response in a semiconductor,” Phys. Rev. B **85**(20), 201304 (2012). [CrossRef]

**8. **A. T. Tarekegne, K. Iwaszczuk, M. Zalkovskij, A. C. Strikwerda, and P. U. Jepsen, “Impact ionization in high resistivity silicon induced by an intense terahertz field enhanced by an antenna array,” New J. Phys. **17**(4), 043002 (2015). [CrossRef]

**9. **S. Tani, F. Blanchard, and K. Tanaka, “Ultrafast Carrier Dynamics in Graphene under a High Electric Field,” Phys. Rev. Lett. **109**(16), 166603 (2012). [CrossRef]

**10. **C. Lange, T. Maag, M. Hohenleutner, S. Baierl, O. Schubert, E. R. J. Edwards, D. Bougeard, G. Woltersdorf, and R. Huber, “Extremely nonperturbative nonlinearities in gaas driven by atomically strong terahertz fields in gold metamaterials,” Phys. Rev. Lett. **113**(22), 227401 (2014). [CrossRef]

**11. **P. Gaal, K. Reimann, M. Woerner, T. Elsaesser, R. Hey, and K. H. Ploog, “Nonlinear terahertz response of n-type gaas,” Phys. Rev. Lett. **96**(18), 187402 (2006). [CrossRef]

**12. **X. Chai, X. Ropagnol, S. M. Raeis-Zadeh, M. Reid, S. Safavi-Naeini, and T. Ozaki, “Subcycle terahertz nonlinear optics,” Phys. Rev. Lett. **121**(14), 143901 (2018). [CrossRef]

**13. **H. Y. Hwang, N. C. Brandt, H. Farhat, A. L. Hsu, J. Kong, and K. A. Nelson, “Nonlinear thz conductivity dynamics in p-type cvd-grown graphene,” J. Phys. Chem. B **117**(49), 15819–15824 (2013). [CrossRef]

**14. **L. Razzari, F. H. Su, G. Sharma, F. Blanchard, A. Ayesheshim, H.-C. Bandulet, R. Morandotti, J.-C. Kieffer, T. Ozaki, M. Reid, and F. A. Hegmann, “Nonlinear ultrafast modulation of the optical absorption of intense few-cycle terahertz pulses in n-doped semiconductors,” Phys. Rev. B **79**(19), 193204 (2009). [CrossRef]

**15. **H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green, N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Teichert, Z. Wang, K.-J. Tielrooij, Z. Liu, Z. Chen, A. Narita, K. Müllen, M. Bonn, M. Gensch, and D. Turchinovich, “Extremely efficient terahertz high-harmonic generation in graphene by hot dirac fermions,” Nature **561**(7724), 507–511 (2018). [CrossRef]

**16. **O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, “Sub-cycle control of terahertz high-harmonic generation by dynamical Bloch oscillations,” Nat. Photonics **8**(2), 119–123 (2014). [CrossRef]

**17. **J. Hebling, M. C. Hoffmann, K.-L. Yeh, G. Tóth, and K. A. Nelson, “Nonlinear lattice response observed through terahertz spm,” in * Ultrafast Phenomena XVI*, P. Corkum, S. Silvestri, K. A. Nelson, E. Riedle, and R. W. Schoenlein, eds. (Springer, Berlin, Heidelberg, 2009), pp. 651–653.

**18. **J. Hebling, K. Yeh, M. C. Hoffmann, and K. A. Nelson, “High-Power THz Generation, THz Nonlinear Optics, and THz Nonlinear Spectroscopy,” IEEE J. Sel. Top. Quantum Electron. **14**(2), 345–353 (2008). [CrossRef]

**19. **K. Dolgaleva, D. V. Materikina, R. W. Boyd, and S. A. Kozlov, “Prediction of an extremely large nonlinear refractive index for crystals at terahertz frequencies,” Phys. Rev. A **92**(2), 023809 (2015). [CrossRef]

**20. **K. J. Blow and D. Wood, “Theoretical description of transient stimulated raman scattering in optical fibers,” IEEE J. Quantum Electron. **25**(12), 2665–2673 (1989). [CrossRef]

**21. **R. W. Boyd, * Nonlinear Optics* (Elsevier, 2003).

**22. **T. Brabec and F. Krausz, “Nonlinear optical pulse propagation in the single-cycle regime,” Phys. Rev. Lett. **78**(17), 3282–3285 (1997). [CrossRef]

**23. **M. Kolesik and J. V. Moloney, “Nonlinear optical pulse propagation simulation: From maxwell’s to unidirectional equations,” Phys. Rev. E **70**(3), 036604 (2004). [CrossRef]

**24. **J. Andreasen and M. Kolesik, “Nonlinear propagation of light in structured media: Generalized unidirectional pulse propagation equations,” Phys. Rev. E **86**(3), 036706 (2012). [CrossRef]

**25. **A. Couairon, E. Brambilla, T. Corti, D. Majus, O. de J. Ramírez-Góngora, and M. Kolesik, “Practitioner’s guide to laser pulse propagation models and simulation,” Eur. Phys. J.: Spec. Top. **199**(1), 5–76 (2011). [CrossRef]

**26. **M. Kolesik, J. V. Moloney, and M. Mlejnek, “Unidirectional optical pulse propagation equation,” Phys. Rev. Lett. **89**(28), 283902 (2002). [CrossRef]

**27. **M. D. Feit and J. A. Fleck, “Beam nonparaxiality, filament formation, and beam breakup in the self-focusing of optical beams,” J. Opt. Soc. Am. B **5**(3), 633–640 (1988). [CrossRef]

**28. **L. M. Kovachev, “The light filament as vector solitary wave,” AIP Conf. Proc. **1629**(1), 167–171 (2014). [CrossRef]

**29. **P. Saari, “Evolution of subcycle pulses in nonparaxial gaussian beams,” Opt. Express **8**(11), 590–598 (2001). [CrossRef]

**30. **M. Kolesik, “gUPPEcore simulation framework” (2019). http://acms.arizona.edu/FemtoTheory/MK_personal/guppelab/

**31. **J. Dai, J. Zhang, W. Zhang, and D. Grischkowsky, “Terahertz time-domain spectroscopy characterization of the far-infrared absorption and index of refraction of high-resistivity, float-zone silicon,” J. Opt. Soc. Am. B **21**(7), 1379–1386 (2004). [CrossRef]

**32. **A. Deneuville, D. Tanner, and P. H. Holloway, “Optical constants of znse in the far infrared,” Phys. Rev. B **43**(8), 6544–6550 (1991). [CrossRef]

**33. **G. P. Agrawal, * Nonlinear Fiber Optics* (Academic, 2001).