An analytical approach to high harmonic generation

With the recently introduced concept of dominant interaction Hamiltonians, we construct numerically as well as analytically the spectrum of high harmonics (HH) generated in electron–ion scattering under an intense laser field. This is achieved by switching the interaction along classical electron trajectories between the integrable cases of the dipole coupled laser field or the ionic potential, depending on which potential is stronger. As a side effect, the intrinsically chaotic character of the dynamics is mapped onto the potential switching sequence while the trajectories become regular. As a consequence, only a few per cent of the trajectories needed in full semiclassical calculations are sufficient to construct the HH spectrum.


Introduction
High harmonic generation (HHG) is one of the fundamental phenomena in intense laser matter interaction. Since the initial experimental observations [1,2] and theoretical explanations of it [3,4] (for an early review, see, e.g., [5]), much experimental work has been devoted to HHG with laser-driven atoms. More recently, also molecules [6,7] and clusters [8,9] have come into the focus of interest. This is partially due to the very simple yet accurate description that was given in [3,4] and that is nowadays called the three-step model. According to this model, the electron tunnels out of the combined nuclear plus laser potential and is then accelerated in the field. Upon recombination with the ion an HHG photon is generated with a maximum energy of 3.17U p + I p , where U p is the ponderomotive energy and I p is the ionization potential. This follows from a purely classical strong field consideration [3] and predicts the cutoff observed in HH spectra. Extended to include the ionization potential in the form of the phase exp(iI p t) [4], the HH spectrum with its plateau could be qualitatively reproduced.
The tunneling step is not decisive for the phenomenon and, upon using scattering initial conditions, the maximal energy is slightly reduced to 2U p + I p [10]. A semiclassical approach based on the Herman-Kluk (HK) propagator [11] has reproduced the quantum HH spectrum in the collision scenario very well and has revealed that the plateau formation is a quantum interference effect [10]. This implies that (a small) part of the continuum wavefunction has to recombine with the ion and assumes ground-state character in order for harmonic generation to occur.
While the continuum as well as the bound wavefunction can be approximated analytically in the spirit of the three-step model [4,12], the fraction of the continuum admixture could be 3 determined only numerically so far [12]. Here, we will calculate this fraction also analytically, with the help of the so-called dominant interaction Hamiltonians (DIH). We have introduced the DIH approach in [13] with the general motivation to facilitate the computational treatment of complicated dynamics. As an example, we successfully generated numerically HH spectra with DIH. In the following, we give a more complete account of this work and demonstrate, in addition, that with DIH the HH spectrum can be obtained analytically, as indicated above.
The paper is organized as follows. In sections 2 and 3, we present brief reviews of the theoretical approach to HHG as well as of the HK and the linearized semiclassical initial value representation (LSC-IVR) formalism. Section 4 explains in detail the implementation of the switching mechanism between dominant interactions for HH, where far away from the nucleus the laser potential is dominant (which is always assumed in the strong field approximation (SFA) [14][15][16]), while close to the nucleus the ionic binding potential is dominant. In this sense, the concept of DIH applied to HH dynamics is a natural extension of the SFA. It also completes the three-step model by providing analytically the time-dependent amplitudes in addition to the phases for the bound and continuum wavefunctions involved. In section 5, we compare accurately calculated quantum dipole accelerations and harmonic spectra with full semiclassical and DIH spectra, thereby demonstrating that the DIH spectra require considerably fewer trajectories compared to the full semiclassical approach yet retaining quantitative agreement. In section 6, we determine analytically the relative amplitudes of the Volkov wavefunction versus the ground state in the total wavefunction from the DIH. In section 7, we conclude and give an outlook on future developments.

High harmonic generation
As is well known, all essential facets of HHG are contained in a one-dimensional one-electron atom described by the so-called soft-core potential driven by a laser field, which in the length gauge is described by the Hamiltonian (all quantities are given in atomic units (au) if not stated otherwise) where E x f (t) cos(ωt) represents the interaction potential with a linearly polarized laser field in dipole approximation. We use here a pulse length of three and a half laser cycles with the parameters E = 0.1 and ω = 0.0378 and a step-like envelope function f (t) [10]. This implies a ponderomotive energy of In the soft-core potential [17] V sc ( the parameter a removes the singularity in the Coulomb potential, which causes numerical instabilities in the one-dimensional description. It has been chosen as a = 2, which yields the ground state energy E = −1/2 = −I p of the hydrogen atom. The central quantity behind theoretical calculations of the high harmonic spectrum is the dipole acceleration [18]. Using Ehrenfest's theorem twice, it reads gives the power spectrum which is proportional to the spectrum of emitted harmonics.

Quantum, classical and semiclassical approximations for the high harmonic spectrum
We will supply several approximate solutions (x, t) of the non-relativistic time-dependent Schrödinger equation to construct the high harmonic spectrum according to equation (4). These spectra will serve as a comparison with our new DIH approach whose numerical and analytical results will be presented in sections 5 and 6, respectively. Firstly, we generate the full quantum solution, obtained numerically with the split-operator fast Fourier transform (FFT) method. Secondly, we provide classical Wigner-type calculations, and thirdly, we have redone some of the semiclassical calculations presented in [10] for comparison. The classical and semiclassical calculations are based on the LSC-IVR developed by the Miller group [19][20][21] and the underlying concept will be briefly reviewed in the following.

The Herman-Kluk propagator
In the time domain a wavefunction (x, t) at time t is obtained by applying the quantum mechanical propagator with the time-ordering operatorT to an initial wavefunction (x, 0) according to A semiclassical approximation for the propagator based on multiple frozen Gaussians is based on the work by Heller [22]. It was justified later by Herman and Kluk [11], and was revitalized by the seminal work of Kay [23][24][25]. For a system with one degree of freedom whose coordinate is denoted by x, the HK propagator is given by Its main ingredients are normalized Gaussian wavepackets with a real and positive width parameter γ . The integration is performed over phase space points which serve as initial conditions of classical trajectories ( p t = p( p 0 , q 0 , t), q t = q( p 0 , q 0 , t)). Furthermore, classical mechanics comes into play via the classical action S = S( p 0 , q 0 , t) = t 0 L dt , with the Lagrangian L = T − V . The prefactor, devised by HK, is given by It ensures that the propagator is unitary in the stationary phase sense [26] and consists of elements of the so-called monodromy (or stability) matrix, This matrix describes the time evolution of small deviations in the initial conditions of a classical trajectory and can be obtained by solving linearized Hamilton equations. An even more general expression for the prefactor in equation (9) can be found in [27]. We note in passing that it has recently been shown that the HK propagator is the lowest order term in a series expression for the full quantum mechanical propagator [28,29]. The semiclassical approach is based solely on classical dynamical input. It allows, however, for interference effects due to the phase differences of the different trajectories in the sum of (7). This fact has been used in [10] to explain the plateau formation in HHG. It has been shown that 10 6 trajectories are needed to converge the results [30] and that a trajectory removal procedure [25,31,32] is not appropriate for the HHG problem.
If the HK propagator is to be applied to a Gaussian initial state according to equation (6), the x integration can be performed analytically and the integral over initial phase space contains a weighting function with a smooth decrease when leaving the wavepacket center (q α , p α ). This weighting function allows the application of a Monte Carlo procedure for the phase space integration [33].

The linearized semiclassical initial value representation
In order to compare our results to a purely classical dynamical approach to HHG, we now briefly sketch the LSC-IVR that was developed by the Miller group [19][20][21].
The LSC-IVR that we will employ here is based on the SC-IVR discussed in detail in [34]. It is used for both propagators in the correlation function expression After a linear expansion of the action difference in the exponent of equation (12), this leads to a purely classical expression with no quantum interference effects. The final result is given in terms of a single phase space integral, where A w is the Wigner transform of the operatorÂ, and an analogous relation for the operatorB. The method is therefore also referred to as the classical Wigner method. It has a long history of analysis and applications [35][36][37][38]. For our purposes,Â = | α (0) α (0)| is the operator representing the initial state, whereas for the dipole acceleration,B has to be chosen according toB = −|x x| ∂ V sc /∂ x. The LSC-IVR 6 expression for the dipole acceleration appearing in (4) is then given by The final expression contains quantum mechanical initial conditions due to the Gaussian nature of the initial state but is not plagued by an oscillatory integrand and thus is easy to determine numerically. However, interference effects are not contained in the LSC-IVR description.

The concept of dominant interaction Hamiltonians for high harmonics
This concept is based on the idea to simplify the dynamics of a system by partitioning the phase space into different regions where each region is defined by an interaction which is an approximation to the full interaction but which is dominant in this region. Ideally, the respective interactions lead even to integrable DIH. This is indeed the case for an atomic electron interacting with a laser field as described by equation (1). Far away from the nucleus the influence of the atomic potential V sc on the electron is negligible. In the vicinity of the nucleus it can occur that |V sc (x)| > |E x cos(ωt)|. In this case, the atomic potential V sc dominantly governs the dynamics such that the laser potential can be neglected. On the other hand, far away from the nucleus the influence of the atomic potential V sc on the electron is negligible. Both cases of dominant interaction lead to an integrable Hamiltonian, where the latter is for all times assumed by the SFA with the corresponding Volkov wavefunctions [16]. Clearly, the identification of a dominant interaction relies on local properties in phase space, and therefore our entire approach is based on a semiclassical information in which exclusively the parameters of the local classical dynamics enter. Hence, for DIH the underlying classical trajectories are propagated under the dominant interaction Hamiltonian until the specific trajectory reaches the phase space boundary where a switch occurs to another interaction which becomes dominant.
For the specific case of HH, we have devised the following switching scheme (see also figure 1) which was proven by a first test [13] to yield quantitative results.
If the trajectory starts out under the influence of the dominant laser potential V L , it will experience a switch to the soft-core potential V sc only if p(t c ) = 0 and |q (t c ) | x c = 3.0083 at a particular time t c . For |x| x c the laser potential will never become larger than the soft-core potential, where x c obviously depends on the properties of V sc as well as on the peak intensity of the laser light. In addition, the trajectory has to fulfil the condition that, immediately before the switching, the momentum points inward. We note in passing that the former conditions realize the recombination step in the three-step model, which assumes that only the accelerated electrons that return back to the nucleus can emit harmonics by recombining to the ground state [3,4,39]. Hence, the switching trajectories are responsible for HH in the DIH description and their initial conditions lie inside stripes in phase space, as worked out analytically in more detail in section 6.

Numerical results
Before showing the HH spectra using full quantum, full semiclassical (using the full Hamiltonian) as well as classical results (using LSC-IVR) and comparing them to the DIH results with the switching SFA Hamiltonian, we perform quite a stringent test for the DIH

Fidelity amplitude
A convenient quantity to compare both semiclassical approaches is the so-called fidelity amplitude, which is the overlap of the two wavefunctions, each propagated in a different potential Here,Û is calculated with the Hamiltonian given in (1) containing the full potential, whereasŨ uses DIH. Initially, the electron is described by the Gaussian wavepacket with q α = 70 au, p α = 0 au and the width parameter γ = 0.05 au.
In figure 2, the absolute square of the fidelity amplitude is plotted for a time scale of three and a half laser periods (after which the laser pulse has died out). Within the first period, the agreement between both wavefunctions is good, although temporarily strong decreases occur. Toward the end of the third period, the overall curve decreases to less than 0.2. Hence, we expect that the fine details in the HH spectrum requiring long-time propagation of the Fourierrelated dipole acceleration equation (3) will not be represented accurately. Nevertheless, will it be possible to reproduce the characteristics of the HH spectrum, i.e. the height and extension of the plateau as well as the overall shape with DIH, eventually even analytically? This question will be answered in the following.

Dipole acceleration and high harmonic spectrum
Using equation (3)  different approaches, i.e. with a full quantum propagation, a semiclassical propagation using the full Hamiltonian and DIH, the analytical DIH approach (see below), as well as a classical LSC-IVR time evolution in the full potential. For reference, also the Volkov result, derived from the analytical expression for the Volkov wavefunction [40], with q(t) = q α + p α t + E[cos(ωt) − 1]/ω 2 and p(t) = p α − E sin(ωt)/ω, which is based on the dynamics of the free electron in a laser field (i.e. no Coulomb field is present), is plotted. For the full quantum calculations, we have used the split operator FFT [41] and a grid of 2048 points extending from −300 to +300 au. As can be clearly observed, neither the classical (panel (e)) nor the Volkov result (panel (f) in figure 3) for the acceleration shows the detailed interference pattern, manifesting after the first laser period in the semiclassical and the quantum results. The agreement between all results within the first laser period is excellent, suggesting that the electron is prominently under the influence of the laser field. Furthermore, it is remarkable that although the full semiclassical calculations require 10 6 trajectories for convergence, the DIH semiclassical one does not need more than the LSC-IVR classical one, about 5 × 10 4 trajectories. Now we are in a position to look at the HH spectra, namely the Fourier transforms of the dipole accelerations, obtained with the different propagation schemes. As shown in figure 4, differences become visible due to the observed differences in the long-time behavior of the acceleration. In the classical as well as the Volkov case (e) and (f), the characteristic plateau is not visible due to the absence of interferences as already suggested by the acceleration in figure 3. To be more specific, in the classical case interferences are per se not included even though the full potential has been considered. In turn, the Volkov solution does not even take the soft-core potential into account, so that the laser-driven electron cannot recombine with the atom emitting a HH photon. In contrast, the quantum, the full and the DIH semiclassical results show the typical HHG spectrum (even harmonics are not strictly excluded due to the finite pulse length of the laser) with a plateau of the correct height and width, i.e. the correct cutoff frequency of 2U p + I p . Semiclassically, this can be attributed to the interference of phase-carrying free trajectories with trapped ones, i.e. trajectories being for some time mainly under the influence of the Coulomb potential oscillating around the nucleus [10]. Obviously, the DIH approach approximates this feature quite well.

Analytical wavefunction
Based on the success of the switching scheme in the numerical results shown above, we aim now at deriving an analytical expression for the total time-dependent wavefunction of the electron using the dominant interaction Hamiltonian idea. In the spirit of the stimulated recombination  effect, this wavefunction has to consist of a ground state and a Volkov part in order for harmonic generation to be possible.

Initial condition analysis
For the trajectories switching from the laser potential to the Coulomb potential, we require that at the switching instant t c the momentum p(t c ) is zero and we define q(t c ) = q c . The solution of the equations of motion for a free classical trajectory in a laser field with initial conditions q(0) = q 0 and p(0) = p 0 is given by Setting p(t c ) = 0 leads to switching times where n can take either odd or even values. By inserting this expression into the equation for q(t) we obtain for the initial position as a function of initial momentum By expanding arcsin (ωp 0 /E) and 1 − (ωp 0 /E) 2 in powers of p 0 , we obtain the series For the chosen parameters the prefactors of higher than linear order are small compared to the factor in front of the linear term. This observation suggests a linear approximation which becomes even better when n increases. Furthermore, due to the fact that the initial momenta for all the trajectories are distributed around small p α = 0 (see figure 5), we can neglect higher order terms in the expansion, which gives rise to This is the relationship that the initial conditions, p 0 and q 0 , have to fulfil in order for a trajectory to experience a switch at t c , which now simplifies to t c ≈ nπ/ω. Resolving (24) for the initial momentum leads to the linear relationship We stress that, for n even, q c lies in the interval [−x c , 0], whereas for n odd, q c lies in the interval [0, x c ]. For the different values of n this leads to origins of switching trajectories that are (narrow) stripes of width ω nπ x c in phase space which are depicted in figure 5 (see the inset). It is worthwhile to note that the initial conditions of the switching trajectories are close to those of the trapped and stranded ones that play a crucial role in the semiclassical explanation of the plateau [10]. A comparison of the complete respective trajectories is shown in figure 6. In the following, we will evaluate the integral for the case of n even. By using the following change of variables: p 0 → u = p 0 + ω nπ q 0 , the integral is transformed into The integral over q 0 can be performed analytically, yielding In order to evaluate this remaining integral, it is important to note that the integrand varies smoothly in the whole integration interval, so that we can obtain a very good approximation by considering it to be linear. Hence, the integral is well approximated by applying the mean value theorem of analysis, which gives the result Finally, we note that there is no need to implement the inward condition discussed in section 4, because the restriction of the initial conditions to a straight strip in phase space already cuts off the contributions of outgoing trajectories.

Comparison of the results
We now again compare full quantum results with wavefunction splitting results. In contrast to section 5, we calculate the switching wavefunction analytically according to with a time-dependent coefficient c(t) and a time-dependent normalization factor N (t) that can also be calculated analytically due to the fact that both the Volkov and the soft-core potential wavefunctions are of Gaussian type. In principle, we could have used the part of the Volkov wavefunction remaining after subtracting the contribution of the trajectories that perform a switch, but owing to the fact that only 5% of the trajectories do perform a switch, using the full Volkov solution is a reasonable approximation. In order to find an analytical expression for c(t), we need to consider the 2N possible instants, during the total propagation time, at which the trajectories experience a jump (once per half laser period and where N is the number of those half peiods). This can be managed by adding up the possible contributions c n calculated above to c(t) every half laser period, which can be expressed as  In order to understand what the overlap means in terms of trajectories, we should start by calculating the expression for the overlap between the Volkov and the ground state wavefunctions. The resulting expression is Clearly, we can see that the first exponential term vanishes as t → ∞, while the second one tends asymptotically to exp (−λp 2 t /4). For the Volkov dynamics the momentum as a function of time is given by p t = −E sin(ωt)/ω, from which we see that the momentum vanishes every half laser cycle and therefore we can see the maxima in the overlap every half laser period (see figure 8(d)).
Similar to the Volkov case, also in the analytical case displayed in figure 8(c), we observe peaks every half laser period. In addition, between these peaks we observe a step structure, which is related in a straightforward manner to the percentage of switched trajectories, namely the step function is given by c(t) 2 . For our full semiclassical overlap the essence is the same as the analytical. We still see the peaks every half laser period; however, between these peaks, now there is a more complex structure. In contrast to the overlap obtained with the analytical wavefunction, the structure between the peaks is no longer a step function but a more complex function. The difference comes from the fact that during a half laser period in the full semiclassical scheme several events can occur. This is reflected in the shape of the overlap; as soon as the laser field oscillation reaches it maximum, a number of trajectories start to leave the nucleus, resulting in a decay of the overlap during the first quarter laser period threafter. During the next quarter laser period other trajectories get trapped into the nucleus and this gives an increment in the overlap. Although the structure obtained by the analytical and the full semiclassical wavefunctions is not the same, it is noteworthy to observe the qualitative agreement between the number of switched trajectories and the number of trapped trajectories.
A question remains open: why can we not observe an increment in the number of trapped trajectories between 2.5 and 3.0 laser periods as predicted by the quantum result? We conjecture that the presence of tunneling in the quantum scheme is responsible for the disagreement with the semiclassical result. In principle, tunneling effects can be incorporated also semiclassically by taking into account higher order terms in the series expression of the propagator [28,29], which is beyond the scope of this presentation, however.
With the final result for c(t) plugged into the approximate wavefunction (30), we are now able to calculate an 'analytical' result for the dipole acceleration. In panel (d) of figure 3, it is compared with the full quantum results. The corresponding spectra are shown in panel (d) of figure 4 and, for both quantities, we observe excellent agreement.

Conclusions and outlook
We have shown that the concept of DIH together with the semiclassical propagation of the wavefunction according to the HK propagator reproduces the crucial features of the HH spectrum of a laser-driven atomic system. This is remarkable due to the fact that the much simpler classical dynamics allows for a calculation with two orders of magnitude fewer trajectories. Furthermore, the potential switching scheme we have devised to realize DIH for laser-driven atomic electrons leads to an analytical time-dependent wavefunction, where the crucial missing part was the determination of the relative amplitude between the ground state and the Volkov part of the total wavefunction that has previously been chosen as a constant, leading to good numerical agreement with the full quantum result [12].
For future calculations of larger systems (molecules, clusters), the absence of laser-induced chaos suggests that the DIH approach is a realistic possibility to obtain numerical results if appropriate switching schemes can be identified. Classical calculations show no plateau at all and full quantum solutions are at present not available for systems with more than five to six degrees of freedom. Finally, fundamental questions of electron-electron correlation appearing, e.g., in non-sequential double ionization of helium [42] could be tackled in the spirit of DIH.