New Journal of Physics

The open access journal at the forefront of physics

Deutsche PhysikalischeGeseUschaft DPG IOP Institute Of PhySjCS

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

C Zagoya1, J Wu1, M Ronto2, D V Shalashilin2 and C Figueira de Morisson Faria1

1 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

2 School of Chemistry, University of Leeds, Leeds LS2 9JT, UK E-mail: c.faria@ucl.ac.uk

Received 12 May 2014, revised 6 August 2014 Accepted for publication 9 September 2014 Published 28 October 2014

New Journal of Physics 16 (2014) 103040

doi:10.1088/1367-2630/16/10/103040

Abstract

We assess the suitability of quantum and semiclassical initial-value representations (IVRs), 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 phasespace trajectories to model such features. Overall, we find good agreement with the full solution of the time-dependent Schrödinger equation (TDSE) for Wigner distributions constructed with both IVRs. Our results indicate that the HK propagator does not fully account for tunneling and over-the-barrier reflections. This leads to a dephasing in the time-dependent wave function, which becomes more pronounced for longer times. 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 CCS method represents a fully quantum initial value representation and accurately reproduces the results of a standard TDSE solver. Finally, we show that the HK propagator may be

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal

citation and DOI.

New Journal of Physics 16 (2014) 103040 1367-2630/14/103040+26$33.00

©2014 IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

successfully employed to compute the time-dependent dipole acceleration and high-harmonic spectra. Nevertheless, the outcome of the semiclassical computation exhibits disagreements with the TDSE, as a consequence of the previously mentioned dephasing.

Keywords: phase space dynamics, strong laser fields, semiclassical propagators, ionization, high-harmonic generation

1. Introduction

Initial-value representations (IVRs) such as the coupled coherent states (CCS) method [1] and the Herman-Kluk (HK) 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 (for discussions related to the CCS and HK methods applied to multidimensional systems, see [3-5], respectively). This efficiency may be increased by employing several strategies, such as dominant Hamiltonians in specific phase-space regions [6, 7] or quantum-state reprojection [8, 9]. For that reason, IVRs have been successfully used in the modeling of the short-time dynamics of quantum systems with many degrees of freedom, such as polyatomic molecules and small clusters. For instance, the CCS method has been applied to the intramolecular vibrational energy exchange [10], tunneling in large molecules [11], and calculations of Franck-Condon [12] spectra. Further examples include applications of the multiconfiguration Ehrenfest method, which is closely related to the CCS approach, to pyrazine [13], and to the spin-boson model [14].

In principle, all these features make IVRs very attractive to strong-field and attosecond science. In fact, it has been well known for the past two decades that strong-field phenomena, such as high-order harmonic generation (HHG), above-threshold 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 [15]. This implies that electron orbits play a very important role in understanding these phenomena. This has led not only to a myriad of applications, such as the attosecond imaging of dynamic processes in matter (see [1618]), but to the extensive use of orbit-based approaches (see [19]).

In particular, classical and semiclassical methods are very popular. Classically, ensembles of electrons that behave according to the previously mentioned recollision picture are constructed 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 [20] and the V-shaped structure observed in NSDI [21, 22]. These methods, however, cannot account for quantum-interference effects or tunnel ionization. Most quantum mechanical and semiclassical methods in strong-field physics 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 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 [20, 2325] and fan-shaped interference patterns [26, 27] in ATI and strongly influences NSDI in circularly polarized fields [28, 29]. Because this interplay is important, in the past few years, Coulomb-corrected analytic approaches have been developed and successfully applied to strong-field phenomena [24, 30-32]. These approaches, however, require the external field to be dominant.

On the other hand, IVRs have only been employed in strong-field physics in relatively few publications. Specifically, in [6, 7, 33] HHG spectra have been computed using the HK propagator and in [34-36] NSDI ion and electron momentum distributions have been calculated employing the CCS method. 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 what extent (see [37-41]). To circumvent this problem, in [6, 7, 33] 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 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 quasistatic 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 [42].

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 IVRs. This analysis will be performed in phase space, for reduced dimensionality models, under the assumption that the electronic wave packet is initially bound. These reduced-dimensionality models are a first step toward more realistic, multidimensional computations and will serve to assess the ionization dynamics in a simplified setting. A further key advantage is that the full solution of the time-dependent Schrödinger equation (TDSE) is straightforward for one-dimensional (1-D) models, whereas for more degrees of freedom this does not hold. For this reason, the TDSE will be used as a benchmark throughout this paper.

This article is organized as follows. In section 2, we provide the necessary theoretical background to understand our results. Subsequently, in section 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-the-barrier ionization, in this case we will also perform a comparison

with the CCS method. In section 4, we present approximate estimates, in which the transmission through the potential barrier is computed analytically using uniform Wentzel-Kramers-Brillouin (WKB) approximations and a comparison with an inverted harmonic oscillator is made. In section 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 section 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.

2. Background

2.1. Model

For the sake of simplicity, we consider a one-electron, 1-D atom. The Hamiltonian of such a system can be expressed as

. p2 „ „

h = Y + Va + Vs, (1)

where Va and V represent 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 length gauge has been chosen, as it leads to the physical picture of an electron tunneling through an effective potential barrier or undergoing over-the-barrier ionization. This picture will be very convenient for the phase-space analysis performed in this work.

The external field exhibits the same temporal dependence as in [6, 33], i.e.,

VS(x, t) = Sx cos (mt), (2)

where S and m denote the field strength and frequency, respectively, of the laser field. The binding potential is taken either as the soft-core potential

Vsc(x) = - *—, (3)

Vx2 + x

or as the short-range Gaussian potential

Vg(x) = -exp (-Xx2), (4)

where X is an arbitrary parameter. As a benchmark, we employ the full solution of the TDSE

idtl Y (t)> = H\Y (t )>, (5)

which is solved in coordinate space for Y (x, t) = (x IY (t)) using a split operator method. The initial wave function is defined as the Gaussian wave packet

<x| Y(0)> = ( £ exp { -2(x - qa)2 + ipa (x - qa) } (6)

centered at pa, qa, with qa = 0. Furthermore, the parameter y defines the width of the wave packet. This initial choice facilitates the implementation of IVRs.

The initial energy expectation value, computed for the system in the absence of the external driving field, is given as

If the wave packet is initially centered at q

(E(t = 0)> = 4 + f + (r(0)JVa(*)| r(0)).

= 0, this expectation value reads as

r(0)| VG(x)\ r(0)} =

r(0)|VSC(*)| r(0)} = -((()1/2^0((xjexp((j

for the Gaussian and the soft-core potentials, respectively, where Kn(z) is the modified Bessel function of the second kind. Equation (7) depends on the width and the initial momentum of the initial wave packet and on the potential parameter A. Throughout, this parameter has been chosen as A = 1 for the soft-core potential and as A = 1/2 for the Gaussian potential.

The width and the position of the initial wave packet were chosen to minimize its initial energy. This procedure yields (pa, qa) = (0, 0) and y ~ 0.46 and (pa, qa) = (0, 0) and y ~ 0.65 for the soft-core and the Gaussian potentials, respectively. The ground-state energies associated with the potentials Vsc(x) and VG(x) are Esc = -0.67 a.u. and EG = -0.594 a.u., respectively.

2.2. Quantum and semiclassical IVRs

Although quantum mechanics is usually formulated in coordinate or momentum space, it can also be formulated in phase space. In this section, we present the IVRs in phase space employed in this work. These representations are the CCS method, which is a formally exact approach for describing quantum mechanics, and the HK propagator, which is semiclassical. More details can be found in the review article [1].

The CCS method represents a time-dependent wave function as a superposition of nonorthogonal Gaussian coherent states (CS), which can either be chosen as static, or as time dependent and trajectory-guided. In this work, we will adopt the latter formulation, i.e., Iz> = Iz (t) >.

Gaussian CS are eigenstates of the creation and annihilation operators

a\z) = z |z) and = (z|z*,

where the eigenvalue z (q, p) is a complex number in phase space parametrized as

2q + T2Y"p'

z* = §■

which is defined in terms of the position q and the momentum p of the particle. CS are not orthogonal; their overlap is given as

(zi\zj) = üij = exp

The identity operator in the CS representation reads as

f\zXd ^,

where d2z = dqdp/2 is the notation for phase-space integration. In this case, the time-dependent wave function is explicitly given by

\Y (t)> = J \z)Dz (t)ei

where the expression S = J

¡s, d2z

*dz dz* z*— - z— dt dt

denotes the classical action along the trajectory defined with regard to the matrix element Hord(z*, z) = (z IHHord(iit, a)I z). This represents the diagonal elements of the ordered Hamiltonian matrix HHord(iit, a). In general,

(z|z')Hord(z*, z'),

where Iz), Iz') denotes two arbitrary CS.

By inserting the identity operator in the TDSE (5) and closing it with (z I, one obtains the integro-differential equation

dDz' (t) dV

r dDz'(t) r, d2z'

J expLl(Sz'- Sz=

- i J (z|z')52Ho*,(z*, z') exp[i(Sz'- Sz)'(t)-for the amplitudes Dz' (t), where the coupling kernel is

52Hor*d(z*, z') = Hord(z*, z') - Hord(z *, z') - idt ((* - z'*)

with the matrix elements of the ordered Hamiltonian defined in equation (16). An important property of (18) is that the kernel 52Ho*j(z*, z') is always small for z and z close to each other and their overlap (zIz) vanishes when they are far away.

The exact form of the kernel depends on the choice of the trajectories, which are used to guide the basis, and it becomes particularly simple if they are determined by Hamilton's equations with a quantum-averaged Hamiltonian

dz ' . dHord (z *, z')

In coordinate space

-—(x - q)2 + ip (x - q) + ipq

2 2 is a Gaussian wave packet centered at the phase-space coordinates q and p.

In practice, a finite basis of Gaussian CS is used so that the integral in (14) becomes a finite sum. The identity operator is discretized as

I = /lV z>, (21)

where Q-1 is the inverse of the overlap matrix Q, = {n\z,). This specific discretization reduces the integro-differential equation (17) to a system of linear equations for the derivatives Dz (t), which are solved numerically together with equation (19) for the trajectories.

For the HK propagator, the time-dependent wave function is once more expressed in terms of Gaussians in phase space (for seminal papers, see [2, 43-45]). Following the original work [2], they are usually labeled, not with a single complex number z, but with two real numbers p and q. This specific formulation is considered in our implementation and in the results that follow. Explicitly,

|SHK(t)) = //1q, P)R(t, qo, Po))qo, Po\T(°))e*d—(22)

where \q, p) represents a coherent state whose expression in coordinate representation is given by

<x|q, p> = ( n) exp[ -- q)2 + ip(x - q)], (23)

which differs from equation (20) employed in the CCS approach by a phase factor, and

1 I i Y'2

R ( qo, Po) = 1 mpp + mqq- iymqp + Ympq) (24)

is given in terms of the elements muv = du/dvo of the monodromy matrix, which is composed of the derivatives of the final variables q, p with regard to their initial values qo, po. In equation (22), the action reads as

Sci(q, p) = / (pq - Hi)dt, (25)

where Hcl is the classical Hamiltonian, in which the operators x, p in (1) have been replaced by the phase-space variables q, p. For the initial wave packet (6) considered here,

^ po(o)) = exp j - 4(qa - qo)2 - 4Y p - po)2 + ^(a + po)o - qa ) }• (26)

The integral is carried out over phase-space coordinates, which are used as initial conditions of the classical solutions (q, p). In the p, q-notation, the phase of CS in (23) differs from that of z-notations by ipq/2, which is compensated by the different form of action (25) used in [2] as opposed to that of equation (15).

Apart from trivial notations, the HK propagator and the CCS method differ in three main ways:

(1) The trajectories of the Gaussian CS in the HK method are purely classical, whereas in CCS, the trajectory of a CS \z (t)) is driven by the Hamiltonian {z \H\ z). This latter Hamiltonian is the average of the quantum Hamiltonian with regard to the CS and

considers 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 potentials.

(2) The CCS representation is a formally exact basis set technique. It uses coupled quantum equations for the coefficients Dz(t), obtained simply by substitution of (14) into the Schrödinger equation. 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 considers the first and second terms in the Taylor expansion of the potential energy around a specific trajectory, and assume that only the coupling of CS near this trajectory is important. Physically, this implies that, while in the CS representation the trajectories are coupled through the amplitudes Dz(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 basis sets allow us 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 [8, 9].

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 IVRs, see also [46-48].

3. Phase-space dynamics

In the following, we analyze different ionization mechanisms in phase space. According to the quasistatic tunneling picture, a low-frequency, time-dependent field and the binding potential determine an effective potential barrier

Veff(x, t) = Va (x) + xE(t), (27)

whose maximum will be given by Vfffe, t) at the coordinate xs such that dVff (x, t)/dx I 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.

3.1. 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 equation (22), and the wave functions ¥hk (x, t) and W (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 (27) if the frequency is low enough. In this case, the time-dependent field E (t) is replaced by S in equation (2).

Figure 1. Phase portrait of the system defined by equation (28) for a static electric field. The solid lines represent the separatrix in phase space for driving-field amplitudes of £ = 0.075 a.u. (left) and £ = 0.05 a.u. (right) according to equation (29). Dashed and dotted lines show solutions for energies E = -0.67 a.u. and E = -0.3 a.u., respectively. The dashed-dotted lines illustrate the evolution of some sample trajectories from t = 0 to t = 20 a.u.

In phase space, the classical static Hamiltonian then reads as

q) = p2 + va (q) + £q. (28)

In these studies, we will consider the soft-core potential (3).

In figure 1, we present the phase portrait of the system for the Hamiltonian (28) and static fields of different amplitudes. The figure shows that the condition

p2/2 + Va (q) + £q = Emin , (29)

where Emin is the minimal energy necessary for the electron to undergo over-the-barrier ionization, defines a separatrix, which crosses at (q, p) = (qs, 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 (qs, 0). This means that an orbit with E < Emin will remain unbound if its initial spatial coordinate q0 lies on the left-hand side of the saddle point, and that it will remain bound if q0 lies on the right-hand side of qs (figure 1). On the other hand, orbits with E > Emin 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 Emin always remain bound and propagate along closed orbits bounded by the separatrix. Those that follow the separatrix from below (Emin < 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 phase-space 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 vice versa. This is shown in the upper panels of figure 2, where we display the phase-space representation (26) of two Gaussian wave packets (6) centered at (qa, pa) = (0, 0) with different widths y = 0.5 and y = 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

Figure 2. Classical trajectory distributions, together with the initial quantum mechanical wave packet, considering a soft-core potential and a static field of amplitude £ = o.o75 in equation (28). Upper panels: initial quantum mechanical distributions in phase space for a Gaussian wave packet centered at (qa, pa) = (o, o), y = o.5 (bound-state energy E ^ -o.67 a.u.; left), and y = o.o5 (bound state energy E ^ -o.5 a.u.; right). Middle panels: initial conditions in phase space for an ensemble of classical trajectories corresponding to these distributions, represented by dots. Lower panels: classical distributions at t =20 atomic units for an ensemble whose initial conditions are given by those in the middle panels [(qa, pa) = (o, o) and, from left to right, y = o.5 = and Y = o.o5]. In our computations using the HK method, we have considered around 107 trajectories. The thick black lines indicate the separatrix defined according to equation (29) and the constant-energy curve for E =0 a.u.

-1.5 ....................

-15 -10 -5 0 5

Figure 3. Time-dependent separatrix corresponding to the laser field represented by equation (2) (solid lines), together with a specific electron trajectory (dashed line), with frequency m = 0.05 a.u. and amplitude £ = 0.075 a.u. The gradient color shows the time variation from t =0 (blue) to t =30 a.u. (red), which spans slightly less than the first quarter cycle of the field. The dot depicts the initial condition of an initially unbound trajectory that 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 > T14, where T = 2nlm is the field cycle, the region bounded by the separatrix starts to decrease, until the trapped trajectory is eventually able to escape.

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 figure 2. These are the starting points (q0, p0) 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 figure 2), the distributions follow the constant-energy curves in phase space. For the wave packet with y = 0.5, the distributions are concentrated in the bound region or, for q < qs, 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 (y = 0.05), there are phase-space events in this region, but closely following the separatrix given previously.

These features are determined by the initial position and momentum spreads of the wave packet and the corresponding trajectory ensemble. For y = 0.5, the wave packet 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 < Emin 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 (y = 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 previously determined constraints vary in time may a classical bound trajectory become unbound or vice versa. This is what happens if the external field is time dependent. In figure 3, we illustrate this behavior for the laser field associated with equation (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.

3.2. Wigner quasiprobabilities

We will now employ the Wigner quasiprobability distribution (also known as the Wigner function) to relate the preceding picture to the wave-packet propagation in phase space [49]. This function is defined as

W (q, p) = - I dy¥* (q + y) ¥ (q - y) exp [2ipy]. (30)

If equation (30) 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 to assess the dynamics of ionization [50], rescattering [51, 52], double ionization [53], and HHG [54]. A common feature identified in [50] and [54] 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, W2(q, p) as it makes the previously mentioned tail slightly clearer.

In the left and right panels of figure 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 (qa, pa) = (0, 0) and width y = 0.5 (see figure 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-the-barrier ionization: the electronic wave packet does not leave the core with vanishing momentum, but, rather, with the minimum necessary momentum to overcome Veff(qs). As time progresses, interference fringes start to build up on the left-hand side just above the separatrix. These fringes have been identified in [50], for a delta-potential model in a static field, and in [51, 54] for long- and short-range potentials in time-dependent 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 = qs, the momentum at this tail is even approximately vanishing, in agreement with the model in [50]. 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 [39].

All previously mentioned features are present both for the semiclassical and quantum mechanical computations. In fact, the agreement between the outcomes of both approaches is

-15 -10 -5 0 5 -15 -10 -5 0 5

Figure 4. Square of the Wigner function calculated for a wave packet in a soft-core potential (3) and a static field of amplitude £ = 0.075 a.u. by using the full quantum wave function (left column) and the semiclassical wave function (right column) for t =10 a.u. (first row) and t =20 a.u. (second row). The thicker lines in the figure show the separatrix and the phase-space trajectory for E = 0. The width and initial momentum of the initial wave packet are y = 0.5 and pa = 0, respectively.

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 wave packet to cross classically forbidden regions. In other words, classical trajectories with over-the-barrier energy appear to mimic transmission through a barrier for E < Emin . This is not obvious, as these trajectories do not violate any phase-space constraints, or cross classically forbidden regions.

However, we have verified that the trajectories whose initial coordinates (q0, p0) lie outside the bound phase-space region lead to the tail in the Wigner distributions. For clarity, in figure 5 we show the Wigner function computed from the HK propagator leaving out these trajectories. In the figure, both the previously mentioned tail and the interference fringes are absent. These findings are in agreement with the results in [55], in which a nonlocal behavior of the Wigner function around separatrices has been observed for an inverted harmonic oscillator, and with those in [56], 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 < qs also supports the assumption that the Wigner function behaves nonlocally.

For a short-range potential, the separatrix becomes much steeper and the tail in the Wigner function starts to form along the separatrix, i.e., for lpI ~ Ipmin I. This is explicitly shown in

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 figure 4 and times t =10 a.u. and t = 20 a.u. (left and right panels, respectively).

Figure 6. Square of the Wigner function calculated using the TDSE and the HK propagator (left and right panels, respectively), for the same field parameters as in figure 4, time t =10 a.u., and a truncated soft-core potential according to equation (31), with a0 = 1 and L = 7.55. The solid and dashed lines indicate the separatrices for the truncated and the long-range soft-core potential, respectively.

figure 6, in which the soft-core potential has been truncated by multiplying equation (3) with the function

cos 0,

n |x| - a0 2 L - a 0

|x| < a0 a0 < |x| < L, \x\ > L

in which the parameters a0 and L are chosen so that the ground-state energy remains practically the same but the effective potential barrier changes. This type of potential has been employed in our previous publications [57, 58] to study the influence of the Coulomb tail on HHG enhancements, and HHG spectra from specific Bohmian trajectories, respectively. Physically, this may be understood as follows: although for a long-range potential the electron may escape via highly excited states and over-the-barrier ionization, for short-range potentials these

Figure 7. Square of the Wigner function computed using the TDSE and the HK propagator (left and right panels, respectively) employing a soft-core potential (3) and a time-dependent field associated to equation (2), with frequency m = o.o5 a.u. and amplitude £ = o.o75 a.u. The top, middle, and bottom panels have been calculated for t = 12.6 a.u., t = T/4, and t = T/2, respectively, where T = 2n/m is the field cycle.

pathways are strongly suppressed. Hence, the wave packet is expected to escape through phasespace regions associated with tunnel ionization.

For the time-dependent field in equation (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 [51, 52]). Snapshots of W2(q, p) in the time-dependent case are presented in figure 7. 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

-15 -10 -5 0 5 -15 -10 -5 0 5 -15 -10 -5 0 5

Figure 8. Modulus square of the Wigner quasiprobability distributions computed using a static field of the same amplitude as in figure 4 and the Gaussian potential (4), using the TDSE, the HK propagator, and the CCS method (left, middle, and right panels, respectively). The upper and lower panels correspond to t =10 a.u. and t =20 a.u., respectively. The separatrix and the curve in phase space for the energy E = 0 are illustrated by the thick lines in the figure. For the HK propagator and the CCS method, we employ 107 and 1600 trajectories, respectively.

and the TDSE worsens with time. This is related to the fact that nonclassical effects such as tunnel ionization and over-the-barrier reflections cannot be fully accounted for, and become dominant for longer times [39].

3.3. Comparison with the CCS method

We will now investigate the Gaussian potential (4) and an initial wave packet with y = 1. As discussed in the previous section, 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. A Gaussian potential is very convenient for the subsequent comparison of the ordered and the classical Hamiltonians, as it may be computed analytically in a coherent-state basis.

In figure 8, we display snapshots of W2(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 given previously. This behavior holds even for the results obtained with the semiclassical propagator and is different from that observed for the soft-core potential (see figure 4). Physically, this suggests that there will be a nonvanishing probability density leaving the core with Ip I < Ipminl. 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. 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 [59] (see also [8, 9] for details).

We have also verified that, although the trajectories employed in the HK propagator never cross the separatrix, defined with regard to the classical Hamiltonian those used in the CCS method do. This is a consequence of the fact that the ordered Hamiltonian Hord (z*, z) is different from the classical Hamiltonian Hci(p, q). In fact, as a function of the phase-space coordinates (q, p), for a static field Hord (z*, z) is given by

2 / \l/2

HoldCp, q) = J + y - I y+j I exp [-nq2] + q£, (32)

with n = (Jy/(y + J)). If the separatrix is however defined with regard to Hord(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 Hos^d(p, q) and HCi1t(p, q) shows an effective energy shift y/4 and a shift

\1/2 f 32„2 I

AVg (q) = Vg (q)

y 1 exp- 1

y + J I I y + J

in the binding potential, where VG(q) is given by equation (4). For discussions of these shifts, see [47, 48].

4. 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 [55]. In [39], 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 Veff(x) around the saddle xs, and, using the uniform WKB approximation, we compute the transmission coefficient through this barrier for a static field. This coefficient reads as

P (E) =

1 + exp j 2 ^ ^ J 2(Veff(x) - E) dxj

where xi and xr are the left and right turning points, respectively, for which Veff(x;) = V;ff(xr) = E .In figure 9, 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.

S 10"2 h

10"4 -0.75

exact ..........

second order -

fourth order

sixth order o

„ ,..-'**"'0.025

i ° 0.015 -

0.005 ---------- , . , -

-0.7 -0.68 -0.66

-0.65 E

Figure 9. Transmission coefficients for different order approximations to the barrier in the soft-core (left) and the inverted Gaussian potentials (right) for a static field of amplitude £ = 0.075 a.u. The insets show the behavior of such coefficients within the range in which the energy of the respective ground state lie.

A better agreement is observed for the soft-core potential. This is consistent with the results in the previous section (see figures 4 and 8).

5. 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. HHG spectra are a particularly critical test for the present IVRs as HHG is a coherent process, for which the phase of the wavefunction is important and which builds up over at least a few cycles. Furthermore, it occurs close to the core, for which the interplay between the binding potential and the external laser field is important.

The HHG spectra are given by

J dta (t )exp (iQt )

dx¥*(x, t)^^¥(x, t), dx

denotes the dipole acceleration. The time-dependent wave function is either given by the solution of equation (5) in coordinate space or by ¥Hk (x, t). The acceleration form of the dipole operator has been taken as it emphasizes regions near the core, and it is numerically much less expensive than the dipole length for the HK method. However, we have verified that the length form of the dipole operator leads to similar results (not shown). In this latter case, it was necessary to employ a Hanning filter to remove the existing numerical background. The high-frequency part of the spectrum from the dipole length is underestimated by a factor QA, in comparison to that of the dipole acceleration. All these features are expected and consistent with what is known about the different forms of the dipole operator (see [60, 61] for discussions of these issues).

Laser periods

Harmonic Order

Figure 10. Dipole acceleration over two field cycles, together with the HHG spectra (left and right panels, respectively), computed using a laser field of frequency m = 0.05 a.u. and amplitude £ = 0.075 a.u., the soft-core potential (3) and an initial wave packet (6) of width y = 0.5, centered at qa = 0 and with vanishing momentum pa = 0. Panels (a) and (b) were computed using the full TDSE computation, and the HK propagator, respectively. The dashed lines show the energy position of the cutoff, which in this case is located at Ip + 3.17 Up = 49 m.

In figure 10, we display a(t), together with its power spectrum, for an initial wave packet (6) centered at qa = 0 and with vanishing momentum pa = 0, whose width is y = 0.5 a.u. The figure shows a qualitative agreement between the fully quantum and semiclassical acceleration [panels (a) and (b), respectively], which roughly follow the field and exhibit a series of 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 [62] and the HK propagator [6, 7, 33], for an initial wave packet far from the core. They have also been studied in a different context, namely the adiabatic approximation [63] and Bohmian trajectories [58, 64]. The agreement between the outcome of the TDSE and the HK

propagator is particularly good for times below half a cycle. However, phase differences arise between 0.5 T and 1.5 T, with T = 2n/m.

This agreement persists for the spectra, which exhibit a long plateau followed by a sharp cutoff at Ip + 3.17Up 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 previously.

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 wave packet is increased, for instance, by changing its width or initial momentum pa. This, together with an overall decrease in the plateau height, is consistent with the fact that, for pa # 0 or y ^ 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 [38]. These issues may constitute a problem for long times [39]. We have observed that our results are reasonably accurate for t ^ 3.5 T.

One should also note that the spectrum from the HK computation exhibits a second, much less intense plateau. A similar plateau, extending to Ip + 8Up, has been identified in [65], in the context of an improved SFA method, and has been attributed to ionization and recombination at different points in space and/or ionization with initial nonvanishing velocities. Although we cannot state with certainty whether the mechanism in [65] is behind the second plateau observed in this work, our model does account for the spatial extension of the core and an initial velocity spread. Other possibilities for this second plateau are the temporal degradation of the HK computation and classical trajectories becoming trapped in orbits whose energies are lower than that of the ground state (in fact until the energy — 1/VX defined by the potential minimum). We have ruled out the latter possibility, as the removal of such trajectories does not eliminate the second plateau. Furthermore, we have found no obvious connection between the degradation of the HK propagator and this plateau.

Finally, in figure 11, we display a distribution of trajectories employed in the HK propagator as functions of the kinetic energy obtained from the field at the recombination time. This kinetic energy is expected to be converted in harmonic radiation. In the figure, we have assumed that the electron returns to the core when the condition \x I ^ r0, where r0 is the Bohr radius, is met. This implies that recombination to the ground state has occurred. Because our model incorporates an initial velocity spread, for each orbit we have subtracted the initial kinetic energy from its value upon return to obtain the kinetic energy transfer from the field. A similar procedure has been employed in [66] and in our previous publication [58] for more simplified, classical trajectory models.

To discard the trajectories that remain in the bound phase-space region, we only considered the kinetic energies starting from the times in which the trajectories reach the saddle, i.e., at x (ts) = xs. The figure exhibits a clear cutoff around Ip + 3.17Up, in agreement with the three-step model. We have verified that, if this initial assumption is relaxed and we consider all trajectories starting at xi (0) according to the distribution in figure 2, there is a substantial increase in the low-energy trajectories and a more uniform distribution of trajectories across the plateau (not shown). This is expected according to the preceding discussion.

A Ek/u

Figure 11. Histogram constructed with the trajectories that go beyond the saddle xs as functions of the net kinetic energy AEk/m = [Ek (t) - Ek (ts)]/m gained by the electron from the field in units of the field frequency m, where Ek(t) and Ek (ts) denote the kinetic energy at the return time and at the time ts for which the electrons reach the saddle, respectively. We consider the same field and atomic parameters as in the previous figure. A total sample of 5 X 105 trajectories has been used, and the return condition has been chosen as Ix (t) I < 1 a.u. The dashed line in the figure indicates the cutoff energy.

6. Conclusions

The results presented in this paper strongly suggest that semiclassical IVRs, a concrete example of which is the HK 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 TDSE, or the 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 zero-range potential [50, 52] or the TDSE [54]. 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 [15]. 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-the-barrier ionization, whereas 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, although the trajectories in the semiclassical method used in this paper are classical and will never cross a separatrix, the initial position-momentum 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 [55] for an inverted harmonic oscillator.

On the other hand, if the classical trajectories that start outside the bound phase-space 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 IVR. Furthermore, a semiclassical approximation would only be exact for a parabolic barrier, whereas 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 [39, 67] as a criticism to the findings in [55]. 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 [39, 67], the potential barrier is flat, i.e., limx^xVa(x) = 0, whereas, for the effective potential Veff(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 [39, 67]. Because 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 of [40, 68] or by adding higher-order correction terms to the HK propagator as in [69, 70]. Still, the results in section 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.

Nevertheless, the HK propagator did reproduce low-energy regions in the HHG spectra to a reasonable extent (e.g., the hump around the frequency Q = 10 m), for which the application of standard strong-field methods such as the SFA is questionable. Furthermore, a very good agreement in the cutoff region has been obtained. This implies that, at the interpretational level, semiclassical IVRs provide a tool for investigating quantum interference between different types of trajectories returning to the core, in which the binding potential is incorporated. This holds even if the electronic wave packet is initially bound. This is a noteworthy fact in view of the wide range of applications for quantum interference in strong-field physics, especially in the attosecond imaging of matter.

Acknowledgements

We are grateful to A Fring, A Olaya-Castro, A Serafini, J Guo, X S Liu, and in particular to F Grossmann and J M Rost for very useful discussions. CFMF, CZM, and JW would like to thank the University of Leeds and the Max Planck Institute for Physics of Complex Systems, Dresden, where the final stages of this work were carried out, for their kind hospitality. This work was funded by the UK EPSRC (Grant EP/J019143/1) and the CSC/BIS (China-UK Studentship for Excellence).

Appendix. Relation between the CCS and the HK IVRs

For the reader's convenience, in this appendix we briefly focus on how the CS representation and the HK propagator share a common origin. The details of the derivation of CCS working equations and the HK formula can be found in [1]. Originally, the semiclassical phase-space HK method has been derived differently [2]. Nonetheless, the current derivation shows that it also can be obtained as an approximation of the exact quantum dynamics in phase space.

The HK method can be derived from equation (17) by employing the local quadratic approximation, which only considers the first and second terms in the Taylor expansion of the potential energy around a specific trajectory. This allows us to use classical trajectories instead of those driven by the quantum averaged Hamiltonian Hord (z*, z). Furthermore, the local quadratic approximation assumes that only the CS that are very close to each other, and therefore are driven by the same quadratic potential, are coupled. Under these assumptions, the integrals in (17) may be calculated analytically.

The resulting wave function may be then represented as

(t)> = / |z>VM^ei5zcl(zo| Y (0))

where ^Mzz denotes the HK prefactor in the z notation. In this notation Mzz is a single element of the monodromy (or stability) matrix

' Mzz Mzz* ' Mz*z Mz*z*

which describes how stable the dynamics around a specific trajectory are. The monodromy matrix elements in this representation are related to those in the p, q representation by

2-1(mqq + mpp - iymqp + iy~lMpq), 2-1(mqq - mpp + iymqp + iy-1mpq), 2-1(mqq - mpp - iymqp - iy-1mpq),

1 (mqq + mpp + iymqp - iy lmpq),

where the elements Mv^ are defined in section 2.2. The evolution of the monodromy matrix is given by the matrix of second derivatives of the Hamiltonian

d2Hcl/dz*dz d2Hcl/dz* -d2Hcl/dz2 - d 2Hcl/dz*dz

zz* Mz*z*

In this work, only a 1D case has been considered. However, a generalization to more than one degree of freedom is straightforward, as a multidimensional CS is simply a product of 1-D CSs. The CCS method can also be generalized for multielectronic states by introducing Fermionic CS with proper permutational symmetry [35].

References

[1] Shalashilin D V and Child M S 2004 The phase space CCS approach to quantum and semiclassical molecular

dynamics for high-dimensional systems Chem. Phys. 304 103

[2] Herman M F and Kluk E 1984 A semiclasical justification for the use of non-spreading wavepackets in

dynamics calculations Chem. Phys. 91 27

[3] Shalashilin D V and Child M S 2008 Basis set sampling in the method of coupled coherent states: coherent

state swarms, trains, and pancakes J. Chem. Phys. 128 054102

[4] Ronto M and Shalashilin D V 2013 Numerical implementation and test of the modified variational

multiconfigurational gaussian method for high-dimensional quantum dynamics J. Phys. Chem. A 117 6948

[5] Brewer M L 1999 On the scaling of semiclassical initial value methods J. Chem. Phys. 111 6168

[6] Zagoya C, Goletz C-M, Grossmann F and Rost J M 2012 Dominant-interaction Hamiltonians for high-order

harmonic generation in laser-assisted collisions Phys. Rev. A 85 041401(R)

[7] Zagoya C, Goletz C-M, Grossmann F and Rost J M 2012 An analytical approach to high harmonic generation

New J. Phys. 14 093050

[8] Shalashilin D V and Jackson B 2000 Guiding paths and time-dependent basis sets for wavefunction

propagation Chem. Phys. Lett. 318 305

[9] Burant J C and Batista V S 2002 Real time path integrals using the Herman-Kluk propagator J. Chem. Phys.

116 2748

[10] Shalashilin D V and Child M S 2003 Nine-dimensional quantum molecular dynamics simulation of

intramolecular vibrational energy redistribution in the CHD3 molecule with the help of coupled coherent states J. Chem. Phys. 119 1961

[11] Sherratt P A J, Shalashilin D V and Child M S 2006 Description of multidimensional tunnelling with the help

of coupled coherent states guided by classical hamiltonians with quantum corrections Chem. Phys. 322 127

[12] Shalashilin D V and Child M S 2004 Real time quantum propagation on a monte carlo trajectory guided grids

of coupled coherent states: 26d simulation of pyrazine absorption spectrum J. Chem. Phys 121 3563

[13] Shalashilin D V 2010 Nonadiabatic dynamics with the help of multiconfigurational ehrenfest method:

improved theory and fully quantum 24d simulation of pyrazine J. Chem. Phys. 132 244111

[14] Shalashilin D V 2009 Quantum mechanics with the basis set guided by Ehrenfest trajectories: theory and

application to spin-boson model J. Chem. Phys. 130 244101

[15] Corkum P B 1993 Plasma perspective on strong field multiphoton ionization Phys. Rev. Lett. 71 1994

[16] Lein M 2007 Molecular imaging using recolliding electrons J. Phys. B: At. Mol. Opt. Phys. 40 135

[17] Krausz F and Ivanov M 2009 Attosecond physics Rev. Mod. Opt 81 163

[18] Haessler S, Caillat J and Salières P 2011 Self-probing of molecules with high harmonic generation J. Phys. B:

At. Mol. Opt. Phys. 44 203001

[19] Salières P et al 2001 Feynman's path-integral approach for intense-laser-atom interactions Science 292 902

[20] Quan W et al 2009 Classical aspects in above-threshold ionization with a midinfrared strong laser field Phys.

Rev. Lett. 103 093001

[21] Emmanouilidou A 2008 Recoil collisions as a portal to field-assisted ionization at near-UV frequencies in the

strong-field double ionization of helium Phys. Rev. A 78 023411

[22] Ye D F, Liu X and Liu J 2008 Classical trajectory diagnosis of a fingerlike pattern in the correlated electron

momentum distribution in strong field double ionization of Helium Phys. Rev. Lett. 101 233003

[23] Blaga C I, Catoire F, Colosimo P, Paulus G G, Muller H G, Agostini P and DiMauro L F 2009 Strong-field

photoionization revisited Nat. Phys. 5 335

[24] Yan T-M, Popruzhenko S V, Vrakking M J J and Bauer D 2010 Low-energy structures in strong field

ionization revealed by quantum orbits Phys. Rev. Lett. 105 253002

[25] Telnov D A and Chu Sh-I 2011 Low-energy structure of above-threshold-ionization electron spectra: role of

the Coulomb threshold effect Phys. Rev. A 83 063406

[26] Rudenko A, Zrost K, Schröter C D, de Jesus VLB, Feuerstein B, Moshammer R and Ullrich J 2004

Resonant structures in the low-energy electron continuum for single ionization of atoms in the tunnelling regime J. Phys. B: At. Mol. Opt. Phys. 37 L407

[27] Arbo D G, Yoshida S, Persson E, Dimitriou KI and Burgdörfer J 2006 Interference oscillations in the angular

distribution of laser-ionized electrons near Ionization Threshold Phys. Rev. Lett. 96 143003

[28] Mauger F, Chandre C and Uzer T 2010 Recollisions and correlated double ionization with circularly

polarized light Phys. Rev. Lett. 105 083002

[29] Kamor A, Mauger F, Chandre C and Uzer T 2013 How key periodic orbits drive recollisions in a circularly

polarized laser field Phys. Rev. Lett. 110 253002

[30] Popruzhenko S V, Paulus G G and Bauer D 2008 Coulomb-corrected quantum trajectories in strong-field

ionization Phys. Rev. A 77 053409

[31] Smirnova O, Spanner M and Ivanov M 2006 Coulomb and polarization effects in sub-cycle dynamics of

strong-field ionization J. Phys. B: At. Mol. Opt. Phys. 39 307-21

[32] Smirnova O, Spanner M and Ivanov M 2008 Analytical solutions for strong field driven atomic and

molecular one- and two-electron continua and applications to strong-field problems Phys. Rev. A 77 033407

[33] van de Sand G and Rost J-M 1999 Irregular orbits generate higher harmonics Phys. Rev. Lett. 83 524

[34] Shalashilin D V, Child M S and Kirrander A 2007 Mechanisms of double ionisation in strong laser field from

simulation with coupled coherent states: beyond reduced dimensionality models Chem. Phys. 347 257

[35] Kirrander A and Shalashilin D V 2011 Quantum dynamics with fermion coupled coherent states: theory and

application to electron dynamics in laser fields Phys. Rev. A 84 033406

[36] Guo J, Liu X-S and Chu Shih-I 2010 Exploration of strong-field multiphoton double ionization, rescattering,

and electron angular distribution of He atoms in intense long-wavelength laser fields: the coupled coherent-state approach Phys. Rev. A 82 023402

[37] Keshavamurthy S and Miller W H 1994 Semi-classical correction for quantum-mechanical scattering Chem.

Phys. Lett. 218 189

[38] Grossmann F and Heller E J 1995 A semiclassical correlation function approach to barrier tunneling Chem.

Phys. Lett. 241 45

[39] Maitra N T and Heller E J 1997 Barrier tunneling and reflection in the time and energy domains: the battle of

the exponentials Phys. Rev. Lett. 78 3035

[40] Kay K 1997 Semiclassical tunneling in the initial value representation J. Chem. Phys. 107 2313

[41] Spanner M 2003 Strong field tunnel ionization by real-valued classical trajectories Phys. Rev. Lett. 90 233005

[42] Figueira de Morisson Faria C and Liu X 2011 Electron-electron correlation in strong laser fields J. Mod. Opt.

58 1076

[43] Kluk E, Herman M F and Davis H L 1986 Comparison of the propagation of semiclassical frozen Gaussian

wave functions with quantum propagation for a highly excited anharmonic oscillator J. Chem. Phys. 84 326

[44] Herman M F 1986 Time reversal and unitarity in the frozen Gaussian approximation for semiclassical

scattering J. Chem. Phys. 85 2069

[45] Kay K G 1994 Integral expressions for the semiclassical time-dependent propagator J. Chem. Phys. 100 4377

[46] Grossmann F and Xavier A L Jr 1998 From the coherent state path integral to a semiclassical initial value

representation of the quantum mechanical propagator Phys. Lett. A 243 243

[47] Miller W H 2002 On the relation between the semiclassical initial value representation and an exact quantum

expansion in time-dependent coherent states J. Phys. Chem. B 106 8132

[48] Child M S and Shalashilin D V 2003 Locally coupled coherent states and Herman-Kluk dynamics J. Chem.

Phys. 118 2061

[49] Wigner E 1932 On the quantum correction for thermodynamic equilibrium Phys. Rev. 40 749

[50] Czirjak A, Kopold R, Becker W, Kleber M and Schleich W P 2000 The Wigner function for tunneling in a

uniform static electric field Opt. Comm. 179 29

[51] Kull H-J 2012 Position-momentum correlations in electron-ion scattering in strong laser fields New J. Phys.

14 055013

[52] Czirjak A, Majorosi S, Kovács J and Benedict M G 2013 Emergence of oscillations in quantum entanglement

during rescattering Phys. Scr. 153 014013

[53] Lein M, Engel V and Gross E K U 2001 Phase-space analysis of double ionization Opt. Express 8 411

[54] Gräfe S, Doose J and Burgdörfer J 2012 Quantum phase-space analysis of electronic rescattering dynamics in

intense few-cycle laser fields J. Phys. B: At. Mol. Opt. Phys. 45 055002

[55] Balazs N L and Voros A 1990 Wigner's function and tunneling Ann. Phys. 199 123

[56] Heim D M, Schleich W P, Alsing P M, Dahl J P and Varro S 2013 Tunneling of an energy eigenstate through

a parabolic barrier viewed from Wigner phase space Phys. Lett. A 377 1822

[57] Figueira de Morisson Faria C, Kopold R, Becker W and Rost J M 2002 Resonant enhancements of high-order

harmonic generation Phys. Rev. A 65 023404

[58] Wu J, Augstein B B and Figueira de Morisson Faria C 2013 Local dynamics in high-order-harmonic

generation using Bohmian trajectories Phys. Rev. A 88 023415

[59] Symonds C, Wu J, Ronto M, Zagoya C, Figueira de Morisson Faria C and Shalashilin D V 2014 Coupled

coherent state approach for high-order harmonic generation, submitted

[60] Burnett K, Reed V C, Cooper J and Knight P L 1992 Calculation of the background emitted during high-

harmonic generation Phys. Rev. A 45 3347

[61] Krause J L, Schafer K J and Kulander K C 1992 Calculation of photoemission from atoms subject to intense

laser fields Phys. Rev. A 45 4998

[62] Protopapas M, Lappas D G, Keitel C H and Knight P L 1996 Recollisions, bremsstrahlung, and attosecond

pulses from intense laser fields Phys. Rev. A 53 R2933

[63] Okajima Y, Tolstikhin O I and Morishita T 2012 Adiabatic theory of high-order harmonic generation: one-

dimensional zero-range-potential model Phys. Rev. A 85 063406

[64] Wu J, Augstein B B and Figueira de Morisson Faria C 2013 Bohmian-trajectory analysis of high-order-

harmonic generation: ensemble averages, nonlocality, and quantitative aspects Phys. Rev. A 88 063416

[65] Plaja L and Pérez-Hernández J A 2007 A quantitative s-matrix approach to high-order harmonic generation

from multiphoton to tunneling regimes Opt. Express 7 3629

[66] Hostetter J A, Tate J L, Schafer K J and Gaarde M B 2010 Semiclassical approaches to below-threshold

harmonics Phys. Rev. A 82 023401

[67] Maitra N T and Heller E J 1997 Tunneling and the semiclassical propagator: a new perspective Classical,

Semiclassical and Quantum Dynamics in Atoms (Lecture Notes in Physics) ed H Friedrich and B Eckhardt (Berlin: Springer) p 94

[68] Ankerhold J and Saltzer M 2002 Semiclassical wave packet tunneling in real-time Phys. Lett. A 305 251

[69] Zhang S and Pollak E 2003 Monte Carlo method for evaluating the quantum real time propagator Phys. Rev.

Lett. 91 190201

[70] Kay K G 2006 The Herman-Kluk approximation: derivation and semiclassical corrections Chem. Phys. 322 3

Copyright of New Journal of Physics is the property of IOP Publishing and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.