Semiclassical two-step model for ionization by a strong laser pulse: further developments and applications

We review the semiclassical two-step model for strong-field ionization. The semiclassical two-step model describes quantum interference and accounts for the ionic potential beyond the semiclassical perturbation theory. We discuss formulation and implementation of this model, its further developments, as well as some of the applications. The reviewed applications of the model include strong-field holography with photoelectrons, multielectron polarization effects in ionization by an intense laser pulse, and strong-field ionization of the hydrogen molecule.


Introduction
Strong-field physics studies phenomena arising from the interaction of strong laser pulses with atoms and molecules. The most well-known examples of these highly nonlinear phenomena are above-threshold ionization (ATI), formation of the high-energy plateau in the electron energy spectrum (High-order ATI), generation of high-order harmonics (HHG) and nonsequential double ionization (NSDI), see Refs. [1][2][3][4][5] for reviews. Both experimental and theoretical approaches used to analyze these processes are constantly being improved. The vast majority of the modern theoretical methods used in strong-field physics are based on the strong-field approximation SFA [6][7][8], the direct numerical solution of the time-dependent Schrödinger equation (see Refs. [9][10][11][12] and references therein), and the semiclassical models applying classical mechanics to describe the electron motion in the continuum. The widely known examples of the semiclassical models are the two-step model [13][14][15] and the three-step model [16,17].
In the SFA ionization is described as a transition from an initial state unaffected by the laser field to a Volkov state, i.e., the wave function of an electron in an electromagnetic field. Therefore, the SFA neglects the intermediate bound states and the Coulomb interaction in the final state. The SFA provides the illustrative physical picture of many strong-field phenomena and often allows for the analytic solutions. Nevertheless, the approximations used in the SFA are strong enough and may sometimes lead to wrong results. The widely known example is the fourfold symmetry of the photoelectron angular distributions in the elliptically polarized field predicted by the SFA [18]. In contrast to this, a e-mail: nikolay.shvetsov@itp.uni-hannover.de (corresponding author) the experimental angular distributions show only the inversion symmetry: They are asymmetric in any half of the polarization plane [19]. The theoretical studies [20][21][22][23][24][25] have shown that the fourfold symmetry of the angular distributions is a direct consequence of neglecting the effect of the Coulomb potential on the electron motion in the continuum.
In most cases the direct numerical solution of the TDSE provides a good agreement with the experimental results. However, it is often difficult to understand the physical mechanism of the phenomena under study with only the numerical wave function. What is also important, the capabilities of modern computers are not unlimited. One of the most prominent examples is the strong-field ionization of molecules. The solution of the TDSE in three spatial dimensions is possible only for the simplest molecules and with selection of the most relevant degrees of freedom [26,27]. Indeed, ionization of a molecule by an intense laser pulse is much more complicated than ionization of an atom. This is because of the existence of additional degrees of freedom (nuclear motion), the associated time scales, and the complex shape of the electronic orbitals. For typical laser paremeters used in experiments nuclear motion should be treated on an equal footing with the processes induced by a strong laser field. Simultaneously, the rich nuclear structure of molecules results in orbitals of diverse symmetries.
Although the first semiclassical model (i.e., the two-step model) was formulated in 1988-1989 [13][14][15], the trajectory-based models are still widely used for description of various strong-field phenomena. This is due to a number of important advantages characteristic to the semiclassical approaches. The semiclassical models provide a great insight into strong-field processes. They allow to reveal the specific mechanism responsible for the process under investigation, as well as visualize it using classical trajectories. This point needs to be discussed in more details.
In the ATI an electron absorbs more photons than necessary for ionization. The studies of the ATI have revealed that the majority of the ionized electrons do not experience hard recollisions with their parent ions. These electrons are referred to as direct electrons. They contribute to the low energy part of the ATI energy spectrum E < 2U p , where U p = F 2 0 /4ω 2 is the ponderomotive energy. Here, in turn, F 0 and ω are the amplitude and the frequency of the laser field (atomic units are used throughout the paper). The two-step model allows to describe the spectrum of the direct electrons. In the first step of the this model an electron tunnels out of an atom. In the second step it moves along a classical trajectory in the laser field towards a detector.
There are also rescattered electrons that are driven back by the laser field to their parent ions. Upon their returns the rescattered electrons scatter from the parent ions by large angles close to 180 • . These electrons form the high-energy plateau of the ATI spectrum. The rescattering scenario provides the basis for an understanding of the HHG and NSDI. Indeed, the returning electron can recombine to the parent ion and as the result of the recombination a high-frequency photon (harmonic) radiation is emitted. Alternatively, if the energy of the scattered electron is sufficient enough, it can release the second electron from the ion, e.g., by impact ionization. The three-step model comprises the interaction of the rescattered electron with the parent ion as the third step. As the result, the three-step model provides the qualitative description of the rescatteringinduced processes.
The three-step model explained a number of features revealed in the studies of the high-order ATI, HHG, and NSDI: the cutoffs in high-order ATI spectrum [28] and HHG [16,29], the maximum angles of the angular distributions of ionized electrons [30], the characteristic recoil ion momenta in NSDI [31,32], etc. Originally the two-step and the three-step models did not account for the effect of the ionic potential on the electron motion in the continuum. The inclusion of the ionic force in the Newton's equation of motions allowed to uncover the Coulomb focusing effect [33], study the Coulomb cusp in the angular distributions of the photoelectrons [34], investigate the low-energy structures in ionization by the strong midinfrared pulses [35][36][37][38][39][40][41][42][43] (the so-called ionization surprise observed for the first time in experiment [44]), explore the nonadiabatic effects in ionization by intense laser pulses (see, e.g., Refs. [45][46][47]), etc.
The trajectory-based simulations are often (although not always) computationally less expensive than the solution of the TDSE. Furthermore, for some strongfield processes the semiclassical simulations are presently the only feasible approach. The most well-known example of such process is the NSDI of atoms by circularly [48] or elliptically polarized pulses [49][50][51], as well as the NSDI in molecules [52]. Therefore, further development of the semiclassical approaches to strong-field phenomena is an important objective.
Until recently the trajectory-based models were not able to describe quantum interference effects. However, a significant progress along these lines has been made in the last decade. The trajectory-based Coulomb SFA (TCSFA) [37,53], the quantum trajectory Monte Carlo model (QTMC) [54], the semiclassical two-step model (SCTS) [55], and the Coulomb quantum orbit strongfield approximation (CQSFA) [56][57][58][59][60] (see Ref. [61] for the foundations of the CQSFA approach) are recent trajectory-based models that are capable to reproduce interference structures in photoelectron momentum distributions of the ATI process. These models assign certain phases to classical trajectories, and the contributions of different trajectories leading to the same final momentum are added coherently.
The TCSFA is an extension of the CCSFA [62,63] that on an equal footing accounts the laser field and the Coulomb force in the Newton's equation for electron motion in the continuum. The TCSFA applies the firstorder semiclassical perturbation theory [64] to account the Coulomb potential in the phase associated with every trajectory. The same first-order semiclassical perturbation theory was used in the phase of the QTMC model. In contrast to this, the SCTS and the CQSFA approaches go beyond the perturbation theory.
The SCTS model operates with large ensembles of classical trajectories that are propagated in the continuum to find the final asymptotic momenta and bin them (and, therefore, the corresponding contributions assigned to these trajectories) in bins in momentum space. This approach is often referred to as "shooting method" (see, e.g., Ref. [37]). Instead, the CQSFA model solves the so-called inverse problem, i.e., finds all the trajectories leading to a given final momentum. This allows to avoid large ensembles of trajectories and establish a better control over cusps and caustics that are inevitable in trajectory-based simulations. The price that is to be paid is that the solution of the inverse problem is a difficult task. In addition to this, the approach with the inverse problem can often be less versatile.
In this paper we review the SCTS model, as well as two recent implementations of this model. We also discuss some of the applications of the SCTS. The SCTS model has been applied to the investigation of the intrahalf-cycle interference of photoelectrons with low energies [65], to the studies of the interference patterns arising in the strong-field photoelectron holography [66][67][68], to the analysis of the sub-cycle interference in ionization by counter-rotating two-color fields [69], to the investigation of sideband modulation by subcycle interference in ionization by circularly polarized two-color laser fields [70], etc. Here we focus on the applications of the SCTS to the strong-field photoelectron holography, study of the multielectron polarization effects, and the ionization of the H 2 molecule.
The paper is organized as follows. In Sect. 2 we review the SCTS and discuss different approaches used to implement this model numerically. In Sect. 3 we discuss the further modifications of the SCTS model: the semiclassical two-step model with quantum input and the SCTS model accounting for the preexponential factor of the semiclassical propagator. In Sect. 4 we briefly review applications of the SCTS model to the strongfield photoelectron holography. The application of the SCTS to the study of the multielectron polarization effects in the ATI are discussed in Sect. 5. In Sect. 6 we review the usage of the SCTS model to describe the strong-field ionization of the H 2 molecule. The conclusions of this colloquia paper are given in Sect. 7.
2 Semiclassical two-step model

Formulation of the semiclassical two-step model
As any semiclassical model, the electron trajectory in SCTS is calculated using classical equation of motion: where F (t) is the laser field and V (r, t) is the ionic potential. In order to find the trajectory from Eq. (1), we need to specify the initial conditions, i.e., the initial velocity of the departing electron and the starting point. In the original version of the SCTS model it is assumed that the electron starts with zero initial velocity along the laser field v 0,z = 0, but it can have a nonzero initial velocity v 0,⊥ in the perpendicular direction. We note that the application of the SFA to describe the electron motion under the potential barrier leads to a nonzero initial longitudinal velocity v 0,z = 0. The effect of the nonzero v 0,z will be discussed later. Let us first assume that the interaction of the ionized electron with the ion is modelled by the Coulomb potential. Then the starting point of the trajectory, i.e., the tunnel exit point, can be obtained using the separation of the static tunneling problem in parabolic coordinates. For the static field polarized along the z-axis we define the parabolic coordinates as ξ = r + z, η = r − z, and ϕ = arctan (y/x) and find the tunnel exit coordinate η e from the following equation: Here m is the magnetic quantum number of the initial state, I p (F ) is the Stark-shifted ionization potential, and The tunnel exit point is given by z e = −η e /2. In the general case, the ionization potential I p (F ) in Eq. (2) is given by Here I p (0) is the ionization potential in the absence of the field, and μ N,I and α N,I are the dipole moments and static polarizabilities, respectively. The index N refers to the neutral atom (molecule), and the index I stands for its ion. We note that for atom the term linear with respect to F is absent in Eq. (4). The static field F in Eqs. (2), (3), and (4) should be replaced by the instantaneous value of the laser field at the time of ionization t 0 . The instants of ionization and the initial transverse velocities are distributed in accord with the static ionization rate [71]: where κ = 2I p . Following the original formulation of the SCTS model we omit the preexponential factor in Eq. (5). For atoms it only slightly affects the shape of the electron momentum distributions. After the laser pulse terminates an electron moves in the Coulomb field only. If the electron energy at the time t = t f at which the laser pulse terminates is negative E < 0, the electron moves along the elliptical orbit, and it should be treated as captured into a Rydberg state [72,73]. The corresponding process is often referred to as frustrated tunnel ionization, see, e.g., Refs. [74][75][76][77]. It is clear that the trajectories with E < 0 should be excluded from consideration, if we are interested in ionized electrons. The latter obviously correspond to the hyperbolic trajectories (E > 0). The asymptotic momentum k of the electron is determined by its position r (t f ) and momentum p (t f ) at the time t = t f : see Refs. [73,78]. In Eq.
are the angular momentum and the Runge-Lenz vector, respectively. The magnitude of the momentum k is determined by the energy conservation: The key ingredient of the SCTS model is the expression for the phase associated with every trajectory. This phase corresponds to the phase of the matrix element of the semiclassical propagator U SC (t 2 , t 1 ) between the initial state at time t 1 and the final state at time t 2 [79][80][81] (for a text-book treatment see Refs. [82,83]). Depending on the variables used to describe the initial and final states there exist four equivalent forms of the semiclassical propagator U SC : Here r 1 (r 2 ) and p 1 (p 2 ) are the initial (final) coordinates and momenta, respectively. The phase φ 1 that corresponds to the transition from the initial state to the final state, which are both described by the position, is determined by the classical action: where H [r (t) , p (t)] is the classical Hamiltonian function that depends on the canonical coordinates r (t) and momenta p (t). The other three phases φ 2 , φ 3 , and φ 4 are related to φ 1 by the canonical transformations: Then the question arises: Which of these phases should be chosen for description of the strong-field ionization process? On the assumption that for a given ionization time the starting-point of the electron trajectory is localized in space [see Eq. (2)] and the final state is characterized by the asymptotic momentum k, the phase φ 3 is used in the SCTS model. Indeed, the strong-field ionization can be viewed as a half-scattering process of an electron that is initially localized near the atom (molecule) and detected with the final momentum k. We note that if the initial longitudinal velocity is equal to zero, the initial electron momentum p 1 is orthogonal to the initial position vector r 1 (i.e., p 1 · r 1 = 0), and therefore, the phases φ 3 and φ 4 coincide with each other. For nonzero v 0,z the term p 1 ·r 1 is to be accounted in the phase. However, in most cases this term almost does not affect the resulting electron momentum distributions.
As the result, after a partial integration, the phase corresponding to a given trajectory in the SCTS model is given by: (11) where it is assumed that the trajectory has also the initial phase exp (iI p t 0 ) that describes the time evolution of the ground state. The expression (11) can be also written as follows: To arrive at the expression (12), we use the explicit form of the Hamiltonian for an arbitrary effective potential (13) and employed Newton's equation of motion (1). This formula is applicable for any single-active-electron potential used to describe the multielectron system (atom or molecule), including pseudopotentials (see, e.g., Ref. [84] and references therein). For the specific case of the Coulomb potential, the phase (12) reads as This formula should be compared with the phase used in the QTMC model: . (15) It is seen that the QTMC phase can be obtained from Eq. (14) by neglecting the term r(t) · ∇V [r(t)] in the integrand. This term leads to the double weight of the Coulomb term in the SCTS compared to the QTMC. Therefore, the QTMC phase can be considered as an approximation to the SCTS one. The double weight of the Coulomb contribution leads to a better agreement with the TDSE results [55]. The SCTS phase (14) is divergent at t → ∞, and therefore, it is to be regularized. The regularization (see Ref. [55]) can be accomplished by decomposing the SCTS phase as and separating the time-independent part of the integrand in the term ∞ t f dt {E − Z/r (t)}. Although this time-independent part leads to the contribution lim t→∞ E (t − t f ) that diverges linearly when t → ∞, it does not produce a phase difference for electron trajectories ending up in the same bin. Indeed, the final momenta of such trajectories (and, therefore, their energies) should be considered as equal. Using the solution of the Kepler problem (see, e.g., Ref. [85]) we calculate the divergent integral analytically The parameter ξ is used to parametrize the time t and the distance r from the Coulomb center: Here, in turn, b = 1/ (2E) and g = √ 1 + 2EL 2 . The constant C in Eq. (18) is to be found using the initial conditions, i.e, r (t f ) and p (t f ). It is easy to verify that see Ref. [55]. Therefore, for trajectories arriving at the same bin, we can discriminate between the common divergent part ln 2t/ √ b 3 and the finite contributions determined by − ln (g). We note that the latter depends not only on the energy, but also on the angular momentum L, and thus is different for different trajectories interfering in a given bin. Since we obtain the following contribution to the phase accumulated in the time interval [t f , ∞] due to the Coulomb potential (see Ref. [55]): This asymptotic correction of the phase which we call post-pulse phase is missing in the QTMC model.

Implementation of the semiclassical two-step model
The expression for the phase can be conveniently treated as an additional equation in the system of the first-order ordinary differential equations for electron coordinates and velocity components following from (1). This system can be solved using the fourth-order Runge-Kutta method with adaptive step size [86]. The ability of the numerical method to change the integration step is particularly important at small distances from the Coulomb center.
It is clear that the convergence of the results must be controlled with respect to both the size of the bin in the momentum space and the number of trajectories. It is particularly convenient to control convergence by using the energy spectra. In contrast to the three-dimensional (3D) differential momentum distributions or their twodimensional (2D) cuts, the spectra are functions of only one variable. They can be easily compared to each other in, e.g., logarithmic scale.
Already the first practical application of the SCTS model has shown that a large number of trajectories is needed for convergence (see Ref. [55] for details). Typically, for the same laser parameters a thousand times more trajectories are needed for the simulations with the phase compared to a semiclassical model disregarding the interference effect. This can be expected taking into account the fine interference details of electron momentum distributions generated in ionization by strong laser pulses. For this reason, it is important to consider optimization of the codes implementing the SCTS model. The most obvious way to speed up the SCTS calculations is to use parallelization. Indeed, any trajectory-based simulation can be very easily and efficiently implemented on a computer cluster by paralleling the loop over the number of trajectories.
Another approach consists in an efficient sampling of the initial conditions, i.e., times of ionization t j 0 and initial velocities v j 0 , where index j enumerates the trajectories of an ensemble. In a standard trajectorybased approach the initial conditions are chosen either randomly or from a certain uniform grid. Neglecting interference effect the ionization probability R (k) for the final momentum k that corresponds to the bin while the similar formula for the SCTS model reads as The sums in Eqs. (22) and (23) are calculated over all n p trajectories arriving at the given bin. However, the approach sketched here is not the only possible one. Importance sampling widely used in Monte-Carlo integration (see, e.g., Ref. [86]) can be used to implement the SCTS model. We turn first to the semiclassical simulations disregarding interference. In the important sampling approach the weights (importance) of classical trajectories are accounted already at the sampling stage. More specifically, the sets of initial conditions t j 0 , v j 0 are distributed in accord with the tunneling rate w t j 0 , v j 0 and the ionization probability R (k) is given by a number of trajectories reaching the bin corresponding to the final momentum k. It is easy to see that the ionization probability in the SCTS model based on the importance sampling reads as with the initial conditions distributed in accord to the square root of the ionization probability. In many situations the important sampling technique provides faster convergence compared to the standard approach of Eqs. (22)(23). Its performance, however, depends on the laser-atom parameters and the specific part of photoelectron momentum distribution under study.

Benchmark case: ionization of the H atom
The SCTS model was compared with the QTMC approach and direct numerical solution of the TDSE for ionization of the hydrogen atom (see Ref. [55]). The 2D electron momentum distributions calculated in accord to the all three approaches are shown in Fig. 2a-c. The simulations are done for ionization by a few-cycle laser pulse linearly polarized along the z-axis and defined through the vector-potential: Here n is the number of optical cycles within the pulse, and e z is the unit vector in the polarization direction. The pulse (25) is present between t = 0 and t = t f = (2π/ω) · n. The laser field is to be calculated from Eq. (25) as F (t) = −dA (t) /dt. The factor (−1) n in (25) ensures that F z (t) for ϕ = 0 has its maximum (equal to F 0 ) for ωt = πn, i.e., at the center of the pulse. We do calculations for ϕ = 0.
It is seen that the most important features of the TDSE result are reproduced by the semiclassical models (see Fig. 1a, c, e). Indeed, the electron momentum distributions are stretched along the z-axis and show clear ATI rings as well as the pronounced interference structure in their low-energy parts. The width of the momentum distributions along the polarization direction is obviously underestimated by both semiclassical models. This is due to the initial condition v 0,z = 0 (see Ref. [55] for details).
However, a closer examination of the low-energy parts of the distributions reveals remarkable deviations Indeed, for |k| < 0.3 a.u. the photoelectron momentum distributions demonstrate pronounced fanlike interference structures, see Fig 1b, d, f. These structures are similar to the ones of Ramsauer-Townsend diffraction oscillations, see Refs. [87][88][89][90]. It is seen that the SCTS model reproduces the interference pattern of the TDSE, whereas the QTMC model predicts fewer nodal lines. This fact was attributed to the underestimate of the Coulomb potential in the expression for the phase used in the QTMC model. The comparison of the photoelectron energy spectra dR/dE shows that the QTMC and the SCTS qualitatively reproduce the ATI peaks, see Fig. 2a-c. However, both semiclassical approaches can quantitatively reproduce the amplitude of interference oscillations only for a few low-order peaks. This is related to the fact that due to the initial conditions [Eq. (5)] used in both semiclassical models too few trajectories with large initial momenta in the polarization direction are launched. This also explains why the semiclassical energy spectra fall off too rapidly with the increase of energy. In order to test this hypothesis, the initial longitudinal velocity for every ionization time is set to the value predicted by the SFA, see, e.g. Ref. [37]. This change in initial conditions leads to a better agreement between the SCTS model and the TDSE, see Fig. 3 and Ref. [55]. Therefore, the main reason of deviations of the SCTS results from the TDSE solutions is not the semiclassical treating of the electron motion in the continuum, but the fact that the SCTS model does not describe the tunneling step accurately enough.
Here we compare the semiclassical simulations with the direct numerical solution of the TDSE assuming that the latter approach is exact. However, it should be noted that the numerical solution of the TDSE can sometimes have its own limitations. This is also true for the extraction of the photoelectron momentum distributions from the time-dependent wave function, which is a non-trivial problem. Some methods used for this purpose can cause wrong results (see Ref. [91] for details).

Modifications of semiclassical two-step model
Substantial efforts have been recently made to modify the SCTS model. These modifications are aimed at providing not only a qualitative, but also a quantitative agreement with the TDSE. To achieve this goal, it is necessary to overcome the deficiencies of the SCTS model (as well as of any other semiclassical model) in description of the ionization step. The simplest way is to use the SFA formulas to distribute the initial conditions of classical trajectories. This approach dates back to the studies of Refs. [92,93]. It is used in the various semiclassical models (see, e.g., Refs. [37,[45][46][47]), as well as in the implementations of the SCTS model developed in Refs. [70,94]. We note, however, that the validity of the SFA formulas used as initial conditions for classical trajectories requires a systematic study. To the best of our knowledge, such a study has not been accomplished so far. Here we discuss two modifications of the SCTS model: The semiclassical two-step model with quantum input (SCTSQI) [95] and the SCTS model with the prefactor [94].

Semiclassical two-step model with quantum input
The SCTSQI model combines the SCTS with initial conditions obtained from the solution of the TDSE. Such a combination leads to a novel quantum-classical approach. The SCTSQI model is formulated for ionization of a one-dimensional (1D) model atom. Therefore, before reviewing the SCTSQI, we briefly discuss the solution of the 1D TDSE, as well as the application of the SCTS model in 1D case.
For the 1D model, the TDSE in the velocity gauge is given by where Ψ (x, t) is the wave function in coordinate space. The 1D soft-core Coulomb potential with a = 1.0 (see Ref. [96]) is used in Ref. [95]. The corresponding time-independent Schrödinger equation reads as: Equation (28) can be easily solved on a grid using, e.g., the well-known three-step formula for approximation of the second derivative and subsequent diagonalization. In Ref. [95] the TDSE (26) is solved using slit-operator method [97]. In the regions x b ≤ |x| ≤ x max the wave function is multiplied by a mask where x = ±x b correspond to the internal boundaries of the absorbing regions, and x max is the size of the computational box. The mask prevents unphysical reflections of the propagating wave function from the grid boundary and allows to calculate the electron momentum distributions using the mask method [98].
In the 1D case the Newton's equation for an electron moving in the laser field and the field of the potential (27) reads as The corresponding SCTS phase is given by (see Ref. [95]): We note that the ionization rate (5) in the 1D case is to be replaced by where E 0 = −0.6698 a.u. is the ground-state energy in the potential (27). Equation (30) is to be numerically integrated up to the end of the laser pulse at t = t f . The asymptotic momentum of the photoelectron can be found from x (t f ) and p x (t f ) using the energy conservation law. We note that after the end of the pulse the unbound electron cannot change its direction of motion, and, therefore, k x has the same sign as that of p x (t f ). In order to correctly apply the SCTS model in the 1D case, the post-pulse phase is to be calculated. This calculation can be performed as follows (see Ref. [95]). At first, we decompose the phase as: As in the 3D case, we separate the post-pulse phase into parts with time-dependent and time-independent integrands and disregard the linearly divergent contribution from the first part. As the result, the post-pulse phase is determined by: The divergent part of this integral can be efficiently isolated. Indeed, Eq. (34) can be equivalently rewritten as follows: Since the second divergent term in Eq. (35) depends on the electron energy E and the parameter a, it is the same for every trajectory that arrives at a given bin [k x − Δk x , k x + Δk x ]. Therefore, it does not affect the resulting interference pattern and can be omitted [95]. The post-pulse phase is determined by the first term in Eq. (35). This converging integral is easily calculated numerically. It depends on the position x (t f ) and velocity p x (t f ) at the end of the laser pulse what suggests an efficient way to calculate it by interpolation [95]. This is not a simple task to unify the direct solution of the TDSE and the trajectory-based approach in one single model. The main problem of such combination has a fundamental origin. Indeed, both the starting point and the initial velocity are needed to uniquely determine the classical trajectory. On the other hand, the Heisenberg's uncertainty principle imposes a limit to the precision with which position and momentum (as other canonically conjugated variables) can be simultaneously known. The application of quasiprobability distribution allows to extract the information from the wave function about both the coordinate and momentum.
The most widely-known examples of the quasiprobability distributions are the Wigner function and Husimi distribution [99] (see Ref. [100] for a textbook treatment). The latter can be obtained by smoothing of the Wigner function with a Gaussian weight. The Gabor transformation [101] was used in Ref. [95]. The Gabor transformation is presently widely used in studies of the ATI (see, e.g., Ref. [102]) and, especially, the HHG (see, e.g., Refs. [103][104][105]). The Gabor transform of the wave functionΨ (x, t) near the point x 0 is given by: where δ 0 is the width of the Gaussian window. The square modulus |G (x 0 , p x , t)| 2 corresponds to the momentum distribution of the particle in the vicinity of x = x 0 at time t and is just the Husimi distribution [99]. The Husimi distribution is a positive semidefinite function, which helps to interpret it as a quasiprobability distribution.
The SCTSQI model employs the solution of the TDSE in the length gauge: Two additional spatial grids containing N points are introduced the absorbing regions |x| ≥ x b : Here j = 0, ..., N and Δx = (x max − x b ) /N . In the SCTSQI the Gabor transforms of the absorbed part of the wave functionΨ ( are calculated at every time at the points x j 0,− and x j 0,+ of the grids (38). The value of the Gabor transformation at an arbitrary point belonging to D 1 or D 2 can be obtained by interpolation (see Ref. [95] for a details of the implementation of the SCTSQI model). Hence, at every time t the Gabor transform G (x, p x , t) is known on the grids in the phase-space ]. An example of the Husimi distribution obtained in the domains D 1 and D 2 at t = 3t f /2 is shown in Fig. 4. It should be stressed that the size of the computational box x max used in the SCTSQI can be much smaller than the one required to obtain accurate momentum distributions by using the mask method.  Fig. 4 The Husimi quasiprobability distribution |G (x, px, t)| 2 at t = 3t f /2 calculated for ionization of 1D model atom by a laser pulse with a duration of n = 4 cycles, wavelength of 800 nm, and peak intensity of 2.0 × 10 14 W/cm 2 . The distribution is calculated in the phase space domains D1 and D2 (see text). The points S1, S2, and S3 depicted by a magenta square, cyan triangle, and green circle, respectively show the three main maxima of the Husimi distribution. A logarithmic color scale is used This phase coincides with the phase of the semiclassical propagator describing a transition from an initial state characterized by the momentum to a final state, which is also described by the momentum value. The ionization probability R (k x ) is calculated as: where N T is the number of steps that is used in the TDSE propagation and n k is the number of trajectories arriving at the same bin centered at k x . It is important to stress that G t m 0 , x j 0 , p j x,0 is a complex function having both modulus and the phase.
The SCTSQI model was tested by comparing its predictions with the numerical solution of the TDSE and the SCTS model, see Fig. 5a, b. It is seen that the SCTSQI provides not only qualitative, but also quantitative agreement with the TDSE result. This is true for both the width of the electron momentum distributions and the positions of the interference maxima and minima. The small discrepancy in the heights of some interference peaks (see Fig. 5a) is attributed to the fact that the SCTSQI model does not account for the preexponential factor of the semiclassical matrix element. We note that as in the 3D case (see Sect. 2.3) the 1D SCTS model shows only a qualitative agreement with the fully quantum results, see Fig. 5b. Specifically, the SCTS model underestimates the width of the momentum distributions. The electron energy spectra calculated within the SCTSQI model and from the solution of the TDSE are in almost perfect agreement. Simultaneously, the spectrum calculated using the SCTS model falls off too rapidly with the increase of the energy. This is caused by the underestimation of the width of the electron momentum distributions in the SCTS model. It was shown that the phase of the Gabor transform is very important in the SCTSQI [95]. Without this phase the SCTSQI model does not provide even a qualitative agreement with the TDSE result. This could be expected, since the amplitude G (t, x, p x ) contains all the information about the quantum dynamics of the absorbed part of the wave function before it was transformed in an ensemble of trajectories. In a way the term I p t 0 in the phase (11) of the SCTS model plays the same role as the phase of G (t, x, p x ) in the SCTQI approach.
As any semiclassical approach, the SCTSQI model can visualize the physical mechanism responsible for the strong-field process under study using classical trajectories (see Ref. [95] for details). Since the initial conditions in the SCTSQI model are determined from the direct solution of the TDSE, we expect that this model will be able to provide more accurate trajectory-based pictures of strong-field phenomena compared to the standard semiclassical approaches. This advantage of the SCTSQI model should be used in studies of complicated strong-field processes. In addition to this, after some modification the SCTSQI model can be applied to studies of the rescattering-induced phenomena, especially the high-order ATI and the HHG. The ways of this modification are suggested in Ref. [95]. Finally, the extension of the SCTSQI model to the threedimensional (3D) case is straightforward and developments in this direction are on the way. Most importantly, the model solves the non-trivial problem how to choose initial conditions for classical trajectories. In the SCTSQI model these initial conditions are determined by the exact quantum dynamics.

SCTS model with preexponential factor
An efficient modification and extension of the SCTS model was proposed recently in Ref. [94]. This study for the first time investigates systematically the influence of the preexponential factor of the semiclassical matrix element (8c) (see Refs. [106,107]) that was not explicitly considered in all other versions of the SCTS. This preexponential factor for strong-field processes was calculated in Appendix B of Ref. [108]. The modulus of this prefactor that corresponds to the mapping from initial conditions to the final momentum components influences the weights of the classical trajectories. Its phase known as the Maslov phase can be identified as a case of Gouy's phase anomaly and modifies the interference structures [94]. In addition, the authors propose a novel way of solving the so-called inverse problem based on a clustering algorithm.
Since the SCTS implementation of Ref. [94] employs the SFA and the saddle-point approximation to calculate the ionization weight of the classical trajectories and their initial positions, the ionization time t 0 for each initial electron momentum k is determined by the real part of the corresponding saddle-point time t s = t 0 + it 1 . The saddle point t s satisfies the equation: The ionization probability is calculated as: Here, the summation is over all the initial momenta k leading to the final momentum k. D is the matrix element emerging when the saddle-point method is applied to calculate the SFA ionization amplitude and C Coul is the Coulomb correction of the ionization rate [64]. The phase associated with every trajectory is decomposed in Eq. (42) corresponds to the ionization step (motion under the potential barrier), and accounts for the electron motion in the continuum. We note that the phase S → coincides with the third term of the SCTS phase [see Eq. (12)]. The Jacobian J is calculated as The Maslov index ν changes at focal points, i.e, at times T when the Jacobian is zero J (T ) = 0. The change (jump) of the Maslov index when the trajectory passes through a focal point is calculated as: where the m × m matrix g is given by Here, in turn, m is the number of linearly independent directions d (i) (i = 1, ..., m), which can be found at the focal points, such that infinitesimal changes of the initial momenta in these directions k → k + d (i) do not affect k (T ) in the first order of . These changes of the initial momenta correspond to the changes of the position see Eq. (47). The Hessian Hesse r,r (H) of the Hamiltonian function is calculated with respect to the position vector r.
The inverse problem is solved in Ref. [94] by using clustering algorithms. More specifically, density-based spatial clustering of applications with noise algorithms was applied. The solution of inverse problem with clustering shows an example of the application of machine learning (see Ref. [109] for a text-book treatment) to strong-field phenomena. Other recent applications of the machine learning in strong-field physics are discussed in, e.g., Refs. [110,111].
The fact that the Jacobian is explicitly taken into account in Eq. (42) along with the solution of the inverse problem ensures the correct preexponential weight of every trajectory, namely, 1/ |J|. It should be emphasized that this weight cannot be reproduced in "shooting method", since the distribution of the trajectories over the cells in accord with their final momenta automatically creates a factor of 1/ |J| instead of the 1/ |J|. This problem was ignored in the implementation of the SCTS [55], since the implementation of Ref. [55] accounts only for the exponential factors in the trajectories weights.
The simple relation between the Jacobian in the 3D case and the corresponding Jacobian for two spatial dimensions was derived in [94] for systems (ionic potential and the laser field) with cylindrical symmetry: where k ⊥ = k 2 x + k 2 y (the field is polarized along the z-axis). This correction weight allows to obtain the results for the 3D system performing only the 2D simulations, and, by doing so, reduce the computational costs of the SCTS model significantly. We note that Eq. (50) has been already used in the SCTS simulations of Ref. [55].
The modified version of the SCTS is in excellent agreement with solution of the TDSE. This applies for both electron momentum distributions and energy spectra [94]. It is shown that the inclusion of the preexponential factors is crucial for quantitative agreement with the TDSE results. The extended version of the SCTS can be applied not only to the linearly polarized pulses, but also to non-cylindrically-symmetric laser fields, e.g., bicircular ones, see Ref. [94]. Undoubtedly the version of the SCTS developed in [94] is a valuable tool that is extremely useful in studies of strong-field ionization.

Semiclassical two-step model and the strong-field holography with photoelectrons
Development of the techniques capable to image the atomic positions that change in time in a chemical reaction will lead to a revolution in chemistry, biology, nanoscience, etc. At present there are many methods for time-resolved molecular imaging (see Ref. [112] for a review). These methods have been developed due to the prominent progress in laser technologies. This applies above all to the development of the technology for pulse compression and the emergence of free-electron lasers. Moreover, the availability of table-top intense femtosecond lasers, which led to the emergence of strong-field, ultrafast, and attosecond physics, gave a strong impulse to the development of new techniques for time-resolved molecular imaging. Among these techniques are: laserinduced Coulomb-explosion imaging [113][114][115][116], laserassisted electron diffraction [117,118], high-order har-monic orbital tomography [119,120], laser-induced electron diffraction (see, e.g., Refs. [121][122][123]), and strongfield photoelectron holography (SFPH) [124].
The SFPH method implements the widely-known idea of holography (1971 Nobel Prize in Physics awarded to Dennis Gabor, see Ref. [125]) in strong-field physics. It was for the first time shown in 2011 by Y. Huismans et al. [124] that a holographic pattern can be clearly recorded in experiment. This pattern in the electron momentum distributions is created by the signal (rescattered) and reference (direct) electrons. The SFPH can be implemented in a table-top experiment. It was shown that the holographic patterns encode a lot of spatio-temporal information about both the parent ion and the recolliding electron [124]. Last but not least, the electron dynamics can be imaged with subcycle (i.e., attosecond) time resolution. These advantages have triggered extensive studies of the SFPH, both experimental [66,[126][127][128][129] and theoretical [56][57][58][59][60]124,126,127,[130][131][132][133][134][135].
However, the first SFPH experiments [124,[126][127][128]130] investigated the ionization process and the dynamics of the electron wave packet rather than molecular structure or dynamics. This is because of the fact that for diatomic and small molecules the holographic structures are mostly determined by the long-range and the alignment-independent Coulomb potential. As the result, the short-range effect reflecting the molecular structure cannot be observed on the background of the more intense Coulomb contribution. This problem was elegantly solved in experiment of Ref. [129] by considering the difference between the normalized photoelectron holograms for aligned and antialigned molecules. This approach is based on the fact that for large scattering angles the differential cross section deviates from the Coulomb one and depends on the alignment of the molecule at the ionization instant. A similar method was also used in Ref. [66]. Various approaches were used for theoretical analysis of the SFPH: the three-step model [131][132][133]135], the SFA version that accounts for rescattering [124,130], the Coulomb-corrected strongfield approximation [124,130], the CQSFA [56,57,59,60], etc. (see Ref. [136] for recent review).

SCTS model and experimental holographic patterns
The SCTS model was applied to the simulations of the holographic interference patterns observed in the experiment [66]. In the study [66] the electron momentum distribution produced in ionization of the NO molecule were calculated for two different cases. In the first case the electron density of the highest occupied molecular orbital (HOMO) is aligned along the polarization direction, whereas in the second case this density is orthogonal to it. These distributions, as well as their normalized difference are shown in Fig. 6. To apply the SCTS model, the distributions over the initial transverse velocities are needed for both these cases. These distributions were determined using the approach based on partial Fourier transform generalized to molecules (MO-PFT) [137][138][139]. The MO-PFT approach works with the electron wave function in mixed (coordinatemomentum) representation and uses the Wentzel-Kramers-Brillouin (WKB) approximation. The MO-PFT requires the corresponding HOMO's that were obtained using the GAMESS package [140]. The semiclassical simulations are in a perfect agreement with the experimental results [66]. The simulations within the SCTS model reproduce all characteristic features of the holographic patterns. The regions of constructive and destructive interference predicted by the model of Ref. [131] that neglects the Coulomb potential are shown in Fig. 6 with white and black color, respectively. It is seen that the three-step model overestimates the spacing between the holographic fringes in the direction perpendicular to laser polarization. Therefore, the account of the Coulomb potential leads to the improved agreement between the experiment and the semiclassical simulations.

Effects of the Coulomb potential and the strong-field photoelectron holography
The three-step semiclassical model predicts different types of subcycle interferometric structures, see Ref. [131]. Various types of the holographic structures arise due to the fact that the reference and signal electrons can start from different quarter cycles of the laser field. In Ref. [67] the various types of the subcycle interference patterns revealed in [131] were calculated accounting for the Coulomb potential of the ion with the adapted version of the SCTS model. Here we sketch the main points of this adapted SCTS.
First, an ensemble of classical trajectories is launched only from the central period of a long (8 optical cycles) laser pulse. Second, the simple formula entirely neglecting the Coulomb potential, i.e., considering triangular potential barrier formed by the laser field and the ground state energy, is used for the tunnel exit point: where the sign of z e (t 0 ) is to be chosen to ensure the electron tunnels in the direction opposite to the instantaneous field F (t 0 ). This makes it possible to directly compare the resulting interference patterns with the patterns of the three-step model. Third, the weights (5) of classical trajectories were not taken into account, and the trajectories were distributed uniformly, which is justified by the fact that holographic patterns and not electron momentum distributions were calculated in Ref. [67]. Finally, a special approach instead of Eq. (23) has to be used in the semiclassical model to obtain the phase difference between the signal and reference electrons. Indeed, to calculate the phase difference we need to isolate only one kind of rescattered trajectories and only one kind of the direct ones. This is a complicated task if the Newton's equation of motion (1) is solved treating the laser field and the Coulomb force on equal footing. First of all, it is necessary to answer the question: How to distinguish between the direct and rescattered electron trajectories in the presence of the Coulomb field? Indeed, all the trajectories are, to some extent, affected by the Coulomb potential.
The following simple recipe is used in Ref. [67]. The reference trajectories were defined as those passing the ionic core at large distances and thus experiencing small-angle scattering only. More precisely, the reference electrons obey the condition v 0,⊥ k y ≥ 0. In contrast to them, the signal trajectories come close to the parent ion and undergo large-angle scattering that changes direction of the k y component compared to the initial one. Therefore, the signal trajectories can be defined as obeying the condition v 0,⊥ k y ≤ 0. However, these conditions are not sufficient to calculate the holographic structures correctly. The fact is that in the presence of the Coulomb field the mapping from the plane of initial conditions (t 0 , v 0,⊥ ) to the (k x , k y ) plane is a complicated function. For example, in the domain where the condition v 0,⊥ k y ≤ 0 defining the signal trajectories is fulfilled, this mapping is not one-to-one: Different sets of initial conditions lead to the same momentum k, see Ref. [67] for details. The separation of trajectories of different kinds can be efficiently done by using the clusterization algorithms. In Ref. [67] this trajectory separation was accomplished manually by careful inspection of the mapping (t 0 , v 0,⊥ ) → (k x , k y ). It was found that the Coulomb potential changes interference patterns significantly. Three main effects of the Coulomb field in the holographic patterns were identified in Ref. [67]. These are: shift of the interference pattern as a whole, filling of the parts of the pattern that are unfilled when the Coulomb potential is disregarded, and the characteristic kink of the interference pattern in the vicinity of k y = 0 (cf. Fig. 7a, b). This kink at zero transverse momenta was attributed to the Coulomb focusing effect [141]. However, the question remains, how sensitive are the predicted Coulomb effects to focal averaging. Therefore, further studies are required to understand which of these effects can be observed in experiment.

Semiclassical two-step model and multielectron polarization effects
The theoretical methods used in strong-field physics usually employ the single-active electron approximation (SAE). In the SAE an atom or molecule interacting with the laser pulse is replaced by a single electron. This single electron moves in the laser field and in the field of an effective potential. Therefore, the ionization is treated as a one-electron process. The SAE is a basis for understanding of many strong-field processes, including ATI and HHG [2,83]. Nevertheless, the role of the multielectron effects (ME) in strongfield and ultrafast physics has been attracting particular attention (see, e.g., Refs. [142,143] and references therein). By now many theoretical approaches aimed at the description of the ME effects have been developed. The most well-known and widely used of them are: the time-dependent density-functional theory [144] (see Refs. [145,146] for a text-book treatment), multiconfiguration time-dependent Hartree-Fock theory [147,148], time-dependent restricted-active-space [149] and timedependent complete-active-space self-consistent field theory [150], time-dependent R-matrix theory [151,152] and R-matrix theory with time-dependence [153,154], time-dependent analytical R-matrix theory [155], etc. (see Ref. [156]). There are also many semiclassical approaches capable to account for the ME effects, see, e.g. Refs. [43,78,142,[157][158][159]. The advantages of the trajectory-based models discussed in Sect. 1 are particularly valuable in studies of complex ME effects.
One of the most well-known ME effects in strong-field ionization is laser-induced polarization of the parent ion. Recently the polarization effects in the ATI have been actively studied, see, e.g., Refs. [43,78,142,[157][158][159][160]. In Refs. [161,162] and [160] the effective potential for the outer electron that accounts for the external laser field, the Coulomb interaction, and the polarization effects of the ionic core, is derived in the adiabatic approximation. It was for the first time found in Ref. [157] that the time-independent Schrodinger equation with this effective potential and accounting for the Stark-shift of the ionization potential can be approximately separated in parabolic coordinates. This separation determines a certain tunneling geometry. The emerging physical picture of the flow of the electron charge associated with the tunneling electron is referred to as tunnel ionization in parabolic coordinates with induced dipole and Stark shift (TIPIS). The semiclassical model based on the TIPIS approach and disregarding the interference effect has shown a good agreement with experimental data (see Refs. [157][158][159]) and the TDSE results [78,157].
The electron momentum distributions generated in ionization of different atoms and molecules, including Ar, Mg, CO, naphthalene, etc., are very sensitive to the ME effects accounted by the induced dipole of the ionic core [43,78,[157][158][159]. These studies consider ionization by circularly or elliptically polarized laser pulses. This is due to the fact that the effective potential derived in Refs. [161,162] and [160] is valid only at large and intermediate distances from the ionic core. In close to circularly polarized laser fields the rescattering-induced processes are suppressed (see Ref. [163]), and, therefore, the vast majority of the ionized electrons do not return to the parent ion. However, this is not true for linearly polarized field, and the applicability of the TIPIS approach in semiclassical simulations in the case of linear polarization raised questions. This problem is addressed in Ref. [156]. Furthermore, the study [156] combines the TIPIS approach with the SCTS model. The resulting two-step semiclassical model for strongfield ionization is capable to describe quantum interference and accounts for the Stark-shift, the Coulomb potential, and the polarization induced dipole potential.

Combination of the TIPIS model and the STCS
The ionic potential derived in Refs. [160][161][162] reads as: where ME effect is accounted through the induced dipole potential α I F (t) · r/r 3 . For the potential of Eq. (52) the starting point of a classical trajectory can be obtained as the tunnel exit in the TIPIS model. More specifically, the tunnel exit point is given by z e ≈ −η e /2, where η e satisfies the equation: It is seen that Eq. (53) has the additional ME term in the left-hand side compared to the equation (2). Since the ME term in the potential (52) is proportional to the laser field F(t), it is absent at t > t f . Therefore, after the laser pulse terminates, the electron moves in the Coulomb field only. This makes it possible to use Eq. (6) for calculation of the asymptotic momentum of the electron from its position and momentum at t = t f . The SCTS phase (12) with the potential V (r, t) defined by Eq. (52) reads as: In order to implement the resulting semiclassical model, the importance sampling method was used in Ref. [156]. We note that in addition to the inapplicability of the potential (52) at small distances, there exist other conditions that restrict the range of applicability of the TIPIS model (see Refs. [78,156] for details).
The study [156] focuses on the cases of Mg (I p = 0.28 a.u., α N = 71.33 a.u., α I = 35.00 a.u.) and Ca (I p = 0.22 a.u., α N = 169.0 a.u., α I = 74.11 a.u.) atoms, which have similar ionization potentials. But their static ionic polarizabilities are different by approximately two times. In order to avoid the application of the potential (52) at small distances, a special cutoff radius r C was introduced in Ref. [156], and all the trajectories entering the sphere r < r C were ignored. The remaining trajectories do not reach the vicinity of the ion. It is clear that the elimination of the whole class of the trajectories (the returning ones) depletes some parts of electron momentum distributions. However, these depleted parts usually correspond to the boundary of the direct ionization spectrum. Therefore, they do not affect the main part of the momentum distributions that provides major contribution to the ionization yield, see Ref. [156] for details.

Application of the combined semiclassical model
The 2D photoelectron momentum distributions calculated in accord with the resulting semiclassical model are shown in Fig. 8a-d. Figure 8a, c correspond to the We first discuss the narrowing of the longitudinal distributions. This effect is further illustrated in Fig. 9a, c that show the longitudinal momentum distributions obtained with and without the ME term for Mg and Ca, respectively. Since the widths of the distributions do not change due to the interference effects, the phase is disregarded in the calculations of Fig. 9a, c. The corresponding electron energy spectra are shown in Fig. 9b, d. It is seen that the spectra calculated accounting for the ME term fall off more rapidly with increase of the energy than the ones obtained neglecting the ME effects. This is a direct consequence of the narrowing of the corresponding 2D electron momentum distributions.
The mechanism underlying the narrowing effect has a kinematic origin [156]. The analysis of classical trajec- tories has shown that there is a certain class of trajectories strongly affected by the induced polarization of the ionic core. The trajectories of this class start closer to the parent ion that other trajectories and their initial transverse velocities are not too large (see Ref. [156] for details). Indeed, the force acting on the electron due to the ME polarization effect (the ME force) decays as 1/r 2 with increasing r. Therefore, this force can change the electron motion only at the initial part of the trajectory adjacent to the tunnel exit. The ME force reduces both longitudinal and transverse components of the electron final momentum, and, as the result, the trajectories belonging to this class lead to the bins with smaller k. We note that for close to circularly polarized laser pulses the ME effects result in the rotation of the 2D electron momentum distributions towards the small axis of polarization ellipse [157].
It is seen that the presence of the ME term in the equations of motion and the phase does not dramatically change the interference patterns. The interference structure is modified only in the first and the second ATI peaks and also in the vicinity of the k z axis. The analysis of the mechanism behind the polarizationinduced interference effect showed that the changes in interference patterns are mostly caused by the ME term in the equation of motion, whereas the presence of the term −3α I F · r/r 3 in the phase (54) does not play a substantial role [156].
It was found that the trajectories interfering in a given bin often have similar ME contributions to the phase and therefore, the difference of these contributions is small. This difference is the only important quantity for the interference effect. The ME contributions to the phase are similar due to the combination of the following reasons: (i) the tunneling probability is a sharp function of the laser field F (t 0 ) at time of ionization, (ii) the tunnel exit depends only on F (t 0 ) and the parameters of the atom (molecule), and (iii) only the initial part of the electron trajectory is relevant in the integral (55). Nevertheless, for atoms and molecules with large values of the ionic polarizability α I , the difference of the ME contributions to the phase is essential. As the result, the changes in the interference patters due to the ME effect can be significant. This is illustrated in Fig. 10a-d. It is seen that the number of radial nodal lines in the fanlike interference pattern at low energies is different when calculated with and without the ME term in the phase of Eq. (54). In Fig. 10c there are six nodal lines for positive k ⊥ , while in the presence of the ME term only five such lines are visible [see Fig. 10d].

Semiclassical two-step model for H 2 molecule
To the best of our knowledge, there are only a few works that apply semiclassical models accounting for the quantum interference effect to describe strong-field ionization of molecules, see Refs. [139,164,165]. The studies [139,164] extend the QTMC model to the molecular case. The SCTS model was applied to the hydrogen molecule in Ref. [165]. Two-dimensional electron momentum distributions, energy spectra, and angular distributions were compared to the ones calculated for ionization of the atomic hydrogen. The study [165] revealed substantial differences in electron momentum distributions and energy spectra as compared to the atomic case.

SCTS model for hydrogen molecule
The ionic potential experienced by a single-activeelectron in the H 2 molecule is given by Here R is the vector pointing from one nucleus to another. It is assumed that the origin of the coordinate system is located in the center of the molecule. The effective charges Z 1 and Z 2 are chosen to be equal to 0.5 a.u. [139,164]. It is obvious that the question how to distribute initial conditions of classical trajectories is more complicated for a molecule than for an atom. Presently there are two well-known approaches to this problem: Molecular quantum-trajectory Monte-Carlo model (MO-QTMC) [139,164] and MO-PFT [137][138][139]. The MO-QTMC model applies expressions of the molecular strong-field approximation (MO-SFA) [166,167]. The MO-SFA is a generalization of the SFA that was initially developed for atoms to the case of molecules. The MO-PFT model was used in Ref. [165]. The bound state orbital in the H 2 molecule is the bonding superposition of the two 1s atomic orbitals located at the centers of the atoms: The corresponding partial Fourier transform is given by (see Ref. [138]): Here θ m and ϕ m are the polar and azimuthal angles of the molecular axis, respectively, and Π atom (p x , p y , z) is the partial Fourier transform of the 1s orbital. Substituting the expression for Π atom (p x , p y , z) (see Ref. [137]) in Eq. (58) we obtain the following formula for the mixed-representation wave function of the H 2 molecule applicable just beyond the tunnel exit: As in Ref. [139], this expression (without prefactor) was used in [165] as a complex amplitude describing ionization at time t 0 with initial transverse velocity v 0,⊥ = p 0,⊥ . In the simplest case analyzed in Ref. [165] the molecule is oriented along the laser polarization direction (θ m = ϕ m = 0), and the factor in brackets in Eq. (59) is constant for a fixed internuclear distance R. This allows to use only the exponential factor of Eq. (59). Different approaches can be used to find the tunnel exit point, i.e., the starting point of the trajectory, in the molecular case. The simplest one consists in neglecting the molecular potential, i.e., considering triangular potential barrier (51). An alternative approach, the so-called field direction model (FDM) (see Ref. [78]), accounts for the molecular potential. The potential barrier in the FDM model is formed by the molecular potential and the laser field in a 1D cut along the field direction. Therefore, the tunnel exit point in the FDM model is defined by the equation: To finalize the generalization of the SCTS model to the case of the H 2 molecule, we need to obtain the phase, which is assigned to a classical trajectory. This phase is derived by substituting the potential (56) in Eq. (12): see Ref. [165]. It is seen that for r R/2 this phase corresponds to the SCTS phase for the Coulomb potential −Z/r with the effective charge Z = Z 1 +Z 2 . In contrast to this, the QTMC phase for the H 2 molecule is given by: The expression (61) can be simplified at large distances and, as the result, the SCTS phase for the H atom is reproduced. Finally, it is assumed in Ref. [165] that at the end of the laser pulse the ionized electron is far enough from both nuclei, i.e., r (t f ) R. If this condition is met, after the end of the pulse the electron moves in the Coulomb field with the effective charge Z. Therefore, its asymptotic momentum can be calculated from Eq. (6), and the post-pulse phase is determined by Eq. (21).

Application of the SCTS model to H 2 molecule
In Fig. 11 we compare the photoelectron momentum distributions calculated within the SCTS model for the hydrogen atom (Fig. 11a) and hydrogen molecule (Fig.  11b, c), see Ref. [165]. The starting point of the trajectory for H is calculated using the triangular potential barrier (51). The distribution of Fig. 11b for H 2 is also obtained for the exit point calculated from Eq (51). We note that the molecular potential is fully taken into account in the classical equations of motion (1) and in the phase (12) when calculating Fig. 11b. The electron momentum distribution of Fig. 11c corresponds to the tunnel exit obtained by using the FDM model. The electron momentum distributions shown in Fig. 11a, b are similar to each other. Therefore, it can be concluded that if the molecular potential is not accounted in calculating the starting point, the effects of the molecular structure are not visible in electron momentum distributions. This result can be expected bearing in mind that r 0 = I p /F 0 R/2 for the parameters of Fig. 11, and the distance between the ionized electron and the molecular ion increases further when the electron moves along the trajectory. As the result, the departing electron feels only the Coulomb asymptotic instead of the full molecular potential. The FDM model predicts smaller exit points as compared to the triangular barrier formula (51), see Ref. [165]. For this reason, the effects of the molecular potential (56) are visible in Fig. 11c. First, the photoelectron momentum distribution is more extended in the polarization direction. As the result, the energy spectra for the hydrogen molecule falls off slower with the increase of energy than the ones for the H atom. Simultaneously, the angular distributions in the molecular case are more aligned along the polarization direction. Second, at the same parameters of the laser pulse the holographic interference fringes are more pronounced for H 2 than for H (see Fig. 11). The comparison of the distributions calculated using the SCTS and QTMC models for ionization of the H 2 molecule is presented in Fig. 12a-d. Two different pulse envelopes were used in Fig. 12a-d. Figure 12a, b corresponds to the sine squared pulse, whereas Figure  12c, d shows the distributions obtained for the trapezoidal pulse (see Ref. [165] for details). Figure 12a, c display the momentum distributions calculated within the QTMC model, and Fig. 12b,d show the corresponding SCTS results. For the sine squared pulse these distributions have a pronounced fan-like structure in their low-energy part. For the trapezoidal envelope the fans are substituted by the characteristic blobs (see Fig. 12c, d) lying on a circle with the radius k = 0.30 a.u. Similar to the atomic case, the QTMC predicts fewer nodal lines in the interference structure at low energies than the SCTS model. This fact can be again attributed to the underestimation of the Coulomb potential in the QTMC phase [165].

Conclusions
The semiclassical models using classical mechanics to describe the electron motion after it has been released from an atom or molecule are one of the powerful methods of strong-field, ultrafast, and attosecond physics. The standard formulation of the trajectory-based models does not allow to describe the effects of quantum interference. Nevertheless, a substantial progress in simulations of the interference effects using the semiclassical models has been achieved recently. By present several trajectory-based models capable to describe the interference effects have been developed and successfully applied to the studies of the ATI. Here we discuss one of these models, namely, the SCTS.
The SCTS model allows to reproduce interference patterns of the ATI process and accounts for the ionic potential beyond the semiclassical perturbation theory. In the SCTS the phase assigned to every classical trajectory is calculated using the semiclassical expression for the matrix element of the quantum mechanical propagator [79][80][81]. As the result, the SCTS model yields a good agreement with the direct numerical solution of the TDSE, better than, e.g., the QTMC model applying the first-order semiclassical perturbation theory to account for the Coulomb potential in the phase. Here we review further developments and applications of the SCTS. At first, we review the formulation of the SCTS and its numerical implementation. The application of the model was illustrated in the case of the H atom. We next turn to the further developments of the SCTS: the SCTSQI model [95] and the SCTS model with the prefactor [94]. In the SCTSQI model the initial conditions for classical trajectories are determined from the exact quantum dynamics of the wavepacket. For ionization of the 1D atom the SCTSQI model yields not only qualitative, but also quantitative agreement with the numerical solution of the TDSE. Further work is needed to accomplish the generalization of the SCTSQI on the 3D case. The developments in this direction have already begun. The quantitative agreement with the TDSE was also achieved by the extension of the SCTS model that accounts for the pref-  [165]). The molecule is oriented along the laser polarization direction (z-axis). A logarithmic color scale in arbitrary units is used actor of the semiclassical matrix element. Furthermore, the 3D implementation of the SCTS [94] has a number of other important modifications.
We discuss the application of the SCTS approach to the SFPH. The semiclassical simulations within the SCTS model are in perfect agreement with the results of the recent experiment [66]. The model is able to reproduce all characteristic features of the observed holographic patterns. The SCTS model also allows to investigate the effect of the Coulomb potential on the holographic structures. Three main Coulomb effects in the interference patterns were predicted [67]. However, it should be investigated how sensitive are these Coulomb effects to focal averaging. This further work will allow to understand, which of the predicted effects can be observed.
We also present a quick review of the application of the SCTS to study of the multielectron polarization effects. We discuss the modification of the SCTS model accounting for the multielectron polarizationinduced dipole potential. The semiclassical simulations predict narrowing of the electron momentum distributions along the polarization direction. This narrowing arises due to the focusing of the ionized electrons by the induced dipole potential. Furthermore, the polarization of the ionic core can also modify the interference patterns in electron momentum distributions.
Finally, we briefly reviewed the extension of the SCTS model to ionization of the hydrogen molecule. The SCTS model for the H 2 can be generalized to an arbitrary laser polarization and orientation of the molecule, as well as to heteronuclear and polyatomic molecules. We believe that these generalizations being combined with the extended versions of the SCTS will result to an emergence of powerfool tools for studies of the strong-field processes.
Funding Open Access funding enabled and organized by Projekt DEAL.

Data Availibility Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All data included in this manuscript can be reproduced by performing simulations in accord with the approaches described here, or are available upon request by contacting the corresponding author.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.