{'fraction_non_alphanumeric': 0.057414670783610115, 'fraction_numerical': 0.01922167515090823, 'mean_word_length': 3.9381868131868134, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
First-order strong field approximation for high-order harmonic generation

Recently it was shown [A. Gordon and F. X. Kaertner, Phys. Rev. Lett. 95, 223901 (2005)] that the strong field approximation (SFA) for high-order harmonic generation (HHG) is significantly improved when the SFA wave function is used with the acceleration rather than the length form of the dipole operator. In this work it is shown that using the acceleration form upgrades the SFA from zeroth-order to first-order accuracy in the binding potential. The first-order correct three-step model (1st-order TSM) obtained thereby is systematically compared to its standard zeroth-order counterpart (0th-order TSM) and it is found that they differ significantly even for energetic electrons. For molecules (in the single-electron approximation), the 0th-order and the 1st-order TSMs in general disagree about the connection between the orbital symmetry and the positions of the minima in the HHG spectrum. At last, we briefly comment on gauge and translation invariance issues of the SFA.


I. INTRODUCTION
The strong field approximation (SFA) [1,2,3] is a key technique in the study of interactions of matter with intense laser fields. In particular, the SFA is used to describe high-order harmonic generation (HHG) [4,5,6,7]. The three-step model (TSM) [5,6,7], which is based on the SFA, has proven very successful describing much of the experimental behavior. The best known examples are the cutoff formula [5,6,7] and time-frequency structure of the HHG signal [8], which led to the prediction of attosecond pulses [9]. Reviews can be found in Refs. [10,11].
It has been long known however that the TSM is incorrect by 1-2 orders of magnitude predicting the spectral intensity of HHG in atomic hydrogen, as found from comparisons with numerically-exact results [11,12]. For the H + 2 ion the TSM gives a spectrum which is many orders of magnitude away from numerically-exact results [13], and the shape of the spectrum can be heavily distorted [13,14].
HHG experiments are now coming to the point where more accurate theory is needed. The race towards harnessing HHG for a coherent short wavelength source [15,16,17] can benefit from quantitative theoretical estimates of the HHG efficiency, which to date are very scarce [18]. The recent orbital imaging experiment [19] uses the HHG spectrum to infer the structure of a molecular orbital, which also requires quantitatively reliable theory, capable of giving a precise description of the shape of the HHG spectrum.
In a recent theoretical work [13] we proposed a modified version of the TSM, where the SFA wavefunction is used with the dipole operator in the acceleration rather than the length form. Comparison of the TSM with a numerical solution of the time dependent Schrödinger * Electronic address: gariel@mit.edu equation (TDSE) then demonstrated excellent quantitative agreement for atomic hydrogen, and significantly improved agreement for the H + 2 ion. The argument why using the acceleration form improves the TSM so much was that the TSM then becomes correct to first order in the binding potential (we shall henceforth refer to the TSM obtained this way as 1 st -order TSM), whereas the standard TSM [5,6,7] (0 th -order TSM) is correct to zeroth order in the binding potential.
This work is a followup on Ref. [13], and it has three goals. First, the first-order accuracy of the 1 st -order TSM is established. Second, a detailed and general comparison between the 1 st -order TSM and the 0 th -order TSM is given. It is shown that the two always disagree with respect to the connection between the orbital symmetry and the positions of the minima in the HHG spectrum, as demonstrated for H + 2 [13,14]. Third, the opportunity of re-deriving the SFA is used to give a perspective on some delicate issues in the derivation, such as gauge and translation invariance, which have been under debate [20,21,22,23,24]. This paper is organized as follows. In Sections II and III the derivation of the SFA and 0 th -order TSM is reviewed, in an attempt to illuminate some subtleties in the derivation and in order to set the stage for deriving the 1 st -order TSM in Sec. IV. Section V compares the 0 th -order and 1 st -order TSMs, and Sec. VI is dedicated to the discussion of gauge and translation invariance issues in the SFA. Section VII gives a brief summary.
Atomic units are adopted throughout the paper.

II. STRONG FIELD APPROXIMATION
This work is restricted to the single-electron approximation, where an atom or a molecule are modeled by an electron in an effective (local) potential V (r): I p is the binding energy of the ground state, and is added to Eq. (1) for convenience reasons, such that the ground state of H 0 has zero energy. The atom or molecule is placed in a linearly polarized electric field E(t). The x axis is chosen along the direction of polarization, the wavelength is assumed sufficiently long such that the dipole approximation holds, and the length gauge is chosen, to give the Hamiltonian The SFA is usually presented [see, e. g. Ref.
[10] and references therein] as a perturbative expansion in V (r), where the unperturbed Hamiltonian is the Volkov Hamiltonian The evolution operator is known exactly [6]. Using the Lippmann-Schwinger equation one can approximate U (t, t ′ ), the evolution operator associated with H(t) (defined through Eq. (4) with all V subscripts omitted) to arbitrary order in V (r) [25]. The zeroth and first order would be With the operator U (t, t ′ ) at hand, any problem can be solved. If at time t = 0 the electron is in the ground state |0 of H 0 , the wavefunction |ψ(t) at any time would be given by One therefore could suggest approximating U (t, t ′ ) in Eq. (7) by the perturbative expansion outlined in Eqs. (5,6). However this is not the way the SFA is usually performed. The route is rather by making the ansatz [6] |ψ = a(t)|0 + |ϕ(t) .
a(t) has the initial conditions a(0) = 1 and is determined later. In other words, the way we split |ψ into a(t)|0 and |ϕ(t) is not yet defined. Let now |ϕ(t) be the exact solution of which in terms of U reads Let us assume for the moment that a(t) is specified. The SFA approximates U in Eq. (10) by a perturbation series in V upon U V . One can now ask why the SFA approximates U in Eq. (10) rather than in Eq. (7), if both lead to an exact wavefunction when U is exact? In other words, why one needs the ansatz (8) and the rather non-straightforward Eq. (9)? To our understanding, the answer to these questions is that the perturbed U V (t, t ′ ) is a very bad approximation for describing the evolution of (quasi) bound states.
Without the laser field, perturbation theory for the evolution operator (the Born series) diverges at bound states of the full Hamiltonian [26]. The presence of a laser field may change the situation, since formally there are no bound states at all. We are not in a position to make a rigorous mathematical statement about the convergence of the SFA, but since the ground state remains quasibound, one may expect convergence difficulties. Obviously, at least for the the zeroth order theory, U V (t, t ′ ) itself gives a very bad approximation for the evolution of a quasi-bound state. In fact, in the Keldysh ionization rate [1] is found rather indirectly, by calculating the rate at which the norm of the continuum increases, and the rate at which the ground-state decays is inferred solely by conservation of the total norm.
The last term in Eq. (9) is therefore added because we know in advance that we are going to use an approximate propagator, which will give a very bad description for the evolution of the ground state. Yet, Eq. (9) provides a clear route for systematically improving the SFA: a(t) is first assumed to be known, and Eq. (9) is solved to in principle an arbitrary order in V (r). Then a(t) is found by demanding the conservation of norm, or by means of other approximations, as the perturbed U V (t, t ′ ) alone is not well suited for this purpose. Note that in the limit where U is exact, one easily finds from Eq. (9) and Eq. (8) that a(t) ≡ 1.

III. ZEROTH ORDER TSM
In this section we re-derive the standard 0 th -order TSM. The main reason is that our slightly different derivation lays the technical foundations for establishing the first order accuracy of the 1 st -order TSM in Sec. IV.
The 0 th -order TSM is obtained by simply replacing U (t, t ′ ) by U V (t, t ′ ) in Eq. (9). Since U V (t, t ′ ) is known exactly, the solution can be written in a closed form [6]: dt ′ is a vector potential that describes the electric field E(t), and At this point we slightly deviate from the standard derivation of the TSM [6,7]: We perform the saddlepoint integration in Eq. (11) right now, before proceeding. This calculation was pioneered by Keldysh [1] and is further discussed in a vast number of SFA studies [see e. g. Ref. [27] and references therein]. The discussion in this work is limited to the tunneling regime, which is defined by the requirements [27]: In Eq. (13a) E 0 is the amplitude of the driving field and ω is its frequency. These parameters are precisely defined for a sinusoidal driving field, and more loosely for a general field, such as a sinusoidal field with an envelope.
In the latter case, ω can be thought of as the parameter characterizing the timescale over which E(t) varies. γ is the well known Keldysh parameter [1]. Note that the requirement (13c) follows from (13a) and (13b). Under the conditions (13), the saddle-point integration of Eq. (11) is carried out to give √ 2Ip e −iS(p,t,tn(px)) . (14) w(E) is the static Stark ionization rate associated with the ground state, p 2 ⊥ ≡ p 2 y + p 2 z , and the function t n (p x ) is the set of positive real solutions of the equation The number of solutions depends on E(t). a(t) can now be found by requiring conservation of the norm in Eq. (8) [neglecting 0|ϕ(t) ] to give the well-known expression [11] Using the ansatz (8) with the wavefunction (14), one can now compute the expectation value of the dipole moment: The origin of coordinates can always be chosen such that the first term in Eq. (17) vanishes. The second term has no high harmonics [28], since matrix elements of x between two different Volkov states vanish. The high harmonics come from the cross term ξ 0 (t). Using Eq. (14) we find: The p integration is now carried out in the stationary phase approximation as in Ref. [6]. The stationary phase is attained at p ⊥ = 0 and at all p x values such that t n (p x ) satisfy [see Appendix A for a more careful discussion] For any given t,t n (t) is defined as the set of solutions to Eq. (20), i. e. birth times of trajectories that end up at the origin at time t. Then in the stationary phase approximation Eq. (19) gives whereS n (t) ≡ S(A(t n (t)), t,t n (t)). Eq. (21) is identical to the original TSM expression derived in Ref. [7]. The only difference is that Eq. (21) takes into account the depletion of the ground state, and that all numerical prefactors (2 3/2 π) were calculated. The rederivation of the 0 th -order TSM is now concluded.

IV. FIRST ORDER TSM
A straight forward way to upgrade the TSM from zeroth to first order accuracy in V would be improving the TSM wavefunction (11) by adding the first-order correction to the evolution operator, Eq. (6). This method has been employed in Ref. [29], for improving the SFA theoretical description of above threshold ionization. For HHG we propose another method, which is significantly simpler, and does not require correcting the wavefunction.
If |ψ(t) is the exact solution of the time-dependent Schrödinger equation with the Hamiltonian (2), then the Ehrenfest theorem holds: Therefore computing the time-dependent dipole expectation value in all three forms of Eq. (22), length (22a), velocity (22b), or acceleration (22c), is equivalent if accompanied by the appropriate differentiation or integration in time.
However with the approximate wavefunction given by Eqs. (8) and (14), the results will be in general different for each form. One therefore has to make a choice which form to use for this particular wavefunction. We argue that the acceleration form gives in general the best results, since even if |ϕ(t) is only correct to zeroth order in V (r), the expectation value is automatically correct, at least formally, to first order in V (r). One therefore needs to evaluate whereξ The last term of Eq. (22c) has been dropped, since it is the driving field itself and thus contains no harmonics. It is easy to show, exploiting the commutator [p, H 0 ], that 0|∂ x V (r)|0 = 0 [irrespectively of the position of origin], which is why the corresponding term in Eq. (23) is missing. Note that the mere appearance of V in Eq. (23) is not sufficient to warrant first order accuracy. One can see that, for example, from the fact that zeroth order SFA expressions in V can be transformed such that they appear linear or quadratic in V [30,31]. The (formal) first-order accuracy of Eq. (23) is established by the fact that a first-order correction in V to the wavefunction in Eq. (11) will result in only a second order correction in V in Eq. (23). This is why we dedicated Sec. II to carefully defining the procedure by which the SFA is corrected order by order in V .
The cross termξ 1 has the same structure as ξ 0 [Eq. (18)]. Going through the same steps that have led to Eq. (21), one arrives at the same expressionξ 1 , with the only difference that the matrix element 0|x|k is replaced by − 0|∂ x V (r)|k : Eq. (26) is the improved version of the TSM presented in Ref. [13]. The change in the expression compared to Eq. (21) is very small and easy to implement, and yet results in a very large difference, especially in the case of molecules.
Eq. (26) is an expression forξ 1 , but Eq. (23) contains also ξ c . It turns out that the contribution ofξ c to HHG is smaller by at least O(ω 3/2 ) than that ofξ 1 .ξ c is therefore negligible for ω ≪ 1, which coincides with Eq. (13c) when I p is of O(1). The latter holds for all neutral (or not highly charged) atoms and molecules. The evaluation of ξ c is rather lengthy, especially for potentials with a longranged Coulomb tail, and is given in Appendix B.
The 0 th -order TSM suggests Eq. (18) as an approximation to the dipole moment, whereas the 1 st -order TSM suggests Eq. (24). In order to compare them conveniently, we now differentiate Eq. (18) twice in time and compareξ 0 withξ 1 . This is done in detail in Appendix C. The result is that under the condition (13a), in order to obtainξ 0 (t), one has to replace the 1 st -order TSM recombination amplitude in Eq. (24) by Comparing the 1 st -order TSM and the 0 th -order TSM is thus reduced to comparing a new rec and a old rec respectively. The two expressions look different, and indeed they are. A detailed comparison requires the knowledge of V . However in order to gain some general insight, in what follows we study the k → ∞ asymptotic behavior of a new rec and a old rec . Some pretty general statements can be made about k → ∞, which turn out to provide important insights also for k-s of O(1).
A. High-momentum asymptotic evaluation a old rec (k) can be written as (and similarly for a new rec (k) and χ new (r)), where with ψ 0 (r) ≡ r|0 . Using H 0 |0 = 0, Eq. (30) can be transformed into The k → ∞ asymptotic behavior of a old rec (k) and a new rec (k) is dictated by the most singular part χ old (r) and χ new (r) respectively, since they are connected by a Fourier transform [40].
We now assume that V (r) is an effective potential that represents a molecule (or, as a special case, an atom) with nuclei at positions r j with charges Z j , j = 1...N . Common choices for an effective potential to model such a system have a Z j /|r − r j | singularity near the j-th nucleus, which is partially screened off away from the nucleus. Then the last term in Eq. (32) can be written as Now one has to distinguish two cases. The first case is a molecule. Whenever there is more than one nucleus, Eq. (33) does not vanish [apart from possibly a very few orientations of the molecule]. Then the term in Eq. (33) is the most singular one among all terms on the right hand side of Eq. (32), and thus governs a old rec (k) at k → ∞. In this case where the notation '∼' means 'equal in its most singular part'. The second case is an atom, which has only one nucleus, whose position can be always chosen at r = 0 [see Sec. VI for a discussion about translation invariance]. In this case Eq. (33) obviously vanishes. Then the most singular term on the right hand side of Eq. (32) is the first one, and we find that

B. Atoms
Eq. (35) implies that for atoms where '→' denotes 'approaches asymptotically at k → ∞'. Thus a new rec (k) and a old rec (k) do not agree even asymptotically for k → ∞, the limit were plane waves approximate the exact continuum states increasingly well. For atoms, the 0 th -order TSM predicts a 4 times larger HHG yield than the (1 st -order TSM). For a Z/r singularity in the potential, the asymptotic behavior of a new rec can be easily worked out to give When ψ 0 (0) vanishes Eq. (37) is modified. For finite k the disagreement between a new rec (k) and a old rec (k) can be even greater. For the hydrogen atom one finds [32] a new rec (k for k lying along the x axis. Fig. 1 compares the two expressions in Eq. (38).

C. Molecules
The differences between the 0 th -order TSM and the 1 storder TSM for molecules is particularly striking. The two models give opposite relations between orbital symmetry and the positions of the minima in the HHG spectrum, a subject of many recent theoretical and experimental studies [33,34,35,36,37,38,39]. A set of minima is predicted by the 0 th -order TSM with an odd orbital, corresponds to an even orbital in the 1 st -order TSM, and vice versa.
From Eq. (34) we find that for molecules whereas from Eq. (27) and Eq. (37) we find  40) is the Fourier transform of xψ 0 (r), whereas Eq. (39) is the Fourier transform of ψ 0 (r) (both sampled at the singular points). xψ 0 (r) and ψ 0 (r) obviously have the opposite symmetry with respect to a x → −x reflection. In order to illustrate this we consider the H + 2 ion, where the nuclei are positioned along the x axis, at x = ±R/2. In this case Eq. (39) gives where C ≡ ψ 0 (− 1 2 Rx) = ψ 0 ( 1 2 Rx). In contrast, Eq. (40) gives A zero of a rec at a given k x corresponds to a minimum in the HHG spectral intensity at the frequency I p + 1 2 k 2 x . The 0 th -order TSM therefore predicts minima at energies given by I p + (2πn/R) 2 (n is an integer), whereas the 1 st -order TSM predicts minima at energies given by I p + (2π(n + 1 2 )/R) 2 . The latter condition well agrees with numerically-exact results [13,39].
It is easy to see from Eq. (39) and Eq. (40) that if ψ 0 (r) were an odd rather than even wavefunction with respect to x → −x, the cosine [sine] in Eq. (41) [Eq. (42)] is replaced by a sine [cosine]. Therefore if the 0 th -order TSM is used for reconstructing ψ 0 (r) from the HHG spectrum, it flips the symmetry of ψ 0 (r) from even to odd and vice versa. This statement is, of course, based only on the k x → ∞ asymptotic behavior of the a rec -s. Yet it is interesting to note that for H + 2 , Eq. (41) and Eq. (42) approximate Eq. (27) and Eq. (28) reasonably well for low momenta, as Fig. 2 shows. Note that the k x /k 2 envelope, which is present in Eq. (40) but absent in Eq. (39), results in an ordersof-magnitude difference in the magnitudes of the a rec -s, in addition to the distorted shape.

VI. GAUGE AND TRANSLATION INVARIANCE ISSUES
The dependence of the SFA on gauge has been the subject of lively discussions [3,20,21,22,23,24]. Here we do not attempt a comprehensive discussion of this subject. However since some steps in Sec. V seem to rely on placing the origin of the laser potential and the atom at the same point, we briefly address the related issue of translation invariance.
Consider thus the Hamiltonian  (43) can give rise to the generation of even harmonics [21], which is an artifact.
There is an obvious way by which this problem can be cured, namely, by replacing the ansatz (8) by Substituted Eq. (44) into the TDSE with the Hamiltonian (43), it is easy to convince oneself that the results no longer depend on x 0 . This is because wavefunction (8) undergoes exactly the gauge transformation that connects Eq. (43) to Eq. (2), guaranteeing gauge invariance. The SFA depends on the choice of gauge simply because it crucially depends on the choice of the initial ansatz.
One can now ask how to obtain an ansatz for a given gauge a-priori, without e. g. gauge-transforming Eq. (8), and how to obtain the ansatz (8) itself a-priori. Our answer is that the ansatz should provide the best approximation for the evolution of the quasi-bound ground state, since the (perturbed) Volkov evolution operator cannot do it [see Sec. II].
In the low frequency regime [Eq. (13c)], a good approximation for the evolution of the ground state can be found using the adiabatic theorem [41]: It is approximately given by e −i ε(E(t))dt |φ(E) , where the |φ(E) and ε(E) are the field-dependent ground state wavefunction and energy, defined through the eigenvalue equation The adiabatic theorem holds even though |φ(E) becomes a resonance [42]. If we neglect the Stark shift of the energy, we obtain ε(E) = Ex 0 . Going further and approximating |φ(E) by |φ(0) = |0 , we obtain exactly the ansatz (44), a-priori, without referring to Eq. (8). In fact, Eq. (8) itself is obtained for x 0 = 0.
The discussion can be extended to at least one common gauge, namely the velocity gauge. Then the adiabatic argument would again give a modified ansatz, similar to the one used in Ref. [43]. It is likely that the differences between the length and velocity gauges reported in Ref. [23] would then disappear.

VII. DISCUSSION
The large discrepancy between the 0 th -order TSM HHG spectra and numerically-exact calculations has been attributed long ago [11] to the SFA wavefunction being inaccurate. It is often argued that this SFA wavefunction is especially inaccurate for molecules [45]. This work strongly supports these ideas: The length and acceleration forms of the dipole operator used on the same SFA wavefunction are shown to give dramatically different results, especially for molecules. For an exact wavefunction the results would be identical, and thus the large discrepancy is evidence for the inaccuracy of the wavefunction.
Since the wavefunction is so far from being exact, it is especially important to use it properly. In this work we have argued that together with the acceleration form of the dipole, the zeroth-order SFA wavefunction leads to an approximation (the 1 st -order TSM) with first-order overall accuracy in the binding potential. Since the wavefunction itself is so inaccurate, it can lead to very large errors if used otherwise than in the specially-suited way introduced in this work. The accuracy of the 1 st -order TSM should be tested in terms of the expectation values it generates and not in terms of the wavefunction itself.
We have shown that the continuum-continuum term in the 1 st -order TSM is negligible compared to the boundcontinuum term. Therefore, the widely-used 0 th -order TSM is upgraded to first-order accuracy by simply replacing x by ∂ x V (r) in the recombination amplitude, which is simple to implement. The 1 st -order TSM shows excellent agreement with numerically-exact results for atoms [13,46] and good agreement for H + 2 [13]. The SFA decomposes the wavefunction into a continuum, and a quasi-bound ground state. This is justified, since the strong laser field smears the excited states of the unperturbed system and turns them into a continuum [6]. The dynamics of the continuum is approximated by the Volkov propagator, which can be corrected order by order in the binding potential. The quasi-bound state however should be treated separately.
It should be noted that there could be more than one quasi-bound state. For example, the initial state of the electron can be an excited state of the laser-free Hamiltonian H 0 , a case that was not treated in this work. This is commonly the case in single-electron models of multielectron atoms. Then the states lying energetically below the initial state are also quasi-bound, and should also be treated separately. One simple way to do that is to project these states out of the SFA wavefunction, as briefly discussed at the end of Ref. [46].
The second term in Eq. (A1) is of the order of E 0 /ω 2 , while the first one is of the order of I p /E 0 . Therefore the first term in Eq. (A1) is a second-order correction in γ to Eq. (20). It is interesting to note however that the subleading term term is meaningful, since the electron is released by tunneling roughly at the turning point (x = I p /E) [27] rather than at the origin. Note that the denominator (t − t n (t)) 3/2 in Eq. (21) cannot result in a divergence. If the first term in Eq. (20) is retained, the reason is obvious -the duration of a trajectory is never zero, since it begins and ends at different locations. If the first term in Eq. (20) is neglected, as often happens, one can easily show that trajectories for which t approaches t n (t) can only occur at points where E(t) = 0, and then w(E) vanishes exponentially. Performing the t ′ integration first thus spares the need for the regularizing parameter ǫ used in, e. g. Ref. [6].
Note also that the denominator (t − t n (t)) 3/2 is correct only to zeroth order in γ. The leading correction is obtained if the e −p 2 ⊥ /(|E(tn(px))| √ 2Ip factor in Eq. (19) is included in the Gaussian in integration and if Eq. (15) is taken into account when the second derivative of S with respect to p x is computed.

APPENDIX B: EVALUATION OF THE CONTINUUM-CONTINUUM TERM
For completeness, in this Appendix the continuumcontinuum termξ c (t) [Eq. (25)] is evaluated. To this aim we define the Fourier transform of the force field derived from the potential V . It turns out that a long-range Coulomb behavior of V requires more careful evaluation ofξ c than a short-ranged V . Since once |ϕ(t) is givenξ c (t) is linear in V , one can write V as a sum of the pure Coulomb potential plus a short-ranged potential, evaluate each term separately, and sum up the results at the end. In what follows we therefore treat the two cases separately.

Short-range potentials
In order to keep the expressions from becoming too cluttered, we introduce the notatioñ Substituting Eq. (14) in Eq. (25) we obtain We now carry out the integration in Eq. (B3) in the stationary-phase approximation. To this aim we assume that apart from the exponentials in Eq. (B3), the rest of the integrand is slowly varying in p. Here it is where we use the short-range property of V , which assures that Q is indeed slowly varying in p.
The stationary phase condition in the six-dimensional p-p ′ momentum space is identical to the one of the integral (19): p ⊥ = p ′ ⊥ = 0, whereas p x and p ′ x are obtained by finding all solutionst n (t) that satisfy Eq. (20) and using Eq. (15) to find the corresponding p x -s. The result of the integration is where the dependence oft n on t has been suppressed. Eq. (B4) confirms our claim thatξ c (t) is small compared toξ 1 (t).ξ 1 (t) [Eq. (26)] has only a (2π/(t − t n )) 3/2 prefactor (which is O(ω 3/2 )), whereas forξ c (t) it is O(ω 3 ). Moreover, ξ 1 (t) is proportional to w(E), whereasξ c (t) is linear in w(E). HHG experiments typically operate under the condition of small ionization per cycle, which means w(E) ≪ ω. Therefore althoughξ c (t) is formally of the same order in V (r) asξ 1 (t),ξ c (t) is more than O(ω 2 ) smaller thanξ 1 (t), and is thus negligible under the conditions (13).
It should be noted that since ξ 1 (t) is proportional to a * (t), at very high fields, where the ground state is almost completely ionized in one cycle, ξ 1 (t) becomes exponentially small in the field amplitude [44]. In this caseξ c (t) may become significant. However this regime is of little interest form the point of view of HHG, since HHG basically disappears under these operating conditions [44].

Coulomb potential -different trajectories
Now we consider the case V (r) = −1/r, which leads to Due to the singularity at k = 0, the stationary phase approximation should be now used more carefully when integrating Eq. (B3). To this aim we now consider an integral of the form (B6) I represents the third line of Eq. (B3), and the rest of the integrand is slowly varying and can be added later. Eq. (B6) has two oscillating phase factors, centered at p 0 and p ′ 0 , and the Coulomb potential. τ and τ ′ are real and positive, and by comparison with Eq. (B3) one can see that they represent the traveling times of the two trajectories.
Using the convolution and Plancharel's theorems, Eq. (B6) is transformed to whereĨ where s ≡ 1/τ ′ − 1/τ and ∆p 0 ≡ p 0x − p ′ 0x . We assumed that the p 0 − p ′ 0 is parallel to the x axis, since this is the case of interest in Eq. (B3).
The integral in Eq. (B8) can be carried out analytically in spherical coordinates and expressed in a closed form using the error function: Eq. (B9) and Eq. (B7) give an exact expression for I in Eq. (B6). Let us now perform the integration in Eq. (B6) using the stationary-phase approximation instead. The result is Since each exponential picks up only an environment of radius τ −1/2 (or τ ′−1/2 ) around its center p 0 (or p ′ 0 ), one expects that when 1/τ, 1/τ ′ ≪ ∆p 2 0 , the integration does not reach the singularity and Eq. (B10) holds. Eq. (B7) and Eq. (B9) verify this expectation. Figure 3 visualizes this statement, and shows thatĨ(∆p 0 , s) basically gives a smoothed version of the 1/∆p 0 singularity near ∆p 0 = 0.
The discussion in Sec. B 1 therefore holds as it is as long as the condition 1/τ, 1/τ ′ ≪ ∆p 2 0 is met. This condition is violated when ∆p 0 approaches zero, and this case will be our concern in what follows. By observing Eq. (B4) one can see that the latter always happens when n = n ′ , and can also accidentally happen if n = n ′ . We begin with the second case. By Taylor-expanding the modulus of Eq. (B9) keeping the two leading orders in ∆p 0 , one can find the maximal value of |Ĩ(∆p 0 , s)| with respect to ∆p 0 , and see that it is proportional to s −1/2 . Using this result for an upper bound on I, one obtains |I| < 2 7/2 4 700 80919 Note that we are considering two different trajectories, which means that τ = τ ′ . Moreover, it is easy to show that two different trajectories that return at the same time must have their birth times separated by more than a quarter of a driving cycle [ω|τ − τ ′ | > π/2]. It follows therefore that Eq. (B11) is O(ω 5/2 ). For two different trajectories, n = n ′ , the Coulomb singularity therefore enhances the integral Eq. (B3) by at most a factor of ω −1/2 . This is seen by comparing Eq. (B7) and Eq. (B11). In Sec. B 1 we have shown thaẗ ξ c (t) is smaller thanξ 1 (t) by a factor of O(ω 2 ). Now we arrive at the conclusion that in the long range case, for two different trajectories,ξ c (t) can be sometimes smaller than ξ 1 (t) by factor of O(ω 3/2 ) only. Yet,ξ c (t) remains negligible for ω ≪ 1.

Coulomb potential -same trajectory
If n = n ′ , which means τ = τ ′ and ∆p 0 = 0, the integral (B6) diverges. This divergence is an artifact of the stationary phase approximation and can be avoided if one uses the fact that the first two lines in Eq. (B3) vanish exponentially at p, p ′ → ∞. Doing so is however technically cumbersome, and we adopt a simpler approach.
We start with having another look at Eq. (B4). For a given pair n and n ′ , the main HHG frequency to be generated is given by the derivative of the phase in Eq. (B4): Eq. (B12) has a simple intuitive meaning.ξ 1 (t) gives the beat frequency between the ground state (at frequency I p ) and a continuum electron with a frequency corresponding to the kinetic energy upon return. In contrast, ξ c (t) is obtained from two different electron trajectories, which return to the parent ion at the same time with two different kinetic energies. The emitted radiation is at the beat frequency corresponding to the difference in the kinetic energies upon return. Obviously, if n = n ′ , that is, if the two trajectories originated at the same birth time, they are identical. Therefore the beat frequency is zero and no high harmonics are generated. In other words, for n = n ′ Eq. (B3) varies slowly in time. There is yet a subtle point to be checked: loosely speaking, if Eq. (B3) is very large, then even if it varies very slowly in time, its high-frequency tail may be comparable toξ 1 (t).
In order to show that this is not the case, we introduce an infrared cutoff to the Coulomb potential, replacing it by the Yukawa potential V (r) = −e −r/r0 /r. After changing the variables of integration to p + ≡ 1 2 (p + p ′ ) and p − ≡ p − p ′ , Eq. (B3) has the form i 2 π d 3 p + d 3 p − e iτ p+·p− p −x g(p + , p − ) where g represents the slowly-varying part of the integrand. The integral (B13) is governed by the vicinity of p + = p − = 0, and we therefore Taylor-expand g around the origin. The zeroth order term, where g(p + , p − ) is replaced by g(0, 0), vanishes since the integrand is odd in (p + , p − ). Since the first two lines of the integrand in Eq. (B4) consist of two identical functions of p and p ′ , the derivative of g with respect to p − vanishes at the origin. The leading non-vanishing term in Eq. (B13) comes from the derivative of g with respect to p +x , giving: where g + ≡ ∂ p+x g(0, 0) = −ȧ(t n (t))/E(t n (t)). The range of the Yukawa potential, r 0 , should be set to a value that is large enough such that the wavefunction is contained within. Classically, the electron travels to distances of the order of E 0 /ω 2 , whereas τ is of the order of 1/ω. Therefore the modulus of expression in Eq. (B14) is of the order ofȧ(t n (t))E 0 .
We shall not delve into the analysis ofȧ(t n (t))E 0 any further, being content with the fact that it is much smaller than one, and is a smooth function of time that varies over a timescale of ω −1 . These two properties assure that frequencies of O(1) have amplitudes which are exponentially small in ω −1 , and thus can be neglected. Our statement, that the continuum-continuum termξ c (t) is negligible compared to the recombination termξ 1 (t) as far as HHG is concerned, is now established.
APPENDIX C: DERIVATION OF EQ. (28) The zeroth order SFA wavefunction is defined through the equation We now differentiate ξ 0 twice in time, using Eq. (C1). We obtain: where terms including derivatives of a(t) have been dropped, as they are negligible. Terms that do not con-tain |ϕ have been dropped as well, as they cannot contribute to HHG. If we substitute Eq. (14) into Eq. (C2), we obtain an expression forξ 0 (t) which is identical to Eq. (26), with the only difference that the matrix element 0|∂ x V (r)|p is replaced by Of the terms on the right hand side of Eq. (C3), we argue that the first one is the most dominant in the regime defined by Eq. (13a). In order to show this we note that 0|p has a pole at p 2 = −2I p . This follows simply from the asymptotic behavior of 0|r at r → ∞, which is an exponential falloff at the rate of 2I p . Therefore for p 2 ≫ 2I p , 0|p typically falls off like p 2 + 2I p to some negative power. This fact can be used to estimate the modulus of the second term on the right hand side of Eq. (C3) as |E(t) p 2 + 2I p 0|x|p | up to a factor of order one. Using |E| ≪ (2I p ) 3/2 from Eq. (13a) one arrives at the conclusion that the second term is negligible compared to the first one. Using similar argumentation and using Eq. (13a) one shows that the first term on the right hand side of Eq. (C3) is indeed the dominant one. Eq. (28) is thus established.