Site-specific tunnel-ionization in high harmonic generation in molecules

We demonstrate that the standard picture of strong-field tunnel-ionization from molecules should be reformulated. The extended nature of the molecular potential implies the separation of some of the molecular sites from the edge of the ionization barrier. We show that the dependence of the tunnel probability with the distance to the barrier is translated into the ionized wavepacket, modifying substantially the high-order harmonic emission. The introduction of the dependence of tunnel ionization with the molecular site significantly improves the theoretical description of high-order harmonic generation in molecules, which is used as a cornerstone in high-harmonic spectroscopy and attosecond imaging.


Introduction
High-order harmonic generation (HHG) is an extreme non-linear process induced by intense fields. In atomic or molecular gases, it can be described as a three-step process [1,2]: first, near the maxima of the driving field's amplitude, an electronic wavepacket is tunnel-ionized from the parent atom; in the second step, the electronic wavepacket is accelerated and, after reversal of the sign of the electric field, it is redirected to the parent ion; finally, upon recollision, the electron's kinetic energy is emitted in the form of high-frequency radiation. The harmonic spectra encode information of the target structure and dynamics, that can be disentangled using high-harmonic spectroscopy (HHS) and time resolved attosecond spectroscopy techniques [3]. These procedures have been successful in retrieving the information from the HHG spectra about molecular structure [4][5][6], nuclear dynamics [7], molecular orbitals [8,9], energy dispersion in solids [10], dynamics in strongly correlated systems [11], tunneling times [12], and orbital tomography [13][14][15].
The study of molecular systems interacting with strong laser fields has a main tool in computationally solving the time-dependent Schrödinger equation (TDSE). The detailed physics, however, is frequently hidden by the complexity of the processes involved. To this end, approximated models are needed to describe the problem in terms of fundamental physical mechanisms. Among them, those based in the strong-field approximation (SFA) [16][17][18] allow to establish a link between the harmonic spectral signatures and the structural details of the radiating matter system. The SFA has successfully demonstrated to reproduce the main characteristics of the harmonic spectra in atoms [19][20][21], and molecules [13,[22][23][24][25][26], and their electron and nuclear dynamics [7,[27][28][29][30][31].
A key ingredient in HHS is the mapping of the molecular orbital matrix elements into the harmonic spectrum and phase, as predicted from the SFA. This relation allows for the reconstruction of the molecular geometry, as well as the molecular orbital structure. The quality of the retrieval of the molecular features relies, therefore, on the fairness of the SFA in describing the HHG process. The first step of HHG is the tunnel ionization, which constitutes a fundamental process in strong-field interactions and has been extensively studied in atoms [16,[32][33][34] and extended to molecules [35][36][37][38][39] in the SFA. Notably, the SFA neglects the shape of the molecular potential barrier at ionization, thus, describing the ionization from each molecular site on equal foot. However, for the actual potential shapes, the elimination of the ion potential barriers inside the molecule leads to extended molecular orbitals, where part of the electron wavefunction is effectively separated from the edge of the tunnel ionization barrier. Recently, Liu and Liu have suggested treating the ionization from each atomic site differently, though by considering different effective barrier widths [39]. Finally, Labeye et al have studied the implications of the internal molecular barrier in the semiclassical description of HHG from molecules [40].
In this work, we demonstrate that the standard tunnel ionization picture has to be revised in the molecular case. As we will show, the tunnel ionization probability is affected by the separation of part of the electronic wavefunction from the potential barrier, introducing a molecular-site specificity in the electron ionization. As tunnel ionization is a key ingredient in HHG, the signature of its site-specificity is evidenced in the details of the high-order harmonic and photoelectron spectra. The total ionization rates, as given by the molecular Ammosov-Delone-Krainov theory [35] and the weak-field asymptotic theory [37], are little affected. As a practical consequence, the orbital information encoded in the harmonic spectrum is substantially distorted from the prediction of the standard SFA models. We, however, show that the proper description of electron tunneling from molecular orbitals can be incorporated into these existing models with the introduction of a modified molecular form factor. The consequences of this study are two-fold: on the one side, it modifies the interpretation of the harmonic signal in HHS, since the retrieved orbital is not a raw image of the actual molecular wavefunction; on the other side, it demonstrates that the HHG spectrum shows well-resolved signatures of the tunnel-ionization site-specificity. Finally, our considerations lead to a derivation of a molecular SFA HHG model with quantitative accuracy.
The article is organized as follows. First, we describe our theoretical approach to the site-depending tunneling ionization and its incorporation in the strong field approximation formalism. Second, we present the results for HHG and photoelectrum spectra obtained from our model compared to those from the exact TDSE calculations, and finally we conclude.

Theoretical approach to site-dependent tunnel ionization in molecules
We consider as the subject of study the hydrogen molecular ion (H + 2 ) in a linearly polarized laser field, since the exact TDSE calculations are affordable for this simple case [41]. The molecule is assumed to lie parallel to the polarization direction, a feasible scenario after using laser alignment techniques [42]. We also neglect nuclei dynamics, which is reasonable for time scales of few femtoseconds [43]. We consider an 800 nm wavelength laser pulse described as E(t) = E 0 sin 2 πt/τ sin (ω 0 t), ω 0 being the laser carrier frequency, τ = 10.67 fs, which leads to 3.88 fs full-width at half maximum (FWHM) in intensity, and a peak intensity of 3.5 × 10 14 W cm −2 . The Hamiltonian governing the interaction is where m is the electron mass, and is the Coulomb potential in cylindrical coordinates (q is the electron charge, and R is the internuclear distance). p is the kinetic momentum operator, p = −i ∇ − (q/c)A(t)e z , A(t) being the vector potential of the laser field, polarized in the z direction. As shown in appendix A, the ground-state of H + 2 is a symmetric orbital that can be described in terms of atomic orbitals centered at the molecular ion sites, as where χ 0 (r) is a localized wavepacket, described as a linear combination of Gaussian orbitals [44,45], and C is the normalization factor. The harmonic spectrum is calculated from the Fourier transform of the mean dipole acceleration along the z axis, a z (t) = â z = −(1/m)∂V M /∂z . The TDSE is integrated using the Crank-Nicholson algorithm in the finite differences scheme. Molecular HHG can be described using an extension of the atomic SFA models. We use the SFA + approach, which provides a quantitative accurate reproduction of the harmonic spectra in atoms [21]. The extension to molecules leads to the following expression for the dipole acceleration (see appendix A): where 0 is the bound state energy, η(P, t 1 ) is the transition matrix element (see equation (A15) in is the molecular form factor, arising from the application of the translation operator to the different molecular ion sites [46]. P is the electron's canonical momentum, S(P, t, t 1 ) is the ionized electron's action, |P are free-electron wavefunctions [47] and V F (P, t) = −(q/mc)A(t) · P + (q 2 /2mc 2 )A 2 (t). We shall refer to equation (4) using the form factor F s [p (t)] as the standard SFA (s-SFA). According to reference [4], p is interpreted as the kinetic momentum of the recolliding electron, as seen from the bottom of the molecular potential well, i.e. p z (t) = sign {p z (t)} p 2 z (t) + 2mI P , p z (t) being the kinetic momentum of the recolliding electron outside the molecular well, and I P being the bound-state ionization energy. This description has been validated by exact calculations using Coulomb-Volkov wavefunctions [48]. Finally, the term α + (P) = (P 2 /2m − 0 )/Δ s in equation (4) is a prefactor necessary to describe the harmonic yield with quantitative accuracy [21], being Δ s the ground-level Stark shift at the instant of recollision. As it is well-known, the form factor F s [p (t)] in equation (4) is responsible for the interference pattern in the harmonic spectrum, dubbed as structural minimum [4].

Introduction of the site-depending tunneling
As pointed out above, a main assumption in SFA is to neglect the form of the Coulomb potential upon ionization and, therefore, it oversimplifies the nature of the tunnel barrier [33]. In figure 1 the molecular potential is depicted showing that the inner potential barrier is nonexistent. As a consequence, one of the wavepackets χ 0 in equation (3) is effectively separated from the outer tunneling barrier, reducing its ionization probability. According to the derivation in appendix B, the SFA description of molecular ionization should be modified to include this site-dependent tunnel probability. To incorporate it, we reformulate the decomposition of the molecular orbital, equation (3), as where T(z c , z b ) is the ratio between the tunnel amplitude probability of the wavepacket placed at the z c molecular site, P(z c , z b ) 1/2 , to the amplitude probability near the z b edge of the potential barrier, P(z b , z b ) 1/2 . The positive (negative) sign of z b corresponds to ionization through a barrier located at the right (left) side of the molecule. According to appendix B, T(z c , z b ) in equation (5) modulates the SFA ionization, so that the wavepackets that are located further away from the barrier have less probability to be ionized, as depicted schematically in figure 1.
Thus, due to this different probability of tunneling, the wavepackets ionized from each molecular site have a different amplitude. The site-specific tunnel probability, T(z c , z b ), can be included in equation (4) redefining the s-SFA form factor, F s , as a new site-dependent tunneling (SDT) version. This new form factor not only accounts for the application of the translation operator but also for the consequences of site-dependent tunneling in the ionized wavepacket: Note, however, that the modulus of the F ± SDT form factor is the same regardless the position ±R/2 of the barrier. We can, therefore, drop the superindices (±) when substituting F s by F SDT in equation (4). After this substitution, we shall refer to equation (4) as the site-dependent tunneling SFA (SDT-SFA).
The interpretation of the site-depending tunneling is probabilistic, as given in appendix B, right after equation (B6). Basically, equation (B5) is the compound probability of two independent events: (i) the particle being located at the edge of the barrier and (ii) the tunnel ionization of a particle in contact with the internal edge of the barrier. The final probability is the product of the probabilities these two events.
The compound event appears intrinsically in the TDSE calculations, but it is not well described within the strong field approximation models, as they neglect the shape of the potential, and, therefore the information of the location of the inner edge of the tunneling barrier.

2
In this section, we show the comparison between the HHG spectra obtained from the TDSE, the standard SFA and the site-dependent tunneling SFA using the form factor F SDT in equation (7). Second, we show the comparison of the HHG spectra calculated from the s-SFA, the SDT-SFA and the TDSE, for different laser parameters and different molecular internuclear distances, showing a much better agreement of SDT-SFA with the TDSE. Next, we extract F SDT form factor directly from the TDSE calculation, and compare it with equation (7), showing that the TDSE results are consistent with our interpretation. Afterwards, we present results where macroscopic propagation in a molecular target is taken into account, demonstrating that the spectral signature of the site-dependent tunneling is resilient to propagation. Finally, we show the suitability of the SDT-SFA in the photoelectron spectrum emitted in H + 2 . Figure 2(a) shows the quantitative comparison between the HHG spectrum obtained with the exact TDSE (blue), the s-SFA (red), and the SDT-SFA (green) as raw data (no rescaling is done). The spectra show the typical feature of non-perturbative harmonic generation: a plateau followed by a cut-off. The s-SFA shows a qualitative recovery of the harmonics close to the cut-off (though, one order of magnitude quantitative error), but a serious qualitative and quantitative departure (up to four orders of magnitude) in the plateau, most evident for harmonic orders below the 45th. These trends are also present in calculations for different laser parameters, as it will be shown below. In contrast, this clear departure is not found in HHG from atoms. To illustrate this, we show in figure 2(b) the same comparison as figure 2(a) but for the helium atom, with an ionization potential similar to H + 2 . The helium case corresponds to equation (4), replacing F s [p (t)] by the atomic form factor, which is equal to one, and substituting the matrix elements and orbital energies accordingly. Figure 2(b) shows that the agreement between atomic SFA and TDSE is excellent even for harmonics well into the cutoff (in this case, for harmonic orders above the 30th). The departure of the SFA description for the lowest frequencies (below 30th harmonic in this case) is a known artifact of the strong field approximation. Thus, while the SFA offers an accurate description of HHG in atoms, it substantially fails to describe HHG in H + 2 . Such departure reveals that some relevant information is missing in the s-SFA formulation, equation (4).

High harmonic spectra: single-molecule response
The presence of the sharp minimum in the s-SFA spectra has been reported before [49,50], and it corresponds to the molecular structural minimum mentioned above. For H + 2 at equilibrium internuclear distance, interacting with an 800 nm-wavelength laser pulse, the structural minimum is centered at the 22nd harmonic, extending up to the 49th harmonic (as found when imposing the condition λ dB = 4R/3, where λ dB is the deBroglie wavelength of the recombining electron). However, in contrast to the s-SFA results, the exact integration of the TDSE shown in figure 2(a) shows a weak trace of this interference. This attenuation of the structural minimum in H + 2 at the equilibrium internuclear distance, parallel to the laser polarization, has been also evidenced in previous works [4,6,50]. On the other hand, it is also known that HHG spectra of the atom of helium for the same laser parameters. In the case of atoms, the SFA approach (red) gives an accurate quantitative reproduction of the exact TDSE (blue) for harmonics above 30th. Both panels show raw data since no rescaling has been done. the minimum shows up for tilted molecular orientations [6,25], where the projection of the molecular axis onto the field polarization results into an effective internuclear distance smaller than the equilibrium one. We will later show results of the TDSE for such case, where the structural minimum is clearly observed.
The comparisons of the results of SDT-SFA with s-SFA and the exact TDSE are shown in figure 2(a). In agreement with the TDSE, the SDT-SFA spectrum shows a weaker signature of the structural minimum. Remarkably, also the overall correspondence of SDT-SFA with the TDSE-both quantitative and qualitative-is substantially improved in comparison with the s-SFA. It should be reminded that the departure at the lowest part of the spectral plateau (harmonic orders <30th) is also found in the atomic case (figure 2(b)). As discussed before, it reflects a fundamental inaccuracy of the SFA approach, not connected to the atomic or molecular nature of the species.
We validate the SDT-SFA model for different laser parameters and internuclear distances. First, in figure 3 we present the HHG spectra considering different laser peak intensities and wavelengths. As a general conclusion, the SDT-SFA (green lines) reproduces very satisfactorily the exact TDSE results (blue lines).
On the other hand, tilted molecules present an effective internuclear distance smaller than the equilibrium one: R eff = R eq cos θ. Under such configurations the ratio of the tunnel probability of a distant wavepacket (equation (6)) increases, and the structural minimum shows up. Thus, we have chosen a smaller internuclear distance in order to study the capability of the SDT-SFA to reproduce the well-known structural minimum [4] in a configuration where it is evidenced. In figure 4 we present the HHG spectra and its corresponding time-frequency analysis for a different internuclear distance, R = 0.6Å, which would correspond to a rotation of θ = 53 • . The minimum appears in the TDSE spectrum at the harmonic order q = 61, approximately. While the standard SFA fails in predicting the minimum's depth, the SDT-SFA provides a better approach. Note that the time-frequency description of the SDT-SFA has also a better coincidence with the TDSE results. As pointed out before, the departures for harmonic orders below 30th are a consequence of the SFA, also present in the atomic case ( figure 2(b)).

Extraction of the form factor from the harmonic spectra
A strong evidence of site-specificity in molecular tunnel ionization can be found directly from the exact TDSE using high-harmonic spectroscopy techniques to retrieve the form factor from the TDSE and compare it with our SDT proposal, equation (7). For that purpose, we follow the same philosophy used in tomographic studies and we compare the HHG molecular emission to that from an atom with similar   ionization potential. Figures 5(a) and (b) show the time-frequency analysis corresponding to the TDSE harmonic spectra of figures 2(a) and (b), for H + 2 and He respectively. These maps reveal that each harmonic in the plateau is emitted in a discrete series of bursts, corresponding to the rescattering of the different electron trajectories, the well-known short and long trajectories [19,51]. We have selected the harmonic emission corresponding to a single rescattering event emitted from the short trajectories during a half-cycle of the driving field (as indicated with white dashed circles in figures 5(a) and (b)) and calculated the associated spectral amplitudes for both species: a z,He (ω) and a z,H + 2 (ω). Then, the form factor can be retrieved as |F TDSE (ω)| 2 a z,H + 2 (ω)/a z,He (ω).
This assumption follows from the SFA saddle point method, which applied to the integrals in equation (4) leads to a simple identification of the harmonic spectral amplitudes and the quasiclassical electron trajectories [19]. After the saddle-point analysis, the Fourier transform of equation (4) can be cast into a simple expression [52]: where η includes the terms containing the electronic wavefunctions, χ 0 , and ξ includes the remaining terms, except the form factor F TDSE , which is to be determined from the TDSE. The summation in equation (9) extends over all saddle points, st, each representing a recolliding electron trajectory responsible of the harmonic emission at frequency ω = p 2 st /2m . Each rescattering corresponds to a term in the summation in equation (9), therefore, in selecting a single rescattering event from the TDSE time-frequency maps, we are isolating a single term in the summation. Thus, in this case, equation (8) follows directly from equation (9). Figure 5(c) shows the comparison of the form factor extracted from the TDSE, F TDSE , using equation (8), against the s-SFA form factor, F s , and the modified SDT-SFA version, F SDT . The excellent agreement between F SDT and F TDSE strongly supports equation (7) and, therefore, our interpretation of distant tunneling as well as molecular site-dependent ionization. We note that evidences of the site-dependent ionization due to the distant tunneling may also be found implicity in photoionization studies where electron localization leads to an asymmetric molecular dissociation [53][54][55].

Signature of the site-dependent tunnel-ionization in the propagated high-harmonic signal
We introduce the macroscopic picture in order to gain insight about the survival of the features of the single-molecule spectrum in an experimental situation. It is known that macrosocopic phase-matching can strongly influence the HHG spectrum measured in an experiment [56,57]. In previous works, the experimental HHG spectrum has been used to retrieve the molecular orbital using tomographic techniques [13]. Therefore, it is important to know if the SDT features of the single-molecule HHG spectrum are also found in the macroscopic picture.
The macroscopic simulation of HHG in H + 2 is based on the electromagnetic field propagator [58], in which we discretize the target (molecular gas jet) into elementary radiators. The dipole acceleration of each elementary source is computed using the SFA models described by equations (A21) and (B10). We assume that the harmonic radiation propagates with the vacuum phase velocity, which is a reasonable assumption for high-order harmonics. The low-density molecular gas jet, flowing along the perpendicular direction to the beam propagation is modelled as a Gaussian distribution of 200 μm at full width half maximum, and with a peak pressure of 5 torr. The driver field has a Gaussian profile with a beam waist at the focus position of 120 μm. All molecules are assumed to be oriented parallel to the polarization direction.
In figure 6, we show a qualitative comparison between the single-molecule and macroscopic HHG spectra for the SDT-SFA and the s-SFA models. The main features of the single-molecule spectrum calculated using both models are conserved in the macroscopic case and, although the shape of SDT-SFA spectra suffers modifications upon macroscopic propagation, it remains clearly distinct from the s-SFA results. Our results show that the main properties of the single-molecule HHG spectrum shown in figure 2 survive after macroscopic propagation.

Signatures of the site-dependent tunnel-ionization in the photoelectron spectra
In this section, we calculate the photoelectron spectra from the H + 2 molecule along the z coordinate, |δΨ(P z )| 2 , after the interaction with the pulse. In figure 7, we present the photoelectron spectra calculated from the TDSE (a), standard SFA (b) and SDT-SFA (c) models. We show that, while the SDT-SFA spectrum exhibits a better agreement with the TDSE spectrum, the standard SFA spectrum presents two deep minima [59]. The origin of those minima is the interference of the ionization from each atomic orbital. The SDT-SFA, on the other hand, includes the side-dependency of the ionization, which sharply decreases the structural interference.
Note that, in order to obtain the photoelectron spectrum from the TDSE, we have filtered the fundamental state from the final wavefunction. As a result, an artificial minimum at the lowest kinetic energies appears in the spectrum.

Conclusions
In conclusion, we have demonstrated that the standard picture of tunnel ionization needs to be modified for the non-atomic case. Our exact computations of the TDSE in H + 2 reveal that wavepacket portions located at the ion sites separated from the barrier ionize with lower probability. We propose a corrected molecular form factor to implement into the existing strong-field models. The new form factor agrees extremely well with the one extracted from the exact TDSE solution, and it improves both the HHG and photoelectron spectra. In addition, the signatures of the site-dependent tunneling are present in the HHG spectra for different laser parameters and molecular internuclear distances. We show that those spectral signatures are also present when the macroscopic phase-matching is taken into account. We believe that the implementation of the site-dependent corrections in the retrieval algorithms will improve substantially the accuracy of high-harmonic spectroscopy measurements, as well as tomographic orbital image reconstruction.

Appendix A. Standard SFA description of molecular HHG
In this appendix, we present the standard (s-SFA) description of molecular high-order harmonic generation (HHG) within the strong-field approximation (SFA).
The interaction of a one-electron system with electromagnetic radiation, in the dipole approximation, is governed by the Hamiltonian where V C (r) is the Coulombic potential and V F (P, t) = −(q/mc)A(t) · P + (q 2 /2mc 2 )A 2 (t), being A(t) the vector potential of the laser field, P the canonical momentum and m and q the mass and the charge of the electron, respectively. The standard SFA approach is based in the assumption that, once ionized, the electron dynamics is governed by the field interaction, neglecting the effect of the ion Coulombic potential. The SFA propagator is written as [18] G SFA (t, t 0 ) = G 0 (t, t 0 ) where G 0 is the propagator of the non-interacting Hamiltonian, H 0 = p 2 2m + V C (r), and G F is the propagator of the free electron in the electromagnetic field, described by H F (t) = p 2 2m + V F (t). Within SFA, the electron's wavefunction is given by |ψ(t) = iG SFA (t, t 0 )|φ 0 , with |φ 0 an initial bound-state of the system. We can split the wavefunction into two terms, |ψ(t) = |φ 0 (t) + |δψ(t) , with |φ 0 (t) = iG 0 (t, t 0 )|φ 0 the bound electron evolving in the absence of the field, and the electron in the continum. Defining 0 as the bound-state energy, we have G 0 (t 1 , t 0 )|φ 0 = −ie −i 0 (t 1 −t 0 )/ |φ 0 . For the electron in the continuum, we resort to the Volkov basis [47]: |P(t) = e i 1 S(P,t,t 0 ) |P , with S(P, t, t 0 ) is the action defined as S(P, t, t 0 ) = − 1 2m t t 0 p 2 (τ )dτ , with p(t) = P − (q/c)A(t), the kinetic momentum. Accordingly, the free electron propagator can be expressed in the Volkov basis as G F (t, t 1 ) = −i e i 1 S(P,t,t 1 ) |P P|dP. Using these definitions in (A3) the ionized electron wavefunction in momentum space reads as where η(P) is the transition matrix element where the factor C F /r n is a Coulomb correction [33] that improves the quantitative accuracy of the SFA description [21,60], with C F = [4| 0 |/(|q|E 0 )] 2 and n = Zq 2 / m/2| 0 |, Z being the charge of the atomic or ionic core (Z = 1 for the hydrogen atom and Z = 2 for H + 2 ). The coherent radiation spectrum is proportional to the Fourier transform of the mean acceleration. Since we are interested only in high-harmonics, we compute the complex acceleration amplitude as [21] being the contribution of each Volkov wave to the total acceleration. The prefactor of the integrand in (A6) accounts for the boundstate dressing at the instant of recollision, necessary to describe the harmonic yield with quantitative accuracy, 0 being the energy of the bound orbital and Δ s the Stark shift at the instant of recollision [21,61]. This formulation has been previously applied to atoms successfully, and an example is shown in figure 2(b).
Let us now introduce the extension of this SFA description to molecules, named as standard SFA (s-SFA) in this paper. We describe the molecular hydrogen ion ground state as a linear combination of atomic orbitals (LCAO), with two centers at ±R/2, R = 2 a.u being the internuclear distance. The LCAO are described with a 6-311G Pople basis, and determined from a variational calculation. For the case of H + being g s (α, r) = (2α/π) 3/4 e −αr 2 . By using Hartree-Fock, the molecular orbital for the ground state is found to be: where C is the normalization factor and χ 0 (r) = 0.1937 Φ 1s (r) + 0.3990 Φ 2s (r) + 0.0484 Φ 3s (r).
Thus, the time-dependent molecular orbital can be expressed as where |χ 0 (t) = e −i 0 (t−t 0 )/ |χ 0 and the binding energy is calculated to be 0 = −29.65 eV (−1.09 a.u.). Consequently, the ionized wavefunction can be written as Figure 8. Scheme of the dipole energy shift of the potential barrier's maximum due to the molecular configuration. An atomic potential (blue) is displaced R/2 from the origin (purple) under the presence of an external field (red). As a consequence, the height of the potential barrier (green) seen by the electronic wavepacket (whose ionization potential is depicted in black) decreases in an amount of E(t 1 )R/2.
wherep z (t) is the projection of the kinetic momentum operator, as seen from the molecular well (see below). The expression for |δψ(t) in equation (A4) should be corrected to take into account the dipole energy shift of the barrier's maximum, as the molecule is an extended object. For a symmetric aligned molecule of length R, the change in the barrier height affects the ionization probability by a fraction W ADK [ (R/2, t 1 )]/W ADK ( 0 ), where (R/2, t 1 ) = 0 − q|E(t 1 )|R/2 is the ionization potential as seen from the top of the barrier at the ionization time t 1 (see figure 8). W ADK is the Ammosov-Delone-Krainov the ionization rate [34] where e is the Euler number, m and are the quantum numbers of the atomic orbital, E 0 is the amplitude of the electric field and We implement this correction including the probability ratio into the bound-to-continuum amplitude probability, thus, redefining η in equation (A4) as where R is the equilibrium internuclear distance (R = 1.055 Å = 2 a.u.). Equations (A7), (A11) and (A12), lead to the following expression for the dipole acceleration for the H + 2 molecule, V M where is the Coulomb molecular potential. For H + 2 , V M can be written as a superposition of the hydrogen potential at the ionic sites, V M (r) = V C r − R 2 e z + V C r + R 2 e z . The dipole acceleration a(P, t), therefore, results from the added contributions of eight different physical paths (see figure 9). Naming {α, β, γ} the sign of the displacements (+1 for right and −1 for left) of the recombination wavefunction χ 0 (r − αR/2), the rescattering potential V C r − βR/2 , and the ionizing wavefunction χ 0 (r − γR/2), respectively, we rewrite (A16) as the sum over the different paths a(P, t) = α,β,γ a α,β,γ (P, t), where Our calculations show that the main contributions come from terms with α = β, meaning that the scattering on a potential site is most probably followed by a recombination to the same site's bound wavefunction. Thus, using we can approximate equation (A17) to For a free electron p (t) = p(t) = P − (q/c)A(t), with A(t) the electromagnetic vector potential. In molecules, however, p (t) describes the kinetic momentum of the free electron at the instant of recombination, t, as seen from the molecular sites, so it includes the acceleration by the potential well, i.e.
(1/2m)p 2 (t) = (1/2m)p 2 (t) + I P , I P being the molecular ionization potential. Therefore p z (t) = sign{p z (t)} p 2 (t) + 2mI P [48]. This correction to the free electron's kinetic momentum is found necessary to recover the correct position of the molecular structural minimum in the harmonic spectrum [4]. Summing over all the relevant paths, the acceleration a(P,t) can be written as where is the molecular form-factor. Substituting in equation (A6), we finally obtain the total acceleration

Appendix B. Tunnel of a distant wavepacket
In this appendix, we compute the tunnel probability of a wavepacket located at a finite distance from the barrier's edge. In a description using localized atomic orbitals (LAO), the molecular orbital is decomposed into a basis of wavepackets, each centered at a different ion site. In molecules at equilibrium nuclear distances, these localized wavefunctions are generally not separated by internal potential barriers, therefore Figure 10. Scheme of the tunneling of a distant wavepacket. The wavefunction χ(z, t) (red) is centered at z c and its energy is 0 , while the blue triangle is the potential barrier V(z) whose edges are placed at z b and z f . dwelling in a multi-ion Coulomb potential well. The presence of a strong field modulates the molecular potential, forming an external potential barrier, that separates the bound orbitals from the continuum. The extended nature of the molecular well, therefore, results in some of the LAO positioned at finite, non zero, distances from the edge of the barrier (see figure 1 in the main text). Let us consider a stationary wavepacket with energy 0 centered at the coordinate z c , at the left of a potential barrier V(z). The barrier's edges z b and z f are defined so that V(z) = 0 if z < z b , and V(z) > 0 if z b < z < z f (see figure 10).
We express the wavepacket as a planewave decomposition, χ(z, t) = e −i 0 t/ χ 0 (z), with χ 0 (z) = ∞ −∞g (p)e iS(p,z)/ e ip(z−z c )/ dp ∞ −∞g (p)e iS WKB (z)/ e −ipz c / dp, ( B 1 ) where we have defined S WKB (z) = z −∞ √ 2m [ − V(z )]dz as the semiclassical approximation ( → 0) to the planewave phase S(p, z) + pz, according to the 0th-order WKB approximation [62]. Note that it is crucial to preserve the localized nature of the wavepacket, therefore the 0th order WKB approximation is used at the level of the individual planewaves composing it, rather that to the total wavefunction. We can rewrite equation (B1) as The tunnel transmission probability of the particle located at the wavepacket's mean position z c , is therefore given by where W WKB corresponds to the WKB probability of tunneling of a particle located at the left edge of the barrier [62], Therefore, equation (B5) can be interpreted as the probability of tunneling times the probability of the particle being near the edge of the barrier. The quotient in equation (B5) shows that the tunnel probability is reduced when the wavepacket's center is at a distance z b − z c from the edge of the barrier. We define the ratio of the tunnel amplitude probability of a distant wavepacket to the amplitude probability of the particle being at the barrier's edge is given by In the case of the H + 2 molecule, the wavepackets χ 0 are atomic orbitals localized at the molecular ion sites, z c = ±R/2, where R is the internuclear distance. The coordinate of the barrier's edge, z b , corresponds approximately to −R/2 (R/2) for a positive (negative) field amplitude. To introduce the site-dependent probabilities in the SFA formalism it is sufficient to redefine the molecular form factor, F s p z (t) , in the s-SFA formula (A21) as where the positive (negative) index corresponds to ionization through a barrier at the right (left) side of the molecule. Note however that the modulus of the form factor is the same, regardless the right/left position of the barrier, and therefore we can drop the superindex, Therefore equation (A21) finally becomes e i 1 S(P,t,t 1 ) e −i 0 (t 1 −t 0 ) V F (P, t 1 ) P| C F r n |χ 0 dt 1 dP.
The new form factor, F SDT , defines the difference between the standard SFA (s-SFA) and the site-dependent-tunneling SFA (SDT-SFA) proposed in this work.