Semiclassical two-step model with quantum input: Quantum-classical approach to strong-field ionization

Strong-field physics is a fascinating field of research resulting from remarkable progress in laser technologies during the last three decades. The interaction of strong laser radiation with atoms and molecules leads to many highly nonlinear phenomena, including above-threshold ionization (ATI) along with the formation of the plateau in the energy spectrum of the photoelectrons (high-order ATI), generation of high-order harmonics of the incident field (HHG), and nonsequential double ionization (NSDI) (see, e.g., Refs. DeloneKrainovBook2000 ; BeckerRev2002 ; MilosevicRev2003 ; FaisalRev2005 ; FariaRev2011 for reviews). The main theoretical approaches used to study all these phenomena are based on the strong-field approximation (SFA) Keldysh1964 ; Faisal1973 ; Reiss1980 , direct numerical solution of the time-dependent Schrödinger equation (TDSE) (see, e.g, Refs. Muller1999 ; Bauer2006 ; Madsen2007 ; Grum2010 ; Patchkovskii2016 ; Tong2017 and references therein), and the semiclassical models.

The semiclassical models apply classical mechanics to describe the motion of an electron after it has been released from an atom or molecule by a strong laser field. The most widely known examples of the semiclassical approaches are the two-step Linden1988 ; Gallagher1988 ; Corkum1989 and the three-step Kulander_Schafer1993 ; Corkum1993 models. The two-step model corresponds to the following picture of ionization process. In the first step an electron is promoted into the continuum, typically by tunneling ionization Dau ; PPT ; ADK . In the second step the electron moves in the laser field towards a detector along a classical trajectory. In addition to these two steps, the three-step model involves the interaction of the returning electron with the parent ion. Accounting for this interaction allows the three-step model to qualitatively describe high-order ATI, HHG, and NSDI.

The semiclassical approaches have important advantages. First, the trajectory-based models, including those that take into account both the laser field and the ionic potential, are often computationally simpler than the numerical solution of the TDSE. What is even more important, the analysis of the classical trajectories helps to understand the physical picture of the strong-field phenomenon under study.

In order to calculate the classical trajectory, it is necessary to specify the corresponding initial conditions, i.e., the starting point and the initial velocity of the electron. To obtain the former, i.e., the tunnel exit point, the separation of the tunneling problem for the Coulomb potential in parabolic coordinates can be used, see, e.g., Ref. Dau . In trajectory-based models it is often assumed that the electron starts with zero initial velocity along the laser field. Simultaneously, it can have a nonzero initial velocity in the direction perpendicular to the field. The initial transverse momenta, as well as the instants of ionization, are usually distributed in accord with the static ionization rate DeloneKrainov1991 with the field strength equal to the instantaneous field at the time of ionization.

In the standard formulation, the trajectory models used in strong-field physics are not able to describe quantum interference effects. Accounting for interference effects in trajectory-based simulations have attracted considerable interest (see, e.g., Refs. SandRost2000 ; Spanner2003 ; Zagoya2012 ; Faria2014 ). The recently developed quantum trajectory Monte-Carlo (QTMC) Li2016 and semiclassical two-step (SCTS) models Shvetsov2016 describe interference structures in photoelectron momentum distributions of the ATI process. These models assign a certain phase to each classical trajectory, and the corresponding contributions of all the trajectories leading to a given asymptotic (final) momentum are added coherently. The QTMC model accounts for the Coulomb potential within the semiclassical perturbation theory. In contrast to this, in the SCTS the phase associated with every trajectory is obtained using the semiclassical expression for the matrix element of the quantum-mechanical propagator (see Ref. Miller1974 ). Therefore, the SCTS model accounts for the binding potential beyond the semiclassical perturbation theory. This explains why for identical initial conditions after the ionization step the SCTS model shows closer agreement with solution of the TDSE than the QTMC model (see Ref. Shvetsov2016 ).

The analysis of the photoelectron momentum distributions and energy spectra calculated within both the QTMC and the SCTS models showed that the ATI peaks are qualitatively reproduced by the semiclassical approaches Shvetsov2016 . However, the semiclassical approximation does not quantitatively reproduce the amplitude of the oscillations. The photoelectron spectra calculated within the semiclassical models fall off too rapidly for energies exceeding $U_{p}$, where $U_{p}=F_{0}^{2}/4\omega^{2}$ is the ponderomotive energy, i.e., the cycle-averaged quiver energy of a free electron in an electromagnetic field (atomic units are used throughout the paper unless indicated otherwise). Here, $F_{0}$ and $\omega$ are the amplitude and the frequency of the field, respectively. This deficiency is closely related to the fact that the initial conditions usually employed in semiclassical models provide too few trajectories with large longitudinal momenta Shvetsov2016 .

Recently several approaches to improving the quality of the initial conditions in semiclassical models have been proposed. The simplest method is to distribute the initial conditions for electron trajectories using the SFA formulas, see, e.g., Refs. Yan2010 ; Popruzhenko2008 ; Boge2013 ; Hofmann2014 ; Geng2014 ; Li2016dop . We note that this method dates back to Refs. Yudin2001 ; Bondar2009 . In most cases it leads to closer agreement with the TDSE. However, to the best of our knowledge, the validity of the usage of the SFA expressions in trajectory-based simulations has not been systematically analyzed so far. In this paper we consider an alternative approach: the combination of the semiclassical models with direct numerical solution of the TDSE.

A significant step in this direction has been taken with the development of the backpropagation method (see Refs. Ni2016 ; Ni2018 ). In this method, the wave packet of the outgoing electron obtained from the TDSE is transformed into classical trajectories. These trajectories are then propagated backwards in time, which makes it possible to retrieve the information about the tunnel exit point and the initial electron velocity. Various approximations to the distributions of the starting points and the initial velocities were analyzed by choosing different criteria to stop the backpropagating trajectories Ni2018 ; Niwe2018 . However, the backpropagation method requires the numerical solution of the TDSE up to some point in time after the end of the laser pulse in the whole space. This restricts its applicability in the case of computationally difficult strong-field problems.

A promising approach would be a combination of the SCTS model with extended virtual detector theory (EVDT), see Refs. Wang2013 ; Wang2018 . For the first time the concept of virtual detector (VD) was proposed in Ref. Thumm2003 as a method for calculating momentum distributions from the time-dependent wave function. The EVDT approach combines the VD method with semiclassical simulations. The EVDT employs a network of virtual detectors that encloses an atom interacting with the external laser field. Each detector detects the wave function $\psi\left(\vec{r},t\right)=A\left(\vec{r},t\right)\exp\left[i\phi\left(\vec{r},t\right)\right]$ obtained by solving the TDSE and generates a classical trajectory at the same position with the initial momentum $\vec{k}$ determined from the gradient of the phase, $\vec{k}\left(\vec{r}_{d},t\right)\equiv\nabla\cdot\phi\left(\vec{r}_{d},t\right)=\vec{j}\left(\vec{r}_{d},t\right)/\left|A\left(\vec{r}_{d},t\right)\right|^{2}$. Here $\vec{r}_{d}$ is the position of a virtual detector and $\vec{j}\left(\vec{r}_{d},t\right)$ is the probability flux at this position. The latter determines the relative weight of the generated trajectory. The subsequent motion of an electron is found from the solution of Newton’s equations. The final photoelectron momentum distribution is obtained by summing over all classical trajectories with their relative weights. It should be stressed that EVDT solves the TDSE only within some restricted region centered at the atom. A network of virtual detectors is placed at the boundary of this region. This reduces the computational load of numerically difficult strong-field problems. Recently the VD approach was used for study of tunneling times Teeny2016 and longitudinal momentum distributions Tian2017 in strong-field ionization.

Leaving the combination of the SCTS with the EVDT for future studies, in this paper we formulate an alternative quantum-classical approach: the semiclassical two-step model with quantum input (SCTSQI). To this end, we combine the SCTS model with initial conditions obtained from TDSE solutions using Gabor transforms. For simplicity, we consider ionization of a one-dimensional (1D) atom. The generalization to the real three-dimensional case is straightforward. The benefit of the 1D model, however, is that potential deficiencies of trajectory models are exposed better and, therefore, it makes the comparison with the fully quantum simulations more valuable.

The paper is organized as follows. In Sec. II we sketch our approach to solve the TDSE, we briefly review the SCTS model, and we formulate our SCTSQI approach. In Sec. III we apply our model to the ionization of a 1D model atom and present comparison with the TDSE results. The conclusions and outlook are given in Sec. IV.

We benchmark our semiclassical model against the results obtained by direct numerical solution of the 1D TDSE and by using the SCTS model. For this reason, before formulating the SCTSQI model and discussing its outcomes, we briefly review the technique used to solve the TDSE and sketch the SCTS model. We define a few-cycle laser pulse linearly polarized along the $x$-axis in terms of a vector-potential:

Here $n$ is the number of the cycles within the pulse, $\varphi$ is the carrier envelope phase, and $\vec{e}_{x}$ is a unit vector. The laser pulse is present between $t=0$ and $t_{f}=\left(2\pi/\omega\right)n$, and its electric field $\vec{F}$ can be obtained from Eq. (1) by $\vec{F}=-\frac{d\vec{A}}{dt}$.

For our numerical examples we use the intensity of $2.01\cdot 10^{14}$ W/cm² ($F_{0}=0.0757$ a.u.) and the wavelength 800 nm ($\omega=0.057$ a.u.). This corresponds to the Keldysh parameter $\gamma=\omega\sqrt{2I_{p}}/F_{0}$ (see Ref. Keldysh1964 ) equal to 0.87. For simplicity, we set the absolute phase of the pulse (1) equal to zero: $\varphi=0$.

We benchmark our SCTSQI approach against the SCTS model and the exact numerical solution of the TDSE. We implement the SCTS by solving Newton’s equation of motion using a fourth-order Runge-Kutta method with adaptive step size Nuref . In order to fully resolve the rich interference structure, we need to use the a momentum-space bin size of $\Delta k_{x}=0.0019$ a.u. For this value of $\Delta k_{x}$ the convergence of the interference oscillations is achieved for an ensemble consisting of $1.2\times 10^{7}$ trajectories. At first, we consider photoelectron momentum distributions. In Fig. 2 (a) we compare the SCTS model with the solution of the TDSE. The TDSE photoelectron momentum distribution has a rather complicated structure. This is due to the fact that the laser pulse used in calculations is neither long nor very short. The side maxima at $k_{x}=-1.35$ a.u. and $k_{x}=1.33$ a.u. are created due to the interference of contributions from times near the central maximum and minimum of the vector potential, respectively, see Fig. 1(c). The central minimum of the vector potential is also responsible for the formation of the maximum at $k_{x}=1.0$ a.u. On the other hand, the ATI peaks in the electron momentum distributions are most pronounced in the range of $k_{x}$ from $-1.0$ a.u. to $-0.25$ a.u. The SCTS model predicts a caustic of the momentum distribution around $k_{x}=0.38$ a.u. For this reason, we normalize the distributions of Fig. 2 (a) to the total ionization yield. Fig. 2 (a) shows that there is only a qualitative agreement between the SCTS approach and the TDSE result. Indeed, the SCTS model underestimates the width of the momentum distribution.

In Fig. 2 (b) we compare the SCTSQI model with the TDSE. In our SCTSQI simulations we have used $N=50$, $x_{\text{max}}=500$ a.u., and $x_{b}=70$ a.u. In order to achieve convergence of the momentum distribution, the bin size was chosen to be $1.5\times 10^{-4}$ a.u., and $n_{p}=10^{6}$ trajectories were launched at every time step of the TDSE propagation. We note that in the mask method it is difficult to achieve full convergence of the TDSE momentum distribution for small momenta. The distribution in the vicinity of $k_{x}=0$ is formed by the slow parts of the electron wave packet. A long propagation time is needed, in order to let these parts reach the absorbing mask, and, therefore, to obtain converged distribution for small $k_{x}$. Thus we do not consider the region of small $k_{x}$ when comparing the SCTSQI with the TDSE. It is clearly seen from Fig. 2 (b) that for $\left|k_{x}\right|\gtrsim 0.15$ a.u. the SCTSQI model provides quantitative agreement with fully quantum-mechanical result. This applies to both the width of the momentum distribution and the positions of the interference maxima (minima). The small remaining discrepancy in the heights of some of the interference maxima is caused by the fact that similar to the SCTS, the SCTSQI model does not account for the preexponential factor of the semiclassical matrix element Miller1974 .

In Fig. 2 (c) we present the photoelecton energy spectra obtained from the SCTS, the solution of the TDSE, and the present SCTSQI model. It is seen that the SCTSQI and the TDSE spectra are almost identical, while the spectrum predicted by the SCTS model falls off to rapidly with the increase of the electron energy. This is a direct consequence of the fact that the SCTS model underestimates the width of the electron momentum distribution, see Fig. 2 (a).

In order to further test the SCTSQI model, we calculate the electron momentum distributions for different positions of the mask $x_{b}$ and fixed $x_{\text{max}}$ of the computational box, see Fig. 3 (a). The distributions corresponding to different values of $x_{b}$ are in good quantitative agreement with each other. The same is also true for momentum distributions obtained for fixed $x_{b}$ and different values of $x_{\text{max}}$, see Fig. 3 (b). Here, we have used the two values $x_{\text{max}}=500$ a.u. and $x_{\text{max}}=200$ a.u. It should be stressed that it is impossible to obtain accurate electron momentum distributions for the small value $x_{\text{max}}=200$ a.u. using the mask method. We also note that for the 1D soft-core Coulomb potential used in this work, the smallest allowed $x_{b}$ should exceed 30-40 a.u., to be outside of the region where the bound-state wave function is localized. Indeed, due to the large number of time steps, even the absorption of a small fraction of the bound-state wave function at each step will result in a severe distortion of the final momentum distribution.

Finally, we check how important is the phase of the factor $G\left(x,p_{x},t\right)$ in Eq. (24). To this end, in Fig. 4 we compare photoelectron momentum distribution calculated using the formula

instead of the Eq. (24). We find that neglecting the phase of the Gabor transform is severe: The SCTSQI distribution cannot even be qualitatively reproduced when using Eq. (25). This result could be expected. Indeed, the factor $G\left(x,p_{x},t\right)$ contains all the information about the quantum dynamics of the absorbed part of the wave packet prior its conversion to the ensemble of classical trajectories. In a sense the $I_{p}t_{0}$ term in the SCTS phase [see Eq. (11)] plays the role of the phase $G\left(t,x,p_{x}\right)$ of the Gabor transform in Eq. (24).

In conclusion, we have developed a trajectory-based approach to strong-field ionization: the semiclassical two-step model with quantum input. In the SCTSQI every trajectory is associated with the SCTS phase and, therefore, the SCTSQI model allows us to describe quantum interference and account for the ionic potential beyond the semiclassical perturbation theory. Furthermore, the SCTSQI corrects the inaccuracies of the SCTS model in treating the tunneling step. This has been achieved by the numerical solution of the TDSE with absorbing boundary conditions in a restricted area of space, applying the Gabor transform to the part of the wave function that is absorbed at each time step, and transforming this absorbed part into classical trajectories. The Gabor transform determines quantum amplitudes assigned to trajectories of the ensemble. Therefore, in the SCTSQI model the initial conditions of classical trajectories are governed by the exact quantum dynamics rather than by the quasistatic or SFA-based expressions as in other semiclassical approaches.

We have tested our SCTSQI model by comparing its predictions with the numerical solution of the 1D TDSE. We have shown that the SCTSQI model yields quantitative agreement with the fully quantum results. This is true not only for the widths of the electron momentum distributions, but also for the positions of the interference maxima and minima. The model can be straightforwardly extended to the three-dimensional case. Most importantly, the SCTSQI circumvents the non-trivial problem of choosing the initial conditions for classical trajectories. This makes the SCTSQI model extremely useful for study of strong-field ionization of molecules.

We are grateful to Professor Lars Bojer Madsen (Aarhus University), as well as to Nicolas Eicke and Simon Brennecke (Leibniz Universität Hannover) for stimulating discussions. This work was supported by the Deutsche Forschungsgemeinschaft (Grant No. SH 1145/1-1).