It is all about phases: ultrafast holographic photoelectron imaging
C. Figueira de Morisson Faria A. S. Maxwell

	
	Ultrafast holographic photoelectron imaging]It is all about phases: ultrafast holographic photoelectron imaging	
	
	
	Department of Physics & Astronomy, University College London
 Gower Street, London WC1E 6BT, United Kingdom
	c.faria@ucl.ac.uk, andrew.maxwell@ucl.ac.uk
	
	
	
	
		Photoelectron holography constitutes a powerful tool for the ultrafast imaging of matter, as it combines high electron currents with subfemtosecond resolution, and gives  information about transition amplitudes and phase shifts. Similarly to light holography,  it uses the phase difference between the probe and the reference waves associated with qualitatively different ionization events for the reconstruction of the target and for ascertaining any changes that may occur. These are major advantages over other attosecond imaging techniques, which require elaborate interferometric schemes in order to extract phase differences. For that reason, ultrafast photoelectron holography has experienced a huge growth in activity, which has led to a vast, but fragmented landscape. The present review is an organizational effort towards unifying this landscape. This includes a historic account in which a connection with laser-induced electron diffraction (LIED) is established, a summary of the main holographic structures encountered and their underlying physical mechanisms, a broad discussion of the theoretical methods employed, and of the key challenges and future possibilities. We delve deeper in our own work, and place a strong emphasis on quantum interference, and on the residual Coulomb potential.
	
	
[
    C. Figueira de Morisson Faria and 
	A. S. Maxwell

    December 30, 2023
======================================================




§ OVERVIEW

Attosecond (10^-18s) science deals with some of the shortest time scales in nature, and has emerged from the study of the interaction of matter with strong laser fields, typically of intensities around I=10^13W/cm^2 or higher. The key idea behind it is to probe and ultimately steer  electron dynamics in real time (for reviews see, e.g., <cit.>). Among other applications, electrons create or destroy chemical bonds, carry energy in biomolecules, and information in human-made devices, so that the above-mentioned control has the potential to revolutionize many areas of knowledge. For that reason, strong-field phenomena such as high-order harmonic generation (HHG) or above-threshold ionization (ATI) have established themselves as powerful imaging tools.


This is made possible because their underlying physical mechanisms, laser-induced recollision or recombination <cit.>, take place within a fraction of a field cycle.  
In ATI, an electron is released close to the peak of the field and, if it reaches the detector without further interaction, it is called a direct ATI electron. If it is driven back by the field and rescatters with the parent ion, it is known as a rescattered electron  <cit.>. For HHG, instead of a laser-induced collision, recombination with a bound state of its parent ion will occur, and the kinetic energy acquired in the continuum will be released as high-frequency radiation <cit.>. The time for which the electron will return will be close to a field crossing. Hence, the driving field will dictate the time window for which ionization and recollision occur. This has led to a wide range of applications such as attosecond pulses <cit.>, high-order harmonic spectroscopy <cit.> and photoelectron holography <cit.>.

The information about the target obtained upon recollision is imprinted in the high-harmonic or photoelectron spectra. Particularly in the case of 
photoelectrons, light-induced electron diffraction (LIED) has called a great deal of attention since it has been first proposed <cit.>, due to the high electron recollision currents of around 10^11A/cm^2 incident on the target <cit.>. This is several orders of magnitude higher than what is available in conventional time-resolved electron microscopy <cit.>.
The physical picture outlined above also makes it possible to associate specific transition amplitudes with electron trajectories. 
Since there may be many pathways for an electron to reach the detector with the same final energy, the corresponding transition amplitudes will interfere. Due to allowing an intuitive physical interpretation, this trajectory-based picture has been incorporated in a myriad of theoretical approaches that allow for quantum interference. Examples are the Strong-Field Approximation (SFA), the Eikonal Volkov Approximation (EVA), Analytical R-Matrix theory (ARM), the Quantum Trajectory Monte Carlo method (QTMC), the Semiclassical Two-Step model (SCTS), the Coulomb Volkov Approximation (CVA), Quantitative Rescattering Theory (QRS), the Coulomb-Corrected Strong-Field Approximation (CCSFA), and the Coulomb Quantum-orbit Strong-Field Approximation (CQSFA). For a review of many of these methods see, e.g., <cit.> and Sec. <ref>.

Throughout the years, ultrafast imaging has moved from purely structural questions such as the reconstruction of molecular orbitals <cit.>, to dynamic effects such as electron migration in larger molecules <cit.>. In all these studies, quantum interference plays a huge role. A vital and highly non-trivial issue has been how to retrieve relative phases from the high-order harmonic or photoelectron signal. This is exemplified by the fact that it took six years between the seminal paper reconstructing the highest molecular orbital (HOMO) of N_2 using HHG, and the experimental retrieval of the molecular phase shifts using the Reconstruction of Attosecond Burst By Interferences of Two-photons Transitions (RABBIT) technique <cit.>. 
In fact, in <cit.> it was necessary to postulate such phase shifts in order to achieve a successful bound-state reconstruction. 

The wish to record the magnitude and the phase of photoelectron scattering amplitudes has led to the inception of ultrafast photoelectron holography. In conventional holography <cit.>, interference patterns between a reference and a probe wave are recorded and the phase differences are used to construct the hologram.  Similarly, in ATI one may employ the quantum interference between different types of electron trajectories to obtain this information. This has been first proposed in <cit.>.
Originally, the definition “holographic structures” referred to those 
generated by 
the interference of direct trajectories, for which there is no interaction with the core subsequent to ionization, i.e., “the reference”, with those recolliding elastically with the core, i.e., “the probe” <cit.>. As the field progressed this idea has been relaxed and generalized by some research groups, which refer to holographic patterns also as resulting from  the interference between trajectories undergoing different types of rescattering. 

Within this generalized framework, several holographic patterns have been predicted theoretically, some of which are easily identified in experiments, and some of which are fairly obscure.
Well known examples 
are the spider-like structure caused by the interference of two types of forward-scattered trajectories <cit.>, and fork-  <cit.> or fishbone-like <cit.> structures resulting from backscattered trajectories interfering with forward-scattered orbits. The fishbone holographic pattern is particularly sensitive to the target structure and has been employed in <cit.> to probe diatomic molecules. 
A largely unexplored and to a certain extent contentious issue is how the residual long-range potential influences the contributing trajectories and the holographic patterns. While earlier models incorporate the binding potential only upon rescattering and employ Coulomb-free orbits during the propagation in the continuum <cit.>, in the past few years there have been detailed studies of Coulomb-distorted orbits <cit.>, whose topology is significantly altered by the residual potential <cit.>.  A recent example  is the fan-shaped structure that forms near the ionization threshold. It is a known fact that this pattern is caused by the joint influence of the external field and the long range of Coulomb potential, but the precise mechanism behind it has raised considerable debate. For instance, in <cit.> it has been attributed to resonances with specific bound states of the Coulomb potential with a well-defined orbital angular momentum. In <cit.>, it has been related to the interference of several types of trajectories. In <cit.> we have shown that this structure is a type of holographic pattern resulting from the interference of direct and forward-deflected trajectories.  The residual Coulomb potential causes an angle-dependent, radial distortion in a pattern known as the “temporal double slit”. If Coulomb distortions are not incorporated, the temporal double slit corresponds to the interference between two reference, non scattering events and it is thus not classified as holographic <cit.>. 

There are currently three main research trends related to photoelectron holography: more complex targets, longer wavelengths, in the mid- or far-IR regime, and tailored fields. These trends are not mutually exclusive and may be combined in subfemtosecond imaging and control of electron dynamics.  Longer wavelengths λ lead to larger ponderomotive energies U_p ∼ Iλ^2, and thus provide the electron with larger kinetic energies upon return <cit.>. This implies that the holographic patterns will be clearer and closer to those predicted by simplified models, as Coulomb distortions become less prominent <cit.>. Furthermore, large ponderomotive energies may be achieved with lower driving-field intensities, which means that ionization of the target may be substantially reduced. This makes the mid-IR regime more favourable for probing complex systems, as for larger molecules the ionization potentials are typically lower. Longer wavelengths have been successfully used in LIED of organic molecules, in order to retrieve bond lengths and track bond dynamics in real time <cit.>.  They have also been used in key publications on photoelectron holography, such as those in which spider-like patterns were first identified <cit.>. Targets such as molecules or multielectron atoms allow for assessing/probing dynamic changes such as those involving the coupling of electronic and nuclear degrees of freedom, core polarization, resonances, electron-electron correlation or charge migration. 
Finally, external driving fields can be tailored with the purpose of manipulating specific types of trajectories, and thus expose holographic patterns which are obfuscated by more prominent structures. Examples are orthogonally polarized <cit.> or bicircular <cit.> fields.

However, this hugely popular research field has evolved in a highly fragmented way. This lack of consensus includes the main holographic patterns encountered, the theoretical approaches developed and the current trends. The aim of the present review is to compare, unify and organize the vast research landscape around photoelectron holography, from its early days to the current trends, with special emphasis on our own work on the subject. This includes a summary of the main interference patterns identified in experiments (Sec. <ref>), and the theoretical models used to model photoelectron diffraction and holographic structures (Sec. <ref>). In  Sec. <ref>, we delve deeper into the Coulomb Quantum-orbit Strong-Field Approximation, which is the method derived and employed in our publications. Subsequently, in Sec. <ref>, we discuss the key results obtained with the CQSFA. This includes the types of orbits encountered (Sec. <ref>), single-orbit photoelectron momentum distributions (Sec. <ref>) and several types of holographic structures (Sec.  <ref>). Finally, in Sec. <ref> we conclude the review with a summary of the existing trends. We use atomic units throughout.



§ HOLOGRAPHIC PATTERNS AND TYPES OF RESCATTERING


In order to understand photoelectron holography, one must first discuss laser-induced photoelectron diffraction (LIED).  In LIED, photoelectrons emitted in above-threshold ionization (ATI) are employed to image specific targets. As an imaging tool, LIED exhibits three key advantages: (i) It allows for high electron currents, which are orders of magnitude higher than those in standard photoelectron microscopy; (ii) The electron emission is highly coherent, as it is controlled by the external laser field; (iii) It allows for resolving dynamic changes in the target with subfemtosecond resolution.
Since its early days in the late 1990s <cit.>, LIED has led to a wide range of applications such as the tracking of the coupling between electronic and nuclear degrees of freedom, giant resonances, multi-orbital effects <cit.> and bond dynamics in real time <cit.>. 
One should note, however, that the concept of LIED is broader than that of photoelectron holography, as it does not necessarily involve rescattered trajectories. In fact, even direct ATI trajectories, in which the electron reaches the detector  without rescattering, carry information about the electron's parent ion. Rescattered trajectories are however expected to be far more sensitive with regard to the target, so that they have been widely studied in the context of LIED and photoelectron holography.  In this section, we will commence by discussing both and will subsequently focus on specific holographic patterns that have been predicted or measured. Examples of laser-induced electron diffraction patterns that are not holographic structures are ATI rings, the spatial and the temporal double slit.


 §.§ ATI rings


Ionization events associated to ATI happen at specific times, which repeat with the periodicity of the laser field. Thus, the transition amplitudes related to events separated by an integer number of cycles will interfere. The periodicity of the field will lead to discrete peaks of energy

    𝐩_f^2/2=nω-U_p-I_p,

where 𝐩_f is the electron momentum at the detector, ω the field frequency, U_p is the ponderomotive energy, n is an integer and I_p the ionization potential. The remaining frequencies will be averaged out.  In photoelectron momentum distributions, this inter-cycle interference will manifest itself as concentric rings centered around the origin of the plane spanned by the electron momentum components parallel and perpendicular to the driving-field polarization. In the limit of an infinitely long pulse N_c →∞, intercycle interference is described by a Dirac-delta comb <cit.>, i.e., the peaks given by (<ref>) are infinitely sharp. In <cit.> we show that, for a  monochromatic field with a finite number N_c of cycles, the probability modulation associated with inter-cycle interference reads as 

    Ω_N_c(𝐩_f)=cos[2π i N_c/ω(I_p+U_p+1/2𝐩_f^2)]-1/cos[2π i /ω(I_p+U_p+1/2𝐩_f^2)]-1

and that the ATI rings remain invariant if the Coulomb potential is incorporated.



 §.§ Temporal and spatial double slits


Temporal and spatial double slits are interference patterns that were first defined for direct ATI, in a theoretical framework for which the residual Coulomb potential is neglected in the electron propagation. A well-known approach that allows for tunneling and quantum interference and falls within this category is the strong-field approximation (SFA). The SFA neglects the electric field when the electron is bound and the binding potential when the electron is in the continuum. Formally, the SFA may be viewed as a Born series with a field-dressed  basis, for which there is a clear-cut definition of whether rescattering events are present or absent. Direct ATI electrons correspond to the zeroth order term of such a series, and do not undergo any act of rescattering.
The energy of the direct ATI electrons extends up to 2U_p. This corresponds to the maximal kinetic energy that a classical electron may acquire in the absence of the Coulomb potential. Rescattered ATI electrons correspond to the events determined by the first-order term of the abovementioned series, for which a single act of rescattering is present. 


A spatial double-slit type of interference is present for direct ATI in aligned diatomic molecules. There will be electron emission at different centers in the molecule, so that the associated transition amplitudes will carry different phases. This will leave quantum-interference  imprints in the ATI spectra and photoelectron angular distributions, which will provide structural information about the molecule, such as 
the internuclear distance <cit.>. This idea has been initially proposed in <cit.>, and subsequently explored by several research groups both theoretically <cit.> and experimentally <cit.>. Since the mid 2000s, it has also been generalized to account for more complex types of molecules, with many scattering centers, multi-electron dynamics, coupling of different degrees of freedom, and different orbital shapes (for a review see, e.g., <cit.>). Spatial double slits are also present for high-order harmonic generation <cit.>.

The temporal double slit involves the interference of the two different types of trajectories that occur in direct ATI, i.e.,  the long and the short orbits. An electron along the short orbit will be released in the continuum towards the detector, while the electronic wave packet following the long orbit will start on the opposite side of the ion, return and eventually reach the detector <cit.>. The electric field at both instants will have opposite signs and the same amplitude.   Diffraction patterns coming from ATI direct trajectories also carry information about the geometry of the bound state from which the electron has been released, and can be manipulated using tailored fields.   This is exemplified in our previous work for both the temporal and spatial double slits in elliptically polarized fields <cit.>. It is important to stress that, for the long direct ATI orbit, return is allowed but not rescattering. This is a particularity of the methods used in these studies, which neglect the residual binding potential in the continuum. If the residual Coulomb potential is incorporated, the distinction between direct and rescattered electrons become blurred.   More details will be provided in Sec. <ref>.

The role of rescattering in laser-induced electron diffraction in diatomics has been first proposed in <cit.>. Therein, the emphasis was on how recollision at spatially separated centers would affect the above-threshold ionization photoelectron angular distributions. Nonetheless, key ideas of how to construct simplified classical models in order to compute phase shifts are already to be found. Similar ideas have been used in <cit.> in order to construct holographic patterns.



 §.§ Rescattered trajectories and photoelectron holography


Photoelectron holography uses the phase difference between a recolliding (“the probe”) and a direct (“the reference”) ATI electronic wave packet to map a specific target, as briefly mentioned in Ref. <cit.>. These wavepackets may be associated with different types of trajectories, employing the recollision physical picture. This seminal paper touches upon the possibility of photoelectron holography as an imaging tool and hints at it being able to resolve amplitudes and phases. Its main focus is however on images obtained by electron-self diffraction of atoms and diatomic molecules, and on how to maximize their resolution. A detailed account of the possible distortions that may be caused by the external driving field is provided, together with prescriptions for avoiding them. 
Nonetheless, there is an extensive discussion of the types of interference that may occur, and on rescattered trajectories, which are classified as forward-scattered and backscattered according to their scattering angle.  Backscattered trajectories may lead to very high photoelectron energies and have been widely studied since the 1990s, in the context of the ATI plateau[The plateau is a well-known structure with ATI peaks of comparable intensities, which may extend up to the energy of 10U_p. For a review see <cit.>, and for seminal plateau measurements see <cit.>.]. Forward-scattered trajectories, in contrast, typically yield a much lower photoelectron energy range, which is comparable to that of the direct ATI trajectories. They have only called the attention of the strong-field community much more recently. This sudden shift of focus may be attributed to two main reasons:

	
  *  The existence of a myriad of low-energy structures, which have been identified in experiments and theoretical studies. Within the past decade, there have been several measurements of  near-threshold enhancements in ATI spectra for long-wavelength driving fields. An example of such enhancements is the so-called Low-Energy Structure (LES) reported in <cit.> and attributed to the interplay between the binding potential and the external field.  Other examples are the very low-energy structure (VLES)<cit.> and zero energy structure (ZES)[Sometimes called the near-zero energy structure (NZES)<cit.>] <cit.>. Although classical in nature, these features require a strong interaction with the core, which may be viewed as scattering.  Indeed, in <cit.> it has been shown that the energy bunching of neighboring low-energy trajectories that turn around the core lead to a series of peaks whose existence neither depends on the range nor on the shape of the potential <cit.>. The binding potential however serves to focus these trajectories and determines the strength and absolute position of these peaks <cit.>. Support for this interpretation has been provided in subsequent work, in which a detailed investigation of the LES has been performed using classical trajectories, analytical models and a systematic mapping of the initial conditions on the plane of the final momentum <cit.>.  Therein, it has been shown that the binding potential considerably influences the transverse and longitudinal momentum components, and the role of Coulomb focusing and Coulomb defocusing is investigated. The residual Coulomb potential leads to looplike orbits that give rise to a caustic, and  a bunching similar to that described in <cit.> occurs. The changes in the electron transverse momenta caused by the potential also play an important role. The existence of caustics, orbits whose transverse momenta are modified by the potential, and their relation to the LES have also been discussed                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 in <cit.>. Within the context of the analytical R-Matrix theory (ARM), the LES and the near-zero energy structure are associated with steep changes in the imaginary part of the action, which may be related to the presence of cusps <cit.>.
	Multiple forward-scattering events in the LES have been addressed in <cit.>. For seminal work on Coulomb focusing see, e.g., <cit.>.
	
	
  * Recent SFA results for Coulomb-type potentials show that the contributions from the rescattered trajectories are not obfuscated by those of the direct trajectories.  It was widely believed that, for photoelectrons whose energies extend up to the direct-ATI cutoff energy 2U_p, the direct electron contributions would prevail and thus obfuscate those from scattered trajectories. A largely overlooked issue is however the  large scattering cross section of the Coulomb potential, which enhances the contributions from rescattered electrons by several orders of magnitude. Extensive studies of such types of rescattering have been done in the context of photoelectron holography. Features such as a three-pronged fork-type structure and rhombi that occur for scattering angles almost orthogonal to the polarization axis have been associated with such trajectories <cit.>.


The interference between direct 
and rescattered 
trajectories and its relation to photoelectron holography was first explored in detail in <cit.>. Therein, four types of holographic structures, provided for clarity in Fig. <ref>, have been obtained assuming several types of intra-cycle interference. The phase shifts between the probe and the reference depend on the instant of ionization and recollision.  Explicitly, holographic patterns are attributed to the interference between (i) long direct orbits and forward scattered trajectories; (ii) short direct orbits and forward scattered trajectories; (iii) long direct orbits and backscattered trajectories; (iv)  short direct orbits and backscattered trajectories. Therein, it is argued that only pattern (i) would be easily identified in experiments, as it is located along the polarization axis. 
One should note, however, that the patterns in <cit.>  markedly differ from those in experiments and in ab-initio computations. An exception is perhaps pattern (i), which resembles the spider-like structure identified in <cit.>, and is attributed to the same type of interference (see discussion below). Nonetheless, the model in <cit.> does not reproduce the converging fringes near the origin. 

This follows from the fact that, in <cit.>,  simplifying assumptions have been performed. First, the influence of the residual Coulomb potential has been neglected, which means that the real topology of the orbits and the Coulomb phases are not taken into consideration. Second, the distance from the origin at which the electron is born and the ionization rates are  kept fixed. This means that contributions from ionization times close to a field crossing are overestimated. Recent work <cit.> has shown that, if the residual Coulomb potential is taken into consideration, one must considerably change the time intervals used in order to reproduce the features reported in <cit.>.



 §.§ Most common holographic structures
 



  §.§.§ The spider.


The so-called “spider” is possibly the best known  ultrafast holographic structure. It consists of broad fringes nearly parallel to the polarization  axis extending up to very high photoelectron energies. It has been first measured in ATI from a metastable state of Xenon in fields of intensities of the order of 10^11W/cm^2 and wavelengths around 7μ m <cit.> (for even lower frequencies see <cit.>). Shortly thereafter, this structure has been identified experimentally for ATI in rare gases starting from the ground state, with near- and mid-IR driving fields <cit.> of much higher intensities (∼ 10^13W/cm^2) and has been reported by several groups in atoms <cit.> and molecules <cit.>. Traces of the spider can also be seen in <cit.>.


In <cit.>, the phase differences observed have also been associated with the interference of two types of forward-scattered trajectories starting from the same half cycle of the driving field. A simple model of this structure that considers the interference of a plane wave with a Coulomb scattering wave has been proposed in <cit.>. It reproduced not only the spider, but also inner spider-like structures that occur near the ionization threshold and are associated with multiply recolliding orbits. Some features observed in this structure, such as the fringes parallel to the driving-field polarization, as well as the contributing orbits, were predicted qualitatively and compared with the experiments <cit.>. This early model was however too simple to account for the diverging patterns close to the threshold and the finer details of the spider-like pattern. Trajectory based models that allow for the residual Coulomb potential and carry quantum phases, such as the Coulomb Corrected SFA (CCSFA) and its variations <cit.>, the Coulomb Quantum-orbit Strong-Field Approximation (CQSFA) <cit.> and the Quantum Trajectory Monte Carlo method (QTMC) <cit.> reproduce the spider in striking detail, in agreement with experiments and ab-initio methods.  An alternative explanation for the spider in terms of glory rescattering has been proposed in <cit.>.
The type of interference that leads to the spider has also been associated with side lobes encountered in angle-resolved ATI photoelectron distributions <cit.>. Studies of how the spider behaves with regard to the field parameters have been reported in <cit.>.

It has been initially argued that the spider-like structure is not very sensitive to the target due to being observed for many different species, and stemming from the interference of two distinct types of forward scattered trajectories, whose interaction with the core is brief.  Further studies however reveal that the spider-like fringes exhibit an offset for aligned molecules <cit.>, and carry information about nodal planes  and rotational degrees of freedom <cit.>. An ab-initio computation of photoelectron momentum distributions for diatomic molecules has in fact found that the spider and other structures are very sensitive with regard to the molecular orientation, nodal planes and the coupling of different continua <cit.>.  For traces of the spider in the context of a multielectron computation for CO_2 see <cit.>. Recently, it has been proposed that phase differences in the spider-like pattern may be used to resolve   electron motion in molecules by preparing the system in a non-stationary coherent superposition of states <cit.>.



  §.§.§ The fishbone structure.

A less prominent holographic pattern is a fishbone structure that forms for moderate to high photoelectron energies, with fringes nearly perpendicular to the polarization axis. This structure stems from the interference of direct and backscattered trajectories. Since backscattered trajectories spend in principle a longer time near the core region, it is expected that they will be more sensitive to the target structure than  the forward-scattered trajectories that lead to the spider <cit.>.  Typically, however, the fishbone structure is obfuscated by the spider-like fringes, so that the latter must artificially be removed. A simple SFA computation proposes several possible types of interfering trajectories for this structure <cit.>. Therein, however, several disagreements with ab-initio methods are pointed out and related to the absence of the residual Coulomb potential in the model. Recently, we have shown that the presence of the Coulomb potential will lead to spiral-like structures that are also usually obfuscated by the spider, but may be identified for large scattering angles <cit.>.  Within the SFA framework, it has also been proposed that structures involving backscattered trajectories could be retrieved in a realistic setting by using orthogonally polarized fields <cit.>. For an early experimental example of the fishbone structure see Fig. <ref>.



  §.§.§ The fan.




The “fan” is a widely known pattern consisting of radial interference fringes, which forms near the ionization threshold <cit.>. This pattern does not occur for short-range potentials, so that it requires the interplay between the laser field and the Coulomb potential to be able to form. In fact, it is absent both in SFA <cit.> computations and in the photodetachment of negative ions <cit.>, which show fringes nearly perpendicular to the driving-field polarization in this region. In the literature, the fan-shaped structure has been interpreted as a resonant process involving intermediate bound states with a specific angular momentum, both in the seminal experimental work <cit.> and in subsequent theoretical papers <cit.>. Specifically, in <cit.> the fan's degradation for few-cycle pulses suggested a resonance-like character, and in 
<cit.> the fan has been related to Ramsauer-Townsend fringes, by comparing it to ab-initio methods and with CMTC computations. Therein, laser-dressed hyperbolic orbits with neighboring angular momenta have been identified as those responsible for the fan. One should note, however, that classical-trajectory methods do not allow for quantum interference. This means that only indirect statements about how the fan forms can be inferred from this method. 

In <cit.>, the number of fringes in the fan have been related to a specific angular momentum, which gives the minimal number of photons necessary for the electron to reach the continuum, and empirical rules for predicting the number of fringes have been provided. These rules, however, work well only in the multiphoton regime[Roughly speaking, the multiphotoen regine is characterized by a Keldysh parameter γ=√(I_p/(2U_p))>1, where I_p is the target's ionization potential. In this regime, tunneling is not the prevalent ionization mechanism, but, rather, the electron is freed by multiphoton absorption. In the multiphoton regime, resonances are extremely important. For a a more rigorous definition see <cit.>.]. This, together with the fact that he fan-shaped structure is also present for a wide range of species, much lower frequencies and higher intensities <cit.>, indicates that resonances cannot be the sole mechanism behind it. 

It is worth mentioning that the fan-shaped structure has also been identified using orbit-based methods that allow for quantum interference, such as the Coulomb-Volkov Approximation (CVA)<cit.>, the Coulomb-corrected strong-field approximation (CCSFA) <cit.>, the Quantum Trajectory Monte Carlo (QTMC) method  <cit.> the semiclassical two-step model (SCTS) for strong-field ionization <cit.>, and the Coulomb Quantum-orbit Strong-Field Approximation (CQSFA) <cit.>. However, if such methods employ Coulomb-distorted trajectories and solve the direct problem (see discussion in Sec. <ref>), they require a huge number of contributed orbits to obtain converged photoelectron momentum distributions. This makes it difficult to ascertain how the fan forms. The CQSFA is an exception as it only needs a few contributed trajectories for each value of the final momentum. This means that they can be switched on and off at will. 
For that reason, in our previous publications <cit.> we have focused on the question of what types of orbits are responsible for specific patterns, including the fan.  We have found that the fan  may be viewed as a holographic pattern resulting from the interference of direct with forward deflected trajectories. The latter undergo an angular-dependent distortion due to the presence of the Coulomb potential, which is maximal for momenta parallel to the laser-field polarization, and cancels out for momenta perpendicular to it. These distortions take place in a fraction of the laser field, which are much shorter than typical timescales for which resonances occur. For short-range potentials or the SFA, these distortions are absent. In this case, one may speak of return, but not of deflection, and the temporal double slit is recovered. The CVA obtains the fan by considering Coulomb distortions in the final continuum state, but keeps the same orbits as the SFA. Thus, it does not assess how the Coulomb potential modifies the orbits <cit.>. 
One should also note that the right number of fringes has only been obtained with the CQSFA and the SCTC methods. This is made possible due to an extra phase, which was overlooked in the remaining orbit-based methods mentioned in this section.  More details about this phase will be provided in Sec. <ref> and can be found in <cit.>. The CQSFA, in conjunction with analytic approximations, is also an excellent means to disentangle sub-barrier and continuum distortions and the associated phase differences <cit.>.
Finally, near-threshold fan-shaped fringes have also been observed in molecules <cit.>. 
Interesting features encountered therein are a marked enhancement in the fan-like fringes, which has been related to electronic-nuclear coupling in H_2 <cit.>, and a suppression associated to the population of Rydberg states in dimers <cit.>.



  §.§.§ Fork and Rhombi-like structures.

In <cit.>, a three-pronged fork-like structure was observed experimentally for photoelectron momentum distributions using xenon in long, mid-IR driving pulses (see Fig. <ref> for a depiction of this structure). This structure is nearly orthogonal to the laser polarization axis and is very sensitive to the pulse length and frequency. In fact, it is markedly enhanced with increasing wavelength and disappears for few-cycle pulses. Physically, it has been explained in terms of low-energy forward scattered trajectories. An SFA computation has shown that this structure is universal and does not necessarily require Coulomb-type potentials to form. However, the divergent Coulomb scattering cross section makes it visible by allowing it to rise above the contributions of the direct electrons.  

A key difference between the fork and the previous holographic patterns is that its existence is of a classical nature, similar to that leading to the low-energy structures (LES). The fork is not caused by quantum-interference effects, but, rather, stems from the kinematical constraints imposed upon the electron's returning energies. This is particularly important for orbits whose excursion times are much longer than a field cycle. This explains their absence for few-cycle pulses. Specifically, SFA studies provide a detailed analysis of the role of forward- and backscattered trajectories in the energy ranges for which the fork appears, and show that it is closely related to the LES <cit.>. This raises the question of whether the presence of the Coulomb potential is a necessary condition for the LES to arise. Therein, rhombi-like structures of similar origin have also been identified.  To the present date, there is no study of the influence of the residual Coulomb potential on the fork- or rhombi-shaped structure. Because, however, the fork occurs for momenta nearly perpendicular to the field-polarization axis, the phase distortions introduced by the Coulomb potential when the electron is in the continuum are expected to cancel out <cit.>.  This very likely explains the excellent agreement of the SFA with existing experiments in which the fork is observed.




§ THEORETICAL METHODS


In order to compute ATI spectra or electron momentum distributions, one must determine the evolution of an electronic wave packet from a bound state to a continuum state associated to a final momentum 𝐩_f when the electron reaches the detector. This evolution is given by the time-dependent Schrödinger equation, which can either be solved numerically, or approximately. For one-electron atoms, the ab-initio solution of the TDSE is meanwhile  straightforward, and there are well-established Schrödinger solvers available to the strong-field community <cit.>. For multielectron systems, however, ab-initio solutions are a formidable task. Key issues are how to model electron-electron correlation, the core dynamics, and the coupling of different degrees of freedom. A key difficulty of a full ab-initio solution is the so called “exponential wall”, i.e., the fact that the numerical effort increases exponentially with regard to the degrees of freedom. This implies that a series of approximations must be made in order to render the problem tractable. Examples are the use of the time-dependent density functional theory <cit.>, trajectory-based grids <cit.>, the time-dependent restricted-active-space self-consistent-field theory with space partition (TD-RASSCF-SP) <cit.>, and the R Matrix with time dependence <cit.>. For discussions of classical, semiclassical, quantum and ab-initio methods for correlated two-electron systems in the context of laser-induced nonsequential double ionization (NSDI) see the review <cit.>.
Here, we will consider single-electron systems and, unless otherwise stated, use atomic units throughout.  

Our starting point will be the Hamiltonian H(t)=H_a+H_I(t), of an electron under the influence of an external laser field and a binding potential. The field-free atomic Hamiltonian is given by 

    H_a=𝐩̂^2/2+V(𝐫̂)

and H_I(t) gives the interaction with the external field. In Eq. (<ref>), 𝐫̂ and 𝐩̂ denote the position and momentum operators, respectively. The binding potential  V(𝐫̂) is chosen to be of Coulomb type, i.e.,  

    V(𝐫̂)=-C/√(𝐫̂·𝐫̂),
 
where 0≤ C≤ 1 is an effective coupling. The interaction Hamiltonian in the length and velocity gauge read

    H_I(t)=-𝐫̂·𝐄(t)

and

    H_I(t)=𝐩̂·𝐀(t)+𝐀^2/2,

respectively, where 𝐄(t)=-d𝐀(t)/dt
is the electric field of the external laser field and 𝐀(t) the corresponding vector potential.  The evolution of the system is described by the time-dependent Schrödinger equation

    i∂_t|ψ(t)⟩=H(t)|ψ(t)⟩ ,

which is either solved fully numerically in a specific basis set or computed semi-analytically using a series of approximations. Some of them will be discussed below, but our emphasis will be the strong field approximation and beyond, more specifically the Coulomb Quantum-orbit Strong-Field Approximation (CQSFA). 
It is instructive to write the time-dependent Schrödinger equation (<ref>) in integral form using time evolution operators. This leads to 
    U(t,t_0)=U_a(t,t_0)-i∫^t_t_0U(t,t^')H_I(t^')U_a(t^',t_0)dt^' 
    ,

where the time evolution operators U_a(t,t_0)=exp[iH_a (t-t_0)] and 

    U(t,t_0)=𝒯exp [i ∫^t_t_0H(t^')dt^'],

where 𝒯 denotes time-ordering, are associated with the field-free and full Hamiltonians, respectively, evolving from an initial time t_0 to a final time t.  Above-threshold ionization requires computing the transition amplitude ⟨ψ_p_f(t) |U(t,t_0) |ψ_0⟩ from a bound state |ψ _0⟩ to a final continuum state |ψ_p_f(t)⟩ with momentum 𝐩_f, which may be written as 

    M(𝐩_f)=-i lim_t→∞∫_-∞^t d
    t^'⟨ψ_p_f(t)
    |U(t,t^')H_I(t^')| ψ _0(t^')⟩ 
    ,

with | ψ _0(t^')⟩=exp[iI_pt']|ψ _0⟩, where I_p is the ionization potential. One should note that Eq. (<ref>) is formally exact, and it is used to construct many of the approaches presented below.  Perturbation theory with regard to the laser field is obtained by iterating Eq. (<ref>). This specific series works well for weak fields, but breaks down in the parameter range of interest. Below we will briefly discuss the theoretical orbit-based methods whose starting point is Eq. (<ref>) approximating the time-evolution operator (<ref>). We will examine the SFA and CQSFA in more depth, as these are the approaches that lead to our main results and to our publications that are revised in this work. 




 §.§ Strong-field approximation


The strong-field approximation, also known as the Keldysh-Faisal-Reiss (KFR) <cit.> theory, is the most widespread semi-analytical approach employed in the modelling of strong-field ionization. In its standard form, it consists in replacing the full time evolution operator (<ref>) by the Gordon-Volkov <cit.> time-evolution operator U^(GV) related to the Hamiltonian 

    H^(GV)(t)=𝐩̂^2/2+H_I(t)

of a free particle in the presence of the laser field.  This brings a key advantage as the time dependent Schrödinger equation using the Hamiltonian (<ref>) can be solved analytically.  This procedure will lead to the SFA transition amplitude 
for the ATI direct electrons.

The SFA may also be generalized to account for rescattering by proceeding as follows. One commences by noticing that the full time-evolution operator may also be written as 

    U(t,t_0)=U^(GV)(t,t_0)-i∫^t_t_0U(t,t^')VU^(GV)(t^',t_0)dt^' 
    ,

in terms of the Gordon-Volkov time evolution operator. If one now iterates Eq. (<ref>) to first order in V and inserts the resulting expression in Eq. (<ref>), the zeroth and the first-order term in the series will give the direct and the rescattered electrons when the resulting time evolution operator is employed in (<ref>). It is commonly accepted that the “rescattered electrons” contain one act of rescattering. One may also write the direct ATI transition amplitudes in terms of the binding potential V. Explicitly, this yields

    M_d(𝐩)=-i ∫_-∞^∞dt' ⟨𝐩+𝐀(t')|V|ψ_0⟩ e^i S_d(𝐩,t,t'),

where

    S_d(𝐩,t,t')=-lim_t→∞1/2∫^t_t'[𝐩+𝐀(τ)]^2 dτ + I_pt'

is the action describing a process in which an electron is freed in the continuum by tunnel ionization and reaches the detector without further interaction.  For a thorough explanation and detailed derivations we refer to <cit.>. 
One may also combine the direct and rescattered ATI transition amplitudes in order to obtain the expression

    M_r(𝐩)=- ∫^∞_-∞ dt ∫^t_-∞ dt'∫ d^3kexp[iS_r(𝐩,𝐤,t,t')]V_𝐤0V_𝐩𝐤,

which allows for up to one act of rescattering and contains, apart from the ionization time t' and the final time t, the intermediate momentum 𝐤 as integration variables. Thereby, the action reads

    S_r(𝐩,𝐤,t,t')=
       -1/2∫^∞_t[𝐩+𝐀(τ)]^2 dτ -1/2∫^t_t'[𝐤+𝐀(τ)]^2 dτ + I_pt'.

In the SFA, the influence of the core is incorporated in the ionization prefactor

    V_𝐤0=⟨𝐤+𝐀(t)|V | ψ_0 ⟩

and in the rescattering prefactor 

    V_𝐩𝐤=⟨𝐩+𝐀(t)|V | 𝐤+𝐀(t) ⟩.

Although Eq. (<ref>) also contains direct electrons, it is more convenient to employ the transition amplitude (<ref>) in this case. Calculating the direct electrons from Eq. (<ref>) would require the computation of limits and is thus more cumbersome. This is referred to by some groups as ISFA (improved SFA).

A very intuitive description can be obtained if Eqs. (<ref>) and (<ref>) are solved using saddle-point methods. This requires seeking the values of the ionization times t', the rescattering times t and of the intermediate momentum 𝐤 so that the actions (<ref>) and (<ref>) are stationary, i.e., their derivatives with regard to such variables must vanish. The saddle-point equation obtained from the direct action S_d(𝐩,t') reads

    [𝐩+𝐀(t')]^2 = - 2I_pt',

which gives the kinetic energy conservation at the ionization time. Eq. (<ref>) has no real solutions, as a direct consequence of tunnel ionization having no classical counterpart. If the condition ∂ S_r(𝐩,𝐤,t,t')/∂ t'=0 is imposed upon (<ref>), a formally identical equation will follow, with the final momentum 𝐩 being replaced by the intermediate momentum 𝐤. This gives

    [𝐤+𝐀(t')]^2 = - 2I_pt'.

The condition ∂ S_r(𝐩,𝐤,t,t')/∂ t=0  leads to the conservation of energy

    [𝐩+𝐀(t)]^2=[𝐤+𝐀(t)]^2

upon recollision, and ∂ S_r(𝐩,𝐤,t,t')/∂𝐤=0 restricts the electron's intermediate momentum  𝐤, such that it must return to the site of its release, i.e., the origin. Explicitly, 

    𝐤=-1/t-t'∫_t^t'𝐀(τ)dτ.

Physically, the approximations introduced above have the following consequences:

	
  * The influence of the laser field is neglected when the electron is bound. This means that the SFA does not take into consideration processes involving bound states only, such as excitation, the core dynamics, or field-induced distortions such as Stark shifts. It is also quite common to neglect bound-state depletion. One may however modify the SFA in order to include such effects to a certain extent. Bound-state polarization has for instance been introduced in <cit.>. Furthermore, electron-electron correlation and excitation have been incorporated in our previous work <cit.> and <cit.>, respectively. For a thorough review on these issues see <cit.>.
	
  * The influence of the residual binding potential is neglected when the electron is in the continuum.  This implies that the field-dressed momentum is conserved, and that the influence of the core is incorporated via prefactors at a single point, namely the origin. Relaxing this approximation leads to several Coulomb-distorted approaches, whose overview will be provided here.
	
  * Formally, the SFA may be viewed as the Born series with a modified, field-dressed basis. This allows a clear distinction between direct and rescattered events.  This becomes even clearer if one looks at the saddle-point equations stated above and the corresponding actions. Eqs. (<ref>) and (<ref>), in principle, allow return, but not rescattering, while Eqs. (<ref>), (<ref>), (<ref>) and (<ref>) allow a rescattering event to occur. 
	
  * If the steepest descent method is used, one gains an intuitive orbit-based interpretation but loses the spatial and momentum widths at the instant of ionization and recombination that would stem from the uncertainty relation. Relaxing these approximations has shown some softening on phase jumps in molecular HHG <cit.> and quantitative differences in the HHG spectra <cit.>.


One should note, however, that there exists a more general formulation of the SFA which takes into consideration exact scattering waves instead of Volkov states. This formulation has been developed in <cit.> and is reviewed in <cit.>.



 §.§ Beyond the strong-field approximation




In this section, we summarize several orbit-based theoretical methods and models that incorporate Coulomb effects in the continuum.  Some of them, such as the Coulomb Volkov approximation (Sec. <ref>) and the Quantitative Rescattering Theory (Sec. <ref>), rely on factorization and modified continuum states, but employs orbits from the SFA, either direct or rescattered. The low-frequency approximation (Sec. <ref>) is developed along similar lines, but employs a higher-order Born expansion.
Other methods either include Coulomb distortions in the orbits approximately or perturbatively. For instance, the Coulomb-corrected Strong-Field Approximation (CCSFA) (Sec. <ref>) expand the orbits perturbatively around the SFA, while the Eikonal Volkov Approximation (EVA) (Sec. <ref>), employ a field-dressed Wentzel-Krammers-Brillouin (WKB) and further approximations in the scattering angles. This makes it possible to compute the orbits recursively starting from the SFA, and makes the EVA suitable for the outer region in the Analytical R-Matrix (ARM) theory (Sec. <ref>). One may also treat the Coulomb potential and the laser field on equal footing by solving classical equations of motion in which both are present in order to describe the continuum propagation, but incorporating these equations in the action, so that quantum interference is included. This is the procedure taken in the Trajectory-Based Coulomb SFA (TCSFA) (Sec. <ref>), which is similar to the CCSFA, but non-perturbative with regard to the residual potential, and in path-integral strong-field approaches such as the Quantum Trajectory Monte Carlo (QTMC) method (Sec. <ref>), the semiclassical two-step model (SCTS), the semiclassical approximation for strong-field processes, and the Coulomb-Quantum Orbit Strong-Field Approximation (CQSFA). The latter is the method employed by us and will be discussed in more detail in Sec. <ref>. In order to facilitate the discussion, a brief timeline with these methods is shown in Fig. <ref>. 

Further approximations and differences in implementation exist in all the above-mentioned methods, the most important of which are: 

	
  * Real orbits in the continuum. All methods with full Coulomb distortion in the orbits  (see dark blue boxes in Fig. <ref>), except a very recent version of the CQSFA, employ real orbits in order to describe the electron's continuum propagation. This is due to the fact that complex orbits lead to branch cuts, which are tractable if the orbits are Coulomb free or if approximate methods are employed around the SFA, but are much more difficult to tackle otherwise. 
	
  *  Adiabatic tunneling rates. Because sub-barrier dynamics are not easy to address if Coulomb-corrected or Coulomb-distorted orbits are taken (sky blue and dark blue boxes in Fig. <ref>, respectively), some methods, such as the QTMC and the SCTS employ adiabatic, cycle-averaged tunneling rates to weight the initial orbit distributions (see dark blue boxes in Fig. <ref> outlined in green).  Methods that go beyond that must make several approximations in order to deal with this part of the problem, which includes cusps and singularities. A proper treatment of sub-barrier dynamics also requires complex orbits, and the above-mentioned tunneling rates provide a way to avoid them. 
	
  * Direct vs inverse problem. The full presence of the Coulomb potential in the continuum significantly alters the topology of the orbits. This makes it difficult to guess their shapes, relevance or main features. For that reason, methods such as the TCSFA, QTMC and SCTS solve the direct problem, i.e., given the initial momentum distribution, they employ a shooting method and a large number of contributed orbits (typically 10^8-10^9). Once these orbits have been propagated in time, they are binned according to the resulting final momentum. An advantage of this type of implementation is that it is versatile, as it can be easily adapted to different targets and driving-field shapes. Furthermore, the initial momentum distribution is also known. Two major disadvantages are however the presence of cusps and caustics if sub-barrier corrections are included, and the huge number of contributed orbits. Another possibility is to solve the inverse problem, i.e., given a final momentum, find the corresponding initial momentum assuming a specific type of orbit and sub-barrier corrections. This approach is typically used if the Coulomb potential is only incorporated approximately around the SFA (see sky blue boxes in Fig. <ref>), and, for fully Coulomb-distorted orbits, in the CQSFA. Advantages of solving the inverse problem are: (i) a much smaller number of contributed orbits, typically four for fixed final momentum components; (ii) “cleaner” interference patterns and a better control over cusps and caustics, which facilitates the inclusion of sub-barrier corrections and complex trajectories. A drawbacks is the loss of versatility, i.e., prior knowledge of the approximate behavior of the orbits and precise initial guess conditions are required. This is because the inverse problem is harder to solve.

One should bear in mind that this is only a short summary.
A discussion of these methods, focused on their advantages, shortcomings, approximations, implementation and on their success in reproducing specific holographic structures, is provided below. For a summary of the the features included in many of these models see Table <ref> presented subsequent to the discussion.



  §.§.§ Coulomb-Volkov approximation (CVA)

	

In order to account for Coulomb effects in the continuum one approach one may replace the Gordon-Volkov states with so called Coulomb-Volkov states, which approximately account for both the laser field and the Coulomb potential. To make the problem tractable it is assumed that the laser distorts the Coulomb continuum adiabatically. Then, the resultant Coulomb-Volkov wavefunction is the product of a Coulomb continuum scattering state and a Gordon-Volkov state. Explicitly, the final state of the system reads

    ⟨𝐫|ψ^(CV)_𝐩(t)|=⟩ψ^(CV)_𝐩(𝐫,t)=ψ^(C)_𝐩(𝐫)ψ^(GV)_𝐩(𝐫,t),
 
where ψ^(C)_𝐩 is the field-free Coulomb scattering wave for the (final) continuum state and 

    ⟨𝐫|ψ^(GV)_𝐩(t)|   ⟩=ψ^(GV)_𝐩(𝐫,t)
       =exp[i [𝐩+𝐀(t)]·𝐫]/(2π)^3/2expiS_d(t,t'),

is the Gordon-Volkov wavefunction. In (<ref>),  S_d(t,t') is given by the SFA action for direct ATI [Eq. (<ref>)]. This ansatz is then used in the direct ATI transition amplitude. The above-stated equations show that, in the CVA, the propagation from the instant of ionization to the instant of detection is Coulomb-free, as it is dictated by the Volkov propagator. 
A summary of the history and main results of the CVA can be found in the review <cit.>.
This method was originally proposed in a strong field context in Refs. <cit.> and since then there has been a wealth of publications further developing it <cit.>.

One of the early successes of the CVA was that it was found to break a four-fold symmetry in the momentum distributions for photoelectrons ionized via elliptically polarized fields <cit.> that was exhibited by SFA-like approaches but was not observed in experiments <cit.>. In the case of relatively high frequency and low intensity linearly polarised laser fields it was found <cit.> that the CVA had markedly better agreement with ab-initio solutions than the SFA (left panels in Fig. <ref>) for low photoelectron momenta in the multiphoton regime. 

On the other hand, this agreement worsens for higher photoelectron momenta or lower-frequency fields (right panels in Fig. <ref>). In fact, because the continuum propagation uses Volkov states, the CVA is only suitable for a parameter range in which significant Coulomb distortion of the continuum propagation is not required. Thus, it is not surprising that the CVA reproduces the fan in the multiphoton regime, as the Coulomb factor in Eq. (<ref>) can describe a resonant behavior appropriately. However, the CVA does not satisfactorily reproduce interference effects in the low-frequency regime, as in this case the Coulomb potential would exert a strong influence on the electron continuum propagation. In the right panels of Fig. <ref>, one clearly sees that, in this latter caase, the CVA neither reproduces the correct number of fringes in the near-threshold fan, nor the spider-like structure, which is almost absent. 	

Recently however, an alternative version of the CVA <cit.>, known as the second order CVA or CVA^2 has been used to account for rescattering by using the Born expansion, introduced in the above section, replacing the scattering states in the SFA with Coulomb-Volkov states. This enables the CVA to account for hard-rescattering in the same way as the SFA and extends its description to the high energy photoelectron momentum distributions. Care should be taken with this method, as Coulomb-Volkov states are not solutions of the Schrödinger equation and their use will introduce errors in a non-standard way. 
Additionally, there is some inconsistency using Coulomb-Volkov states in direct terms if a Born series about the potential has been applied. The direct terms will contain only a Gordon-Volkov time evolution operator so should not contain any Coulomb effects, only the rescattered term will contain the full time evolution operator that will warrant this. 	
Hence, it is not clear what effects are neglected and some parts may be “doubly”/ over accounted for. However, as the original CVA did not reproduce effects due to rescattering electrons there is an empirical reason to construct such a model and “doubly” counted effects will most likely be small.

In its original form the CVA was not suited for use in electron holography, adding no real improvement over the SFA for the typical parameter range. Furthermore, the lack of quantum trajectories/orbits made it difficult to interpret any interference holographically. In its current form, the CVA has been used to study some holographic effects <cit.> in some cases exploiting the new-found ability to account for rescattering. However, the CVA still tends to closely mimic the SFA for the low frequency fields used in holography (see Fig. 5 in Ref. <cit.>). This prevents a good description of effects that deviate non-perturbatively from the `traditional' SFA such as the LES[In fact, it has been shown <cit.> that, by considering previously neglected forward scattered trajectories, the LES can be reproduced within the SFA. Thus, if this were applied to the CVA, presumably it would be able to reproduce the LES as well.] <cit.>, which the CVA does not correctly reproduce. Thus, the CVA can be used to interpret and understand holographic interference, but mainly by exploiting and improving the mechanisms that already qualitatively exist in the SFA.





  §.§.§ Quantitative rescattering theory (QRS).



The QRS, also referred to as scattering-wave SFA (SW-SFA) in earlier papers   <cit.>, uses a particular factorization of the momentum dependent transition amplitude. This can be employed to describe processes in which electron recollision or recombination occurs, such as HHG, high-order above-threshold ionization (HATI) and NSDI. For more information on the QRS method and its current status see the review <cit.> and the tutorial <cit.>.  With this factorization the momentum dependent probability distribution is then written as

    M(𝐩)=W(p^(r))σ(p^(r),θ),

where σ(p^(r),θ) is the differential elastic scattering cross section between free electrons with momentum  𝐩^(r)  at the instant of rescattering, and scattering angle θ and the ion, while W(p^(r)) describes the flux of a returning wavepacket. It is shown that the flux term can be treated accurately by the SFA, while the dipole term uses a Coulomb scattering wave to account for effect of the Coulomb attraction of the ion. The final momentum 𝐩 is given in terms of the momentum of the returning wavepacket by 

    𝐩=𝐩^(r)-𝐀(t^(r)),

where t^(r) denotes the rescattering time.
This factorization was tested both theoretically <cit.> and experimentally <cit.> and found to work well in key cases for atoms <cit.> and molecules <cit.>, in comparison with ab-intio solutions. For an example of a photoelectron spectra for ATI calculated with the QRS see Fig. <ref>, where it is compared with ab-initio methods. Following this success it was used to explain the species dependence of the high-energy plateau in HHG <cit.> and HATI <cit.> and was applied to HHG of small <cit.> and even polyatomic <cit.> molecules  where ab-intio solutions would be too intensive. Despite this, it has been stated that <cit.> using the QRS to model HHG for complex polyatomic molecules is troublesome as it is a single active electron model that does not account for the electron or nuclear correlations. Recently, in Ref. <cit.>, it was argued that for CF_3I there are contributions from many orbitals for the photoelectron distributions. Thus, single active electron models, where only the highest-occupied molecular orbitals are considered (such as the QRS) will not be successful.

The QRS method works well for the range of intensities of interest for photoelectron holography. However, in contrast to CVA it does not describe the direct electrons in ATI and thus is not suited to study most holographic interference patterns as it will not account for the reference signal. Despite employing a trajectory-based approach for the returning wavepacket this description uses SFA orbits. Thus the Coulomb potential is neglected in the continuum propagation besides hard scattering events. Then, although this method can certainly be used for imaging purposes such in HHG spectroscopy (e.g. <cit.>) or using a LIED technique in ATI <cit.>, the target information will not be encoded in the trajectories but in the scattering cross-section term. 




  §.§.§ Low-Frequency Approximation (LFA).
	


The idea behind the LFA is to factorize the transition amplitude in a similar fashion to the QRS. However, this is  achieved by extending the Born series expansion about the potential to beyond first order. A early form of the LFA was proposed in <cit.> but it was developed and used in a strong field context in the following publications <cit.>. Eqs. (<ref>) and (<ref>) give the zeroth and first-order terms in the Born expansion, which are usually associated with direct and rescattered electrons, respectively. The inclusion of the first-order term is sometimes referred to as the improved SFA (ISFA) or SFA2. 
The Born series is obtained by iteratively substituting Eq. (<ref>)[Note that an alternative form of Eq. (<ref>) is used for LFA, where the full and Gordon-Volkov time evolution operators are swapped in the integrand.] into Eq. (<ref>). Both 
equations are formally exact and an approximation is only made when, after iteration, the remaining full time-evolution operators are replaced by their Gordon-Volkov counterparts. 
Repeating this process twice and retaining the full time time evolution operator adds an additional “correction” term to the standard rescattering transition amplitude, namely


    M^(2)_r(𝐩)   =(-i)^2lim_t→∞∫_-∞^t dt'∫ d^3𝐤∫_t'^t dt”[V_𝐩𝐤V_𝐤0e^iS_r(𝐩,𝐤,t,t')
           -i∫_t”^t dt”'⟨ψ^(GV)_𝐩(t”')|VU(t”',t”)V|ψ^(GV)_𝐤(t”)|⟩_M^(c)(𝐩,𝐤)⟨ψ^(GV)_𝐤(t')|H_I(t')|ψ_0(t')|⟩],


where the Gordon-Volkov states are defined in Eq. (<ref>).
The extra term can be interpreted as describing two rescattering events. Including it makes the total transition amplitude formally exact. Now as in <cit.>, one of the seminal papers, the LFA is used to further evaluate M^(c)(𝐩,𝐤). The assumption is made that, if the laser intensity is high and its  frequency is low, the laser-dressed momentum will change only a little between rescattering events.  This means that one may use a stationary approach in the second-order term in Eq (<ref>). Explicitly,  the time evolution operator U(t”',t”) is replaced by the time-independent Green's operator G_V(E_𝐤+𝐀(t”'))δ(t”'-t”) calculated at time  t”' for the field-dressed instantaneous energy  

    E_𝐤+𝐀(t”')=1/2[𝐤+𝐀(t”')]^2.


Using this approximation, it is possible to write the whole rescattered contribution, in a similar fashion to the QRS, as a product of a differential cross section of the laser-free continuum electron back to the initial bound state and a returning electron wavepacket part. This procedure changes the scattering states but leads to the same orbits as the improved SFA as given by Eq. (<ref>), and the resulting transition amplitude can be evaluated using the saddle-point or uniform approximation.

In <cit.>, it was argued that the QRS is equivalent to taking only a single (short) orbit in the above expression. It has even been stated <cit.> that the QRS can be derived entirely from the LFA by making some additional approximations. However, an exact correspondence between these two methods is not completely clear as the cross-section in the QRS is between a continuum-scattering state (accounting for the Coulomb potential) and the ion, which is evaluated exactly using ab-initio techniques. However, in the LFA the cross-section is between Gordon-Volkov states and the initial atomic state. A further difference with regard to the QRS is that the transition amplitude in the LFA is in integral form, due to applying the relations for the time-evolution operator. This  integral is then evaluated using the saddle-point approximation.


There has been extensive comparison between the LFA and the ISFA with improvements found both in the high-energy ATI spectra <cit.> and in the low-energy part when studying the LES <cit.>. Reasonable agreement has been found when compared to experiment for the ATI spectra of atoms <cit.>. A comparison with ab-initio methods has been performed for the case of photodetachment of F^- <cit.>, see Fig. <ref>, and in the extension of the LFA to HHG <cit.> with bicircular laser fields. In both cases very good agreement with the TDSE was found. The LFA has also been extended to beyond the dipole approximation for ATI of H_2^+, where good agreement was found along the momentum components parallel to the laser field polarization when compared to ab-initio methods <cit.>.

The LFA is a fully fledged quantum-orbit method and unlike the QRS it has been shown to give improvements in the low and high-energy part of the ATI spectrum. For these reasons, it seems appropriate to apply this method to electron holography. Despite this, the LFA does not satisfactorily reproduce the holographic interference of interest that arises for linearly polarized fields. This may be because only hard-scattering Coulomb effects are considered and Coulomb-distortion of the orbits is not incorporated. These distortions include deflections and soft collisions, which may happen away from the origin. For linearly polarized fields these “soft” Coulomb effects have been shown to be crucial for the correct formation of interference structures present in the distributions <cit.>. The Born series will in theory, if taken to high enough order, describe all Coulomb effects. However, the direct orbits in Born-series approaches neglect Coulomb distortion and deflection. Instead, these effects will be described in terms of the higher order terms that relate to orbits undergoing multiple hard recollisions. This is an inappropriate basis with which to describe “soft” Coulomb effects and as such should not be expected to converge quickly for cases where this plays an important role.



  §.§.§ Coulomb-Corrected SFA (CCSFA).



The CCSFA <cit.> introduces Coulomb effects perturbatively to the SFA in the action. For a detailed introduction to the CCSFA see the review <cit.>. The peturbative expansion is implemented  after the application of the saddle-point approximation to Eq. (<ref>). Both the action and the trajectories can be expanded in the following way

    𝐫(𝐩_f,t)   =𝐫^(0)(𝐩_f,t)+𝐫^(1)(𝐩_f,t)+…
    𝐩(𝐩_f,t)   =𝐩^(0)(𝐩_f,t)+𝐩^(1)(𝐩_f,t)+…
    
    S(𝐩_f,t)   =S^(0)(𝐩_f,t)+S^(1)(𝐩_f,t)+…,

where the zeroth order gives the Coulomb-free trajectories and the direct SFA action [Eq. (<ref>)].
The first order correction to the action is then given by

    S^(1)(𝐩_f,t)=∫_t_s^t dτC/|𝐫^(0)(𝐩_f,τ) |,

where a hydrogenic Coulomb potential has been used and C is the Coulomb correction factor, set to unity for singly charged ions. The corrections are all dependent on the final momentum, thus the CCSFA uses an inverse solution approach. The final transition amplitude can be written as

    M(𝐩_f)=∑_s 𝒞(t_s)exp(i S(𝐩_f,t_s) + i F(𝐩_f,t_s) ),
 where t_s are the corrected times of ionisation given by inserting a corrected initial momentum 𝐩_0 into to the solutions of the saddle-point equation (<ref>),  𝒞(t_s) is the SFA prefactor, 

    F(𝐩_f,t)=𝐫·[𝐩_f+𝐀(τ)]|_t_s^t≈ F^(0)(𝐩_f,t)± C/√(2 I_p),
 
The Coulomb phase is integrated along the whole temporal contour, even up to the singularity. A regularization procedure has been introduced to remove the singularity, in which the asymptotic form of the initial bound state wavefunction is matched to the action for a matching time t_*, where the laser is dominant over the Coulomb force but only a fraction of a laser cycle has taken place since ionization. A detailed formulation of the procedure is given in <cit.>. Bound-state regularization was also applied to other methods that use a Coulomb phase such as the EVA <cit.> and the CQSFA <cit.> discussed below. 

This method is in theory very quick to calculate. However, a particularly debilitating issue that affects it (and all methods that employ a Coulomb phase) is branch cuts, which arise due to the complex trajectories in the square root of Eq. (<ref>). This causes defects in the momentum distributions if crossed in the temporal integration contour. Recently, algorithms that avoid the branch cuts have been introduced <cit.>, these work well for such perturbative methods, enabling 2D photoelectron momentum distributions to be produced without defects or additional approximations. (For an example see Fig. <ref>).  This will be discussed in more detail later in this review in Secs. <ref>, <ref> and <ref>.
The CCSFA, together with branch-cut corrections, was successfully used to explain orders of magnitude difference between ab-initio methods and the SFA around the direct ATI cutoff, in terms of soft recollisions <cit.>. This description required complex trajectories so that the overall probability could be modified during recollision.

The CCSFA provides a consistent and clear extension of the SFA but as a perturbative extension of the direct SFA there are many effects it can not account for. One obvious extension is to solve Newton's equations for the trajectories, including both the laser and Coulomb interaction on an equal footing, instead of using a perturbative expansion. This extension of the CCSFA is referred to as the TCSFA and will be discussed next. Another option could be to, similar to the CVA, include the SFA rescattered electrons. However, using two alternative perturbative series could make the method lose consistency and it may be unclear what effects are included or left out. Recently, the CCSFA has been extended to describe HHG <cit.>, which shows it is possible to incorporate the returning wavepacket in this method. 



  §.§.§ Trajectory-Based Coulomb SFA (TCSFA)



In the TCSFA <cit.> the CCSFA is built on by including equations of motion that include the Coulomb potential for continuum trajectories. Thus the trajectories are now described by

    𝐫̇(τ)=𝐩+𝐀(τ)       𝐩̇(τ)=-∇ V(𝐫(τ)).

Solving these equations of motion in conjunction with the tunnel ionization [Eq. (<ref>)] led to a new categorization of quantum orbits <cit.> (later also used in the CQSFA) of four possible contributing trajectories for each final momentum point. Methods that solve the full Newton's equations of motion including the Coulomb potential are some times called Coulomb-distorted trajectory approaches. This approach led to unprecedented qualitative agreement with ab-initio methods for the direct ATI 2D momentum distributions. In an additional publication the TCSFA was used to show the importance of including the Coulomb phase in the action for the tunnelling part of the trajectory, known as sub-barrier Coulomb-corrections (sub-CC), see Fig. <ref>. It was shown that as well as influencing the overall tunnelling probability, the interference patterns in momentum distributions were also affected by sub-CC <cit.>. Recently, a version of the TCSFA was formulated that goes beyond the dipole approximation <cit.>.




In contrast to the CCSFA and CQSFA, the TCSFA solves the direct problem by specifying the initial rather than final momentum in a shooting method. The resulting final momenta are binned  similar to some Monte Carlo approaches such as the QTMC and SCTS methods, see Table <ref>. Because the momentum is not conserved, an extra 𝐩̇·𝐫 term must be incorporated in the action, as it is written in phase space. This term is not included in the TCSFA. It has been conjectured <cit.> that this results in some discrepancies when compared with ab-initio solutions, such as the wrong number of fringes appearing on the inner ATI ring for the fan-like structure <cit.>. Branch cuts are an even bigger problem for Coulomb-distorted trajectory models because of the complex square root in equations of motion. There are no full solutions to this branch cut problem, but a partial one has been suggested in <cit.>. In the TCSFA (and most other models of this type), branch cuts are avoided by taking real trajectories in the continuum. This is only a partial remedy  because, as shown in Refs. <cit.>, complex trajectories are essential for capturing certain physical effects, such as the effective deceleration of the electron wavepacket due to the Coulomb potential. In principle, the TCSFA could be used for holographic imaging. However, the caustics in the 2D momentum distributions may cause some difficulties. Additionally, the binning of trajectories into final momenta makes it is more difficult to isolate particular interference effects than other methods. On the other hand, the binning approach means that classes of solutions for the trajectories do not have to be known in advance. This makes it much easier to generalize the TCSFA to more complex systems and/or driving fields, for which the classes of solutions for the trajectories may change significantly.



  §.§.§ The Eikonal-Volkov Approximation (EVA).
 


The EVA <cit.> uses a totally different approach to most quantum trajectory methods. The Coulomb field is accounted for by corrections to Gordon-Volkov states using
the eikonal approximation. The wavefunction is written using the WKB (Wentzel-Kramers-Brilliouin) ansatz

    ψ(𝐫,t)=P(𝐫,t)exp(i S(𝐫,t)/ħ).

Here S(𝐩,t)=S^GV_𝐤(𝐫,t)+G_k(𝐫,t), where the first term is the phase of a Gordon-Volkov and the second is a correction/ distortion due to the Coulomb potential. This is then inserted into the Schrödinger equation, expanded about ħ and terms of the same order are collected. This results in an equation that determines G_k, which can be solved to integral form 

    G_k(𝐫,t)=-∫_T^t dτ V(𝐫_0(τ))+G_0𝐤(𝐫_0(τ)),

where 𝐫_0 denotes the Coulomb-free electronic coordinate in the laser field. 
Note that, in order to get to Eq. (<ref>), the approximation that the change in electron momentum during scattering is small i.e. |∇ G_𝐤(𝐫,t)|≪|𝐤| is used. When an electron travels close to the core, in many cases, the large change in momentum will cause this condition to fail. However, in conjunction with this method the analytical R-Matrix (ARM) method was developed, which enables the use of an alternative approach for electrons that come close to the core. This is discussed in the next section.

Using the EVA one can construct continuum states |𝐩^EVA_T(t)⟩, that are defined by backpropagating the field-free eikonal states from after the laser is off at time T to the desired time t. These continuum states in position representation are given by


    ⟨𝐫|𝐩^EVA_T(t)|=⟩1/(2π)^3/2   exp[i(𝐩+𝐀(t))·𝐫
    -i/2∫_T^tdτ(𝐩+𝐀)^2
    .
       .
    -i∫_T^tdτ V(𝐫_0(τ))+iG_0𝐩(𝐫_0(T))
    ].





  §.§.§ Analytical R-Matrix (ARM) Method



The ARM method <cit.> splits the space into an inner and outer region, near and away from the atom/ molecule, respectively. In the inner region, exact ab-initio methods are used, while in the outer region the EVA is employed. This has the capacity to give very accurate results, while retaining the concept of quantum trajectories.
The Hamiltonian can be split into two parts that describe the dynamics in  the inner and outer regions, denoted + and -, respectively, given by ℋ̂^±=Ĥ+L̂^±(a), where L̂^± is the Bloch operator, which preserves the hermiticity of the both the regions and is defined in position notation as

    L̂^±=±δ(r-a)(d/d r+1-b/r),

where a is the radius of the circular boundary which separates the two region and b is an arbitrary constant. In Ref. <cit.>  b is ultimately set to b=C/√(2I_p), which is unity in the case of hydrogen. 
A solution of the TDSE for the inner and outer region may be written  integral form

    ψ^±(𝐫,t)=-∫_-∞^t d t' ∫ d^3 𝐫' G^±(𝐫,t,𝐫',t')L^±ψ(𝐫',t'),

where G^±(𝐫,t,𝐫',t')=⟨𝐫|U^±(t,t')|𝐫'|$⟩ is the Green's function in the inner or outer region. The initial wavefunctionψ(𝐫',t')can be taken as the atomic ground state before the laser is turned on. This can be inserted into the outer equation to describe the ionization process

    ψ^-_(1)(𝐫,t)   =-∫_-∞^t d t' ∫ d^3 𝐫' G^-(𝐫,t,𝐫',t')
    
       ×
    L^-ψ_0(𝐫')exp(iI_p t').

To describe a wavepacket returning to the region of the coreψ^-_(1)may be substituted into the inner equation to getψ^+_(2), then finally the scattered wavepacket is described by insertingψ^+_(2)into the outer equation to getψ^-_(3). In theory this process could be repeated for further recollision processes but the parameterais chosen such that return to the inner region results in a “hard” recollision and thus the electron wavepacket is unlikely to return. The saddle point approximation is applied to the integrals stemming from Eq. (<ref>), resulting in a complex time of ionization in common with other quantum orbit approaches. For the propagation in the outer region the EVA is exploited to make this region easy to solve. The Green's function for the outer region is replaced with one constructed via the EVA

    G^-(𝐫,t,𝐫',t')   =⟨𝐫|U^EVA(t,t')|𝐫'|
    ⟩   =θ(t-t')∫ d𝐩⟨𝐫|𝐩^EVA_T(t)|⟨%s|%s⟩⟩𝐩^EVA_T(t')|𝐫' ,

where the states|𝐩^EVA_T(t)⟩are given by Eq. (<ref>).
One of the main benefits of the ARM method is that it incorporates a rigorous approach to “hard” scattering, something that is missing from many approaches such as the CCSFA, TCSFA and CQSFA. Additionally, it is easy to extend to many electron systems, shown in <cit.>, as the R-matrix method is already designed for this.

The combination of the EVA and ARM has given very good results. For example it has been used in the case of nearly circularly polarized fields leading to the interpretation of ionization times in the attoclock[The attoclock exploits nearly circularly polarized fields to use the offset in photoelectron ejection angle away from the peak of the field to give a measure of time taken for the electron to tunnel.]<cit.>. This demonstrates that, for hydrogen, there is no tunnelling delay once effects due to the Coulomb potential have been accounted for. A high level of accuracy was achieved for the ARM computations for nearly circularly polarized light, as, in comparison with three different ab-initio solutions, deviations did not exceed 2%<cit.>. Furthermore, it was used to show that including the imaginary part of the complex trajectories improves agreement with ab-initio solutions of the TDSE and can be related to wave-packet deceleration by the Coulomb potential <cit.>. In a recent set of three companion papers <cit.>, the ARM method is extended to look at the effect of fields with non-vanishing ellipticity upon initial bound electronic states with co- and counter-rotating orbital angular momentum. Thereby, differences in ionization yields, the interplay of the Coulomb potential, its effect on the tunnel exit velocities and coordinate, and the high level of spin polarization possible in the photoelectrons were examined.

Despite the high level of success of the ARM method, it has not been used to calculate 2D photoelectron momentum distributions in order to understand the holographic interference structures present for linearly polarised fields. 
Given that the EVA is ultimately perturbative and assumes small changes in momentum, the strong deflections by the Coulomb potential that are not quite classified as “hard” collision by the ARM method may not be fully accounted for. These intermediate interactions occur much more commonly for linear polarized fields. 
Since the EVA employs complex trajectories, it will also suffer from branch cuts. Unlike the TCSFA and CQSFA, the incomplete solution to take the trajectories to be real in the continuum does not seem to have been used, perhaps to preserve consistency of the method. It seems likely that a proper treatment of branch cuts is the main stumbling block for modelling holographic patterns from linear polarized fields for the ARM method.



  §.§.§ Quantum-Trajectory Monte Carlo (QTMC).


The QTMC <cit.> is an extension upon the classical-trajectory Monte Carlo (CTMC) method <cit.>, which uses the quasi-static Ammosov-Delone-Krainov (ADK) tunnelling theory <cit.> to describe the first step of ionization, while continuum propagation is described by classical trajectories using Newton's equations of motion accounting for both the Coulomb and laser fields. The final momentum distribution is then constructed, via a direct approach, by randomly choosing a spread of initial conditions, in the Monte Carlo fashion. These initial conditions are used to both provide a weight via the ADK tunnelling rate and to compute trajectories in order to calculate the final momenta. Then all the trajectories are binned via their final momentum, such that ADK probabilities are summed if they correspond to trajectories with momenta that lies in a particular bin. 

In the classical case, the probability distribution in thei-th bin corresponding to a momentum𝐩_ifreads

    Ω_ϵ(𝐩_if)   =∑_𝐩_jf∈𝕊^bin_iW(t_0,p_j0⊥),
    
    𝕊^bin_i:   ={𝐩_f : |p_fk-p_ifk|<ϵ  ∀ k∈{∥,⊥}},

whereW(t_0,p_j0⊥)is the ADK rate andϵ≪1. Note that computed trajectories are denoted by the indexj, with final momentum𝐩_jf, initial transverse momentump_j0⊥and tunnelling timet_0. The probability distribution is calculated for a discrete set of momenta𝐩_if, each centred on a bin. Square bins are shown here but in theory any shapes could be used.

The CTMC is improved on by the QTMC by including a phase (the action) with each trajectory, allowing for interference effects. The phase is given by the classical action and no additional fluctuation factor/prefactor term is included. Hence, this corresponds to a zeroth-order expansion of the action in the continuum propagator given by Feynman's path integral around classical saddle-points (a.k.a. the classical limit).

    ⟨𝐩_f|U(t,t_0)|𝐫_0|≈⟩(2π)^-3/2exp(iS[𝐫]-i𝐩_f·𝐫_f).

The probability in the QMTC is given by

    Ω_ϵ(𝐩_if)   =|∑_𝐩_jf∈ S^bin_i√(W(t_j0,p_j0⊥))exp(i S[𝐫_j])|^2
    
    
    S[𝐫_j]   =∫_t_j0^∞d t
    (1/2𝐫̇_j(t)^2+I_p+V(𝐫_j(t))).

Some approximations have been used to arrive at this formula. Firstly, only a real probability is provided by the ADK tunnelling rate formula, so that any phases accumulated in the tunnelling process will be lost. It has been shown in various studies that important phases are picked up in this process (e.g. see <cit.> for a small sample of such work). Secondly, the Monte Carlo approach only approximately solves the saddle-point equation and there is an error dependent on the bin size, something shared by all direct approaches. Given that the use of the saddle-point approximation is motivated by a fast varying action, the bin size may have to be very small to account for this. Also the initial momenta may vary considerably from the final, thus a large spread of initial momenta should be taken. All this means that many trajectories must be used in order to populate the bins and resolve interference features. Thirdly, as shown in <cit.>, there should be a factor of two in front of the potential. In the QTMC, this does not appear because the term𝐩̇·𝐫that should appear in a phase-space action has been neglected. The use of the ADK rate does remove the caustics that were present in the TCSFA. This is because the ionization time and initial conditions will all be real. Thus, the problem of branch cuts is avoided.

The 2D-momentum dependent probability distributions produced by the QTMC have been compared with ab-initio methods in <cit.> and many qualitative features agree. For an example see Fig <ref>. 

However, in the same study <cit.> it was shown that the QTMC does not reproduce many qualitative features in the 2D-momentum dependent probability distributions, if fields with very high ellipticity are used. For instance, the peak of the momentum distribution over the radial momentum component differed by more 50% over a range of intensities. This is due to the adiabatic/quasi-static approximation used in the ADK tunnelling rate. On the other hand, it was found the CCSFA (which is a non-adiabatic method) was able to give very good agreement with solutions of the TDSE. In response to this, in <cit.>  an extension of the QTMC was proposed, where the non-adiabatic features of the SFA were introduced to replace the ADK tunnelling rate. In this new version of the theory, the SFA is used to calculate a tunnel-ionization rate, while the the Monte Carlo approach is used for continuum propagation. This makes it more similar to the TCSFA, but the phases and sub-barrier Coulomb-corrections in the ionization step are neglected <cit.>. Given that the QTMC now uses a complex time, the tunnel exit is also complex. Thus, its imaginary part must be neglected to circumvent the difficulty of branch cuts.

The QTMC has been used to model the fan-like structure <cit.>, the spider-like structure <cit.> and distributions that contain features from both of these <cit.>. Because the QTMC is a trajectory-based description with a phase associated to each path, it is well placed to analyze holographic interference. However, given that a tunnelling rate is used, phases that come from the ionization step or the initial bound state of the atom or molecule are ignored. Such phases may be very important to capture some interference features <cit.>.





  §.§.§ Semiclassical Two-Step Model (SCTS).


In the SCTS <cit.> a similar approach to the QTMC is used, with a tunnelling rate calculated by ADK theory weighting the initial conditions. 
In order to improve the speed of the calculation, this weight is included via importance sampling, where initial conditions are chosen with a biased random distribution instead of a uniform one. This means that less computational time is wasted calculating trajectories with very low probabilities.
Additionally, the SFA is used to provide an initial condition for the initial momentum coordinate parallel to the laser field. Just as in the QTMC, the trajectories are propagated using Newton's equations of motion, the action is included for each trajectory via the classical limit and they are combined via a binning method is of the final momenta. This gives the final momentum-dependent probability distribution as

    Ω_ϵ(𝐩_if)   =|∑_𝐩_jf∈𝕊^bin_iexp(i S[𝐫_j])|^2
    
    
    S[𝐫_j]=   -𝐩_j0·𝐫_j0+I_p t_j0
       +∫_t_j0^∞d t
    (1/2𝐫̇_j(t)^2+2V(𝐫_j(t))),

with tunnelling timest_j0andp_j0⊥,𝐫_jdistributed via√(W(t_j0,p_j0,⊥))in order to compute the trajectories𝐫_j(t). In this work <cit.> a careful analysis of the action in different representations (e.g. position space and phase-space) was presented, which leads to the inclusion of the𝐩̇·𝐫term that was neglected in other approaches such as the TCSFA and QMTC. This is what leads to the factor of two in front of the potential, which also is used in the case of the CQSFA <cit.>. The SCTS also calculates the asymptotic momentum by modelling the electron as a Kepler hyperbola once the laser field has been turned off att_f. Using this the action fromt_fto∞may be calculated analytically, enabling a more efficient and accurate computation.



In Refs. <cit.> a direct comparison of the SCTS, QTMC and ab-initio methods is performed, see Fig. <ref>. In the low-energy region it was found the STMC has an improved angular distribution, that is closer to the TDSE results than the QTMC. This was traced back to the improved initial conditions given by the SFA.  The SCTS, like the QTMC, can be applied to understand holographic interference patterns. 
The SCTS was also extended to model ionization ofH_2, and again discrepancies with regard to the QTMC were identified in the low energy region <cit.>. With this extension, holographic imaging of small molecules using the SCTS is not far off. 



§ THE COULOMB QUANTUM-ORBIT STRONG-FIELD APPROXIMATION

The CQSFA combines the power of the path integral description of the photoelectron with the semi-classical approximation. This enables the use of quantum trajectories that put the Coulomb and laser fields on the same footing, with neither being treated perturbatively. The CQSFA bears some resemblance to the TCSFA method and contains the same type of trajectories but there are a few key differences: 1) The CQSFA solves the inverse problem enabling a description with much fewer trajectories. 2) The TCSFA neglects the phase space action termṗ·𝐫, included in the CQSFA. 3) The CQSFA considers perturbations in the action around the the saddle point solutions to second order, which leads to a stability prefactor term not included in any of the other models listed above. This term leads to important features that are visible in the final probability distributions, such as the “holographic” sidelobes <cit.>. A new model that also uses a similar path integral approach, which is still currently under development <cit.>, referred to as the semi-classical approximation, perforns an in-depth derivation of this prefactor.

In order to derive the Coulomb-quantum orbit strong-field approximation, we insert the closure relation∫d
𝐩̃_0 |𝐩̃_0⟩⟨𝐩̃_0 | =1in Eq. (<ref>). This yields

    M(𝐩_f)=
       -i lim_t→∞∫_-∞^t d t'
    ∫ d 𝐩̃_0 ⟨𝐩̃_f(t)
    |U(t,t') |𝐩̃_0⟩⟨𝐩̃_0 | H_I(t')| ψ
    _0(t')⟩  ,

where|𝐩̃_f(t)⟩=|ψ_𝐩(t)
⟩. The variables𝐩̃_0=𝐩_0+𝐀(t')and𝐩̃_f(t)=𝐩_f+𝐀(t)give the electron's initial and final field-dressed momenta, respectively.  The matrix element ⟨𝐩̃_f(t)
|U(t,t') |𝐩̃_0⟩is then calculated using time-slicing techniques  <cit.>. This leads to the CQSFA transition amplitude

    M(𝐩_f)   =-ilim_t→∞∫_-∞^tdt' ∫ d𝐩̃_0
    ∫_𝐩̃_0^𝐩̃_f(t)𝒟'
    𝐩̃∫𝒟𝐫/(2π)^3
       ×  
    e^i S(𝐩̃,𝐫,t,t')⟨𝐩̃_0 | H_I(t')| ψ _0  ⟩  ,

where𝒟'𝐩and𝒟𝐫are the integration measures for the path integrals <cit.>, and the prime indicates a restriction. The action in Eq. (<ref>) is given by

    S(𝐩̃,𝐫,t,t')=I_pt'-∫^t_t'[
    𝐩̇(τ)·𝐫(τ)
    +H(𝐫(τ),𝐩(τ),τ]dτ,

and the Hamiltonian reads as

    H(𝐫(τ),𝐩(τ),τ)=1/2[𝐩(τ)+𝐀(τ)]^2
    +V(𝐫(τ)),

where both the intermediate momentum𝐩and coordinate𝐫have been parameterized in terms of the timeτand the tilde indicates laser dressing, i.e.,𝐩̃=𝐩(τ)+𝐀(τ). 


Eq. (<ref>) is then solved using saddle-point methods. This means that we must seek variables𝐩,𝐫andt'such that the action is stationary, i.e., its partial derivatives with regard to these variables vanish. This gives

    [𝐩(t')+𝐀(t')]^2/2+V(𝐫(t'))+𝐩̇(t')·𝐫(t')=-I_p,

    ṗ=-∇_rV(𝐫(τ))

and

    ṙ= 𝐩+𝐀(τ).

Eq. (<ref>) gives the tunneling probability att', and Eqs. (<ref>) and (<ref>) give the subsequent electron propagation. These equations are quite different from their SFA counterparts in the sense that, according to (<ref>) the field-dressed momentum is no longer conserved due to the presence of the Coulomb potential. 

Within the saddle-point approximation, the ATI transition amplitude reads

    M(𝐩_f)∝-i lim_t→∞∑_s{[  ∂𝐩_s(t)/∂𝐫_s(t_s)] }^-1/2𝒞(t_s) e^i
    	S(𝐩̃_s,r_s,t,t_s) ,
wheret_s,𝐩_sand𝐫_sare determined by Eqs. (<ref>)-(<ref>), the term in brackets is associated with the stability of the orbit, and𝒞(t_s)is given by

    𝒞(t_s)=√(2 π i/∂^2	S(𝐩̃_s,r_s,t,t_s) / ∂ t^2_s)⟨𝐩+𝐀(t_s)|H_I(t_s)|Ψ_0⟩.

In practice, we use the stability factor∂𝐩_s(t)/∂𝐩_s(t_s)instead of that stated in Eq. (<ref>), which may be obtained employing a Legendre transformation. As long as the electron starts from the origin, this choice will not affect the action. We normalize  Eq. (<ref>) so that the SFA transition amplitude is obtained in the limit of vanishing binding potential, and, unless otherwise stated, take the electron to be initially in a1sstate <cit.>. We will refer to the product of the stability factor with Eq. (<ref>) as “the prefactor”.

One should note that, although the initial momentum states|𝐩̃_0⟩form a complete basis, by inserting the above closure relation we pay a price. While the CQSFA is very good for reproducing the dynamics of continuum and possibly high-lying bound states, it is not accurate for describing processes involving deeply bound excited states, such as excitation. For that purpose, a closure relation∫d
ϕ̃_c |ϕ̃_c⟩⟨ϕ̃_c|+∑_n |ψ_n⟩⟨ψ_n|=1involving the exact continuum states and the system's bound states would be more appropriate. In this case, however, how to apply the time slicing in order to obtain path integrals would be much less straightforward. 
Furthermore, the presence of the Coulomb potential, together with the fact that Eq. (<ref>) has complex solutions, brings a great deal of difficulty. First, since the initial conditions are taken at the origin, a singularity is expected atV(𝐫(t')). Second, strictly speaking, Eqs. (<ref>)-(<ref>) should be solved in the complex plane. This means that such equations themselves will contain poles or branch cuts. In order to solve these equations, we have performed a series of approximations, which have been improved from publication to publication. They will be briefly outlined below. For details, we refer to the original papers <cit.>.


  §.§.§ Contours and explicit expressions

Unless otherwise stated, we compute the ATI transition amplitude along a two-pronged contour, whose first and second parts are parallel to the imaginary and real-time axis, respectively. This is the most widespread contour used in the implementation of Coulomb-corrected approaches <cit.>, and it is very convenient for subsequent approximations. The contour starts at the complex timet'=t'_r+it'_i. Its first arm, fromt'tot'_r, may be related to the sub-barrier dynamics, and the second arm of the contour, from the real ionization timet'_rto the final timet, can be associated with continuum propagation.  The total action is then given by 
    S(𝐩̃,𝐫,t,t')=S^tun(𝐩̃,𝐫,t'_r,t')+S^prop(𝐩̃,𝐫,t,t_r').


In the first arm of the contour, the electron momentum is assumed to be constant and is kept as𝐩_0.  Physically, this means that the acceleration caused by the Coulomb potential is neglected in the under-the-barrier dynamics.
The action in the first arm of the contour then reads

    S^tun(𝐩̃,𝐫,t'_r,t')   =I_p(it'_i)-1/2∫_t'^t'_r[ 𝐩_0+𝐀 (τ)]^2dτ
       -∫_t'^t'_rV(𝐫_0(τ))dτ,
where𝐫_0is defined by 

    𝐫_0(τ)=∫_t'^τ(𝐩_0+𝐀(τ'))dτ',
 
which will be referred to as “the tunnel trajectory”. This type of trajectory has also been defined in many publications, in which the EVA <cit.> or the CCSFA <cit.> have been used.  
For the Coulomb potential, this approximation renders Eq. (<ref>) formally identical to its SFA counterpart (<ref>). This is highly advantageous as in some cases, such as monochromatic linearly polarized fields, this can be solved analytically.
In practice, however, the CQSFA solution for the tunnel ionization equation will be different from that obtained with the SFA. It must be matched to the saddle-point equations (<ref>) and (<ref>) giving the continuum propagation in the presence of the Coulomb potential. Thus, we must use the initial rather the final momentum as the independent variable. The actionS^tun(𝐩̃,𝐫,t'_r,t')inside the barrier is then calculated from the origin until the tunnel exit𝐫_0(t'_r). 

A widely used approximation is to take a real tunnel exit 

    z_0=Re[r_0z(t'_r)].

Although this will considerably simplify the computations by rendering the trajectories in the continuum real, this makes the problem contour dependent.  The main issue is that, while the actual problem should be contour independent and satisfy the Cauchy Riemann conditions, with this approximation analyticity is lost. Furthermore, from a physical point of view, this  assumption cannot be justified, as the only  variable that must be real is the final momentum at the detector. 
Still, as we will discuss next this will lead to an overall good agreement with ab-initio computations and avoids the issue of branching points and branch cuts.
Under  assumption (<ref>), the imaginary parts of the variables that will determine the single-orbit ionization probabilities, and probability densities, will be “picked up” in the sub-barrier part of the contour. The subsequent propagation will mainly lead to phase shifts in the real part of the action, which will influence the holographic patterns. Inclusion of the imaginary part of the tunnel exit will lead to complex trajectories, which may cause further suppression or enhancement in the ATI transition amplitudes during the continuum propagation. Complex trajectories may also lead to effective, further deceleration or acceleration, with regard to their real counterparts. This latter effect has been highlighted in <cit.>, and is related to the shifting probabilities of the trajectories in the continuum.

Another issue that must be addressed is the fact that there will be a logarithmic divergency at the endpointt'of the contour, as in this case the tunnel trajectory (<ref>) vanishes. We have applied several strategies in order to deal with this issue, such as introducing a very small imaginary time, close tot'and computing the limit of this extra term tending to zero <cit.>, or using the regularization method outlined in <cit.>, in which the contribution to the action coming from the Coulomb potential is matched to the asymptotic value of the bound-state Coulomb wave function. Further restrictions to the contour will be dictated by the condition𝐫^2=0, which is not necessarily at the origin. 
For more detailed discussions we refer to our recent publication <cit.>.

In the second part of the contour, the action reads

    S^prop(𝐩̃,𝐫,t,t'_r)   =I_p(t'_r)-1/2∫_t'_r^t[ 𝐩(τ)+𝐀 (τ)]^2dτ
        - 2∫_t'_r^tV(𝐫(τ))dτ,

where the factor 2 in front of the Coulomb integral is obtained as, for the specific, Coulomb-type potential used in our publications,𝐫 ·𝐩̇=V(𝐫(τ)). This extra factor is essential in order to obtain a good agreement with ab-initio methods, such as the right number of fringes in the near-threshold fan-shaped structure <cit.>. 
It is important to bear in mind that concepts such as “the tunnel exit” have raised considerable debate.  Indeed, if analyticity is preserved, the problem should be contour independent. This would render the tunnel exit arbitrary and hence of no physical significance  <cit.>. 
It is convenient to write  the action (<ref>) in terms of𝐩(τ)=𝒫(τ)+𝐩_f, such that𝒫→0whenτ→t.
This leads to 

    S^prop(𝐩̃,𝐫,t,t'_r)   =I_p(t'_r)-1/2∫_t'_r^t[ 𝐩_f+𝐀 (τ)]^2dτ
       +S_𝒫(𝒫,t,t'_r) +S^prop_C(𝐫,t,t'_r),

where 

    S_𝒫(𝒫,t,t'_r)=-1/2∫_t'_r^t𝒫(τ) ·[ 𝒫(τ) +2 𝐩_f +2 𝐀(τ)
    ]

is the phase difference due to the acceleration caused by the Coulomb potential, and 

    S^prop_C(𝐫,t,t'_r)=∫_t'_r^t2C/√(𝐫(τ)·𝐫(τ))dτ

is the Coulomb phase acquired during the electron's propagation in the continuum. These two phases, together with the sub-barrier Coulomb phase 

    S^tun_C(𝐫,t'_r,t')=∫_t'^t'_rC/√(𝐫_0(τ)·𝐫_0(τ))dτ,

vanish if the Coulomb coupling is zero, i.e., in the limitC→0. ThenV(𝐫)→0,𝐩_f→𝐩_0→𝐩and𝒫→0, so that the SFA action (<ref>) is recovered.



  §.§.§ Analytic approximations

For a monochromatic, linearly polarized driving field, Eqs. (<ref>) and (<ref>) describing the action in the sub-barrier and continuum parts of the contour can be computed analytically using a series of approximations. Thereby, a key issue is how to calculate the phases (<ref>), (<ref>) and (<ref>). This was the main topic of our previous publication <cit.>, to which we refer for details. Below we briefly sketch the main points to be considered. 

  1)Sub-barrier dynamics. In order to compute the sub-barrier action analytically, we apply the long wavelength approximation to the imaginary part of the ionization time to first order. This leads to convenient approximations for the tunnel trajectory, and to an analytic solution for the sub-barrier Coulomb phase (<ref>). It is worth stressing that this is an excellent approximation to the full CQSFA solution. This is justified by the fact that, in the full CQSFA, the momentum along the tunnel trajectory is already approximated by a constant value 𝐩_0. Thus, the only difference between both models is the long wavelength approximation applied to Im[t'] and to hyperbolic functions that depend on it.  The resulting transition amplitudes may then be factorized into a Coulomb- and an SFA-like part, which are very useful for interpretational purposes <cit.>. 
	
  2)Continuum propagation.  To find approximations that render the action S^prop_C(𝐫,t,t'_r) analytically solvable is more involved, as one cannot assume the intermediate momentum is fixed. It is however possible to find analytic expressions that make physical sense by (i) using the long-wavelength approximation around relevant times (not necessarily the ionization time); (ii) finding simple approximate functions for the intermediate momentum when dealing with the acceleration integral (<ref>); (iii) assuming that the field-dressed momentum is constant or piecewise constant when computing the Coulomb phase (<ref>). Typically, we take the difference 𝒫(τ) between a generic momentum 𝐩_j at a  relevant time and the final momenta 𝐩_f to be exponentially decaying and to vanish at the final time t. The relevant times for the long-wavelength approximation in (i) and the number of subintervals chosen in (iii) will depend on the orbit considered and will carry different physical meanings. One may construct a single or more subintervals using  the time of ionization, recollision, etc. If the CQSFA orbit is relatively close to its Coulomb free, SFA counterpart, it suffices to assume a single interval from t'_r to t for the continuum propagation, use the final momentum 𝐩_f to compute the Coulomb phase S^prop_C(𝐫,t,t'_r) and impose that, in this time interval, the momentum will tend monotonically from 𝐩_0 to 𝐩_f. 
	
  3)Incorporating soft recollision. If, for a particular orbit, there is pronounced interaction with the core, apart from the approximations employed in 1) and 2), one must incorporate one act of rescattering in order to obtain a good approximation for the full CQSFA solution. Specifically, the continuum propagation will require two subintervals: (i) from the ionization time t'_r to the time t_c of closest approach to the core, for which the electron coordinate parallel to the core vanishes, and (ii) from t_c to the final time t. We also assume that, in the first subinterval, the electron momentum perpendicular to the laser field polarization will remain the same, while the parallel momentum can be computed using 
	
    ∫_t'_r^t_c(p_oc∥+A(τ))dτ +z_o0=0,

	where p_oc∥ is the parallel component of the momentum at the instant of recollision, and z_o0 is the tunnel exit for a generic orbit o that interacts strongly with the core. 
	Upon recollision, the energy of the electron is conserved, the momentum changes instantaneously 𝐩_oc to 𝐩_o  and the rescattering angle does not change until the end of the pulse. These assumptions were “borrowed” from rescattered ATI, with the difference that, in our analytic approach, rescattering does not occur at the origin, but at a distance r_c perpendicular to the polarization axis. In the last subinterval, we assume the same monotonic behavior as for the previous orbits for the phase S_𝒫_o(𝒫_𝐨,t,t'_r) associated with the acceleration caused by the Coulomb potential in the continuum. 

An advantage of these analytical expressions is that they allow a detailed assessment of how the phases (<ref>), (<ref>), and (<ref>) influence the continuum dynamics and the holographic patterns. A drawback of the present analytic approach is that it is not a stand-alone model and requires the times of closest approacht_cto be obtained from the full CQSFA. 


  §.§.§ Treatment of branch cuts and singularities

The ideal scenario and most rigorous way to proceed would be to relax the assumptions in the previous section and solve the full complex problem. One should note, however, that the Coulomb phases (<ref>) and (<ref>) 
present in CQSFA the transition amplitude contain a square root of the complex variableu=𝐫(τ)·𝐫(τ). This means that, if saddle-point methods are used,  there will be branch cuts and two Riemann sheets for nonvanishing perpendicular final momenta. If the standard convention is applied, the branch cuts will lie in the negative real half axis of𝐫(τ)·𝐫(τ), i.e., for

    Re( 𝐫(τ)·𝐫(τ) )<0     and    Im( 𝐫(τ)·𝐫(τ) )=0.

Thereby, an important issue is that only one of the two Riemann sheets leads to physically meaningful results. Hence, it is vital to define a contour so that branch cuts are not crossed and the physical Riemann sheet is employed. This implies that the branch cuts must be mapped into the complex time domain, asτis the variable with regard to which𝐫·𝐫is parametrized. 
This is a non-trivial task if the residual Coulomb potential is incorporated in the electron's (complex) equations of motion for𝐫(τ)and𝐩(τ), as they will have branch cuts themselves. A first approximation is to employ Coulomb-free trajectories. Physically, this approximation works if 
the Coulomb acceleration is not substantial, such as for driving fields of large ellipticity, but is inaccurate for linearly polarized fields <cit.>. Still, it illustrates the branch cuts' overall behavior even in the latter case.  


A temporal mapping using Coulomb-free trajectories for linearly polarized fields shows that the branch cuts occur in pairs and are separated by gaps <cit.>. 
In <cit.>, a prescription for building contours that make use of this structure has been proposed. First, one calculates the middle point of each of the above-mentioned gaps. Subsequently, for a specific momentum value, a contour that passes through the middle point of each such gates is sought. The branching point pairs are solutions of𝐫(τ)·𝐫(τ)=0, while the condition𝐫(τ)·𝐯(τ)=0give the times of closest approach, which are located at the center of the gate, as well as the classical turning points. 

This method has been applied in <cit.> to above-threshold ionization in the context of soft recollisions under the Coulomb barrier. Therein, the trajectories were taken to be Coulomb-free.  It was shown to improve the agreement of the SFA and ab-initio computations in orders of magnitude, especially near the direct ATI cutoff energy of2U_p. This was a counter-intuitive  and rather intriguing result related to  under-the-barrier collisions. 
Nonetheless, its implementation poses a series of difficulties. First, one must find the midpoints in advance, which may not always be possible. For instance, this could prove particularly difficult for Coulomb-distorted trajectories.
Second, to calculate photoelectron momentum distributions using this method, it would be necessary to construct and apply a different contour for each pair of final electron-momentum components.

Therefore, we apply a different method: we keep the standard two-pronged contour discussed in Sec. <ref>, whose first and second arms are parallel to the imaginary and the real time axis, respectively. After computing the trajectory, we test whether it crosses a branch cut. If it does, we deform the contour in such a way that it goes around the branch cut.  If a specific timet_balong the second arm of the contour satisfies the branch-cut conditions (<ref>), it may be employed to find the branching pointt_kdefined  by

    𝐫(t_k)·𝐫(t_k) = 0+0 i.


The intersectiont_bof the branch cut with the real axis and the complex timet_kassociated with the branching point can then be used to build a three-pronged contour along which the Coulomb phase is computed. Explicitly,

    S^(cut)_C=∫_c_1V[𝐫(τ)]d τ+∫_c_2V[𝐫(τ)]d τ+∫_c_3V[𝐫(τ)]d τ,

wherec_1goes fromt_btot_kfollowing the upper side of the branch cut,c_2returns fromt_ktot_balong the other side of the cut andc_3connects them together with a semi-circle aroundt_k. This approach is applicable if the transverse momentum component is non-vanishing. For vanishing transverse momentum the branching points will merge into a pole and the integrals alongc_1andc_2become divergent.  An example of this contour is presented in Fig. <ref>. Physically, a branch cut may be associated with a collision (see a detailed discussion in <cit.> and a brief explanation in Sec. <ref>). 

In order to make the problem tractable for Coulomb-distorted complex trajectories, additional approximations must be made. First, we assume that their imaginary parts behave as their Coulomb-free counterparts during the continuum propagation. This means that they will be kept constant when the time contour is chosen to be along the real time axis. This approximation may be loosely justified as follows: (i) In the sub-barrier part of the contour, we neglect the Coulomb potential for the CQSFA, thus describing the dynamics in the same way as in the standard, Coulomb-free SFA; (ii) The final momentum𝐩_fmust be real, as it is an observable, which means that thelim_t→∞ Im[𝐫(t)] =const.If a branch cut is met, we move vertically in the complex time plane by assuming that the influence of the slowly varying Coulomb field on the electron intermediate momentum can be neglected. This is consistent with the physical picture for which collisions happen in much shorter time scales than the driving-field cycle. The approximations described in (ii) simplify the problem by  rendering all momenta real in the second part of the contour.

One should note that the treatment of branch cuts and singularities in Coulomb-distorted strong-field approaches is a contentious issue. For instance in <cit.> it has been suggested that, when the electron is close to the core, such features are absent. Therein, an effective potential has been derived with no singularities or branch cuts, and it has been argued that the classical action employed in most Coulomb-distorted approaches is a limit that is only valid for negligible quantum spreading of the electronic wave packet. A softened effective potential is also obtained in the coupled coherent states method, in which the full Hamiltonian is expanded into a trajectory-guided, coupled coherent state basis <cit.>. Rigorous ways of treating branch cuts, caustics and singularities are more established in other areas of research (see, e.g., <cit.> for in depth discussions).



§ ORBITS, TRANSITION AMPLITUDES AND HOLOGRAPHIC PATTERNS

In this section, we highlight the main results of our previous publications, but provide a more unified, broader discussion. Such results have been obtained with linearly polarized monochromatic fields of frequencyωand amplitudeE_0=2ω√(U_p), whereU_pis the ponderomotive energy. The corresponding vector potential reads

    𝐀(t)=2√(U_p)cosω t ê_z.

This simplifying assumption has several advantages. First, the ionization times can be written analytically as functions of the momentum components parallel and perpendicular to the driving-field polarization. Second, the periodicity of the field facilitates the interpretation of many holographic patterns and can be used in analytical estimates. 
For this particular case, the CQSFA action obtained combining the sub-barrier and continuum contributions reads


    S(𝐩̃, 𝐫,t,t')
        =   (I_p + U_p ) t'
    +1/2𝐩^2_f t'_r+i/2𝐩^2_0  t'_i
    +U_p/2ωsin(2ω t')
    
             
    + 2√(U_p)/ω[ p_ 0∥sin(ω t')-(p_0∥-p_f∥)sin(ω t'_r) ] 
    -∫^t'_r_t' V(𝐫_0(τ))dτ
          
    -1/2∫_t'_r^t𝒫(τ)· (𝒫(τ)+2𝐩_f+2𝐀(τ))dτ
    -2 ∫^t_t'_rV(𝐫(τ))dτ,

wherep_ν∥andp_ν⊥, withν=0,f, give the momentum components parallel and perpendicular to the laser-field polarization. The explicit expression for the direct ATI transition amplitude within the SFA is given by

    S_d(𝐩,t')   =(p^2_∥+p^2_⊥/2+I_p+U_p)t'
       +2 p_∥√(U_p)/ωsin[ω t']+U_p/2 ωsin[2ω t'],

where the electron momentum remains constant throughout (i.e.,𝐩_0=𝐩_f=𝐩). 
A direct inspection of the SFA transition amplitude (<ref>) shows that their dominant term is proportional tot'. The remaining terms are either bounded, or increase more slowly. This implies that the rough behavior of the action is expected to mirror that of the solutions of the saddle-point equations. 
For the CQSFA transition amplitude (<ref>), this will depend on the orbit in question. If the orbits remain in regions for which the laser field is dominant, i.e., if it is a direct orbit or if the Coulomb potential merely deflects the electron, then the dominant term will be proportional tot'. If the electron crosses a region in which the residual binding potential becomes dominant, then the Coulomb phase may prevail and further investigation is necessary. 



 §.§ Trajectory classification

The solutions of the CQSFA saddle-point equations and of their SFA counterparts can be associated with electron trajectories, whose quantum interference leads to the  ultrafast holographic patterns. For the SFA transition amplitude(<ref>) leading to direct ATI, the only integration variable is the ionization timet'. Within the steepest descent method,t'can be obtained by solving Eq. (<ref>). Specifically for linearly polarized monochromatic fields,  this equation has an analytical solution. Thus,  one may parametrize the ionization times in terms of the electron momentum components parallel and perpendicular to the laser-field polarization. 
The general SFA solution related to an eventewithin a cyclecis 

    t'_ec = 2π n/ω±1/ωarccos(-p_∥∓ i√(2 I_p+p^2_⊥)/2√(U_p)),

wherenis an integer.  The times predicted by Eq. (<ref>) occur in conjugate pairs, and only those in the upper half complex plane, i.e., fulfilling the conditionIm[t'_ec]>0, have physical significance. 

Within a single cycle, the SFA exhibits two solutions: the electron may either be released in the direction of the detector, or half a cycle later in the opposite direction. In the latter case, the electron will be freed in the “wrong” direction, brought back to its parent ion and then reach the detector. These solutions are known as the short orbit and the long orbit, respectively, and have been briefly mentioned in the context of the temporal double slit (Sec. <ref>).
Specifically, we use the solutions

    t'_1c   = 2π/ω-1/ωarccos(-p_∥+ i√(2 I_p+p^2_⊥)/2√(U_p),) 
    
    t'_2c   = 1/ωarccos(-p_∥- i√(2 I_p+p^2_⊥)/2√(U_p))
 
for each cycle. Solutiont'_1c(t'_2c) corresponds to the short (long) orbit forp_∥>0, with the scenario reversing forp_∥<0. Their imaginary parts will be the same, and their real parts will differ. Fig. <ref> shows the behavior of the real parts of the two SFA solutions and of the corresponding actions, as functions of the parallel momentum component.  Therein, analytic expressions for the SFA case can also be found.  In general, the expression for the coordinate𝐫obtained from the SFA is complex. This will also be the case in the CQSFA as the tunnel exit𝐫(t'_r)is not real if no further approximations are used. 


Once the Coulomb potential has been included, the number and the topology of the orbits will change. The central Coulomb force will act upon the transverse momentum component of the electron and thus considerably alter the dynamics.  We will first analyze these changes for real trajectories by neglecting the imaginary part of the tunnel exit. This is the approach employed in most Coulomb-distorted strong-field approaches, and significantly simplifies the computations. Subsequently,  we will discuss the complex-trajectory case.


  §.§.§ Real trajectories

The following classification has been first proposed in <cit.> and has been successfully employed by us. It makes use of the initial and final momentum components, and of the direction of the tunnel exit with regard to the detector.  It may be summarized as follows:

  * Type 1 orbits: the tunnel exit is located on the same side as the detector, i.e., z_0p_f∥>0, and the transverse momentum component does not change sign during the continuum propagation, i.e., p_0⊥p_f⊥>0. They can be related to the short SFA orbit in direct ATI. The residual binding potential decelerates an electron along this orbit.
	
  * Type 2 orbits are related to those events in which the electron is ionized from the opposite direction to the detector,  i.e., z_0p_f∥<0 but eventually changes its direction. It is deflected by the combined action of the field and the residual binding potential. Its Coulomb-free counterpart would be the long SFA orbit in direct ATI. Nonetheless, the presence of the Coulomb potential distorts this orbit and leads to field-dressed Kepler hyperbolae. The net effect of the Coulomb potential for an electron along this orbit is to accelerate it <cit.>.
	
  * Type 3 orbits are similar to type 2 orbits with regard to the direction of the tunnel exit (z_0p_f∥<0), and, in typical scenarios, its hyperbolic shape. The key difference is that, for type 3 orbits, the interaction with the core is stronger. This means that (i) the momentum component perpendicular to the laser-field polarization will change sign, i.e., p_0⊥p_f⊥<0, and (ii) this type of orbit has no counterpart in the direct SFA. Our recent work has in fact shown that rescattering is an important assumption with regard to type 3 orbits, if one is willing to recover the spider-like structures <cit.>.
	
  * Type 4 orbits start in the same direction as the detector, but, unlike orbit 1, they go around the core. This means that z_0p_f∥>0, but  p_0⊥p_f⊥<0, and that these orbits are very strongly influenced by the Coulomb potential. These orbits can be loosely associated with backscattered trajectories in high-order ATI, although, strictly speaking, the rescattered SFA does not allow for events encircling the core.

For clarity, we provide an illustration of these orbits (Fig. <ref>), together with a summary of their classification criteria (Table <ref>). In the language employed in photoelectron holography, one may say that orbits type 2 and 3 are forward-scattered trajectories, and this scattering may in principle range from a light deflection to a hard collision. Type 4 orbits are backscattered trajectories. On a more rigorous note, the two laser-dressed hyperbolae in Fig. <ref> actually form a torus in the three-dimensional space for ionization events parallel to the driving-field polarization, which becomes tilted for non-vanishing scattering angle. For the SFA, this torus decreases to a point, so that orbits type 2 and 3 become degenerate. For details see our previous publication <cit.>. 


The orbit classification presented above invites the following questions: 

  Q1:Can one systematically relate the CQSFA orbits to those obtained with the SFA, be it direct or rescattered? The main difficulty here lies on the fact that, while the SFA is constructed as a field-dressed Born series, the CQSFA is not. It is difficult to establish when and at what point in configuration space the electron “returns” or  “rescatters” if the binding potential and the driving field are considered on equal footing. In contrast, for the SFA the electron may only return to the origin, and its very structure ensures that either no, or a single act of rescattering occurs. 
	
  Q2:Does the classification provided suffice to describe the dynamics? This is an important question as it could be that the conditions provided in Table <ref> are satisfied for orbits of different topology.


In order to address the question Q1, it was necessary to devise ad-hoc criteria for when the electron recollides and whether this collision is “soft” or "hard”, based on (i) the Bohr radiusr_B; (ii) the perimeterr_Tdefined by the (real) tunnel exit; (iii) the distance of closest approachr_cof the electron subsequently to rescattering. We then assume that a hard collision will happen if an electron orbit enters the region defined by the Bohr radius. In this region, the Coulomb forces will be dominant. A soft collision will occur ifr_cis larger than the Bohr radius, but smaller thanr_T. Finally, if the electron's orbit does not cross the region delimited byr_T, it may be related to the ATI direct orbits, as the laser field will be dominant. A schematic representation of this region is provided in Fig. <ref>. A boundary-based approach is also employed in the ARM to separate hard collisions from deflection.


The real parts of the ionization times as functions of the final parallel momentum and scattering angle are good  indicators of how the CQSFA orbits behave. For instance, Fig. <ref>(a) shows that the real partst'_rof the CQSFA ionization times for orbits type 1, 2 and 3 follow the direct SFA solutions from below, but, depending on the orbit, behave differently in the high-energy limit, i.e., when|p_f∥|is large. For orbit 1,t'_rtends to the SFA solution, because increasing the photoelectron energy renders the influence of the Coulomb potential less and less relevant. In contrast, the real partst'_rassociated with the CQSFA orbits 2 and 3 never tend to the direct SFA solutions. For orbit 4, this behavior is particularly extreme and never follows that of the SFA. Instead, the ionization times are practically constant. The marked discrepancies for orbits 3 and 4 are a consequence of the fact that these types of orbit interact much more with the core. This means that to relate them to a direct ATI orbit is inaccurate for a broad parameter range. Fig. <ref>(b) also shows that the imaginary parts of such orbits do not approach those predicted by the SFA.  This issue will be discussed in more detail in Sec. <ref>.

Other good indicators are the time and distance of closest approach, denotedt_candr_cin our previous publications <cit.>. These quantities are plotted in Fig. <ref>, as functions of the final scattering angle [panels a) and c)] and of the absolute value of the final momentum|𝐩_f|[panel (b)].  The figure shows that the distances of closest approach determined by the direct and rescattered SFA are limits to what is obtained from  the CQSFA. This is particularly clear if one looks at the behavior with regard to the final scattering angle, which, within the CQSFA framework, should encompass from direct orbits to hard scattering.  If one considers an orbit with small scattering angles, the distance of closest approach is the tunnel exit and follows what can be obtained if the Coulomb coupling is taken to zero from above, i.e., what one would expect from the direct SFA. Once the angleθ_f=1.1πis reached, the distance of closest approach drops drastically, crossing the region delimited by the Bohr radius and tending towards the origin. That would be precisely the expected return coordinate for high-order ATI. An even clearer pattern can be observed for the CQSFA times of closest approach, which follow the solution of the saddle-point equation determined for direct and rescattered ATI very closely.


With regard to question Q2, i.e., of the classification  discussed in this work being sufficient, Fig. <ref> provides a very good example of topologically very different trajectories, which however fall under the umbrella of type 3 and type 4 orbits, according to the classification in Table <ref>. The standard case for such orbits, shown in Fig. <ref> and described above, is only represented in panel a). Panel b) shows multi-pass orbits, which are driven past the core many times before recollision, and panel c) displays orbits that undergo hard scattering  before the laser field changes sign. These orbits, and in particular their influence on holographic patterns have not yet been investigated in detail.  There is however a strong indication that multipass orbits are responsible for near-threshold effects, and patterns such as the “inner spider”, both from simplified methods <cit.> and from the CQSFA <cit.>. A proper study of these orbits will however require novel asymptotic expansions, as there are radical changes in the number and nature of the saddles <cit.>. 

It is important to stress that this is not the only classification employed in the literature. For instance, in <cit.>, orbits leading to holographic patterns have been classified as direct (D_1), and forward scattered (R_1andR_2). TrajectoriesR_1have large, initial positive transverse momenta, while for trajectoriesR_2the initial transverse momenta are small and negative. Although we cannot directly state, as the direction of the final transverse momenta are not given, it seems thatR_1andR_2may be associated with type 2 and type 3 orbits, respectively.  This impression is supported by the stronger influence of the core and higher ionization rate that exists for theR_2trajectories in <cit.> (see also supplemental material therein).



  §.§.§ Complex trajectories and branch cuts


If complex trajectories are taken, two legitimate questions are: (i) Will any of the trajectories discussed above, either Coulomb-free or Coulomb-distorted, meet a branch cut? (ii) If so, what is the physical meaning of this encounter? 
In order to understand the problem and ultimately answer these questions, we will map the branch cuts in the complex momentum plane and match this behavior with that of qualitatively different trajectories.

In Fig. <ref>, such a mapping is presented for Coulomb-free trajectories, i.e., the SFA orbits 1 and 2. The figure shows two sets of branch cuts separated by gaps that may get closer or distance themselves from the real axis depending on the parameters used. For orbit 1, the branch cuts do not cross the real time axis in the time range of interest, i.e, forRe[τ]>t'_r. However, they do for orbit 2. This implies that the contour must be altered for this latter orbit.  

If the Coulomb potential is incorporated, as shown in Fig. <ref>, we also see that orbit 1 does not encounter a branch cut, but the other three types of orbits do depending on the choice of the electron momentum components. In Fig. <ref>(a), orbit 2 crosses the branch cut, while for the other types of orbits, branch cuts are avoided. In Fig. <ref>(b), the final momentum has been carefully chosen as to lead to a final momentum almost perpendicular to the driving-field polarization direction. In this case, the branch cuts intersect the real time axis for orbits 3 and 4. The situation is particularly complex in  Fig. <ref>(c), for which the branch cuts associated with orbit 4 cuts through the real time axis twice within a field cycle. A closer inspection [Fig. <ref>(d)] even shows that, in this case, the branch cuts overlap, which would rule out the applicability of the present procedure. 
Physically, a trajectory encountering a branch cut is strongly associated with an act of return or rescattering. For that reason, this may occur for all orbits except orbit 1.  The fact that orbit 4 may cross twice a branch cut within a field cycle reflects its behavior of going around the core.  This physical meaning has already been mentioned in <cit.>, and in <cit.>, which state that, every time a contour pass through the gates, soft recollisions do occur. Because the variables are complex, such recollisions influence not only the phases, but also the amplitudes associated with specific orbits. One should note that, due to the restrictions associated with the method in <cit.>, in such studies only Coulomb-free trajectories were taken into consideration. In contrast, our method allow us to treat Coulomb-distorted trajectories. 


 §.§ Sub-barrier dynamics and single-orbit distributions

In the standard formulation of the CQSFA, understanding the contributions of each orbit to the overall probability density means having a closer look at the sub-barrier dynamics and the imaginary parts of the ionization times. In the CQSFA framework, the ionization probability is proportional toexp[-2Im[S(𝐩̃, 𝐫,t,t')]], where the variables upon which the action depends are determined using the saddle-point equations in the first part of the integration contour discussed in Sec. <ref>.


Fig. <ref> presents the single-orbit distributions computed for orbits type 1, 2 and 3, without and with prefactors. For the contributions of orbits 1 and 2 (left and middle panels), one may easily identify two peaks along the parallel polarization axis, forp_f∥ ≠0and elongated shapes near thep_f∥axis.  This behavior mirrors that of the imaginary parts of the CQSFA action, which are plotted in the upper panels of Fig. <ref>, and that observed forIm[t'], plotted in Fig. <ref>(b). The figure clearly shows two minima in both cases, whose positions correspond to the maxima in this distributions. 
Qualitatively, this doubly-peaked shape is very different from that observed for the SFA. In the SFA case,Im[S_d(𝐩,t')]andIm[t']exhibit a single minimum for vanishing momentum, which leads to a single peak in the corresponding electron-momentum distribution. This may be understood as a result of the fact that the barrier to be overcome by the electron is narrowest at peak field. For peak-field times, the electron tunnels with vanishing momentum. Such events, however, are highly suppressed in the CQSFA due to the residual Coulomb attraction. This may be easily understood as an electron with vanishing momentum at the tunnel exit will be pulled back by the core. 
For orbit 3, one observes a much broader momentum range delimited by caustics. This is consistent with the fact that the imaginary part of the ionization time is very flat and much closer to zero for this orbit. Once the prefactor is included, however, the probability density associated with this orbit is restricted to a narrow range along the polarization axis. For orbit 4 the significant contributions are located perpendicular to the polarization axis, so that they only start to play a role for nearly right scattering angles. It is noteworthy that, for orbits 3 and 4, the shapes of the electron momentum distributions are strongly influenced by the prefactors, which also causes a drastic reduction of the yield. This is a strong indicator that these orbits are far less stable than orbits 1 and 2. 

Hence, we conclude that, for orbits 1 and 2, the key influence on these distributions comes from the Coulomb potential modifying the SFA-like part of the sub-barrier actionS^tun(𝐩̃,𝐫,t'_r,t')via the ionization times, while for orbits 3 and 4 the prevalent contribution is less clear. Thus, it is legitimate to ask what role the sub-barrier Coulomb phase plays. Analytical estimates for each term can be used to disentangle both contributions. The imaginary parts of the action with and without the Coulomb integral are shown in Fig. <ref>(b), while the corresponding single-orbits distributions for orbit 2 are shown in the lower panels of the figure. The main effect of the Coulomb integral is to broaden the distributions in the direction perpendicular to the polarization axis, and lead to secondary maxima. 

More detailed information is given  in Fig. <ref>, which shows the analytic estimates together with the full CQSFA electron momentum distributions for a wide parameter range, computed using orbits 1 and 2. Throughout, the SFA-like terms lead to elongated shapes along the polarization axis, while the contributions from the sub-barrier Coulomb phase have non-trivial behavior. In most cases, the latter are located at or around thep_f⊥axis. This is consistent with the momentum regions occupied by the probability distributions associated with orbits 1 and 2, along the polarization axis, and for orbit 4, in the direction perpendicular to it (see Fig. <ref>). For orbits 1 and 2, the coupling with the laser field dominates, while for orbit 4 the binding potential prevails. 
For very long wavelengths, the influence ofS^tun_C(𝐫,t'_r,t')is minimal. This may be attributed to the long electron excursion amplitudes for the electron in this regime, which reduce the influence of the core.  As the wavelength decreases, the sub-barrier Coulomb phase becomes more important by modifying the shape of the electron momentum distributions. 

Our studies also show that sidelobes are already present for single-orbit distributions, so that they are primarily caused by sub-barrier contributions. In the seminal work <cit.>, they were attributed to quantum-interference effects. Quantum interference does enhance the sidelobes, but it is only a contributory cause. This assertion was only made possible because the CQSFA allows one to separate and visualize the contributions of specific types of orbits in the presence of the binding potential and assess the influence of sub-barrier dynamics in each. Interestingly, features that could be associated with sidelobes have also been reproduced by a classical model and associated with rescattering in ATI <cit.>.

Finally, it is noteworthy that many semi-classical approaches that include the residual Coulomb potential weight the electron's release in the continuum with the quasi-static ADK tunnelling rate, such as the QTMC <cit.>, and the SCTS  <cit.> methods. It is a well-known fact, however, that these rates are maximal at peak field. However, the results presented in this section do not seem to support this picture. 
Still, both the QTMC and the SCTS exhibit a good agreement with ab-initio  methods and reproduce many of the holographic features encountered in experiments. This apparent contradiction may be understood if one considers that an electron freed with vanishing momentum would be recaptured by the core into a high-lying excited state. This mechanism is known as frustrated tunnel ionization (FTI) <cit.>, and has been recently incorporated in the SFA in <cit.>. For early studies of transient recapture of an electron in a high-lying state see <cit.>.  This mechanism is explicitly mentioned in the context of the QTMC <cit.>, and of an in-depth study using classical-trajectory methods <cit.>. It is important to note the resemblance between our single-orbit distributions for orbits 1 and 2 and those in <cit.>, with two  peaks stretched along the direction of the laser-field polarization. Frustrated tunneling ionization has also been studied in conjunction with the attoclock <cit.>, with the very low energy structure <cit.> and with the fan-shaped holographic fringes <cit.>. Nonetheless, the issue is not fully settled, as recent studies indicate that non-adiabatic ionization conditions play an important role in enhancements that occur for the ATI signal near photoelectron energy2U_p<cit.>. For a thorough discussion of several tunneling criteria see <cit.>, a phase-space analysis of tunneling using initial value representations see <cit.>, and non-adiabatic effects in circularly polarized fields see <cit.>.



 §.§ Holographic patterns: well-known and overlooked

The objective of this section is to discuss the key types of interfering trajectories contributing to specific patterns encountered in ATI photoelectron momentum distributions.  This will include not only the fan and the spider-like structures, but also patterns not yet reported in experiments. 


  §.§.§ Preamble: several types of temporal double slits.

Thereby, a good place to start is the direct SFA. An example of several possible types of inter- and intracycle interference  within the SFA is provided in Fig. <ref> (left and right panels, respectively). Intercycle interference leads to the ATI rings, which  become finer if more and more cycles are added, while intracycle interference leads to double-slit patterns.  A striking feature is that, even if only two solutions per cycle exist, several intra-cycle patterns can be constructed. This follows from the variety of phase differences that may be obtained depending on the time difference between both solutions is chosen. One may for instance choose the interfering orbits so thatRe[t'_2c-t'_1c]is smaller (panel (d)) or greater (panel (e)) than half a cycle, respectively, or if no restriction is imposed (panel (f)). We denote the interference stemming from the conditionsRe[t'_2c-t'_1c]<π/ωas type A interference, while that coming fromπ/ω≤Re[t'_2c-t'_1c]≤2π/ωis called type B interference. The former restriction has been widely used in the literature, while the latter case has only been addressed systematically in our previous publication <cit.>. A methodical  study of the type A double-slit interference is presented in <cit.>. The present  classification is not related to that used in the seminal paper <cit.>. 
Type A interference, when symmetrized with regard top_∥=0,leads to the double slit pattern reported in  <cit.>. Fig. <ref>(d) shows such fringes, which are nearly vertical and slightly divergent at the origin. Type B interference, shown in Fig. <ref>(e), leads to much finer, slightly convergent fringes. The differences in the widths may be understood by inspecting the differences in the real parts of the ionization times and actions of the two contributing orbits (see Fig. <ref>). The further apart they are, the finer the fringes they will be. If no restriction is imposed uponRe[t'_2c-t'_1c], the fringes move from finer to thicker in the example presented. In Fig. <ref>, types A and B intracycle interference are marked by circles and rectangles, respectively. Analytic interference conditions for all double-slit cases were derived in our previous publication <cit.>.



  §.§.§ The fan and the spider.





We will now view the fan-shaped pattern discussed in Sec. <ref> in the light of the CQSFA.  In Fig. <ref>, we plot photoelectron momentum spectra computed in the polarization plane, using only orbits 1 and 2 over five cycles. Apart from the ATI rings caused by inter-cycle interference, the figure shows the fan-shaped structure, both in the CQSFA and ab-initio results. This is in striking disagreement with the SFA results, which exhibit the nearly vertical fringes that are characteristic of the type A interference. One should note as well that there is excellent agreement between the number of fringes in the fan for the CQSFA and the TDSE. This is not the case in many Coulomb corrected approaches, such as the CCSFA and CVA (see Sec. <ref>).
If only type A intra-cycle interference is allowed, the fan becomes even clearer (see Fig. <ref>).  A  comparison with the SFA results, displayed on the right panels of Figs. <ref> and <ref>, supports the statement that the fan is due to unequal angular distortions in the temporal double-slit pattern caused by the Coulomb potential. The influence of the binding potential decreases with photoelectron energy and scattering angle, being very pronounced near the polarization axis, but practically cancelling out close toθ=90^∘. Thus, resonances with excited bound states are not necessary to reproduce the fan.



In Fig. <ref>, we plot the spider-like pattern calculated with the CQSFA.  In agreement with previous studies in the literature <cit.>, our results show that its main contributors are trajectories 2 and 3, which are forward-scattered orbits starting in the same half cycle.  However, the CQSFA allows a much deeper analysis as we can use it to probe directly how the spider forms. An interesting feature is that there is a caustic around the spider, which is very visible if the prefactor is not included. If the prefactor is inserted, the fringes become much more localized near thep_f∥axis. One should note that there are no type A or B interferences for the spider. 

One may also use the analytic approximations mentioned in Sec. <ref> to dissect what each of the phasesS_𝒫(𝒫,t,t'_r)andS^prop_C(𝐫,t,t'_r)given by Eqs. (<ref>) and (<ref>), as well as soft recollision, does to the spider-like fringes. This can be seen in Fig. <ref>, in which we also include the outcome of the full CQSFA and that of its analytic approximation including all phases [Fig. <ref>(a) and (b), respectively]. The sub-barrier Coulomb phase in taken into consideration in all cases.
If soft recollision is neglected [Fig. <ref>(c)], the central part of the spider and its fringes are poorly reproduced. This shows that it is an essential part of the physics leading to the spider. The role of the Coulomb phase is to modify the slope of the spider-like fringes. If it is removed, they are “bent down” beyond the polarization axis [Fig. <ref>(d)], but are still present as they are caused by rescattering and the phaseS_𝒫(𝒫,t,t'_r)related to the acceleration caused by the Coulomb potential. A remarkable example of this is provided in Fig. <ref>(e), for which the acceleration phase has been removed. This leads to the spider-like fringes being truncated at the direct-ATI cutoff energy2U_pprescribed by the SFA. If rescattering is removed [Fig. <ref>(f)], the electron momentum distribution bears a striking resemblance to those obtained for direct ATI using the SFA if only the long orbits are taken. This is expected as the ionization times for orbits 2 and 3 are relatively close. 

It is important to mention that there exist  interpretations for the spider that are alternatives to the assumption that it is caused by the interference of two types of trajectories. For instance, in <cit.>, spider-like fringes are obtained by considering the interference of a spherical wave with an incoming plane wave, and, in <cit.>, the spider is attributed to an axial caustic singularity and to glory effects in rescattering. This is justified by computing the transverse width of the electron momentum distributions and the positions of the fringes. Nonetheless, the issue is not yet settled as the glory trajectories in <cit.> resemble orbits 2 and 3, which form a torus in the three-dimensional space (see discussion in Sec. <ref> and <cit.>). Furthermore, the trajectory-based methods used in the comparisons therein were less sophisticated than many of those reported here as it used either the plain SFA or trajectories in which the binding potential has been added perturbatively. Still, it raises important questions with regard to the role of caustic singularities along the polarization direction.



  §.§.§ Overlooked holographic patterns.
  
Apart from the fan and the spider, there are many overlooked patterns that may be obtained if (i) type 4 orbits are taken into consideration; (ii) if one considers ionization times displaced by more than one half cycle and different types of orbits. Below we provide a few examples, taken from <cit.>. In Fig. <ref>, we show spiral-like patterns that are obtained if orbit 4 is incorporated. In experiments and ab-initio computations, these patterns can be easily obfuscated and/or mistaken for ATI rings.  They are most visible close to the perpendicular momentum axis and relatively high photoelectron energies, because in this region the spider-like fringes and the fan-shaped structure are strongly suppressed. Caustics associated with this orbit and with orbit 3 are also visible in the figure.


By restricting the real parts of the ionization times to specific intervals, one may obtain a myriad of patterns. One such pattern is the fan-shaped structure given in Fig. <ref> and <ref>, which results from differences smaller than half a cycle and symmetrization with regard top_f∥=0[see also Fig. <ref>(a)]. Another possibility is to take this difference to lie between half a cycle and a whole cycle (type B interference according to the notation in <cit.>). Examples of these patterns are provided in Fig. <ref>.  Depending on the choice of interfering orbits, the interference can be divergent [Figs. <ref>(a) and (f)], convergent [Figs. <ref>(b), (c) and (d)], spiral-like [Figs. <ref>(e) and (g)] or borderline cases. Clearly, if no symmetrization or restriction is employed, we will find quite different fringes on the left and right sides ofp_f∥=0. 


One may also shift[t']in a whole cycle for one of the interfering orbits in a pair, and this will lead to strikingly different patterns. However, because the types of orbits in the interfering pairs are different, the patterns are very distinct from ATI rings. 
In Fig. <ref>, we display results obtained using orbits starting from the same side of the core, i.e., 2 and 3 (upper panels), and 1 and 4 (lower panels). For forward scattered trajectories, apart from the spider [<ref>(a)], if one of the orbits are shifted of a full cycle, the fringes will become much finer, and either convergent [Fig. <ref>(b)] or spiral-like [Fig. <ref>(c)]. Finer fringes indicate that the phase shift between orbits 2 and 3 has been increased. 
If one considers the interference of backscattered and direct trajectories (lower panels of Fig. <ref>), it is clear that, if orbit 1 starts one cycle later [Fig. <ref>(e)], the fringes become much broader and there is a fan-shaped structure. This is due to the fact that it will be more “in phase” with a backscattered orbit that started much earlier but went around the core. Clearly, if orbit 4 starts later, or within the same cycle as orbit 1, the fringes will be much finer [Fig. <ref>(d) and (f)]. 

The message to be taken out of this section is that the interplay between the Coulomb potential and the driving field leads to a far more complex scenario than anticipated in the seminal paper <cit.>.  Those pictures have been drawn by assuming that the tunnel exit remained fixed, the amplitudes associated with each orbit were equal and that rescattering took place at one single point: the origin. The CQSFA made it possible for us to blur many of these distinctions. First, both the weight of the orbit and the tunnel exit are dictated by the sub-barrier dynamics and the first part of our contour, as discussed in Sec. <ref>. This means that the tunnel exits will vary with regard to the field and that the probability amplitudes associated with each interfering orbits will be unequal, which, in some instances, will lead to the blurring of the fringes. Second, the phase shifts acquired in the continuum will be determined by the interplay between the laser field and the long-range binding potential, which means that they will differ a lot from those predicted in <cit.>. Finally, it is not trivial to define whether an electron has returned to the core, or suffered rescattering, and what type of rescattering it underwent. Indeed, to make this kind of statement we followed two strategies: (i) to employ analytic approximations in which individual contributions to the overall action were switched on and off at will; (ii) to compare the CQSFA orbits with those in the direct or rescattered SFA. Thus, some of the patterns predicted in <cit.> that have been sought in experiments are an over simplification and may need to be reassessed.



 §.§ Comparison with ab-initio methods

Apart from allowing the in-depth analysis presented above, the CQSFA exhibits a good agreement with the ab-initio solution of the time-dependent Schrödinger equation over a wide momentum range.  This is shown explicitly in the following three figures, which also provide information about what types of interference would be prominent in a realistic setting, and how this is influenced by the field parameters. The CQSFA electron momentum distributions were computed without any restriction upon the real parts of the ionization times and using no symmetrization. 

Fig. <ref> displays the results from the CQSFA taking into consideration the contributions of orbits 1, 2 and 3. Throughout, both the near-threshold fan and the spider-like pattern are very visible, and become more and more symmetric as the pulse length increases. For the shortest pulses, one may see traces of the converging fringes that result from the interference of orbits 1 and 3. Once the number of cycles increases, these fringes become washed out and there is the presence of very clear ATI rings.
The effect of the prefactor is mainly to suppress contributions of high-energy photoelectrons for large scattering angles, improving the resemblance of the CQSFA with the TDSE computation. 

It is noteworthy that the number and approximate location of the spider-like and fan-shaped fringes agree both in the CQSFA and TDSE. This agreement is particularly good for moderate and near-threshold photoelectron energy. The slope of the spider-like structure is however slightly steeper for the ab-initio computation. In the CQSFA, this slope is caused by the interplay of the phasesS_𝒫(𝒫,t,t'_r)andS^prop_C(𝐫,t,t'_r), which undergo some approximations as the continuum trajectories are taken to be real. Furthermore, all CQSFA results reveal a caustic for high-energy photoelectrons, which marks a boundary within which the present asymptotic expansions are valid for orbit 3. This caustic is also exacerbated by the continuum trajectories being real. Finally, one may see that, close to the perpendicular momentum axis, the ATI rings appear to be “broken” for the TDSE distributions, while they remain intact for the CQSFA. This caustic has also been reported in <cit.>.

Fig. <ref> shows that this effect is caused by the contributions of orbit 4. If this orbit is included in the computations, the spiral-like structures resulting from its interference with orbit 3 will cause the ATI rings to break close to thep_f⊥axis. In <cit.>, the broken rings have been attributed to the constructive and destructive interference of the direct SFA short and long orbits. However, in the energy range for which they occur, one expects rescattering to play an important role. Therein, discrepancies with regard to the experimental results have also been reported. This means that the problem is still open to a great extent.  Orbit 4 will also lead to an additional caustic in the distribution [see Figs. <ref>(b) and (e)]. This caustic has been identified in <cit.> using the adiabatic theory, and has been rightly attributed to backscattered trajectories. In <cit.>, this shape has been obtained from backscattered SFA orbits with maximum photoelectron momentum, and was explained analytically using the mapping (<ref>).

Finally, in Fig. <ref>, we show that, depending on the driving field parameters, the agreement between the CQSFA and the TDSE may improve or worsen. For instance, in the limit of long wavelengths (see upper panels in the figure),  the slope of the spider obtained in the TDSE computation decreases, so that it becomes more similar to the outcome of the CQSFA. Overall, in this wavelength regime the similarity between the result of both approaches is striking. This is possibly due to the long electron excursion amplitudes involved, which reduce the influence of the core.  For the very same reason, the agreement worsens with decreasing wavelength (see middle and lower panels in the figure). The figure also shows that the analytic approximations employed in order to analyze the patterns in greater depth work well together and bear a strong resemblance to the results  obtained with the full CQSFA (see far left and left panels in the figure). 



 §.§ Branch-cut corrections

In the following, we incorporate complex trajectories in the CQSFA and analyze the effect of branch cuts in the photoelectron momentum distributions. The simplest scenario is given in Fig. <ref>, in which the Coulomb phase has been incorporated in the electron propagation but Coulomb-free, SFA orbits were used. The ovals in Fig. <ref> and in the lower panels of the figure delimit the region for which orbit 2 does not cross a branch cut. In the absence of branch-cut corrections, fringe discontinuities and a trumpet-shaped structure near thep_f∥axis closely following the oval boundary are very visible [Figs. <ref>(a) and (e)]. These features are eliminated once the branch-cut corrections have been applied [Figs. <ref>(c) and (f)]. Overall, we also see that the distributions computed with the strong-field approximation [Fig. <ref>(d)], or with real trajectories [Fig. <ref>(e)] are much broader and decay more slowly. This supports the statement that complex trajectories may lead to an effective deceleration of the electronic wave packet <cit.>. One should bear in mind, however, that their imaginary parts may also lead to enhancements in the photoelectron yield, as reported in <cit.>. Throughout, the Coulomb phase considerably modifies the shape of the quantum-interference fringes, although the holographic structures discussed in the previous section are absent. The figure also makes it clear that our procedure does not work well on thep_f∥axis, as previously stated [see Fig. <ref>(c)].
 

If Coulomb-distorted orbits are used, similar features are observed. An example is provided in Fig. <ref>, in which specific holographic patterns are shown. For the interference between orbit 1 and 2, one clearly sees that the distributions computed with complex orbits [Fig. <ref>(b) and (c)] decay much faster than that in which real orbits are taken into the continuum propagation [Fig. <ref>(a)]. Furthermore, orbit 2 crossing a branch cut also leads to discontinuities, which our method eliminates. Similar features can also be observed for the spider, and the spiral (middle and lower rows in the figure, respectively). Furthermore, employing complex trajectories alters the slope of the spider and changes the spacing and the contrast of the spiral. Finally, the fringes near thep_f⊥axis that exist for the spider if real trajectories are taken are softened in the complex-trajectory case. This is due to a faster decay in the contributions of orbit 2.

In Fig. <ref>, we show a direct comparison of all the versions of the CQSFA developed and employed by us with the ab-initio solution of the time-dependent Schrödinger equation. If branch-cut corrections are taken into consideration, complex trajectories make the electron momentum distribution decay faster for increasing transverse momentum componentsp_f⊥. This means that they become more localized near the polarization axis, improving the agreement with the ab-initio solution. Furthermore, the slope of the spider-like fringes is nearly horizontal if real trajectories are taken, but becomes closer to the TDSE outcome in the complex-trajectory case.
There is however a drawback in the present procedure for correcting branch cuts, which tends to worsen the overall agreement between the CQSFA and ab-initio computations. The approximations employed in the Coulomb-distorted case to render the problem tractable (see Sec. <ref>) work well for orbit 1, reasonably for orbit 2, but less so for orbits 3 and 4. This is expected, as those two latter orbits interact very strongly with the core, and are therefore quite different from their direct SFA counterparts (see upper panels of Fig. <ref>). An immediate consequence is that the contributions of orbits 3 and 4 to the overall probability densities are over-enhanced. This is caused by the imaginary parts of orbits 3 and 4 being  overestimated. In order to counteract this problem, we reduce the contributions of orbits 3 and 4 by employing an empirical factor. This considerably improves the agreement, as  a direct comparison of Figs. <ref>(d) and (f) shows. 




§ TRENDS


Overall, the existing trends for photoelectron holography reflect those one may observe for ultrafast imaging as a whole. There is a shift from merely structural questions towards the aim of steering electron and core dynamics in real time, which requires more sophisticated modeling, more complex targets and tailored fields, which are chosen in such a way as to probe specific properties of these targets. Thereby, an important issue is to move from an unstructured continuum, which can be approximated by field-dressed plane waves, towards a structured continuum.  This, together with an appropriate treatment of multielectron dynamics, is a key requirement for photoelectron holography in polyatomic molecules, solids and nanostructures. A schematic representation is presented in Fig. <ref>.


The simplest possible example of modeling a structured continuum, which is to include the residual Coulomb potential and compute photoelectron spectra for a one-electron atom, has been extensively discussed in this review, with particular emphasis on the CQSFA. Our results show that the interplay between the residual potential and the laser field plays a huge role. The very fact that the Coulomb potential alters the ionization dynamics and the number and topology of the interfering orbits means that the resulting holographic patterns are strikingly different from those predicted with the strong-field approximation and in the seminal papers <cit.>. Clearly, the present studies have only paved the way for the understanding and control of this interplay in more complex systems. 
Below we elaborate on some of the  present trends:

  * Molecules and multielectron atoms.  These are the quintessential ultrafast imaging targets, and have called a great deal of attention since the mid 2000s. Apart from the geometry of the bound states involved, that may not be spherically symmetric, removing one electron will alter the dynamics of the the core. This removal may for instance create a hole that will migrate and/or induce further rearrangement processes in the remaining charge, such as relaxation, excitation and electron migration. The external field may also polarize the core electrons, or even induce double and multiple ionization. In case of molecules, there may be also the coupling of electronic and nuclear degrees of freedom, such as vibration and rotation, or even dissociation and fragmentation. The latter effects will be particularly important for large systems, or highly excited states. To be able to probe such systems using photoelectrons, one must ensure that a strong interaction with the core takes place. This can be achieved by focusing on the near-threshold energy regime, for which the electron orbits are substantially distorted by the core, or on patterns involving orbits that undergo appreciable acceleration by the residual ion.
	Well-known examples of the former and the latter situation are the fan and the spider, respectively. Enhancements in the fan have for instance been related to resonances with Rydberg states in dimers <cit.>, nuclear-electronic coupling of degrees of freedom <cit.>, and suppression of the fan-shaped fringes has been attributed to frustrated tunnel ionization <cit.>.
	Phase shifts related to molecular alignment <cit.>, as well as imprints of nodal planes  and evidence of coupling with nuclear degrees of freedom  have been identified in the spider <cit.>. Recent computations indicate that polarization effects influence both the fan and the spider <cit.>.  Backscattered trajectories and the resulting quantum interference may also be employed. For instance, in <cit.> backscattered trajectories have been used to retrieve information about multiple orbitals in LIED of an aligned organic molecule, and required modelling beyond the single-active orbital and plane rescattering wave approximation. However, special care must be taken that they are not obfuscated by more prominent patterns. An example is the fishbone structure in  <cit.>, which could only be identified by using differential holograms in which the spider was removed. This fishbone pattern is highly dependent on the internuclear separation, and thus could be used for tracking nuclear dynamics. 
	
	
  * Tailored fields, such as linearly <cit.> or orthogonally polarized <cit.> two-color fields, few-cycle pulses <cit.>, elliptically or circularly polarized fields <cit.>, or bicircular fields <cit.>, allow a greater deal of control in the electron ionization, propagation and deflection/recollision by the potential. Hence, they may be used to probe a myriad of properties such as tunneling times <cit.>, the influence of the Coulomb potential <cit.> and holographic patterns that are not easily resolved <cit.>. In principle, obfuscated holographic patterns could be revealed by (a) moving more prominent holographic patterns to other momentum regions or suppressing them altogether; (b) changing the curvature and shapes of specific holographic fringes. 
	For instance, in <cit.> it has been proposed, within the SFA framework, that orthogonal two-color (OTC) fields may be employed to distinguish the fishbone holographic patterns related to backscattered orbits from those stemming from the temporal double slit. The SFA predictions, however, are very distinct from the outcome of ab-initio computations in which the residual potential is taken into consideration. Still, OTC fields may move overlapping holographic patterns into distinct momentum regions.
	Furthermore, orthogonally polarized fields  provide a good testing ground for the influence of the Coulomb potential on photoelectron momentum distributions. In <cit.>, asymmetries in photoelectron momentum distributions for neon in orthogonally polarized two-color fields could be traced back to the influence of the Coulomb potential using classical-trajectory Monte Carlo (CMTC) methods. Thereby, a strong deviation from the mapping 𝐩_0=-𝐀(t')	has been encountered, which has been related to the Coulomb potential holding back departing electrons, i.e., to frustrated tunnel ionization.	
	Bicircular fields are of paramount interest as they may be used for probing chiral molecules, and also allow for rescattering. Above-threshold ionization in such fields is expected to reveal a myriad of LES and interference patterns <cit.>, some of which exist only due to nonvanishing initial momenta and the presence of the Coulomb potential <cit.>. 
	
  * Plasmonically enhanced fields. It has been suggested that plasmon-induced resonant enhancements may amplify modest fields in several orders of magnitudes <cit.>. This amplification does however lead to driving fields with an inhomogeneous spatial profile. Effective models in which this inhomogeneity is approximated within a one-electron framework show that this causes significant distortions in electron-momentum distributions <cit.>. A quantum-orbit analysis using the SFA also shows that the photoelectron orbits significantly change in this case <cit.>. It is noteworthy that even very simplified models as those used in both papers may already reveal a very rich, non-trivial dynamics. For instance, in <cit.>, in the context of HHG, we have found that even very small inhomogeneities lead to spatial confinement, bifurcations, symmetry breaking, two very distinct time scales and different stability regimes. This all plays a huge role in the continuum electron dynamics, and is expected to influence above-threshold ionization as well.
	
  * Longer wavelengths. Typically, the complex systems mentioned above have lower ionization potentials. Furthermore, in specific cases, such as in the seminal work of <cit.>, it is advantageous to prepare the target in an excited metastable state. 
	This means that, to probe such systems, it is desirable (i) to reduce the driving-field intensity and keep the ponderomotive energy fixed; (ii) to ensure tunneling is still the dominant ionization mechanism. Condition (i) implies that high-energy photoelectrons may be created, as the ATI cutoffs  both for direct and rescattered electrons are functions of U_p, and that one is still operating below the saturation intensity. Condition (ii) means that all the formalism and the physical interpretations developed for the quasi-static regime are valid, and that ionization is very much dependent on the instantaneous field strength. Both conditions can be achieved by using fields of longer wavelengths, which have become increasingly popular since the beginning of decade <cit.>. For that reason, many of the studies reported in this review have been performed in the mid-IR regime. Thereby, another difficulty arises: if the
	magnetic field associated with the laser field causes the electron motion perpendicular to the
	driving-field polarization to become greater than a Bohr radius, then the dipole approximation breaks down <cit.>. This will alter several features described in this work, such as Coulomb focusing and defocusing <cit.>, and the topology of the electron orbits <cit.>. For studies of radiation pressure in the context of LIED see, e.g.,  <cit.>, and for non-dipole effects with elliptically polarized fields see <cit.>.

The above-stated list is obviously non-exhaustive, and ultrafast photoelectron holography will bring with itself many challenges and questions, both theoretical and experimental. With the present review, we hope to have provided a thorough and multifaceted contribution to the understanding of photoelectron holography's past, present and future.



§ FUNDING

This research was is part funded by the UK Engineering and Physical Sciences Research Council (EPSRC) (grants EP/J019143/1 and EP/P510270/1). The latter grant is within the remit of the  InQuBATE Skills Hub for Quantum Systems Engineering.



§ ACKNOWLEDGEMENTS

We would like to thank  X. Liu, X. Y. Lai, T. Das, A. Al-Jawahiry, and  S. Popruzhenko for their collaboration, which led to some of the publications discussed in this review. We are also grateful to D. B. Milošević, D. Bauer, M. Yu. Ivanov, W. Becker, M. Lewenstein, E. Pisanty, L. Torlina, J. Z. Kaminski, N. I. Shvetsov-Shilovski, W. Koch and D. J. Tannor for useful discussions.



§ ABBREVIATIONS

The following abbreviations are used in this manuscript:


	ADK     Ammosov-Delone Krainov

	ATI     Above-threshold ionization

	ARM     Analytical R-Matrix theory

	CCSFA    Coulomb-corrected strong-field approximation

	CMTC    Classical-trajectory Monte Carlo method

	CQSFA     Coulomb Quantum-orbit Strong-Field
    approximation

	CVA     Coulomb-Volkov approximation

	DATI     Direct above-threshold ionization

	EVA    Eikonal Volkov approximation

	HATI     High-order above-threshold ionization

	HHG     High-order harmonic generation

	HOMO     Highest-occupied molecular orbital 

	ISFA     Improved strong field approximation

	KFR     Keldysh-Faisal-Reiss

	LIED    Laser-induced electron diffraction

	LES     Low-energy structure

	LFA     Low-frequency approximation

	NSDI     Nonsequential double ionization

	OTC     Orthogonally polarized two-colour

	QRS     Quantitative Rescattering Theory

	QMTC     Quantum-trajectory Monte Carlo method

	SCTS     Semi-classical two-step model

	SFA     Strong-field approximation

	RABITT     Attosecond Burst By Interferences of
    Two-photons Transitions

	TCSFA     Trajectory-based Coulomb strong-field
    approximation

	TDSE     Time-dependent Schrödinger equation

	VLES     Very low energy structure

	WKB     Wentzel-Krammers-Brillouin

	ZES     Zero energy structure


§ REFERENCES
iopart-num