Quantum and semiclassical phase-space dynamics of a wave packet in strong fields using initial-value representations

We assess the suitability of quantum and semiclassical initial value representations, exemplified by the coupled coherent states (CCS) method and the Herman Kluk (HK) propagator, respectively, for modeling the dynamics of an electronic wave packet in a strong laser field, if this wave packet is initially bound. Using Wigner quasiprobability distributions and ensembles of classical trajectories, we identify signatures of over-the-barrier and tunnel ionization in phase space for static and time-dependent fields and the relevant sets of phase-space trajectories in order to model such features. Overall, we find good agreement with the full solution of the time-dependent Schr\"odinger equation (TDSE) for Wigner distributions constructed with both initial-value representations. Our results indicate that the HK propagator does not fully account for tunneling and over-the-barrier reflections. However, it is able to partly reproduce features associated with the wave packet crossing classically forbidden regions, although the trajectories employed in its construction always obey classical phase-space constraints. We also show that the Coupled Coherent States (CCS) method represents a fully quantum initial value representation and accurately reproduces the results of a standard TDSE solver. Furthermore, we sow that both the HK propagator and the CCS approach may be successfully employed to compute the time-dependent dipole acceleration and high-harmonic spectra. Nevertheless, the semiclassical propagator exhibits a worse agreement with the TDSE than the outcome of the CCS method, as it neither fully accounts for tunneling nor for over-the-barrier reflections. This leads to a dephasing in the time-dependent wave function which becomes more pronounced for longer times.


Introduction
Initial-value representations (IVRs) such as the coupled coherent states method [1] and the Herman Kluk propagator [2] are widely used in many areas of science. These approaches allow an intuitive interpretation of a time-dependent wave packet in terms of trajectories in phase space, and account for binding potentials, external fields and quantum-interference effects. Furthermore, the numerical effort in IVRs does not scale exponentially with the degrees of freedom involved. This efficiency may be increased by employing several strategies, such as dominant Hamiltonians in specific phase-space regions [3,4], or quantum-state reprojection [5,6].
In principle, all these features make initial-value representations very attractive to strong-field and attosecond science. In fact, it is well known since two decades that strong-field phenomena such as high-order harmonic generation (HHG), abovethreshold ionization (ATI) or nonsequential double ionization (NSDI), may be described as the result of the laser-induced scattering or recombination of an electron with its parent ion [7]. This implies that electron orbits play a very important role in the understanding of these phenomena. This has led not only to a myriad of applications, such as the attosecond imaging of dynamic processes in matter (see, e.g., [8,9,10]), but to the extensive use of orbit-based approaches (see, e.g., [11]).
In particular, classical and semiclassical methods are very popular. Classically, ensembles of electrons that behave according to the above-mentioned recollision picture are constructed in order to mimic the behavior of the quantum mechanical wave packet and both the external laser field and the binding potentials are fully incorporated. This is the key idea behind classical-trajectory methods, which have reproduced key features such as the low-energy structure in ATI [12] and the Vshaped structure observed in NSDI [13,14]. These methods, however, cannot account for quantum interference effects, or tunnel ionization. Most quantum-mechanical and semi-classical methods in strong-field physics, on the other hand, incorporate such effects, but make drastic approximations on the residual binding potential. A typical example is the strong-field approximation (SFA), which, in conjunction with the steepest descent method, is the most used approach in this field. In the SFA, the continuum is approximated by field-dressed plane waves, i.e., the potential is not accounted for in the electron propagation. Nonetheless, the interplay between the external laser field and the binding potential is important and has revealed itself in many ways. For instance, this interplay leads to a prominent low-frequency structure [15,12,16,17] and fan-shaped interference patterns [18,19] in ATI, and strongly influences NSDI in circularly polarized fields [20,21]. Because this interplay is important, in the past few years, Coulomb-corrected analytic approaches have been developed and successfully applied to strong-field phenomena [22,16,23,24]. These approaches, however, require the external field to be dominant.
On the other hand, initial-value representations have only been employed in strong-field physics in relatively few publications. Specifically, in [25,3,4] HHG spectra have been computed using the Herman Kluk (HK) propagator, and in [26,27,28] NSDI ion and electron momentum distributions have been calculated employing the Coupled Coherent States (CCS) method. In order to be able to apply these methods widely, major challenges must be overcome. Concretely, the trajectories employed in the construction of the time-dependent wave function with the HK propagator are real, i.e., they cannot cross classically forbidden regions in phase space. Hence, tunnel ionization may not be properly accounted for. Since the 1990s, there has been considerable debate about whether one may model tunneling employing semiclassical IVRs such as the HK propagator, and if so, to which extent (see, e.g., [29,30,31,32]). In order to circumvent this problem, in [25,3,4] the initial electronic wave packet has been placed far away from the core. Unfortunately, these initial conditions leave out many strong-field problems, for which tunneling is expected to be the dominant ionization mechanism. On the other hand, tunneling is present in the CCS method, because it is a basis-set method. However, special effort must be made to choose a trajectory-guided basis which is suitable for tunneling.
Nevertheless, one may wonder whether the over-the-barrier dynamics, per se, would not be sufficient for the modeling of strong field wave-packet dynamics in the presence of the Coulomb potential. This is a legitimate question, especially if one considers that classical models, for which tunnel ionization does not occur, have been hugely successful in reproducing a number of features in ATI and NSDI. In some of these methods, tunnel ionization has been mimicked by employing the quasi-static Ammosov-Delone-Krainov (ADK) tunneling rate, which may explain this success. However, there exist also purely classical models in which the electron ensemble is left to propagate without the need for any ad-hoc quantum mechanical ingredient. For a discussion of these models in the NSDI context see our review article [33].
For that reason, in this work we perform a systematic analysis of the dynamics of an electronic wave packet in a strong field, with particular focus on ionization, and on what is left out by initial value representations. This analysis will be performed in phase space, for reduced dimensionality models, under the assumption that the electronic wave packet is initially bound. As a benchmark, we will employ the full solution of the time-dependent Schrödinger equation (TDSE).
This article is organized as follows. In Sec. 2, we provide the necessary theoretical background in order to understand our results. Subsequently, in Sec. 3, we will have a closer look at how the time-dependent wave packet overcomes the potential barrier resulting from the combined action of the external field and the binding potential. For that purpose, we will employ quasiprobability distributions in phase space, for a static and a time dependent field, and for long-and short-range potentials. As, for short-range potentials, tunneling is expected to be more prevalent than over-thebarrier ionization, in this case we will also perform a comparison with the CCS method. In Sec. 4, we present approximate estimates, in which the transmission through the potential barrier is computed analytically using uniform WKB approximations, and a comparison with an inverted harmonic oscillator is made. In Sec. 5, we will show how the phase differences between the semiclassical and quantum mechanical computations lead to discrepancies in HHG spectra and the time-dependent dipole acceleration. In Sec. 6, we will provide a summary of the main results and the conclusions to be drawn from this work. Finally, in the Appendix we discuss the fact that the HK and the CCS methods share a common origin.

Model
For the sake of simplicity, we consider a one-electron, one-dimensional atom. The Hamiltonian of such a system can be expressed aŝ whereV a andV E represents the binding potential and the interaction with the laser field in the length gauge, respectively, and the hats denote operators. Atomic units are used throughout. The external field exhibits the same temporal dependence as in [25,3], i.e., with E and ω the intensity and frequency, respectively, of the laser field. The binding potential is taken either as the soft-core potential with λ = 1 constant, or as the short-range Gaussian potential where λ = 1/2. As a benchmark, we employ the full solution of the time-dependent Schrödinger equation which is solved in coordinate space for Ψ(x, t) = x|Ψ(t) using a split operator method.
The initial wave function is defined as the Gaussian wave packet centered at p α , q α , with q α = 0. Furthermore, the parameter γ defines the width of the wave packet. This initial choice facilitates the implementation of initial-value representations. The initial energy expectation value, computed for the system in the absence of the external driving field, is given as If the wavepacket is initially centered at q α = 0, this expectation value reads as and for the Gaussian and the soft-core potentials, respectively, where K n (z) is the modified Bessel function of the second kind. Eq. (7) depends on the width and the initial momentum of the initial wave packet, and on the potential parameter λ. The width and the position of the initial wave packet were chosen to minimize its initial energy. This procedure yields (p α , q α ) = (0, 0) and γ ≃ 0.46 and (p α , q α ) = (0, 0) and γ ≃ 0.65 for the soft-core and the Gaussian potentials, respectively. The groundstate energies associated with the potentials V sc (x) and V G (x) are E sc = −0.67 a.u. and E G = −0.594, respectively.

Quantum and semiclassical initial value representations
Although quantum mechanics is usually in coordinate or momentum space, it can also be formulated in phase space. In this section we present the initial-value representations in phase space employed in this work. These representations describe the system dynamics in phase space, and are the Coupled Coherent States (CCS) method, which is a formally exact approach for describing quantum mechanics, and the Herman-Kluk (HK) propagator, which is semiclassical. Both methods can be derived from the same source, and this derivation is provided in the Appendix. More details can be found in the review article [1].
The CCS method represents a time-dependent wave function as a superposition of time-dependent, nonorthogonal Gaussian coherent states (CS) |z = |z(t) guided along the trajectory determined by the Hamiltonian averaged within such a basis. Explicitly, where the coherent state |z is labeled by a single complex number which is defined in terms of the position q = q(t) and the momentum p = p(t) of the particle. The expression denotes the classical action along the trajectory defined with regard to the matrix element H ord (z * , z) = z|Ĥ ord (â † ,â)|z . This represents the diagonal elements of the ordered Hamiltonian matrixĤ ord (â † ,â). In general, where |z , |z ′ denotes two arbitrary coherent states. In coordinate space, is a Gaussian wave packet centered at the phase-space coordinates q and p. For the Herman Kluk propagator, the time-dependent wave function is once more expressed in terms of Gaussians in phase space (for seminal papers see [2,34,35,36]). Following the original work [2], they are usually labeled not with a single complex number z but with two real numbers p and q. In the p, q-notation the phase of CS in (16) differs from that of z-notations by ipq/2 which is compensated by the different form of action (18) used in [2] as opposed to that of Eq. (12). This specific formulation is considered in our implementation and in the results that follow. Explicitly, where |q, p represents a coherent state whose expression in coordinate representation is given by which differs from Eq. (14) employed in the CCS approach by an insignificant phase factor, and is given in terms of the elements m uv = ∂u/∂v 0 of the monodromy matrix. In Eq. (15), the action reads as where H cl is the classical Hamiltonian, in which the operatorsx,p in (1) have been replaced by the phase-space variables q, p. For the initial wave packet (6) considered here, The integral is carried out over phase space coordinates which are used as initial conditions of the classical solutions (q, p). Apart from trivial notations, the HK propagator and the CCS method differ in three main ways: (1) The trajectories of the Gaussian Coherent States in the HK method are purely classical, while in CCS the trajectory of a CS |z(t) is driven by the Hamiltonian z|Ĥ|z . This latter Hamiltonian is the average of the quantum Hamiltonian with regard to the coherent states, and takes into account the local zero-point energy and further corrections due to commutators. This makes all wells more shallow and lowers all potential barriers. The CCS and HK trajectories are identical only in the case of harmonic potential. (2) The CCS representation is a formally exact basis set technique. It uses coupled quantum equations for the coefficients D z (t), obtained simply by substitution of (10) into the Schrödinger equation, to propagate the wave function (10). In contrast, in the HK theory the coefficients are obtained by an analytical semiclassical formula, which includes the elements of the monodromy or stability matrix. These expressions result from the so called local quadratic approximation, which only takes into account the first and second term in the Taylor expansion of the potential energy around a specific trajectory, and assume that only the coupling of coherent states near this trajectory is important. Physically, this implies that, while in the CS representation the trajectories are coupled through the amplitudes D z (t), each trajectory in the swarm employed in the HK method contributes independently. Indeed, for each trajectory there is one prefactor, which depends only on the information carried out by that specific trajectory. (3) On a more technical level, the CCS is often used in conjunction with various algorithms of basis set expansion, which generate additional basis functions and reproject the wave function on the new basis set. These adoptive basis sets allow to follow complicated features of the dynamics more efficiently. Although reprojection has been used in conjunction with the HK method as well, in the latter case it is less common [5,6].
In the Appendix a sketch is presented of how both CCS coupled equations and the HK formula are obtained from the same source, which is the integro differential form of the Schrödinger equation in the continuous CS representation, closely following the review article [1]. For a detailed account of the similarities and differences between both initial-value representations, see also [37,38,39].

Phase-space dynamics
In the following, we wish to analyze different ionization mechanisms in phase space. According to the quasi-static tunneling picture, a low-frequency, time-dependent field and the binding potential determine an effective potential barrier whose maximum will be given by V eff (x s , t) at the coordinate x s such that ∂V eff (x, t)/∂x| x=xs = 0. If the total electron energy is larger than this value, it may escape via over-the-barrier ionization. If the total energy is smaller, tunneling is expected to be the dominant ionization mechanism.

Phase portraits and classical-trajectory analysis
With this aim in mind, we will perform a phase space analysis of both the trajectory ensemble used to construct the semiclassical wave function in Eq. (15), and the wave functions Ψ HK (x, t) and Ψ(x, t). For simplicity, we will first study what happens if the driving field is static. This is a good approximation for the instantaneous barrier (20) if the frequency is low enough. In this case, the time-dependent field E(t) is replaced by E in Eq. (2). In phase space, the classical static Hamiltonian then reads as In these studies, we will consider the soft-core potential (3). In Fig. 1, we present the phase portrait of the system for the Hamiltonian (21) and static fields of different amplitudes. The figure shows that the condition p 2 /2 + V a (q) + Eq = E min , where E min is the minimal energy necessary for the electron to undergo over-the-barrier ionization defines a separatrix, which crosses at (q, p) = (q s , 0). For energies E below the separatrix energy the dynamics can be either bounded or unbounded depending on whether the spatial coordinate of the orbit is larger or smaller than that of the saddle point (q s , 0). This means that an orbit with E < E min will remain unbounded if its initial spatial coordinate q 0 lies on the left-hand side of the saddle point, ant that it will remain bound if q 0 lies on the right-hand side of q s (Figure 1). On the other hand, orbits with E > E min will remain unbounded regardless of the initial value of their spatial coordinate. Furthermore, the trajectories in the ensemble always respect the constraints dictated by classical dynamics. This is explicitly shown by the thin lines in the figure, which illustrate the time evolution of some sample trajectories. In fact, the trajectories whose energy lie below E min always remain bound and propagate along closed orbits bounded by the separatrix. Those that follow the separatrix from below (E min < E < 0.3 a.u.) go around the bound region in phase space, but never cross this region.
Classically, if a trajectory has a specific energy it will occupy a well-defined phasespace orbit. Quantum mechanically, however, an initial wave packet of a specific  energy may occupy many regions in phase space. In fact, due to the uncertainty relation, a strong spatial localization will lead to a larger momentum spread and viceversa. This is shown in the upper panels of Fig. 2, where we display the phase-space representation (19) of two Gaussian wave packets (6) centered at (q α , p α ) = (0, 0) with different widths γ = 0.5 and γ = 0.05. These widths give bound-state energies of roughly E = −0.5 a.u. or E = −0.67 a.u., respectively. If now a set of initial conditions in phase space is generated to match these distributions, these conditions will spread over several regions, bound and unbound, as shown in the middle panels of Fig. 2. These are the starting points (q 0 , p 0 )) for an ensemble of classical trajectories, which will be propagated in time and used in the construction of the semiclassical wave function Ψ HK (x, t).
After some time has elapsed [lower panels in Fig. 2], the distributions follow the constant-energy curves in phase space. For the wave packet with γ = 0.5, the distributions are concentrated in the bound region or, for q < q s , in the momentum region below the separatrix. The region on the left-hand side of the saddle around the axis p = 0 is practically unpopulated. In contrast, for the wider wave-packet in position space (γ = 0.05), there are phase-space events in this region, but closely following the separatrix from above.
These features are determined by the initial position and momentum spreads of the wavepacket and the corresponding trajectory ensemble. For γ = 0.5, the wavepacket is more localized in position space, so that the associated classical ensemble practically does not occupy the region on the left-hand side of the saddle. The trajectories that lie between the separatrix and the curve associated with the energy E = −0.3 a.u. then follow the separatrix from below. The remaining trajectories, for E < E min or E > −0.3 a.u. either remain bound or lead to the distribution below the latter curve, for high negative momenta. In contrast, for a larger momentum spread (γ = 0.05) the region on the left-hand side of the saddle is reasonably populated from the start.
According to this analysis, only if the above-stated constraints vary in time may a classical bound trajectory become unbound or vice-versa. This is what happens if  The gradient color shows the time variation from t = 0 (blue) to t = 30 (red), which spans slightly less than the first quarter cycle of the the field. The dot depicts the initial condition of an initially unbound trajectory which becomes trapped because of the time dependent field, and the remaining symbols illustrate the phase-space coordinates (pt, qt) at specific times t. For t > T /4, the region bounded by the separatrix starts to decrease, until the trapped trajectory is eventually able to escape.
the external field is time dependent. In Fig. 3, we illustrate this behavior for the laser field associated with Eq. (2). In this case, the bound region first increases in time, so that an initially unbound trajectory may become trapped. After a field crossing, the field amplitude starts to increase, and, consequently, the bound region will shrink. This will lead to a bound trajectory becoming unbound. Hence, the electron will be able to escape. According to the classical constraints, however, its momentum must be such that it only may go over the barrier.

Wigner quasiprobabilities
We will now employ the Wigner quasiprobability distribution (also known as the Wigner function) in order to relate the above picture to the wave-packet propagation in phase space [40]. This function is defined as If Eq. (22) is integrated over the momentum or position space, the corresponding probability densities are recovered. One should note, however, that the Wigner distribution may exhibit positive as well as negative values. Hence, strictly speaking, it cannot be associated with a probability density. Nonetheless, it does give an intuitive picture of the wave-packet dynamics in phase space, and provide valuable information about momentum-position correlation. Wigner distributions are widely employed in quantum optics, and have also been used in strong-field physics in order to access the dynamics of ionization [41], rescattering [42,43], double ionization [44] and HHG [45].
A common feature identified in [41] and [45] was a tail in the Wigner distribution,  which has been associated with tunnel ionization. In the following, we will plot the square of the Wigner function, W 2 (q, p) as it makes the above-mentioned tail slightly clearer.
In the left and right panels of Fig. 4, we depict the squares of the Wigner distributions obtained using the TDSE and the HK propagator, respectively, for a static field and the soft-core potential (3). This quasiprobability has evolved from an initial Gaussian state centered at (q α , p α ) = (0, 0) and width γ = 0.5 (see Fig. 2 for details). The time propagation of the wave packet, and in particular the shape of the Wigner function, are strongly influenced by the separatrix. Throughout, the figure shows a distinct tail in the Wigner functions leaving the bound phase-space region. For short times, this tail follows the separatrix from below, as shown in the upper panels of the figure. This strongly suggests that the continuum is reached by over-thebarrier ionization: the electronic wave packet does not leave the core with vanishing momentum, but, rather, with the minimum necessary momentum to overcome V eff (q s ). As time progresses, interference fringes start to build up on the left-hand side just above the separatrix. These fringes have been identified in [41], for a delta-potential model in a static field, and in [42,45] for long-and short-range potentials in timedependent fields, and have been associated with the quantum interference of ionization processes occurring at different times. Apart from that, there is a pronounced tail now following the separatrix from above. This implies that parts of the wave packet are reaching the continuum with lower energy than that required for over-the-barrier ionization to occur. In other words, tunnel ionization may be taking place, both for the TDSE and the HK wave packets. Near q = q s , the momentum at this tail is even approximately vanishing, in agreement with the model in [41]. A decrease in the momentum associated with the tail is intuitively expected as, physically, the components of the wave packet with lower energies will take longer to reach the barrier [31].  Figure 5. Square of the Wigner function calculated using the HK propagator leaving out the trajectories starting outside the bound region, for the same parameters as in Fig. 4 and times t = 10 a.u. and t = 20 a.u. (left and right panels, respectively).
All the above-mentioned features are present both for the semiclassical and quantum mechanical computations. In fact, the agreement between the outcomes of both approaches is quite striking. The presence of interference patterns in both cases is not surprising, as it is well known that the HK propagator is capable of reproducing such effects. It seems, however, that the semiclassical propagator allows the wavepacket to cross classically forbidden regions. In other words, classical trajectories corresponding to over-the-barrier energy appear to mimic transmission through a barrier for E < E min . This is not obvious, as these trajectories do not violate any phase-space constraints, or cross classically forbidden regions.
We have however verified that the trajectories whose initial coordinates (q 0 , p 0 ) lie outside the bound phase-space region lead to the tail in the Wigner distributions. For clarity, in Fig. 5 we show the Wigner function computed from the HK propagator leaving out these trajectories. In the figure, both the above-mentioned tail and the interference fringes are absent. These findings are in agreement with the results in [46], in which a nonlocal behavior of the Wigner function around separatrices has been observed for an inverted harmonic oscillator, and with those in [47], which find that the transmission coefficient associated with a parabolic barrier is related to the quantum-mechanical weight of all classical trajectories with enough energy to go over the barrier. The disappearance of the fringes in the Wigner functions around p = 0 for q < q s also support the assumption that the Wigner function is non-local.
For the time-dependent field in Eq. (2), we observe that the tails with momenta following the separatrix from below prevail in the Wigner function. Furthermore, there are many more fringes associated with quantum interference and rescattering events (for discussions of these fringes see [42,43]). Snapshots of W 2 (q, p) in the time dependent case are presented in Fig. 6. The separatrix, which is now time dependent, is also displayed in the figure as the thick black lines. This suggests that, for this type of potential, there will be enough over-the-barrier dynamics for the approximate modeling of strong-field ionization and rescattering. We see, however, that the agreement between the HK propagator computation and the TDSE worsens with time. This is related to the fact that non-classical effects such as tunnel ionization and over-the-barrier reflections cannot be fully accounted for, and become dominant for longer times [31].

Comparison with the CCS method
We will now investigate the Gaussian potential (4), and an initial wave packet with γ = 1. For a short-range potential, the effective potential barrier is much steeper and there is no Rydberg series. Hence, tunneling is expected to be dominant for the parameter range in question. Furthermore, we will also perform a direct comparison with the CCS method, which is well known to incorporate tunneling.
In Fig. 7, we display snapshots of W 2 (q, p) computed with the TDSE, the HK propagator, and the CCS method (left, middle and right columns, respectively) using a static field. For t = 10 a.u. (upper panels) there is a very good agreement between the outcome of all approaches, with a distinct tail in the quasiprobability distributions following the separatrix from above. This behavior holds even for the results obtained with the semiclassical propagator, and is different from that observed for the softcore potential (see Fig. 4). Physically, this suggests that there will be a nonvanishing probability density leaving the core with |p| < |p min |. For t = 20, a longer tail extending to far beyond the core region and interference fringes are present in all cases. For these longer times, the agreement between the HK propagator computation and the TDSE worsens. For the standard CCS, this is also the case, as, in this region, anharmonicities associated with the binding potential become more important. In order to overcome this problem, it is necessary to perform a periodic reprojection of the trajectories onto a static grid. A detailed discussion of this method in the strong-field context will be given elsewhere [48] (see also Refs. [5,6] for details).
We have also verified that, while the trajectories employed in the HK propagator never cross the separatrix, those used in the CCS method do. This is a consequence of the fact that the ordered Hamiltonian H ord (z * , z) is different from the classical Hamiltonian H cl (p, q). In fact, as a function of the phase-space coordinates (q, p), for a static field H ord (z * , z) is given by with η = (λγ/(γ + λ)). If the separatrix is however defined with regard to H ord (p, q) instead of the classical Hamiltonian, the trajectories defined by the CCS method will not cross. In other words, averaging the potential over a Gaussian CS basis leads to the lowering of the barrier, which partially takes tunneling into account. A direct comparison of H st ord (p, q) and H st cl (p, q) shows an effective energy shift γ/4 and a shift in the binding potential, where V G (q) is given by Eq. (4). For discussions of these shifts see [38,39].

Approximate estimates
It is a well known fact that, for an inverted harmonic oscillator, IVRs lead to exact descriptions of the time-dependent wave packet dynamics and Wigner quasiprobability distributions exhibit nonlocal behavior [46]. In [31], however, it has been argued that this does not hold for a general barrier, unless transmission occurs close enough to its top. In this case, the potential barrier may be approximated by an inverted harmonic oscillator, and IVRs give reasonable, albeit not exact, results. It is thus our objective to assess whether, for the potential barriers employed here, the results obtained in the previous section may be justified in this way.
For that reason, we expand V eff (x) around the saddle x s , and, using the uniform WKB approximation, we compute the transmission coefficient through this barrier. This coefficient reads as where x l and x r are the left and right turning points, respectively, for which V eff (x l ) = V eff (x r ) = E. In Fig. 8, these results are displayed as a function of the energy of the initial wave packet. The figure confirms that the inverted harmonic oscillator is a reasonable approximation, for the parameter range employed in this work. This approximation becomes increasingly more accurate as the energy approaches the threshold and over-the-barrier regime. A better agreement is observed for the soft-core potential. This is consistent with the results in the previous section (see Figs. 4 and 7)

High-harmonic generation
We will now establish a connection with previous work in the literature, and compute high-harmonic spectra with the HK propagator. We will use the soft-core potential (3) and the interaction Hamiltonian (2), but assume that the electronic wave packet is initially localized at the core. The HHG spectra are given by where denotes the dipole acceleration. The time-dependent wave function is either given by the solution of Eq. (5) in coordinate space or by Ψ HK (x, t). In Fig. 9, we display d(t), together with its power spectrum, for an initial wave packet (6) centered at q α = 0 and with vanishing momentum p α = 0, whose width is γ = 0.5 a.u. The figure shows a good qualitative agreement between the fully quantum and semiclassical acceleration, which roughly follow the field and exhibit a series of  Figure 9. Dipole acceleration over two field cycles, together with the HHG spectra (left and right panels, respectively), computed using a laser field of frequency ω = 0.05 a.u. and amplitude E = 0.075, the soft-core potential (3) and an initial wavepacket (6) of width γ = 0.5, centered at qα = 0 and with vanishing momentum pα = 0. Panels (a), (b) and (c) were computed using the full TDSE computation, the Herman Kluk propagator, and the CCS method, respectively. The dashed lines show the energy position of the cutoff, which in this case is located at Ip + 3.17Up = 49ω.
high-frequency oscillations. These oscillations, together with spatial localization, are responsible for the HHG plateau. They have been identified and discussed in previous publications employing the TDSE [49] and the HK propagator [25,3,4], for an initial wave packet far from the core. They have also been studied in a different context, namely the adiabatic approximation [50] and Bohmian trajectories [51,52]. The agreement between the outcome of the TDSE and the HK propagator is particularly good for times below half a cycle. Phase differences however arise between 0.5T and 1.5T , with T = 2π/ω 0 . This good agreement persists for the spectra, which exhibit a long plateau followed by a sharp cutoff at I p + 3.17U p and a reasonably similar substructure, such as the intensity modulations near the cutoff (harmonic orders 45 ≤ N ≤ 59) and in the below-threshold harmonic region (harmonic orders N ≤ 19). Discrepancies, however, exist in the overall intensity of the plateau, which is slightly higher for the semiclassical spectrum, and in this substructure. These discrepancies are associated with the phase differences mentioned above.
Physically, this dephasing is a consequence of the fact that the semiclassical IVRs do not fully account for processes in which classically forbidden regions in phase space are crossed, which affect the overall phase of the wave packet. Examples of such processes are tunneling ionization and over-the-barrier reflections. We have indeed found that this dephasing decreases if the energy of the initial wavepacket is increased, for instance, by changing its width or initial momentum p α . This, together with an overall decrease in the plateau height, is consistent with the fact that, for p α = 0 or γ ≪ 1, there will be an enhancement in the over-the barrier pathways. A similar dephasing has also been observed in the context of the transmission of a wave packet through an Eckart barrier [30]. These issues may constitute a problem for long times [31]. We have observed that our results are reasonably accurate for t ≤ 3.5T . In contrast, the results from the CCS method are once more practically identical to those of the TDSE.

Conclusions
The results presented in this paper strongly suggest that semiclassical initial-value representations, a concrete example of which is the Herman Kluk propagator, may be employed for describing strong-field wave-packet dynamics, even if this wave packet is initially bound and located within the core region. We have observed a reasonably good agreement with fully quantum mechanical methods, such as the numerical solution of the time-dependent Schrödinger equation (TDSE), or the Coupled Coherent States (CCS) representation, for static and time-dependent fields, and different types of binding potentials. This agreement manifests itself in the time evolution of Wigner quasiprobability distributions, and in the computation of time-dependent quantities such as the dipole acceleration. A noteworthy feature is the presence of tails in the Wigner distributions, which leave the core region following the separatrices very closely. Depending on the momenta associated with this tail, it may be related to over-the-barrier or tunneling ionization. Similar tails have been identified in the literature, using either a zerorange potential [41,43] or the TDSE [45]. Therein, however, focus has been placed on how this tail behaves outside the core region, and on its agreement with classical trajectories as defined by the strong-field recollision model [7]. In this work, we place more emphasis on the position-momentum correlation in the vicinity of the core as evidence for different ionization mechanisms. For the soft-core, long-range potential employed here, the momentum in this tail indicates substantial over-thebarrier ionization, while for the Gaussian, short-range potential, tunnel ionization seems to be dominant. It therefore appears that the Wigner function is crossing a classically forbidden region.
Thus, one may argue that, while the trajectories in the semiclassical method used in this paper are classical and will never cross a separatrix, the initial positionmomentum spread will provide the wave function with access to classically forbidden regions. This probability density, and hence the Wigner function, may exhibit nonlocal behavior around a separatrix. This specific behavior has been shown in [46] for an inverted harmonic oscillator.
On the other hand, if the classical trajectories that start outside the bound phasespace region are removed, the tail in the Wigner functions disappears. This implies that they are an essential, and fully classical ingredient for reproducing this tail in the context of a semiclassical initial-value representation. Furthermore, a semiclassical approximation would only be exact for a parabolic barrier, while for the potentials employed here it will be only approximate. The presence of anharmonicity in more realistic barriers, together with the existence of tunneling loop structures in phase space for anharmonic potentials, was in fact employed in [31,53] as a criticism to the findings in [46]. These structures become dominant for long times. Nevertheless, any barrier in the vicinity of its maximum may be approximated by an inverted harmonic oscillator. Specifically, our approximate estimates for the transmission coefficient show that this is a reasonable approximation for the parameter range employed in this work. Another noteworthy aspect is that, in [31,53], the potential barrier is flat, i.e., lim x→∞ V a (x) = 0, while, for the effective potential V eff (x) employed in this work, this condition does not hold. In such references, it was repeatedly emphasized that this condition led to the tunneling contributions becoming dominant for longer times.
One should bear in mind, however, that the trajectories do need to cross classically forbidden regions for the phase of the wave function to build up correctly. If this does not occur, there will be a degradation of this phase for longer times [31,53]. Since the dipole acceleration and the HHG spectra are strongly dependent on this phase, only a qualitative agreement with the full quantum mechanical result may be reached if a standard semiclassical IVR is employed. A quantitative agreement would require more sophisticated approaches, possibly along the lines in [32,54] or by adding higher order correction terms to the HK propagator as in [55,56]. Still, the results in Sec. 5 show that the HK propagator may be quite useful in modeling strong-field phenomena and understanding quantum-interference effects in, for instance, few-cycle laser pulses, for which this critical regime has not been reached.