Ab-initio calculations of laser-atom interactions reveal harmonics feedback during macroscopic propagation

TL;DR

This paper presents a coupled quantum-classical model to simulate ultrashort laser pulse propagation in atomic media, revealing how self-generated harmonics influence ionization and harmonic yield during macroscopic propagation.

Contribution

It introduces a novel ab initio coupled model that captures the feedback of generated harmonics on laser propagation and ionization processes.

Findings

Harmonics up to order 11 significantly affect pulse propagation.

Generated harmonics influence ionization probability.

Feedback mechanisms alter harmonic yields during propagation.

Abstract

We couple the full 3D ab initio quantum evolution of the light pulse polarization in interaction with an atom with a propagation model to simulate the propagation of ultrashort laser pulses over macroscopic dimensions, in the presence of self-generated harmonics up to order 11. We evidence a clear feedback of the generated harmonics on propagation, with an influence on the ionization probability as well as the yield of the harmonic generation itself.

Tables

Table 1. Table 1: The Lyman series

Notation	Transition	$\lambda$ (nm)	Order ($\omega /{\omega }_0$)
$\alpha$-Lyman	$n=2$ to $n=1$	121.5	6.58
$\beta$-Lyman	$n=3$ to $n=1$	102.5	7.80
$\gamma$-Lyman	$n=4$ to $n=1$	97.2	8.23
Limit	$n\to \infty$ to $n=1$	91.15	8.78

Equations

$$\Delta \mathbf{E}-\frac{1}{c^2}{\partial }_t^2\mathbf{E}={\mu }_0({\partial }_t\mathbf{J}+{\partial }_t^2\mathbf{P}).$$

$$\Delta \mathbf{E}-\frac{1}{c^2}{\partial }_t^2\mathbf{E}={\mu }_0{\partial }_t^2\text{mathbfcal}\mathbf{P}.$$

$$k_x^2{\tilde{E}}_z-\frac{{\omega }^2}{c^2}{\tilde{E}}_z={\mu }_0{\omega }^2{\tilde{\mathcal{P}}}_z.$$

$${\tilde{E}}_z=\frac{{\mu }_0{\omega }^2}{k_x^2-\frac{{\omega }^2}{c^2}}{\tilde{\mathcal{P}}}_z.$$

$$\dfrac{1}{k_{x}^{2}-\frac{\omega^{2}}{c^{2}}}=\dfrac{1}{2\frac{\omega}{c}}\Bigg{(}\dfrac{1}{k_{x}-\frac{\omega}{c}}-\dfrac{1}{k_{x}+\frac{\omega}{c}}\Bigg{)}.$$

$$-k_x{\tilde{E}}_z=\frac{\omega }{c}{\tilde{E}}_z+\frac{\omega }{2{ϵ}_{0}c}{\tilde{\mathcal{P}}}_z,$$

$$-\text{i}{k}_x{\tilde{A}}_z={\partial }_x{\tilde{A}}_z=\text{i}\frac{\omega }{c}{\tilde{A}}_z+\frac{1}{2{ϵ}_{0}c}{\tilde{\mathcal{P}}}_z$$

$$\partial_{x}\widetilde{A}_{z}=\textrm{i}\,\Big{(}\dfrac{\omega}{c}-\dfrac{\omega}{v}\Big{)}\,\widetilde{A}_{z}+\dfrac{1}{2\epsilon_{0}c}\,\widetilde{\mathcal{P}}_{z}\ .$$

$$\partial_{x}\widetilde{A}_{z}=\textrm{i}\Big{(}\dfrac{\omega}{c}-\dfrac{\omega}{v}\Big{)}\,\widetilde{A}_{z}+\dfrac{\widetilde{\varmathbb{P}_{z}}+\textrm{i}\omega\epsilon_{0}\rho_{\textrm{at}}\alpha_{0}\widetilde{A}_{z}}{2\epsilon_{0}c},$$

$$\partial_{x}\widetilde{A}_{z}=\textrm{i}\,\Bigg{[}\dfrac{\omega}{c}\Big{(}1+\dfrac{\rho_{\textrm{at}}\alpha_{0}}{2}\Big{)}-\dfrac{\omega}{v}\Bigg{]}\,\widetilde{A}_{z}+\dfrac{1}{2\epsilon_{0}c}\,\widetilde{\varmathbb{P}}_{z}.$$

$$\partial_{x}\widetilde{A}_{z}=\textrm{i}\,\Big{(}\frac{n_{0}}{c}\omega-\dfrac{\omega}{v}\Big{)}\,\widetilde{A}_{z}+\dfrac{1}{2\epsilon_{0}c}\,\widetilde{\varmathbb{P}}_{z}.$$

$${\partial }_x{\tilde{A}}_z(\omega )=\frac{1}{2{ϵ}_{0}c}(\widetilde{{\mathcal{P}}_z}(\omega )-\text{i}\omega {ϵ}_0{\rho }_{\text{at}}\alpha {\tilde{A}}_z(\omega )),$$

$$\text{i}\hbar\frac{\partial}{\partial t}\Psi(\mathbf{r},t)=\Big{[}\frac{\bf{p}^{2}}{2m_{e}}+V(r)-\frac{q}{m_{e}}\mathbf{A}(t)\cdot\text{\bf{p}}\Big{]}\Psi(\mathbf{r},t),$$

$$\text{i}\frac{\partial }{\partial t}\Psi (\mathbf{r},t)=[{\text{H}}_l+\text{D}(t)]\Psi (\mathbf{r},t),$$

$$\Psi (\mathbf{r},t)=\sum _{l=0}^{l_{\text{max}}}\sum _{i=0}^N{c}_i^l(t)\frac{B_i^k(r)}{r}{Y}_l^0(\theta ,\varphi ).$$

$$\text{i}\text{S}\cdot \frac{\partial }{\partial t}\text{c}(t)=({\text{H}}_l+\text{D}(t))\cdot \text{c}(t),$$

$$s_{i,j}=⟨B_j^k|B_i^k⟩={\int }_0^{\infty }{B}_i(r)B_j(r)dr.$$

Full text

I Introduction

The propagation of high-intensity, ultrashort laser pulses gives rise to a wide variety of nonlinear processes, from higher-harmonic generation (HHG) Li et al. (1989); Brabec and Krausz (2000); Agostini and DiMauro (2004); Krausz and Ivanov (2009) to filamentation Chin et al. (2005); Bergé et al. (2007); Couairon and Mysyrowicz (2007); Wolf (2017). Most descriptions of the processes rely on the single atom picture as in the so-called 3-step model Corkum (1993); Schafer et al. (1993). Propagation effects through the target are often ignored or modeled as a macroscopic phase matching between the fundamental and the harmonics Rundquist et al. (1998); Popmintchev et al. (2009), neglecting phenomena like the interplay between the generated harmonics and the non-linear propagation of the high intensity driving pulse.

An exact treatment combining the time-dependent ab-initio Schrödinger equation (TDSE) description of the interaction of the electric field with atoms and the Maxwell solving of the non-linear propagation is extremely complex. For this reason, simplified microscopic models like Strong-Field Approximation (SFA) Lewenstein et al. (1994); Priori et al. (2000); Gaarde et al. (2008), quantum re-scattering Le et al. (2009); Jin et al. (2011), or the metastable state theory Kolesik and Moloney (2013); Bahl et al. (2017) were developed to make the entire numerical processing tractable.

Filamentation and harmonic generation have long been treated independently although they occur at similar intensities and rely on nonlinear atom-field interactions. Still, efficient third harmonic (TH) generation with a conversion efficiency up to a few percent was measured in laser filaments Liu et al. (2011); Gaarde and Couairon (2009); Bergé et al. (2005); Aközbek et al. (2002); Vockerodt et al. (2012); Kolesik et al. (2006), even on distances as short as a few mm Rundquist et al. (1998). Harmonics up to the 23${}^{\text{rd}}$ order have been observed in filaments in argon Steingrube et al. (2011); Vockerodt et al. (2012), at intensities in the $100\text{\,}\mathrm{T}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ range at $800\text{\,}\mathrm{n}\mathrm{m}$.

Recently, the TH and its relative phase were shown to significantly alter the ionization probability in laser filaments at $800\text{\,}\mathrm{n}\mathrm{m}$ Béjot et al. (2014, 2015); Doussot et al. (2016). The TH also affects the nonlinear refractive index Perelomov and Popov (1967), possibly saturating the Kerr effect Béjot et al. (2011); Bache et al. (2012) and affecting phase matching that is crucial to the HHG yield Paul et al. (2003); Gibson et al. (2003); Weerawarne et al. (2015). Harmonics therefore seem to have a non-negligible feedback on the filamentation dynamics. Despite this potential impact, they are often omitted in propagation codes. To the best of our knowledge, no work to date has included harmonics beyond the $5^{th}$.

Contradictory results about the higher-order Kerr effect (HOKE) in air Loriot et al. (2009, 2010), that raised both positive Kolesik et al. (2010); Bache et al. (2012); Petrarca et al. (2012); Richter et al. (2013); Weerawarne et al. (2015) and negative reactions Odhner et al. (2012); Wahlstrand et al. (2012); Köhler et al. (2013); Kosareva et al. (2011), showed the limitations of modeling of filamentation on the basis of a simple empirical parametrization of the medium polarization and ionization Mitrofanov et al. (2015). A comprehensive picture requires to describe the full quantum polarization response of the medium to the electric field Rensink et al. (2014); Köhler et al. (2013); Volkova et al. (2013); Schuh et al. (2016); Lu and Bandrauk (2001); Bandrauk and Shen (1993); Bandrauk et al. (2009), including a saturation or an inversion of the Kerr effect Nurhuda et al. (2002, 2008); Kano et al. (2006); Köhler et al. (2013) as well as a closure of ionization channels due to the Stark effect Béjot et al. (2013); Wiehle et al. (2003); Cormier et al. (2001). Similarly, the impact of various harmonics, including low-order ones, on the generation of the HHG themselves, has recently been evidenced Krüger and Dudovich (2016); Orenstein et al. (2017); Hofmann et al. (2018).

So far, very few ab initio microscopic models are fully coupled with macroscopic propagation codes. In fact, global descriptions of ultrashort pulse propagation including harmonic generation is necessary to overcome the limitations of the above-mentioned simplified models. Such global descriptions have been proposed along the years Lorin et al. (2007, 2011, 2012); Schuh et al. (2014); Lorin et al. (2015). However, the computing requirements limited them to $10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ propagation distance Lorin et al. (2011). Here, we have implemented an ab-initio description of the atom and a semi-classical representation of the interaction, and coupled it to a propagation model over $5\text{\,}\mathrm{m}\mathrm{m}$, a distance in line with standard gas jet interaction lengths in HHG experiments Li et al. (1989); Brabec and Krausz (2000); Agostini and DiMauro (2004); Krausz and Ivanov (2009). By filtering out selected spectral ranges, we showed that harmonics up to the 9${}^{\text{th}}$ impact on both the ionization process and the lower-order harmonics generation and therefore have a significant feedback on the pulse propagation.

II Methods

We performed one-dimensional, time-resolved simulations (1D+1) of the propagation of a laser pulse initially centered at $800\text{\,}\mathrm{n}\mathrm{m}$, with a peak intensity of $50\text{\,}\mathrm{T}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ and a full width at half maximum (FWHM) of $20\text{\,}\mathrm{f}\mathrm{s}$, in atomic hydrogen at atmospheric pressure.

To include the quantum calculation of the polarization in the propagation equation, we start with a real-field propagation equation

[TABLE]

where $\mathbf{E}$ is the electric field, $c$ is the speed of light, the polarization $\mathbf{P}$ is obtained from the bound states of the atomic wavefunction, and the current density $\mathbf{J}$ from its continuum states. We then define $\mathbfcal{P}$ to gather them in a single term in Eq. (1):

[TABLE]

Considering a field linearly polarized along $\mathbf{z}$, propagating along $\mathbf{x}$, we can express this equation in the momentum space as :

[TABLE]

where $k_{x}$ is the wave vector and $\omega$ is the frequency. Then, the field can be described as :

[TABLE]

Considering only the forward-propagating field, we can keep only the right term of:

[TABLE]

In this condition, the propagation equation become:

[TABLE]

$\epsilon_{0}$ being the permittivity of vacuum. We define the vector potential $\mathbf{A}$ such that $\mathbf{E}=-\partial_{t}\mathbf{A}$. In the Fourier space, $\widetilde{E}_{z}=\textrm{i}\omega\widetilde{A}_{z}$, so that:

[TABLE]

Moving to the referential moving at the pulse velocity:

[TABLE]

If we define the polarization as $\widetilde{\mathcal{P}}_{z}=\widetilde{\varmathbb{P}}_{z}+\textrm{i}\omega\epsilon_{0}\rho_{\textrm{at}}\alpha_{0}\widetilde{A}_{z}$, ($\rho_{\textrm{at}}$ being the density of atoms in the propagation medium and $\alpha\simeq 4.593\,\text{a.u.}$ the hydrogen polarizability at 800 nm) we can write the previous equation as:

[TABLE]

then

[TABLE]

By definition $\rho_{\textrm{at}}\alpha_{0}=\chi_{0}^{(1)}$ and $1+\dfrac{\chi_{0}^{(1)}}{2}\approx\sqrt{1+\chi_{0}^{(1)}}=n_{0}$, if $\chi_{0}^{(1)}\leavevmode\nobreak\ \ll\leavevmode\nobreak\ 1$. Therefore, the propagation equation rewrites:

[TABLE]

Since the electric field is linearly polarized along $\mathbf{z}$, the propagation equation in the referential moving at the phase velocity $v_{\text{p}}=c/(1+\frac{\rho_{\text{at}}\alpha}{2})$ reduces to the unidirectional pulse propagation equation (UPPE) Kolesik et al. (2002); Berti et al. (2015):

[TABLE]

The polarization $\mathbfcal{P}$ is deduced from the wavefunction $\Psi(\mathbf{r},t)$ of the atom under the influence of an electric field $\mathbf{E}(t)$, linearly polarized along $\mathbf{z}$. We use the 3D TDSE in the velocity gauge Cohen-Tannoudji et al. (1973):

[TABLE]

where p is the momentum. In atomic units, this equation displays as:

[TABLE]

where $\text{{H}}_{l}={\text{\bf{p}}^{2}}/{2}+V(r)$ is the free Hamiltonian of the system, and $\text{D}(t)=\mathbf{A}(t)\cdot\text{\bf{p}}$. The angular dependence of the wavefunction is expressed in terms of spherical harmonics, while its radial part is discretized on a B-spline basis:

[TABLE]

Injecting this form into Eq. 14, we get the matrix form of the coupled equation system:

[TABLE]

where the matrix S stems from the non-orthogonality of the B-spline functions:

[TABLE]

and $\text{\bf{c}}(t)$ is the projection of $\Psi$ on $Nl_{\textrm{max}}$ dimensions, that will be propagated thanks to a Crank-Nicholson scheme. This scheme involves the matrix $\text{\bf{M}}(t)=\text{\bf{S}}+\text{i}\left(\text{\bf{H}}_{0}+\text{\bf{D}}(t)\right){\delta t}/{2}$, that is highly sparse and bloc-tridiagonal, so that the bi-conjugate gradient method with preconditioning Sleijpen and Fokkema (1993); Sleijpen et al. (1994) is highly efficient on a single processor.

We then solve the full TDSE at each propagation distance and temporal point, using a highly efficient $B$-spline basis Bachau et al. (2001); Cormier and Lambropoulos (1997). In this basis the interaction Hamiltonian matrix is several orders of magnitude sparser as compared to the eigen-basis. Together with the biconjugate stabilized gradient (BiCGStab) algorithm Sleijpen and Fokkema (1993), this reduces the calculation of polarization to a few hours on a single processor, allowing coupling with propagation. Eq. 12 was integrated by using the split-step Fourier method Weideman and Herbst (1986), with a $4^{\textrm{th}}$-order Runge-Kutta scheme. The spatial propagation step is $\textrm{d}x=$15\text{,}\mathrm{\SIUnitSymbolMicro m}\text{/}$$. We used 2048 grid points per optical cycle at the fundamental wavelength ($800\text{\,}\mathrm{n}\mathrm{m}$), corresponding to a time step of $\simeq$ $1\text{\,}\mathrm{a}\mathrm{s}$ and a box size of $450\text{\,}\mathrm{f}\mathrm{s}$. To keep adequate sampling of the wave carrier, we limited the frequency range to $12\ \omega_{0}$.

The knowledge of the atomic wavefunction $\Psi$ then allows to calculate the polarization $\mathcal{P}_{z}$ as $q_{\textrm{e}}\,\rho_{\textrm{at}}\langle\Psi|z|\Psi\rangle$, with $q_{\textrm{e}}$ the electron charge. Note that this polarization includes both free and bound electrons, as their respective contributions are calculated in the same term, without introducing the usual arbitrary distinction between the current of free charges $\mathbf{J}$ and polarization $\mathbf{P}$ (See Eqs. 1 and 2). Unlike phenomenological models, it also intrinsically includes quantum effects like the Stark effect, that strongly affect the ionization via channel closure Béjot et al. (2013).

III Results

Figure 1 shows the electric field spectrum after propagating over $5\text{\,}\mathrm{m}\mathrm{m}$. As expected at $50\text{\,}\mathrm{T}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ Bandrauk et al. (2009), the harmonic yield decays by 1 to 2 orders of magnitude between two consecutive odd orders: 2%, 0.08%, and 0.008% for the 3${}^{\textrm{rd}}$, 5${}^{\textrm{th}}$, and 9${}^{\textrm{th}}$ harmonics, respectively, as also evidenced by the ranges of colorbars in Figs. 2 and 3(c)-(e). The propagated pulse spectrum also includes the well-known Lyman emission lines between $6\ \omega_{0}$ and $9\ \omega_{0}$, especially where the $7^{\textrm{th}}$ harmonic is expected.

These Lyman lines correspond to transitions from excited sates with levels $n>1$ to the ground state $n=1$ (See Table 1). At intensities of several tens of TW/cm2, they are much more intense than the harmonics, and only one order of magnitude below the fundamental. Since the decoherence between excited and fundamental states is not considered in our model, the system continues to oscillate even after the field is switched off. As they travel slower than the fundamental frequency, the Lyman lines would be transferred from the positive (pulse trail) end to the negative (pulse front) end of the temporal box during the split-step calculation, be caught up by the fundamental pulse, and resonantly amplified. To prevent this, the Lyman frequencies are absorbed on the edge of the temporal numerical box with a soft filter (erf function). This artificial damping forbids quantitative analyses of the Lyman lines evolution, and of their interaction with the remaining of the pulse, especially the $7^{th}$ and $9^{th}$ harmonics.

Figure 2 displays the evolution of the temporal shape of the harmonics along the propagation. The fundamental intensity decays linearly till $2.5\text{\,}\mathrm{m}\mathrm{m}$ and then re-grows (Figure 3b). Conversely, harmonics 3, 5, and 7 first grow over approximately $2.5\text{\,}\mathrm{m}\mathrm{m}$, and then decay. In contrast, harmonic 9 reaches a peak after as little as $0.2\text{\,}\mathrm{m}\mathrm{m}$, then decays almost fully till $3\text{\,}\mathrm{m}\mathrm{m}$, before entering in a new growth cycle. The faster rise and decay of the $7^{th}$, and even more of the $9^{th}$ harmonics, can be expected to stem from their partial overlap with the Lymann resonances, that strongly affect both their absorption, and, via Kramers-Kronig relationships, their dispersion.

The evolution of the harmonic temporal pulse shape is more marked for the higher-order harmonics. The $3^{\textrm{rd}}$ harmonic remains rather bell-shaped over more than $3\text{\,}\mathrm{m}\mathrm{m}$, before a temporal shift accompanied with a slight pulse splitting. The $5^{\textrm{th}}$ harmonic self-shortens downto as little as less than $10\text{\,}\mathrm{f}\mathrm{s}$, before splitting during the first mm of propagation and around $4\text{\,}\mathrm{m}\mathrm{m}$. Pulse splitting is even more obvious in the $7^{\textrm{th}}$ harmonic, both around $1\text{\,}\mathrm{m}\mathrm{m}$ and from $4\text{\,}\mathrm{m}\mathrm{m}$ on, with even three sub-pulses at the end of the simulation. The $9^{\textrm{th}}$ harmonic splits from the very first steps of the propagation, and the two sub-pulses merge back after $2.5\text{\,}\mathrm{m}\mathrm{m}$. We attribute the complex evolution of the harmonics pulse shape to their wide band associated with their short duration. As a consequence, group-velocity dispersion is significant within each harmonic. The phase matching and the associated (back-) conversion efficiency are therefore different among spectral slices within each harmonic.

To get more insight into the impact of each harmonics on the propagation and ionization, we applied a low-pass filter with different cut-off frequencies to the polarization, preventing the generation and propagation of wavelengths beyond $2\omega_{0}$, $4\omega_{0}$, $6\omega_{0}$, $10\omega_{0}$, $12\omega_{0}$, i.e., respectively allowing propagation of the fundamental only, and of harmonics up to the $3^{\text{rd}}$, $5^{\text{th}}$, $9^{\text{th}}$, and $11^{\text{th}}$. Note that the effect of the $7^{\textrm{th}}$ harmonics on the propagation has not been investigated due to potential interferences, which cannot be adequately quantified, with the overlapping Lyman transitions lines (See Fig. 1).

The ionization probability strongly depends on the harmonics taken into account in the propagation (Fig. 3(a)). For example, considering harmonics up to the $11^{\text{th}}$ reduces the ionization probability by a factor 2, as compared with simulations considering the fundamental frequency only. The harmonics intensity is insufficient to yield any significant ionization by themselves. They contribute to the ionization via quantum interferences between pathways implying only photons at the fundamental wavelength, and pathways implying, e.g., one harmonic photon, of higher energy, beside fundamental photons, as was evidenced experimentally in the case of the third harmonic Doussot et al. (2015). This effect therefore depends on both the intensity of the harmonics (Fig. 3(c)-(e)) and their phases relative to each other and to the fundamental (Fig. 3(f)), as discussed below in more details. Note that these phases do not depend much on the different spectral truncations. The TH contribution explains the oscillation of the ionization probability over the propagation, as the $3^{\text{rd}}$ harmonic is generated up to $0.8\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ and subsequently decays (Fig. 3(c)).

The fast oscillation of the ionization yield over the first $600\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$ when considering the $9^{\text{th}}$ and $11^{\text{th}}$ harmonics also clearly correlates with that of the intensity of the $9^{\text{th}}$ harmonic, as well as the relative phase of the $9^{\text{th}}$ and $11^{\text{th}}$ harmonic during the propagation (Fig. 1(f)). All in all, depending on the range of harmonics considered in the calculation, the ionization yield can vary by a factor of more than 2 after $2\text{\,}\mathrm{m}\mathrm{m}$ propagation, clearly evidencing that the harmonics are not a byproduct of the propagation, but rather actively contribute to its dynamics, at least up to the $9^{\text{th}}$. In contrast, the $11^{\text{th}}$ harmonics with only $0.002\text{\,}\%$ conversion efficiency has no influence on the ionization yield, although a single 11${}^{\textrm{th}}$ harmonic photon has an energy of $17\text{\,}\mathrm{e}\mathrm{V}$, sufficient to photoionize the hydrogen atom.

The presence of harmonics also affects the evolution of the fundamental (Fig. 3(b)). Without harmonics, the linear decay of the fundamental intensity is due to dispersion, as well as to losses due to the energy spent to excite and/or ionize the atoms. Accounting for harmonics reduces the ionization yield, but as expected it accelerates the intensity decay of the fundamental along the propagation, due to its conversion into harmonics.

Harmonics up to the 9${}^{\textrm{th}}$ therefore substantially impact the propagation and deform the carrier wave, despite the very low conversion efficiencies. As shown in Fig. 4(a), the carrier wave deformation is mostly explained by the contribution of the third harmonic. This is confirmed by Fig. 4(c), where the envelope deformation is strong at the propagation distance ($2\text{\,}\mathrm{m}\mathrm{m}$) where the TH intensity is close to its maximum. In spite of a minimal impact on the carrier wave, harmonics from the 5${}^{\textrm{th}}$ have a substantial influence on ionization.

The feedback of the harmonics on the pulse propagation is also clearly visible on the harmonic generation itself, as well as on their temporal pulse shape (Fig. 3(c)-(e)). Considering a broader spectral range covering the $5^{\text{th}}$ harmonic increases the yield of the $3^{\text{rd}}$ harmonic (Fig. 3(c)). In contrast, considering the $9^{\text{th}}$ harmonics boosts the $3^{rd}$ harmonics while reducing the $5^{\text{th}}$ one (Fig. 3(d)), that saturates already after $0.5\text{\,}\mathrm{m}\mathrm{m}$ of propagation and even decay beyond $3\text{\,}\mathrm{m}\mathrm{m}$. Finally, the $11^{\text{th}}$ harmonic has little influence, if any, on the generation of the harmonics (Fig. 3(c)-(e)), as was the case for ionization. This suggests that, in the considered conditions, the spectral range of the simulation does not need to be extended beyond the $9^{\text{th}}$ harmonics.

The ionization yield also depends on the relative phases $\Delta\phi_{\omega_{n}}=\phi_{\omega_{n}}-n\phi_{\omega_{0}}$ between the $n^{\textrm{th}}$ harmonics and the fundamental field, which modify the interferences between the above-mentioned ionization pathways. This can also be seen in terms of the influence of the harmonics on the carrier wave shape. In-phase (resp. out-of-phase) harmonics lead to a triangular (resp. square) optical cycles as expected from the Fourier series decomposition of triangle- and square-like periodic signals. Figure 3(f) displays $\Delta\phi_{\omega_{n}}$ as a function of the propagation distance. At the beginning of the propagation, the third harmonic is weak (below 1%) and in phase quadrature with the fundamental as expected from a three-photon process. As determined by Doussot et al. Doussot et al. (2016), this results in a lower ionization probability. However, when the third harmonic rises beyond 1%, its contribution to ionization becomes positive, regardless of the phase Doussot et al. (2016). But the monotonous increase of the phase during the propagation due to group-velocity dispersion finally induces back-conversion of the third harmonic, reducing its intensity, and, consequently resulting in a drop in the ionization yield.

The fifth harmonic is out-of-phase with respect to the fundamental field at the beginning of the propagation. This implies that the former is generated by the cascading process combining two photons of the fundamental field with a third harmonic photon: $5\omega_{0}=3\omega_{0}+\omega_{0}+\omega_{0}$.

The relative phases display wiggles that translate into fluctuations of the ionization probability (See Fig. 3(f)). This complex evolution of the phases cannot be described by group-velocity dispersion only (See dotted lines for the 3${}^{\textrm{rd}}$ and 5${}^{\textrm{th}}$ harmonics) Karplus and Kolker (1963). It emphasises that harmonics follow a complex dispersion relationship influenced by the phase with which they are generated along the propagation, by the non-linear refractive index they encounter, as well as by the free electrons, in a process comparable to macroscopic quasi-phase matching relying on dispersion of a mixture of the neutrals and plasma Paul et al. (2003); Gibson et al. (2003).

Our simulations are consistent with experimental results displaying efficient harmonic generation in the propagation of high-intensity laser pulses Liu et al. (2011); Gaarde and Couairon (2009); Bergé et al. (2005); Aközbek et al. (2002); Vockerodt et al. (2012); Kolesik et al. (2006); Steingrube et al. (2011). Furthermore, the observed strong effect of the harmonics on ionization widely generalizes and extends previous experimental and numerical results on the third harmonic Béjot et al. (2014); Doussot et al. (2016). Note that the same approach could be used for, e.g., rare gases, by using a single-active electron potential instead of the purely Coulombian potential of the atomic hydrogen.

IV Conclusion

As a conclusion, a full coupled TDSE-UPPE propagation model has been developed and used to investigate the propagation of an ultrashort, high-intensity laser pulse over $5\text{\,}\mathrm{m}\mathrm{m}$. Simulations over such distance has been made possible by an efficient ab-initio TDSE solver based on $B$-spline basis set. The consideration of harmonics and their relative dephasing up to the $9^{th}$ strongly influences the propagation of ultrashort pulses, and especially their harmonic generation and ionization yield. Our results illustrate that empirical models based on a complex envelope and phenomenological coefficients for ionization and the Kerr effect cannot adequately reproduce the complex behavior of the high-field – atom interaction. The emergence of single attosecond pulses may allow to probe the propagation at the sub-fs scale, ie, at the sub-optical-cycle scale Sommer et al. (2016), therefore providing direct experimental comparison with our calculations, especially if the propagation is considered in two dimensions. Such $2D+z$ extension, implementing a transverse resolution, together with radial symmetry and efficient parallelization of the propagation steps, will allow its application to filamentation.

Acknowledgements.

We warmly thank M. Moret and the Computer resource center (CRI) of the University of Burgundy for technical and numerical support, as well as J. Doussot for fruitful discussion. This work was supported by the ERC advanced grant "filatmo" and the Conseil Régional de Bourgogne (PARI program).