Strong-field ionization of complex molecules

Strong-field photoelectron momentum imaging of the prototypical biomolecule indole was disentangled in a combined experimental and computational approach. Experimentally, strong control over the molecules enabled the acquisition of photoelectron momentum distributions in the molecular frame for a well-defined, narrow range of incident intensities. A novel, highly efficient semiclassical simulation setup based on the adiabatic tunneling theory quantitatively reproduced these results. Jointly, experiment and computations revealed holographic structures in the asymptotic momentum distributions, which were found to sensitively depend on the alignment of the molecular frame. We identified the essential molecular properties that shape the photoelectron wavepacket in the first step of the ionization process and employ a quantum-chemically exact description of the cation during the subsequent continuum dynamics. The detailed modeling of the molecular ion, which accounts for its polarization by the laser-electric field, enables the simulation of laser-induced electron diffraction off large and complex molecules and provides full insight into the photoelectron's dynamics in terms of semiclassical trajectories. This provides the computational means to unravel strong-field diffractive imaging of biomolecular systems on femtosecond time scales.


I. INTRODUCTION
Strong-field ionization is a versatile and powerful tool for the imaging of molecular structure, allowing simultaneous access to valence-electronic topology and atomic positions with sub-atomic-unit resolution in time and space [1]. A generalization of the underlying imaging techniques, which were born in the field of atomic physics, to complex biomolecular targets promises deep insights into the biochemical machinery of life and this extension is the aim of the current work.
In strong-field-ionization experiments the desired structural information can be obtained in two complementary ways: through the observation of the photoelectron itself, which is discussed in detail below, or through the burst of light that is emitted when the electron recombines with the cation, referred to as high-harmonic generation. The latter was employed in many studies on atoms and molecules unraveling valence-orbital [2] and atomic structure [3], but is not discussed further in our work.
The current article focuses on the photoelectron wavepacket as it picks up molecular structural imprints during its formation and motion in the continuum. The entire strong-field process may descriptively be resolved into three stages, which are linked through the photoelectron's propagation in the combined time-and positiondependent field that is exerted by the laser and the cation: (i) The photoelectron is born in the continuum. Its wavefunction is shaped by the ionization potential [4] and the nodal geometry [5][6][7][8] of the initial bound electronic state.
(ii) Upon propagation within close vicinity to the ionized target the photoelectron wavepacket picks up additional information about the cation through interaction with its electric potential. For circularly polarized laser pulses this encounter is essentially limited to early times right after birth in the continuum and, if the target molecule is chiral, results in photoelectron circular dichroism [9][10][11]. For linear polarization, however, the photoelectron is, in first instance, quickly driven away from the molecule yet might partly return at a later time.
If the momentum at return is rather small, i. e., comparable to the initial momentum distribution at birth, the wavepacket may holographically interfere with another portion of itself that does not return to the cation. Those holographic interferences encode the tunnel exit position and the initial momentum distribution of the photoelectron [12][13][14]. However, if the re-collision occurs at the maximum of the kinetic energy, namely about three times the ponderomotive potential [15], the photoelectron may intrude deeply into the ionic field probing the instantaneous positions of the nuclei. This process includes only a temporally narrow slice of the initial wavepacket that returns after less than one optical cycle near a zero-crossing of the laser field. It is referred to as laser-induced electron diffraction [16][17][18][19][20].
(iii) After both, the decay of the laser field and the escape from the cation's Coulomb potential, are completed, the photoelectron wavepacket assumes its asymptotic form in momentum space. This is the stage one can readily access experimentally. It is highly desirable to extract unambiguous information about the target molecule, that was imprinted onto the photoelectron, from this detectable final momentum distribution. Although attempts of inversion recipes were formulated already, building for example upon the quantitative rescattering theory [21] for the evaluation of laser-induced electron diffraction experi-ments, they rely on severe approximations regarding the initial photoelectron wavepacket or its continuum motion. For instance, the returning electron's momentum distribution is frequently assumed to be free of any angular modulation. But since molecular systems may give rise to rather structured initial electron waves [7,[22][23][24], the target-dependent formation and continuum motion of the photoelectron needs to be well understood in order to come up with a general inversion technique.
During the last twenty years numerous models have been brought into being, which elegantly tackle various different aspects of strong-field physics. The numerical evaluation of the time-dependent Schrödinger equation (TDSE) [25,26] represents the most rigorous approach among them, but remains mostly restricted to atomic systems. Further approaches to tackle strongfield ionization of complex molecules rely on the timedependent density-functional theory [8,27]. They can be efficiently applied for wavelengths up to the mid-infrared regime, but so far do not allow to track the photoelectron's dynamics in the continuum.
A broad class of models building upon the strong-field approximation (SFA) in the single-active electron picture allows for the treatment of more complex targets. Those SFA-based theories commonly divide the whole process into two steps corresponding to stage (i) and stages (ii+iii) described above. In the low-frequency limit, the initial continuum wavefunction can be described by means of quasistatic tunneling theories [4], which can be step-wise supplemented to grasp non-adiabatic behavior during tunneling [28] and allow for the consideration of imprints from the initial bound electronic state [6,7]. An alternative, wavelength-independent description of the early photoelectron wavefunction can be obtained through the sub-barrier Coulomb-corrected SFA [29].
The very classical nature of the electron's motion in continuum has stimulated several studies, which rely on its decomposition into many classical trajectories. A careful analysis of these trajectories can then help to gain an illustrative understanding of intricate strong-field phenomena [30]. The two most important approaches to describe the electron's phase in the continuum are the frequently used Coulomb-corrected SFA [12,29] and semiclassical time-dependent propagators [31].
We note that in the vast majority of published models the cation is approximated by a point charge during the continuum propagation of the freed electron.
Here, we introduce a novel and very efficient semiclassical simulation setup. Our approach builds upon the adiabatic tunneling theory and allows for the consideration of structural imprints from the initially bound orbital. The continuum propagation of the photoelectron wavepacket is treated through semiclassical trajectories that experience the detailed, quantum-chemically exact, field of the cation.
Due to its efficiency our model is ideally equipped for disentangling the structures of even highly complex molecules through strong-field ionization. In conjunction FIG. 1. Indole molecules were three-dimensionally fixed in the laboratory frame according to two different alignments of their principal-polarizability frame, depicted on the left, with respect to the polarization axis of the ionizing laser field, depicted as double-headed red arrow. In the experiment, this was achieved through irradiation with a nonresonant laser field carrying elliptical polarization, sketched by the blue ellipses. Both laser beams propagated collinearly along the X axis. The molecule's HOMO is shown in an isosurface representation for all four degenerate orientations.
with a strong-control experimental scheme, which enabled a well-defined ionization study in the molecular frame, we present a combined experimental and computational survey on the strong-field ionization of the prototypical biomolecular building block indole (C 8 H 7 N), which is the major ultraviolet-absorbing chromophore of proteins.
Unless stated otherwise, atomic units (at. u.) are used throughout this article.

II. SEMICLASSICAL MODEL
Our description of the strong-field ionization process employs a two-step approach: Treating the photoelectron wavefunction at birth in the continuum through quasistatic tunneling theory allows for the analytical description of the first step, which is explained in section II A. Subsequently, the final wavefunction is found by sampling its dynamics in the combined electric field of the laser and the polarizable cation on plane-wave quantum orbits, which is depicted in section II B, while section II C illustrates the procedure to compose the probability density of the final wavefunction from the asymptotic trajectory properties. We tested our model by simulating the photoelectron momentum distribution from metastable xenon at a wavelength of 7 µm, which yielded excellent agreement with the published experimental and TDSE results [12].
Throughout this manuscript the reference coordinate frame will be set through the polarization of the laser field, i. e., laboratory coordinate frame, with the Z axis corresponding to the polarization axis. The molecular frame is assumed to be linked to the laboratory frame via the indole molecule's principal axes of polarizability, α xx > α yy > α zz . Two different molecular alignment configurations will be used, x-axis alignment with (z, ±y, ±x) → (X, Y, Z) and y-axis alignment with (z, ±x, ±y) → (X, Y, Z). Note, that the signs of x and y axes are not fixed, giving rise to four simultaneous orientations, while the sign of the z coordinate is insignificant as a result of the molecule's C s point group. See Fig. 1 for a scheme of all four orientations for both alignment cases in the laboratory coordinate frame.

A. Initial photoelectron wavefunction
The quasistatic Ammosov-Delone-Kraǐnov (ADK) tunneling theory [4], empirically extended to the barriersuppression regime [32], yields the instantaneous ionization rate w. In conjunction with the initial momentum distribution in the plane orthogonal to the polarization axis, the probability density of the initial photoelectron wavefunction [33] can be obtained: ξ = 2I p / is the ratio of the initial state's characteristic momentum, i. e., the square root of twice the ionization potential I p , and the instantaneous field strength . η represents the electronic-structure imprint of the initial bound state onto the momentum distribution of ψ 0 and can be acquired by linking the momentum in the initial state k to the momentum at the tunnel exit p via the temporal saddle point equation [33]: Here, A Z (t s ) is the Z component of the vector potential at the complex saddle-point time t s . Subsequently, the corresponding molecular electronic structure factor can be reorganized to the form η · e i∆φ0 , with η providing a momentum imprint onto ψ 0 and ∆φ 0 introducing an initial phase offset. A simplified description of the highest occupied molecular orbital (HOMO) as a p X orbital, results in a constant phase shift and a momentum imprint Within the scope of this model description it is sufficient to draw tuples of (t 0 , p 0X , p 0Y ) from |ψ 0 | 2 to get complete sets of initial phase-space coordinates, since the tunnel exit r 0 and p 0Z can be unambiguously inferred.
r 0 is deduced from the field-direction model [34] and the corresponding momentum component along the polarization axis, p 0Z , can be accessed with the aid of the adiabatic tunneling theory embodying the first nonadiabatic correction [28]. Within the barrier-suppression regime the immanent excess energy is transformed into additional longitudinal momentum [34].
Finally, starting from the instant of birth, the planewave phases along the quantum orbits, which are employed to trace the wavefunction's dynamics, are given by [31] The second term of the integrand in the action integral of (5), −2/r = V ( r) − r · ∇V ( r), originates from a pure far-field description of the potential energy, V = −1/r. In order to also properly track the phases of those quantum trajectories that intrude deeply into the cationic potential, one would need to model both, V ( r) and ∇V ( r), in quantum-chemical detail. However, due to the large de Broglie wavelength of the photoelectron wavepacket under the present experimental conditions, vide infra, we approximate V ( r) by this coarse description without losing agreement with the experimental data.

B. Continuum propagation
In order to acquire the wavefunction's asymptotic probability density in momentum space, the classical equations of motion are numerically solved within the total combined electric field of the laser and the cation. For the sake of numerical efficiency, the ionic field − ∇V (t, r) is computed at varying levels of accuracy, in three different spatial domains: At large distances from the center of nuclear charge a multipole description of the electric field, including monopole, permanent dipole moment, and induced dipole moment, is employed, in close proximity to the nuclei the total electric field is approximated by an electric monopole with an atom-dependent effective charge, fully neglecting the field of the laser, and elsewhere the quantum-chemically exact electric field is utilized. See section III for the spatial boundaries of these domains.
In the close-proximity regime the application of Kepler's laws of orbital mechanics [35] to the resulting effective two-body problem allows for an analytical description of the orbital motion and provides the means to sort out trajectories that end up on stationary orbits around the ion. If the total energy of such a two-body system is negative, recombination of electron and cation is assumed and the corresponding quantum orbit does not contribute to the asymptotic photoelectron wavefunction.
Starting from a quantum orbit's individual time of birth in the continuum, t 0 , the ordinary differential equations for position, momentum, and phase are numerically integrated until the time when the laser field is fully decayed. The asymptotic momentum and the relative phase can then be obtained through Kepler's laws [31].

C. Asymptotic wavefunction
From the asymptotic momenta and relative phases of all available quantum orbits the probability density of the final photoelectron wavefunction in momentum space, |ψ ∞ ( p)| 2 , is then composed as follows. First, the momentum-space volume of interest is divided into bins in spherical energy coordinates: Subsequently, for each bin U j , θ k , ϕ l the modulus of the coherent sum over all associated M jkl asymptotic plane waves is computed: Before generating the modulus of the bin-wise coherent sum, the complex-phase histogram may be coherently symmetrized to exploit the intrinsic symmetry of the photoelectron wavefunction; see section III for the detailed treatment of the three-dimensionally aligned indole molecule. After interpolation on a Cartesian momentum-space grid the asymptotic momentum map is then obtained through summation along the Y axis, which mimics the momentum projection in a velocity-map imaging (VMI) spectrometer. The use of spherical energy coordinates comes with three major advantages: (i) The curvatures of the bin edges resemble the shapes of the interference structures of |ψ ∞ | 2 . (ii) The volume element scales with p · sin(θ), giving rise to larger bins with increasing radial momentum, which partly makes up for the smaller number of trajectories. (iii) The inherent energy coordinate allows for a straightforward emulation of focal-volume averaging. While (i) and (ii) result in a highly increased convergence of |ψ ∞ | 2 as a function of the number of samples, (iii) enables the efficient handling of a distribution of incident intensities as it is often encountered in strong-field experiments. Within a narrow range of incident intensities and in absence of electronic resonances the change of |ψ ∞ | 2 with respect to intensity is dominated by the ponderomotive shift of the nonresonant above-threshold ionization (ATI) interferences [36]. Those ATI structures represent a purely radial pattern with an energy spacing equal to the photon energy, ω, between adjacent maxima. Thus, a small uncertainty in the ponderomotive energy due to focal-volume averaging can easily be emulated by introducing this uncertainty to the energy axis of |ψ ∞ | 2 , which smears out the ATI pattern but leaves the angular structures unchanged.

III. COMPUTATIONAL SETUP
Ionization of the indole molecule was described as occurring from a single p X orbital, which resembles HOMO and HOMO−1 both in energy and symmetry; see Table I. Already for HOMO−2, 2.7 eV below the HOMO, the corresponding ADK tunneling rate is roughly three orders of magnitude smaller. So contributions from all lower-lying orbitals were readily neglected.
The electric field of the ionizing laser was defined as with a cosine-square envelope that for n = 44 reproduced the full width at half maximum (FWHM) of the intensity envelope of the experimental pulses. A peak electric field of 0 = 2.9×10 −2 was utilized, corresponding to the larger of the two pulse energies employed in the experiment, vide infra.
The ionization potential of the molecule was computed considering Stark shifts up to second order: with the field-free ionization potential I between cationic and neutral species. These properties were computed using the GAMESS quantum chemistry software [37,38] at the MP2/aug-cc-pVTZ level of theory, see section A for details. Throughout the numerical treatment, the cation was frozen in the nuclear equilibrium geometry of the neutral molecule. Samples of initial phase-space coordinates were generated through rejection sampling [39] from the probability density (1). The tunnel-exit positions and possible additional momentum offsets along the polarization axis in the barrier-suppression regime were obtained based on quantum-chemistry calculations using Psi4 [40] at the HF/aug-cc-pVTZ level of theory, following the fielddirection model, see section B for details. For the laser intensity utilized, the ionization process was found to enter the barrier-suppression regime close to the peak electric field for the x-axis alignment case. In contrast, for y-axis alignment the tunnel exit tightly follows the simple description r 0 = − I p / 2 [41], which considers the laser field to be quasistatic and ignores the Coulomb distortion of the barrier.
Within a (60 a 0 ) 3 cube, centered at the molecule's center of nuclear charge, ab initio quantum chemistry calculations were used to describe the electric field of the cation, see section C. This detailed modeling of the ion's field is necessary, because the indole cation experiences strong polarization in response to the laser-electric field. Such behavior is expected to be general for molecules with π-conjugated electrons and, thus, makes this external-field dependent description indispensable for biomolecular systems. For both, the application of the field-direction model and the evaluation of the quasistatic field of the ion, quantum-chemical calculations were conducted beforehand for discrete positions and external fields, see section C for details. Further values were obtained through interpolation in the semiclassical propagation setup.
At distances < 0.5 a 0 with respect to any of the molecule's nuclei the total electric field was modeled using an electric monopole with an atom-dependent effective charge. For the H, C, and N atoms in indole effective charges of +1, +4 and +5 were used. These values were retrieved by finding the integer-charged monopoles that fitted the quantum-chemically exact field best on spherical surfaces with 0.5 a 0 radius around the atomic centers. Fig. 2 shows a slice through the total electric field that is exerted by the indole cation.
The equations of motion (6) were integrated following an embedded Runge-Kutta scheme with adaptive step size (Cash-Karp) [42]. For each alignment case roughly 5 × 10 9 non-recombining quantum orbits were traced. Due to the use of the spherical energy coordinate system introduced in section II C roughly 5 × 10 8 plane-wave samples sufficiently reconstructed all relevant radial and angular features of |ψ ∞ | 2 .
To account for the coherent superposition of all four degenerate orientations in ψ ∞ the symmetry of ψ ∞ with respect to inversion of the X and Y axes was exploited a posteriori by coherently symmetrizing the histogram of bin-wise coherent sums in (8). The X-axis symmetry is a result of the molecule's C s point group and the Y -axis inversion transforms two orientations into each other. The Z-axis symmetry of the aligned molecular ensemble could not be exploited in the same way, as the corresponding axis inversion would also flip the sign of the laser-electric field. As a consequence, for each of the two alignment scenarios two orientation cases were independently simulated and coherently added.

IV. EXPERIMENTAL SETUP
The experimental setup was described in detail elsewhere [8,43,44]. In brief, a cold molecular beam was created by supersonically expanding indole seeded in 95 bar of helium through an Even-Lavie valve [45] operated at a repetition rate of 100 Hz and heated to 110 • C. Using an electrostatic deflector an ultracold high-purity sample was obtained [46]. Employing elliptically polarized laser pulses at a peak intensity of 1 × 10 12 W/cm 2 , which carried a saw-tooth temporal shape with a rise time of 500 ps, the molecular ensembles were aligned within the center of a VMI spectrometer [47]. These alignment laser pulses, which were spectrally centered at 800 nm and had a polarization ellipticity of 3:1 in intensity, allowed for the nonresonant, quasi-adiabatic fixation [43,48] of the indole molecule's principal axes of polarizability in the laboratory frame. Two different alignment scenarios, Fig. 1, could be realized by rotating the polarization with a half-wave plate. The resulting degree of alignment was estimated to be cos 2 δ ≈ 0.9 [49].
A second laser pulse, spectrally centered around 1300 nm, with a duration of 70 fs (FWHM) and linear polarization (ellipticity 200:1), singly ionized the molecules and photoelectron momentum maps were recorded using a high-energy VMI spectrometer. A position-sensitive detector, consisting of a microchannel-plate stack, a phosphor screen, and a high-frame-rate camera, was used for counting and two-dimensional momentum-mapping of individual electrons. To lower the impact of focal-volume averaging onto the incident-intensity distribution, intensitydifference photoelectron-momentum maps [50] were obtained using peak intensities of 2.5 and 3.0 × 10 13 W/cm 2 . The resulting range of the ponderomotive potential ∆U p = 0.83 ω = 0.79 eV between these two peak intensities is still too large to resolve nonresonant ATI structures but represents a compromise between incident-intensity resolution and the amount of signal that is left in the intensitydifference momentum maps. The maximum kinetic energy at which the photoelectrons could re-encounter the cation [15], 3.17 U p ≈ 15 eV, resulted in a minimum de Broglie wavelength of ∼6 a 0 , justifying the coarse description of the ion's potential energy in section II A. Fig. 3 shows the calculated and measured asymptotic momentum maps. Each map exhibits an intense part in the center (intensity ≥ 0.4, yellow to orange color), that overlaps with two relatively faint circular patterns (intensity < 0.4, purple color) centered on the polarization axis at the peak vector potential, |p Z | = 0 /ω ≈ 0.83. Fig. 3 e, f show the experimentally obtained intensitydifference data sets, which clearly differ for the two alignment scenarios employed. This difference mainly manifests itself in the star-shaped angular structure of the intense central parts of the images. For the x-axis alignment case there is an angular local minimum along the polarization axis, whereas the momentum map for y-axis alignment exhibits a maximum. Generally, the angular modulations are a lot less pronounced for y-axis alignment. Furthermore, both momentum maps show a few sharp radial maxima along the Z axis at momenta around ±0.4, which are easier to see in the linear-scale representation in Fig. 4. The projected out-of-plane angles,

V. RESULTS AND DISCUSSION
are 25 • for x-and 24 • for y-axis alignment. The simulated asymptotic momentum distributions following irradiation through a laser pulse with a distinct temporal peak intensity are depicted in Fig. 3 a, b. The momentum maps for both alignment cases show rich radial and angular structure in the intense central parts. The radial ATI interference patterns are clearly dominant and are overall shifted in momentum between the two cases. In the high-intensity central parts of the images a star-like pattern is visible with an hourglass shape imprinted on it. Here, the projected out-of-plane angles for x-and y-axis alignment are 28 • and 26 • . Fig. 3 c, d show the modeled asymptotic momentum maps with focal-volume averaging, assuming a range of the ponderomotive potential of 0.83 ω. Here, the ATI interferences are almost fully extinguished [50], which results in the star-like angular pattern becoming the dominant feature. The angles Ω are the same as those for the momentum maps at fixed laser intensity.
These phenomena in the momentum maps are linked to their respective individual sources within the complete molecular strong-field ionization process: The inner, intense part of the asymptotic momentum distribution, which constitutes the majority of the total probability density, represents the fraction of the photoelectron wavefunction that leaves the molecule directly, i. e., without close interaction with the cationic potential. Its continuum motion is dominated by the electric field of the laser, giving rise to ATI interferences between adjacent optical cycles. However, the two rather faint, disc-like shapes in the momentum maps result from strong Coulomb interaction of a small fraction of the photoelectron wavefunction with the ionic potential upon return to the cation [51]. The outermost rings of these discs correspond to electroncation recollisions at maximum momentum and carry information on the molecular geometry [52]. Where the inner ellipse and the outer discs overlap in momentum space, holographic interferences can occur, which give rise to the star-like pattern observed here. This holographic fingerprint encodes both, the tunnel exit position and the initial momentum distribution of the photoelectron wave-function [12]. The embedded hourglass, which is especially pronounced in the x-axis alignment case, results from the nodal plane imprinting onto the momentum distribution at birth and thus onto the holographic pattern.
Both the out-of-plane angle Ω and the absolute kinetic energies of the ATI interferences provide a sensitive probe of the alignment-dependent ionization potential. While Ω is, for identical laser parameters, mainly determined by the initial transverse momentum distribution in (1), the energy of an ATI ring, N ω − I p − U p [36], represents a direct link to the in-field ionization threshold for a given number N of absorbed photons.
Due to their invariance towards laser-intensity averaging the sharp maxima in the radial distributions of the experimental momentum maps are ascribed to Freeman resonances, resonance-enhanced multiphoton ionization through Rydberg states [36,50,53,54]. Conceptually, these electronic resonances cannot be grasped by our current semiclassical model.
The predictive strength of the present model, in comparison to the experimental data is very good. This is seen, for instance, by the agreement of the momentum maps from experiment and the semiclassical model incorporating the emulation of focal-volume averaging. All features encountered in the experimental data are qualitatively reproduced with comparable relative intensities. Furthermore, the trend in the out-of-plane angle is reflected by the simulations, predicting a larger angle for the x-axis alignment case, which can be ascribed to a stronger depression of the ionization potential.
The hourglass shape in the x-axis alignment case represents the only structural deviation in a modeled momentum map with respect to its experimental counterpart: From the semiclassical simulation a local maximum in the angular distribution along the polarization axis is predicted, while in experiment a local minimum is observed. This discrepancy is attributed to the nonperfect linear polarization of the experimentally used laser pulses: A decreasing degree of linearity lowers the return probability to the cation and thus the peculiarity of the holographic pattern, which carries a maximum along the polarization axis. However, the imprint of the HOMO's nodal plane, giving rise to a minimum along the polarization axis, remains at least as pronounced. In the limit of circular polarization the node will be maximally distinct [7,22], while the holographic pattern will be fully suppressed. Although the description of initial in-polarization-plane momenta by means of the adiabatic tunneling theory, including the first-order nonadiabatic correction, was established for elliptical polarization shapes close to circular [55], so far there is no equivalent theoretical framework that tackles just slightly elliptical polarization shapes. In future experimental studies special attention should be paid to optimizing the linearity of the laser field. This would facilitate the comparison with simulation results and moreover maximize the re-collision probability with the cation, which would increase the quality of holographic structures and imprints from laser-induced electron diffraction.

VI. CONCLUSIONS
We unraveled the strong-field photoelectron imaging of the prototypical biomolecule indole using a combined experimental and computational approach. Strongly controlled molecules and an experimental technique suppressing laser-intensity-volume averaging enabled the recording of photoelectron-momentum distributions directly in the molecular frame and for a well-defined, narrow spectrum of incident intensities. As a numerical counterpart we developed a novel, highly efficient semiclassical model that builds upon the adiabatic tunneling theory. Both procedures revealed holographic structures in the asymptotic momentum distributions that were found to sensitively depend on the direction of the ionizing field's polarization axis in the molecular frame.
Based on the very good agreement between experiment and theory we are confident to have identified all essential molecular properties that shape the photoelectron wavepacket as it is born at the tunnel exit. Owing to the quantum-chemically exact treatment of the cation during the subsequent continuum dynamics our model is ideally suited for studies of highly complex molecular structures through strong-field ionization and laser-induced electron diffraction. It allows to describe electron-diffractive imaging of biomolecules on femtosecond time scales while offering the opportunity to fully follow the photoelectron's motion along semiclassical trajectories. Our model description takes account of the cation's laser-induced polarization and gives rise to a significantly faster convergence of the asymptotic photoelectron wavefunction than previously described models [31] due to the use of a spherical energy coordinate system. Furthermore, we sketched a clear path that leads to the simplified emulation of focal-volume averaging, enabling a direct comparison with common experimental data.
The applicability of our model is only limited by the feasibility of the quantum-chemical computation of the molecular ion and the validity of the single-active electron approximation. However, the simultaneous ionization from multiple orbitals could be added easily and would come without any additional numerical expense. Furthermore, our model could be used to simulate experiments involving circular polarization, for example measurements of strong-field photoelectron circular dichroism, employing the corresponding framework of the adiabatic tunneling theory with the first-order nonadiabatic correction [55].
The scattering of electrons at the quantum-chemically exact potential with the Kepler-law approximation would also be useful for trajectory-based descriptions of conventional electron diffraction off molecules, especially in the low-energy-electron diffraction (LEED) regime [56]. The equilibrium geometry of the neutral indole molecule in its electronic ground state, the field-free ionization potential, and the permanent dipole moments and polarizability tensors of the neutral and singly charged indole species were computed using GAMESS [37,38] at the MP2/aug-cc-pVTZ level of theory. For the sake of comparability all molecular-frame dependent quantities in this section are given in the principal-axes-of-inertia frame, marked with a prime. However, all computations were performed in the principal-axes-of-polarizability frame of the neutral species, which relates to the inertial coordinate system by a rotation of ∼1 • within the molecular plane. Table II lists the atomic coordinates in the resulting equilibrium geometry, which is in excellent agreement with previously performed calculations at CC2/cc-pVTZ level of theory [57]. The field-free ionization potential was calculated as the difference of the field-free total energies of cation and neutral molecule I The isotropic polarizability of the neutral species, tr(α n )/3 = 14.75 · 10 −3 nm 3 , is in good agreement with previous quantum chemistry results [61].

Appendix B: Field-direction model
The tunnel-exit positions were obtained in the fielddirection model at the HF/aug-cc-pVTZ level of theory using Psi4 [40]. The quasistatic electric potential of the neutral indole molecule, V ( ), was computed along the alignment-dependent molecular axis that coincided with the polarization axis of the ionizing laser field. That is,  TABLE II. Calculated equilibrium geometry of the neutral indole molecule in its electronic ground state. The coordinates are given in the principal-axes-of-inertia frame. All atoms lie in the z = z = 0 plane. The same atomic coordinates were used for the description of the indole cation during the semiclassical computations.
FIG. 5. The sum of quasistatic molecular potential, laser potential, and ionization energy is shown for a) x-axis and b) y-axis alignment; each panel depicts a single orientation. Only positions on the one principal polarizability axis were evaluated that coincided with the polarization axis, which is depicted as white line within the corresponding molecule sketch. The deep potential wells near the center of nuclear charge reflect the proximity to the closest nuclei. For comparison, the respective tunnel exit positions for a quasistatic laser-electric field and neglecting the Coulomb distortion of the barrier [41] are shown as dashed black lines.
only positions on the x axis were evaluated for the x-axis alignment case and only positions on the y axis for y-axis alignment. In both cases the sampling line was chosen to cross the molecule's center of nuclear charge. Position and external electric field were sampled in intervals of ∆r = 0.5 a 0 and ∆ = 1.4 × 10 −3 . Subsequently, the sum of quasistatic potential, laser potential, and fielddependent ionization energy, V ( ) + Z · Z + I p ( ), was examined, which assumes a value of 0 at the tunnel exit. Fig. 5 illustrates the corresponding results for the two alignment scenarios. The actual tunnel-exit position used in the semiclassical computations was eventually determined by finding the root of V ( ) + Z · Z + I p ( ) at the given instantaneous laser-electric field, which appears as a yellow seam in Fig. 5. For the x-axis alignment case and Z ≈ −3 × 10 −2 the maximum of the potential barrier was found to be suppressed below 0, i. e., in the barrier-suppression regime. For these cases the barrier maximum was used as tunnel exit position and the energy difference to 0 was assumed to be transformed into additional longitudinal momentum opposing the direction of the instantaneous laser-electric field.

Appendix C: Cationic electric field
The external-field dependent cationic field of the indole molecule, − ∇V ( , r), was calculated at the HF/aug-cc-pVTZ level of theory using Psi4 [40]. Positions were sampled in intervals of ∆r = 0.2 a 0 and the laser-electric field was sampled in steps of ∆ = 2.8 × 10 −3 over a (60 a 0 ) 3 cube centered at the cation's center of nuclear charge.