Symphony on Strong Field Approximation

TL;DR

This paper reviews and extends the Strong Field Approximation (SFA) for intense laser-matter interactions, focusing on multi-electron atoms and molecules, and introduces new time-dependent models for complex systems.

Contribution

The paper presents novel generalizations of SFA to time-dependent multi-electron and molecular systems, enhancing understanding of strong field phenomena.

Findings

Extended SFA to two-electron atoms with EII and RESI dynamics

Formulated SFA for large molecules with ongoing electronic dynamics

Developed analytical models for molecular electronic wave packets

Abstract

This paper has been prepared by the Symphony collaboration (University of Warsaw, Uniwersytet Jagiello\'nski, DESY/CNR and ICFO) on the occasion of the 25th anniversary of the "simple man's models" which underlie most of the phenomena that occur when intense ultrashort laser pulses interact with matter. The phenomena in question include High-Harmonic Generation, Above-Threshold Ionization, and Non-Sequential Multielectron Ionization. "Simple man's models" provide, both an intuitive basis for understanding the numerical solutions of the time-dependent Schr\"odinger equation, and the motivation for the powerful analytic approximations generally known as the Strong Field Approximation (SFA). In this paper we first review the SFA in the form developed by us in the last 25 years. In this approach SFA is a method to solve the TDSE using a systematic perturbation theory in a part of the Hamiltonian describing continuum-continuum transitions in the presence of the laser field. In this review we focus on recent applications of SFA to HHG, ATI and NSMI from multi-electron atoms and from multi-atom. The main novel part of the presented theory concerns generalizations of SFA to: (i) time-dependent treatment of two-electron atoms, allowing for studies of an interplay between Electron Impact Ionization (EII) and Resonant Excitation with Subsequent Ionization (RESI); (ii) time-dependent treatment in the single active electron (SAE) approximation of "large" molecules and targets which are themselves undergoing dynamics during the HHG or ATI process. In particular, we formulate the general expressions for the case of arbitrary molecules, combining input from quantum chemistry and quantum dynamics. We formulate also theory of time-dependent separable molecular potentials to model analytically the dynamics of realistic electronic wave packets for molecules in strong laser fields.

Tables

Table 1. Table 1: The operating conditions in terms of photon energy, ℏ ​ ω Planck-constant-over-2-pi 𝜔 \hbar\omega , ionization potential, I p subscript 𝐼 𝑝 I_{p} , and ponderomotive energy, U p subscript 𝑈 𝑝 U_{p} are presented for single-photon ionization (SPI), multi-photon ionization (MPI), above-threshold ionization (ATI) and tunnel ionization (TI).

Ionization regime		Operating condition
Single-photon ionization (SPI)		$\hslash \omega >I_p\gg U_p$
Multi-photon ionization (MPI)		$I_p>\hslash \omega \gg U_p$
Above-threshold ionization (ATI)		$I_p>U_p>\hslash \omega$
Tunnel ionization (TI)		$U_p>I_p>\hslash \omega$

Equations

$$\gamma =\sqrt{\frac{I_p}{2U_p}},$$

$$U{}_{\mathrm{p}}=\frac{e^2{\mathcal{E}}_0^2}{4m_e{\omega }_0^2}=\frac{I_0{e}^2{\lambda }_0^2}{8{\pi }^2{m}_{\mathrm{e}}{ϵ}_0{c}^3}=9.337\times 10^{-20}{I}_0{\lambda }_0^2\frac{\mathrm{eV}}{\mathrm{W}\mathrm{c}{\mathrm{m}}^{-2}\mathrm{n}{\mathrm{m}}^2},$$

$$\text{E}(t)={\mathcal{E}}_{0}f(t)\sin ({\omega }_{0}t+{\varphi }_0){\mathbf{e}}_z.$$

$$i\hslash \frac{\partial }{\partial t}|\Psi (t)⟩=\hat{H}|\Psi (t)⟩,$$

$$\hat{H}$$

$${\hat{H}}_0$$

$${\hat{H}}_0|\text{p}⟩=\frac{1}{2m}{\text{p}}^2|\text{p}⟩$$

$$|\Psi(t)\rangle=e^{\textit{i}I_{p}\textit{t}/\hbar}\bigg{(}a(t)|0\rangle+\>\int{\textit{d}^{3}\textbf{p}\>\textit{b}(\textbf{p},t)|\textbf{p}\rangle}\bigg{)}.$$

$$e⟨\mathbf{p}|\hat{\mathbf{r}}|{\mathbf{p}}^{\prime }⟩=ie\hslash {\nabla }_{\mathbf{p}}\delta (\mathbf{p}-{\mathbf{p}}^{\prime })+\hslash \mathbf{g}(\mathbf{p},{\mathbf{p}}^{\prime }),$$

$$\dot{a}(t)$$

$$\dot{b}(\text{p},t)$$

$$-e\mathbf{E}(t)\cdot {\nabla }_{\mathbf{p}}b(\mathbf{p},t)+i\text{E}(t)\cdot \int {\text{d}}^3{\text{p}}^{\prime }\text{b}({\text{p}}^{\prime },t)\mathbf{g}(\text{p},{\text{p}}^{\prime }).$$

$$e⟨\text{p}|\hat{\text{r}}|0⟩=\text{d}(\text{p}).$$

$$e⟨\mathbf{p}|\hat{\mathbf{r}}|{\mathbf{p}}^{\prime }⟩=ie\hslash {\nabla }_{\mathbf{p}}\delta (\mathbf{p}-{\mathbf{p}}^{\prime })+\hslash \mathbf{g}(\mathbf{p},{\mathbf{p}}^{\prime }),$$

$${\partial }_t{b}_0(\text{p},t)=-\frac{i}{\hslash }\left(\frac{{\text{p}}^2}{2m}+I_p\right)b_0(\text{p},t)+\frac{\text{i}}{\hslash }\text{E}(t)\cdot \text{d}(\text{p})a(t)+\text{E}(t)\cdot {\nabla }_{\mathbf{p}}{\text{b}}_0(\text{p},t).$$

$$\begin{array}{rl}{b}_0(\text{p},t)=& \frac{i}{\hslash }{\int }_0^t\text{d}{\text{t}}^{\prime }\text{E}(t^{\prime })\cdot \text{d}\left(\text{p}+e\text{A}(t)/c-e\text{A}(t^{\prime })/c\right)\\ & \times \exp \left(-\text{i}{\int }_{t^{\prime }}^{t}d\tilde{t}\left[\frac{1}{2m}(\text{p}+e\text{A}(t)/c-e\text{A}(\tilde{t})/c{)}^2+I_p\right]/\hslash \right)a(t^{\prime }).\end{array}$$

$$b_0(\text{p},t)=\frac{i}{\hslash }{\int }_0^t\text{d}{\text{t}}^{\prime }\text{E}(t^{\prime })\cdot \text{d}\left(\text{p}-e\text{A}(t^{\prime })/c\right)a(t^{\prime })\exp \left(-\text{i}{\int }_{t^{\prime }}^{t}d\tilde{t}\left[\frac{1}{2m}(\text{p}-e\text{A}(\tilde{t})/c{)}^2+I_p\right]/\hslash \right).$$

$$S(\mathbf{p},t,t^{\prime })={\int }_{t^{\prime }}^{t}d\tilde{t}\left[\frac{1}{2m}(\mathbf{p}-e\text{A}(\tilde{t})/c{)}^2+I_p\right].$$

$$\begin{array}{rl}{b}_1(\text{p},t)=& {\left(\frac{i}{\hslash }\right)}^2{\int }_0^t\text{d}{t}^{\prime }\exp \left[-\text{i}S(\mathbf{p},t,t^{\prime })/\hslash \right]\text{E}(t^{\prime })\cdot {\int }_0^{t^{\prime }}\text{d}{\text{t}}^{\prime \prime }\int {\text{d}}^3{\text{p}}^{\prime }\text{g}\left(\text{p}-e\text{A}(t^{\prime })/c,{\text{p}}^{\prime }-e\text{A}(t^{\prime })/c\right)\\ & \times \text{E}(t^{\prime \prime })\cdot \text{d}\left({\text{p}}^{\prime }-e\text{A}(t^{\prime \prime })/c\right)a(t^{\prime \prime })\exp \left[-\text{i}S({\mathbf{p}}^{\prime },t^{\prime },t^{\prime \prime })/\hslash \right].\end{array}$$

$$\text{E}(t^{\prime \prime })\cdot \text{d}\left({\text{p}}^{\prime }-e\text{A}(t^{\prime \prime })/c\right)a(t^{\prime \prime }).$$

$$|b(\mathbf{p},t_{\mathrm{F}}){|}^2$$

$$=|b_0(\mathbf{p},t_{\mathrm{F}}){|}^2+|b_1(\mathbf{p},t_{\mathrm{F}}){|}^2+b_0(\mathbf{p},t_{\mathrm{F}})b_1^{\ast }(\mathbf{p},t_{\mathrm{F}})+c.c.$$

$$\begin{array}{rl}⟨\mathbf{r}(t)⟩& =\mathrm{Re}[\frac{i}{\hslash }{\int }_0^t\text{d}{\text{t}}^{\prime }\int d^3\mathbf{p}{a}^{\ast }(t){\text{d}}^{\ast }\left(\text{p}-e\text{A}(t)/c\right)\text{E}(t^{\prime })\cdot \text{d}\left(\text{p}-e\text{A}(t^{\prime })/c\right)a(t^{\prime })\\ & \times \exp \left(-\text{i}{\int }_{t^{\prime }}^{t}d\tilde{t}\left[\frac{1}{2m}(\text{p}-e\text{A}(\tilde{t})/c{)}^2+I_p\right]/\hslash \right)],\end{array}$$

$$\tilde{\mathbf{r}}(\Omega )$$

$$=\mathrm{Re}[\frac{i}{\hslash }{\int }_{-\infty }^{\infty }dt{\int }_0^t\text{d}{\text{t}}^{\prime }\int d^3\mathbf{p}{a}^{\ast }(t){\text{d}}^{\ast }\left(\text{p}-e\text{A}(t)/c\right)\text{E}(t^{\prime })\cdot \text{d}\left(\text{p}-e\text{A}(t^{\prime })/c\right)a(t^{\prime })$$

$$\times \exp \left(-\text{i}{\int }_{t^{\prime }}^{t}d\tilde{t}\left[\frac{1}{2m}(\text{p}-e\text{A}(\tilde{t})/c{)}^2+I_p\right]/\hslash +i\Omega t\right)],$$

$$\exp \left(-i{S}_{\Omega }(\text{p},t,t^{\prime })/\hslash \right)=\exp \left(-\frac{\text{i}}{\hslash }{\int }_{t^{\prime }}^t\left[\frac{1}{2m}(\text{p}-e\text{A}(\tilde{t})/c{)}^2+I_p\right]d\tilde{t}+i\Omega t\right),$$

$$\frac{\partial }{\partial t^{\prime }}{S}_{\Omega }(\text{p},t,t^{\prime })=\frac{1}{2m}(\text{p}-e\text{A}(t^{\prime })/c{)}^2+I_p$$

$$\frac{\partial }{\partial t}{S}_{\Omega }(\text{p},t,t^{\prime })+\Omega =\frac{1}{2m}(\text{p}-e\text{A}(t)/c{)}^2+I_p$$

$${\nabla }_{\text{p}}{S}_{\Omega }(\text{p},t,t^{\prime })=\frac{1}{m}{\int }_{t^{\prime }}^t\left[\text{p}-e\text{A}(\tilde{t})/c\right]d\tilde{t}$$

$$\tilde{\mathbf{r}}(\Omega )$$

$$|\psi (t)⟩=\exp (i{E}_{0}t/\hslash )\left(a(t)|0⟩+\int d^3\mathbf{p}b(\mathbf{p},t)|\mathbf{p},0⟩+\sum _{\eta }\int d^3\mathbf{p}c(\mathbf{p},\eta ,t)|\mathbf{p},\eta ⟩+\iint d^3\mathbf{p}{d}^3{\mathbf{p}}^{\prime }d(\mathbf{p},{\mathbf{p}}^{\prime },t)|\mathbf{p},{\mathbf{p}}^{\prime }⟩\right),$$

$${\hat{H}}_0=\sum _{i=1}^2\left(\frac{{\hat{\mathbf{p}}}_i^2}{2m}+V({\hat{\mathbf{r}}}_i)\right)+V({\hat{\mathbf{r}}}_1-{\hat{\mathbf{r}}}_2),$$

$${\hat{H}}_0|0⟩$$

Full text

We dedicate this work to the memory of Bertrand Carré, who passed away in March 2018 at the age of 60.

Accepted Manuscript for Rep. Prog. Phys. 82 no. 11, 116001 (2019), arXiv:1812.11447, available under the CC BY-NC-ND license.

I Introduction

I.1 Motivation

In the last three decades, we have witnessed incredible advances in laser technology and in the understanding of nonlinear laser-matter interactions, crowned recently by the award of the Nobel prize to Gérard Mourou and Donna Strickland Strickland1985 ; CPA-Nobel-2018 . It is now routinely possible to produce few-cycle femtosecond (1\text{\,}\mathrm{f}\mathrm{s}$=${10}^{-15}\text{\,}\mathrm{s}) laser pulses in the visible and mid-infrared regimes Nisoli1997 ; Hemmer2013 . By focusing such ultrashort laser pulses on gas or solid targets, possibly in the presence of nano-structures Ciappina2017 , the targets are subjected to an ultra-intense electric field, with peak field strengths approaching the binding field inside the atoms themselves. Such fields permit the exploration of the interaction between strong electromagnetic coherent radiation and an atomic or molecular system with unprecedented spatial and temporal resolution Krausz2009 . On one hand, HHG nowadays can be used to generate attosecond pulses in the extreme ultraviolet Drescher2001 ; Chang2012 , or even in the soft X-ray regime Silva2015 . Such pulses themselves may be used for dynamical spectroscopy of matter; despite carrying modest pulse energies, they exhibit excellent coherence properties Pascalcoh ; AnneML . Combined with femtosecond pulses they can also be used to extract information about the laser pulse electric field itself Paulus2001 . HHG sources therefore offer an important alternative to other sources of XUV and X-ray radiation: synchrotrons, free electron lasers, X-ray lasers, and laser plasma sources. Moreover, HHG pulses can provide information about the structure of the target atom, molecule or solid Itatani2004 ; Blaga2012 ; Pullen2015 . Of course, to decode such information from a highly nonlinear HHG signal is a challenge, and that is why a possibly perfect, and possibly “as analytical as possible” theoretical understanding of these processes is in high demand. Here is the first instance where SFA offers its basic services.

In HHG, an electron is liberated from an atom or molecule through ionization, which occurs close to the maximum of the electric field. Within the oscillating field, the electron can thus accelerate along oscillating trajectories, which may result in recollision with the parent ion, roughly when the laser field approaches a zero value. Since electronic motion is governed by the waveform of the laser electric field, an important quantity to describe the electric field shape is the so-called absolute phase or carrier-envelope phase (CEP). Control over the CEP is paramount for extracting information about electron dynamics, and to retrieve structural information from atoms and molecules Baltuska2003 ; Itatani2004 ; Wolter2014 . Control over the CEP is particularly important for HHG, when targets are driven by laser pulses comprising only one or two optical cycles. In that situation the CEP determines the relevant electron trajectories, i.e. the CEP determines whether emission results in a single or in multiple attosecond bursts of radiation Baltuska2003 ; Sola2006 .

The influence of the CEP on electron emission is also extremely important. It was demonstrated for instance in an anti-correlation experiment, in which the number of ATI electrons emitted in opposite directions was measured Paulus2001 ; Milosevic2006 . Since the first proof-of-principle experiment Paulus2001 , the stereo ATI technique has established itself as a direct measure of the CEP, and demonstrated its ability for single-shot measurements even at multi-kHz laser repetition rates. Both bound-free and the rescattering continuum-continuum transitions are CEP sensitive; hence, the photoelectron distribution of ATI can also be used to extract structural information about the target. Again, theoretical understanding of these processes and the SFA, which are “as analytical as possible”, are more than welcome in this line of research.

Laser-induced electron diffraction (LIED) is a technique that uses the doubly-differential elastic scattering cross-sections to extract structural information Zuo1996 ; Lein2002 ; Lin2010 . Meeting the requirements to extract structural information has, however, proven difficult due to the stringent prerequisites on the laser parameters. During recent years, the development of new laser sources has dramatically advanced, leading to the first demonstrations of the technique Xu2014 ; Meckel2008 ; Morishita ; Xu2012 ; Pullen2015 , and the successful retrieval of the bond distances in simple diatomic molecules with fixed-angle broadband electron scattering Xu2014 . Recently, Pullen et al. Pullen2015 have exploited the full double differential cross section to image the entire structure of a polyatomic molecule for the first time. Again, to exploit the full potential of the recollision physics and the intrinsic time resolution of LIED (and eventually HHG), we need the comprehensive and complete understanding of the ATI and HHG processes and its theoretical, possibly analytical description Milosevic2006 ; Faisal1973 ; Reiss1980 ; Lewenstein1995 ; Starace2014 ; Becker2002 ; Morishita2010 . Here is the third instance where SFA is indispensable.

The key for gaining dynamical structural information and for realizing nonlinear dynamical spectroscopy with HHG, ATI and, last but not least, NSMI, however, consists in generalizations of the existing theories to the case, when the target in question itself undergoes dynamics beyond its SAE electronic structure. For NSMI this requires including two, or even more electrons in the SFA theory. For molecules that would mean, for instance, developing theory that takes into account vibration or dissociation processes, occurring on the time scale comparable with the laser pulse duration. Ab initio theory of this kind is generally too computationally intensive for numerical simulation. Generalizations of the SFA to include many electrons and/or nuclear motion are thus more than welcome.

I.2 History of the SFA

The history of the SFA has at least five beginnings, many intertwined and entangled paths, and many interfering endings. The point of view presented here is dominated by the opinion of one of us (M.L.), who started to work on the subject in 1980. The five beginnings of SFA, according to this point of view, are the following:

• •

Keldysh ionization theory. In 1964 L.V. Keldysh Keldysh1965 attempted and solved analytically a very fundamental problem of tunnelling ionization of atoms in strong low-frequency oscillating electric fields. In the epochal paper entitled “Ionization in the Field of a Strong Electromagnetic Wave” he derived his famous Bessel function formulae for ionization rates. F.H.M. Faisal Faisal1973 and H.R. Reiss Reiss1980 , building also on earlier work by Perelomov, Popov and Terent’ev Perelomov1966 ; Perelomov1967 ; Perelomov1967b , extended this theory to calculate electron spectra in the process to be named Above-Threshold ionization (ATI). These results were the first instances of the Strong Field Approximation, named “KFR theory” after the authors. Note that in this initial phase, SFA was exclusively applied to ionization problems. The first experiments on ATI date from 1979 Agostini1979 ; in fact, the name of the observed phenomenon slowly shifted from “continuum-continuum transitions” to ATI Kruit1981 ; Fabre1982 ; Kruit1983 . The first observations of HHG date from 1988 Ferray1988 ; McPherson1987 , and of non-sequential two-electron ionization originate from 1983 [Anne1983a, ; Anne1983b, ; for more recent reviews see Mainfray1991, ; Krausz2007, ]. Keldysh theory was further developed to calculate total ionization rates for various atomic species and states by M.V. Ammosov, N.B. Delone and V.P. Krainov Ammosov1986 —the resulting expressions are known as ADK rates. More broadly, Keldysh theory has been an inspiration for many years for many scholars. M.L. first learned the Keldysh theory from a preprint by L. Davidovich et al., ICTP Trieste. In fact, Davidovich later published several interesting papers on Keldysh theory Brandi1982 ; Neto1990 ; Becker1994 . In the beginning of the 1980s several non-perturbative, quantum-optics-inspired models of “continuum-continuum transitions” and ATI were introduced by Z. Białynicka-Birula Bialynicka1983 , J.H. EberlyDeng1984 , K. Rzążewski Rzazewski1985 and others (for a review see Ref. EberlyRev1991, ). These models stimulated M.L. to try to combine them with the Keldysh theory.

• •

Kroll-Watson theory. Another inspiration for the contemporary SFA came from the seminal papers of N.M. Kroll and K.M. Watson on electron-atom Kroll1973 and atom-atom Kroll1976 scattering in the presence of a strong electromagnetic wave. The earlier paper clearly dealt with the problem of “continuum-continuum transitions”, dressed by the laser field, leading to the expected Bessel-function dependence of the corresponding transition amplitudes. This observation led to the formulations of the ATI theory as a theory of multichannel decay and continuum-continuum transitions, dressed by the laser field Lewenstein1985 . This approach, employing the relation between ATI and electron scattering in the intense laser field has been deepened and developed further in Ref. Grochmalicki1986, . In the contemporary language, the results of these studies described ATI as a combination of direct tunnelling processes, and rescattering processes occurring in the laser dressed continuum. At that time, however, the underlying quasi-classical theory and the simple man’s model was yet not known. It is worth noting that this approach was also applied to two-electron ionization in Ref. Lewenstein1986a, , where the direct two-electron tunnelling processes were analysed.

• •

Numerical studies of TDSE. Numerical simulations always played, play and will play a fundamental role in our understanding of physics of matter in intense laser fields. A particularly important role was played here by 1D models of one- and two-electron systems, initiated by J.H. Eberly on “Eberlonium”, also known as the Rochester atom model. This series of studies allowed the description of several qualitatively important results but, more importantly, it allowed—by appropriate tuning of the parameters—the finding of accurate quantitative predictions concerning ATI Javanainen1988 , HHG Eberly1989 ; Eberly1989a , stabilization of a 1D atom in the presence of a strong field of high frequency Su1990 —all of that optimizing the “model atom for multiphoton physics” Su1991 . This approach was very successfully generalized to 1D two-electron models Grobe1992 ; Grobe1993 , which in turn stimulated the development of other quasi-1D approaches to the two-electron problems in intense laser fields. These developments were very important, especially in a view of the difficulties and computational cost of solving TDSE for helium in 3D [Dundas1999, ; Taylor2003, ; Parker2000, ; Emma2011, ; see also Feist2008, ; Pazourek2012, ]. Analysis of the classical pathways for simultaneous escape of two electrons showed that there are two saddle points in the effective potential landscape located symmetrically with respect to the field polarization axis Sacha2001 . This led to a modification of the 1D model, where electrons move along axes inclined symmetrically with respect to the polarization direction Eckhardt2006 ; Prauzner2007 ; Prauzner2008 . Within this model the ionization for three active electrons was also recently considered Thiede2018 . Another model, in which the center of mass movement was restricted to the polarization axis was introduced by Ruiz et al. Ruiz2006 , and successfully applied to momentum distributions Staudte2007 —for a comparison of various quasi-1D approaches, see Ref. Efimov2018, .

Nevertheless, the most important were investigations of the TDSE in 3D, led in those years by K. Kulander, who not only developed codes for solving the TDSE, but also propagation and phase matching in HHG, and collaborated intensively with the top experimental groups. Several seminal papers were written on ATI Schafer1993 , double ionization of helium Fittinghoff1992 ; Walker1994 and phase matching in HHG Lhuillier1991 . The one that was the most important for the formulation of the simple man’s models was the theory paper on HHG from atoms and ions in the high-intensity regime Krause1992 , in which the famous cut-off law for HHG was discovered: the high harmonics cease to exist above the photon energy $I_{p}+3.2U_{p}$, where $I_{p}$ is the ionization potential, and $U_{p}$ is the ponderomotive energy.

• •

Classical phase-space methods. The key idea in these approaches was to mimic the evolution of the electronic system in terms of classical trajectories, governed by the completely classical Hamiltonians, but satisfying an initial phase-space distribution. In fact, that is why some versions of this approach were termed the “truncated Wigner function” approach. Initially, these methods were proposed for microwave perturbations of atomic systems by J. Leopold and I.C. Percival Leopold ; Percival . They gained a lot of attention and popularity in the studies of quantum chaos and quantum dynamical localization Casati1987 . Later, these methods were extended to the regime of mid-infrared, and even high laser frequencies and two-electron systems. There are two variants of these methods that can be distinguished: those where the initial distribution is calculated classically Groch1991 ; Gajda1992 ; Rzazewski1994 ; Wojcik1995 ; Ho2007 , and those which use the below-barrier tunneling approximations for the calculation of the initial distribution Chen2000 . The quality of classical phase-space averaging methods can be checked by comparing their results with those from the corresponding quantum-mechanical models Shvetsov2016 . Considerable progress was made with these methods on the study of the ionization yields Ho2005 ; Mauger2009 and momentum distributions Panfili2001 . Still, getting HHG spectra within classical methods exclusively seems to be a complicated task—for a recent discussion see Ref. Berman2018, . These approaches are related to the so-called Initial Value Representation methods used in quantum chemistry VanDeSand1999 ; Zagoya2012 ; Zagoya2014 ; Symonds2015 , though those methods are not purely classical, since they also account for phase shifts.

• •

Simple man’s models. There were simple man’s models before the simple man’s model. An extended discussion of precursors is contained in a recent review in Ref. Ivanov2014, , where the earlier quantum formulation of “Atomic Antennas” of M.Y. Kuchiev is discussed Kuchiev1987 , as well as early attempts by F. Brunel and P. Corkum himself Brunel1987 ; Brunel1990 ; Corkum1989 . Essentially the same formulae as the ones derived later in the framework of the SFA for HHG were obtained by F. Ehlotzky Ehlotzky1992 , but without the underlying quasi-classical picture. The history of science chooses its own heroes. Nowadays, the simple man’s model, also known as the “three-step” or “recollision” model, is usually attributed to P. Corkum Corkum1993 , K. Kulander Krause1992 ; Kulander1993 , and to a conference contribution of H. Muller. These formulations were done in the right place in the right time, and were truly seminal—they have revolutionized the whole area! M.L. learned about the simple man’s model for HHG during and immediately after the famous NATO Workshop on Super Intense Laser Atom Physics (SILAP) in Han-sur-Lesse, Belgium, in January of 1993 footnote . After the workshop M. Yu. Ivanov visited Saclay and stayed at M.L.’s house—that is how our first paper on the SFA for HHG, based on simple man’s model, and co-authored by A. L’Huillier and other colleagues, was written Anne1993 . After a long fight with Phys. Rev. Lett., this paper was published as a Rapid Communication in Phys. Rev. A, entitled “High-order harmonic-generation cutoff”. We termed the formula we used to evaluate the time-dependent dipole moment a dynamical Landau-Dykhne formula.

M.L. went to JILA in February 1993 and started to work on the long version of the theory, including a detailed discussion of the relation of the simple man’s model with the quasi-classical (better termed quantum-orbit) saddle point approximations, along with concrete calculations for what we called Gaussian models, i.e. models in which the ground state of the atom of interest was approximated by a Gaussian function. The paper on the theory of HHG by low-frequency laser fields appeared in Phys. Rev. A in 1994, and soon became a reference paper for theorists and experimentalists working in the field Lewenstein1994 . During his stay in JILA M.L. worked with K. Kulander on the extension of the newly-formulated version of the SFA to ATI, stimulated by the observation of the intensity-dependent rings in the high order ATI Yang1993 . The paper that extended a quasi-classical analysis (i.e. one based on simple man’s model) of rescattering processes in ATI appeared in Phys. Rev. A in 1995 Lewenstein1995 . On one hand, it explicitly demonstrated in which sense the SFA is a systematic perturbation theory in part of the Hamiltonian describing the continuum-continuum transitions. On the other hand, we introduced here for the first time the model atom involving a separable (non-local) potential. This kind of approach was a generalization of the so-called zero-range Becker’s model Becker1994a . Recently, it turned out to be extremely useful in modelling HHG and ATI in molecular dimers, trimers and quadrimers Noslen1 ; Noslen2 ; Noslen3 ; Noslen4 ; Noslen-diss . Other uses of separable potentials in the literature are discussed in Refs. Faisal1987, ; Faisal1987b, ; Faisal1989, ; Faisal1990, ; Faisal1991, ; Faisal1993, ; Faisal1993b, ; Radozycki1993, ; Kaminski1989, ; Klaiber2005, ; TetchouNganso2007, ; TetchouNganso2011, ; TetchouNganso2013, ; Galstyan2017, ; Galstyan2018, .

The quantum simple man’s models, as the novel SFA was termed sometimes, proved very useful in explaining the relation of quantum orbits to phase matching in HHG. It allowed thus to realise coherence control in high-order harmonics Pascal1995 , and understand the behaviour of the phase of the atomic polarization in high-order harmonic generation Maciek1995 , which in turn allowed the construction of schemes for generation of attosecond pulse trains using HHG Antoine1996 . The first such trains were observed in by P. Agostini et consortes in 2001 Paul2001 . Equally well, the quantum simple man’s model allowed for explanations, both intuitive and quantitatively accurate, of the generation of a single isolated attosecond pulse by an ultrashort, few-cycle laser pulse Hentschel2001 ; Drescher2001 ; Kienberger2004 ; Sansone2006 .

The crowning of the SFA applications in the 1990s was perhaps the Science paper Pascal2001 , in which theory was confronted with the experimental results of the groups of the late B. Carré and P. Salières at Saclay, on quantum orbits in HHG, and of G. Paulus and H. Walther at MPI Garching, on quantum orbits in ATI, driven by elliptically-polarized laser fields.

I.3 The SFA today

Since the formulation of various versions of the SFA, starting from Keldysh in 1964-1965 until the approaches based on simple man’s models, formulated in 1993-1994, the SFA kept being one of the most important theoretical tools of the physics of matter in intense laser fields. There are several important review articles and books that either describe or include description of developments and applications of the SFA Krausz2009 ; Milosevic2006 ; Becker2012 ; Reiss2018 . There were various directions in which the theory of the SFA developed in the recent 20 years, with many authors working on technical issues of improving the accuracy of the theory, while others generalised the theory to novel regions, as we move forward into the future Lindroth2019 . We summarize here some of the main trends in this body of work:

• •

Coulomb corrections. A number of strong-field phenomena, particularly in ionization experiments, show features caused by the ion’s Coulomb potential that evade the SFA, ranging from Coulomb focusing Rudenko2005 and the asymmetric photoelectron spectra produced by elliptical polarizations Bashkansky1988 to the more recent “ionization surprise” of the so-called Low-Energy Structures Blaga2009 . Early work focused on including the Coulomb potential through a Born series, often with a single rescattering used very successfully for the ATI plateau and non-sequential double ionization (NSDI) Milosevic1998 , but this is generally insufficient for Coulomb-dependent phenomena. Current approaches include the use of an oscillating Coulomb-wave basis for the continuum (the Coulomb-Volkov approximation Arbo2008 ), the eikonal inclusion of the Coulomb Hamiltonian to solve the continuum TDSE (the analytical $R$-matrix theory Torlina2012 ) and the direct modification of the SFA’s trajectory language to include Coulomb potential influence on the action and the trajectories (the Coulomb-Corrected SFA Popruzhenko2008 ). In addition to more “phenomenological” approaches (cf. Refs. Yan2010, ; Smirnova2008, ) there have been very elegant approaches based on the Feynman path integral formulation Lai2015 . The most recent results based on the path-integral approach are discussed in Refs. Maxwell2017, ; Milosevic2017, ; Maxwell2018, ; Maxwell2018a, .

• •

Saddle point techniques. The simple man’s model’s classical trajectories are encoded in the SFA as the quantum orbits obtained as the saddle-point contributions to the oscillatory integrals. Understanding the nature of these saddle points in the complex time and momentum planes Kopold2000 ; Kopold2002 ; Morisson2000 allows for a clear understanding of the coherent contribution of each pathway Pascal2001 , and it also paves the way for experiments showing the contributions of other orbits Sansone2006b ; Zair2008 . Technical improvements include the regularization of discontinuities at the cutoff via uniform approximations Morisson2002 ; Milosevic2002 , and the extension to multi-electron configurations Ivanov2014 . On the other hand, some problems, such as the inclusion of field dressing of the ground state PerezHernandez2009 , are less amenable to saddle-point analysis.

• •

Applications to novel systems: Two-electron systems. In the last two decades the SFA has been successfully applied to two-electron systems. Most of these approaches used SFA formulations based on $S$-matrix theory à la W. Becker Kopold2000 ; Becker2012 , while the use of the physics of strong laser fields for imaging goes back to the seminal Refs. Itatani2004, (for HHG) and Huismans2011, (for ATI). In fact, one could argue that two-electron experiments on cold target recoil ion momentum spectroscopy (COLTRIMS) Weber2000 pioneered the imaging methods using strong-field physics. The phenomenon of interest here is NSDI, which occurs in accordance with P. Corkum’s idea of a recollision-driven model Corkum1993 . Still, NSDI has two faces. If the ionization potential of the target ion is smaller than $3.17U_{p}$, the recolliding electron may directly cause stripping of another electron, since there is enough energy for that to occur; this scenario is called Electron Impact ionization (EII). On the other hand, if the recolliding electron does not have enough energy, it may still excite the target ion to an excited state, from which direct tunneling might easily take place; this scenario is known as Recollision Excitation with Subsequent ionization (RESI). In EII, electrons are typically ejected step by step, with most of the quantum interference effects getting washed out, and the standard SFA and quasi-classical trajectory models work very well Figueira2011 ; Becker2012 ; Staudte2007 ; Parker2006 ; Ye2008 ; Emmanouilidou2008 . The early studies of the RESI observed that time delay leads to back-to-back electron ejection, and it was assumed that interferences between different channels (different intermediate excited states, etc.) were irrelevant Eremina2003 ; Haan2010 . Pretty soon, however, a myriad of shapes in the electron momentum distributions were observed in RESI and, moreover, experiments were in clear contrast to the simplest SFA theories Bergues2012 ; Sun2014 ; Kubel2014 ; Ye2010 ; Emmanouilidou2009 . It was then realised that the interference must be accounted for in RESI Shaaran2010 ; Chen2010 ; Hao2014 . A lot of insight was gained by the sophisticated analysis of the saddle-point approach Shaaran2012 . A more complete understanding of the RESI phenomenon, taking into account interference effects, was only achieved recently Maxwell2015 ; Kubel2014 ; Maxwell2016 . All these results allow, in principle, to work backwards toward the experimental data to reconstruct the states of the excited electron involved in RESI. Amazingly, the channel interference in RESI seems to have been observed recently Quan2017 .

• •

Applications to novel systems: Attosecond streaking. The ideas of SFA have also been successfully applied to the single-photon ionization of atoms by high-frequency XUV pulses in the presence of a moderately strong infrared (IR) field, where the IR dresses the continuum and affects the final momentum. This configuration, known as attosecond streaking Itatani2002 ; Yakovlev2005 , is extremely sensitive to the details of both fields and the ionization process. Because of this, it is routinely used as a benchmark for ultrafast light sources, providing direct measurements of optical waveforms Goulielmakis2004 , as well as high-temporal-resolution observations of ionization processes Schultze2010 .

• •

Applications to novel systems: Large molecules. One of the most natural arenas for application of the SFA and related theories is to the strong-field dynamics of molecules. The earliest approaches used the intense-field many-body $S$-matrix theory (IMST) Faisal1999 ; Muth-Bohm2000 ; Muth-Bohm2001 ; Jaron-Becker2003 ; Jaron-Becker2004 later referred to as the molecular SFA Kjeldsen2004 ; Milosevic2006b . This is a Bessel-function expansion of SFA theory, which was used to model ATI ionization rates Faisal1999 ; Muth-Bohm2000 ; Muth-Bohm2001 and photoelectron angular distributions Jaron-Becker2003 ; Jaron-Becker2004 ; Kjeldsen2004 of molecular targets. Further work used saddle-point methods to link back to the semi-classical dynamics, such as the recollision physics analysed in Refs. Lein2002, , Spanner2004, , for HHG Chirila2006a ; Chirila2006b ; Faria2007 , including multi-electron effects Smirnova2009 , photoelectron angular distributions Milosevic2006b and spectra Hetzheim2007 for ATI, as well as NSDI Faria2008 ; for reviews on these methods see Refs. Lein2007, ; Augstein2012, . Two related alternatives to SFA for molecules are the so-called molecular ADK (MO-ADK) Tong2002 and quantitative rescattering (QRS) Le2008 theories. This theoretical work has also been mirrored and enhanced by substantial advances in the experimental control of strong-field processes for molecules, including work on molecular structure Krausz2009 ; Xu2016 ; Meckel2008 ; Blaga2012 ; Wolter2016 ; Burt2018 and dynamics Krausz2009 ; Lein2007 ; Baker2006 ; WenLi2008 ; Worner2010 ; Wolter2016 ; Amini2019 , electronic excitation in molecules in the XUV and X-ray energy domains Villeneuve2018 ; Chatterley2016 , and the ability to follow the dynamics of an attosecond electron wave packet ejected from molecules Itatani2004 ; Krausz2009 .

In a series of recent papers we revisited the SFA model on ATI for few-cycle infrared laser pulses Noslen1 ; Noslen2 ; Noslen3 ; Noslen4 ; Noslen-diss . We compared it first with the numerical solution of the TDSE in one (1D) and two (2D) spatial dimensions for an atomic system Lewenstein1995 . We developed and generalized there an analytical atomic model based on a non-local (short range) potential. In the first paper Noslen1 we analysed ATI for an atom, followed by ATI Noslen2 and HHG Noslen3 for diatomic molecules. Here we paid special attention to non-physical terms which arise in the theory if plane waves are used instead of the exact continuum scattering states of the system. Finally, in Refs. Noslen4, and Noslen-diss, , we generalized our approach to molecular trimers and quadrimers, and attempted to describe Laser-Induced Electron Diffraction for such targets. The ultimate goal of this theory is to characterize the time evolution of the target (its size, configuration, its molecular orbital, its dynamical configuration, etc.) by looking at the ATI spectrum and angular distributions, especially in the region of high energies, corresponding to rescattering processes.

• •

Applications to novel systems: Solid state. In the past 5 years, the strong field community has increasingly turned its attention to condensed matter systems, following the observation of high-order harmonics from bulk crystals subjected to strong laser fields Ghimire2011 ; Ghimire2019 . Our understanding of dynamics in gases, based essentially on the SFA, has recently been extended to yield crucial insights into microscopic attosecond phenomena in condensed matter Vampa2015 ; Vampa2014 ; Osika2017 . The merger of strong field physics with solids has the potential to revolutionize contemporary electronics Krausz2014 , as well as yield crucial fundamental insights into long-standing problems in condensed matter physics. For instance, the first direct measurements of the Berry phase Liu2017 ; Luu2018 were accomplished using HHG. These recent measurements promise to be of great interest to the condensed-matter community, due to the important role played by Berry phase in the anomalous Hall effect and in topological insulators, among other fields.

• •

Applications to novel systems: Atto-nanophysics. In the last decade the SFA has been successfully applied to situations in which HHG, ATI or NSMI come directly from nano-structures, for instance through plasmonic excitations, or from atomic/molecular targets located close to nano-structures. In the latter case, plasmonic enhanced electromagnetic fields close to the structures serve to excite the targets. Recent developments of this area, termed atto-nanophysics, are extensively described in the review Ciappina2017 .

• •

Applications to novel systems: Quantum simulators. The strong-field dynamics described by the SFA have very close analogues in the motion of cold atoms in optical traps Dum1998 , particularly via the Kramers-Henneberger correspondence between a dipole coupling and a ‘shaken’ atomic potential Morales2011 . In the decades since, cold trapped atoms have become one of the primary platforms for quantum simulation LewensteinQuantumSimulators , and several works have explored the possibility of using cold-atom simulators to probe the strong-field dynamics described by SFA Arlinghaus2010 ; Luhman2015 ; Sala2017 ; Rajagopal2017 , thereby allowing a complementary look at observables (like e.g. the instantaneous wavefunction, or a full quantum state tomography on the outgoing particles) that are inaccessible to conventional experiments. While some dedicated experimental efforts to implement this are still at the tool-building stage Ramos2018 , there are already functioning quantum simulation platforms for ultrafast physics Senaratne2018 ; Ramos2019 , which should provide growing opportunities to test SFA physics in new ways. Similar complementary views on strong-field dynamics should also be available via photonic simulation, using the natural Schrödinger-equation correspondence for optical fibers Kahn2017 .

I.4 The present paper

The present article is organized as follows. Section I, Introduction, covers the motivations for the paper (Subsection I.1), and the past (Subsection I.2) and present (Subsection I.3) of the SFA. In Section II we present a short explainer of the basic phenomena and processes: HHG, ATI and NSMI, including representative and explanatory figures. In Section III we review in more detail the theory, which describes the HHG and ATI processes within the version of Strong Field Approximation following Ref. Lewenstein1994, . In particular, we present the derivation of the transition amplitudes for both the direct and the rescattered electrons, as well as for the time-dependent dipole moment. We develop in detail the mathematical foundations towards the final results by starting from the Hamiltonian, which describes the atomic system and the TDSE associated to it. Section IV is devoted entirely to two-electron processes. We derive the SFA for this case in an explicit time-dependent approach, using exact eigenfunctions for the states in the one- and two-electron continua and perturbation theory with respect to the “less-singular” parts of the continuum-continuum matrix elements. We then analyze here the interplay between the EII and RESI processes.

In Section V, we formulate our theory for the case when the single active electron (SAE) approximation is applied to a molecule undergoing temporal evolution of its nuclear configuration. This is done using the Born-Oppenheimer approximation and classical equations of motion for the nuclei. We consider first the simplest case in which the molecule’s dynamics (vibrations, dissociations) do not affect the SAE electronic dynamics. Even in this simple case, novel effects appear in the SFA dynamics, such as the appearance of a temporal Berry phase. In Subsection V.2 we go beyond this approximation and consider the self-consistent dynamics, in which SAE dynamics affect nuclei and vice versa. We conclude in Subsection V.3 where we present an outlook on extending our quasi-analytical model to more complex atomic and molecular systems. Finally, in Section VI we briefly review the recent application of SFA theory to the generation of harmonics in solid-state systems.

Appendix A discusses time dependent ADK rates. In Appendix B we sketch calculations of $a(t)$—the amplitude of the ground state within our generalized SFA theory. In Appendix C we introduce the model for our atomic and molecular systems that uses a particular form of a non-local short-range separable potential. The matrix elements to describe the direct ionization and the re-scattering processes are then computed analytically.

We then offer additional material regarding the two-electron theory: in Appendix D we discuss the properties of the dipole matrix elements involved, and in Appendix E we pose full forms for the two-electron integro-differential equations derived from the TDSE. In Appendix F we analyze a single-electron toy model for HHG from a quenched molecule. Finally, in Appendix G, we present solutions of the RESI and EII equations using additional approximations for the dipole matrix elements, neglecting electron-electron interaction effects for those elements.

II Strong-field processes in an atomic gas

Over the lifetime of strong-field physics, SFA theory has accounted for a broad variety of physical phenomena which could not be explained by traditional perturbation theory. These phenomena involve light-matter interaction using lasers whose field strengths are comparable to the Coulomb force of attraction between electrons and protons. Consequently, this can lead to the distortion of the Coulomb potential and the subsequent lowering of the barrier to ionization in the strong-field regime. Many nonlinear ionization processes can be initiated in this regime, such as multi-photon ionization (MPI), above-threshold ionization (ATI), tunnel ionization (TI), and over-the-barrier (OTB) ionization. These processes are shown in Fig. 1 with their operating conditions summarized in Table 1. The ionization regime of operation can be identified by the Keldysh parameter, $\gamma$, given by

[TABLE]

where $I_{p}$ is the ionization potential (i.e. the energy required to eject an electron from the ground state to the ionization continuum), and $U_{p}$ is the ponderomotive energy (i.e. the average kinetic energy of the oscillations of a free electron in a laser field) given by

[TABLE]

where $e$ is the elementary charge, $\mathcal{E}_{0}$ is the electric field amplitude, $m_{\rm e}$ is the mass of an electron, $\omega_{0}$, $\lambda_{0}$ and $I_{0}$ are the central laser frequency, wavelength and intensity, respectively, $\epsilon_{0}$ is the vacuum permittivity, and $c$ is the speed of light. We give, in the rest of this section, a brief summary of strong-field processes that are based on ATI and TI.

II.1 Above-Threshold Ionization (ATI)

Above-threshold ionization (ATI) is an extension of multi-photon ionization where multiple photons are absorbed to not only access the ionization continuum but to surpass the $I_{p}$ by more than one photon, $\hbar\omega$ Agostini1979 . In the initial studies of ATI there was substantial interest in determining the nature of this process, and discerning the relative contributions of perturbative multi-photon ionization and tunnel ionization. Nowadays, ATI is considered holistically as a process where non-perturbative phenomena coexist with perturbative ones (see, for instance, Ref. Paulus1994, ).

In a typical ATI photoelectron spectrum, as shown in Fig. 2, a series of peaks are observed that correspond to each photon absorbed above the $I_{p}$, each of which is separated by a single photon energy, $\hbar\omega$. More strongly, ATI can be observed in the high energy range of the photoelectron spectrum ($2U_{p}\leq E_{r}\leq 10U_{p}$), referred to as high-order ATI (HATI), where recollision-based strong-field physics can appear, giving rise to elastic and inelastic scattering.

II.2 High-Harmonic Generation (HHG)

Attosecond laser pulses of high-photon energies in the extreme ultraviolet (XUV) and X-ray energy region can be produced by high-harmonic generation Krause1992 ; Corkum1993 ; Krausz2009 . Pulse trains of attosecond radiation are generated using a multi-cycle femtosecond driving laser pulse, as presented in Fig. 3a with the $27^{\rm th}-85^{\rm th}$ harmonics shown Jin2011 . Similarly, a single attosecond pulse with a broadband spectrum can be generated using a near-single cycle driving laser pulse with a continuous broadband spectrum, as shown in Fig. 3b Buades2018 , though other so-called ‘gating’ schemes are also possible CorkumBurnett1994 ; Kovacev2003 ; Mashiko2008 ; Pfeifer2007 ; Vincenti2012 . The maximum HHG cut-off energy, $E_{\rm max}$, that can be generated is given by the simple-man’s-model-like $E_{\rm max}=3.17U_{p}+I_{p}$, and is observed as the abrupt end to the HHG plateau.

These microscopic aspects aside, it is also important to remark that HHG is a macroscopic nonlinear optical process that requires the coherent combination of a large number of emitters to be observed experimentally, and this requires that specific attention be paid to the phase-matching conditions Heyl2016 , which are often the determining limitation in the production of harmonics.

II.3 Inelastic Scattering: Non-Sequential Double Ionization (NSDI)

If the recollision of the returning electron ($e_{1}$) with the parent ion is inelastic, then it can transfer enough energy to eject a second electron ($e_{2}$). This process is known as non-sequential double ionization (NSDI) Walker1994 , and it can proceed through two ionization pathways upon the recollision of $e_{1}$ Pullen2017 , as shown in Fig. 4a: (i) direct ionization of $e_{2}$, known as electron-impact ionization (EII) Corkum1993 ; or (ii) impact excitation of $e_{2}$ subsequently followed by its delayed tunnel ionization, known as recollision-excitation with subsequent ionization (RESI) Feuerstein2001 . A typical signature of NSDI, particularly in the EII regime, is the correlated detection of two electrons ($e_{1}$ and $e_{2}$) in the same emission direction within the two-dimensional momentum map $(p_{\parallel,1},p_{\parallel,2})$ of the longitudinal momenta of the two electrons, as shown in Fig. 4b.

II.4 Elastic Scattering: Laser-Induced Electron Diffraction (LIED)

The highly-energetic returning electron can collide elastically and scatter on the target ion, leading to a momentum transfer between the electron and parent ion. This is known as laser-induced electron diffraction (LIED) Meckel2008 ; Blaga2012 ; Wolter2016 ; Pullen2016 ; Ito2016 ; Amini2019 and it can be explained in the framework of laser-driven electron recollision Krause1992 ; Corkum1993 . Structural information is embedded in the photoelectron momentum distribution, appearing as oscillations in the high-energy part of spectrum corresponding to recollision-based physics ($2U_{p}\leq E_{r}\leq 10U_{p}$) as a function of the emission angle, as shown in Fig. 5 with a zoomed-in view of these oscillations given in the inset. It should be noted that, as opposed to HHG, phase-matching is irrelevant for inelastic and elastic scattering processes (i.e. NSDI and LIED), since the macroscopic observable is an incoherent combination of the emission from the different atoms in the laser focus.

III Strong field approximation for single-electron processes

Strictly speaking, neglecting nuclear motion, an atomic or molecular system interacting with a strong electric field pulse is described by the time-dependent Schrödinger equation (TDSE) that captures both the evolution of the (electronic) wave function and the time evolution of the physical observables. The numerical solution of the TDSE offers a full quantum mechanical description of the laser-matter interaction processes; it has been used extensively to study HHG KulanderGaarde1998 ; TateDiMauro2007 ; Marcelo2012 ; Scrinzi2014 and ATI Kulander1987 ; Muller1999 ; Bauer2005classical ; Blaga2009 ; RufKeitel2009 in atomic and molecular systems. However, the full numerical integration of the TDSE in all the degrees of freedom of the system is computationally very demanding, when it is at all possible. Moreover, a physical interpretation of the numerical results is highly nontrivial, as always for an ab initio technique. Within this framework, then, approximate methods are welcome, and the SFA has consistently been shown over the years to be the workhorse tool for that role.

III.1 Hamiltonian and TDSE

Let us consider an atom or molecule under the influence of an intense laser field in the so-called single active electron (SAE) approximation. In the limit when the wavelength of the laser $\lambda_{0}$ is large compared with the Bohr radius, $a_{0}$ ($5.29\text{\times}{10}^{-11}\text{\,}\mathrm{m}$), the electric field of the laser beam around the interaction region can be considered spatially homogeneous. Consequently, the interacting atoms will not experience the spatial dependence of the laser electric field and, hence, only its time variation is taken into account—this is the so-called dipole approximation. Note, on the other hand, that certain dynamical effects, even in the long-wavelength limit, can break this approximation Reiss2018 . Within this framework, the laser electric field can be written as:

[TABLE]

The field of Eq. (3) has a carrier frequency $\omega_{0}=\frac{2\pi c}{\lambda_{0}}$, where $c$ is the speed of light, and a peak amplitude $\mathcal{E}_{0}$. We consider here that the laser field is linearly polarized along the $z$ direction, with a pulse envelope ${f(t)}$ and a carrier-envelope phase $\phi_{0}$. More generally, we could consider time-dependent polarization, i.e. replace $\mathcal{E}_{0}{\bf e}_{z}\to{\bf E}_{0}(t)$.

The TDSE reads:

[TABLE]

where the Hamiltonian, $\hat{H}$, describes the laser-target system in SAE approximation, and is the sum of two terms, i.e.

[TABLE]

where $\hat{H}_{0}$ is the laser-free Hamiltonian of the atomic or molecular system

[TABLE]

with $V(\hat{\textbf{r}})$ the effective SAE atomic or molecular potential, $m$ the electron mass, and $\hat{U}=-e{\bf E}(t)\cdot\hat{\bf r}$ the dipole coupling, which describes the interaction of the atomic or molecular system with the laser radiation, written in the length gauge Bauer2005 ; Galstyan2016 and under the dipole approximation. Note that in atomic units, the electron charge, denoted by $e$, is $e=-1\>$a.u., and the Planck constant and electron mass are both set to unity, $\hbar=m=1\>$a.u. In this work, however, we keep the explicit constants.

III.2 “Standard” SFA à la Lewenstein

We shall restrict ourselves to the regime of low laser frequency and relatively high intensity, where the SFA is expected to be valid Keldysh1965 ; Faisal1973 ; Reiss1980 ; Lewenstein1994 ; Lewenstein1995 ; Grochmalicki1986 and to describe well the laser-matter interactions. This corresponds to the tunnelling regime, where the Keldysh parameter $\gamma$ is less than one, $\gamma<1$. In this regime the effects of atomic effective potential on the dynamics of electrons in the continuum are assumed to be small, and they can be treated using perturbation theory. These observations suggest to formulate the “standard SFA” as follows:

1. (i)

The strong field laser does not couple with any bound state beyond the ground state, $|0\rangle$, so that only it and the continuum (scattering) states, $|\textbf{p}\rangle$, are taken into account in the dynamics; 2. (ii)

The amplitude of the ground state, $a(t)$, is considered to be known. 3. (iii)

The continuum states are taken from the basis of exact scattering states, which are eigenstates via

[TABLE]

of the atomic Hamiltonian with a fixed outgoing (kinetic) momentum ${\bf p}$. The continuum-continuum matrix element from ${\bf p}$ to ${\bf p}^{\prime}$ are then decomposed into their most singular part, proportional to $i\hbar{\nabla}_{\bf p}\delta({\bf p}-{\bf p}^{\prime})$, and the “rest” Grochmalicki1986 ; Lewenstein1994 ; Lewenstein1995 . The “rest” is then treated in a perturbative manner Lewenstein1995 .

The following comments are necessary in order to specify more precisely the above points.

1. Ad (i)

Based on the statement (i), the electronic state $|\Psi(t)\rangle$ that represents the time evolution of the system is a coherent superposition of the ground $|0\rangle$ and the continuum $|\textbf{p}\rangle$ states Lewenstein1994 ; Lewenstein1995 :

[TABLE]

The factor $a(t)$, representing the amplitude of the ground state, is assumed to be known (see below for the ways to evaluate or estimate it). The prefactor $e^{\textit{i}I_{p}\textit{t}/\hbar}$ represents the phase oscillations which describe the accumulated electron energy in the ground state ($I_{p}=-E_{0}$ is the ionization potential, with $E_{0}$ the ground-state energy of the target system). Furthermore, the transition amplitude to the continuum states is denoted by $\textit{b}(\textbf{p},t)$, and it depends both on the kinetic momentum of the outgoing electron and the laser pulse. Note that, if needed, other (relevant) bound states may be taken into account in the expression (8) (cf. Refs. Ivanov2014, ; PerezHernandez2009, ; Sanpera1996, ). 2. Ad (ii)

There are several ways of evaluating or estimating $a(t)$, depending on the regime of parameters.

• •

First, one can use ab initio TDSE of the target system to determine the amplitude $a(t)$. This is obviously quite costly numerically, but it is much less costly than a full solution of the TDSE which is also required to calculate photoelectron momentum spectra or angular distributions which would need much higher precision, memory and disk storage, and higher computation times.

• •

Second, one can use any “cheap” approximate method to calculate $a(t)$, such as phase-space averaging or the truncated Wigner approximation Berman2018 .

• •

Third, a broadly-used method is to calculate $a(t)$ analytically using the ionization rates according to the Ammosov-Delone-Krainov theory (ADK rates Ammosov1986 ). To this end, one generalizes these rates to depend on time locally through the time dependence of the laser electric field (also known as the quasi-static approximation), which is generally a rather straightforward task (see Appendix A). This approach is valid in the quasi-static regime, when not only the laser frequency, but also the rate of change of the pulse envelope function $f(t)$ are small—meaning that the laser pulse is longer, so that it includes several optical periods.

• •

Fourth, when the pulse is very short, or it is long but not too strong, there is practically no depletion of the ground state, i.e. $a(t)\simeq 1$. This happens, for instance, for moderately long pulses when the ponderomotive energy is lower than the saturation energy of the system $(U_{p}<U_{\mathrm{sat}})$.

• •

Fifth, one can calculate $a(t)$ within our SFA self-consistently. This approach was already discussed in Ref. Lewenstein1994, , but it turned out not to be very precise for the longer pulses—the ADK rates were giving much better agreement with the exact solutions of the TDSE and with the experimental data. This approach seems to be, however, much more adequate and precise for ultrashort, few-cycle pulses. We describe it in detail in Appendix B. 3. Ad (iii)

The continuum-continuum matrix element, independently of the fact whether the effective SAE potential is short-range (as it is for model atoms and negative ions) or Coulomb-like, has the general form:

[TABLE]

where the part $\hbar{\bf g}({\bf p},{\bf p}^{\prime})$ is less singular—typically the strongest singularity it contains corresponds to the on-energy-shell gradient of the Dirac delta of ${\bf p}^{2}-({\bf p}^{\prime})^{2}$. This part is responsible for rescattering effects in ATI and recollisions in NSMI. Note that since we insist on using the exact scattering states, the dipole matrix element $\langle{\bf p}|\hat{\bf r}|0\rangle$ (together with the rescattering continuum-continuum matrix elements) does include the full effects of the effective SAE potential, comprising both the short-range effects as well as any long-range Coulomb effects (if present).

Note also that the SFA in the present formulation (actually equivalent to that of Ref. Lewenstein1995, ) does not involve plane waves or Volkov solutions! The majority of authors, including ourselves, “erroneously” (in the view of the present formulation) claim that SFA corresponds to the use of Volkov states in the continuum. This is, in principle, false and dangerous. One can use additional approximations, and approximate the exact scattering states by plane waves in the calculations, but this is an additional approximation! It does simplify life and allows one to calculate many things more easily, but it also leads to problems, especially in the case of molecules and other extended targets.

These problems are due to the fact that a plane wave $|{\bf p}\rangle$ is not orthogonal to the ground state $|0\rangle$, so that the matrix element $\langle{\bf p}|{\bf R}_{0}|0\rangle\neq 0$, where ${\bf R}_{0}$ is the typical internuclear distance, which is just a constant vector. The lack of orthogonality of $\langle{\bf p}|0\rangle\neq 0$ leads to various non-physical and misleading results in applications of, say, “primitive” SFA to molecules (for remedies see Refs. Noslen2, ; Noslen3, ). We stress: no remedies are needed, on the other hand, if the exact scattering states are used, since then the orthogonality is assured by construction.

Why, then, do the plane waves and Volkov solutions appear at all? Clearly, this is due to the fact that in the zeroth approximation of SFA we neglect the contribution of $\hbar{\bf g}({\bf p},{\bf p}^{\prime})$. In this case, the full continuum-continuum matrix element becomes $\langle{\bf p}|\hat{\bf r}|{\bf p}^{\prime}\rangle=i\hbar{\nabla}_{\bf p}\delta({\bf p}-{\bf p}^{\prime})$, and is exactly equal to that obtained for plane waves and Volkov solutions. That means that the quasi-classical action, describing the propagation of electrons in the continuum, does indeed have a free electron form. For short-range effective potentials this is acceptable, but not for the Coulomb-like ones. That is why the so-called Coulomb corrections are easily included in $\langle{\bf p}|\hat{\bf r}|0\rangle$ or $\hbar{\bf g}({\bf p},{\bf p}^{\prime})$, but much effort has been devoted to find Coulomb corrections to the action—see the Introduction for the relevant references.

III.3 Solutions of the SFA equations

Our main task in this subsection will be to derive a general expression for the amplitude ${b}({\bf p},t)$, which then will be used to calculate ATI spectra and angular distributions, as well as HHG spectra. After some algebra, the time variation of the ground state amplitude, $a(t)$, and the transition amplitude $\textit{b}(\textbf{p},t)$ read:

[TABLE]

The first term on the right-hand side of Eq. (11) represents the free phase evolution of the electron in the absence if the oscillating laser field. In the second term we have defined the bound-free transition dipole matrix element as

[TABLE]

Finally, the last two terms describe the continuum-continuum transition, $\nabla_{\bf p}b({\bf p},t)$, without the influence of the scattering center, and by considering the core potential, $\int{\textit{d}^{3}e\textbf{p}^{\prime}\>\textit{b}(\textbf{p}^{\prime},t){\bf g}({\bf p},{\bf p}^{\prime})}$. Here ${\bf g}({\bf p},{\bf p}^{\prime})$ denotes the rescattering transition matrix element, where the potential core plays an essential role:

[TABLE]

Note that (11) is a linear integro-diffential equation for $\textit{b}(\textbf{p},t)$. In the following, we shall describe how it is possible to compute the transition amplitude, $b({\bf p},t)$, by applying zeroth- and first-order perturbation theory to the solution of the partial differential equation (11). We will split the solution of the transition amplitude, $b({\bf p},t)$, into two parts: $b_{0}({\bf p},t)$ and $b_{1}({\bf p},t)$, i.e. $b({\bf p},t)=b_{0}({\bf p},t)+b_{1}({\bf p},t)$. The zeroth order of our perturbation theory $b_{0}({\bf p},t)$ will be called the direct term. It describes the transition amplitude for a laser-ionized electron that will never rescatter with the remaining ion. On the other hand, the first-order term, which we call the rescattered term, $b_{1}({\bf p},t)$, refers to the electrons that, once ionized, will have a certain probability of rescattering with the potential of the parent ion.

III.4 Direct-ionization amplitude

Let us consider the process in which the electron is ionized and does not return to its parent ion. This process is modelled by the direct photoelectron transition amplitude $b_{0}({\bf p},t)$. As the direct ionization process should have a larger probability compared with the rescattering one Lewenstein1995 , one can neglect the last term in Eq. (11). This is what we refer to as the zeroth order solution:

[TABLE]

The above equation is a first-order inhomogeneous differential equation, which is easily solved by conventional integration methods (see e.g. Ref. Lev, ). Therefore, the solution can be written as

[TABLE]

Here, we have considered that the electron appears in the continuum with kinetic momentum ${\bf p}(t^{\prime})={\bf p}+e{\bf A}(t)/c-e{\bf A}(t^{\prime})/c$ at the time $t^{\prime}$, where p is the final kinetic momentum, and $\textbf{A}(t)=-c\int^{t}{\textbf{E}(t^{\prime})dt^{\prime}}$ is the vector potential of the electromagnetic field, with $c$ the speed of light. In particular, the vector potential at the time when the electron appears at the continuum $t^{\prime}$ is ${\bf A}(t^{\prime})$, and at a certain detection time $t$, the vector potential reads ${\bf A}(t)$. In addition, it is possible to write Eq. (15) as a function of the canonical momentum $\textbf{p}_{c}$, defined by $\textbf{p}_{c}={\bf p}+e{\bf A}(t)/c$, and therefore the probability transition amplitude for the direct electrons simplifies to that from Ref. Lewenstein1994, , where we have eliminated the subscript ${c}$:

[TABLE]

This expression is understood as the sum of all the ionization events which occur from the time $t^{\prime}$ to $t$ Pascal2001 . Then, the instantaneous transition probability amplitude of an electron at a time $t^{\prime}$, at which it appears into the continuum with momentum ${\bf p}(t^{\prime})=\textbf{p}-e\textbf{A}(t^{\prime})/c$, is defined by the argument of the integral in Eq. (16). Furthermore, the exponent phase factor in Eq. (16) denotes the “semi-classical action”, ${S}({\bf p},t,t^{\prime})$, that defines a possible electron trajectory from the birth time $t^{\prime}$ until the “detection” time $t$ Lewenstein1995 :

[TABLE]

As our purpose is to obtain the final transition amplitude $b_{0}({\bf p},t)$, the time $t$ will be fixed at the end of the laser field, $t=t_{\rm F}$. For our calculations, we thus define the integration time window as $t\in[0,t_{\rm F}]$. Therefore, we set ${\bf E}(0)={\bf E}(t_{\rm F})={\bf 0}$, in such a way to make sure that the electromagnetic field is a time oscillating wave and does not have static components. The same arguments are applied to the vector potential ${\bf A}(t)$. In concrete calculations we have defined the laser pulse envelope as $f(t)=\sin^{2}(\frac{\omega_{0}t}{2N_{c}})$ where $N_{c}$ denotes the number of total cycles.

Note that, for an arbitrary electromagnetic field, it is possible for ${\bf A}_{\rm F}(t)\neq 0$, i.e. the vector potential does not necessarily vanish at the end of the pulse, implying that the kinetic momentum at $t_{\rm F}$ is ${\bf p}_{\mathrm{kin}}={\bf p}-e{\bf A}(t_{\rm F})/c$; if that is the case then it should be considered carefully, since it is ${\bf p}_{\mathrm{kin}}$ which is detected in experiments. However, for laser pulses that are focused away from their source and in the paraxial approximation, nonzero-area pulses of this form are not possible, and the vector potential can be taken as zero on both sides of the pulse.

III.5 Rescattering transition amplitude

In order to find a solution for the transition amplitude of the rescattered photoelectrons, $b_{1}({\bf p},t)$, we have considered the rescattering core matrix element $\textbf{g}(\textbf{p},\textbf{p}^{\prime})$ term of Eq. (11) different than zero, i.e. $\textbf{g}(\textbf{p},\textbf{p}^{\prime})\not=\textbf{0}$. In addition, the first-order perturbation theory is applied to obtain $b_{1}({\bf p},t)$ by inserting the zeroth-order solution $b_{0}({\bf p},t)$ in the right-hand side of Eq. (11). Then, we obtain $b_{1}({\bf p},t)$ as a function of the canonical momentum ${\bf p}$ (neglecting the subscript $c$) as follows:

[TABLE]

This last equation contains all the information about the rescattering process. In particular, it refers to the probability amplitude of an emitted electron at the time $t^{\prime\prime}$, with an amplitude given by

[TABLE]

In this step the electron has a kinetic momentum of ${\bf v}^{\prime}(t^{\prime\prime})=\textbf{p}^{\prime}-e\textbf{A}(t^{{\prime}{\prime}})/c$. The last factor, $\exp{\left[-\textit{i}S({\bf p}^{\prime},t{{}^{\prime}},t^{\prime\prime})\right]}$, is the accumulated phase of an electron born at the time $t^{\prime\prime}$ until it rescatters at time $t^{\prime}$. The intervening term, ${\bf g}({\bf p}-e{\bf A}(t^{\prime})/c,{\bf p}^{\prime}-e{\bf A}(t^{\prime})/c)$, contains the structural matrix element of the continuum-continuum transition at the re-scattering time $t^{\prime}$. At this particular moment in time, the electron changes its kinetic momentum from ${\bf p}^{\prime}-e{\bf A}(t^{\prime})/c$ to ${\bf p}-e{\bf A}(t^{\prime})/c$. We stress, however, that the term ${\bf g}({\bf v},{\bf v}^{\prime})$ does not necessarily imply that the electron returns to the ion core.

In addition to this, the phase term $\exp\left[-\textit{i}S({\bf p},t,t^{\prime})\right]$ defines the accumulated phase of the electron after the rescattering from the time $t^{\prime}$ to the “final” one $t$ when the electron is “measured” at the detector with momentum p. In particular, note that the photoelectron spectrum, $|b({\bf p},t_{\rm F})|^{2}$, is a coherent superposition of both solutions, $b_{0}({\bf p},t_{\rm F})$ and $b_{1}({\bf p},t_{\rm F})$, together with an interference term:

[TABLE]

So far we have formulated a model, which describes the photoionization process leading to two main terms, namely, a direct $b_{0}({\bf p},t_{\rm F})$ and a rescattering $b_{1}({\bf p},t_{\rm F})$ one. As the complex transition amplitude, Eq. (16), is a single time integral, it can be integrated numerically without major problems. However, the multiple time (“2D”) and momentum (“3D”) integrals of the re-scattering term, Eq. (18), present an increasingly difficult and demanding task from a computational perspective. In order to reduce the computational difficulties, and to obtain a physical meaning of the ATI process, one may employ saddle-point methods to evaluate these highly-oscillatory integrals (see Section III.7 below for details).

The main challenge to calculate the ATI spectrum is then the computation of the bound-free transition dipole matrix element, ${\bf d}({\bf p})$, and the continuum-continuum transition re-scattering matrix element ${\bf g}({\bf p},{\bf p}^{\prime})$ for a given atomic or molecular system. In the Appendices, we discuss how to do this analytically for a model atom or molecule with a short-range separable potential.

III.6 Time-dependent dipole moment

Finally, to analyse the HHG we need to know the electron acceleration, or at least the time dependent electron dipole moment. This is dominantly given by the zeroth order solution of the SFA equations. It is then given by the dynamical version of the celebrated Landau-Dykhne formula,

[TABLE]

which is then generally compared to experiment via its frequency-domain version, the Fourier transform

[TABLE]

where $\Omega$ is the frequency of the emitted harmonic.

III.7 Saddle-point methods and quantum orbits

The SFA results as we have obtained them thus far, i.e. Eqs. (16) and (18) for the direct- and rescattered-electron momentum wavefunctions and Eq. (III.6) for the frequency-domain harmonic dipole, in what is known as their time-integrated versions with explicit integrals over the times of ionization and recollision or recombination, are often perfectly sufficient for the evaluation of the relevant physical observables via a direct numerical integration. However, they generally involve the integration of highly oscillatory terms, such as those contained in the phase factor of the harmonic dipole in (III.6),

[TABLE]

where the phase of the exponential can vary wildly, introducing extreme cancellations in the integrand that require increased precision in the numerical integration to calculate correctly. Generally speaking, the phase in this factor can be estimated by considering its scaling once the pulse amplitude $\mathcal{E}_{0}$ and frequency $\omega$ are factored out by de-dimensionalizing the integral in the exponent, giving contributions which scale with the so-called strong-field parameter $z=U_{p}/\omega$ and with $I_{p}/\omega$. For experiments with a strong low-frequency field, both of these parameters are large, and the exponent in Eq. (23) will quickly cover many radians without giving the rest of the integrand time to change, giving rise to cancellations in the integral.

As mentioned earlier, this problem can be overcome by approximating any relevant oscillatory integrals using the method of steepest descents, which approximates the integrals using the values of the integrand at stationary points of the action – in exactly the same way as the classical trajectories emerge as the stationary-action points of the Feynman path integral Pascal2001 . Using the paradigmatic case of HHG as an example, this requires us to find solutions to the stationary-point equations over the three integration variables,

[TABLE]

which are often termed the tunnelling, recollision, and return equations, respectively. (For other phenomena, these should be adjusted accordingly, by e.g. dropping (25, 26) for the direct-electron ionization amplitude.) The tunnelling equation here, Eq. (24), is the central, determining structure, both because of its prevalence over all SFA applications and because both of the terms on its left-hand side, $\frac{1}{2}(\textbf{p}-e\textbf{A}(t^{\prime})/c)^{2}$ and $I_{p}$, are ostensibly positive (for real $t^{\prime}$), which means that solutions will only be possible if $t^{\prime}$ (and, with it, all the other variables) are complex-valued.

Within that steepest-descent approximation, then, SFA amplitudes are given by a sum over all the relevant saddle-point roots that contribute to the deformed integration contour,

[TABLE]

with an additional Hessian factor $H(t_{s},t^{\prime}_{s},\textbf{p}_{s})$ that accounts for the width of the complex-integration Gaussians being approximated Bleistein1975 . (Alternatively, it is also possible to perform a partial saddle-point approximation over momentum only, keeping the unique root of the return equation (26) $\textbf{p}_{s}=\textbf{p}_{s}(t,t^{\prime})$ as a function of the ionization and recollision times, and then integrate numerically.) For the full saddle-point method, the ionization time typically has a large imaginary part and it is confined to a small window shortly after the peak of the field, while the recollision time comes in a series of so-called quantum orbits that span the following periods, as shown in Fig. 6.

Typically, the only quantum orbits that contribute significantly to harmonic generation are the so-called short and long trajectories, shown in black and blue (resp.) in Fig. 6(b). The long trajectories, which ionize closer to the peak of the field, have a higher single-atom harmonic yield, but the phase-matching conditions are typically chosen to select the contribution of the short trajectories, which are easier to phase-match. The higher-order returns (shown in red, pink, green and purple), which recollide more than one period after the ionization, typically spend too much time in the continuum accumulating an intensity-dependent phase for them to form a macroscopically-coherent emission, but under dedicated circumstances it is still possible to observe signatures of their presence Sansone2006b ; Zair2008 .

III.8 Polarization effects

The above analysis of quantum orbits becomes, obviously, more complex when the laser fields have more complicated patterns of polarization beyond the simple linear one. Elliptical polarization was considered in the context of HHG already in the pioneering papers by P. Corkum Corkum1993 : the electron trajectories in such situation form ellipses, and essentially miss the parent ion, leading to a rapid decrease of the HHG efficiency with increasing ellipticity. These trajectory-based predictions were first confirmed in experiments by Budil et al. Budil1993 , and they can be used to produce ‘gating’ schemes CorkumBurnett1994 ; Kovacev2003 ; Mashiko2008 to produce isolated attosecond pulses by using a time-dependent polarization that changes across the pulse from circular to linear and back.

The late Bertrand Carré was also among the pioneers of polarization studies in HHG. The first experimental results of ellipticity dependence of the harmonic yield were published in Ref. Antoine1996elliptical, , in which the detailed SFA theory of HHG by an elliptically polarized laser field was investigated. The following seminal paper Antoine1997 was the first one to discuss measurements of the (partial) polarization of high harmonics generated by elliptically polarized laser fields, with careful comparison to SFA-based theory including propagation. This paper stimulated many researchers to search for ways to control polarization of harmonics. The Holy Grail was to generate high-order harmonics with left- and right-circular polarization to be able to use them to study circular dichroism in absorption—to distinguish, for instance, chiral molecules—, or to study chiral effects in magnetism.

The rapid decrease of HHG efficiency with ellipticity suggested looking for scenarios based around a linearly-polarized IR driver. Pioneering ideas were formulated by P.-M. Paul in his doctoral thesis Paul2001 , and developed further in the group of B. Carré by Y. Mairesse, first employing two-photon absorption of one XUV and one IR photon, and later resonant HHG Ferre2015a ; Ferre2015b , as well as HHG generated by linearly-polarized light pulses applied to aligned molecules Levesque2009 ; Zhou2012 ; Shafir2012 .

A breakthrough method was proposed by D. B. Milošević Milosevic2000 and later implemented by O. Cohen, using two circularly-polarized beams with a frequency ratio of $1{:}2$ and opposite helicity Fleischer2014 , which permits the generation of bright phase-matched circularly-polarized extreme ultraviolet high harmonics [Kfir2015, ; for a review see PisantyThesis2016, ]. The original method and results of Ref. Antoine1997, was developed further to completely characterize the state of elliptically-polarized light by electron-ion vector correlations Veyrinas2013 , and finally to realise the complete polarimetry of high harmonics Veyrinas2016 . These methods have recently been applied to HHG generated under O. Cohen’s ‘bicircular’ fields Barreau2018 , providing a clear evidence for depolarization of high harmonics.

It is worth mentioning that more laser fields with “exotic” polarization (spin) and orbital angular momentum have been proposed recently (see Ref. Pisanty2018conservation, and references therein). These so-called polarization torus knots, proposed in Ref. Pisanty2018knotting, , when applied to atoms in a form of ultrashort and ultraintense pulse, generate “exotic” harmonics that conserve torus-knot angular momentum, a topologically-nontrivial mixture between spin and orbital angular momentum.

One should also say that ellipticity and polarization effects play an important role in ATI and multielectron ionization. A nice example of these ATI results is included in the Science paper of 2001 Pascal2001 , where G. Paulus was able to characterize a whole plethora of trajectories corresponding to rescattering of electrons in elliptically polarized laser fields. In the same paper, Carré and Salières present spectra of high harmonics that allow one to identify directly the contribution of the “short” and “long” electronic trajectories.

IV Two-Electron Processes

IV.1 Hamiltonian and states

In order to describe higher order processes we must extend this formalism to include more active electrons. We start by formulating the SFA equations for two active electrons, which allows us to model higher-order strong field ionization process such as non-sequential double ionization (NSDI). We also include excited states in the wave function to allow for the recollision excitation with subsequent ionization (RESI) pathway of NSDI; otherwise only the direct electron impact ionization (EII) pathway would be present. Following a similar procedure to the one electron case, the ansatz for the wavefunction can be written as

[TABLE]

where $\ket{0}$ is the two-electron ground state, $\ket{\mathbf{p},0}$ gives the two-electron state, where the first electron has been promoted to a continuum scattering state with momentum $\mathbf{p}$, $\ket{\mathbf{p},\eta}$ is similar but the second electron is in an excited state with a principal quantum number $\eta$, while $\ket{\mathbf{p},\mathbf{p}^{\prime}}$ denotes both electrons in continuum scattering states. The state $\ket{\eta,0}$ could be included to allow for some additional effects such as single electron frustrated tunnelling, but we will neglect it for now as we are interested in two-electron effects, where this state will play almost no role. These states are all eigenstates of the two-particle Hamiltonian

[TABLE]

where $\hat{\mathbf{p}}_{i}$ are single particle momentum operators, $V(\mathbf{r}_{i})$ gives the interaction of each particle with the atomic or molecular core and $V(\mathbf{r}_{1}-\mathbf{r}_{2})$ gives the interaction between the two electrons. Note that including the interaction between electrons means that none of the two-particle states introduced above can be written as products of one-particle states, e.g. $\ket{\mathbf{p},0}\neq\ket{\mathbf{p}}\ket{0}$ . However, we can write the energy eigenvalue equations for each sector:

[TABLE]

Given that we are accounting for electron correlation, $E_{0}$ will generally be different from $2E_{10}$; however, often the correlation is weak and then this is a good approximation to make. Note that here the eigenstates $|\mathbf{p},0\rangle$, $|\mathbf{p},\eta\rangle$ and $|\mathbf{p},\mathbf{p}^{\prime}\rangle$ do not denote plane-wave states in the continuum, but instead the full one- and two-electron scattering eigenstates with asymptotic outgoing momenta $\mathbf{p}$ and $\mathbf{p}^{\prime}$. This explains the particular form of the eigenenergies.

If one assumes non-interacting electrons, this amounts to dropping the last term in Eq. (29), then the following substitution can be made for each of the two particle states

[TABLE]

Note we do not use anti-symmetric superpositions as we consider two electrons from a singlet spin state, so that the spins will already be anti-symmetric. These will be eigenstates of the Hamiltonian $\hat{H}_{0}$ without the electron-electron interaction term.

As in the one-particle case the full Hamiltonian is given by laser-free and dipole coupling Hamiltonians as in Eq. (5) but now in the two particle case the dipole coupling is given by

[TABLE]

We will proceed as before and derive the integro-differential equations for $a(t)$, $b(\mathbf{p},t)$, $c(\mathbf{p},\eta,t)$ and $d(\mathbf{p},\mathbf{p}^{\prime},t)$. However, first we will introduce the dipole matrix elements required for each possible kind of transition between the two particle states. The matrix elements will follow the convention that the left-hand state will have a lower or equal energy to the right-hand state. Then the dipole matrix elements can be defined in the following way

[TABLE]

where the subindices refer to the two electrons. Note that, due to symmetry, the dipole matrix element from $\ket{0}$ to $\ket{0}$ will be zero. Each matrix element has an important physical meaning, which we will discuss in some detail for both the interacting and non-interacting cases in the Appendix D.

IV.1.1 Example for RESI

Here we use this formulation to recover the equations for the RESI mechanism of NSDI. The process of RESI goes through each states in our two-electron wavefunction ansatz given by Eq. (28). Hence, it goes through the ‘chain’

[TABLE]

where the dipole matrix elements underneath are essential for the transitions between states, so must be included to describe RESI. In Fig. 7, the complete Feynman diagram for RESI is shown. In addition to these dipole matrix elements we include those responsible for self propagating the states

[TABLE]

The remaining dipole matrix elements, $\mathbf{d}(\mathbf{p},\eta),\ \mathbf{d}(\mathbf{p},\mathbf{p}^{\prime})$, and $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\mathbf{p}^{\prime\prime})$, will not contribute significantly to RESI and thus we can neglect their contributions to the corresponding part of the time-dependent Schrödinger equation.

In this example we will take the simplest case where electron interaction is only considered in the necessary step of the excitation of the second electron. As in the case of SFA in the SAE approximation, the crucial point is to determine the most singular parts of the relevant matrix elements and the less singular “rest”. The SFA will then correspond to the systematic expansion of the “rest”.

Below we provide the necessary decompositions of the relevant matrix at two levels: (a) at the level of exact matrix elements, calculated for two-electron Hamiltonian; (b) approximated matrix elements, calculated, neglecting the electron-electron interactions. In this way we will be able to compare with previous SFA results for RESI. We list both expressions (a) and (b) in their most explicit form,

[TABLE]

The convention that we use above is that the less singular parts of the matrix elements with tilde are calculated “exactly”, taking into account electron-electron interaction, while the matrix with subscript 1 stem from an approximate calculation, in which we neglect the electron-electron interactions, so that these matrix elements can be obtained from the corresponding single electron dipole moments, calculated in the SAE approximation. Thus $\hbar{\mathbf{g}}_{1}(\mathbf{p},\mathbf{p}^{\prime})=\hbar{\mathbf{g}}(\mathbf{p},\mathbf{p}^{\prime})$ from the previous chapters.

In addition we have perhaps the most important matrix element that describes re-scattering of the electron accompanied by the excitation of the remaining electron,

[TABLE]

We will keep this matrix element in its entirety, but we will treat it as a perturbation.

Since we treat $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\eta)$ as a first order perturbation, we can then keep only the most singular parts of the remaining matrix elements, neglecting the less singular parts such as $\tilde{\mathbf{g}}(\mathbf{p},\mathbf{p}^{\prime})$. These less singular contributions can be very interesting, leading to Coulomb effects such as distortion of interference structures in ATI Maharjan2006 ; Meckel2008 , and should certainly receive some attention, but in the present instance we will take the most basic form of RESI and neglect them. Thus, we will take, for instance, the single electron re-scattering $\tilde{\mathbf{g}}(\mathbf{p},\mathbf{p}^{\prime})=0$ in the above equations. Substituting these into the integro-differential equations (see Eq. (136) in Appendix E) leads to a much simplified form:

[TABLE]

The above equations contain only the terms relevant for the perturbative solution in the first order in $\tilde{\mathbf{g}}(\mathbf{p},\mathbf{p}^{\prime},\eta)$ – they have thus reduced to a very simple form. Now integral solutions of each of these equations can be formulated, where $\dot{d}(\mathbf{p},\mathbf{p}^{\prime},t)$ is expressed in terms of $\dot{c}(\mathbf{p},\eta,t)$, while $\dot{c}(\mathbf{p},\eta,t)$ is in terms of $\dot{b}(\mathbf{p},t)$, and $\dot{b}(\mathbf{p},t)$ is in terms of $a(t)$, which we assume to know (or we set to unity for the not-too-strong and not-too-long driving pulses). The solutions are as follows:

[TABLE]

where

[TABLE]

The S-Matrix transition amplitude of this process can be related to the above expression in the following way

[TABLE]

We get thus the final result for the RESI amplitude

[TABLE]

In this expression, the three actions will be combined to give the SFA action for the RESI processes. The three integrals with the times $t^{\prime}$, $t^{\prime\prime}$ and $t^{\prime\prime\prime}$ can be associated with ionization, recollision excitation, and final ionization, respectively. The integral over momentum can be related to the intermediate momentum. These integrals can be solved by the saddle point approximation, which makes this problem computationally tractable, as was done in Maxwell2015 ; Maxwell2016 . In that work, the probability distributions for monochromatic fields were calculated for the momentum components parallel to the laser field polarisation, where the components perpendicular to the laser field polarisation are integrated over, as shown in Fig. 8. This shows that different pathways for the RESI process will interfere, by plotting coherent and incoherent sums. These are pathways related to ionization via different excited states; for argon there are six pathways that contribute, as well as pathways related to those via symmetries such as the indistinguishability of electrons.

IV.1.2 Example for EII

Using the same logic it is easy to do the same for the EII of NSDI. Considering similar restrictions as before only the integro-differential equation for $d(\mathbf{p},\mathbf{p}^{\prime},t)$ needs to be changed. This can easily be seen in the EII “chain” (see Fig. 9)

[TABLE]

Here, we neglect the dipole matrix element for recollision excitation RESI contribution given by $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\eta)$, and instead include in a similar way the matrix element $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\mathbf{p}^{\prime\prime})$. We can proceed as before, and now the integro-differential equation for $d(\mathbf{p},\mathbf{p}^{\prime},t)$ can be written out, this time it only depends on $b(\mathbf{p}^{\prime\prime},t^{\prime\prime})$ and is given by,

[TABLE]

Then the solution can be written as

[TABLE]

Ultimately, both the EII and RESI mechanism can be included in $d(\mathbf{p},\mathbf{p}^{\prime},t)$ and it will still be integrable and, as expected, will simply be equal to the sum of these two solutions. Note that EII is much more “classical” than RESI, i.e., that interference effects get washed out upon integration over momentum components perpendicular to the driving field polarization Figueira2004 ; Figueira2004b .

V Large molecules and targets in strong laser fields

In this section we formulate the general problem we want to attack with the help of SFA. The Hamiltonian describing a multi-atomic molecule or atomic cluster has the following general form:

[TABLE]

where the nuclear hamiltonian reads

[TABLE]

with the inter-nuclear potential

[TABLE]

where $N$ is the number of nuclei involved, $i$ enumerates the nuclei, and ${\bf R}_{i}$, ${\bf P}_{i}$ are their positions and momenta, respectively, $Z_{i}$ are the nuclear charges, and $M_{i}$ the nuclear masses. In principle, we could include more complex nucleus-nucleus interactions, taking into account deeply-bound electrons via effective potentials and similar methods. We neglect here the influence of the laser electric field on the nucleus—they are simply too heavy to be affected by the short laser pulses.

The electronic Hamiltonian depends parametrically on the positions of the nuclei, via

[TABLE]

Here $M$ is the number of electrons involved, $i$ enumerates them, and ${\bf r}_{i}$, ${\bf p}_{i}$ are their positions and momenta, respectively. Again, we could replace bare Coulomb potentials by the dressed effective ones. Also, we assume that the target is large, but still smaller than the wavelength, so that a global dipole approximation holds.

Born-Oppenheimer Approximation

In the following we assume that the nuclear motion is slower than that of the electrons, so we use the Born-Oppenheimer approach. To this is end, we first determine the electronic wave function, $\Psi_{\mathrm{el}}(\{{\bf r}_{i}\}_{i=1}^{M},t;\{{\bf R}_{i}\}_{i=1}^{N})$, in Dirac’s notation denoted as $|\Psi_{\mathrm{el}}(t)\rangle$, that fulfils TDSE with fixed nuclear positions

[TABLE]

We define then instantaneous electronic potential for the nuclei

[TABLE]

and treat the motion of the nuclei classically by solving the resulting Newton equations

[TABLE]

The solutions of these classical equations are then inserted into the electronic TDSE, and so on. In general, it has to be done self-consistently. We will discuss below a couple of cases when some simplifications are possible.

Single Active Electron Approximation

The theory formulated above can be reduced to a single-electron TDSE using the SAE approximation. In principle, it can be done in the same way as it is done in the static case for atoms or molecules. The only difference is that we now have to consider the fact that the effective potential must now be time-dependent through the parametric dependence on the nuclear coordinates,

[TABLE]

While calculations of $V_{\mathrm{eff}}$ for atoms belong still to the domain of the atomic physics, calculations of the effective potential for molecules, especially in dynamical situations, clearly require the use of methods from molecular physics and quantum theoretical chemistry.

V.1 Strong Field Approximation for quenched molecules

The equations of the above section are very complex. There are some situations, however, when they can be radically simplified. One example of such a situation is the instant quench, in which the molecule is suddenly stripped of, say, one of the electrons, or photoexcited to a certain well-defined state. This can be achieved, for instance, applying an ultrashort attosecond XUV or soft X-ray pulse to the molecule. Right after the pulse, the molecule will find itself in the ground state corresponding to one missing electron, or in the well-defined excited state. In both situations, the nuclei configuration will be by no means stable. The molecule will start to vibrate, rotate, and maybe even dissociate.

If the excitation occurs to a weakly bound molecular state, the following vibrations or dissociation will occur on a rather slow time scale of $100\text{\,}\mathrm{f}\mathrm{s}$ to $1\text{\,}\mathrm{p}\mathrm{s}$. In that case, the HHG or ATI caused by an intense few-femtosecond pulse may be used for an instant imaging of the dynamically changing molecular structure (for seminal experiments see Refs. WenLi2008, ; Worner2010, ). If the electron removal or excitation occurs to a strongly-bound state, the resulting dynamics might be much faster: stripping of electrons, for instance, might lead to dissociation completely controlled by the Coulomb forces, and occurring on the timescales of an atomic unit (fractions of a femtosecond). These are the situations we want to consider in this section.

SFA and molecular dynamics

If we then apply a short femtosecond laser pulse in the mid-infrared range, we may expect that, similarly to the standard HHG or ATI processes, the femtosecond laser induced electronic dynamics will not affect the intrinsic molecular dynamics. That means that, from the point of view of nuclei, we can replace the electronic Hamiltonian

[TABLE]

with

[TABLE]

The Born-Oppenheimer Newton equations for the nuclei can then be solved self-consistently, as

[TABLE]

where

[TABLE]

Assuming the that the ionization during the process is weak, the contribution of the continuum part of the electronic wavefunction will give negligible contribution to the electronic energy, so that

[TABLE]

where $|\Psi_{0}(t)\rangle$ is the time-dependent ground electronic state, and $a(t)$ is the probability amplitude of being in this state. If we know $|a(t)|^{2}$ (for instance, if we can assume that $|a(t)|^{2}\simeq 1$), then the solutions can be simply introduced into Eq. (75), and one can then calculate explicitly—i.e. without self-consistency conditions—both $|\Psi_{0}(t)\rangle$ and the corresponding continuum functions that fulfill

[TABLE]

Both of these functions depend explicitly on time through the time dependence of the positions of the nuclei. Note, that the equations can be even more simplified if we can simplify the effects of $E_{0}$ in Eq. (74)—the equations will not even require self-consistency! For instance, in the case of stripping of, say, $K$ electrons, for large internuclei distance, the only effect of $E_{0}$ in Eq. (74) will be to screen the nuclei charges, that is replace $Z_{i}$s by $\tilde{Z}_{i}$s, where $\sum_{i}Z_{i}=\sum_{i}\tilde{Z}_{i}-K$.

SFA for a quenched molecule

The expression derived above implicitly assumes that we proceed in fact as in Section III. That is, we write the full electronic wave function as

[TABLE]

where we set $I_{p}(t)=-E_{0}(t)$. The new effect here is $\phi_{B}(t)=\langle\Psi_{0}(t)|\partial_{t}\Psi_{0}(t)\rangle$—the Berry phase arising from projecting/expanding the electronic wave function in the time dependent basis. The equations still have practically the same form as before; for instance the direct transition amplitude fulfills:

[TABLE]

where we have included now the Berry phase in $\tilde{I}_{p}(t)=I_{p}(t)+\hbar\phi_{B}(t)$. Note that the Berry phase is nonzero if and only if the ground-state wavefunction is complex. This typically happens if the time-reversal symmetry is broken, i.e. for instance in the presence of a magnetic field or a so-called “artificial” gauge field. Also, the matrix element now depends explicitly on time, through the time dependence of the positions of the nuclei.

Generalized SFA expression for a quenched molecule

The above equations for the electronic dynamics, together with the Newton equations (73) and (74), as well as the expressions for the electronic energy (75-77), allow us to derive thus:

• •

The direct ATI amplitude:

[TABLE]

• •

The re-scattering amplitude:

[TABLE]

where

[TABLE]

• •

The time-dependent dipole moment

[TABLE]

Note that if $|a(t)|^{2}$ is “known”, then the solutions do not require self-consistency. Otherwise, they have to be obtained in the manner discussed below. In the Appendix F we discuss a toy model of 1D quenched $\rm H_{2}^{+}$ molecule, looking at qualitative and even semi-quantitative effects in HHG.

V.2 SFA for large targets

Here we consider another situation: the molecule (a large target) is initially in the ground state, and is impinged by an intense, short (few-cycle) laser pulse in the mid-infrared range. This pulse causes the ionization of the single active electron, and induces thus structural dynamics of the target, i.e. the motion of the nuclei. Amazingly, the expressions describing the quantities of interest are exactly the same as in the previous section. The way to obtain them, however, is much more complex: now we have to determine the evolution of ${\bf R}_{i}(t)$ and ${\bf P}_{i}(t)$ simultaneously and self-consistently with the dynamics of the electronic wave function, $|\Psi_{\mathrm{el}}(t)\rangle$.

The protocol to follow is thus:

Calculate the electronic state (the ground state of $\hat{H}_{\mathrm{eff}}$),

[TABLE]

for the initial positions of the nuclei $\pm\Delta{\bf R}_{i}$ (to be able to calculate gradients). 2. 2.

Propagate the equations for nuclei,

[TABLE]

to the next time instant, $t$. Calculate the new $\{{\bf R}_{i}(t)\}_{i=1}^{N}$. 3. 3.

Calculate the new electronic state $|\Psi_{\mathrm{el}}(t)\rangle$. This is the state propagated using $\hat{H}_{\mathrm{eff}}$,

[TABLE]

for the actual positions of the nuclei $\pm\Delta{\bf R}_{i}$ (to be able to calculate gradients). Note that this propagation should be done using the SFA ansatz (78)). Calculate then

[TABLE]

for the actual positions of the nuclei $\pm\Delta{\bf R}_{i}$ (to be able to calculate gradients). 4. 4.

Calculate

[TABLE] 5. 5.

Go to 2.

Obviously, the above procedure is quite complex, but it does not present giant numerical challenges, and it is relatively straightforward to implement. Evidently, it is much easier than solving the TDSE involving classical (Born-Oppenheimer) dynamics of the nuclei, the feasibility of which is not entirely obvious.

V.3 SFA and quantum molecular dynamics

The use of classical Newton equations for molecules in dissociation or vibration processes might be questionable. There is a simple method of including certain aspects of the quantum motion of molecules that we describe now. Our starting assumption is the generalization of the SFA ansatz to the full wave function:

[TABLE]

where $\xi_{0}$, $\xi_{1}$ are the normalized wave functions of the nuclei for the molecule with $M$, $M-1$ electrons correspondingly.

As before $\hat{H}_{\mathrm{eff}}=\hat{H}_{0}-e{\bf E}(t){\bf r}$, whereas

[TABLE]

The electronic ground state is now time-independent, but it does explicitly depend on the nuclear positions via

[TABLE]

Similarly, the states in the continuum do not depend on time, but on the nuclear positions, entering via $\hat{H}_{0}$ as

[TABLE]

We still use the Born-Oppenheimer approximation, but in the quantum version. Also, we use different Hamiltonians for the non-ionized and ionized part of the molecular electronic dynamics. Thus, for $\xi_{0}({\bf R}_{1},\ldots,{\bf R}_{N})$ we use

[TABLE]

while for $\xi_{1}({\bf R}_{1},\ldots,{\bf R}_{N})$ we use simply

[TABLE]

We neglect here the laser part of the electronic energy, as well as the kinetic energy of electrons in the continuum. Note that the equation (95) can be solved without any self-consistency conditions. As in the previous sections, equation (94) can also be solved that way, provided that the time dependence of $|a(t)|^{2}$ is known.

The last point is the derivation of the SFA equation. To this end we assume that the quantum fluctuations of the nuclear positions are small, and replace the ${\bf R}_{i}$ dependence in $|\Psi_{0}(t)\rangle$ by the average $\bar{\bf R}_{i}(t)=\int d^{3}{\bf R}\>{\bf R}_{i}|\xi_{0}({\bf R},t)|^{2}$. Similarly, we replace the ${\bf R}_{i}$ dependence in the continuum part by the average $\bar{\bf R}_{i}(t)=\int d^{3}{\bf R}\>{\bf R}_{i}|\xi_{1}({\bf R},t)|^{2}$; after that trick, the SFA equations can be projected on the normalized functions $\xi_{0,1}(t)$. This leads to the following modified equations:

[TABLE]

and

[TABLE]

As we see, the final equation depends only on the overlap $\langle\xi_{1}(t)|\xi_{0}(t)\rangle$, generally called the nuclear autocorrelation function, which, despite the fact that the positions of the nuclei in each branch of the process are quite “classical”, might become very small as the positions of the nuclei in the two channels change. This can then seriously limit the HHG and LIED signals from the process, both for the direct and re-scattering parts Lein2005 ; Baker2006 ; Patchkovskii2009 .

VI SFA for solids

The most important challenges of strong laser field physics concern using strong laser fields, be it the ones used in atto-science, or in free electron lasers (FELs), or in X-ray lasers, to characterize the statics and dynamics of the targets. The complexity of those targets grows with their “size”: in a sense, atoms are the simplest, molecules are much more complex, and solid state and condensed matter systems are perhaps the most challenging, especially if we consider strongly-correlated systems.

The question of the optical response of solids has been a long-lasting subject of investigation, from the basics of nonlinear optics Boyd2003 , through studies of femtosecond X-ray sicence and the optical properties of semiconductors Pfeifer2006 ; Haug2008 ; it is directly related to the fundamentals of attosecond optics Chang2011 and contemporary photoelectron spectroscopy Hufner2013 . Before we focus on solids, let us summarize shortly the visions with atoms and molecules over the past decade.

• •

Visions with atoms. In a sense, the basis for all of these applications is Ref. Lewenstein1994, . However, to use this theory in practice, one has to include important physical and quantitative effects, such as Coulomb focusing Brabec1996 . In the last ten years or so, the primary focus of the work on atoms was on dynamics, in particular the dynamics of the tunnelling processes, including the birth time of HHG photons Dudovich2006 , collective many-electron dynamics Shiner2011 , explicit probing of the atomic wave function Shafir2009 , or time resolved auto-ionization Wang2010 ; the Lund group proposed to study the birth time of attosecond pulses using two-color fields He2010 ; Dahlstrom2011 . Recent highlights of attacking atomic challenges with ultra-strong and ultrashort laser pulses include studies of the time when an electron exits a tunnelling barrier Shafir2012b , observations of electron propagation and screening on the atomic length scale Neppl2015 , and attosecond tunnelling interferometry Pedatzur2015 . For a very recent discussion of these areas, see the beautiful review on ultrafast holographic photoelectron imaging in Ref. Figueira2019, .

• •

Visions with molecules. These visions, obviously, come back to the seminal work of Itatani et al. Itatani2004, , which was focused on the reconstruction of the static molecular wave functions from measured HHG spectra. In recent years these visions moved very strongly forward. A few important examples (though by no means totally representative) concern interferometry of multi-electron dynamics in molecules with HHG Smirnova2009 , following time-resolved chemical reactions  Worner2010 , performing generalized molecular orbital tomography Vozzi2011 , or observing the bending of molecules in real time – the Renner-Teller effect using laser-induced electron diffraction Amini2019 .

• •

Visions with solids. This is probably the most challenging area, especially if one thinks about strongly correlated systems or systems with topological order.

The pioneering works on theory and experiments concerning HHG in solids go back to the works of P.C. Becker et al. and, independently, G. Petite and Ph. Martin in CEA Saclay Becker1988 ; Martin2000 ; Pisani2001 , who developed and analyzed a “happy electron model”. A more rigorous approach to SFA in solids was followed by F. Faisal’s group few years later Faisal2005a ; Faisal2006 . Of course, the basis to analyze HHG or LIED in solids must involve precise and detailed static structure calculations Goano2007

The whole area started to grow exponentially roughly ten years ago. Important early theory developments were done by the S. Koch’s group Golde2008 . Around 2010, the “volcano” erupts: Goulielmakis and coworkers report real-time observation of valence electron motion Goulielmakis2010 , while S. Ghimire et al. publish their seminal paper on observation of high-order harmonic generation in a bulk crystal in 2011 Ghimire2011 . This is followed by a true explosion of interest in this area: HHG generation and propagation in crystals Ghimire2012 , controlling dielectrics with light Schultze2013 , attosecond band gap dynamics in silicon Schultze2014 , electron propagation and dielectric screening on the atomic length scale Neppl2015 , and efficient HHG in liquids Heissler2014 ; Luu2018-liquids .

While the disorder of the liquid phase makes spectroscopic studies of its structure more complicated, the rigid structure of solids allows for a much deeper understanding of the emission mechanisms as well as, potentially, broad and detailed high-harmonic spectroscopy Luu2015 , real-time observation of interfering crystal electrons in HHG Hohenleutner2015 , HHG in silicon Vampa2016 , and, more recently, interferometric measurements of the dipole phase in high harmonics from solids Lu2019 . This is because, when an electron is released in a solid by a strong-field excitation from the valence to the conduction band, it explores a continuum which contains much more structure than the quadratic band of a free electron (as shown in Fig. 10 and explained in depth in the review in Ref. Vampa2017, ), and the dispersion induced by this structured continuum produces non-harmonic motion which leads to the emission of so-called “intraband” emission (cf. Ref. Schubert2014, ).

Moreover, in atoms, the ionized electron leaves behind a stationary hole which is bound to the parent ion and cannot be displaced in space, but in solids this is no longer the case, and the motion of the hole in the valence band also needs to be considered. Nevertheless, when the hole and electron meet, they can recombine and emit so-called interband harmonics, exactly as in the atomic case (cf. Refs. McDonald2015, ; Wu2015, ). However, despite that similarity, there are important differences, since that recombination can happen away from the origin, the electron and hole trajectories are subject to more complicated dynamics in their dispersive bands, and the bands themselves contain nontrivial parallel-transport effects that produce, through a nonzero geometrical phase, additional ‘anomalous’ velocity terms that also contribute to the harmonic emission.

The theory of HHG in solids is based on the “classical” work of Blount Blount1962 . To a good approximation, the electronic dynamics in a solid driven by a strong low-frequency laser pulse is governed by the semiconductor Bloch equations (SBE) HaugKoch ; Picon2019 , as derived for instance in Ref. Chacon2018, :

[TABLE]

where $n_{m}$ is the population in band $m$, $\pi$ is the inter-band coherence (where we assume for simplicity a two-band model), $w=n_{c}-n_{v}$ is the population difference between the valence and conduction bands, $s_{c}=-s_{v}=1$, $\mathbf{d}_{mm^{\prime}}(\mathbf{k})$ is the inter-band dipole moment, and $\bm{\xi}_{g}(\mathbf{k})=\bm{\xi}_{c}(\mathbf{k})-\bm{\xi}_{v}(\mathbf{k})$ is the difference in the Berry connections of the two bands (Chacon2018, , supplemental material). These variables give rise to the harmonic emission via the intra and inter-band components of the total current,

[TABLE]

Let us discuss with a little more detail the theoretical description of HHG in solids, and its similarities to the SFA for atoms. The basis for the SFA is the Keldysh theory of tunnelling ionization Keldysh1965 . Keldysh considered the ionization of atoms in low frequency fields, and thus could use a quasi-static approximation to describe the tunnelling processes. It is worth noting that in the orginal paper from 1965 Keldysh already could generalise his theory to the solid state, and estimate tunneling ionization in solids. As detailed in Section I, Keldysh theory was generalized, on one hand, by Ammosov, Delone and Krainov Ammosov1986 to describe ionization rates from complex atoms, and, on the other hand, by Faisal and Reiss to describe ATI Faisal1990 ; Reiss1980 ; it was also a stimulus for the present formulation of the SFA Lewenstein1994 .

The optical response of solids (insulators and semiconductors, and to some extent even metals) to strong laser fields is described by the semiconductor Bloch equations (98, 99). In general, these equations are too complex to solve directly other than numerically, but, as in the atomic case, in some situations they allow us to calculate an approximate expression for the electronic current.

In the atomic case, approximate analytic solutions of the TDSE are possible in the SFA, where the dynamics is reduced to the ground state and continuum states (with bound states other than the ground state getting neglected), and in the simplest approach Coulomb interactions are ignored. For solids, on the other hand, the situation is different from the very beginning, since we start with a many-body system of electrons in the valence band of the crystal. In some cases, it is possible to use effective-electron descriptions of the system (and even then we deal with effectively non-interacting quasi-particles) and we have to include all of the particles that fill the valence band. This is analogous to considering many filled bound states in an atom. Perhaps the best analogy to the atomic case is achieved by considering the solid-state dynamics using the Wannier description in the valence bond, coupled to Bloch waves in the conduction band Osika2017 .

If the effective Fermi-liquid non-interacting quasi-electron description is a good approximation, then the response of the solid to a strong laser pulse is well-described by a theory similar to the Keldysh framework and to the SFA. On the other hand, if we deal with strongly-correlated systems that cannot be described by the effective mean-field-like description, then one has to rely on much more rigid theoretical tools, including the effects of Coulomb interactions. This is a challenging area, with only a few ab initio numerical studies Silva2018 ; Takayoshi2019 , though some additional headway can be made using DFT Bauer2018 ; Drueke2019 .

Within the effective single-active-electron SFA theory, however, there is substantially more that one can say analytically, at least on an approximate footing. Here, the current can be expressed Chacon2018 in the form

[TABLE]

where $S({\bf K},t,t^{\prime})$ is the so called quasi-classical action for the electron-hole and is defined according to:

[TABLE]

This forms the heart of the SFA description of high-harmonic emission in a crystalline solid, and it contains all of the semiclassical dynamics for the electron trajectories in the conduction band as well as the hole’s trajectory in the valence band Vampa2014 ; Vampa2015 ; Vampa2015a ; Vampa2017 . More recent theoretical developments concern the effects of multiple conduction bands Hawkins2015 .

Traditional treatments of HHG in solids have worked in one-dimensional configurations where the geometrical-phase element of (103) can be set to zero using an appropriate choice of gauge for the Bloch-function basis, in which case the kinematics of the electron and hole wavepackets are fully determined by the band structure Vampa2015 ; Hawkins2013 ; Higuchi2014 . However, there are materials where this gauge transformation is not possible, due to the existence of a Berry curvature on one or more of the bands, and the geometric-phase terms in (103) cannot be neglected. This Berry curvature is crucial for a wide array of solid-state effects vonKlitzing1986 ; Haldane2017 ; Xiao2010 , and it is the driving ingredient of many nontrivial phenomena. Its presence in the SFA harmonic-emission current (102) means that it can in principle be measured via HHG observables, and indeed recent experimental Liu2017 ; Luu2018 and theoretical Chacon2018 ; Silva2018a works show that this is the case; we showcase in Fig. 11 some recent results on that front. A review on high-harmonic generation from solids by Ghimire and Reis was published very recently in Ref. Ghimire2019, .

Acknowledgements.

We thank Sougato Bose, Paul Corkum, Anne L’Huillier, Yann Mairesse and Pascal Salières for enlightening discussions. This work was supported by National Science Centre Poland-Symfonia Grant No. 2016/20/W/ST4/00314. We acknowledge also the project Advanced Research Using High Intensity Laser Produced Photons and Particles (No. CZ. 02.1.01/0.0/0.0/16_019/0000789) from the European Regional Development Fund (ADONIS), Ministerio de Economía y Competitividad through Plan Nacional (Grant No. FIS2016-79508-P FISICATEAMO, and Severo Ochoa Excellence Grant No. SEV-2015-0522), funding from the European Union Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Grant Agreement No. 641272 Laserlab-Europe (Grant No. EU-H2020 654148), Fundació Privada Cellex, and Generalitat de Catalunya (Grant No. SGR 1341 and CERCA Programme). N.S. was supported by the Erasmus Mundus Doctorate Program Europhotonics (Grant No. 159224-1-2009-1-FR-ERAMUNDUS-EMJD). N.S., A.C., A.D., E.P., and M.L. acknowledge ERC AdG OSYRIS and NOQIA, and EU FETPRO QUIC. K.A. and J.B. additionally acknowledge ERC AdG Transformer (788218), Ministerio de Economía y Competitividad through Plan Nacional FIS2017-89536-P and Agencia de Gestió d’Ajuts Universitaris i de Recerca (AGAUR) with 2017 SGR 1639. E.P. acknowledges Cellex-ICFO-MPQ Fellowship funding. A.D. is financed by a “Juan de la Cierva” fellowship (IJCI-2017-33180). A.P. acknowledges funding from Comunidad de Madrid through TALENTO grant ref. 2017-T1/IND-5432. J.A.P.-H. acknowledges Spanish Ministerio de Economía, Industria y Competitividad (MINECO) through the PALMA Grant No. FIS2016-81056-R, from LaserLab Europe IV Grant No. 654148, and from the Regional Government, Junta de Castilla y León, Grant No. CLP087U16 and the Consolidated Research Unit 167, recognized by the Junta de Castilla y León. C.F.M.F. and A.S.M. acknowledge funding by the UK Engineering and Physical Sciences Research Council (EPSRC) under grants EP/J019143/1 and EP/P510270/1, the latter within the remit of the InQuBATE Skills Hub for Quantum Systems Engineering. A.Z. acknowledges funding from UK EPSRC under the grant EP/J002348/1.

Appendix A Time dependent ADK rates

This is the first section in which we remind the reader about a possible way to evaluate/estimate the amplitude of the ground state, $a(t)$, which is the essential ingredient of our formulation of the SFA. The ADK rates apply in principle for arbitrary atoms, but the theory can be extended to molecules Tong2002 . This approach can be used for pulses which are not too short: the rate has to have time to “define itself”. The static expressions for the ADK rates are

[TABLE]

where $F$ is the peak value of the laser electric field, $n^{*}=Z/\sqrt{I_{p}}$, and $Z$ being the charge of the atomic/ionic core. The other symbols are

[TABLE]

with $\Gamma(\cdot)$ denoting the Gamma function, $l^{*}=n^{*}-1$, and

[TABLE]

with $l,\,m$ denoting electrons initial orbital and magnetic quantum numbers.

For shorter pulses, we obtain $W_{\mathrm{ADK}}(t)$ replacing $F\to F(t)=\mathcal{E}_{0}\>f(t)$ (the rates change adiabatically with the pulse envelope). If the laser frequency is even smaller, and the pulse shortened we can even replace $F\to F(t)=|\mathbf{E}(t)|$, i.e. the actual value of the electric field.

Appendix B Ground state amplitude according to SFA

Inserting the expression (15) in the equation for $a(t)$, we obtain

[TABLE]

where

[TABLE]

In principle solving equation (107) in the static, quenched or even self-consistent case presents no basic difficulties. If the pulse is longer, and the rate of change of $a(t)$ slow, we may replace $a(t^{\prime})\to a(t)$ and obtain explicitly

[TABLE]

where the SFA rate is

[TABLE]

Appendix C Model atom and molecule

The paradigm examples of separable potentials used in atomic physics are zero-range potentials. In 1D, a Dirac delta potential can be used for this purpose. In 3D, the Dirac delta must be regularized—that is why the celebrated pseudo-potential must be used. It has found multiple applications in the many-body theory of ultracold atomic gases (cf. Refs. StringariPitaevski, ; LewensteinQuantumSimulators, ); in strong-field physics it was elaborated by W. Becker and his collaborators Becker1994a . In 2D the situation is more complex due to the logarithmic divergences (cf. Ref. Busch1998, and references therein).

One should stress, however, that the use of separable potentials has a long history in strong-field physics Faisal1987 ; Faisal1987b ; Faisal1989 ; Faisal1990 ; Faisal1991 ; Faisal1993 ; Faisal1993b ; Radozycki1993 ; Kaminski1989 ; Klaiber2005 ; TetchouNganso2007 ; TetchouNganso2011 ; TetchouNganso2013 ; Galstyan2017 ; Galstyan2018 . Many of these papers deal with zero-range models, some with general separable potentials, but typically using laser fields of constant strength and circular polarization. For such a case, one can transform the problem to a rotating frame in which the Hamiltonian is time-independent. The recent papers by Galstyan et al. Galstyan2017 ; Galstyan2018 are perhaps the closest to the approach developed by us in Refs. Noslen1, ; Noslen2, ; Noslen3, ; Noslen4, ; Noslen-diss, .

In this section we discuss a useful non-local separable potential with the purpose of applications for atoms, but most importantly for large molecules. The idea will be to compute both the direct and the re-scattering transition amplitudes Lewenstein1995 . These terms involve the dipole and the continuum-continuum matrix elements defined by Eqs. (12) and (13). Then, our main task will be devoted to find analytically the wavefunctions for the ground and scattering states of our model potential. The Hamiltonian, $\hat{H}(\textbf{p},\textbf{p}^{\prime})$, of the atomic system in the momentum representation can be written as:

[TABLE]

where the first term on the right-hand side is the kinetic energy operator, and the second one is the non-local potential $\hat{V}(\bf{p},{\bf p}^{\prime})$. We use non-local separable potentials that can be understood as sums of projectors on certain states. They generally have the form

[TABLE]

where $\gamma$ is a coupling constant.

When we model molecules, each of the orbitals is typically centered in real space at the positions of the nuclei, ${\bf R}_{i}$. In the momentum representation it translates to $\phi_{i}(\bf{p})\exp(i\bf{p}\cdot{\bf R}_{i}\tilde{\phi}_{i}({\bf p})$, where $\tilde{\phi}_{i}({\bf p})$ is a “smooth” function, with the Fourier transform centered at ${\bf R}=0$. The above potential has generically $M$ bound states, so one can model with its help not only the ground state of the molecule in question, but even some of its excited states. Alternatively it may be used to model multielectron molecules.

Ground state

In the present Appendix we will consider $M=1$ only, so the models contain a single bound state, but we will explore the full richness of the orbital $\phi({\bf p})$ to model multicentered ground states of, in principle, arbitrary molecules. By using the non-local separable Hamiltonian, we write the stationary Schrödinger equation as follows:

[TABLE]

where $E$ denotes the energy of the wavefunction $\Psi({\bf p})$. Note that we have defined the non-local potential as $\hat{V}(\textbf{p},\textbf{p}^{\prime})=-\gamma\phi(\textbf{p})\>\phi(\textbf{p}^{\prime})$, which describes the attraction between the electron and the nucleus Lewenstein1995 . This potential has been chosen such that it assures analytical solutions of the continuum or scattering states, i.e. for states with energies $E>0$. Note that the ground state can also be calculated analytically. The parameter $\gamma$ is a constant that, as we will see, determines the energy of the ground state. The shape of the ground state, however, can be controlled to a high degree by the choice of a suitable auxiliary function $\phi({\bf p})$, which may correspond to a multicenter molecular orbital of arbitrary shape.

For the ground state, $\Psi_{0}({\bf p})$, we solve the stationary Schrödinger equation in the momentum representation:

[TABLE]

From Eq. (114) we easily determine the ground state:

[TABLE]

where, $\mathcal{N}$ denotes a normalization constant. Multiplying the last formula by $\phi({\bf p})^{*}$, and taking the volume integral on p, we obtain an equation that determines the ground state energy,

[TABLE]

The solution of the last integral in Eq. (116) gives us the relation between the parameters $I_{p}$, and $\gamma$ given $\phi({\bf p})$. The normalization constant fulfils

[TABLE]

Scattering waves

Let us consider the scattering wave, $\Psi_{{\bf p}_{0}}({\bf p})$, with asymptotic momentum ${\bf p}_{0}$, as a coherent superposition of a plane wave and an extra correction $\delta\Psi_{{\bf p}_{0}}({\bf p})$:

[TABLE]

This state has an energy $E={{\bf p}_{0}^{2}}/2m$. Then, the Schrödinger equation in momentum representation reads:

[TABLE]

To solve analytically the last equation, we apply elementary algebra and the following Dirac delta distribution properties: $(\frac{p^{2}}{2m}-\frac{p_{0}^{2}}{2m})\>\delta(\textbf{p}-\textbf{p}_{0})=0$, and $(\frac{p^{2}}{2m}-\frac{p_{0}^{2}}{2m})\>\delta(\frac{p^{2}}{2}-\frac{p_{0}^{2}}{2m}\pm i\epsilon)=0$. Finally, the correction $\delta\Psi_{{\bf p}_{0}}$ reads:

[TABLE]

Here, $\epsilon$, is the (positive) regularization parameter to avoid the divergence at ${\bf p}={\bf p}_{0}$, and $A(ip_{0}+\epsilon)$ is the function $A(k)=2m\gamma\int\frac{|\phi({\bf p})|^{2}}{p^{2}+k^{2}}\textit{d}^{3}\textbf{p}^{\prime}$ evaluated at $k=ip_{0}+\epsilon$. The sign in front of $\epsilon$ defines the asymptotic behaviour of the scattering solutions. For the present choice, the scattered part of the wavefunction in the position representation behaves asymptotically as $\exp(-\textit{i}p_{0}r)$, i.e. as a incoming spherical wave, which is the correct behaviour for the scattering wave functions, describing the states with asymptotic outgoing momentum ${\bf p}_{0}$. The final expression for the scattering wave functions is:

[TABLE]

Dipole matrix element

The dipole matrix element is

[TABLE]

Tedious, but elementary algebra leads to:

[TABLE]

Continuum-continuum transition matrix element

Let us consider the scattering waves obtained in Eq. (121) and evaluate the continuum-continuum transition matrix element of Eq. (13), i.e.

[TABLE]

Again after tedious (though straightforward) calculations we obtain:

[TABLE]

All of the above formulae are directly and easily generalized to the case when the Hamiltonian is time-dependent.

Appendix D Dipole matrix elements

In the remaining appendices, we turn to the objects of the two-electron theory and their properties.

D.1 Dipole Matrix Elements from the ground state

The dipole matrix element given by the function $\mathbf{d}$ deals with dipole transitions from the two electron ground state into all other admissible two electron states. This is the analogue of the function by the same name for one electron case, given by Eq. (12), but now there are three non-zero variants. The first, bearing the strongest resemblance to the one electron case, can be expressed as a sum of two terms

[TABLE]

Each term will have two contributions, (1) direct laser induced ionization of the electron acted on directly by the operator, while the other electron remains bound in a ground state and (2) correlated ionization, where the the action of laser on one electron is transferred to the other via electron-electron interaction, this could be through elastic collision. This last contribution is expected to be small and will be zero in the non-interacting case, where Eq. (126) simplifies to

[TABLE]

The final expression is written in terms of 1-electron states and operators. This dipole matrix element is vital to most strong field processes as it will describe the initial tunnelling step, hence it is necessary to model NSDI for both the EI and RESI mechanism. The subscript denotes whether it is the first or second electron to be ionized, this is significant as the second will have a much larger ionization potential, which will make it less probable and for both mechanisms of NSDI we will neglect this contribution.

The other two variants of the $\mathbf{d}$ dipole matrix element require at least some electron interaction in all contributing processes and so will both be zero in the non-interacting case. The element $\mathbf{d}(\mathbf{p},\eta)$ will again have two contributions, in the first the operator acts on an electron to ionize it and the electron-electron correlation causes the other electron to be excited in a ‘shake-up’ process. This would be the main term involved in the previously proposed shake-off mechanism for NSDI Fittinghoff1992 , which has since fallen out of favour in preference of the re-collision mechanisms EII and RESI. This will be the dominant of the two contributions. In the second (presumably less likely) transition, the laser excites an electron and the electron-electron correlation causes the second electron to be fully ionized. In both these possibilities the energy transfer between electron could be through direct collision or any other electron-electron interaction, however in the first only some of the electron’s energy is transferred, while in the second scenario most of the electron’s energy will be transferred to the other one.

For the dipole matrix element $\mathbf{d}(\mathbf{p},\mathbf{p}^{\prime})$ the only transition is one where the laser-induced ionization of an electron through interactions between the electrons causes the other electron to ionize as well—in fact, this matrix element, when dressed in the laser field, leads to collective tunnelling, as discussed in Ref. Lewenstein1986a, .

D.2 Dipole Matrix Elements from the scattering and ground state

The dipole matrix elements from the two-electron scattering-ground state $\ket{\mathbf{p},0}$ are given by the function $\mathbf{g}$, there are three possible variants. The first, $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime})$ relates to the transition from the two electron continuum-ground state $\ket{\mathbf{p},0}$ to an alternative continuum-ground state $\ket{\mathbf{p}^{\prime},0}$. The leading contribution will typically come from the laser induced change of momentum of the continuum state, this term should play a strong role as we expect strong coupling between the laser and continuum electrons. Continuum-continuum transitions can, as in the one electron case, also involve contributions via interaction through the single electron potential. In Eq. (13) this was described by splitting of the continuum-continuum matrix element into two parts in the one-electron case. Alternatively, there is the strongly correlated and less likely process, where the laser acts on the bound electron which through electron-electron interaction changes the momentum of the scattering state for the other electron, without changing the state of the original bound electron. This would be quite an exotic case and generally it is a reasonable approximation to assume the two electrons in this state are somewhat physically separated. Of course, when an electron recollides with its parent atom or molecule then there will be much overlap and this term could contribute to elastic recollision processes such as high-order above-threshold ionization (HATI). In Fig. 12 the four pathways discussed are shown in the form of Feynman diagrams; the non-interacting cases are given by panels a) and b).

This dipole moment can be considerably simplified if we consider non-interacting electrons

[TABLE]

This can be treated as before by Eq. (13),

The dipole matrix element $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\eta)$ deals with transitions from $\ket{\mathbf{p},0}$ to $\ket{\mathbf{p}^{\prime},\eta}$. This will have contributions from scattering states interacting with the laser to change the moment from $\mathbf{p}$ to $\mathbf{p}^{\prime}$, which simultaneously through the electron interactions also leads to excitation of the other electron from its ground state into an excited state with principle quantum number $\eta$. It is crucial to include this to allow for the excitation step of the second electron in the RESI mechanism of NSDI. Alternatively, there is the contribution that the excitation of an electron via the laser also results in the momentum change of the scattering state from $\mathbf{p}$ to $\mathbf{p}^{\prime}$. In the non-interacting case the first pathway is has no contribution as the ground state electron cannot be excited without electron correlation. The second pathway is non-zero in the case $\mathbf{p}=\mathbf{p}^{\prime}$ and the dipole matrix element can be written as

[TABLE]

The excitation via electron interaction and non-interacting pathways are depicted in panel a) and b) of Fig. 13, respectively.

It is a very similar situation for the dipole matrix element $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\mathbf{p}^{\prime\prime})$, except the excited state is replaced by a scattering state with momentum $\mathbf{p}^{\prime\prime}$. There are the same kind of contributions, one where the change of momentum of the scattering states leads through electron-electron interaction to ionization of the bound state, this is the crucial ionization step of the second electron in the EI mechanism of NSDI. Another possibility is that the laser induced ionization of the bound electron leads to change of momentum of the scattering state, which could happen through the collision of the two ionized electron. In the non-interacting case the first contribution is zero, while the second is non-zero again if $\mathbf{p}=\mathbf{p}^{\prime}$ which means the matrix element can be simplified to

[TABLE]

The subscript two denotes that this is the laser induced ionization of a second electron, i.e. this electron comes from a +1 ion, thus it is much less probable than the first ionization and for our purposes, where we consider NSDI, we will not consider such contributions. In Fig. 14 we show the EI recollision excitation pathway and also the non-interacting pathway, in panels a) and b), respectively.

D.3 Dipole Matrix Elements from the scattering and excited states

The dipole matrix elements from the two electron continuum-excited state $\ket{\mathbf{p},\eta}$ are given by the function $\mathbf{h}$, there are two variants. The function $\mathbf{h}(\mathbf{p},\eta,\mathbf{p}^{\prime},\eta^{\prime})$ deals with transitions to alternate continuum-excited states. This will have contributions such as the laser induced recollision of the scattering state results in the excited electron moving to a another excited state. This will be non-zero in the non-interacting case if either $\mathbf{p}=\mathbf{p}^{\prime}$ or $\eta=\eta^{\prime}$ and the dipole matrix element can be written as

[TABLE]

The other function $\mathbf{h}(\mathbf{p},\eta,\mathbf{p}^{\prime},\mathbf{p}^{\prime\prime})$ deals with transition to continuum-continuum state $\ket{\mathbf{p}^{\prime},\mathbf{p}^{\prime\prime}}$. This will have a contribution where the laser induced recollision of the scattering state changes the momentum from $\mathbf{p}$ to $\mathbf{p}^{\prime}$ and the excited electron is ionized through electron-electron interaction. Alternatively, the laser induced ionization of the excited electron can change the momentum of the continuum electron through their interaction. The first contribution vanishes in the non-interacting case, while the second contribution is non-zero if $\mathbf{p}=\mathbf{p}^{\prime}$. This what leads to the final ionization step in the RESI mechanism of NSDI, where the excited electron tunnel ionizes via the dipole interaction. This dipole matrix element can be written as

[TABLE]

D.4 Dipole Matrix Elements from the scattering states

This dipole matrix element is between two-electron scattering states. This will mostly contribute to the final evolution of the electrons before detection. It will have strong terms, where the laser induces a change in momentum in one electron and the other electron remains unaffected. But it also includes the strongly correlated contribution, where the laser induces recollision and interaction between the two electrons. The dipole matrix element in the non-interacting case is given by

[TABLE]

Appendix E Two-Electron Integro-Differential Equations

The time variations of $a(t)$, $b(\mathbf{p},t)$, $c(\mathbf{p},\eta,t)$ and $d(\mathbf{p},\mathbf{p}^{\prime},t)$ follow from the TDSE for two electrons and read:

[TABLE]

Appendix F HHG from a quenched molecule: A toy model

In this appendix we analyze a toy 1D model of a quenched H${}_{2}^{+}$ molecule, or any similar diatomic molecule, stripped of one or more electrons, and treated in the SAE approximation. The remaining (or SAE) electron is after the rapid stripping process (say by an ultrashort XUV pulse) in the ground state of the model separable Hamiltonian. In the position representation, the electron feels a potential that consists of two Dirac-delta peaks at positions $\pm R/2$; in the momentum (wave vector) representation, that Hamiltonian reads

[TABLE]

where $R$ is the instantaneous distance between the nuclei. Setting $\frac{\hbar^{2}\kappa^{2}}{2m}=E_{0}=-I_{p}$, the ground-state wave function has the form

[TABLE]

By further denoting

[TABLE]

we obtain the equation that determines the energy (through $\kappa$) as

[TABLE]

where $\lambda=\pi\tilde{\lambda}m/\hbar^{2}$. The normalization of the ground-state wave function is therefore fixed by

[TABLE]

Similarly, for this simple model we can calculate the wave functions of the continuum states with the fixed outgoing momentum $p$, and the dipole moment $d(p)=\langle\Phi_{p}|x|\Phi_{0}\rangle$. Here we are interested in approximate qualitative/semi-quantitative description of HHG, so we will calculate the dipole moment $d(p)$ using plane waves, but employing the careful approach of Refs. Noslen1, ; Noslen2, ; Noslen3, ; Noslen4, ; Noslen-diss, . This amounts, essentially, to neglecting “unphysical” terms $\propto R$ that arise because of the non-orthogonality of plane waves and the ground state, $|\Phi_{0}\rangle$; from that approach, we get

[TABLE]

The motion of nuclei in the simplest case can be treated classically with Coulombic repulsion. Denoting the nuclear positions as $x_{1}$, $x_{2}$, with $x_{1}-x_{2}=r$ and $|r|=R$, and neglecting the electronic energy $\hbar^{2}\kappa(R)^{2}/2m$, we obtain

[TABLE]

where $Z$ is a “charge”, and $a$ is the smoothing parameter in the “Rochester” 1D version of the Coulomb potential. All of these parameters can and, in fact, should be adjusted to mimic real molecules.

In order to get some qualitative/semi-quantitative results, let us make some additional simplifying assumptions and observations:

• •

The femtosecond pulse is not too short, so its variation within a laser period can be neglected – we take it into account at the end integrating over the time-dependent $U_{p}(t)$.

• •

During the dissociation of our toy quenched molecule, $\kappa R$ changes from values $\ll 1$ to the ones $\gg 1$. Note that $R$ for the quenched molecule at the time $t=0$ is much smaller the the equilibrium size of the H${}_{2}^{+}$ molecule (compare Fig. 15 and discussion below). Correspondingly, the ionization potential $I_{p}$ changes from $I_{p}\simeq I_{p}(0)$ to $I_{p}\simeq I_{p}(\infty)=I_{p}(0)/4$. In other words, it undergoes a radical reduction by a factor of 4.

• •

The HHG spectrum cutoff for “local” contributions, i.e. electronic trajectories starting at $x_{i}$ and ending at $x_{i}$, changes thus from $3.2U_{p}+I_{p}(0)$ to $3.2U_{p}+I_{p}(0)/4$, which for large $I_{p}$ is quite significant. If dissociation is very fast, only the cutoff $3.2U_{p}+I_{p}(0)/4$ will be observed.

• •

Similarly, the efficiency and cutoff position for the “cross” contribution, i.e. electronic trajectories starting at $x_{1}$ and ending at $x_{2}$ (and vice versa) also changes significantly with $R$.

• •

At small $R$ (in the present terms $\kappa R<1$), the cutoffs for “local” and “cross” contributions are in distinguishable. Their efficiency is comparable, but they contribute with different phases, which leads to destructive interference (see Fig. 15, Refs. Noslen1, ; Noslen2, ; Noslen3, ; Noslen4, ; Noslen-diss, and references therein).

• •

As $R$ grows, we expect the cutoff of “cross” terms to grow, but their efficiency to decrease (see discussion below).

As commented above, from this figure we can clearly see that there is a deep minimum for both terms, and that it is located at the same position as for the total HHG spectra. This means that, for the case of the recombination on $\textbf{R}_{1}$ (dark green circle line), the electron wave packet ionized at $\textbf{R}_{1}$ interferes with the one coming from $\textbf{R}_{2}$ and the other way around. These minima are then generated by the destructive interference of such electron wavepackets. On the other hand, we have seen that our short-range non-local potential is unable to accurately reproduce the positions of the interference minima for some of the molecular orientation angles. We note, however, that these minima are typically washed out when an average over the molecular orientation is considered, a configuration that is commonly used in molecular HHG experiments.

Let us finally discuss more quantitatively the position of the HHG cutoff for the ‘local’ and ‘cross’ terms. For the local one, the quasi-classical theory predicts $3.17U_{P}+I_{p}$ (or $3.17U_{P}+xI_{p}$, with $x\simeq 1.3$ for small $I_{p}$, and $x\simeq 1$ for large $I_{p}$ Lewenstein1994 .

To this aim we repeat the saddle-point calculations from Ref. Lewenstein1994, for the cross trajectories, setting in this last example $e=1$, $c=1$ and $\omega=1$, to make the expressions essentially equal to the ones used in the 1994 paper. The equation for the momentum gives

[TABLE]

which leads to

[TABLE]

where we have fixed the trajectory to the one going from $-R/2$ to $R/2$. We can now calculate the gain of kinetic energy which should be equal to the emitted photon energy

[TABLE]

Taking $E(t)=\cos(t)$, $A(t)=-\sin(t)$, we get

[TABLE]

Obviously, to estimate the cutoff, we need to maximize the above expression first with respect to $t$, or better to say $t-\tau/2$, and then with respect to $\tau$. The first term gives exactly $3.17U_{p}$ as it does for the “local” trajectories. The second can give maxima which are much larger, depending on $R$, but they occur at completely “dephased” instants of $t-\tau/2$ – note that when the second term in maximal, $\cos(t-\tau/2)=\pm 1$, the first term vanishes. This is why, although the cutoff for the the “crossed” trajectories might be very long, their efficiency for large $R$ is not.

Appendix G RESI and EII using perturbative expressions for matrix elements

G.1 RESI

We have left the dipole matrix element $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\eta)$ unevaluated as we will deal with this next. We need to take into account interactions here and to do this we can use first order perturbation theory. In the one-electron case, to include a direct and ‘rescattering’-like contribution a perturbative series was used on the coefficients themselves. In this case in order to get a more explicit form of the electron-electron contribution we will use perturbation theory on the states themselves. We can split the states into a non-interacting part and a first order approximation to the interacting part as shown here

[TABLE]

where $\ket{\psi_{0}(\mathbf{p})}$ is given by the non-interacting states for $\ket{\mathbf{p},0}$ and

[TABLE]

where $\ket{\psi_{0}(\mathbf{p}^{\prime},\eta)}$ is given by the non-interacting states for $\ket{\mathbf{p},\eta}$. In the basis used in the perturbative expansion we have used excited states but also including the grounds state, i.e. $\eta=0$ as a basis to expand in. This is key to including transition from the ground state to excited states via electron interaction. Substituting this into the matrix element gives,

[TABLE]

The first part is just given by the non-interacting form that we previously calculated for this matrix element, while for the second we can insert Eq. (149)

[TABLE]

The last line shows the different pathways we have revealed by this expansion. Each term in the above equation relates to a different physical process: the first is where the ground state is excited by the laser while the scattering state is unaffected; the second term relates to the continuum electron being driven by the laser to interact with the ground-state electron, causing it to be excited – this is the most relevant one for RESI. These first two are shown in Fig. 13. The third term is similar to the second except there is an additional interaction with the core for the continuum electron, this could be from re-scattering after the RESI process has taken place. The final term excites the second electron in two steps through electron interaction and the laser via an intermediate state. An additional pathway can be revealed if the same expansion is applied to the right ket state of Eq. (150). This pathway is the same as the latter, except the order of the interactions is reversed and the laser excites the electron into the intermediate state from the ground state. We did not consider this, given these pathways will not contribute to the core of the RESI process.

Now we make the assumption that there will be no laser induced transitions between bound states and the matrix element simplifies it to

[TABLE]

This assumption also simplifies the integro-differential equations and we can insert the calculation we performed for $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\eta)$, this yields

[TABLE]

Now we select only terms that contribute to the final probability flow into double ionized states and ignore the flow in the opposite direction, as such we will assume $a(t)=1$ is unity throughout. This simplifies the integro-differential equations into a solvable form,

[TABLE]

G.2 EII

Here, we neglect the dipole matrix element for the recollision-excitation (RESI) contribution given by $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\eta)$, and instead include in a similar way the matrix element $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\mathbf{p}^{\prime\prime})$. We can proceed as before splitting the wavefunction into an non-interacting part and an interacting perturbation, where the perturbation is given by

[TABLE]

Then removing dipole transitions for the second electron from the bound state to the continuum to keep only the EI mechanism for the $\mathbf{g}(\mathbf{p},\mathbf{p}^{\prime},\mathbf{p}^{\prime\prime})$ matrix elements yields

[TABLE]

Now the integro-differential equation for $d(\mathbf{p},\mathbf{p}^{\prime},t)$ can be written out, this time it only depends on $b(\mathbf{p}^{\prime\prime},t^{\prime\prime})$ and is given by

[TABLE]

Then the solution can be written as

[TABLE]

The same equations for $b(\mathbf{p}^{\prime\prime},t^{\prime\prime})$ and $S_{d}\left(\mathbf{p},\mathbf{p}^{\prime},t^{\prime\prime},t\right)$ can be used as before. Additionally, both the EII and RESI mechanisms can be included in $d(\mathbf{p},\mathbf{p}^{\prime},t)$, and it will still be integrable and, as expected, will simply be equal to the sum of these two solutions.