Adiabatic Floquet model for the optical response in femtosecond filaments

We briefly describe the methods for simulating the optical response of a one-dimensional hydrogen atom under the influence of a strong laser field. One-dimensional model atoms are popular toy models and give good qualitative predictions of the dynamics of atoms in strong laser fields. In particular, they were successfully used for theoretical studies of high harmonic generation and atomic stabilization in strong fields [24–26]. Starting from first principles, we write down the TDSE in the dipole approximation using atomic units (au)

$$\begin{eqnarray}&&{\rm{i}}{\partial }_{t}\psi (x,t)=\left[-\displaystyle \frac{1}{2}{\partial }_{x}^{2}-\displaystyle \frac{1}{\sqrt{{x}^{2}+{\alpha }^{2}}}+E(t)x\right]\psi (x,t).\end{eqnarray} \tag{ 2 }$$

To circumvent the singularity of the atomic potential at the origin, a soft-core potential is used which yields the correct ionization potential of hydrogen, ${I}_{{\rm{p}}}=13.6\,\mathrm{eV}\equiv 0.5\,\mathrm{au}$ if $\alpha =\sqrt{2}$. For integrating the TDSE numerically we discretize time $t\to {t}_{n}$ and space $x\to {x}_{j}$ with equal spacings $\tau ={\rm{\Delta }}t=0.01\,\mathrm{au}$ and $h={\rm{\Delta }}x=0.1\,\mathrm{au}$, respectively. For the second derivative in the Hamiltonian we can write the central difference

$$\begin{eqnarray}&&{\partial }_{x}^{2}\psi \approx \displaystyle \frac{{\psi }_{j-1}-2{\psi }_{j}+{\psi }_{j+1}}{{h}^{2}},\end{eqnarray} \tag{ 3 }$$

where ${\psi }_{j}^{n}$ is the wave function at time step t_n on the grid point x_j. Thus, the matrix representing H on the grid has a tridiagonal structure. The implicit Crank−Nicolson propagator is commonly chosen to compute the wave function for the next time step

$$\begin{eqnarray}U({t}_{n+1},{t}_{n}) & = & {\left(1+{\rm{i}}\displaystyle \frac{\tau }{2}H({t}_{n+\tfrac{1}{2}})\right)}^{-1}\left(1-{\rm{i}}\displaystyle \frac{\tau }{2}H({t}_{n+\tfrac{1}{2}})\right),\end{eqnarray} \tag{ 4 }$$

which preserves unitarity and is second order accurate in space and time [27]. When the propagator is applied to the wave function, i.e., ${\psi }^{n+1}=U({t}_{n+1},{t}_{n}){\psi }^{n}$, we obtain a set of linear equations

$$\begin{eqnarray}&&{M}_{n+\tfrac{1}{2}}^{+}{\psi }^{n+1}={M}_{n+\tfrac{1}{2}}^{-}{\psi }^{n}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\rm{with}}\quad {M}_{n}^{\pm }\equiv \left(1\pm {\rm{i}}\displaystyle \frac{\tau }{2}H({t}_{n})\right).\end{eqnarray} \tag{ 6 }$$

These equations can be efficiently solved using the tridiagonal matrix algorithm [28]. As an initial condition we provide that ${\psi }^{0}$ is the atom's ground state. Reflections at the spatial ends of the integration domain are suppressed by a 40 au wide absorbing boundary layer [29].

We consider a laser pulse with a central wavelength of 800 nm, which corresponds to an angular frequency of ${\omega }_{0}=0.057\,\mathrm{au}$ and an optical cycle length of $T=2.67\,\mathrm{fs}\equiv 110.3\,\mathrm{au}$. To ensure that the electric field $E(t)$ has a vanishing DC component, we introduce the vector potential $A(t)$, from which we reconstruct the electric field according to $E(t)=-{\partial }_{t}A(t)$. The vector potential has a cosine square envelope

$$\begin{eqnarray}&&A(t)={A}_{0}{\cos }^{2}(\pi t/{T}_{{\rm{p}}})\cos ({\omega }_{0}t)\end{eqnarray} \tag{ 7 }$$

for $| t| \leqslant {T}_{{\rm{p}}}/2$, where T_p is the total pulse duration, and the peak amplitudes of E and A are related via ${E}_{0}={\omega }_{0}{A}_{0}$. We already know from our previous work that Freeman resonances between laser-dressed states lead to drops in the refractive index, and that those resonances are more pronounced when a long flat-top pulse is used. We therefore additionally employ a flat-top pulse with a four-cycle ramp up and down, respectively, according to (7) but with an 80-cycle long central region of constant amplitude in between.

Having solved the TDSE for a single electron, the polarization density for a gas of hydrogen atoms at standard condition, with number density $\sigma =4\times {10}^{-6}\,\mathrm{au}$, can be obtained from the atomic dipole moment

$$\begin{eqnarray}&&{P}({t}_{n})=-\sigma \langle {\psi }^{n}| x| {\psi }^{n}\rangle .\end{eqnarray} \tag{ 8 }$$

Denoting $\widehat{X}(\omega )$ as the Fourier transform of $X(t)$, we can define a time-averaged susceptibility and thus a refractive index of the hydrogen gas interacting with the laser pulse

$$\begin{eqnarray}&&\chi =\displaystyle \frac{\widehat{P}({\omega }_{0})}{{\epsilon }_{0}\widehat{E}({\omega }_{0})},\quad n=\mathrm{Re}(\sqrt{1+\chi }).\end{eqnarray} \tag{ 9 }$$

By repeating the simulation with variable peak intensity I up to $40\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$, we obtain an intensity-dependent, nonlinear refractive index

$$\begin{eqnarray}&&{\rm{\Delta }}n(I)=n(I)-\mathop{\mathrm{lim}}\limits_{I\to 0}n(I).\end{eqnarray} \tag{ 10 }$$

To ensure that $n(I)$ correctly captures the response from ionized electrons, the simulation volume is constantly enlarged with increasing I. Thus, less than 0.1% of the total ionization is absorbed at the boundaries.

Our next objective is to seek for correlations between the nonlinear refractive index and the probability of ionization, the population of excited states and the ground state depletion. Thus, we compute the probability that the electron occupies a specific field-free bound state after the pulse is gone. This can be extracted directly from the wave function at the end of the pulse by projecting onto those eigenvectors of the discrete Hamiltonian which correspond to bound atomic states, denoted as ${\phi }_{k},$

$$\begin{eqnarray}&&{\beta }_{k}=| \langle {\phi }_{k}| \psi ({T}_{{\rm{p}}}/2)\rangle {| }^{2}\quad {\rm{and}}\quad {\beta }_{{\rm{bound}}}=\displaystyle \sum _{k}{\beta }_{k}.\end{eqnarray} \tag{ 11 }$$

In practice, the sum in (11) includes the first 100 bound states which is enough to ensure convergence of ${\beta }_{{\rm{bound}}}$. Obviously, the ground state depletion is $1-{\beta }_{0}$, and the probability for the electron to be in an excited state is ${\beta }_{{\rm{bound}}}-{\beta }_{0}$. The ionization probability is defined as the non-bound part of the wave function, ${\beta }_{{\rm{free}}}=1-{\beta }_{{\rm{bound}}}$.

Results for the 80-cycle flat-top pulse with ${T}_{{\rm{FWHM}}}=215.8\,\mathrm{fs}$ are shown in figure 1. The nonlinear refractive index (purple) exhibits local drops around peak intensities of 15, 22 and $30\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ as well as above $31\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$. Apparently, these drops coincide with an increased ground state depletion (red). The depletion is mainly caused by enhanced transitions into the continuum (green), since the probability to end up in an excited bound state (orange) is 1−2 orders of magnitude smaller. The origin of the enhanced ionization probability are Freeman resonances, as discussed in [17]. To consolidate the existence of Freeman resonances in the optical response, we briefly sketch how they occur: since the ground state of an atom is strongly bound, its energy remains almost constant when an external field is applied. With increasing intensity, however, the AC Stark effect shifts the energy levels of the dressed Rydberg states by the ponderomotive potential, ${U}_{{\rm{p}}}={E}^{2}/(4{\omega }_{0}^{2})$. Eventually, a Freeman resonance arises when the energy difference between a Rydberg state and the ground state is a multiple of the photon energy. These resonances are usually preceded by a channel closure, which occurs when the ground state is in multiphoton resonance with the ponderomotively up-shifted continuum threshold. In consequence, we expect an enhanced population of the participating Rydberg state at a Freeman resonance. Indeed, above the 9-photon channel closure at $5.5\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ in figure 1, we observe an increase and successive peaks in the occupation probability of the excited states. Interestingly, this seems to be anticorrelated to the ionization probability, which is decreased at the same intensities, most likely due to interference stabilization [30]. This stabilization vanishes for peak intensities $\geqslant 15\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$, where transitions into both bound and continuum states are resonantly enhanced.

Zoom In Zoom Out Reset image size

We now look closer at individual Rydberg states which should be dominantly populated at a Freeman resonance. As shown in figure 2(a), we resolve transitions up to the 31st excited state just after the 9-photon channel closure. Note that parity conservation prohibits a 9-photon transition from the even ground state to even excited states, thus only transitions to odd-parity states are allowed. After the 10-photon channel closure at $30.2\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ in figure 2(d), we again observe Freeman resonances, now with 10-photon transitions to even-parity excited states. Somehow different are the resonances near 15, 22 and $30\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ in figures 2(b) and (c): increasing the peak intensity beyond the range of figure 2(a), we expect transitions to the 7th, 5th and 3rd excited state, since with increasing intensity lower lying excited states come into multiphoton resonance with the ground state. However, these transitions are not dominant since the population of other excited states is of same magnitude. Obviously, multiple excited states contribute to these resonances. This may be related to the fact that the energy shift of lower lying excited states shows a more complicated behavior than that of Rydberg states. Resulting interactions between dressed excited states may lead to a parity change, which explains that also even-numbered excited states are occupied. Nevertheless, these 'mixed resonances' facilitate enhanced ionization as seen in figure 1.

Zoom In Zoom Out Reset image size

The observations of the last paragraph clearly indicated that resonances in the refractive index are correlated to the ground state depletion. Our results also revealed that the ground state depletion is dominated by transitions into the continuum. They therefore confirm the standard model assumption that the refractive index saturation stems from free electrons. A significant difference to the standard model, however, lies in the physical mechanism underlying the ionization enhancement which we identified as resonances between laser-dressed states. Thus, the standard model should be able to reproduce our simulation results if the resonantly enhanced ionization is taken into account. To test this assumption, we write down a modified, time-independent version of the nonlinear refractive index (1)

$$\begin{eqnarray}&&{\rm{\Delta }}n(I)={n}_{2}I-\eta \displaystyle \frac{{\rho }_{\mathrm{fin}}(I)}{2{\rho }_{c}},\end{eqnarray} \tag{ 12 }$$

where I is the peak intensity of the laser pulse. Furthermore, we have to take into account that the density of free carriers is unambiguously defined only for a vanishing external electric field [31]. We therefore replaced in (1) the Keldysh-like ionization density $\rho (t)$ with the electron density at the end of the pulse ${\rho }_{\mathrm{fin}}(I)$ calculated from the direct numerical simulations, multiplied by a correction factor $\eta =0.5$ as justified in the appendix. We then compare three different cases:

$$\begin{eqnarray}{\rho }_{\mathrm{fin}}(I)=\sigma \left\{\begin{array}{ll}{\beta }_{{\rm{free}}}(I) & {\rm{case}}\,{\rm{A}}\\ {\beta }_{{\rm{free}}}(I)+{\beta }_{{\rm{Ryd}}}(I) & {\rm{case}}\,{\rm{B}}\\ 1-{\beta }_{0}(I) & {\rm{case}}\,{\rm{C.}}\end{array}\right.\end{eqnarray} \tag{ 13 }$$

Case A considers only ionized electrons, whereas case B includes the population in Rydberg states. This is motivated by the fact that their contribution to the refractive index is similar to that of free electrons [32]. Additionally, electrons which are freed during the pulse may be captured back into Rydberg states when the field is gone. Here, we consider all states above the ninth excited state as Rydberg states since their corresponding Floquet energies shift with the ponderomotive potential (cf figure 1 in [17]). The probability to find the atom in a Rydberg state is then given by ${\beta }_{{\rm{Ryd}}}={\sum }_{k\gt 9}{\beta }_{k}$. The overall population of excited states, however, is small compared to the ionization probability as seen in figure 1. Therefore, case C considers the full ground state depletion and may serve as a simple, but accurate approximation.

Results are presented in figure 3 together with the exact TDSE calculation (10) for different pulse shapes and lengths. First, we employed the 80-cycle flat-top pulse. For low peak intensities shown in panel (a), the refractive index for both model (A: green, B: orange, C: red) and simulation (purple, dashed) follows the Kerr-nonlinearity ${n}_{2}I$ (gray, dashed–dotted) until ionization becomes non-negligible. However, beyond $8\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$, the TDSE result is obviously superlinear, indicating that the next coefficient in a perturbative description n₄ would be positive. Actually, a superlinear correction proportional to I² renders the agreement between model and simulation astounding close. For the $15\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ resonance, the drop in the refractive index is significantly larger than predicted by model case A and B (which overlap in this panel). Case C, on the other hand, is closer to the simulation result. This observation also holds for resonances at $22\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ and $29.8\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ in panel (b) which shows the full intensity range on a larger scale for the same flat-top pulse. All in all the deviation remains small. After the 10-photon channel closure, the increased population of Rydberg states (cf figure 2) becomes noticeable in the refractive index. As expected, case A (ionization only) deviates from the simulation while case B and C are in good agreement. Model and simulation astonishingly coincide even for a purely ${\cos }^{2}$-pulse as shown in panel (c), where the total pulse duration is 96 cycles (${T}_{{\rm{FWHM}}}=93.3\,\mathrm{fs}$). Again, we notice deviations in case A for peak intensities larger than $30\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$, but case B and C are close to the numerical result. For the 24-cycle pulse (${T}_{{\rm{FWHM}}}=23.3\,\mathrm{fs}$) in panel (d), we observe increasing deviations from our model as the Kerr nonlinearity and the plasma response assume comparable magnitudes, i.e., for intensities beyond the zero-crossing of the nonlinear refractive index. These deviations, however, remain small and the qualitative behavior is reproduced.

Zoom In Zoom Out Reset image size

Overall, our modification of the standard model is in remarkable agreement with the nonlinear refractive index from pure TDSE calculations, manifesting the conclusion that resonantly enhanced ionization is the driving component in its saturation and sign inversion.

Non-Hermitian Floquet theory describes the atom dynamics in terms of Floquet resonances [33]. These are metastable electronic states in a finite volume which are subject to outgoing wave boundary conditions and thus take ionization effects into account. In this approach, the ionization rate is identical to the inverse lifetime of the Floquet ground state resonance. To ensure L²-integrability of the resonance wave functions, we employed complex rotation $x\to {{x}{\rm{e}}}^{{\rm{i}}\theta }$ to the Hamiltonian

$$\begin{eqnarray}&&H(x)=-\displaystyle \frac{1}{2}{\partial }_{x}^{2}-\displaystyle \frac{1}{\sqrt{{x}^{2}+{\alpha }^{2}}}-{\rm{i}}{A}(t)\displaystyle \frac{{\partial }_{x}}{\partial x}.\end{eqnarray} \tag{ 16 }$$

Here, we opted for the velocity gauge, since in the length gauge, complex rotation would result in a blow-up of the wave function during time propagation. The ionization rate and Floquet resonances can be obtained from the Floquet eigenvalue equation [34–36]

$$\begin{eqnarray}&&(H({{x}{\rm{e}}}^{{\rm{i}}\theta })-{\rm{i}}{\partial }_{t})\phi (x,t)={\epsilon }_{\alpha }\phi (x,t),\end{eqnarray} \tag{ 17 }$$

where ${\epsilon }_{\alpha }$ is the complex quasi-energy and ${{\rm{\Gamma }}}_{\alpha }=-2\,\mathrm{Im}\,{\epsilon }_{\alpha }$ is the inverse lifetime of a resonance. The Floquet wave function ${\rm{\Psi }}=\,{{\rm{e}}}^{-{\rm{i}}{\epsilon }_{\alpha }t}\phi$ [37] is quasi-periodic with

$$\begin{eqnarray}&&{\rm{\Psi }}(x,t+T)=\,{{\rm{e}}}^{-{\rm{i}}{\epsilon }_{\alpha }T}{\rm{\Psi }}(x,t),\end{eqnarray} \tag{ 18 }$$

where $T=2\pi /{\omega }_{0}$ is period of the laser pulse. Ψ decays exponentially due to the negative imaginary part of the quasi-energy ${\epsilon }_{\alpha }$ (with ${{\rm{\Gamma }}}_{\alpha }\gt 0$). Typically, the ground state resonance ${\phi }_{0}$ is the most stable eigensolution with the smallest decay rate ${{\rm{\Gamma }}}_{0}$, which can be interpreted as the rate of ionization from the ground state. This is exactly the quantity we are interested in to build our augmented standard model. Instead of solving (17) directly, we note that by substituting ${\rm{\Psi }}(x,t+T)=U(T)\,{\rm{\Psi }}(x,t)$ into (18), it can be rewritten as an eigenvalue equation for the one-cycle propagator ${ \mathcal U }\equiv U(T)$. Consequently, the Floquet wave functions are eigenfunctions of ${ \mathcal U }$, whereas the Floquet multipliers ${{\rm{e}}}^{-{\rm{i}}{\epsilon }_{\alpha }T}$ are the corresponding eigenvalues [38]. For diagonalization of ${ \mathcal U }$ we employ the eigenvectors ${\chi }_{n}$ of the discrete counterpart of the free particle Hamiltonian ${H}_{0}=-{\partial }_{x}^{2}/2$. By propagating ${\chi }_{n}$ along one optical cycle we calculate the matrix elements

$$\begin{eqnarray}&&{{ \mathcal U }}_{{mn}}=\langle {\chi }_{m}(0)| { \mathcal U }| {\chi }_{n}(0)\rangle =\langle {\chi }_{m}(0)| {\chi }_{n}(T)\rangle .\end{eqnarray} \tag{ 19 }$$

Since the results of this section require an accurate determination of the complex quasienergies, for the numerical propagation of the $| {\chi }_{n}\rangle$, we employed a fourth-order scheme for the spatial derivatives based on the Numerov approximation employed in [39]. In fact, this procedure does not affect the tridiagonal property of the short-distance propagator, such that the Crank-Nicolson scheme of section 2.1 can be maintained here. Diagonalization of ${ \mathcal U }$ then yields the complex Floquet multipliers ${{\rm{e}}}^{-i{\epsilon }_{\alpha }T}$ and the corresponding Floquet states ${\phi }_{\alpha }$. Note that in the analytic theory, according to the Balslev−Combes theorem [40, 41], the quasi-energies of the quasi-bound resonances do not depend on the chosen rotation angle θ. However, this does not necessarily hold true for the numerical treatment, since one usually diagonalizes ${ \mathcal U }$ using either a truncated set of basis functions or a finite space grid. In these cases, the rotation angle has to be chosen such that the quasi-energies become stationary with respect to a variation of θ. In our case, we obtain reasonable results for $\theta =0.3$.

In figure 4, we plot the Floquet multipliers ${{\rm{e}}}^{-{\rm{i}}{\epsilon }_{\alpha }T}$ in the complex plane for two different values of θ. The multipliers corresponding to the discrete quasi-continuum can be observed to spiral into the origin. In accordance with the Balslev−Combes theorem, the continuum energies are transformed as ${\epsilon }_{c}\to {\epsilon }_{c}\,{{\rm{e}}}^{-2{\rm{i}}\theta }$, and the positions of the corresponding Floquet multipliers consequently depend on the rotation angle θ. In contrast, the positions of the multipliers corresponding to Floquet resonances are stationary w.r.t. θ-variation. Therefore, a small change in θ allows for a straightforward identification of the Floquet resonances. Furthermore, the distance of the Floquet multipliers to the unit circle is a measure of the lifetime Γ. While the Floquet multipliers of stable bound states (${\rm{\Gamma }}=0$) are located on the contour ${ \mathcal C }$ of the unit circle, the lifetime $\tau =1/{\rm{\Gamma }}$ of a resonance decreases with increasing distance from ${ \mathcal C }$. Going from $I=5\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ in figure 4(a) to $22\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ in (b), resonances move closer to the origin for higher intensities, i.e., they are less stable and ionization is increased. Additionally, in panel (b), we observe resonances with almost the same complex angle in the polar plot. Since a full rotation corresponds to an energy difference of one photon energy, these resonances facilitate multiphoton transitions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To calculate an intrapulse ionization rate ${{\rm{\Gamma }}}_{0}$ for the ground state resonance, we solve the eigenvalue equation (18) for a set of intensities up to $40\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$. When the intensity increases above a channel closure, ${{\rm{\Gamma }}}_{0}$ exhibits pronounced peaks stemming from multiphoton or Freeman resonances, as shown in figure 6(a). In principle, the obtained ionization rate is valid only in the adiabatic regime of a slowly varying pulse envelope. In this case, resonantly enhanced ionization transfers an electron in the ground state resonance directly into the continuum, without populating the intermediate resonant state, as pointed out in [42]. However, as the pulse approaches the few-cycle regime, population transfer to excited Floquet resonances has to be taken into account, and the pronounced resonance peaks in ${{\rm{\Gamma }}}_{0}(I)$ are expected to smear out for shorter pulses.

Zoom In Zoom Out Reset image size

Having calculated the ionization rate ${{\rm{\Gamma }}}_{0}$, we can compare the predictions from our AFR model (14) and (15) with the results of direct numerical simulations. First, we solve the TDSE to obtain the polarization density $P(t)$ for the given electric field $E(t)$. From these quantities, a time dependent susceptibility $\chi (t)$ is computed by forming the ratio of the corresponding complex analytic signals [43]

$$\begin{eqnarray}&&\chi (t)=\displaystyle \frac{{ \mathcal P }(t)}{{\epsilon }_{0}{ \mathcal E }(t)}\end{eqnarray} \tag{ 20 }$$

with ${ \mathcal P }=P+{\rm{i}}{ \mathcal H }(P)$ (and likewise ${ \mathcal E }$), and ${ \mathcal H }$ is the Hilbert transform. From (20), the (ab initio) nonlinear refractive index is obtained as ${\rm{\Delta }}n(t)=\sqrt{1+\chi (t)}-{n}_{0}$, where n₀ is the field-free background refractive index. In order to obtain the corresponding predictions of the AFR model, denoted as ${\rm{\Delta }}{n}^{{\mathsf{AFR}}}(t)$, we have to determine the nonlinear index n₂: for intensities below $1\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$, i.e., in the perturbative Kerr regime, the refractive index reproduces the intensity profile of the pulse, ${\rm{\Delta }}n(t)={n}_{2}I(t)$. After evaluating ${\rm{\Delta }}n(t)$ from TDSE at the local maxima of $I(t)$ for a series of low peak intensities, a linear fit yields ${n}_{2}=6.5\times {10}^{-7}$ ${\mathrm{cm}}^{2}\,\,{\mathrm{TW}}^{-1}$. This allows the evaluation of the Kerr part of the nonlinear refractive index ${\rm{\Delta }}{n}^{{\mathsf{AFR}}}(t)$. The plasma part in (14) can be calculated by plugging the ionization rate ${{\rm{\Gamma }}}_{0}$ into the rate equation (15) and by integrating the latter for temporal intensity profiles $I(t)$.

Figure 5 shows temporally resolved refractive indices for peak intensities 6.2, 22.2 and $38.2\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$ using a cos²-pulse with 96 cycles. If we compare the AFR results (orange) to that of the direct simulations (purple, dashed), we find that our model gives excellent predictions for the nonlinear refractive index: we clearly see that in the leading edge the pulse profile is recovered in accordance with the Kerr nonlinearity ${\rm{\Delta }}n(t)={n}_{2}I(t)$, panel (a) and (b). However, as plasma accumulates during the laser-atom interaction, the refractive index correction in the trailing edge of the pulse becomes more negative and turns the nonlinearity into a defocusing one, panel (b) and (c). This behavior is also well-known from the standard model and responsible for the dynamic spatial replenishment scenario described in [44].

To conclude our comparison with TDSE results, we derive the intensity dependence of the refractive index as predicted by the AFR model. In particular, if we know the polarization density ${P}^{{\mathsf{AFR}}}(t)$, we can employ (9) to calculate the time-averaged susceptibility ${\chi }^{{\mathsf{AFR}}}({\omega }_{0})$ at the center frequency of the pulse. In order to construct ${P}^{{\mathsf{AFR}}}(t)$, we use the general relation

$$\begin{eqnarray}&&P{(t)={\epsilon }_{0}\chi (t)E(t){\rm{and}}\chi (t)=n}^{2}(t)-1,\end{eqnarray} \tag{ 21 }$$

where, in this case, $n(t)={n}_{0}+{\rm{\Delta }}{n}^{{\mathsf{AFR}}}(t)$. The nonlinear refractive index ${\rm{\Delta }}{n}^{{\mathsf{AFR}}}(I)$ then directly follows from (9) upon Fourier transform of ${P}^{{\mathsf{AFR}}}(t)$ and $E(t)$ with various peak intensities I.

In figures 6(b)–(d), we show the intensity-dependent refractive indices obtained from both model and simulation versus peak intensity I: for a flat-top comprising of 80 optical cycles (b) as well as for cos²-envelopes with ${T}_{{\rm{p}}}=96\,T$ (c) and ${T}_{{\rm{p}}}=24\,T$ (d). For the long flat-top pulse, the agreement between the direct numerical calculation and the augmented standard model are indeed excellent. In fact, for long pulses we can employ the adiabatic approximation with respect to the pulse envelope: here, the wave function remains in the ground state resonance which adiabatically adapts to the instantaneous intensity. In this case, the ionization rate ${{\rm{\Gamma }}}_{0}$ accurately describes the actual ionization from the ground state resonance. Moving to shorter pulses in figures 6(c) and (d), however, we see deviations between the two calculations. Especially, above $30\,\mathrm{TW}\,{\mathrm{cm}}^{-2}$, the AFR refractive index exhibits strong resonances, which appear smeared out in the ${\rm{\Delta }}n(I)$ from the direct simulations. Obviously, for shorter pulses, the evolution of the envelope can no longer be considered adiabatic. In consequence, transitions occur from the ground state resonance to higher lying resonances, and the ionization rate ${{\rm{\Gamma }}}_{0}$ becomes less accurate.

To conclude this section, we point out that we cross checked our results for the one-dimensional model atom with results from non-Hermitian Floquet calculations for a full three-dimensional (3D) model of atomic hydrogen, employing the STRFLO eigensolver of [45], which is available for online download [46]. For the 3D case, we found qualitatively comparable resonance structures in the decay rate of the Floquet ground state resonance. These resonances were also present in full 3D TDSE simulations [39] for finite envelope pulses, giving rise to resonantly enhanced ionization yields at the end of the pulses. However, we observed that finite envelope effects, i.e., non-adiabaticity in the pulse envelope, have a more pronounced impact in 3D: the resonances seen by a 90 fs pulse are slightly more smeared out than in the discussed one-dimensional case. We therefore repeated our 3D simulations for slightly longer 120 fs pulses, in case of which the clear visibility of the resonances in the ionization fraction was restored. This leads us to conclude that the qualitative predictions of our model are equally valid for the three-dimensional model atom.