Advanced momentum sampling and Maslov phases for a precise semiclassical model of strong-field ionization

Our starting point is the S-matrix amplitude [37]

$$\begin{equation} M_{\boldsymbol{p}_f} = -i \lim_{t \rightarrow \infty}\int^t_{-\infty} dt^{^{\prime}} {\langle\psi_{\boldsymbol{p}_f}\left(t\right)\vert}{U\left(t,t^{^{\prime}}\right)H_I\left(t^{^{\prime}}\right)}{\vert\psi_0\left(t^{^{\prime}}\right)\rangle}, \end{equation} \tag{ 1 }$$

describing the transition from initial system state $\vert \psi_0\rangle$, to a final scattering state $\vert \psi_{\boldsymbol{p}_f}\rangle$ with asymptotic momentum p _f . The Hamiltonian has been separated into a field-free part H₀ and time-dependent interaction part $H_I(t)$

$$\begin{equation} H\left(t\right) = H_0 + H_I\left(t\right), \end{equation} \tag{ 2 }$$

and $U(t,t_0)$ is the time-evolution operator for the entire system. We consider atomic systems within the single-active electron approximation (SAE), such that the field-free Hamiltonian takes the form

$$\begin{equation} H_0 = \frac{1}{2}\boldsymbol{p}^2 + V\left(\boldsymbol{r}\right), \end{equation} \tag{ 3 }$$

where p is the canonical momentum of the electron and $V(\boldsymbol{r})$ is the effective SAE potential. For all targets considered in this text, the potential is of the form

$$\begin{equation} V\left(\boldsymbol{r}\right) = -\frac{Z + f\left(r\right)}{r} \end{equation} \tag{ 4 }$$

where Z is the charge of the parent ion and

$$\begin{equation} f\left(r\right) = a_1\text{e}^{-a_2 r} + a_3 r\text{e}^{-a_4 r} + a_5\text{e}^{-a_6 r}, \end{equation} \tag{ 5 }$$

with coefficients from [50].

We work exclusively within the electric dipole approximation, for which the interaction Hamiltonian in the length gauge is given by

$$\begin{equation} H_I\left(t\right) = \boldsymbol{E}\left(t\right)\cdot\boldsymbol{r}, \end{equation} \tag{ 6 }$$

where the electric field is related to the vector potential $\boldsymbol{A}(t)$ as $\boldsymbol{E}(t) = - \partial_t\boldsymbol{A}(t)$.

To proceed, we approximate the final continuum eigenstate as a field-dressed plane wave $\vert \psi_{\boldsymbol{p}_f}(t)\rangle \simeq \vert \boldsymbol{p}_f + \boldsymbol{A}(t)\rangle \equiv \vert \tilde{\boldsymbol{p}}_f(t)\rangle$, and insert the identity

$$\begin{equation} \unicode{x1D7D9} = \int d\tilde{\boldsymbol{p}}_0\left(t^{^{\prime}}\right) \vert \tilde{\boldsymbol{p}}_0\left(t^{^{\prime}}\right)\rangle\langle \tilde{\boldsymbol{p}}_0\left(t^{^{\prime}}\right)\vert \end{equation} \tag{ 7 }$$

in equation (1)

$$\begin{equation} \begin{aligned} M_{\boldsymbol{p}_f} = {} -i \lim_{t \rightarrow \infty}\int_{-\infty}^t \mathrm{d}{t^{^{\prime}}} \int \mathrm{d}^3{\tilde{\boldsymbol{p}}}_0\!\left(t^{^{\prime}}\right) e^{iI_p t^{^{\prime}}} d\left(\tilde{\boldsymbol{p}}_0,t^{^{\prime}}\right) \langle{\tilde{\boldsymbol{p}}_f\left(t\right)}\vert{U\left(t,t^{^{\prime}}\right)\vert}{\tilde{\boldsymbol{p}}_0\left(t^{^{\prime}}\right)\rangle}, \end{aligned} \end{equation} \tag{ 8 }$$

with I_p the ionization potential of the target, and the ionization matrix element

$$\begin{equation} d\left(\tilde{\boldsymbol{p}}_0,t\right) = \langle \tilde{\boldsymbol{p}}_0\left(t\right)\vert H_I\left(t\right)\vert \psi_0\rangle, \end{equation} \tag{ 9 }$$

where $\vert \psi_0\rangle$ is the time-independent initial state. In this way, we capture the time-evolution following ionization using the propagator $\langle{\tilde{\boldsymbol{p}}_f(t)\vert}{U(t,t^{^{\prime}})\vert}{\tilde{\boldsymbol{p}}_0(t^{^{\prime}})\rangle}$ from an initial dressed momentum state $\vert \tilde{\boldsymbol{p}}_0(t^{^{\prime}})\rangle$ to a final dressed momentum state $\vert \tilde{\boldsymbol{p}}_f(t)\rangle$. The semiclassical nature of the CQSFA stems from approximating this propagator by a semiclassical expression.

In the CQSFA, the propagator in equation (8) is written in terms of path integrals, to which semiclassical approximations can be applied. This is done by time slicing, where $U(t,t^{^{\prime}})$ is split into N + 1 operators using its composition property, each propagating a time Δt. Between these operators the identity is resolved in dressed momentum states and position states, allowing one to write [51]

$$\begin{equation*} \langle{\tilde{\boldsymbol{p}}_f\left(t\right)\vert}{U\left(t,t^{^{\prime}}\right)\vert}{\tilde{\boldsymbol{p}}_0\left(t^{^{\prime}}\right)\rangle} \! = \!\prod_{n = 1}^{N} \!\left[\int \!\mathrm{d}^{3} \tilde{\boldsymbol{p}}_n\right] \prod_{n = 1}^{N+1} \!\left[\int \!\frac{\mathrm{d}^{3} \boldsymbol{r}_n}{\left(2\pi\right)^3}\right] e^{i \mathcal{A}_N}, \end{equation*}$$

where

$$\begin{equation} \mathcal{A}_N = -\sum_{n = 1}^{N+1} \left[\boldsymbol{r}_n \cdot \left(\tilde{\boldsymbol{p}}_n - \tilde{\boldsymbol{p}}_{n-1}\right) + \Delta t H\left(\tilde{\boldsymbol{p}}_n, \boldsymbol{r}_n, t_n\right)\right]. \end{equation} \tag{ 10 }$$

The limit $N \rightarrow \infty$, is now taken, in order to turn the product of integrals into path integrals.

The above formulation used the transition matrix element in a momentum representation, propagating from an initial momentum $\tilde{\boldsymbol{p}}_0$ to a final $\tilde{\boldsymbol{p}}_f$. It is also possible to work in a 'mixed' representation, where propagation instead is from an initial position r ₀ to a final momentum [38]. The two forms of the matrix elements can be related as

$$\begin{equation} \langle{\tilde{\boldsymbol{p}}_f\vert}{U\left(t,t^{^{\prime}}\right)\vert}{{\boldsymbol{r}}_0\rangle} = \int \mathrm{d}^{3} \tilde{\boldsymbol{p}}_0 e^{-i\tilde{\boldsymbol{p}}_0 \cdot \boldsymbol{r}_0} \langle{\tilde{\boldsymbol{p}}_f}\vert{U\left(t,t^{^{\prime}}\right)\vert}{\tilde{\boldsymbol{p}}_0}\rangle, \end{equation} \tag{ 11 }$$

which, together with equation (10), shows that the 'sliced' action in a mixed representation includes an additional boundary term

$$\begin{equation} \mathcal{S}_N = \mathcal{A}_N- \tilde{\boldsymbol{p}}_0 \cdot \boldsymbol{r}_0. \end{equation} \tag{ 12 }$$

Returning to equation (8) we see that if we pull out an exponential dependence of d, equation (9),

$$\begin{equation} d\left(\tilde{\boldsymbol{p}}_0,t\right) = e^{-i\tilde{\boldsymbol{p}}_0 \cdot \boldsymbol{r}_0} \tilde{d}\left(\tilde{\boldsymbol{p}}_0,\boldsymbol{r}_0,t\right), \end{equation} \tag{ 13 }$$

we can write the transition amplitude in the mixed representation

$$\begin{equation} \begin{split} M_{\boldsymbol{p}_f} = {} -i \lim_{t \rightarrow \infty} \int_{-\infty}^t \mathrm{d} t^{^{\prime}} \prod_{n = 0}^{N} \left[\int \mathrm{d}^{3} \tilde{\boldsymbol{p}}_n\right] \prod_{n = 1}^{N+1} \left[\int \frac{\mathrm{d}^{3} \boldsymbol{r}_n}{\left(2\pi\right)^3}\right] e^{i \mathcal{S}_N} e^{iI_p t^{^{\prime}}} \tilde{d}\left(\tilde{\boldsymbol{p}}_0, \boldsymbol{r}_0, t^{^{\prime}}\right), \end{split} \end{equation} \tag{ 14 }$$

now including the $\tilde{\boldsymbol{p}}_0$ integral in the product. In the limit $N\rightarrow \infty$, we obtain the path integrals

$$\begin{equation} M_{\boldsymbol{p}_f} = -i \lim_{t \rightarrow \infty} \int_{-\infty}^t \mathrm{d} t^{^{\prime}} \int^{\mathbf{p}_f} \mathcal{D}\tilde{\boldsymbol{p}} \int_{\mathbf{r}_0}\frac{\mathcal{D}\boldsymbol{r}}{\left(2\pi\right)^3} e^{iS} \tilde{d}\left(\tilde{\boldsymbol{p}}_0, \boldsymbol{r}_0, t^{^{\prime}}\right), \end{equation} \tag{ 15 }$$

with the mixed representation action

$$\begin{equation} S\left(\mathbf{p}, \mathbf{r}, t^{^{\prime}}\right) = I_pt^{^{\prime}} -\int_{t^{^{\prime}}}^{\infty} \mathrm{d} \tau \left[\dot{\boldsymbol{p}} \cdot \boldsymbol{r}\left(\tau\right) + H\left(\boldsymbol{p}, \boldsymbol{r}, \tau\right)\right] - \tilde{\boldsymbol{p}}_0 \cdot \boldsymbol{r}_0 \end{equation} \tag{ 16 }$$

and Hamiltonian

$$\begin{equation} H\left(\boldsymbol{p}, \boldsymbol{r}, \tau\right) = \frac{1}{2}\left[\boldsymbol{p}\left(\tau\right) + \boldsymbol{A}\left(\tau\right)\right]^2 + V\left[\boldsymbol{r}\left(\tau\right)\right].\ \end{equation} \tag{ 17 }$$

The fact that we have an action corresponding to a mixed representation is the key point of this section. The effect of this is that when we apply the saddle-point approximation (SPA) to the path integrals, in order to obtain a semiclassical description, we obtain a stability matrix of the form $\partial \boldsymbol{p}(t) / \partial \boldsymbol{p}_0$, instead of the $\partial \boldsymbol{p}(t) / \partial \boldsymbol{r}_0$ that is obtained in a momentum representation [38, 52].

Next, we perform the SPA on the time and path integrals [37, 51], thereby obtaining the desired transition amplitude

$$\begin{equation} M_{\boldsymbol{p}_f} \propto -i \lim_{t \rightarrow \infty} \sum_{s} \left[\sqrt{\frac{2\pi i}{\partial^2 S\left(\boldsymbol{p}_s, \boldsymbol{r}_s, t_s\right)/\partial t_s^2}} \tilde{d}\left(\tilde{\boldsymbol{p}}_0, \boldsymbol{r}_0, t_s\right) \frac{1}{\sqrt{|J_s\left(t\right)|}} e^{i S\left({\boldsymbol{p}}_s, \boldsymbol{r}_s, t_s\right)} e^{-i\frac{\pi \nu_s}{2}} \right], \end{equation} \tag{ 18 }$$

where ν_s is the Maslov index, to be properly introduced in the following section. The sum over s is over all paths connecting the final and initial momentum and fulfilling the classical equations of motion, which follow from stationary action in time, position and momentum [37]

$$\begin{align}\frac{1}{2}\left[\boldsymbol{p}_0 + \boldsymbol{A}\left(t_s\right)\right]^2 = -I_\text{p}, \end{align} \tag{ 19 }$$

$$\begin{align}\dot{\boldsymbol{r}}_s\left(\tau\right) = \boldsymbol{p}_s\left(\tau\right) + \boldsymbol{A}\left(\tau\right), \end{align} \tag{ 20 }$$

$$\begin{align}\dot{\boldsymbol{p}}_s\left(\tau\right) = -\nabla_{\boldsymbol{r}_s} V\left[\boldsymbol{r}_s\left(\tau\right)\right], \end{align} \tag{ 21 }$$

respectively. The subscript s indicates that a quantity is connected to the classical trajectory equations (20) and (21) given by initial momentum and time, p ₀ and t_s , related through equation (19). In this way s functions as a collective index $s = (\boldsymbol{p}_0, t_s)$ of the parameters uniquely identifying a given trajectory (equation (19) can have several solutions for a given p ₀).

Equations (19)–(21) are referred to as the saddle-point equations (SPEs). The t_s is usually related to the time of ionization, with the first SPE, equation (19), expressing energy conservation at ionization. For a vanishing potential $\boldsymbol{p}_s (\tau)$ becomes constant, and equation (19) reduces to the SPE from the lowest order SFA [37, 39], providing a link between the ionization time and final momentum of the electron. This is in agreement with the fact that the potential is neglected in the lowest order SFA. When considering SFA trajectories in the following, we are thus considering solutions to equations (19) and (20) for a constant momentum value, $\boldsymbol{p}_0 = \boldsymbol{p}_f$, with p _f the final momentum value of the photoelectron.

J_s in equation (18) is the Jacobian determinant in the mixed representation

$$\begin{equation} J_s\left(\tau\right) = \text{det}\left(\frac{\partial\boldsymbol{p}_s\left(\tau\right)}{\partial\boldsymbol{p}_0}\right), \end{equation} \tag{ 22 }$$

in general evaluated at some time τ along the classical trajectory, equations (20) and (21). Note that in equation (18) $J_s(t)$ is to be evaluated at asymptotic final time t, where the field is over and the electron far from the parent ion. The Jacobian matrix is also referred to as the stability matrix, since it describes the trajectory's sensitivity to the initial momentum value. At points along a trajectory where $J_s(\tau) = 0$, so-called focal points, there is an insensitivity to changes in initial conditions of the trajectory in some direction. Effectively, this means that several trajectories (with different initial conditions) will pass through the same point, and we thus have a bunching of trajectories at the focal points. If the final momentum of a trajectory happens to be in a focal point, this bunching effect will be apparent in the final transition amplitude equation (18) and cause a breakdown of the semiclassical approximation in those points [53]. Surfaces in the final momentum distribution containing these focal points are known as 'caustics' and are clearly visible in calculated photoelectron momentum distributions (PMDs) as a line of high transition amplitude.

The ionization times in equation (18) are generally complex $t_s = t_s^{\text{Re}} + \text{i} t_s^{\text{Im}}$, so the integral of the action, equation (16), becomes a complex contour integral. We pick a contour composed of two parts; first we integrate from t_s to the real time axis, parallel to the imaginary time axis. This part is interpreted to correspond to the tunneling of the electron out of the potential. During the integration of this first contour, the momentum is assumed to be constant. Thereafter, the action is integrated from $t_s^{\text{Re}}$ to infinity, corresponding to propagation in real time. In total, the action, equation (16), is thus split into tunneling and propagation parts [54]. In the tunneling part of the action, an integral over the potential appears. In practice this will be replaced by a Coulomb correction factor [47, 55, 56]:

$$\begin{equation} \text{exp}\left[-\int_{t_s}^{\mathrm{Re}\left(t_s\right)} \mathrm{d}{\tau}V\left[\boldsymbol{r}_s\left(\tau\right)\right]\right]\longrightarrow \left( \frac{2\kappa^2}{|\boldsymbol{E}\left(t_s\right)|}\right)^{\frac{2}{\kappa}}, \end{equation} \tag{ 23 }$$

where $\kappa = \sqrt{2I_p}$. The reasoning for this substitution comes from Perelomov–Popov–Terent'ev (PPT) theory [55, 56], which is widely used to explain the ionization rates of atoms. Specifically, by regularizing the integral and expanding in the path to second order, one obtains a finite integral, which resembles the expression that has been substituted. A similar result was also recently derived in a CQSFA setting in reference [40]. In order to propagate $t\rightarrow\infty$, after the pulse is over and the electron is sufficiently far from the residual ion, we analytically propagate the equations of motion and compute the action, using the same procedure as outlined in [38].

In order to calculate the ionization matrix element, equation (9), we use the asymptotic expansion of a Coulomb wave function

$$\begin{equation} \psi_0\left(\boldsymbol{r}\right) \xrightarrow{r\rightarrow\infty} \sum_{\ell, m} C_{\ell m} r^{Z/\kappa -1} e^{-\kappa r} Y^{m}_\ell \left(\theta, \phi\right), \end{equation} \tag{ 24 }$$

where $Y^m_\ell(\theta,\phi)$ is a spherical harmonic, $\kappa = \sqrt{2I_p}$ and $C_{\ell m}$ an expansion coefficient. This expansion has previously been shown to work well in the length gauge [57], also when comparing with experimental data for molecules [58]. Following the regularization procedure outlined in [57], we find the matrix element evaluated in the SFA time saddle points, t_s , to be

$$\begin{equation} d\left(\boldsymbol{p}, t_s\right) \propto E\left(t_s\right) \left(S^{^{\prime\prime}}\left(t_s\right)\right)^{-\left(1/2\kappa +1\right)} \sum_{\ell, m} C_{\ell, m} \left(\frac{\tilde{\boldsymbol{p}}\left(t_s\right)}{i\kappa}\right)^\ell \left(\alpha_- Y^m_{\ell-1}\left(\Omega_{\tilde{\boldsymbol{p}}}\right) \frac{\kappa^2}{\tilde{\boldsymbol{p}}\left(t_s\right)} - \alpha_+ Y^m_{\ell + 1}\left(\Omega_{\tilde{\boldsymbol{p}}}\right) \tilde{\boldsymbol{p}}\left(t_s\right)\right) \end{equation} \tag{ 25 }$$

with the abbreviations

$$\begin{align} \alpha_+ & = \left[\frac{\left(\ell-m+1\right)\left(\ell+m+1\right)}{\left(2\ell+1\right)\left(2\ell+3\right)}\right]^{1/2}, \end{align} \tag{ 26 }$$

$$\begin{align} \alpha_- & = \left[\frac{\left(\ell-m\right)\left(\ell+m\right)}{\left(2\ell-1\right)\left(2\ell+1\right)}\right]^{1/2}. \end{align} \tag{ 27 }$$

The constant of proportionality is neglected, since we renormalize the simulation results. For the simulations performed here, we only employ one spherical harmonic in our initial state expansion, choosing $C_{0,0} = 1$ for an s-state (hydrogen and helium) and $C_{1,0} = 1$ for a p-state (argon), aligned along the laser polarization direction.

Finally, we assume that the electron starts from $\boldsymbol{r}_s(t_s) = \boldsymbol{0}$ at the start of the tunneling, in which case we can simply set the phase factor $\boldsymbol{r}_0 \cdot \tilde{\boldsymbol{p}_0}$ in equation (13) to zero, and neglect the boundary term in the action, i.e. neglect the last term in equation (16).

The full CQSFA transition amplitude, equation (18), contains the phase-factor $e^{-i\frac{\pi\nu_s}{2}}$ known as the Maslov phase with ν_s the Maslov index. This phase stems from the phase of the Jacobian determinant, J_s , hence the use of $|J_s|$ in equation (18). It is worth noting that previous iterations of the CQSFA used $1/\sqrt{J_s(t)}$ directly in the expression for the final transition amplitude (see, e.g. [37]), and did not explicitly write out the Maslov phase. However, as noted in [51], evaluating the transition amplitude with $1/\sqrt{J_s(t)}$ directly, for a given final time t, is only correct as long as $1/\sqrt{J_s(t)}$ is regular, i.e. as long as the corresponding classical trajectory does not pass a focal point. Earlier CQSFA iterations used fewer trajectories, with only a small number passing through focal points, and it was assumed the effect would be small. Recently, the Maslov phase was included in a semiclassical strong-field model [47], similar to the CQSFA. In this section, we elaborate on this procedure, based on references [53, 59].

In the focal points, for which $J_s(\tau) = 0$, the Maslov phase will change discontinuously, corresponding to an integer change in the Maslov index. This behavior has an analogy in optics [47]: When focusing a Gaussian beam, it asymptotically obtains a π phase-shift when passing through its focus, known as the Gouy phase shift. Similarly, the ionic potential may act as a focusing lens and focus a manifold of trajectories with infinitesimally small deviations in initial conditions to a single point. The Maslov index is hence a time-dependent function, $\nu_s = \nu_s(\tau)$, and the contributions of each focal point along the trajectory must be calculated, as outlined below.

To calculate $J_s(\tau)$, we require the equations of motion of the stability matrices, which are in the form of a Jacobi initial-value problem. The derivation of the equations of motion and their initial conditions are detailed in appendix A. We find that the time-evolution of the stability matrices in the mixed representation are governed by the following coupled matrix-differential equations

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d \tau}} \frac{\partial\boldsymbol{p}_s\left(\tau\right)}{\partial\boldsymbol{p}_0} & = -\mathbf{H}\left[V\left(\boldsymbol{r}_s\left(\tau\right)\right)\right]\frac{\partial\boldsymbol{r}_s\left(\tau\right)}{\partial\boldsymbol{p}_0} \end{align} \tag{ 28 }$$

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d \tau}}\frac{\partial\boldsymbol{r}_s\left(\tau\right)}{\partial\boldsymbol{p}_0} & = \frac{\partial\boldsymbol{p}_s\left(\tau\right)}{\partial\boldsymbol{p}_0}, \end{align} \tag{ 29 }$$

where $\mathbf{H}[V(\boldsymbol{r}_s(\tau))]$ denotes the Hessian of the potential, which, as noted, must be evaluated along the classical trajectories. The initial conditions of the stability matrices read

$$\begin{equation} {\frac{\partial\boldsymbol{p}_s\left(\tau\right)}{\partial\boldsymbol{p}_0}}\Bigg|_{\tau = \mathrm{Re}\left(t_s\right)} = \unicode{x1D7D9}, \qquad {\frac{\partial \boldsymbol{r}_s\left(\tau\right)}{\partial \boldsymbol{p}_0}}\Bigg|_{\tau = \mathrm{Re}\left(t_s\right)} = \boldsymbol{0}. \end{equation} \tag{ 30 }$$

For equation (4) the Hessian reads

$$\begin{equation} \mathbf{H}\left[V\left(\boldsymbol{r}_s\right)\right] = \frac{1}{r_s}\frac{\mathrm{d}V\left(r_s\right)}{\mathrm{d} r_s}\unicode{x1D7D9} - \frac{1}{r_s^3}\left[\frac{\mathrm{d}^2f\left(r_s\right)}{\mathrm{d} r_s^2} + 3\frac{\mathrm{d}V\left(r_s\right)}{\mathrm{d}r_s}\right] \boldsymbol{r}_s \otimes\boldsymbol{r}_s, \end{equation} \tag{ 31 }$$

where $\unicode{x1D7D9}$ denotes the $3\times3$ identity matrix and ⊗ the tensor product. Propagation of equations (28) and (29) along with the SPEs, equations (19)–(21), then allows determination of $J_s(\tau)$ through equation (22).

With $J_s(\tau)$ obtained for any time τ along the trajectory, we can determine the times τ_i for which $J_s(\tau_i) = 0$, that is, the focal points. Since ν_s only changes in these points, its value at the final time t can be written as the sum of these changes at the focal times

$$\begin{equation} \nu_s\left(t\right) = \nu_{s,0} + \sum_i\Delta\nu_s\left(\tau_i\right), \end{equation} \tag{ 32 }$$

where $\nu_{s,0}$ is the initial value at the start of the trajectory, $\Delta \nu_s(\tau_i)$ the change at a given focal point (to be determined below), and the sum runs over all focal points $\tau_i \in [\mathrm{Re}(t_s), t]$. The initial value $\nu_{s,0}$ is only nonzero if the trajectory happens to start upon a focal point [59]. When $\nu_s(t)$ is to be used in the final transition amplitude, equation (18), t is again to be seen as an asymptotic time, $t\rightarrow\infty$, but practically large enough that the laser field is over, and the photoelectron is far from its parent ion. At the focal points, the determinant of the stability matrix will be zero, $J_s(\tau_i) = 0$, and consequently $\partial \boldsymbol{p}_s(\tau_i) / \partial \boldsymbol{p}_0$ will have some number, m, of linearly independent eigenvectors with corresponding eigenvalues of 0, denoted $\boldsymbol{d}^{(k)}$, that is, m solutions to the eigenvalue problem

$$\begin{equation} \frac{\partial\boldsymbol{p}_s\left(\tau_i\right)}{\partial\boldsymbol{p}_0} \boldsymbol{d}^{\left(k\right)} = \lambda_k \boldsymbol{d}^{\left(k\right)}, \end{equation} \tag{ 33 }$$

with $\lambda_k = 0$ for $k = 1,\ldots ,m$. Given the $\boldsymbol{d}^{(k)}$ eigenvectors at these times, the change in Maslov index at a focal point $\Delta\nu_s(\tau_i)$ is given by the following formula

$$\begin{equation} \Delta\nu_s\left(\tau_i\right) = m - 1 + \text{sgn}\text{det}\boldsymbol{\Pi}\left(\tau_i\right), \end{equation} \tag{ 34 }$$

where $\boldsymbol{\Pi}(\tau_i)$ is a matrix with the following components

$$\begin{equation} \Pi_{nm}\left(\tau_i\right) = \boldsymbol{u}^{\left(n\right)}\left(\tau_i\right)\cdot \mathbf{H}\left[V\left(\boldsymbol{r}_s\left(\tau_i\right)\right)\right]\boldsymbol{u}^{\left(m\right)}\left(\tau_i\right), \end{equation} \tag{ 35 }$$

where we have defined

$$\begin{equation} \boldsymbol{u}^{\left(k\right)}\left(\tau_i\right) = \frac{\partial\boldsymbol{r}_s\left(\tau_i\right)}{\partial \boldsymbol{p}_0} \boldsymbol{d}^{\left(k\right)}. \end{equation} \tag{ 36 }$$

From our experience, setting $\nu_{s,0} = 0$ for all cases will not notably change the result due to the occurrence of focal points at the start of the propagation being extremely rare.

In order to calculate the transition amplitude, equation (18), we sum over all saddle-point solutions fulfilling equations (19)–(21). Given that this transition amplitude is formulated via a mixed representation path integral, this means we have an initial condition $\boldsymbol{r}_s(t_s) = 0$ and a boundary condition $\boldsymbol{p}_s(\tau\rightarrow\infty)\rightarrow\boldsymbol{p}_f$. Thus, to find solutions, one must find all the initial momenta p ₀ that lead to a particular final momentum p _f . This is known as the inversion problem. Due to the difficulty of finding all of these solutions, many models take an alternative approach, where many starting trajectories are forward propagated and the final momentum binned, so that trajectories that end in the same final bin are summed coherently. To resolve interference, this requires very small bins and typically in excess of a billion trajectories. Additionally, it has also been shown in [60] that these 'forward' approaches do not yield the correct sampling weight in terms of the Jacobian J_s , leading to $1/|J_s|$ instead of the correct $1/\sqrt{|J_s|}$ computed by inverse approaches (see equation (18)).

In earlier iterations of the CQSFA [39, 41, 42], not all solutions to the SPEs were included. This method assumed that there were only 4 solution per laser cycle and used an iterative procedure to find them. Later, additional solutions were identified [42] but were not incorporated into the CQSFA. The approach taken here, is to use a Monte-Carlo-style sampling method, to find all the initial momenta that lead to any particular final momentum, see figure 1. This works by generating a partially random guess for an initial momentum $\boldsymbol{p}_{0,g}$, which (propagating via the SPEs) leads to a particular final momentum $\boldsymbol{p}_{f,g}(\boldsymbol{p}_{0,g})$. Then, in analogy to the 'forward' method, we bin this solution, associating it with the closest final momentum grid point $\boldsymbol{p}_{f,j}$ to $\boldsymbol{p}_{f,g}(\boldsymbol{p}_{0,g})$, where j denotes the index of the grid. However, now this method differs from the 'forward' method, as we change the initial momentum to make the resulting final momentum approach that of the grid point. More precisely, the initial guess momentum $\boldsymbol{p}_{0,g}$ is used as a starting point to find a slightly different initial momentum $\boldsymbol{p}_{0,g}^{^{\prime}}$ so that the final guess momentum and the grid point momentum are essentially equal, by solving $|\boldsymbol{p}_{f,g}(\boldsymbol{p}_{0,g}^{^{\prime}})-\boldsymbol{p}_{f,j}|\lt\epsilon$, where ε is taken to be a very small value of the order $\epsilon\sim 10^{-12}$. When sampling several initial guesses, care must be taken that any new initial guesses $\boldsymbol{p}_{0,g}$ do not lead to a solution that has already been binned.

Zoom In Zoom Out Reset image size

If enough initial momenta are provided by a random number generator over a large enough range of p ₀, all unique solutions can be found in this way. However, this algorithm can be very inefficient if solutions are clustered very densely in small areas of p ₀ momentum space, which is the case for rescattered trajectories. An example of the set of all initial momentum obtained for a Helium simulation is given in figure 1. The rescattered trajectories have a negative initial momentum coordinate perpendicular to the laser field polarization $p_{0\perp}\lt0$, which form densely clustered points that are more difficult to sample. To address this, we implemented an algorithm that can detect narrow regions of densely clustered solutions, and then sampling can be more concentrated in these regions. These regions can be detected after performing an initial sampling over the full region depicted in figure 1. Then a grid is placed over the region $p_{0\perp}\lt0$, which is coarse enough so that if a dense cluster of points intersects a box in the grid, the initial sampling will have found at least one solution in this box. The clusters can then be intensively searched by restricting guesses only to boxes in the grid with at least one solution already present.

A similar approach was recently implemented in [40], however, in this case the sampling was done with a Gaussian distribution. This is efficient for the direct trajectories that do not undergo rescattering, but less so for the rescattered trajectories. Rodriguez et al [40] also did not consider a pulse envelope or the Maslov phases that we consider here.

In previous versions of the CQSFA the division of the quantum orbits into 4 types [34] was explicitly used in the simulations [37, 42]. For a laser linearly polarized along the z direction, the classification compares the tunnel exit

$$\begin{equation} z_{0} = \int_{t_s}^{\mathrm{Re}\left(t_s\right)} A\left(\tau\right) d\tau, \end{equation} \tag{ 37 }$$

and the final momentum along the polarization direction, $p_{f\parallel}$, to the final and initial perpendicular momentum $p_{f\perp}$ and $p_{0\perp}$, through the signs of the products $r_{0\parallel} p_{f\parallel}$ and $p_{f\perp} p_{0\perp}$, as summarized in the legend of figure 2. We briefly describe the 'standard' trajectories of each orbit type, also shown on the same figure. Orbit type 1 ($r_{0\parallel} p_{f\parallel} \gt 0$, $p_{f\perp} p_{0\perp}\gt0$) corresponds to direct trajectories, not revisiting the ion core. Type 2 ($r_{0\parallel} p_{f\parallel} \lt 0$, $p_{f\perp} p_{0\perp}\gt0$) describes trajectories driven past the core by the field, but without much interaction. Type 3 ($r_{0\parallel} p_{f\parallel} \lt 0$, $p_{f\perp} p_{0\perp}\lt0$) corresponds to trajectories performing a field-driven forward scattering of the ion core. Last, type 4 ($r_{0\parallel} p_{f\parallel} \gt 0$, $p_{f\perp} p_{0\perp}\lt0$) describes trajectories performing a field-driven backwards scattering of the core. Orbit types 1 and 2 only interact weakly with the ion potential, and can thus be related to the trajectories found in the lowest order SFA (see the text following equation (21)). On the other hand, orbit type 3 and 4 require the presence of the Coulomb potential. Since the semiclassical trajectories here are restricted to a plane we use a 2D definition of the orbit types, but a simple generalization to 3D trajectories is possible [61].

Zoom In Zoom Out Reset image size

We note that other trajectories than these 'standard' versions exist, such as multi-pass trajectories, having multiple close encounters with the ion core, and directly re-colliding trajectories, where the ionization happens at low fields and the scattering is primarily driven by Coulomb forces. While these were neglected in previous CQSFA versions [42], the current version of the CQSFA will include all trajectory dynamics, assuming sufficient sampling.

In figure 3, we compare CQSFA and TDSE simulations performed for hydrogen, helium, and argon. The simulation parameters for helium are the same as used in [47], in order to enable direct comparison. For all targets, we see a very good agreement between the TDSE and the CQSFA with the Maslov phase included, figures 3(a)–(c), in the lower energy region before the presence of caustics. As mentioned in the section following equation (22), the caustics are surfaces containing focal points for the classical trajectories at the final asymptotic time. Because of the breakdown of the semiclassical model in these points, the caustics are generally visible in the PMDs as sharp lines of discontinuity in the signal, e.g. in figure 3(a) at $\boldsymbol{p} \sim (0, -0.8)$ a.u.

Zoom In Zoom Out Reset image size

When we include the Maslov phase into the CQSFA, it is evident that especially the location of the interference fringes along $p_\parallel$ at $p_\perp\sim0$ are very well-reproduced. The fact that the Maslov phase is indeed necessary to obtain this agreement is illustrated in figures 3(d)–(f). Here, the CQSFA with and without the Maslov phase is shown on each half, and a shift of interference fringes is seen along the $p_\perp = 0$ axis, see panel (e) box 1 for He. To allow a detailed comparison of these fringes, an enlarged view of this part of the PMD for all the panels in figure 3 is provided in figure 8 in appendix B.

Another important structure affected by the inclusion of the Maslov phase is the fringes at $p_\parallel \gt 0$ almost parallel to the $p_\parallel$ axis, commonly known as the 'spider-legs' [18]. This interference pattern is found for a wide range of targets, and is visible in all the PMDs in figure 3, most clearly for He. As seen on panels (d)–(f) the Maslov phase results in a 'lift' of the central fringe of the spider (panel (e) box 2 for He), increasing its width to a size in much better agreement with the TDSE results. This effect is further illustrated in figures 4(b), (d) and (f), showing the signal along a slice of the PMDs of figure 3 for $p_\perp = 0.1$ a.u. While the CQSFA results both seem to agree with the TDSE simulations for low $p_\parallel$, the inclusion of the Maslov phase clearly creates a difference for $p_\parallel \gt 0.5$ a.u., increasing agreement with the TDSE results significantly. In figures 4(a), (c) and (e), a slice of the signal along $p_\parallel = 0.4$ a.u. also illustrates how the location of the other spider legs are affected by the Maslov phase. For all targets the agreement is better with the inclusion of the Maslov phase, perhaps most clearly for hydrogen, where quite large shifts in the interference minimum can be seen (compare e.g. the gray dashed line at $p_\perp = 0.25$ a.u. to the blue at $p_\perp = 0.30$ a.u. in figure 4(a)). This illustrates the necessity to take the Maslov phase into account, if one wants to perform analysis based on the spider structure.

Zoom In Zoom Out Reset image size

Still, the CQSFA and TDSE results do not agree completely. One particular region where the results seem to disagree, is for high values along the $p_\parallel$ axis, where the CQSFA results drop off too fast. This is clearly illustrated for both hydrogen and helium in panels (b) and (d) in figure 4, and can also be seen in the PMDs of figure 3. The disagreement seems to be worst for hydrogen. The reason for this could be the fact that the hydrogen simulation is performed for the lowest intensity of all the simulations, and in connection also the simulation with the highest Keldysh parameter ($\gamma_\text{K} = 0.94$). To test the range of validity of the CQSFA, we performed several simulations with values of Keldysh parameters in the range $[0.7, 1.1]$. In general, testing suggests that the CQSFA becomes better for longer wavelengths [61], higher intensities and lower Keldysh parameters, a trend in conformity with the well performing helium simulation which has the highest intensity and lowest Keldysh parameter ($\gamma_\text{K} = 0.72$). The CQSFA also performs well for intensities corresponding to over-the-barrier dynamics.

Another source of discrepancy between the CQSFA and TDSE results is, as mentioned, the caustics. The effect of those are worse for the hydrogen and argon simulations, as seen in panels (a) and (c) of figure 3, where they cut off most of the structure at high energies, an effect also visible as a sharp drop-off on panel(a) of figure 4. Since the caustics represent classical turning points, their position can be pushed outwards from the center by increasing, for example, the intensity or the number of cycles of the pulse.

The results also show that the argon simulation is less converged than the hydrogen and helium simulations, see figures 4(e) and (f), where the blue and gray lines showing the Ar CQSFA results are clearly less smooth than the corresponding curves for He and H, figures 4(a)–(d). The reason for this is that the effective potential used for argon gives rise to many additional (and seemingly irrelevant) solutions to the SPEs, meaning significantly more sampling of orbits is required to reach convergence. It is therefore worth noting, that the inclusion of f(r) in equation (4) is not critical for the current simulations. Just using the plain Coulomb potential (but the correct I_p for the targets) yields results almost identical to those in figure 3, with some small difference in the signal value for the high energy rings, while convergence is reached with less sampling. It could be, however, that the effective potential plays a larger role for other systems or field parameters, such as fields with more cycles where higher energies are reached.

Overall, the large degree of agreement between the CQSFA results and the TDSE simulations shown in figure 3 illustrates the accuracy and validity of the semiclassical model.

The high level of accuracy of the CQSFA adds a new level of credibility to orbit-based analyses. Given how well the semiclassical model reproduces the TDSE results, it seems likely that the sampled trajectories and their properties could shed light onto key aspects of the underlying physics. Studying the different types of trajectories found in the new model, along with their impact on the PMDs, could thus lead to an improved understanding of the ATI process.

To investigate the different types of trajectories and their effects on the transition amplitude, we plot in figure 5(a1) the ionization time as a function of final parallel momentum for orbits sampled in the helium simulation shown in figure 3(b), at $p_\perp = 0.25$ a.u. (the CQSFA result in figure 3(b) is shown for $p_\perp\lt0$, but because of the cylindrical symmetry the result for $p_\perp\gt0$ is identical). Each vertical slice in the plot shows all saddle point times for a particular final momentum point in the PMD. The corresponding SFA solutions (see the text following equation (21) for definition) are shown in dashed black lines. Clearly, the majority of orbit 1 (blue) and 2 (orange) solutions lie along the SFA solutions, indicating their close resemblance to the SFA trajectories. Between each pair of SFA 'lines' one sees additional structure composed of orbit type 3 (green) and 4 (red). These structures contain the rescattering trajectories, which will be studied in detail.

Zoom In Zoom Out Reset image size

The rescattering structure with the most impact on the transition amplitude is the framed structure in figure 5(a1), containing orbits with ionization times near the highest peak of the electric field. A zoom-in of this structure is seen in panel (b1) and a collection of corresponding trajectories displayed in panel (b2), to show the evolution in trajectory dynamics that the structure represents. The lowest orbit 4 line, labeled by box A, contains field-driven scattering trajectories; the trajectory has a turning point at the maximum of the electric field, then it returns to the ionic core and scatters after spending around one laser cycle in the field. These orbits correspond to the 'standard' orbit 4 in figure 2, and were used in previous CQSFA models (see, e.g. figure 16 in [18]). Box B marks a turning point in the orbit 4 structure. Here the field maximum occurs after the trajectories turning points and the return to the ionic core now happens within the same cycle as ionization. In this way, the orbits from boxes A and B are similar to the 'long' and 'short' orbit from the rescattered SFA [63] and HHG [25]. As we move to ionization times further from the electric field maximum, these short orbits behave less like the SFA orbits, as the Coulomb interactions become more dominant, making the trajectory turning point happen sooner compared to the field maximum, see the trajectory from box C. At box D the scattering happens before the next field maximum, while at box E we have rescattering before the laser amplitude changes sign. Interestingly, and quite unlike normal SFA orbits, these orbits are ionized at times approaching 0 electric field. For the same reason, however, their contribution to the transition amplitude is small. It is possible that these trajectories play a more prominent role for more complex targets or laser fields. In such cases, they could be an interesting probe of the system potential with which they interact strongly.

It is worth noting that other trajectory types than the ones just discussed are present in figures 5(a1) and (b1). For example, the orbit 1 and 4 structure found between box C and E on figure 5(b1) represents multipass trajectories, rescattering with the parent ion several times. Even 'chaotic' trajectories, where the electron moves in a complicated and unpredictable manner, can be found in the less connected parts of figure 5(a1), primarily for lower times ($\mathrm{Re}(t_s) \in [0, 90]$ a.u.), since these trajectories needs time interacting with the field to develop. Due to their extreme sensitivity to initial conditions, however, these 'chaotic' trajectories should have a large value of J_s , equation (22) and in turn a small contribution to the final transition amplitude.

We also note that the orbit structure shown in figure 5 can become more complicated for other potentials. One example is the argon potential used to obtain the results in figure 3(c) and (f), where a lot of additional solutions become possible, complicating the sampling as previously mentioned. This however also shows the possibility for additional types of trajectory dynamics to be involved in the description of more complicated targets.

In this section, we present and analyze PMDs made from only one or two orbit types. This is the first time that such PMDs have been produced with the CQSFA for few cycle pulses. As will be seen, the structures in the orbit-resolved PMDs can be directly related to figure 5 panel (a1) and especially (b1). This comparison is aided by the fact that the turning points in the orbit structure (e.g. box B in panel (b1)) correspond to a point on a caustic in the PMDs. This is the case since $\mathrm{Re}(t_s)$ is parameterized by the initial momentum p ₀, see equation (19). Thus, a turning point in $\mathrm{Re}(t_s)$ can be expressed as

$$\begin{align} \frac{\partial \mathrm{Re}\left(t_s\right)}{\partial \boldsymbol{p}_f} & = \frac{\partial \boldsymbol{p}_0}{\partial \boldsymbol{p}_f} \frac{\partial \mathrm{Re}\left(t_s\right)}{\partial \boldsymbol{p}_0}. \end{align} \tag{ 39 }$$

Hence, a focal point, $J_s(t) = 0$, and a caustic, will lead to a divergent turning point in $\mathrm{Re}(t_s)$, as can be observed in figure 5(b1) box B.

In figure 6, the PMDs for orbit types 3 and 4 for a 2-cycle and 4-cycle pulse are shown. As we already saw in figure 5, the somewhat arbitrary definition of the orbit types, mean that orbit types 4 and 3 changes into the other at $p_\parallel = 0$. In order to get a continuous PMD, we stitch together the contributions of the two orbit types on each half of the momentum distribution about $p_\parallel$, as indicated in the panels. Some of the features in figures 6(a) and (b) can be directly related to the distribution of $\text{Re}(t_s)$ and final parallel momentum $p_\parallel$, figure 5(a1), with the high signal in figure 6(b) originating from the orbits in the central 'branch', $\mathrm{Re}(t_s)\in[90, 140]$ a.u. as shown in figure 5(b1), and the lower signal in (a) originating from the branch below, $\mathrm{Re}(t_s) \in [50,90]$ a.u. These are the only two orbit branches expected to correspond to a strong signal, since only the photoelectrons that correspond to these solutions enter the continuum at high field values, while also having adequate time to subsequently propagate in the electric field.

Zoom In Zoom Out Reset image size

We see that the momentum value where the caustic appears on the left side of figure 6(b) ($p_\parallel \sim -1.75$ a.u.), matches very well with the orbit 4 turning point in figure 5(b1) box B, indicating that the circular interference pattern within this caustic originates from the interference of the long and short re-scattered orbit 4 (Box A and C), which is well known. The orbit 3 turning point in figure 6(b1), close to Box D, on the other hand, does not result in a visible caustic in the orbit-resolved PMD in figure 6(b) on the right-hand side (at $p_\perp = 0.25$ a.u.). This indicates that the upper part of this main structure (containing box E), associated with Coulomb-driven recolliding orbits 3 and 4 solutions, correspond to a very small transition amplitude. This situation changes at $p_\perp\sim0.6$ a.u., where a caustic appears on the orbit 3 side of figure 6(b) ($p_\parallel \sim 1.75$ a.u.). This is due to a slight change in the structure from figure 5(b1), where at $p_\perp \gt 0.6$ a.u. the orbit 3 line near box D detaches from the orbit 3 solution above and instead connects to the orbit 3 solution below that lies along the SFA solution, forming a new loop structure together with the orbit 4 solutions in Box A, B, and C. This is shown by the light gray orbits in the same panel, indicating the orbit 3 and 4 structure for $p_\perp = 1$ a.u. The formation of this ring structure increases the importance of the orbit 3 solutions around the right turning point, both leading to the caustic at $\boldsymbol{p} \sim (1.8, 0.6)$ a.u. and the appearance of interference fringes above $p_\perp \sim 0.6$ a.u. on the orbit 3 side of figure 6(b).

The difference in transition probability in panels (a) and (b) of figure 6 is clearly an effect of the two-cycle pulse used. Indeed, in PMDs for a 4-cycle pulse, as can be seen on panels (d) and (e) of the same figure, the value of transition amplitude for a given orbit type is about the same for both $p_\parallel\gt0$ and $p_\parallel\lt0$. The effect of 2 extra field cycles is largest for orbit 4, where the high-energy interference rings are now present for both positive and negative $p_\parallel$. Furthermore, one notices an additional caustic within the orbit 4 signal in panel (d) (starting at $\boldsymbol{p} \sim (1.9,0)$ a.u.) also carrying over to the orbit 3 signal. This caustic is similar to the 2-cycle caustic in panel (b), and is indeed linked to a structure similar to that in figure 5(b1), but for the last parts of the 4-cycle pulse. The additional fine interference fringes found within the boundary of the caustic (compared to the 2-cycle case) thus represents interference between trajectories originating from two different cycles of the pulse. Other than these inter-cycle interference fringes, along with a shift of the caustic to higher $p_\perp$ values, the orbit 3 signal is not much affected by the increase in field cycles.

Figures 6(c) and (f) shows the total signal from orbits 3 and 4, that is, the transition amplitude, equation (18), for orbit 3 and 4 coherently summed before taking the norm square. For the two cycle pulse, panel (c), we only have a relatively small amount of interference between the two orbits, since the difference in signal of the orbit types is large (compare figures 6(a) and (b)). The main interference effect is the emergence of fine interference fringes within the area of the left orbit 4 signal. For the 4 cycle pulse, figure 6(f), interference between orbit type 3 and 4 is more pronounced, since the difference in signal for the two orbit types is less. Especially, we see the emergence of the so-called spiral structure [44] in the central parts of the PMD. The results here show that patterns coming from interference with orbit 4 depends most strongly on the few-cycle nature of the pulse, which was to be expected from the single orbit type PMDs.

The orbits found within the structure in figure 5(b1) capture most of the orbit types 3 and 4 signal in the PMD. It is therefore interesting to investigate where in the initial momentum space these orbits originate from, and to optimize sampling of the initial conditions. In figure 7, we consider the initial momentum values for the orbits found in the central structure (figure 5(b1)) across all values of final momenta. Furthermore, the specific initial conditions for the orbits shown in figure 5(b1), with $p_\perp = 0.25$ a.u., are marked with black dots. We also include the boxes A–E from figure 5(b2), now mapped to initial momentum space. Clearly, judging by the density in the boxes in figure 7, the initial conditions in box A, B, C, and D, that are giving rise to the high signal structure in the PMDs figure 6, are located within a small area, $p_{0\parallel} \in [0;1.5]$ a.u., $p_{0\perp} \in [0;-0.2]$ a.u. The other parts of the figure, such as the orbits for $p_{0\parallel} \gt 2.5$ a.u., containing box E, were seen to have very little influence on the PMDs, even though they are densely populated in initial conditions. In addition, we have plotted the initial momentum values for the $p_\perp = 1$ a.u. orbits, as seen in figure 5(b1) in orange. Here, it is seen how the orbits giving the most signal in the PMDs (the ring structure close to box A, B, C) is even more localized in initial momentum space.

Zoom In Zoom Out Reset image size

Comparing with the total amount of sampled initial conditions, figure 1, one sees that it is only a very small amount of the sampled trajectories that matters for the resulting PMDs. In particular, there are many densely populated areas of initial conditions that are low in probability, meaning a clustering approach to sampling quickly leads to unnecessary oversampling. For more complicated potentials, e.g. for molecular targets (and to some degree even for the Ar potential used here) this could make sufficient sampling very difficult. One way to counter this would be to sample in an iterative manner: First a set of random initial conditions are chosen, their trajectories propagated and the parameters important for their contribution to the transition amplitude, primarily $\Im{S(\tilde{\boldsymbol{p}}_s, \boldsymbol{r}_s, t_s)}$ and the Jacobian of the stability matrix $J_s(t)$, see equation (18), are determined. Based on these parameters, one can estimate which regions of initial conditions should be sampled in greater detail, and prevent sampling in regions of low probability density.