Open Access
Add to BibSonomy
Abstract
We establish a first-principle model for the simulation of spatiotemporal light pulse dynamics based on the combination of the time-dependent Schrödinger equation and the unidirectional propagation equation. The proposed numerical scheme enables computationally efficient simulation while being stable and accurate. We use the new model to examine self-focusing of a short pulse in atomic hydrogen and show that an accurate description of the excited-levels dynamics can only be achieved by a propagation model with an ab-initio description of the light-matter interaction, which accounts for the laser-dressed multilevel structure of the system, including bound and free states, and its sub-cycle response.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Since their discovery, self-focusing [1] and filamentation [2] of light in gases have attracted considerable attention due to the fascinating effects associated with these phenomena. Femtosecond pulse filamentation relies on the balance between self-focusing caused by the "nonlinear lens" formed by the Kerr effect and de-focusing caused by the photo-induced electron plasma [3]. While initially the filamentation mechanism itself was primarily studied in various regimes, attention has recently shifted to potential applications in areas such as THz generation [4], LIDAR atmospheric studies [5], detection and sensing [6], air lasing [7], light bullets [8] and pulse self-compression [9].
To adequately account for the effects relevant to such applications, a model is needed that describes the propagation of an intense, short light pulse in the medium and accurately predicts the spatiotemporal propagation dynamics by considering all important optical processes such as self-focusing, plasma formation and its effects, diffraction and group velocity dispersion. To date, a significant number of such models has been developed and successfully applied in a broad range of studies (for reviews see e.g. [10,11]). Most of them are based on the phenomenological description of the local response of the gas constituents, typically in terms of (linear) group velocity dispersion, Kerr-type third-order optical nonlinearity and plasma contribution augmented by an explicit expression for the ionization rate, such as the Ammosov-Delone-Krainov formula. However, the phenomenological models lack an ab-initio description of the (light-dressed) level structure of the atoms or molecules and thus some important effects, including population dynamics of the excited states, three-step high harmonic generation, Freeman resonances, etc.
Effects where the microscopic structure of the system is important for an accurate description of light propagation and related applications can play an essential role in many situations, e.g. for the description of optical nonlinearities [12] or light amplification enabled by special light-dressed states that are stable against ionization [13]. In such cases, the time-dependent Schrödinger equation (TDSE) provides the most accurate first-principle description of the local response of an atom or molecule. However, due to the formidable computational effort required for an accurate TDSE solution, there are few studies that combine pulse propagation equations with TDSE calculations. In Ref. [14,15], the TDSE was used together with a finite-difference time-domain method for light propagation. Several works have proposed refinements to this technique, e.g. finite-element schemes [16]. However, this approach forfeits the chances of improved numerical efficiency associated with the unidirectional propagation of light in homogeneous or near-homogeneous media, and propagation can only be considered in the range of a few tens of micrometers. In Ref. [17], the propagation of the pump light and that of the light generated by molecular polarization are decoupled, which is a numerically efficient scheme but ignores the feedback of the generated light to the pump and does not allow simulation of stimulated emission processes. Feedback processes were considered in Ref. [18], but in a (1+1)D model that considers one spatial coordinate and time and therefore neglects transverse dynamics and self-focusing. Another approach is based on metastable electronic states [19] and is most suitable for non-resonant systems. A comprehensive, stable, and numerically affordable theoretical model in (3+1)D or (2+1)D (three or two spatial coordinates and time) suitable for accurate simulation of spatiotemporal pulse propagation dynamics has been lacking.
In this paper, we establish such an efficient coupled TDSE-Maxwell approach and elucidate its key and novel ingredients, such as a linear approximation of the forward-propagating nonlinear polarization that provides a numerically efficient spatial step size for unidirectional pulse propagation. We demonstrate its performance using the example of spatiotemporal dynamics in a model system of atomic hydrogen gas. Notably, the applicability of our novel approach goes well beyond self-focusing. It can be applied to any problem where the accurate description of unidirectional, non-relativistic propagation requires a first-principle description of light-matter interaction. It can also be extended to liquids and solids, provided that the TDSE is modified accordingly; in this case, the novel ingredients mentioned above would be particularly useful, since the numerical effort is typically greater in optically dense media.
The paper is organized as follows. In the following Section, the numerical scheme is designed and derived. In Section 3, propagation results for different pulse intensities and models of the local response of the hydrogen atom are discussed, followed by our conclusions.
2. Simulation model
We consider a situation where the change of the (nonlinear) refractive index along the propagation direction $z$ is small and/or slow. In this case, the back-reflection can be neglected and a unidirectional equation can be used. In the case of an optically dilute homogeneous medium such as the gas considered here, it takes a particularly simple form [20],
The polarization $P(z,r_\perp,t)$ contains a nonlinear part $P_{\rm NL}(z,r_\perp,t)$ and a linear part $P_{\rm L}(z,r_\perp,t)$, the latter being related to the electric field by a linear susceptibility $\chi (\omega )$. We have previously shown [21] that for a resonant medium, a propagator operator for the field $E$ provides the most efficient and artifact-free approach for pulse propagation. It consists of subtracting $P_{\rm L}(z,r_\perp,t)$ from the total polarization, so that only $P_{\rm NL}(z,r_\perp,t)$ is substituted into Eq. (1). The linear susceptibility is accounted for separately, using the time-domain propagation operator $\hat {B}$ which connects electric field at positions $z$ and $z+dz$ as $E(z+dz,r_\perp,t)=\hat {B}E(z,r_\perp,t)$. This operator is defined by
The diffraction term in Eq. (1) was implemented in the frequency domain by a standard Crank–Nicolson method, augmented by first-order absorbing boundary conditions. We assume azimuthal symmetry and, as usual for a (2+1)D approach, keep only the radial coordinate, which is a good approximation in many practical situations and reduces the required numerical effort.
To calculate the polarization as a function of the electric field, we have used a 3D and a 1D TDSE solver (not to be confused with the global (2+1) dimensionality of the model defined by the spatial coordinates z, $r_\perp$ and the time $t$, which describe the macroscopic pulse propagation) to calculate the local response of the hydrogen atom and a model hydrogen atom, respectively. The time-dependent dipole moment obtained as the output of the TDSE solution determines the nonlinear polarization $P_{\rm NL}(z,r_\perp,t)$, which is substituted into Eq. (1), providing the coupling between the TDSE and propagation equation. The 1D simulations require less numerical effort and produce semi-quantitative results, while the 3D approach provides an accurate description of the electron dynamics in hydrogen. The hydrogen atom is initially in the ground state. The TDSE was solved using the code described in [22,23]. Convergence has been checked against all discretization parameters. A uniform radial spacing of $0.1$ a.u.$=5.29\times 10^{-12}$ m and $l_{\textrm {max}}=40$ for the angular part were used. For the time discretization, a time step of $dt = 0.01$ a.u.$=2.42\times 10^{-19}$ s was used. For the 1D hydrogen model, we used the following soft-core potential [24]:
where the parameter $a=1.4142$ is chosen to reproduce the ground state of the hydrogen atom with ionization potential $I_p=13.6$ eV.Equation (1) was solved using a second-order predictor-corrector method in the so-called predict-evaluate-correct (PEC) mode [25], which allows us to use only one TDSE call per propagation step,
We discretize the $r_\perp$ coordinate using 128 points, with a $dr_\perp$ of 2 $\mu$m. For each step in z, we run the necessary 128 TDSE calculations in parallel.
We have found that the above model, while rigorously derived, suffers from numerical instability for any $dz$ that would allow efficient numerical calculations. This instability typically manifests itself in the accumulation of radiation at frequencies corresponding to transitions with large dipole moments, followed by a catastrophic growth of the electric field. We have examined different approaches to correct this numerical issue.
One possible solution is to apply at each $z$ step a smooth spectral filter, which completely removes all angular frequencies above 9 fs$^{-1}$, corresponding to 6 eV. Such low-pass filtering has been used previously [20] to achieve numerical stability. The spectral range above 9 fs$^{-1}$ includes the transition frequency of 10.9 fs$^{-1}$ corresponding to the transition between the two lowest states of 1D hydrogen atom, which is characterized by the largest dipole moment of the system. We denote this model with "F" as in "Filtering".
Another approach is to develop a more accurate formalism for the extraction of the forward contribution from the available polarization data. It is related to the fact that the nonlinear part of the polarization contains some small contribution to the backward-propagating wave, which can be projected out as detailed below. For an infinitesimally small propagation step, $dz\ll 2\pi c/\omega$, one could use the polarization at each step directly. However, this approach would be extremely numerically inefficient. To keep the size of the data memory manageable, at each step $z$, we use only $P_{\rm NL}(z,r_\perp,\omega )$ and the polarization of the previous step, $P_{\rm NL}(z-dz,r_\perp,\omega )$, to linearly approximate the polarization at position $z^*$ between the grid points in $z$ direction:
In the co-moving frame, the forward-propagating wave $P_f$ has no $z$ dependence, $P_f\sim 1$, and the backward-propagating wave $P_b$ depends on z as $P_b\sim {\rm exp}[-2i\omega z^*/c]$. We represent the linear approximation given by Eq. (5) in the (generally speaking, non-orthogonal) basis $\{P_f,P_b\}$ to obtain the forward component of the polarization. After a straightforward calculation, we derive the following expression for the forward nonlinear polarization $P_{\rm NL}^{\rm forw}$ to be substituted into Eq. (1) at point $z$:
Fig. 1. Function $K(k(\omega )dz)$ defined in Eq. (7). Note the change of the vertical axis from linear to logarithmic at $K(k(\omega )dz)=0.035$.
The comparison of the results obtained with the above corrections is presented in Fig. 2. In particular for ionization and linear losses the results of the two cases differ significantly. This is to be expected because the suppression of radiation at the frequency of the strongest transition (F model) reduces the ionization rate and prevents linear losses at that frequency. Since the F model has no physical justification and is energy-level-specific, we evaluate it as less suitable and adopt the B model approach in the following.
Fig. 2. Changes of electric field on axis due to ionization (red curves), diffraction (green curves), and linear absorption (blue curves) calculated by the F model (thin dashed curves) and B model (thick solid curves), using an 8 fs, 100 TW/cm$^2$ pulse at 830 nm and with a beam waist of 30 $\mu$m. Change of electric field due to diffraction almost coincides in F and B models.
The correct choice of the step $dz$ ensures a balance between the accuracy of the simulation and a reasonable numerical effort. We have analyzed the dynamics of the system for different $dz$, as shown in Fig. 3 (inset) for parameters given in the caption. We find that the difference of the results for $dz=0.25$ $\mu$m (red curve) and $dz=0.5$ $\mu$m (green curve) is negligible, while a step size of $dz=1$ $\mu$m (blue curve) leads to a large change of loss compared to $dz=0.5$ $\mu$m, resulting in lower accuracy. Thus we consider $dz=0.5$ $\mu$m as the optimum step at an intensity of 100 TW/cm$^2$. The dependence of the optimum step on intensity is presented in Fig. 3. In the code, an adaptive, intensity-dependent (and therefore $z$-dependent) choice of the propagation step $dz$ was implemented.
Fig. 3. Dependence of the optimum step in $z$ on the propagation intensity. In the inset, the dependence of the changes of electric field on axis due to ionization is shown for $dz=0.25$ $\mu$m (red curve), $dz=0.5$ $\mu$m (green curve), and $dz=1$ $\mu$m (blue curve). Propagation of an 8 fs pulses at 830 nm with a beam width of 30 $\mu$m was considered.
To gauge the impact of the TDSE-generated atomic response on the spatiotemporal dynamics of the pulse, we compare our TDSE-based results to those obtained with a much simpler response model based on the strong-field approximation (SFA) for the electron dynamics and photoionization. In this case, the above described numerical scheme replaces the TDSE calculations by the following model for the nonlinear polarization [26]:
3. Numerical results and discussion
Let us first analyze the results of the SFA-based model shown in Fig. 4, with parameters given in the caption. For an intensity of 15 TW/cm$^2$, the input pulse has a peak power of 2.1 GW, well above the critical power for self-focusing ($\sim$ 0.65 GW), resulting in rapid self-focusing, see panel 4(a). Here and hereafter, hydrogen at a pressure of 1 atm is considered. Due to photoionization and associated absorption, no high-order filament is formed and plasma defocusing overtakes self-focusing after 15 mm of propagation. In the spectrum, see panel 4(c), we find that the time-dependent plasma contribution to the refractive index leads to nonlinear phase accumulation and a corresponding blue shift as well as the generation of the third harmonic. This model predicts no additional distinct peaks at other frequencies of the spectrum, in contrast to the 1D-TDSE-based model presented below.
Fig. 4. Propagation results obtained with the SFA-based model. (a,b) Pulse energy distribution during propagation (a,b) and output spectrum (c,d) for an 8 fs input pulse at 830 nm with a waist of 67 $\mu$m and an intensity of 15 TW/cm$^2$ (a,c) and 250 TW/cm$^2$ (b,d). The input pulse has Gaussian profile both in time and $r_\perp$. We consider propagation in atomic hydrogen at standard conditions, using parameters fitted to the 1D TDSE results. Shown are the area-integrated spectra at the end of the propagation.
For a much higher input intensity of 250 TW/cm$^2$, plasma creation by photoionization leads to significant absorption in the beam center, accompanied by defocusing, see panel 4(b). Hence, after a few millimeters of propagation, the beam acquires an almost top-hat profile in the transverse direction. In the corresponding spectrum in panel 4(d), one finds significant spectral broadening around the fundamental frequency due to the plasma contribution, as well as generation of the third harmonic, which is structured due to the structured spectrum of the fundamental.
Figure 5 presents the results obtained with the 1D-TDSE-based model. For the lower intensity case of 15 TW/cm$^2$, we find that after the initial drop in intensity caused by photoionization losses, self-focusing starts to overcome transverse spreading, see panel 5(a). Indeed, panel 5(c) shows that at about 1.8 mm propagation distance, the intensity change on the beam axis becomes positive, since the increase in intensity due to the diffraction term (green curve) is larger than the sum of the nonlinear (red) and linear (blue) losses. The dependence of the photoionization loss (and hence peak intensity) on the propagation coordinate (red curve) appears to be the initial stage of a periodic behavior typical for filamentation in gases and related to spatial replenishment [29].
Fig. 5. Propagation results obtained with the 1D TDSE-based model. (a),(b) Pulse energy distribution during propagation, (c),(d) on-axis energy evolution due to photoionization absorption (red curve), diffraction (green curve) and linear loss (blue curve), and (e),(f) area-integrated spectrum at the end of propagation for an 8 fs input pulse at 830 nm. In (a),(c),(e) the input beam has a peak intensity of 15 TW/cm$^2$ and a waist of 67 $\mu$m, in (b),(d),(f) it has a peak intensity of 250 TW/cm$^2$ and a waist of 30 $\mu$m. Shown are the area-integrated spectra at the end of the propagation.
The corresponding spectrum, shown in panel 5(e), exhibits spectral broadening due to self-phase modulation and generation of the third harmonic. Importantly, the spectrum exhibits another prominent feature - a peak at 4.85$\hbar \omega _0=7.2$ eV, corresponding to the transition energy between the ground state and the first excited state, both coherently populated by the pulse. This transition contributes to the spectrum due to the strong contribution to the polarization oscillating at this frequency. Moreover, lying close to the 5$^{th}$ harmonic of the pump light, this coherence is amplified. We have checked (not shown here) that by choosing a central angular driving frequency of $\omega _0=2.72$ fs$^{-1}$, corresponding to 1.8 eV, this transition gets much less excited as it corresponds to the 4$^{th}$ harmonic, which is not generated by the perturbative nonlinearities due to the point inversion symmetry of the medium.
Panels 5(b),(d),(f) show the 1D-TDSE-based results for the higher input intensity of 250 TW/cm$^2$. The evolution of the pulse energy distribution, see panel 5(b), differs significantly from that predicted by the SFA-model, cf. panel 4(b), developing a ring-like structure after about 50 $\mu$m of propagation. The notable defocusing of the beam is a result of the strong energy absorption due to photoionization, which dominates the propagation and leads to a significant reduction of the pulse energy over a fraction of a millimeter, see panel 5(d). The spectrum in panel 5(f) shows again strong spectral broadening, third harmonic generation, and a peak around 4.85$\omega _0$.
The novel capabilities of the TDSE-based propagation model rely on the TDSE-generated polarization being accounted for at each step of the pulse propagation. To demonstrate this, we have performed simulations using three different models for the microscopic response of the hydrogen atom: (i) SFA-based propagation, (ii) TDSE-based propagation, and (iii) a model that we denote as TDSE$^{1/4}$. In the latter model, the SFA-based polarization is calculated at each step and used for propagation. In addition, the TDSE-based polarization is calculated at every fourth step. After each TDSE call, the fourfold difference between the TDSE-based polarization and the SFA-based polarization is substituted as an additional polarization term into Eq. (1). The TDSE$^{1/4}$ model is motivated by the reduction in numerical effort that results from performing the time-consuming TDSE calculations only at every fourth step in $z$.
The results obtained with the three different models are shown in Fig. 6. Panel 6(a) shows that the TDSE-based model and the TDSE$^{1/4}$ model yield indistinguishable transverse profiles of the propagating beam. However, the SFA-model results deviate notably and the deviation is expected to increase with further propagation. We attribute the stronger photoionization losses in the TDSE-based model to the role of the radiation generated around 4.85$\omega _0$. Although the SFA photoionization rate was corrected to match the TDSE-based ionization rate for the local response, this correction does not account for all the spectral components of the light newly generated during propagation and predicted only by the TDSE-based model. Indeed, while the transverse intensity profiles corresponding to the TDSE-based model and the TDSE$^{1/4}$ model virtually coincide, the electronic state populations are quite different, see panel 6(b). The populations agree only up to a factor of about 3, indicating that the use of the TDSE-generated polarization at each propagation step is critical to obtain accurate population distributions of the system, which, in turn, affect the subsequent pulse propagation. In general, to capture these effects, it is not sufficient to simulate pulse propagation governed by the SFA and use TDSE calculations only at certain propagation coordinates of interest, as it was done previously [12,30]. We note that the differences in the population distributions do not manifest in the results in panel 6(a) as the transverse pulse energy distributions are integrated over the spectrum.
Fig. 6. (a) Transverse pulse energy distribution and (b) populations of excited states after propagation over 250 $\mu$m obtained with the TDSE-based model (red curves), the TDSE$^{1/4}$ model (TDSE is applied only every 4th step, green curves), and the SFA-based model (blue curve). In (a), the initial energy distribution is shown by the grey curve. An 8 fs pulse with a central wavelength of 830 nm, an intensity of 100 TW/cm$^2$ and a waist of 30 $\mu$m was used. Populations of excited states are averaged over spatial coordinates, using pulse fluence at each point in space as weight.
The results obtained for the 1D TDSE hydrogen model provide useful insights into the features of the newly developed numerical scheme, but the 1D model does not describe the realistic energy level structure of the hydrogen atom. In Fig. 7 we present the dependence of the energy level populations on the propagation distance obtained with the 3D-TDSE-based model, which shows periodic features. The population oscillations can be attributed to the drop of intensity and corresponding change of pulse area analog to the periodic dependence of the level populations on the pulse area well known for a (much simpler) two-level system. In principle, the change in pulse area could also be predicted by an SFA-based model and certain population oscillations may be inferred. However, with further propagation, some features appear that cannot be predicted by SFA models, such as the sudden (transient) increase in the population of the high-lying $6f$ state, starting around 75 $\mu$m, which can only be revealed by a first-principles TDSE-based model that takes into account the electronic structure of the (light-dressed) system and the complex light-atom interactions. The behavior of the population of the $6f$ state cannot be described in the framework of the pulse-area concept, since the population is clearly non-periodic in the propagation coordinate.
Fig. 7. Dependence of the on-axis energy-level populations on the propagation distance obtained with the 3D-TDSE-based model, for the levels 2s (red solid curve), 2p (green dashed curve), 3p (blue short-dashed curve), and 6f (magenta dotted curve). An 8 fs pulse with a central wavelength of 830 nm, an intensity of 250 TW/cm$^2$ and a waist of 30 $\mu$m was used.
4. Conclusion
We have shown that a combination of the time-dependent Schrödinger equation and unidirectional propagation equation is a suitable basis for an accurate first-principle model for simulation of spatiotemporal light pulse dynamics. Use of the propagator operator and correction of the polarization to exclude the contribution to the backward-propagating wave are demonstrated to be key ingredients for the numerically efficient and artifact-free simulation of relevant effects such as self-focusing and plasma defocusing. We have shown that the TDSE-based polarization needs to be taken into account at each propagation step in order to accurately address the population distribution across the energy levels as well as the spectrum of the propagating field.
Funding
Deutsche Forschungsgemeinschaft (HU 1593/11-1).
Acknowledgments
We gratefully acknowledge stimulating and fruitful discussions with Misha Ivanov. A. H. acknowledges financial support from German Research Council, grant DFG HU 1593/11-1.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. M. R. Rashidian Vaziri, “Describing the propagation of intense laser pulses in nonlinear Kerr media using the ducting model,” Laser Phys. 23(10), 105401 (2013). [CrossRef]
2. M. Hercher, “Laser-induced damage in transparent media,” J. Opt. Soc. Am. 54, 563 (1964). [CrossRef]
3. N. Aközbek, C. M. Bowden, A. Talebpour, and S. L. Chin, “Femtosecond pulse propagation in air: Variational analysis,” Phys. Rev. E 61(4), 4540–4549 (2000). [CrossRef]
4. H. G. Roskos, M. D. Thomson, M. Kreß, and T. Löffler, “Broadband THz emission from gas plasmas induced by femtosecond optical pulses: from fundamentals to applications,” Laser Photonics Rev. 1(4), 349–368 (2010). [CrossRef]
5. P. Rairoux, H. Schillinger, S. Niedermeier, M. Rodriguez, F. Ronneberger, R. Sauerbrey, B. Stein, D. Waite, C. Wedekind, H. Wille, L. Wöste, and C. Ziener, “Remote sensing of the atmosphere using ultrashort laser pulses,” Appl. Phys. B 71(4), 573–580 (2000). [CrossRef]
6. J. L. Liu, J. M. Dai, and X. C. Zhang, “Broadband terahertz wave remote sensing using coherent manipulation of fluorescence from asymmetrically ionized gases,” Nat. Photonics 4(9), 627–631 (2010). [CrossRef]
7. A. Dogariu, J. B. Michael, M. O. Scully, and R. B. Miles, “High-gain backward lasing in air,” Science 331(6016), 442–445 (2011). [CrossRef]
8. Y. Silberberg, “Collapse of optical pulses,” Opt. Lett. 15(22), 1282–1284 (1990). [CrossRef]
9. S. Skupin, G. Stibenz, L. Berge, F. Lederer, T. Sokollik, M. Schnürer, N. Zhavoronkov, and G. Steinmeyer, “Self-compression by femtosecond pulse filamentation: experiments versus numerical simulations,” Phys. Rev. E 74(5), 056604 (2006). [CrossRef]
10. A. Couairon and A. Myzyrowicz, “Femtosecond filamentation in transparent media,” Phys. Rep. 441(2-4), 47–189 (2007). [CrossRef]
11. R. W. Boyd, S. G. Lukishova, and Y. R. Shen eds., Self-Focusing: Past and Present. Fundamentals and Prospects (Springer, 2009).
12. M. Richter, S. Patchkovskii, F. Morales, O. Smirnova, and M. Ivanov, “The role of the Kramers–Henneberger atom in the higher-order Kerr effect,” New J. Phys. 15(8), 083012 (2013). [CrossRef]
13. M. Matthews, F. Morales, A. Patas, A. Lindinger, J. Gateau, N. Berti, S. Hermelin, J. Kasparian, M. Richter, T. Bredtmann, O. Smirnova, J.-P. Wolf, and M. Ivanov, “Amplification of intense light fields by nearly free electrons,” Nat. Phys. 14(7), 695–700 (2018). [CrossRef]
14. E. Lorin, S. Chelkowski, and A. Bandrauk, “A numerical Maxwell–Schrödinger model for intense laser–matter interaction and propagation,” Comput. Phys. Commun. 177(12), 908–932 (2007). [CrossRef]
15. I. P. Christov, “Enhanced generation of attosecond pulses in dispersion-controlled hollow-core fiber,” Phys. Rev. A 60(4), 3244–3250 (1999). [CrossRef]
16. C. H. Yao, Z. Y. Wang, and Y. M. Zhao, “A leap-frog finite element method for wave propagation of Maxwell–Schrödinger equations with nonlocal effect in metamaterials,” Computers Mathematics with Applications 90, 25–37 (2021). [CrossRef]
17. M. B Gaarde, J. L Tate, and K. J Schafer, “Macroscopic aspects of attosecond pulse generation,” J. Phys. B: At. Mol. Opt. Phys. 41(13), 132001 (2008). [CrossRef]
18. N. Berti, P. Bejot, E. Cormier, J. Kasparian, O. Faucher, and J.-P. Wolf, “Ab initio calculations of laser-atom interactions revealing harmonics feedback during macroscopic propagation,” Phys. Rev. A 99(6), 061401 (2019). [CrossRef]
19. M. Kolesik, J. M. Brown, A. Teleki, P. Jakobsen, J. V. Moloney, and E. M. Wright, “Metastable electronic states and nonlinear response for high-intensity optical pulses,” Optica 1(5), 323 (2014). [CrossRef]
20. A. Husakou and J. Herrmann, “Supercontinuum Generation of Higher-Order Solitons by Fission in Photonic Crystal Fibers,” Phys. Rev. Lett. 87(20), 203901 (2001). [CrossRef]
21. F. Morales, M. Richter, V. Olvo, and A. Husakou, “Propagator operator for pulse propagation in resonant media,” Opt. Express 29(18), 29128 (2021). [CrossRef]
22. S. Patchkovskii and H. G. Muller, “Simple, accurate, and efficient implementation of 1-electron atomic time-dependent Schrödinger equation in spherical coordinates,” Comput. Phys. Commun. 199, 153–169 (2016). [CrossRef]
23. H. G. Muller, “An efficient propagation scheme for the time-dependent Schrödinger equation in the velocity gauge,” Laser Phys. 9, 138–148 (1999).
24. J. Javanainen, J. H. Eberly, and Q. Su, “Numerical simulations of multiphoton ionization and above-threshold electron spectra,” Phys. Rev. A 38(7), 3430–3446 (1988). [CrossRef]
25. J. C. Butcher, Numerical Methods for Ordinary Differential Equations (John Wiley and Sons, New York, 2003).
26. A. Husakou and J. Herrmann, “High-power, high-coherence supercontinuum generation in dielectric-coated metallic hollow waveguides,” Opt. Express 17(15), 12481 (2009). [CrossRef]
27. G. L. Yudin and M. Y. Ivanov, “Nonadiabatic tunnel ionization: Looking inside a laser cycle,” Phys. Rev. A 64(1), 013409 (2001). [CrossRef]
28. F. Morales, M. Richter, M. Ivanov, and A. Husakou, “Non-instantaneous third-order optical response of gases in low-frequency fields,” Opt. Express 30(13), 23579 (2022). [CrossRef]
29. M. Mlejnek, E. M. Wright, and J. V. Moloney, “Power dependence of dynamic spatial replenishment of femtosecond pulses propagating in air,” Opt. Express 4(7), 223 (1999). [CrossRef]
30. B. Clough, N. Karpowicz, and X.-C. Zhang, “Modulation of electron trajectories inside a filament for single-scan coherent terahertz wave detection,” Appl. Phys. Lett. 100(12), 121105 (2012). [CrossRef]
