Adiabatic theory of ionization of atoms by intense laser pulses: One

PHYSICAL REVIEW A 81, 033415 (2010)

Adiabatic theory of ionization of atoms by intense laser pulses: One-dimensional zero-range-potential model Oleg I. Tolstikhin,1 Toru Morishita,2,3 and Shinichi Watanabe2 1

Russian Research Center “Kurchatov Institute,” Kurchatov Square 1, Moscow 123182, Russia 2 Department of Applied Physics and Chemistry, University of Electro-Communications, 1-5-1, Chofu-ga-oka, Chofu-shi, Tokyo, Japan 3 PRESTO, Japan Science and Technology Agency, Kawaguchi, Saitama 332-0012, Japan (Received 15 January 2010; published 22 March 2010) As a first step toward the adiabatic theory of ionization of atoms by intense few-cycle laser pulses, here we consider the simplest one-dimensional zero-range-potential model. The theory for this model is developed. This theory rests on a single parameter giving a ratio of the atomic and laser time scales and becomes exact in the limit → 0 uniformly with respect to the amplitude of the laser pulse. We construct the asymptotics of the solution to the time-dependent Schrödinger equation and photoelectron spectrum. The factorization formula for the photoelectron spectrum near its classical boundary, first suggested by T. Morishita et al. [Phys. Rev. Lett. 100, 013903 (2008)] on the basis of ab initio calculations, is derived. The theory is illustrated by numerical calculations. Its comparison with the strong-field approximation is discussed. DOI: 10.1103/PhysRevA.81.033415

PACS number(s): 32.80.Fb, 32.80.Rm, 34.50.Rk, 42.50.Hz

I. INTRODUCTION

The ionization of atoms by laser pulses is one of the fundamental problems of atomic physics. In recent years, due to the great progress in laser technologies, the problem has became important also for applications. From the point of view of applications, the theory must provide a simple qualitative picture of the ionization dynamics as well as an accurate quantitative description of the resulting photoelectron spectrum. The second goal can be achieved by integrating the time-dependent Schrödinger equation (TDSE) numerically. The role of numerical calculations in advancing the field is difficult to overestimate. They are also very helpful for achieving the first goal. However, a qualitative understanding of phenomena always emerges from analytical constructions based on a suitable physical approximation. We focus on intense (I ∼ 1014 –1016 W/cm2 ) few-cycle infrared (λ ∼ 800 nm) laser pulses which are currently of main interest for applications. A theoretical approach capable of treating this range of the laser parameters was proposed by Keldysh [1], Faisal [2], and Reiss [3] and is known as the strong-field approximation. This theory is based on the assumption that at a certain stage of the ionization process, the interaction of the electron with the laser field becomes more important than that with the atomic potential. An advantage of this theory is that it suggests a simplified picture of the dynamics as consisting of several distinct steps (ionization, rescattering, and, in a more general context, harmonic generation), each of which can be described in the appropriate approximation. This to a large extent classical and therefore appealing picture, often referred to as the “simple man’s theory” [4–6], has gained wide popularity in the laser physics community. An overview of more recent developments of the strong-field approximation and its applications can be found in [7]. Another approach to the problem in the low-frequency regime of interest here can be based on the adiabatic approximation. Indeed, the period of oscillations of the electric field in infrared laser pulses (2.7 fs for λ = 800 nm) essentially exceeds the characteristic time of an electron’s motion in atoms 1050-2947/2010/81(3)/033415(27)

(1 a.u. ≈ 0.024 fs). Thus one can expect that the adiabatic approximation provides a suitable framework for studying the ionization dynamics. The adiabatic theory attacks the problem from a different angle and hence may shed a new light on it. Moreover, the strong-field and adiabatic approximations have different (though overlapping) regions of validity, so there are situations where the adiabatic theory should work better. In this paper, we develop the adiabatic theory for a onedimensional zero-range-potential (ZRP) model. Mathematically, the adiabatic approximation amounts to the asymptotic solution of the problem for → 0, where the adiabatic parameter is the ratio of the atomic and laser time scales. The main consequence of the adiabatic approximation for the present problem is that the solution to the TDSE can be divided into adiabatic and rescattering parts which are characterized by actions of orders O( −1 ) and O( −3 ), respectively. The adiabatic part describes the evolution of the initial bound state as it adiabatically follows a variation of the laser field. The rescattering part describes the ionized photoelectron wave packet. The successive implementation of the asymptotic approach enables us to develop a consistently quantum theory which rests on a single parameter and becomes exact in the limit → 0. This theory is uniform in terms of the amplitude of the electric field, i.e., it applies to weak as well as strong fields, provided that is sufficiently small. The tunneling regime of the strong-field approximation is recovered in this theory in the limit of weak (!) fields and not too long pulses. As far as we know, the present theory has no counterparts in the literature. The simplicity of the model considered facilitates the derivation, but its specific properties are not essential, so the general structure of the theory and its major technical elements remain the same for more realistic models. An extension of the present theory to finite-range potentials in the three-dimensional case is in progress. One of the prospect applications of few-cycle infrared laser pulses is ultrafast imaging of transient atomic species [8–14]. A possibility to accurately retrieve the structural information from laser-induced photoelectron and high-order harmonic spectra was recently demonstrated theoretically [15] and soon

033415-1

©2010 The American Physical Society

TOLSTIKHIN, MORISHITA, AND WATANABE


after that confirmed experimentally [16,17]. These results followed by a number of extensions and applications of the original theoretical procedure [18–27] have attracted much interest and opened a new direction of study [28–44]. The key contribution of the pioneering paper [15], which justifies the procedure proposed there, is in recognizing the fact that the photoelectron momentum distribution factorizes in a certain region of its argument into a product of two factors, one of which depends only on the properties of the atomic target, and the other—on the laser pulse. In [15], the factorization formula was inferred on physical grounds, quite in the spirit of the simple man’s theory, from the results of numerical calculations. Two groups have already succeeded in deriving this formula for certain situations analytically [37,41]. To address the issue from the standpoint of the adiabatic theory was the original motivation of the present work. This paper is intended to lay the groundwork for further development of the present theory. So we try to discuss (within reasonable limits of space) all the essential steps of the derivation. In Sec. II, we formulate basic equations of the problem. In Sec. III, we introduce the notion of a rescattering trajectory. Closed and ionizing rescattering trajectories form a skeleton on which the asymptotic solution of the problem is built. The next two sections present our main results: here we construct the asymptotics of the solution to the TDSE (Sec. IV) and photoelectron spectrum (Sec. V). The theory is illustrated by numerical calculations in Sec. VI. Section VII concludes the paper. Several auxiliary issues, which lie outside of the main line but are needed for the derivation, are discussed in Appendixes A–E. A preliminary report of this work was presented elsewhere [45].

The observables are defined by the coefficients in the expansion of ψ(x, t) in terms of the complete set of atomic states for t → ∞. The survival probability is given by ∞ 2 φ0 (x)ψ(x, t)dx . (5) P0 = −∞

The momentum distribution of the ionized electrons is given by [46] ∞ 2 P (±k) = ϕ± (x, k)ψ(x, t)dx , 0 k < ∞, −∞

t→∞

(6) where ϕ± (x, k), k 0, are the atomic scattering in states, eikx + R(k)e−ikx , x 0, ϕ− (x, k) = (7a) x 0, T (k)eikx , T (k)e−ikx , x 0, (7b) ϕ+ (x, k) = −ikx + R(k)eikx , x 0, e and T (k) and R(k) are the transmission and reflection amplitudes, iα k , R(k) = . (8) T (k) = k − iα k − iα The unitarity condition reads ∞ dk P0 + Pion = 1, Pion = P (k) . (9) 2π −∞ The problem consists in finding P0 and P (k) for a given pulse shape F (t). B. Adiabatic regime

II. BASIC EQUATIONS A. Formulation of the problem

We consider an electron interacting with a ZRP and a homogeneous time-dependent electric field—the simplest model to discuss the interaction of atoms with laser pulses. The TDSE in the length gauge reads (atomic units are used throughout) ∂ψ(x, t) 1 ∂2 i = − − αδ(x) + F (t)x ψ(x, t). (1) ∂t 2 ∂x 2 The coefficient of the ZRP α > 0 can be set to unity by rescaling the coordinate and time in Eq. (1), but we refrain from doing so since it is instructive to trace its appearance in the subsequent formulas. The electric field F (t) is assumed to be an analytic function of t sufficiently rapidly decaying in the remote past and future, F (t → ±∞) = 0.

(2)

The initial condition for Eq. (1) is ψ(x, t → −∞) = e−iE0 t φ0 (x),

(3)

where E0 and φ0 (x) are the energy and normalized wave function of the only bound state of the unperturbed atom, E0 = −

α2 , 2

φ0 (x) = α 1/2 e−α|x| .

t→∞

(4)

The problem contains two independent dimensionless parameters [this can be seen if one does rescale the coordinate and time in Eq. (1), as suggested above] =

1 α2T

∼ 0

1 , |E0 |T0

ξ=

F0 F0 ∼ , 3 α |E0 |3/2

(10)

where T0 and F0 are the characteristic time and magnitude of the electric field F (t), respectively. The solution to Eqs. (1) and (3) is known to have a qualitatively different structure in the different regions of the space of these parameters. Let us specify the adiabatic regime to be considered in this paper. The adiabatic parameter gives the ratio of the characteristic time of an electron’s motion in the atom to that of the variation of the electric field. It plays a key role in the following discussion. The validity of the adiabatic approximation requires 1. To uncover this parameter in Eq. (1), we assume that F (t) depends on t only via a combination τ = t dubbed “slow” time. Our goal is to construct the asymptotic solution of the problem for → 0. We do not introduce into the equations to follow explicitly by switching to slow time τ ; to avoid duplication of notation, we retain the original time variable t. However, it is essential to realize the order of the various quantities in terms of as → 0. We shall use the standard notation of asymptotic analysis [47,48]. For example, the coefficient of the ZRP does not depend on , hence α = O( 0 ) and E0 = O( 0 ). The length

033415-2

ADIABATIC THEORY OF IONIZATION OF ATOMS BY . . .


of the pulse and the intervals between any two characteristic points of the function F (t), e.g., its zeros or maxima, are O( −1 ). The magnitude of the electric field also does not depend on , i.e., F0 = O( 0 ). Then the above assumption concerning F (t) can be expressed in the form d n F (t) = O( n ). (11) dt n The parameter ξ gives the ratio of the energies of interaction of the electron with the electric field and atomic potential; the applicability of perturbation theory requires ξ 1. We have ξ = O( 0 ), which means that the limit → 0 will be taken for a fixed value of ξ ; this specifies the adiabatic regime we are interested in. In other words, the present theory is uniform in terms of ξ , i.e., it applies to weak as well as strong fields, provided that is sufficiently small. We shall see that its region of validity is defined by min(ξ, 1).

It is a property of the motion in a homogeneous electric field that any classical trajectory can be expressed in terms of the reference one, (16a) (16b)

Hence all quantities of interest can be expressed in terms of the functions v(t) and x(t). Consider the trajectory (16) in a finite interval of time ti t tf . Let u(tf ) = uf and y(tf ) = yf . One can use ti , yi , tf , and yf as independent variables identifying the trajectory. The initial ui and final uf velocities can be expressed in terms of these variables. The classical action accumulated between ti and tf is S(yf , tf ; yi , ti ) tf 1 2 u (t) − F (t)y(t) dt = 2 ti

(12)

For many-cycle almost monochromatic pulses, one can introduce the pulse frequency ω, ponderomotive energy Up = F02 /4ω2 , and the Keldysh parameter γ = |E0 |/2Up . In this case, ∼ ω/|E0 |, hence ω ∼ |E0 | = O( 1 ), Up ∼ ξ 2 |E0 |/ 2 = O( −2 ), and γ ∼ /ξ = O( 1 ). Thus the adiabatic regime → 0 corresponds to the tunneling regime γ → 0 in terms of the Keldysh theory [1]. However, we emphasize that the two limits are not equivalent, as is clear from the fact that does not depend on F0 .

u(t) = [ui − v(ti )] + v(t), y(t) = [yi − x(ti )] + [ui − v(ti )](t − ti ) + x(t).

= v(tf )yf − v(ti )yi + −

v(t) ˙ = −F (t),

˙ = v(t), x(t)

v(t → −∞) = x(t → −∞) = 0.

x(t → ∞) = x∞ + v∞ t,

where

v∞ = − x∞ =

0 −∞

(17)

(18a) (18b)

We have uf (tf , ti ) ∂ui (tf , ti ) =− , ∂tf tf − ti ∂ui (tf , ti ) ui (tf , ti ) = −F (ti ) + , ∂ti tf − ti ∂uf (tf , ti ) uf (tf , ti ) = −F (tf ) − , ∂tf tf − ti ∂uf (tf , ti ) ui (tf , ti ) = . ∂ti tf − ti

(13a)

(14a)

(19a) (19b) (19c) (19d)

The classical action for a closed trajectory is F (t)dt,

(14b)

S(tf , ti ) = S(0, tf ; 0, ti ).

−∞

∞

v(t)dt +

[v(t) − v∞ ]dt.

∂S(tf , ti ) 1 = u2i (tf , ti ), ∂ti 2 ∂S(tf , ti ) 1 = − u2f (tf , ti ). ∂tf 2

0

˙ = u(t), y(t) y(ti ) = yi .

(20)

We have (14c)

It was argued [50] that a true laser pulse must satisfy v∞ = x∞ = 0. Even if so, these conditions are not intrinsic to the present problem, so we do not impose them here; our theory is valid for arbitrary values of v∞ and x∞ . The velocity u(t) and coordinate y(t) for a general classical trajectory satisfy the same equations with different initial conditions, u(ti ) = ui ,

v 2 (t)dt.

ti

x(tf ) − x(ti ) , tf − ti x(tf ) − x(ti ) . uf (tf , ti ) = v(tf ) − tf − ti

(13b)

∞

˙ = −F (t), u(t)

tf

ui (tf , ti ) = v(ti ) −

Here and throughout, overdots denote differentiation with respect to time. Using Eq. (2), we obtain v(t → ∞) = v∞ ,

Of special interest for the present analysis are closed trajectories which begin and end at the origin. In this case, yi = yf = 0, hence such a trajectory is identified by its initial ti and final tf moments. Its initial and final velocities are given by

C. Classical dynamics in the electric field

It is well known that quantum dynamics in a homogeneous electric field can be described in purely classical terms [49]. Here we summarize formulas needed for the following. By “classical trajectory,” we mean a solution to the classical equations of motion for an electron interacting only with the electric field F (t). Let us introduce a reference classical trajectory, the velocity v(t) and coordinate x(t), defined by

1 2

[x(tf ) − x(ti ) − (yf − yi )]2 2(tf − ti )

(21a) (21b)

(15a)

Alternatively, one can use the asymptotic velocity u∞ = u(t → ∞) instead of yi or yf as one of the independent variables. Then the initial velocity depends only on u∞ and ti and is given by

(15b)

ui (u∞ , ti ) = u∞ − v∞ + v(ti ).

033415-3

(22)



The action for a closed trajectory can be represented by

III. RESCATTERING TRAJECTORIES

S(tf , ti ) = S(u∞ , ti ) − S(u∞ , tf ), where S(u∞ , t) =

1 2

t 0

(23)

u2i (u∞ , t )dt ,

(24)

and u∞ corresponds to the same trajectory. From Eq. (11) in the case of general position, we have v(t) = O( −1 ),

x(t) = O( −2 ),

(25)

S(u∞ , t) = O( −3 ).

(26)

and S(tf , ti ) = O( −3 ),

These relations establish a foundation for developing the adiabatic approximation. D. Integral form of the time-dependent Schrödinger equation

The retarded Green’s function in a time-dependent electric field is defined by ∂ 1 ∂2 i + − F (t)x G(x, t; x , t ) = δ(t − t )δ(x − x ), ∂t 2 ∂x 2 (27a) G(x, t; x , t )|t 0, we find t f − tF , ti (tf , −) ≈ tF − √ 3 2−1 tf − tF . t1 (tf , −) ≈ tF + 2 √ 3 2−1

(39)

Then its final velocity is given by uf (tf ) = uf (tf , ti (tf )) = v(tf ) − v(ti (tf )).

(40)

Using Eqs. (18) and (19), we obtain t˙i (tf ) =

−uf (tf ) . [tf − ti (tf )]F (ti (tf ))

(41)

The graphical solution of Eq. (39) is illustrated by means of the reference trajectory in Fig. 1: the straight line connecting points x(ti ) and x(tf ) must be tangential to the curve x(t) at t = ti . As is clear from the figure, for the solution to exist, the second derivative of x(t) must change sign, i.e., the electric field must vanish, at some point tF between ti and tf , F (tF ) = 0.

(42)

If both ti and tf are located in a vicinity of tF , where F (t) ≈ F˙ (tF )(t − tF ), in the leading order with respect to δi = ti − tF < 0 and δf = tf − tF > 0, we have ui (tf , ti ) ≈ 16 F˙ (tF )(δf − δi )(δf + 2δi ), uf (tf , ti ) ≈ 1 F˙ (tF )(δi − δf )(δi + 2δf ), 6

(47)

When tf passes through tR , the one-loop CRT, for tF < tf < tR , is continuously transformed into a two-loop CRT with σ1 = +, for tf > tR . From Eq. (45) in the first order with respect to tf − tR > 0, we find (48a)

(43b)

The analysis can be continued. Using Eqs. (43), it can be shown that an infinite series of reflection CRTs (all σl are minuses) with growing number of loops is created at tF . The one-loop CRT and the two-loop CRT with σ1 = − discussed above are just the first two members of the series. These CRTs correspond to the bouncing motion of the electron, which is successively reflected by the electric field, on the one side, and atomic potential, on the other. As tf departs from tF , each of these CRTs may undergo continuous transformations resulting in a change of the number of loops. This happens when either one of ul or uf turns zero. It is clear that each zero of F (t) gives rise to the appearance of a similar infinite series of CRTs. And vice versa, the same consideration suggests that each CRT as a function of its final moment tf can be traced back via a number of transformations to its creation as a reflection CRT at one of the zeros of F (t). Thus the structure of the set CRT(tf ) depends on the number of zeros of F (t) lying to the left of tf , but also on the details of the behavior of this function, and in the general case may be quite complicated. The classical action accumulated along a CRT is given by the sum of actions for each of its loops,

(44b)

tf

t1(tf ,−)

uf (tR ) = 0.

(48b)

t1(tf ,+)

tR t

tF

Thus tF is the moment of creation of a two-loop CRT with σ1 = −. In the transmission case, a solution to Eq. (45) can be found if one considers what happens with the one-loop CRT discussed above as tf departs farther from tF . At some point tf = tR the straight line connecting x(ti (tf )) and x(tf ) may become tangential to the curve x(t) also at t = tf , see Fig. 1, which means that the final velocity of the CRT vanishes too. This may happen if F (t) has another zero in the interval tF < t < tR . The point tR is defined by

ti (tf , +) ≈ ti (tR ),

x(t) ti

(46b)

t1 (tf , +) ≈ tR − (tf − tR ).

(44a)

8

(46a)

(43a)

and hence ti (tf ) ≈ tF − 12 (tf − tF ), uf (tf ) ≈ − 3 F˙ (tF )(tf − tF )2 .

(45)

tf

FIG. 1. (Color online) According to Eq. (16b), within each loop, a CRT is described by y(t) = x(t) − (a + bt). The solid line illustrates the behavior of x(t) for a model F (t). The dashed and dash-dotted lines show the linear part a + bt for one- and two-loop CRTs discussed in the text, respectively. Thin solid lines show linear parts for several higher reflection CRTs created at tF .

033415-5

S(tf , ) =

L−1 l=0

S(tl+1 , tl ).

(49)



Using Eqs. (21), it can be shown that dS(tf , ) 1 = − u2f (tf , ). dtf 2

u∞ (ti,Σ) umax ∞ tim k ti(k)

(50)

Of special interest for us are transmission CRTs (all σl are pluses); we shall see that they determine the leading-order contribution to the rescattering part of the asymptotic solution to Eq. (30). The initial moments of the different transmission CRTs from the set CRT(tf ) are given by the different solutions to Eq. (39); their final velocities are then given by Eq. (40). All loops in a transmission CRT are pieces of the same classical trajectory, hence S(tf , ) = S(tf , ti (tf )). B. Ionizing rescattering trajectories

Semi-infinite rescattering trajectories which begin at the origin with zero initial velocity and, maybe after a number of rescatterings, run away to infinity will be called ionizing rescattering trajectories (IRTs). The last rescattering event divides an IRT into a CRT and a tail—a semi-infinite classical trajectory which begins at the origin and never returns to the origin again. IRTs will be classified by the number of loops L in the CRT. The simplest zero-loop IRT is a tail with zero initial velocity. In addition to the characteristics of the CRT, an IRT is characterized by the symbol σL = ± specifying the result of the last rescattering event and the asymptotic velocity u∞ at t → ∞. It will be convenient to use different sets of variables to describe IRTs in the different situations. One can express u∞ in terms of the initial moment ti and signature = σL of the IRT, where is the signature of the CRT. For a zero-loop IRT, u∞ is defined by (we omit an empty argument ) u∞ (ti ) = −v(ti ) + v∞ .

(51)

k umin (−) ∞ ¯ as a FIG. 2. (Color online) The asymptotic velocity u∞ (ti , ) function of the initial moment ti for several simplest IRTs. An interval of ti near the same zero tF of the same model F (t) as in Fig. 1 is ¯ is empty), Eq. (51). Dashed shown. Solid line: transmission IRT ( ¯ = −), Eq. (53) with σ1 = −. Dash-dotted line: 1-reflection IRT ( ¯ consists of n minuses) with lines: several higher n-reflection IRTs ( n = 2, . . . , 5. The IRTs shown here are obtained by continuing the reflection CRTs shown in Fig. 1. The branches of u∞ (ti , −) lying to the left and right of its minimum at ti = tim correspond to the long and short IRTs discussed in Sec. V C, respectively.

¯ seen that the asymptotic velocities for IRTs with the same are given by the same functions of ti . Thus Eq. (51) defines u∞ ¯ is for any transmission IRT ( contains only pluses, hence empty); indeed, Eq. (53) for σ1 = + coincides with Eq. (51). Similarly, Eq. (53) for σ1 = − defines u∞ for any 1-reflection IRT ( contains only one minus in the first position, hence ¯ = −). In other words, u∞ depends only on ti and ; ¯ this ¯ The behavior of u∞ (ti , ) ¯ is defines a function u∞ (ti , ). illustrated in Fig. 2. This function varies in a finite interval whose boundaries will be denoted by ¯ ¯ umin ∞ () = min u∞ (ti , ), ti

(53)

Alternatively, IRTs can be identified by u∞ and . The initial moment of a zero-loop IRT is then defined by ui (u∞ , ti ) = 0 → ti = ti (u∞ ).

(54)

In the general case, the final moment of the CRT in an IRT can be found by solving Eq. (52) with respect to tf , ui (u∞ , tf ) = σL uf (tf , ) → tf = tf (u∞ , ).

¯ ¯ umax ∞ () = max u∞ (ti , ). ti

(52)

where tf is to be substituted by tf (ti , ), the inverse function to ti (tf , ). For example, for a one-loop IRT with = σ1 = ±, we have u∞ (ti , σ1 ) = σ1 [v(tf (ti )) − v(ti )] − v(tf (ti )) + v∞ .

tF

ti(k,−)

In the general case, u∞ can be found from ui (u∞ , tf ) = σL uf (tf , ) → u∞ = u∞ (ti , ),

ti

(55)

Substituting the solution into the solutions to Eqs. (35), one can express all other characteristics of the IRT in terms of u∞ and ; this defines functions ti (u∞ , ), tl (u∞ , ), ul (u∞ , ), and uf (u∞ , ). We again note that there may exist several different IRTs with the same u∞ and corresponding to the different solutions to Eqs. (35), (54), or (55). The set of all IRTs with a given asymptotic velocity u∞ will be denoted by IRT(u∞ ). This set is central to calculating the spectrum (34). ¯ denote the reduced signature of an IRT obtained from Let by omitting all pluses following the last minus. It can be

(56) ¯ will be omitted. For exFor transmission IRTs, the argument = − min v(t ) + v and umin ample, umax i ∞ ∞ ∞ (−) = min[v(ti ) − 2v(tf (ti ))] + v∞ , see Eqs. (51) and (53) and Fig. 2. The initial moment of an IRT can also be found from the equation (see Fig. 2) ¯ = u∞ → ti = ti (u∞ , ). ¯ u∞ (ti , )

(57)

This shows that ti depends only on u∞ and the reduced ¯ Of special interest for us are transmission IRTs; signature . we will see that they determine the leading-order contribution to the asymptotics of P (k). The initial moments of the different transmission IRTs from the set IRT(u∞ ) are given by the different solutions to Eq. (54). The above discussion may seem to be boring, but it introduces notions and notation needed for the following. Using this notation, two famous results in the field which apply to monochromatic pulses can be formulated as follows: (1) the maximum value of 12 u2f (tf ) as a function of tf is 3.173Up [51], and (2) the maximum value of 12 u2∞ (ti , −) as a function of ti is 10.007Up [52] (strictly speaking, the latter statement implies that the pulse is adiabatically switched on and off in such a way that v∞ = 0).

033415-6


IV. SOLUTION OF THE INTEGRAL EQUATION

t +(t) t

(a)

In this section, we construct the asymptotic solution to Eq. (30) for → 0. Let us introduce two classes of functions. Functions f (t) satisfying

A0

d n f (t) = O( ν+n ) (58) dt n will be called “slow” functions of order ν; we will write f (t) = O( ν ). For example, the electric field F (t) is a slow function of order 0, see Eq. (11). The reference velocity v(t) and coordinate x(t) are slow functions of orders −1 and −2, respectively, see Eqs. (25). The solution to Eq. (30) will be expressed in terms of “fast” functions having the form A(t) = O( ν ),

S(t) = O( µ ), (59)

where the amplitude A(t) and action S(t) are slow functions of orders ν 0 and µ < 0. We will deal with classical actions of order −3 and quantum actions of order −1; the examples are given by the exponents in Eqs. (30) and (31), respectively. The fact that g(t) has the form (59) will be denoted by g(t) ∼ (ν, µ). On physical grounds, the solution to Eq. (30) is sought in the form of a sum of adiabatic and rescattering parts, c(t) = ca (t) + cr (t),

(60)

where the two terms are identified by ca (t) ∼ (0, −1), cr (t) ∼ (ν, −3),

ν > 0,

(61a) (61b)

O(

(b)

)

tF

t A0

t'

t' A−1/2

−1

∋

g(t) = A(t) exp [iS(t)] ,

t −(t)

−1/2

∋


O( )

(c)

ti+(t) ti(t)

tF

t A2

I

t'

ti−(t) FIG. 3. (Color online) Schematic images of the different regions in the complex t plane relevant to the asymptotic solution of Eq. (30), see text. Solid circles: the adiabatic saddle point. Solid squares: tunneling and ionization saddle points.

for CRTs from the set CRT(t), there exists an SP tl+ (t, ) which accounts for the interaction with the atomic potential during the rescattering. Our strategy is to analyze the different contributions on the rhs of Eq. (30) and collect together terms having the same action. In the following analysis, we try to avoid unnecessary details, but some details are needed since the asymptotics we are going to discuss is not a standard one. For brevity, we shall use the notation t A(t ) exp[iS(t )]dt |t ∈Z , (63) −∞

and ca (t → −∞) = e−iE0 t , cr (t → −∞) = 0.

(62a) (62b)

Such a representation of the solution should and will be confirmed by the result. Taking into account the first of Eqs. (26), the integrand in Eq. (30) is a rapidly oscillating function, so the main technique of the derivation is the steepest descent method, see Appendix A. The integral evaluated using this method is given by the sum of contributions from the different saddle points (SPs) whose number and positions in the complex t plane √ depend on t, see Fig. 3. Because of the singular factor 1/ t − t , the upper limit of integration acts as a SP, actually as half of a SP because the integration path ends there. If only the contribution from this SP is retained on the rhs (throughout the paper, lhs and rhs stand for the left- and right-hand side, respectively) of Eq. (30), the equation becomes local in time, which corresponds to adiabatic evolution. Hence we call point t the adiabatic SP; it defines the evolution of the initial state that adiabatically follows a variation of the electric field F (t). We will see that there exists a pair of SPs t ± (t) located in the vicinity of t which are associated with tunneling from the initial state in a static electric field equal to the momentary value of F (t). Next, for each CRT from the set CRT(t), there exists an additional pair of SPs ti± (t, ) located near the initial moment ti (t, ) of the CRT. They are associated with ionization of the electrons which, after traveling along the CRT, return to the origin at the moment t. Finally, near each moment of rescattering tl (t, )

which means that the integral is to be evaluated using the steepest descent method and only the contributions from SPs located in zone Z are to be taken into account. A. Adiabatic part of the wave function

We seek ca (t) in the form ca (t) = Ca (t) exp −iE0 t − i

t

−∞

[E0 (t ) − E0 ]dt

, (64)

where E0 (t) and Ca (t) are slow functions of order 0. The action in Eq. (64) is a slow function of order −1, in accordance with Eq. (61a). To satisfy the initial condition (62a) we require E0 (t → −∞) = E0 ,

Ca (t → −∞) = 1.

(65)

A vicinity of t in the complex t plane defined by zone Aν : |t − t| = O( ν )

(66)

will be called the adiabatic zone. Let us substitute Eq. (64) for c(t) into the rhs of Eq. (30) and consider the contribution to the integral from zone A0 , see Fig. 3(a). Since the exponent in Eq. (30) vanishes at t = t, this contribution is proportional to ca (t), with the coefficient being a slow function. We will see that all other contributions to the integral have actions of order −3. Hence retaining only the specified term on the rhs of Eq. (30), one has to retain only the term ca (t) on the lhs. We

033415-7



thus arrive at the equation eiπ/4 α t Ca (t )dt Ca (t) = √ exp[iSa (t, t )] √ , t − t t ∈A0 2π −∞ where

Sa (t, t ) = S(t, t ) +

(67)

E0 (t )dt .

(68)

Equation (67) defines functions Ca (t) and E0 (t). Below we construct the asymptotic solution to this equation. A simple way to proceed is to find all SPs located in zone A0 and calculate their individual contributions to the integral in Eq. (67). The adiabatic SP t is one of them; we have Sa (t, t) = 0. The other SPs are defined by

∂Sa (t, t ) 1 = u2i (t, t ) − E0 (t ) = 0. (69) ∂t 2 This equation has two solutions t + (t) and t − (t) in the vicinity of t lying in the upper and lower half-plane, respectively; we shall call them the tunneling SPs. Using Eqs. (18) and (19), we find 2iα(t) + O( 1 ), (70) t ± (t) = t ± |F (t)| Re α(t) > 0.

(77)

2iα 3 (t) + O( 1 ). 3|F (t)|

(72)

The integration path in Eq. (67) can be deformed into a steepest descent contour passing through t + (t) and t, but not t − (t); hence only these two SPs contribute to the integral. Their contributions can be treated separately if |Sa (t, t + (t)) − Sa (t, t)| 1. As follows from Eq. (72), this condition is violated for sufficiently large F (t). The overall strength of the field is characterized by the parameter ξ , see Eqs. (10). Thus, in order to obtain the asymptotic solution to Eq. (67) which is uniform in terms of ξ, the complex consisting of the three SPs in zone A0 [see Fig. 3(a)] must be treated as a whole. Such a uniform treatment is one of the achievements of this work. Expanding action (68) in zone A0 in powers of , we obtain 1 F 2 (t)δ 3 24 4 ˙

− 12 E˙ 0 (t)δ 2

1 + 24 F (t)F (t)δ + O( 2 ),

(73)

where δ = t − t . Similarly we have Ca (t ) = Ca (t) − C˙ a (t)δ + O( 2 ).

(74)

The terms containing slow functions F (t), Ca (t), and E0 (t) and their derivatives in these expansions have orders 0 and 1, respectively. Substituting Eqs. (73) and (74) into Eq. (67) and collecting together terms of the same order in , in the zeroth order we obtain dδ eiπ/4 α ∞ i 1= √ exp iE0 (t)δ − F 2 (t)δ 3 1/2 . (75) 24 δ 2π 0

(78)

where E0 (F ) is the Siegert eigenvalue which coincides with E0 for F = 0, see Appendix C. To find Ca (t), one has to proceed to the next order of the expansion. In the first order, we obtain d ˙ 2Ca (t) + Ca (t) N (t) = 0, (79) dt where

i 2 3 exp iE0 (t)δ − F (t)δ δ 1/2 dδ 24 0 (80a) 2 4π α dW± (z) . (80b) = −e±iπ/6 F (t) dz z=z0 (t)

e3iπ/4 α 2 N (t) = √ 2π

(71)

If t is not too close to a zero of F (t), points t (t) indeed lie at a distance O( 0 ) from t, which explains the size of the adiabatic zone A0 included in Eq. (67). We have

Sa (t, t ) = E0 (t)δ −

z0 (t) = −21/3 [F (t)]−2/3 E0 (t).

E0 (t) = E0 (F (t)),

±

Sa (t,t ± (t)) = ±

(76)

The upper and lower signs in Eq. (76) stand for the analytic continuation from positive [arg F (t) = 0] and negative [arg F (t) = −π ] values of F (t) on the real t axis, respectively. Comparing Eq. (76) with Eq. (C3) and selecting the solution satisfying Eqs. (65), we find

where α(t) = [−2E0 (t)]1/2 ,

e∓iπ/6 [2F (t)]1/3 , 4π α

W± (z0 (t)) = where

t t

This equation defines E0 (t). Using Eqs. (B2b) and (B3b), it can be cast into the form

∞

We have N (t → −∞) = α −1 , hence the solution to Eq. (79) satisfying Eqs. (65) is Ca (t) = [αN (t)]−1/2 . Comparing this with Eqs. (C7) and (C8), we find Ca (t) = α −1/2 φ0 (0; F (t)),

(81)

where φ0 (x; F ) is the corresponding Siegert eigenfunction. The procedure can be continued; taking into account the next term in the expansion, one can obtain a correction to Ca (t) of order 1, etc. The structure of this asymptotic series resembles that of the semiclassical solution to the stationary Schrödinger equation [53], but the derivation is different. The above derivation is based on Eq. (73): it was assumed that the dotted terms in this expansion are much smaller than the first two terms near the tunneling SPs t ± (t), i.e., for |δ| ∼ |α(t)/F (t)|. This local condition is certainly satisfied if the following global one holds: E0 1/2 . ξ (82) E0 (F0 ) For weak and moderate fields, ξ 1, we have E0 (F0 ) ∼ E0 , so Eq. (82) reduces to ξ . For strong fields, ξ 1, the magnitude of the rhs of Eq. (82) depends on the behavior of E0 (F ) for large values of F . If E0 (F ) grows not faster than F 2 [the meaning of this condition is that the SPs t ± (t) do not recede from the real axis], as in the present model, then the rhs of Eq. (82) becomes larger than 1 for sufficiently large ξ . The region of validity of the present theory can thus be specified by Eq. (12). Summarizing, the asymptotics of the adiabatic term in Eq. (60) for → 0 is given by Eq. (64), where functions E0 (t) and Ca (t) are expressed in terms of the Siegert state by

033415-8



Eqs. (78) and (81). Thus in the presence of an external electric field, the initial bound atomic state turns into a properly normalized Siegert state which adiabatically follows a variation of the field. The integral term in the exponent in Eq. (64) accounts for a difference between the eigenvalues. It describes the accumulation of an additional phase and decay of the initial state due to its interaction with the field. For weak fields, this can be interpreted in well-known terms of the secondorder Stark shift and semiclassical tunneling probability, see Eqs. (C10) and (C11). The exact Siegert eigenvalue E0 (F ) arising in the above analysis makes the result (64) applicable to fields of arbitrary strength. The preexponential factor in Eq. (64) accounts for a difference between the eigenfunctions; for weak fields, it is equal to unity, see Eqs. (81) and (C12). The asymptotics (64) agrees with Eq. (61a). We emphasize that it is uniform in terms of the parameter ξ . B. Formation zone

Meanwhile, it is not evident whether the asymptotics (64) is uniform in terms of t. Indeed, as t approaches a zero tF of F (t), the tunneling SPs t ± (t) depart from the adiabatic one t, move deeper into the complex plane, and eventually leave zone A0 , see Eq. (70) and Fig. 3(a). The analysis of the previous section based on expansions (73) and (74) does not apply in this case. The validity of Eq. (64) for t located in a vicinity of tF requires a special consideration. Let t belong to the formation zone defined by zone F : |t − tF | = O(

−1/2

).

F˙ 2 (tF ) (t − tF )2 (t − t )3 = O( 1 ). (87) 24 Substituting this into Eq. (84) and retaining on the rhs only terms of orders 0 and 1, we obtain 5F˙ 2 (tF ) i 2 ca (t) = 1 + (t − tF ) ca (t) + 2 c˙a (t). (88) 8α 6 α S(t, t ) ≈ −

The solution to this equation is given by ˙2 5i F (tF ) 3 (t − tF ) , ca (t) = ca (tF ) exp 24α 4

The action in this equation can be expanded in powers of , F˙ 2 (tF ) (δ − δ )3 (4δ 2 + 7δδ + 4δ 2 ) + O( 0 ), 360 (85)

where δ = t − tF = O( −1/2 ) and δ = t − tF = O( −1/2 ), so the leading-order term is O( −1/2 ). Consider the SPs contributing to the integral in Eq. (84). We again have the adiabatic SP t. The other SPs are defined by F˙ 2 (tF ) (δ − δ )2 (δ + 2δ )2 − E0 = O( 1/2 ). (86) 72 This equation (with respect to t ) has four solutions lying at a distance O( −1/2 ) from t, see Fig. 3(b); this explains the size of the adiabatic zone A−1/2 included in Eq. (84). The left pair of solutions coincides with what is obtained from the tunneling SPs t ± (t) as t approaches zone F from the left. The right pair turns into the tunneling SPs t ± (t) as t leaves zone F moving to the right. The analysis of the previous section implies that the distance between the two tunneling SPs is much smaller than that between the left and right pairs of the solutions to Eq. (86), which is not the case when t is inside zone F . To extend the uniform approach of the previous section to zone F, one would have to treat the complex consisting of all five SPs in zone A−1/2 as a whole, see Fig. 3(b).

t ∈ F. (89)

This result applies throughout zone F . On the other hand, in this zone, F (t) ≈ F˙ (tF )(t − tF ) = O( 1/2 ). Hence the energy shift of the adiabatic Siegert state is E0 (t) − E0 ≈ −

(83)

Repeating the argumentation that has led us to Eq. (67), we conclude that ca (t) in zone F satisfies eiπ/4 α t ca (t )dt exp[iS(t, t )] √ . (84) ca (t) = √ t − t t ∈A−1/2 2π −∞

S(t, t ) = −

Fortunately, there are two circumstances which greatly simplify the analysis. First, the difference between the values of the action in Eq. (84) at the different SPs is O( −1/2 ). Therefore all the SPs can be treated separately. Second, the contributions to the integral from all but the adiabatic SP are ∝ exp (−const × −1/2 ), where const > 0. Such exponentially small (in terms of the parameter ) terms are not accounted for by the present analysis. Therefore only the contribution from the adiabatic SP is to be retained, which amounts to reducing the adiabatic zone included in Eq. (84) from A−1/2 to A0 , see Fig. 3(b). In this zone, we have

5F˙ 2 (tF ) (t − tF )2 = O( 1 ), 8α 4

(90)

and its width is ∝ exp(−const × −1/2 ), see Eqs. (C10) and (C11). Substituting Eq. (90) into Eq. (64), one finds that Eqs. (64) and (89) coincide in zone F . Thus the asymptotics (64) remains valid throughout zone F and hence is uniform in terms of t. As t approaches zone F from the left, the contribution from the left pair of SPs in Fig. 3(b) to the integral in Eq. (84) describes decay of the initial state via tunneling. The decay rate is exponentially small in zone F , so we have neglected it in the above analysis. As t leaves zone F moving to the right, the same is true for the right pair of SPs in Fig. 3(b). The left pair in this case turns into the ionization SPs ti± (t), see the next section. Their contribution to Eq. (84) is a seed from which the rescattering term in Eq. (60) grows. Thus zone F is where the formation of the rescattering term takes place, which explains the terminology. C. Rescattering part of the wave function

Substituting ca (t) for c(t) into the rhs of Eq. (30) and calculating the contribution to the integral from the adiabatic zone, one obtains a term with action of order −1. This term is compensated by the adiabatic part ca (t) on the lhs of the equation, see Sec. IV A. In this section, we show that if F (t ) has zeros in the interval −∞ < t < t, then the asymptotics of the integral in Eq. (30) contains terms with actions of order −3. This leads to the appearance of the rescattering part cr (t) on the lhs of the equation. The different terms in the asymptotics of cr (t) for → 0 are associated with the different CRTs from the set CRT(t). Following the logic of Sec. III A, we first discuss the contribution to cr (t) from the simplest one-loop CRT. We

033415-9



then proceed to contributions associated with two-loop CRTs. This will suffice to grasp the structure of cr (t) in the general case.

and

∂ 2 S1 (t, t ) ∂t 2 t =ti±

1. One-loop CRT

Let tF be a zero of F (t), see Eq. (42). Consider Eq. (30) for t > tF assuming that t − tF = O(

−1

),

(91)

which is the case of general position in the interval under consideration. The action obtained upon substituting Eq. (64) for c(t) into the rhs of Eq. (30) is t S1 (t, t ) = S(t, t ) − E0 t − [E0 (t ) − E0 ]dt . (92) −∞

The SPs for this action are defined by ∂S1 (t, t ) 1 = u2i (t, t ) − E0 (t ) = 0. (93) ∂t 2 This equation has a simple physical meaning. In the adiabatic regime, a transition can efficiently occur only when the energies of the initial and final states coincide; transitions associated with a change of the energy of the system are suppressed. In the initial state, the electron is bound to the atom and its energy E0 (t ) adiabatically follows the momentary value of the electric field F (t ). In the final state, the electron is ionized and moves under the influence of the electric field. In order to return to the origin and experience rescattering at the moment t, its energy at the moment of ionization t must be equal to 12 u2i (t, t ). Thus Eq. (93) defines the moment of ionization in the adiabatic regime. One can notice that Eq. (93) coincides with Eq. (69); however, now we are interested in the solutions lying in a different region of the complex t plane. The two terms in Eq. (93) have orders −2 and 0, respectively, hence the equation can be satisfied only near a zero of ui (t, t ). Suppose there exists a one-loop CRT born at tF which ends at t, see Sec. III A. Then ui (t, t ) turns zero at the initial moment t = ti (t) of the CRT, see Eq. (39). From Eq. (91) we have t − ti (t) = O( −1 ). Consider the ionization zone defined by zone I : |t − ti (t)| = O( 0 ).

(94)

Equation (93) has two solutions ti± (t) in this zone to be called the ionization SPs, see Fig. 3(c). They coincide with what is obtained from the left pair of SPs shown in Fig. 3(b) as t leaves zone F moving to the right. Introduce notation [compare with Eq. (40)] ± u± f (t) = uf (t, ti (t)).

(95)

To simplify the equations, in the rest of this section, we omit the argument in ti (t), uf (t), ti± (t), and u± f (t). At the SPs we have ui (t, ti± ) = ∓isi α(ti± ),

(96)

where si is defined by Eq. (38). Using Eqs. (18) and (19) we find iα(ti ) ti± = ti ± + O( 1 ), (97) |F (ti )|

iα 3 (ti ) + O( 1 ), 3|F (ti )| ± −1 ∓isi α(ti± )u± dti f = ± dt t − ti

S1 (t, ti± ) = S1 (t, ti ) ±

= ±iα(ti )|F (ti )| + O( 1 ).

(98a) (98b) (98c)

Here S1 (t, ti ) = O( −3 ); the second terms on the rhs of Eqs. (97) and (98a) and the first term on the rhs of Eq. (98c) have order 0. Let cr(a) (t) denote the contribution to the rhs of Eq. (30), where c(t) is substituted by ca (t), from zone I : eiπ/4 α t (a) Ca (t )dt cr (t) = √ exp[iS1 (t, t )] √ . (99) t − t t ∈I 2π −∞ Introduce a function 2/3 3i − + [S1 (t, ti ) − S1 (t, ti )] ζ (t) = 4 =

α 2 (ti ) + O( 1 ). [2|F (ti )|]2/3

(100a) (100b)

Since |ζ (t)| may become small for sufficiently strong fields, for obtaining the asymptotics of cr(a) (t) which is uniform in terms of ξ, the complex consisting of the two SPs ti± must be treated as a whole. Evaluating the integral in Eq. (99) we find (see Appendix A) +

(101a) cr(a) (t) = Cr(a) (t)eiS1 (t,ti )

1/2 + −dti+ /dt eiS(t,ti )+iπ/4 ca (ti+ ), = αU (ζ (t)) + + si α(ti )uf (101b) where Cr(a) (t)

= . (t − t )∂ 2 S1 (t, t )/∂t 2 t =t + iαU (ζ (t))Ca (t )

(102)

i

Cr(a) (t)

S1 (t, ti+ )

The amplitude and action in Eq. (101a) are slow functions of orders 1/2 and −3, respectively. Hence cr(a) (t) gives a contribution to cr (t); the superscript in the notation indicates that this contribution comes directly from the adiabatic term in Eq. (60). A key moment in the derivation consists in recognizing the fact that there exist other terms on the rhs of Eq. (30) with the same action S1 (t, ti+ ). Indeed, let us substitute cr(a) (t) for c(t) into the rhs of Eq. (30) and calculate the contribution to the integral from the adiabatic SP. Since the exponent in Eq. (30) vanishes at t = t, the result has the same action as in Eq. (101a), but a different amplitude. Substituting it for c(t) into the rhs of Eq. (30) and again calculating the contribution from the adiabatic SP, one obtains yet another term with the same action, etc. The sum of the series of all such terms on the rhs of Eq. (30) will be denoted by c1 (t); the subscript indicates that this is a contribution to cr (t) associated with a one-loop CRT. This function satisfies the inhomogeneous equation eiπ/4 α t c1 (t )dt c1 (t) = cr(a) (t) + √ exp[iS(t, t )] √ . t − t t ∈A2 2π −∞ (103)

033415-10



We seek the solution in the form c1 (t) = C1 (t)e

iS1 (t,ti+ )

(104)

,

where C1 (t) is a slow amplitude. Since t = Eq. (93), we have

ti+

dS1 (t, ti+ ) 1 = − (u+ )2 . dt 2 f

is a root of

(105)

The evaluation of the integral in Eq. (103) depends on the −1 value of u+ f . For t − ti = O( ) we find u+ f = uf +

α 2 (ti ) + O( 2 ), 2(t − ti )F (ti )

(106)

where uf = O( −1 ) and the second term has order 1. This explains the size of the adiabatic zone included in Eq. (103). In this case, only the adiabatic SP t contributes to the integral. Calculating its contribution, we obtain an algebraic equation for C1 (t), C1 (t) = Cr(a) (t) +

iα C1 (t), u˜ + f

(107)

where + u˜ + f = sgn(uf )uf .

(108)

+ For a one-loop CRT u˜ + f = si uf , see Eq. (37). Solving Eq. (107), we finally obtain (a) c1 (t) = T (u˜ + f )cr (t),

(109)

where T (k) is the transmission amplitude defined by Eqs. (8). It remains to specify the interval of t where Eq. (109) holds. In the derivation, we assumed that t − tF is not too small, see 1 Eq. (91). In this case uf = O( −1 ), T (u˜ + f ) = 1 + O( ), and one can use expansions (97), (98a), (98c), (100b), and (106). This is the case of general position. Meanwhile, it can be seen that Eq. (109) remains valid in a wider interval of t. Its region of validity is restricted by our method of solving Eq. (103) and is defined by the condition uf O( 0 ). Thus Eq. (109) holds as t approaches tF up to the right boundary of the formation zone, where t − tF = O( −1/2 ) and uf = O( 0 ). For this to be the case, however, all factors in Eqs. (101) and (109) must be calculated as prescribed, without further approximations. This explains why the transmission amplitude T (u˜ + f ) in Eq. (109) cannot be replaced by unity. For the same reason, one should not rely on the expansions mentioned above. Indeed, the first two terms on the rhs of Eq. (106) become of the same order 0 for t − tF = O( −1/2 ), and similarly with the other expansions. On the other hand, as t departs from tF , Eq. (109) holds up to the left boundary of the recombination zone (see Sec. IV D), where uf turns zero and the one-loop CRT disappears. Let us discuss these results. There is a situation where they become especially transparent and can be interpreted in the spirit of the simple man’s theory. This is the case for weak fields, ξ 1, in the interval of t where uf = O( −1 ). First of all, it is instructive to consider an approximate form of cr(a) (t) under these conditions. In this case, α(ti+ ) ≈ α(ti ) ≈ α, U (ζ (t)) ≈ 1, and expansions (97), (98a), (98c), and (106) can

be used. From Eq. (101b) using Eq. (C11), we obtain 0 (F (ti )) dti 1/2 iS(t,t )+iπ/4 (a) i e ca (ti ). cr (t) ≈ αuf dt

(110)

Each factor in this formula has a simple physical meaning. The electron stays in the bound state until the moment ti . Then it is ionized, travels along a one-loop CRT, and returns to the origin at the moment t. The amplitude to survive in the bound state is ca (ti ). The probability to be ionized via tunneling during the interval dti is 0 (F (ti ))|dti |. The classical action accumulated along the CRT is S(t, ti ). The velocity of the returning electron is uf . The function cr(a) (t) has the meaning of the amplitude of the returning wave packet divided by α 1/2 , see Eq. (29). The electrons ionized during the interval dti are rescattered during the interval dt. Equation (110) ensures that the number of the ionized electrons |ca (ti )|2 0 (F (ti ))|dti | is equal to the number of electrons arriving for rescattering α|cr(a) (t)|2 |uf dt|, as it should be. An additional phase π/4 in Eq. (110) is accumulated during the tunneling. Equation (101b) is less transparent: because of the Airy correcting factor U (ζ (t)), because of a difference between the Siegert eigenvalue E0 (t) and the unperturbed energy of the initial state E0 , and that between a quantum object, the SP ti+ , and its classical counterpart, the initial point of the CRT ti . But it describes the same wave packet and remains valid for arbitrary values of ξ and in a wider interval of t. Turning to Eq. (103), this is a kind of Lippman-Schwinger equation. The source term cr(a) (t) describes the ionized wave packet returning for rescattering. The difference between c1 (t) and cr(a) (t) accounts for the interaction with the atomic potential during the rescattering. The series of terms mentioned above Eq. (103) is the Born series, i.e., the perturbative expansion of the solution in terms of the parameter α/uf . The present method of solving Eq. (103) is based on the adiabatic instead of the Born approximation. This amounts to summing up the whole Born series and results in the exact transmission amplitude T (u˜ + f ) in Eq. (109). Summarizing, a one-loop CRT ending at t gives rise to the appearance in the asymptotics of cr (t) a term c1 (t) given by Eq. (109). This asymptotics is uniform in terms of ξ and applies in the interval of t where uf O( 0 ). In the case of general position, uf = O( −1 ), this result agrees with Eq. (61b), where ν = 1/2. It is clear that each one-loop CRT from the set CRT(t), which is created at tF or one of the earlier zeros of F (t), produces a similar contribution to cr (t) given by the same formulas. 2. Two-loop CRT

We proceed to the discussion of contributions to cr (t) associated with two-loop CRTs. Substituting c1 (t) for c(t) into the rhs of Eq. (30), one obtains the action S2 (t, t ) = S(t, t ) + S1 (t , ti+ (t )).

(111)

Using Eqs. (21) and (105), the SPs for this action are defined by 2 u2i (t, t ) − [u+ f (t )] = 0.

(112)

We can rewrite this equation in the form

033415-11

+ ui (t, t ) = σ1 u+ f (t ) → t = t1 (t, σ1 ),

(113)



where σ1 = ±. The solution t1+ (t, σ1 ) can be found near a solution to Eq. (45), t1+ (t, σ1 ) = t1 (t, σ1 ) + O( 1 ).

(114)

Here t1 (t, σ1 ) is the moment of rescattering for a twoloop CRT which ends at t and has the signature σ1 , see Sec. III A, therefore we call t1+ (t, σ1 ) the rescattering SP. Introduce notation ti± (t, σ1 ) = ti± (t1+ (t, σ1 )), u+ 1 (t, σ1 ) = + uf (t1+ (t, σ1 ), ti+ (t, σ1 )), and u+ (t, σ ) = u (t, t (t, σ1 )). In 1 f f 1 the rest of this section, we omit the arguments in ti (t, σ1 ), t1 (t, σ1 ), u1 (t, σ1 ), uf (t, σ1 ), ti± (t, σ1 ), t1+ (t, σ1 ), u+ 1 (t, σ1 ), and (1) u+ f (t, σ1 ). Let cr (t, σ1 ) denote the contribution to the rhs of Eq. (30), where c(t) is substituted by c1 (t), from the SP t1+ ; the superscript indicates that this contribution comes directly from the one-loop term c1 (t). This function has the action S2 (t, t1+ ). One can observe that there are other terms on the rhs of Eq. (30) with the same action. Similar to the one-loop case, an infinite series of such terms arises from the adiabatic SP. In the transmission case, the above action is equal to S1 (t, ti+ ), because two loops in the CRT are parts of the same classical trajectory, and there exists an additional term of the type cr(a) (t) with the same action. The sum of all such terms on the rhs of Eq. (30) will be denoted by c2 (t, σ1 ); this is a contribution to cr (t) associated with a two-loop CRT. This function satisfies an integral equation of the type (103) with the source term given by cr(a) (t) + cr(1) (t, +), in the transmission case, and cr(1) (t, −), in the reflection case. Solving this equation in the adiabatic approximation by the method used to solve Eq. (103), one can find c2 (t, σ1 ). We omit further details of the derivation. The result reads

1/2 −dti+ /dt + + + c2 (t, σ1 ) = T (u˜ f )A(u˜ 1 )αU (ζ (t1 )) si σ1 α(ti+ )u+ f +

+

+

× eiS(t,t1 )+iS(t1 ,ti

)+iπ/4

ca (ti+ ),

transparency of the result, we give only an approximate form of cL (t, ) in the general case: cL (t, ) ≈ T (|uf |)A(|uL−1 |) · · · A(|u1 |) 0 (F (ti )) dti 1/2 iS(t, )+iπ/4 e × ca (ti ). (117) αuf dt Here ti , ul , and uf are the characteristics of the CRT as functions of its final moment t and signature , see Sec. III A, and A is either T or R, depending on whether the corresponding rescattering event results in transmission or reflection, respectively. Equation (117) applies in the situation mentioned above Eq. (110): it is valid for ξ 1 in the interval of t where ul = O( −1 ) and uf = O( −1 ). The scattering amplitudes in Eq. (117) account for the interaction with the atomic potential at the moments of rescattering and the final moment of the CRT, and the other factors can be interpreted similarly to Eq. (110). The individual contributions to cr (t) given by Eq. (117) satisfy Eq. (61b) with ν = 1/2 + r , where r is the number of reflections in the CRT. Thus the leadingorder term in the asymptotics of cr (t) for → 0 is determined by transmission CRTs. We emphasize that Eq. (117) is only an approximate form: the approach developed above enables one to obtain the asymptotics of cL (t, ) which is uniform in terms of ξ and applies in a wider interval of t where ul O( 0 ) and uf O( 0 ), see Eqs. (109) and (115). Each term in Eq. (116) is localized in the interval of t where the corresponding CRT exists, i.e., where ti and tl defined by Eqs. (35) for tf = t are real. This is because only in this interval are the SPs contributing to cL (t, ) located close to the real axis, see Eqs. (97) and (114). When t leaves this interval, all the SPs move deep into the complex time plane, and the magnitude of cL (t, ) rapidly decays.

(115)

where the symbol A stands for T or R in the transmission (σ1 = +) and reflection (σ1 = −) cases, respectively, u˜ + f is + = −s σ defined by Eq. (108) [for a two-loop CRT u˜ + i 1 uf , f + + see Eq. (37)], and similarly u˜ + 1 = sgn(u1 )u1 = si u1 . This asymptotics is uniform in terms of ξ and applies in the interval of t where u1 O( 0 ) and uf O( 0 ). In the case of general position, u1 = O( −1 ) and uf = O( −1 ), hence T = 1 + O( 1 ) and R = O( 1 ). Then c2 (t, +) ∼ (1/2, −3) and c2 (t, −) ∼ (3/2, −3), which agrees with Eq. (61b). Each two-loop CRT from the set CRT(t) produces a contribution to cr (t) given by the same formula (115). 3. General structure of cr (t)

The analysis can be continued. In a similar way, one can derive formulas describing contributions to cr (t) associated with higher CRTs having larger number of loops. All the essential moments of the derivation have been illustrated above. The final result reads cr (t) = cL (t, ), (116) CRT(t)

where the summation runs over all CRTs from the set CRT(t). The one- and two-loop terms in Eq. (116) are given by Eqs. (109) and (115), respectively. For brevity and

D. Recombination zone

As can be seen from the analysis in Sec. IV C, the asymptotics of cr (t) given by Eq. (116) is not uniform in terms of t. Indeed, Eq. (116) does not apply near zeros of F (t). The interaction between an infinite series of terms in Eq. (116) corresponding to reflection CRTs created at tF should be taken into account [which amounts to treating the complex consisting of all five SPs shown in Fig. 3(b) as a whole] to construct a uniform asymptotics of cr (t) in the formation zone. In addition, the asymptotics of the individual terms in Eq. (116) constructed in Sec. IV C does not apply in small intervals of t around the recombination points, where one of ul or uf turns zero. At such points, a change of the number of loops in the CRT occurs. The corresponding term in Eq. (116) should change appropriately as t passes through a recombination point, which is not accounted for by the above analysis. Technically, the lack of uniformity in the latter case results from our method of solving integral equations of the type Eq. (103) defining cL (t, ), whose applicability requires ul O( 0 ) and uf O( 0 ). Physically, vanishing ul or uf is associated with the phenomenon of recombination, which explains the terminology. Let us discuss the simplest situation illustrated in Fig. 1. The final velocity uf of the one-loop CRT born at tF vanishes at tR , see Eq. (47). Consider the recombination zone

033415-12



defined by zone R:

|t − tR | = O( ). 0

(118)

In this zone 1 u+ f (t) = −F (tR )(t − tR ) + O( ).

(119)

Hence [see Eq. (105)]

process to tunneling ionization. It is localized in zone R because only slow electrons can be efficiently captured by the atomic potential in the presence of the electric field. We mention that the other poles in Eq. (124) also contribute to the integral, which corresponds to population of the other Siegert states as t passes through zone R. A more detailed analysis of this recombination dynamics is postponed to future studies.

S1 (t, ti+ (t)) = S1 (tR , ti+ (tR )) − 16 F 2 (tR )(t − tR )3 + O( 1 ), (120) and the inhomogeneous term in Eq. (103) is given by

where

cr(a) (t) = χ (t − tR )[1 + O( 1 )],

(121)

i χ (t) = cr(a) (tR ) exp − FR2 t 3 , 6

(122)

and FR = |F (tR )|. Substituting Eq. (121) into Eq. (103) and expanding the action in the integral term in powers of , in the leading order, Eq. (103) in zone R takes the form

V. CALCULATION OF THE PHOTOELECTRON SPECTRUM

Having constructed the asymptotic solution to Eq. (30), numerical integration is one option for calculating the photoelectron spectrum (34). However, it is more consistent with the present theory to calculate the integral in Eq. (34) using the steepest descent method. The latter approach does not introduce any additional approximations and is justified by the second of Eqs. (26). The spectrum (34) can be presented in the form P (k) = 2π |I (k)|2 ,

iπ/4

e α c1 (t) = χ (t − tR ) + √ 2π t i 2 3 c1 (t )dt exp − FR (t − t ) √ . × 24 t − t −∞

(123)

This equation can be solved using the Fourier transformation. The exact solution is ∞ dE , (124) c1 (t) = T (k; FR )χ (E)e−iE(t−tR ) 2π −∞ √ where k = 2(E + i0), T (k; F ) is the transmission amplitude in the presence of the electric field, see Appendix D, and χ (E) is the Fourier transform of χ (t) given by −2/3 −2/3 χ (E) = cr(a) (tR )24/3 π FR Ai − 21/3 FR E . (125) −2/3

Equation (124) applies in zone R. For |t − tR | FR , the integral in Eq. (124) can be evaluated using the steepest descent method. For t < tR , the steepest descent contour passes through only one SP located at E = 12 FR2 (t − tR )2 . Calculating its contribution, it can be shown that Eq. (124) coincides with Eq. (109) at the left boundary of zone R. Thus Eq. (124) enables one to continue the asymptotics (109) through zone R. Note that T (k; F ) has poles in the complex E plane at the energies of the Siegert states, see Appendix D. As t passes through tR , the steepest descent contour in Eq. (124) sweeps over the pole at E = E0 (FR ), so that the asymptotics of Eq. (124) at the right boundary of zone R in addition to the contributions from SPs contains a term arising from the pole. This term depends on time ∝ exp[−iE0 (FR )(t − tR )], i.e., is proportional to ca (t), see Eq. (64). It should be combined with ca (t) to the right of zone R, which results in a steplike change of the amplitude of ca (t) as t passes through zone R, ca (t > tR ) = (1 + R )ca (t < tR ),

(126)

where R is a constant that can be calculated from Eq. (124) using Eq. (D5). The appearance of the term R in Eq. (126) accounts for partial recombination of the electrons which arrive for rescattering and are represented by the source term χ (t − tR ) in Eq. (123). Such a recombination is the inverse

I (k) = Ia (k) + Ir (k),

(127)

where the two terms in the ionization amplitude I (k) correspond to the two terms in Eq. (60). In this section we discuss the evaluation of the integrals Ia (k) and Ir (k). The technique of the derivation is similar to that in Sec. IV. The different SPs contributing to I (k) are associated with the different IRTs from the set IRT(k). Their number and positions in the complex t plane depend on k. Our strategy is again to analyze the individual contributions and collect together terms having the same action. A. Adiabatic part of the ionization amplitude

Substituting ca (t) given by Eq. (64) for c(t) into Eq. (34), one obtains the action t S0 (k, t) = S(k, t) − E0 t − [E0 (t ) − E0 ]dt . (128) −∞

The SPs for this action are defined by ∂S0 (k, t) 1 = u2i (k, t) − E0 (t) = 0. ∂t 2

(129)

This equation defines the moment of ionization t in the adiabatic regime, see the argumentation after Eq. (93), now as a function of the asymptotic velocity k of the ionized electron. The two terms in Eq. (129) have orders −2 and 0, respectively, therefore the solution is to be sought near a zero of ui (k, t). Such a zero is the initial moment ti (k) of a transmission IRT with the asymptotic velocity k, see Eq. (54) and Fig. 2. Equation (129) has two solutions near ti (k) to be denoted by ti± (k) and called the ionization SPs. In the rest of this section, we omit the argument in ti (k) and ti± (k). At the SPs we have [compare with Eq. (96)] ui (k, ti± ) = ∓isi α(ti± ).

(130)

Similar to Eqs. (97) and (98), we find

033415-13

ti± = ti ±

iα(ti ) + O( 1 ), |F (ti )|

(131)



and

is iα (ti ) + O( 1 ), 3|F (ti )| ± −1 dti = ±isi α(ti± ) dk

S0 (k, ti± ) = S0 (k, ti ) ± ∂ 2 S0 (k, t) ∂t 2 t=ti±

S1 (k, t) = S(k, t) + S1 (t, ti+ (t)).

3

= ±iα(ti )|F (ti )| + O( 1 ). Introduce a function 2/3 3i ζ (k) = [S0 (k, ti− ) − S0 (k, ti+ )] 4 α 2 (ti ) + O( 1 ). = [2|F (ti )|]2/3

(132a)

(135)

Using Eq. (105), the SPs are defined by 2 u2i (k, t) − [u+ f (t)] = 0.

(132b)

(136)

We can rewrite this equation in the form (132c)

(133a)

+ ui (k, t) = σ1 u+ f (t) → t = tf (k, σ1 ),

where σ1 = ±. The solution tf+ (k, σ1 ) can be found near a solution to Eq. (55), tf+ (k, σ1 ) = tf (k, σ1 ) + O( 1 ).

(133b)

Since |ζ (k)| may become small for sufficiently strong fields, for obtaining the asymptotics of Ia (k) which is uniform in terms of ξ, the two SPs ti± must be treated simultaneously. We thus obtain (see Appendix A) + 1/2 + dti /dk Ia (k) = α 3/2 U (ζ (k)) eiS(k,ti ) ca (ti+ ). (134) + si α(ti ) For each transmission IRT from the set IRT(k), there exists a pair of ionization SPs ti± that produce a contribution to Ia (k) given by Eq. (134). The total amplitude Ia (k) is the sum of all such contributions. As is clear from the derivation, the adiabatic part Ia (k) of the ionization amplitude takes into account only transmission IRTs and only in the approximation that transmission amplitudes for all rescattering events are equal to unity. In other words, Ia (k) accounts for the interaction with both the atomic potential and electric field in the initial state, but completely neglects the atomic potential in the final state. This can be compared with the Keldysh approximation [1] which also neglects the atomic potential in the final state. For the present problem, the photoelectron spectrum in the Keldysh approximation is given by Eq. (34), where c(t) is to be substituted by cK (t) = e−iE0 t , see Appendix E. Thus the Keldysh approximation additionally neglects the interaction with the electric field in the initial state. This interaction reveals itself in the Stark shift and depletion of the initial state and in the present theory is accounted for by a difference between the unperturbed initial state and the adiabatic Siegert state, i.e., between cK (t) and ca (t) given by Eq. (64). We conclude that the Keldysh approximation emerges within the present theory in the limit of weak fields, ξ 1, and not too long pulses, when the Stark shift and depletion effects can be neglected. B. Rescattering part of the ionization amplitude

The contributions to Ir (k) from the different terms in Eq. (116), corresponding to the different CRTs from the set CRT(t), can be analyzed independently. We consider the simplest one-loop term and then discuss the structure of I (k) in the general case. 1. One-loop IRT

Let us substitute c1 (t) given by Eq. (109) for c(t) into Eq. (34) and calculate the integral. The corresponding action

(137)

(138)

Here tf (k, σ1 ) is the moment of rescattering for a one-loop IRT with the asymptotic velocity k and signature σ1 , see Sec. III B, so we call tf+ (k, σ1 ) the rescattering SP. The contribution to Ir (k) from this SP will be denoted by Ir(1) (k, σ1 ). + + Let ti± (k, σ1 ) = ti± (tf+ (k, σ1 )) and u+ f (k, σ1 ) = uf (tf (k, σ1 )). Omitting the arguments in ti+ (k, σ1 ), tf+ (k, σ1 ), and u+ f (k, σ1 ), we obtain Ir(1) (k, σ1 ) = R(u˜ + f )J (k, σ1 ),

(139)

where u˜ + f is defined by Eq. (108) and 1/2 + dti /dk + 3/2 J (k, σ1 ) = α U (ζ (tf )) si σ1 α(ti+ ) +

+

+

× eiS(k,tf )+iS(tf ,ti ) ca (ti+ ).

(140)

This result has only an intermediate value. Indeed, in the transmission case, ti+ (k, +) = ti+ (k) and S(k, tf+ ) + S(tf+ , ti+ ) = S(k, ti+ ), see Eq. (23). Hence there exist two different contributions to I (k) with the same action associated with the same one-loop transmission IRT: one comes from the ionization SP ti+ (k) and is included into Ia (k), let us denote it by Ia(1) (k), and the other comes from the rescattering SP tf+ (k, +) and is given by Ir(1) (k, +). To sum up their amplitudes, we note that ζ (tf+ ) = ζ (k), and hence (1) Ir(1) (k, +) = R(u˜ + f )Ia (k).

(141)

(1) Ia(1) (k) + Ir(1) (k, +) = T (u˜ + f )Ia (k).

(142)

Thus

Let I1 (k, σ1 ) denote the total contribution to I (k) associated with a one-loop IRT characterized by k and σ1 . Then the above results can be summarized as follows: I1 (k, σ1 ) = A(u˜ + f )J (k, σ1 ),

(143)

where A stands for T or R for σ1 = + and −, respectively. This asymptotics is uniform in terms of ξ and applies in the interval of k where uf O( 0 ) [this restriction is dictated by the condition of validity of Eq. (109)]. Each one-loop IRT from the set IRT(k) produces a similar contribution to I (k) given by the same formula (143). As can be seen from the derivation, the contribution Ir(1) (k, σ1 ) to I (k) from the rescattering SP tf+ (k, σ1 ) accounts for the interaction with the atomic potential during the rescattering. In the transmission case, this corrects the amplitude of the corresponding adiabatic term Ia(1) (k); in the

033415-14



reflection case, this gives a term with new action. In both cases, the total contribution (143) associated with a one-loop IRT is proportional to the exact scattering amplitude. This again can be compared with the Keldysh approximation. The approach initiated in [1] was later reformulated in the form of a time-dependent perturbation theory with respect to the interaction with the atomic potential in the final state [2,3]. This approach enables one to treat the rescattering effects, see, e.g., [54–60], but only in the Born approximation. To restore the exact scattering amplitude in this approach, one would have to sum up the whole Born series. The present results show that the summation is possible in the adiabatic approximation. 2. General structure of I(k)

The analysis can be continued. The contributions to Ir (k) from terms associated with higher CRTs in Eq. (116) are associated with the corresponding higher IRTs. Omitting further details of the derivation, the final result reads I (k) = IL (k, ), (144) IRT(k)

where the summation runs over all IRTs from the set IRT(k). The zero- and one-loop terms in Eq. (144) are given by Eqs. (134) and (143), respectively. For brevity and transparency of the result, we again give only an approximate form of IL (k, ) in the general case: IL (k, ) ≈ A(|uf |)A(|uL−1 |) · · · A(|u1 |) dti 1/2 iS(k,tf )+iS(tf , ) × 0 (F (ti )) e ca (ti ). (145) dk Here ti , tf , ul , and uf are the characteristics of the IRT as functions of its asymptotic velocity k and signature , see Sec. III B, and A is either T or R, depending on the result of the rescattering event. Equation (145) is obtained using Eq. (117), and hence is valid for ξ 1 in the interval of k where ul = O( −1 ) and uf = O( −1 ). It describes a contribution to the ionization amplitude from electrons that are ionized during the interval dti and run away to infinity along the corresponding IRT with the asymptotic velocity in the interval dk. The appearance of the derivative dti /dk in Eq. (145) ensures the conservation of flux. The individual contributions to I (k) given by Eq. (145) as functions of k satisfy IL (k, ) ∼ (r, −3), where r is the number of reflections in the IRT. Thus the leading-order term in the asymptotics of I (k) for → 0 is determined by transmission IRTs. We emphasize that Eq. (145) is only an approximate form: using the approach developed above, one can obtain the asymptotics of IL (k, ) which is uniform in terms of ξ and applies in a wider interval of k where ul O( 0 ) and uf O( 0 ), see Eqs. (134) and (143). The spectrum defined by IL (k, ) is localized in the interval max ¯ ¯ umin ∞ () k u∞ (),

(146)

where the corresponding IRT exists, see Eqs. (56). Only in this interval the moments ti , tl , and tf characterizing the IRT are real, and hence the SPs contributing to IL (k, ) are located close to the real axis, see Eqs. (131) and (138). Outside this interval, the SPs move deep into the complex time plane, so max ¯ ¯ IL (k, ) rapidly decays. We shall call umin ∞ () and u∞ ()

the classical boundaries of the spectrum. Localization of the spectrum within these boundaries is a signature of the adiabatic regime. C. The factorization formula

Similar to the situation discussed in Sec. IV D, the asymptotics of I (k) constructed in Secs. V A and V B is not uniform in terms of k. First, the contributions to Eq. (144) from IRTs whose initial moments ti (k, ) lie close to a zero of F (t) are not represented correctly. Second, the asymptotics of the individual terms in Eq. (144) does not apply in small intervals of k around the points where one of ul (k, ) or uf (k, ) turns zero. The lack of uniformity in these two cases results from the limitations of the asymptotics of cr (t) constructed in Sec. IV C. In addition, the asymptotics of the individual terms in Eq. (144) constructed in Secs. V A and V B does not apply near the boundaries of the interval (146), where two different IRTs with the same k and coalesce. Here we show how to resolve the latter deficiency. This will lead us to a factorization formula of the type proposed in [15]. We consider one particular pair of IRTs relevant to the situation discussed in [15], namely, a pair of one-loop reflection IRTs with = − created at the same zero of F (t). The left and right branches of the function u∞ (ti , −) shown by the dashed line in Fig. 2 illustrate the behavior of the asymptotic velocity for the long and short members of the pair, respectively. The initial moments ti (k, −) of the long and short IRTs will be denoted by til (k) and tis (k); the final moments tf (k, −) of the corresponding CRTs will be denoted by tf l (k) and tf s (k). Let km denote the extremum value of u∞ (ti , −) [in Fig. 2, this is the minimum umin ∞ (−) at ti = tim ]; km is the lower or upper classical boundary of the spectrum associated with the IRTs, see Eq. (146). The two IRTs coalesce at k = km , i.e., tf l (km ) = tf s (km ) ≡ tf m and til (km ) = tis (km ) = ti (tf m ) ≡ tim . The values of tf m and km are defined by ui (k, t) + uf (t) = 0 ∂ → t = tf m , k = km . (147) [ui (k, t) + uf (t)] = 0 ∂t Using Eq. (41), it can be shown that tf m can be found from F (ti (t))t˙i (t) = 2F (t) → t = tf m ,

(148)

km = −uf m − v(tf m ) + v∞ ,

(149)

uf m ≡ uf (tf m ) = v(tf m ) − v(tim ).

(150)

and then

where

In the case of general position, tf m − tim = O( −1 ), uf m = O( −1 ), and km = O( −1 ). The classical parameters of the coalescing point tim , tf m , and km are real. We need to introduce also the corresponding quantum objects. The rescattering SPs tf+ (k, −) defined by Eq. (137) for the long and short IRTs + , will be denoted by tf+l (k) and tf+s (k). They coalesce at k = km + + + + + + + i.e, tf l (km ) = tf s (km ) ≡ tf m . The values of tf m and km are

033415-15


I (k) from these IRTs is given by (see Appendix A)

t+

t +(k)

im

is

fl

il

−

im

t −(k) is

− iU− (ζls (k))R(u˜ + f s (k))Js (k),

tfm

tF

tim

Ils (k) = U+ (ζls (k))R(u˜ + f l (k))Jl (k)

t +(k)

t +(k)

t


t+

t

fm

+

t (k)

t −(k)

fs

il

FIG. 4. (Color online) The actual positions of the different points introduced in Sec. V C in complex time plane for the same model F (t) as in Figs. 1 and 2. The lines are traced by the corresponding points as k varies along the real axis in the vicinity of km . The solid and dashed lines correspond to the long and short IRTs, respectively.

defined by ⎫ ui (k, t) + u+ f (t) = 0 ⎬ ∂ ⎭ [ui (k, t) + u+ f (t)] = 0 ∂t

→ t = tf+m ,

+ k = km .

(151)

Using Eq. (106), it can be shown that the solutions to Eqs. (147) and (151) are related by tf+m = tf m + O( 1 ), + km = km +

E0 (tim ) + O( 2 ). (tf m − tim )F (tim )

(152a) (152b)

The moments of ionization ti± (k, −) for the long and short IRTs will be denoted by til± (k) and tis± (k). At the coalescing point, ± + + til± (km ) = tis± (km ) = ti± (tf+m ) ≡ tim . The quantum parameters ± + + are complex. The of the coalescing point tim , tf m , and km distribution of these points in the complex time plane is illustrated in Fig. 4. As k varies along the real axis in the vicinity of km , the two IRTs exist on the one side of km [for k > umin ∞ (−) in Fig. 2], but disappear on the other side + [for k < umin ∞ (−) in Fig. 2]. The rescattering SPs tf l (k) and + tf s (k) are located close to the real axis in the former case, and move deep into the complex plane in the latter case, in accordance with Eq. (138). The individual contributions to the ionization amplitude I (k) from the two IRTs discussed above are given by Eq. (143). Summing them up, as in Eq. (144), amounts to treating the two SPs tf+l (k) and tf+s (k) independently. This is justified if |ζls (k)| 1, where 2/3 3i ζls (k) = [S1 (k, tf+s (k)) − S1 (k, tf+l (k))] . (153) 4 + The value of ζls (k) turns zero at k = km ; this complex point will be called the quantum boundary of the spectrum. From + Eq. (152b) we have |km − km | = O( 1 ), hence ζls (k) may become small for real k in the vicinity of km . In this case, a uniform treatment of the two SPs is needed. Let Jl (k) and Js (k) denote the functions defined by Eq. (140) with σ1 = − for the long and short IRTs, respectively. Then the uniform + with respect to |k − km | asymptotics of the contribution to

(154)

+ + + + + where u+ f l (k) = uf (tf l (k)) and uf s (k) = uf (tf s (k)). Equation (154) resolves (for the particular pair of IRTs) the deficiency of Eq. (144) mentioned above. + As k departs sufficiently far from km along the real axis, the value of |ζls (k)| grows, and Ils (k) turns into the sum of the individual contributions from the two IRTs on the one side of km [inside the interval (146)], and vanishes on the other side of km [outside the interval (146)], in accordance with Eq. (144). In the intermediate region, where |ζls (k)| = O( 0 ), Eq. (144) fails. We now can clarify the structure of the photoelectron spectrum in this boundary region. Consider the interval + | = O( 1 ). |k − km

(155)

In this interval, we have ∂S1 (k, t) + 2 + = −uf m (k − km ) + O( ), ∂t t=tf m ∂ 2 S1 (k, t) + = −F (tf m )(k − km ) + O( 3 ), ∂t 2 t=tf+m ∂ 3 S1 (k, t) = ηuf m + O( 2 ), ∂t 3 t=tf+m

(156a) (156b) (156c)

where η = [2F˙ (t) − F˙ (ti (t))t˙i2 (t) − F (ti (t))tï (t)]t=tf m = O( 1 ). (157) Using this, we find O( 0 ), and

|tf+l (k)

−

tf+m |

= O( 0 ),

|tf+s (k)

−

2u2 1/3 fm + ζls (k) = −sgn[F (tim )] ) + O( 2 ). (k − km η

tf+m |

=

(158)

The first term in this expansion has order 0, thus Eq. (155) specifies the width of the boundary region. In this region, assuming that the contributions from all other IRTs can be neglected, the photoelectron spectrum can be presented in the form (see Appendix A) 2 Pls (k) = 2π |Ils (k)|2 = |R(|uf m |)|2 A Ai(ζls (k)) , (159) where

2 24/3 π u2/3 + + f m (a) + iS1 (km ,tf m ) C (t )e A = α r f m η1/3 ≈

28/3 π 2 |uf m |4/3 0 (F (tim ))|ca (tim )|2 . |η|2/3 (tf m − tim )|F (tim )|

(160a) (160b)

Equation (160b) is an approximation valid for ξ 1 and obtained using Eq. (110). Only the last factor in Eq. (159) containing the Airy function depends on k. This factor describes a typical interference structure of the photoelectron spectrum in the boundary region. In the adiabatic regime, this structure is characterized by the following properties. First, the width of the main peak of Pls (k) is O( 1 ) and its height is O( 1 ). We

033415-16



recall that in the case of general position, |k − km | = O( −1 ), the individual contributions from the one-loop reflections IRTs to I (k) have order 1, and hence Pls (k) = O( 2 ). Second, the + interference phase is reckoned from the complex quantum km instead of the real classical km boundary of the spectrum. + cannot be substituted by km in the argument We note that km of the Airy function in Eq. (159), because the coefficient of + + (k − km ) in Eq. (158) is O( −1 ) and |km − km | = O( 1 ), see + Eq. (152b). A difference between km and km leads to a shift of the interference structure; this seemingly small shift has the same order as the width of the main peak and therefore cannot be neglected. It can be shown that Eq. (159) remains valid in the interval |k − km | |km |, and the number of oscillations of Pls (k) in this interval is O( −3 ). In the monochromatic case, an interference structure of photoelectron spectra near the classical cutoff energy 10.007Up [52] associated with the same pair of long and short trajectories was first discussed in [59]. Equation (159) presents a factorization formula of the type proposed in [15]. In agreement with [15], the first factor in Eq. (159) is a scattering characteristic of the target atom taken at the velocity of the photoelectron at the moment of rescattering. Note that the value of this velocity uf m is determined by the laser pulse. The product of the second and third factors in Eq. (159) corresponds to the wave packet introduced in [15]. According to [15], the wave packet is almost independent of the properties of the target atom—this is the essence of the factorization. Equation (159) enables us to make this statement more precise. The shape of the wave packet defined by the third factor in Eq. (159) is indeed independent of the target, but there is a quantum shift defined by the second term on the rhs of Eq. (152b) which is determined by the target. The amplitude of the wave packet A defined by Eq. (160a) essentially depends on both the target and pulse. In the limit of weak fields, ξ 1, it depends on the target only via the factors 0 (F (tim )) and |ca (tim )|2 in Eq. (160b); the latter factor is close to unity if the depletion effects can be neglected. As was mentioned in the Introduction, two groups of authors have already succeeded in deriving the factorization formula for certain situations, see Eq. (31) in [37] and Eq.(59) in [41]. A brief discussion of the differences between the present theory and approaches used in [37,41] is in order here. First, the present theory treats pulses of arbitrary shape F (t), while Refs. [37,41] consider only monochromatic pulses, and for the Floquet approach used in [41] this restriction is essential. In fact, we believe that the question of whether the results obtained for monochromatic pulses can be applied to realistic pulses of finite length, especially in the case of few-cycle pulses and in situations where the depletion effects may be important, and, more generally, under which conditions such results are physically meaningful, is unjustifiably ignored in the literature. The second difference concerns the appearance of the exact scattering characteristic [the reflection coefficient |R(|uf m |)|2 , in the one-dimensional model] in the factorization formula. In the present theory, this naturally results from the consistent implementation of the adiabatic approximation. In the threedimensional ZRP model considered in [41], this results from the use of a relation [Eq. (20) in [41]] which is specific to

the model and whose physical meaning for more realistic models is obscure. In the strong-field approximation used in [37], the scattering characteristic can appear only in the Born approximation. Being aware of this fact, but, at the same time, being inspired by the success of the approach initiated in [15], see Refs. [16–27], the authors of Ref. [37] enforced the appearance of the exact scattering characteristic (the elastic scattering cross section, in the three-dimensional case) in their factorization formula. The third very important difference between the approaches concerns their applicability to strong laser fields. The present theory, including Eq. (159), is uniform in terms of ξ , i.e., the amplitude of the electric field. As was shown above, the tunneling regime of the Keldysh approximation emerges within the present theory in the limit of weak fields, ξ 1. The same is true for the results of Ref. [37]. This can be seen, e.g., from the fact that the so-called quantum correction to the classical cutoff energy 10.007Up [52] derived in [61], which appears in the formulation of [37] and defines a quantum shift of the position of the back rescattered ridge in the photoelectron momentum distribution, coincides with that following from the second term in Eq. (152b) if one substitutes there the unperturbed bound state energy E0 instead of the exact Siegert eigenvalue E0 (tim ). VI. ILLUSTRATIVE CALCULATIONS

In this section, we illustrate the theory by numerical calculations. We consider pulses with the Gaussian envelope, F (t) = F0 e−τ f (τ ),

(161)

τ = t,

(162)

2

where = 2/T ,

and f (τ ) defines the internal shape of the pulse. The “exact” results reported below are obtained by solving Eq. (30) with the initial condition (31) and calculating the spectrum (34) numerically. The adiabatic results are obtained by finding the saddle points and implementing the formulas derived in Secs. IV and V as prescribed, without further approximations. We first consider some model pulses, to demonstrate that the adiabatic results converge to the exact ones as the pulse length T grows uniformly with respect to the pulse amplitude F0 , and then discuss photoelectron spectra produced by more realistic pulses with the parameters of current interest. All the numerical results are for α = 1, hence E0 = −0.5. A. Half-cycle pulse: Stark shift and depletion effects

The adiabatic part ca (t) of the wave function is a foundation on which the rescattering part cr (t) and then the photoelectron spectrum P (k) are built. The accuracy of the present theory is determined by the accuracy of Eqs. (64), (78), and (81), so our first goal is to scrutinize these equations. To this end, we consider the simplest structureless pulse with f (τ ) = 1.

(163)

The electric field F (t) in this case does not have zeros, hence no rescattering occurs, the set CRT(t) is empty, and cr (t) = 0. Even though this model for F (t) is far from realistic laser

033415-17


c(t)/cK(t)

(a)

1.00

1.00

0.99

0.99

0.98 Re

0.98

Im

0.09

0.04

0.99

exact AA -1

T=10 0

1

T=5

0.8 0.6

Re

0.2 Im 0.1 0.0 -2

0.15 0.10

0.03

1.0

0.4

0.97

0.06

0.00 -2

c(t)/cK(t)

1.00

0.98

0.02

(b)


0.00 2 -2 1.0 0.8 0.6 0.4 0.2 0.0

0.05

T=20

-1

0

0.00 -2 1 2 1.0 T=10 0.8 0.6 0.4 0.2 0.0 0.2

T=40 -1

0

1

2

T=20

0.1

exact AA

0.1 0.0

-1

0

1

2

-2

-1

0

1

2

0.0 -2

-1

0

1

2

t/T FIG. 5. (Color online) The exact numerical solution to Eqs. (30) and (31) and adiabatic approximation (AA), Eq. (64), for half-cycle pulses with (a) F0 = 0.1 and (b) F0 = 1. The results are divided by the Keldysh approximation, cK (t) = e−iE0 t . In all the cases, the upper and lower panels show the real and imaginary parts of the ratio, respectively.

pulses, it is quite suitable for testing the performance of the adiabatic approximation. Let us consider two series of pulses with intermediate, F0 = 0.1, and large, F0 = 1, amplitudes and growing length T . The exact solution to Eqs. (30) and (31) is compared with ca (t) obtained from Eq. (64) in Fig. 5. Both functions are divided by cK (t) = e−iE0 t , so the departure of the exact results from unity illustrates an error incurred by the Keldysh approximation, see Appendix E. This error is caused by the Stark shift and depletion of the initial state, which are not accounted for by the Keldysh approximation. As can be seen from the figure, it grows with T and F0 , and for sufficiently long and/or intense pulses the Keldysh approximation fails. The onset of the adiabatic regime for the present case can be specified by |E0 |T = 2π , which yields T = 4π . Indeed, the adiabatic approximation ca (t) rapidly converges to the exact c(t) as T becomes larger than this value. This convergence is faster for larger F0 , in accordance with the condition (12). The photoelectron spectra for the same pulses are shown in Fig. 6. The results in the adiabatic approximation are obtained from Eq. (134). For the present pulse shape, the reference velocity v(t) monotonically decreases from 0, at t → −∞, to √ v∞ = − 2π F0 T [see Eq. (14b)], at t → ∞, so the classical max boundaries of the spectrum are umin ∞ = v∞ and u∞ = 0, see Eqs. (51) and (56). The equation u∞ (ti ) = k has only one max real solution ti (k) in the interval umin ∞ < k < u∞ , hence for each k in this interval, the set IRT(k) consists of only one zero-loop IRT, and there is only one pair of SPs defined by Eq. (129) which contributes to Ia (k). This establishes a one-to-one correspondence between the moment of ionization and the velocity of the ionized electron. For comparison, we also show the spectra in the Keldysh approximation obtained

by numerical integration from Eq. (34) with c(t) substituted by cK (t). Figure 6 confirms the conclusions drawn from Fig. 5. The adiabatic approximation rapidly converges to the exact results as T grows beyond the onset of the adiabatic regime. It works very well even in situations when almost complete ionization occurs, except for a narrow interval around k = 0, where the exact P (k) turns zero. The Keldysh approximation works well only when the Stark shift and depletion of the initial state can be neglected. This can be most clearly seen from the fact that, in all the cases, it nicely agrees with the exact results in the interval of k near and below the lower classical boundary umin ∞ . The electrons in this part of the spectrum are ionized in the initial part of the pulse, when the difference between c(t) and cK (t) is small, see Fig. 5. For larger t, and hence larger k, this difference grows, and the Keldysh approximation becomes progressively worse. While for weak fields the advantage of the adiabatic approximation reveals itself only for sufficiently long pulses [for F0 = 0.1 and T = 40, the adiabatic and Keldysh results in Fig. 6 differ from the exact P (k) at its maximum by 0.4% and 24%, respectively], for stronger fields the Keldysh approximation becomes qualitatively wrong. We note that the results shown in Fig. 6 are very far from the predictions of the first-order perturbation theory: the spectrum in this approximation is an even function of k, which is obviously not the case even for F0 = 0.1 and T = 10. B. One-cycle pulse: One-loop approximation

The rescattering part cr (t) of the wave function arises only after the electric field F (t) for the first time changes sign. Each next zero of F (t) generates a new infinite series of CRTs contributing to cr (t). Before discussing the general case, it

033415-18



-1

10

(a)

T=10

-2

T=20

T=40

P(k) (a.u.)

10

-3

10

-4

1x10

-5

1x10

exact AA KA

-6

10

-7

10

-8

10

1 -1 10 -2 10 -3 10 -4 1x10 -5 1x10 -6 10 -7 10 -8 10

-1

0

1

-2

-1

0

-4

-3

-2

-1

0

T=5

P(k) (a.u.)

(b)

-2

exact AA KA -6

-4

-2

T=20

T=10 0

2

-10

-8

-6

-4

-2

0

2

-15

-10

-5

0

k (a.u.) FIG. 6. (Color online) The exact results, adiabatic approximation (AA), and the Keldysh approximation (KA) for spectra produced by the same half-cycle pulses as in Fig. 5. (a) F0 = 0.1; the ionization probability Pion in the three cases is 1.34 × 10−2 , 1.21 × 10−2 , and 1.62 × 10−2 , respectively. (b) F0 = 1; the survival probability P0 in the three cases is 1.63 × 10−1 , 2.41 × 10−2 , and 5.60 × 10−4 , respectively. The vertical max dotted lines show the classical boundaries of the spectrum umin ∞ ≈ −0.886F0 T and u∞ = 0.

is instructive to analyze the relative role of the contributions from the different CRTs created at the same zero of F (t). To this end, we consider a pulse in which F (t) has only one zero, √ (164) f (τ ) = − 2eτ. The sign here √ complies with the previous case for t < 0, and the factor 2e is introduced to satisfy max |F (t)| = F0 . This model for F (t) is still far from realistic laser pulses, but is suitable for testing the theory. The function F (t) changes sign from positive to negative at t = 0, like in the model illustrated in Figs. 1 and 2. For t < 0, the set CRT(t) is empty and cr (t) = 0. For t > 0, the set CRT(t) consists of an infinite series of reflection CRTs created at t = 0, see Sec. III A. The corresponding contributions to cr (t) have orders L − 1/2, where L is the number of loops in the CRT, see Sec. IV C3. We restrict our treatment of the rescattering effects to the leading-order approximation. In this case, only the contribution from the one-loop CRT is to be retained in Eq. (116), see the dashed line in Fig. 1. Being continued after the rescattering event, this CRT generates one-loop transmission and reflection IRTs. For the present √ 2 pulse, v(t) = −( 2eF0 T /4)e−τ , hence v∞ = 0, umin ∞ = 0, max and u∞ ≈ 0.583F0 T . The equation u∞ (ti ) = k has two real max solutions ti (k) in the interval umin ∞ < k < u∞ ; the negative and positive solutions correspond to one-loop transmission and zero-loop IRTs, respectively. The equation u∞ (ti , −) = k has two real solutions ti (k, −) in the interval umin ∞ (−) < k < 0, corresponding to the long and short one-loop reflection IRTs discussed in Sec. V C, and one real solution in the interval 0 < k < umax ∞ (−), which is a continuation of the solution corresponding to the short IRT, where umin ∞ (−) ≈ −0.301F0 T max and umax ∞ (−) = u∞ , see the dashed line in Fig. 2. Let us

summarize the contributions to the ionization amplitude I (k) to be taken into account in the adiabatic results reported below: for k < 0, I (k) is represented by the uniform [with respect to |k − umin ∞ (−)|] asymptotics of the contributions from the long and short one-loop reflection IRTs given by Eq. (154); for k > 0, I (k) is represented by a sum of the individual contributions from the zero-loop, one-loop transmission, and short one-loop reflection IRTs given by Eqs. (134) and (143) with σ1 = + and −, respectively. The photoelectron spectra produced by six one-cycle pulses with the same combinations of F0 and T as in Fig. 6 are shown in Fig. 7. The onset of the adiabatic regime for the present case can be estimated by T = 8π . In contrast to the half-cycle case, now the ionized electrons can experience rescattering. As a result, there exist several different IRTs contributing to I (k) for each value of k. An interference of their contributions explains the oscillations of P (k). We take into account only two one-loop reflection IRTs for k < 0. As can be seen from Fig. 7, the uniform asymptotics (154) works very well even below the onset of the adiabatic regime, and its agreement with the exact results becomes perfect for larger T . The three IRTs taken into account for k > 0 also nicely reproduce the spectrum. One can notice that convergence of the adiabatic and exact results as T grows is somewhat slower in the region near and above the upper classical boundary umax ∞ . This is explained by the fact that the zero-loop and one-loop transmission IRTs coalesce at k = umax ∞ , so their contributions in this boundary region require a uniform treatment of the type developed in Sec. V C. This is doable, but falls beyond the limits of the present paper. Another region of slower convergence is near k = 0. Here contributions from higher IRTs, corresponding to higher reflection CRTs created at t = 0 (see Figs. 1 and 2), which are not taken into account in the present treatment, play a more important role. In

033415-19


-1

10 -2 10 -3 10 -4 1x10 -5 1x10 -6 10 -7 10 -8 10

T=10

T=20

T=40

P(k) (a.u.)

(a)


1 -1 10 -2 10 -3 10 -4 1x10 -5 1x10 -6 10 -7 10 -8 10

-1

0

1

2

-1

0

1

2

-1

0

1

2

3

T=20

T=10

T=5

P(k) (a.u.)

(b)

exact AA

exact AA -2

0

2

4

6 -4

-2

0

2

4

6

8

-5

0

5

10

k (a.u.) FIG. 7. (Color online) The exact results and adiabatic approximation (AA) for spectra produced by one-cycle pulses. (a) F0 = 0.1; the ionization probability Pion in the three cases is 3.64 × 10−2 , 2.58 × 10−2 , and 2.72 × 10−2 , respectively. (b) F0 = 1.0; the survival probability P0 in the three cases is 8.03 × 10−2 , 6.27 × 10−3 , and 3.73 × 10−5 , respectively. The vertical dotted lines show the classical boundaries of the min max spectrum umin ∞ (−) ≈ −0.301F0 T , u∞ = 0, and u∞ ≈ 0.583F0 T .

spite of these local difficulties, the adiabatic results certainly converge to the exact ones as T grows, and for larger F0 this convergence is faster, in accordance with Eq. (12). This confirms the eligibility of the leading-order approximation for cr (t), in which only the contributions from one-loop CRTs born at each zero of F (t) are retained in Eq. (116). We shall call it the one-loop approximation.

cr (t) is represented by a sum of the one-loop contributions (109) from each significant zero of F (t). The adiabatic results shown in Fig. 8 are obtained using such a mixed analytical and numerical approach. The adiabatic parameter for the present

(a)

-1

10

c

C. Few-cycle pulse

Having finished preliminary tests, we now demonstrate the performance of the adiabatic approximation for realistic fewcycle infrared laser pulses with

b

-3

10

a

-4

1x10

-5

1x10

exact AA+N

-6

10

-7

10

(165)

where noc = ωT /2π is the number of optical cycles in the pulse. The√full width at half-maximum (FWHM) for such pulses is ln 2T ≈ 0.833T . We consider pulses with T = 248.3 (FWHM = 5 fs) and noc = 2.25 (λ = 800 nm), which coincides with the parameters used in [15], for two values of the field amplitude F0 = 0.1 and 1, as in the previous cases. The exact and adiabatic results for spectra produced by these pulses are shown in Fig. 8. The adiabatic results are obtained in the one-loop approximation. It is sufficient to take into account only CRTs created at a finite number of significant zeros of F (t), which lie within a few widths of the pulse; tunneling ionization at far wings of the pulse is negligible. Even then, the total number of IRTs to be taken into account to obtain the ionization amplitude I (k) for long pulses becomes large. In this case, unless one is interested in individual contributions from some particular IRTs, it is easier to calculate the spectrum by exact numerical integration from Eqs. (34) and (60), where ca (t) is given by Eq. (64) and

(b)

-4

-3

-2

a

1 -1

10

P(k) (a.u.)

f (τ ) = cos(noc π τ ) = cos(ωt),

P(k) (a.u.)

-2

10

-2

10

-3

10

-36.0

a d

-35.6

-35.2

-1

0

1

0.006 0.003 0.000

2 0.002 0.000

3

4

e 34.4

f 34.8

35.2

e

f

-4

1x10

-5

1x10

-7

10

-40

b

exact AA+N

-6

10

-30

-20

-10

0

10

20

30

40

k (a.u.)

FIG. 8. (Color online) The exact results and adiabatic approximation (AA + N) for spectra produced by few-cycle pulses with λ = 800 nm and FWHM = 5 fs. (a) F0 = 0.1 (I = 3.51 × 1014 W/cm2 ), Pion = 2.97 × 10−2 . (b) F0 = 1.0 (I = 3.51 × 1016 W/cm2 ), P0 = 3.35 × 10−4 . The vertical dotted lines show the classical boundaries of the direct (without reflections) contributions to the spectrum umin ∞ = −17.40F0 and umax ∞ = 17.40F0 . The characters indicate the boundary interference structures produced by the pairs of IRTs labeled by the same characters in Fig. 9. The inserts in (b) show (in linear scale!) enlarged parts of the spectrum near structures a and f .

033415-20



40

u∞ (t,Σ)/F0 (a.u.)

30

f

e

b

F(t)/F0

20

1.0

0.5

10 0

0.0

-10

u∞ (t) -0.5 u∞ (t,−) u∞ (t,+ −) u∞ (t,+ + −) -1.0

-20 -30 -40 -1.0

c d

a -0.5

0.0

0.5

1.0

t/T ¯ as a FIG. 9. (Color online) The asymptotic velocity u∞ (t, ) ¯ for several function of the initial moment t and reduced signature IRTs with one reflection for few-cycle pulses defined by Eqs. (161) and (165) (compare with Fig. 2). The characters label the pairs of long and short (the left and right branches of the curve, respectively) IRTs producing the corresponding boundary interference structures in Figs. 8 and 10. The thick dashed line shows the electric field F (t).

-2

c

(a) 10

b

P(k) (a.u.)

-3

10

-4

1x10

a

-5

1x10

exact AA FF

-6

10

-7

10

-3.5

-2

-3.0

-2.5

3.5

4.0

10

(b) -3

10

P(k) (a.u.)

case can be estimated by = ω/|E0 | = 0.11, which is rather deep in the adiabatic regime. One can see that, except for some small intervals of k, where contributions from higher reflection CRTs play a more important role, the adiabatic results nicely agree with the exact ones. We again note that for larger F0 the adiabatic approximation works better; the level of agreement illustrated in the inserts in Fig. 8(b) is typical also for the other parts of this rapidly oscillating spectrum. Such agreement for as vastly different regimes in terms of the field amplitude F0 as shown in Fig. 8 demonstrates the overall performance of the present theory for infrared laser pulses. The classical boundaries of the spectrum for transmission ¯ is empty) in the present IRTs (the reduced signature max case are umin = −17.40F and u 0 ∞ ∞ = 17.40F0 . Beyond these boundaries, the spectra in Fig. 8 are determined by IRTs with reflections. Only IRTs with one reflection are taken into account in the one-loop approximation. The asymptotic ¯ for several IRTs of this type are shown velocities u∞ (t, ) in Fig. 9. The characters in this figure label several pairs of long and short IRTs created at the different zeros of F (t). The same characters in Fig. 8 indicate the boundary interference structures produced by the pairs, see Sec. V C. Figure 9 helps us interpret the spectra shown in Fig. 8. For example, the pairs ¯ = −) define the left of IRTs labeled in Fig. 9 by a and b ( min max u∞ (−) = −35.94F0 and right u∞ (−) = +39.07F0 classical boundaries of the whole spectrum, respectively. The value of F (t) at the moments of ionization for pair a is smaller than that for pair b, which explains the ratio of the amplitudes of the corresponding interference structures in the case when the depletion of the initial state can be neglected, see Fig. 8(a). On the other hand, the ionization in pair b happens later than that in pair a, which may inverse the ratio in the case when a considerable depletion between these two events occurs, see Fig. 8(b). In fact, the amplitude of structure b in Fig. 8(b) is so small that our “exact” numerical procedure fails to converge in this interval of k. For the same reason, (depletion) structure

a

f

-4

1x10

-5

1x10

exact AA FF

-6

10

-7

10

-36.0

-35.6

-35.2

34.4

34.8

35.2

k (a.u.) FIG. 10. (Color online) The exact results, adiabatic approximation (AA), and factorization formula (FF) for the same spectra produced by few-cycle pulses with (a) F0 = 0.1 and (b) F0 = 1.0 as in Fig. 8 in the regions of several most pronounced boundary interference structures. The vertical dotted (dash-dotted) lines show the classical km (quantum Re km+ ) boundaries of the spectrum associated with the corresponding pair of IRTs, see Sec. V C and Fig. 9.

c is not seen at all in Fig. 8(b). Another effect caused by increasing the field amplitude is the appearance in the spectrum of interference structures determined by higher IRTs. Indeed, while in Fig. 8(a) one can see only structures associated with ¯ = − (structures a, b, the simplest one-reflection IRTs with ¯ = +− (structures and c), the contributions from IRTs with ¯ = + + − (structure f ) can be clearly seen in d and e) and Fig. 8(b). We conclude this discussion with an illustration of the performance of the factorization formula whose derivation has provided an impetus for the whole work. Figure 10 shows several of the most pronounced boundary interference structures seen in Fig. 8 in larger scale. The adiabatic results in this figure were obtained from the uniform asymptotics (154) applied to each pair of interfering IRTs independently; in other words, the value of 2π |Ils (k)|2 for each pair is shown. Thus obtained results nicely reproduce the exact spectrum in the intervals of k where the particular pair gives a dominant contribution. For each of the five interference structures shown in Fig. 10, this is the case near the classical boundary of the spectrum. To obtain the spectrum in the intervals of k where two or more pairs of IRTs give comparable contributions, e.g., between structures a and c in Fig. 10(a), one has to sum up the amplitudes Ils (k) corresponding to each pair, which results in another type of interference structure. The agreement between the adiabatic and exact results in Fig. 10 is as good as in Fig. 8, which confirms that calculation of the integral in Eq. (34) by the steepest descent method does not introduce additional approximations. We also show the results obtained from the factorization formula (159) with A defined by Eq. (160a). This formula very accurately reproduces both the amplitude and phase of the interference structure in the vicinity of its

033415-21



classical boundary km , i.e., in the region where Eq. (159) applies. A virtually exact matching between the phase of the + oscillations in this region eventually proves the meaning of km as a true quantum boundary of the spectrum, see Sec. V C. As k departs from km , a difference between the exact interference phase and the adiabatic one given by Eq. (158) grows, and Eq. (159) works worse. The region of validity of Eq. (159) is wider for larger F0 . VII. CONCLUSIONS

The theory developed above is based on the adiabatic approximation. It operates with asymptotic expansions in terms of the adiabatic parameter and becomes exact in the limit → 0 uniformly with respect to the amplitude of the laser pulse; more precisely, its region of validity is defined by Eq. (12). The leading-order terms in the expansions have classical counterparts, so the adiabatic theory inherits, to some extent, transparency of the “simple man’s theory” [4–6]. The present theory has many tools in common (classical trajectories and actions, saddle points) with the tunneling regime of the strong-field approximation [1–3,7]. Therefore it is important to emphasize major differences. First, in contrast to the strong-field approximation, the adiabatic theory takes into account the interaction with the laser field in the initial state. This interaction is described by a difference between the unperturbed initial state and adiabatic Siegert state. It leads to the Stark shift and depletion of the initial state—the effects that become important for sufficiently long and/or intense pulses. Only for weak fields can these effects be approximately described in terms of the second-order Stark shift and semiclassical tunneling probability; for stronger fields, their quantitative description requires the knowledge of the exact Siegert state. Second, the adiabatic theory fully takes into account the interaction with the atomic potential in the final state. This, in particular, leads to the appearance of the exact scattering amplitudes in the asymptotics of the rescattering parts of the wave function and photoelectron spectrum, while in the strong-field approximation the scattering characteristics can appear only in the Born approximation. Third, the adiabatic theory is uniform in terms of the amplitude of the electric field, while the tunneling regime of the strong-field approximation is shown to emerge in this theory in the limit of weak fields. This difference in treating strong-field effects concerns both the adiabatic and rescattering parts of the wave function and spectrum. It reveals itself, e.g., in the expressions for the quantum shift of the classical boundary of the spectrum: the second term in Eq. (152b) reduces to the result derived in [61] if one substitutes there the unperturbed bound state energy E0 instead of the exact Siegert eigenvalue E0 (tim ), which is valid only for weak fields. Summarizing, in spite of the apparent similarity between the adiabatic theory and strong-field approximation, the two theories are based on quite different physical grounds, and the differences between them indicated above may lead to significant differences between their quantitative predictions. The approach used to develop the adiabatic theory of ionization in this work can be characterized as local. Generally speaking, its implementation requires constructing local asymptotic solutions in the different regions of the

complex time plane and then matching them together to obtain a global solution. This approach is most straightforward. More sophisticated approaches to the theory of nonadiabatic transitions to the continuum [62–65] may turn out to be simpler in implementation. It would be interesting to see whether the present theory can be developed in a more elegant way using one of these approaches. However, it should be understood that the adiabatic results do not depend on the method of their derivation because of the uniqueness of the asymptotic expansion [47,48]. This statement also gives an answer to another question often asked in a laser physics audience: Do the present results depend on the gauge? The answer is no. Whatever gauge is used in the derivation, it is only a matter of convenience; the adiabatic results obtained by consistently implementing the asymptotic approach (as is the case in the present study) do not depend on the gauge. Even though the atomic model considered in this work may seem to be oversimplified, the theory developed is general and can be extended to more realistic models. The next goal in the development of the theory is its extension to finite-range potentials in the three-dimensional case. ACKNOWLEDGMENTS

O.I.T. thanks the Russian Science Support Foundation for financial support. This work was supported in part by the PRESTO program of JST, Japan, and by a Grant-in-Aid for Scientific Research (C) from the MEXT, Japan. APPENDIX A: THE AIRY CORRECTING FACTORS

For consistency of the presentation, we recall here one result from the asymptotic analysis [48,66] multiply used in the main text. Consider an integral ∞ −1 I= A(z)ei S(z) dz, (A1) −∞

where the amplitude A(z) and action S(z) are sufficiently good functions so that for → 0 the integral can be evaluated using the steepest descent method. Suppose that there are only two SPs, z+ and z− , relevant to the evaluation of the integral; we have S (z± ) = 0. The integration path can be deformed into a steepest descent contour passing through either only one (we assume that it is z+ ) or both SPs. A simple asymptotics of I for → 0 is defined by the individual contributions from the SPs and given by + | arg ζ | < 2π/3, ISP , I≈ (A2) + − ISP + ISP , 2π/3 < | arg ζ | π, where ± ISP

and

= A(z± )e

i −1 S(z± )+iπ/4

2π S (z± )

1/2

2/3 3i . ζ = [S(z− ) − S(z+ )] 4

,

(A3)

(A4)

The validity of Eqs. (A2) requires |ζ | 1. The action S(z) may depend on some parameters which control the positions

033415-22



of the SPs, and hence the value of ζ . Equations (A2) do not hold near the points in the parameter space where the two SPs coalesce and ζ turns zero. A uniform with respect to ζ asymptotics of I , which remains valid in the coalescing case too, is given by + − I ≈ U+ (ζ )ISP − iU− (ζ )ISP ,

where U± (ζ ) =

√

(A5)

2 π exp ± ζ 3/2 [ζ 1/4 Ai(ζ ) ∓ ζ −1/4 Ai (ζ )]. 3 (A6)

Here Ai(ζ ) is the Airy function [67]. For |ζ | → ∞, we have U+ (ζ ) = 1 and U− (ζ ) = 0(i) in the sector | arg ζ | < 2π/3 (2π/3 < | arg ζ | π ), thus in this limit Eq. (A5) coincides with Eqs. (A2). In contrast to Eqs. (A2), Eq. (A5) holds for any value of ζ . In the coalescing case, |z+ − z− | 1, this equation can be simplified. In this case, z± ≈ z0 ± i(2/S3 )1/3 ζ 1/2 , S (z± ) ≈ ±i(2S32 )1/3 ζ 1/2 , and 2 2 − + (A7) ISP exp − ζ 3/2 ≈ iISP exp + ζ 3/2 , 3 3

On the other hand, Eq. (B1) can be solved using the Laplace contour integral method [53]. Thus we obtain ∞ 1 k 3 dk exp kz − , (B3a) W0 (z) = 4π 3/2 0 12 k 1/2 e±iπ/12 ∞ ik 3 dk exp ∓ikz ∓ . (B3b) W± (z) = 4π 3/2 0 12 k 1/2 Using these formulas, one can obtain an integral representation of a product of any two solutions to the Airy equation. APPENDIX C: SIEGERT STATES IN THE ELECTRIC FIELD

The Siegert states (SSs) are the solutions to the stationary Schrödinger equation having only outgoing waves in the asymptotic region [68]. Applying this concept to the present problem, the SSs for an electron interacting with a ZRP and a time-independent electric field F are defined by 1 d2 − − αδ(x) + F x − E φ(x) = 0, (C1a) 2 dx 2 Ai(e±2iπ/3 z) , Ai(e±2iπ/3 z0 ) Ai(z) φ(x 0) = φ(0) , Ai(z0 )

φ(x 0) = φ(0)

hence the leading-order term in Eq. (A5) is I≈

+ U (ζ )ISP ,

(A8a) 1/3 i 2 Ai(ζ ), ≈ A(z0 ) exp [S(z− ) + S(z+ )] 2π 2 S3 (A8b)

where z0 = (z− + z+ )/2, S3 = S (z0 ), and √ 2 3/2 1/4 ζ Ai(ζ ). ζ U (ζ ) = 2 π exp 3

Equations (A8a) and (A8b) hold in the vicinity of ζ = 0 defined by (/S3 )1/3 ζ 1/2 1. For |ζ | → ∞ and | arg ζ | < π, we have U (ζ ) = 1. Thus, in addition to the specified vicinity of ζ = 0, Eq. (A8a) coincides with the leading-order term in Eq. (A5), and hence gives a uniform with respect to |ζ | asymptotics of I , everywhere in the sector | arg ζ | < 2π/3, in which only one SP contributes to the integral. For evident reasons, we call functions U± (ζ ) and U (ζ ) the Airy correcting factors. The branches of multivalued functions in the above formulas are defined by the analytic continuation in the parameter space from a point where arg ζ = 0 and arg S (z± ) = ±π/2. APPENDIX B: INTEGRAL REPRESENTATION OF A PRODUCT OF TWO SOLUTIONS TO THE AIRY EQUATION

The product of any two solutions to the Airy equation satisfies [67] 3 d d − 4z − 2 W (z) = 0. (B1) dz3 dz As three linearly independent solutions to this equation, we choose W0 (z) = Ai(e−2iπ/3 z)Ai(e+2iπ/3 z), W± (z) = Ai(z)Ai(e

±2iπ/3

z).

(B2a) (B2b)

(C1c)

where z = 21/3 F −2/3 (F x − E),

(A9)

(C1b)

z0 = −21/3 F −2/3 E.

(C2)

We need the SSs for generally complex values of F , so a brief discussion of the analytic properties of the solutions to Eq. (C1a) as functions of F for fixed values of x and E is in order here. The solutions satisfying some decent initial conditions for the Poincaré theorem to apply (e.g., similar to that defining the Jost solutions [69]) are triple-valued functions of F having an essential singular point at F = 0. The physically motivated boundary conditions defining the SSs can be formulated on each of the six rays arg F = nπ , n = 0, . . . , 5, on the Riemann surface of the solutions as functions of F , where F takes pure real values, but their particular form depends on the ray. We have chosen rays arg F = 0, for positive F , and arg F = −π , for negative F , which correspond to the upper and lower sign in Eq. (C1b), respectively. For arg F = 0 (arg F = −π ), the solutions to Eqs. (C1) have only an outgoing wave, as x → −∞ (x → +∞), and exponentially decay, as x → +∞ (x → −∞). Importantly, the boundary conditions (C1b) and (C1c) in the two cases cannot be obtained from each other by the analytic continuation in F , i.e., they define two different solutions. This explains why the two cases require separate treatment. Summarizing, Eqs. (C1) define the SSs for positive (arg F = 0, the upper sign) and negative (arg F = −π , the lower sign) real values of F ; for complex F , the SSs can be obtained by the analytic continuation from a suitable ray. Equations (C1) constitute an eigenvalue problem: the solutions exist only for a discrete set of generally complex values of E. Applying the matching condition at x = 0, this differential eigenvalue problem can be cast into the form of a

033415-23



transcendental equation, (F) (a.u.)

-0.50 ∓iπ/6

1/3

exact PT

-0.55

0

(2F ) , (C3) W± (z0 ) = 4π α where functions W± (z) are defined by Eq. (B2b). This defines the eigenvalue E. The eigenfunction is given by Eqs. (C1b) and (C1c), where the surface amplitude φ(0) is to be determined by the normalization condition. We require ∞ φ 2 (x)dx = 1. (C4)

-0.60 (a)

1 (b)

Γ0(F) (a.u.)

e

−∞

Note that there is no complex conjugation in this formula, which is a general property of the theory of SSs [68,70–72]. Apart from this difference (which is not essential for bound states in the one-dimensional case anyway) Eq. (C4) coincides with the usual normalization condition for a bound state. However, in the presence of the electric field the SS eigenfunctions do not decay (in fact, they grow) in the asymptotic region |x| → ∞, hence the integral in Eq. (C4) does not converge in the usual sense of the word and requires a regularization. This can be done using the relation ∞ Ai2 (z)dz = −Ai2 (z0 )f (z0 ), (C5)

-1

10

-2

10

exact SC

-3

φ0(0;F)−1 (a.u.)

10

0.2 0.1

Reφ0(0;F)−1

0.0 -0.1

Imφ0(0;F)

-0.2 -0.3 (c) 0.0

0.5

where Ai (z) f (z) = , Ai(z)

f (z) = z − f (z). 2

(C6)

Thus we obtain from Eqs. (C1b), (C1c), and (C4) N φ 2 (0) = 1,

(C7)

where the normalization factor N is equivalently given by N = (2F )−1/3 f 2 (z0 ) + e±iπ/3 f 2 (e±2iπ/3 z0 ) (C8a) = −e±iπ/6

4π α 2 dW± (z0 ) . F dz0

(C8b)

This defines φ(0), and hence φ(x), up to a sign. Let E0 (F ) and φ0 (x; F ) denote the SS which for F = 0 coincides with the only bound state supported by the ZRP— only this SS is needed for the discussion in the main text. In the presence of the electric field, this state becomes a resonance. We present E0 (F ) in the form i (C9) E0 (F ) = E0 (F ) − 0 (F ). 2 The energy E0 (F ) and width 0 (F ) obtained by solving Eq. (C3) and the surface amplitude φ0 (0; F ) defined by Eqs. (C7) and (C8) as functions of F are shown in Fig. 11. For real F and |F | → 0, these functions can be expanded in asymptotic series in terms of the parameter ξ = |F |/α 3 [see the second of Eqs. (10)]. We thus obtain 5F 2 α2 − , 2 8α 4 2α 3 , 0 (F ) ≈ α 2 exp − 3|F | E0 (F ) ≈ −

and

5F 2 φ0 (0; F ) ≈ α 1/2 1 − . 4α 6

(C10) (C11)

(C12)

1.0

1.5

2.0

F (a.u.)

z0

FIG. 11. (Color online) The energy (a), width (b), and surface amplitude (c) for the Siegert state which coincides with the only bound state supported by the ZRP for F = 0. Solid lines in (a) and (b): exact results obtained by solving Eq. (C3). Dashed lines in (a) and (b): perturbation theory (PT) and semiclassical (SC) approximations defined by Eqs. (C10) and (C11), respectively. The results are for α = 1.

The leading-order terms in Eqs. (C10) and (C12) agree with Eqs. (4). The second term in Eq. (C10) coincides with the result of the second-order perturbation theory (the quadratic Stark shift), and Eq. (C11) coincides with the result of the semiclassical approximation [53,73]. Equations (C10) and (C11) are compared with the exact results in Fig. 11; the error remains less than 10% up to F = 0.45 and F = 0.05 for E0 (F ) and 0 (F ), respectively. We note an exact relation E0

e±iπ/2 25 π 3 α 3 34 6 (2/3)

= 0,

(C13)

where the two signs correspond to the two signs in Eq. (C3). Let us mention without going into detail that in addition to the solution discussed above, Eq. (C3) has two infinite series of solutions which for real F and |F | → 0 are located near the zeros of the two Airy functions whose product is denoted by W± (z0 ). All these solutions lie in the lower half of the complex E plane, forming two series which run away to infinity along the rays arg E = 0 and arg E = −2π/3, see Fig. 12. These SSs do not have counterparts for F = 0 and correspond to resonances induced by the electric field.

033415-24



The denominator in these formulas turns zero if E and F are related by Eq. (C3), thus the SSs reveal themselves as poles of the coefficients (D4). In particular, using Eqs. (C7) and (C8b), we obtain

Im E (a.u.)

0

-1

−αφ02 (0; F ) . E − E0 (F )

T (k; F )|E→E0 (F ) =

-2

F=0.2 -2

-1

0

1

2

3

(D5)

It is clear that for real E and F, there is no flux associated with the state ϕ(x, k; F ), hence the appropriate scattering characteristic is the elastic phase shift. For example, for arg F = 0 we have

4

Re E (a.u.) FIG. 12. (Color online) The energies of the Siegert states defined by Eq. (C3) for F = 0.2 and α = 1. The arrow indicates the state shown in Fig. 11.

ϕ(x, k; F )|x→−∞ ∝

e−is(x) − ieis(x)+2iδ(k;F ) , |x|1/4

(D6)

where s(x) =

APPENDIX D: SCATTERING STATES IN THE ELECTRIC FIELD

Here we discuss the solutions to Eq. (C1a) corresponding to scattering states. Let us introduce functions [we again use notation (C2); the upper and lower signs stand for the analytic continuation from the rays arg F = 0 and arg F = −π , respectively] √ 3 iπ 2 πk ik − Ai(z), (D1a) e+ (x) = exp (2F )1/6 3F 4 √

iπ iπ 2 πk ik 3 e− (x) = − ± Ai e±2iπ/3 z , exp − (2F )1/6 3F 4 6 (D1b) √ where k = 2E. We have de+ (x) de− (x) − e+ (x) = 2ik. (D2) dx dx These functions satisfy Eq. (C1a) for x = 0 and have a suitable behavior in the asymptotic region: e+ (x) exponentially decays (arg F = 0) or has only outgoing wave (arg F = −π ) as x → +∞, and e− (x) does the same, but in reverse order, as x → −∞. By analogy with the usual scattering theory [69], we define the scattering states ϕ(x, k; F ) in the presence of the electric field as the solutions to Eq. (C1a) satisfying e− (x)

F 1/2 |2x|3/2 E|2x|1/2 , + 3 F 1/2

(D7)

and exp [2iδ(k; F )] =

T (k; F ) . T ∗ (k; F )

(D8)

The first and second terms in Eq. (D6) represent the incoming and reflected waves, respectively. One can “switch off” the ZRP by setting into Eqs. (D4) α = 0; then R(k; F ) = 0, T (k; F ) = 1, and δ(k; F ) = 0. Thus a nonzero phase shift results from the interaction with the atomic potential. The energy dependence of δ(k; F ) is illustrated in Fig. 13. The phase shift acquires a positive increment less than π each time when E passes above a SS pole, see Fig. 12, and decays between the poles. In spite of the apparent similarity between Eqs. (7) and (D3), these two sets of scattering states have quite different physical natures. In particular, coefficients (D4) in contrast to coefficients (8) do not have the meaning of transmission and reflection amplitudes. However, it is clear that states (D3) must converge in some sense to states (7) as the electric field vanishes. For real F , this happens in the limit defined by |F | → 0,

ϕ(x, k; F ) e+ (x) + R(k; F )e− (x), x 0, = , arg F = 0, x 0, T (k; F )e+ (x),

Im k 3 > 0.

(D9)

3.0

F=0.2 2.5

(D3a)

2.0

δ(k;F)

ϕ(x, k; F ) T (k; F )e− (x), x 0, = , arg F = −π. e− (x) + R(k; F )e+ (x), x 0,

1.5 1.0

(D3b) 0.5

Applying the matching condition at x = 0, we find k , k − iαe+ (0)e− (0) 2 iαe± (0) . R(k; F ) = k − iαe+ (0)e− (0)

T (k; F ) =

0.0 -2

(D4a) (D4b)

-1

0

1

2

3

4

E (a.u.) FIG. 13. The scattering phase shift defined by Eq. (D8) for F = 0.2 and α = 1.

033415-25



Importantly, the latter condition is fulfilled at the upper rim of the cut along the positive real semiaxis on the physical energy sheet, i.e., where the usual scattering theory is formulated [69]. In this limit we have e± (x) → e±ikx ,

(D10)

which explains our choice of the normalization factor in Eqs. (D1). We also obtain 5iαF 2 T (k; F ) = T (k) 1 − T (k) + · · · , (D11a) 8k 7 5iF + ··· . (D11b) R(k; F ) = R(k) 1 ± 12k 3

APPENDIX E: THE KELDYSH APPROXIMATION

In the Keldysh approximation [1], the spectrum of the ionized electrons in the present model is given by 2 ∞ ∞ PK (k) = dt dxχk∗ (x, t)F (t)xφ0 (x)e−iE0 t . (E1) −∞

−∞

Here χk (x, t) is the Volkov state defined by ∂χk (x, t) 1 ∂2 + F (t)x χk (x, t), i = − ∂t 2 ∂x 2 and

i t 2 ui (k, t )dt . χk (x, t) = exp iui (k, t)x − 2 0

(E2)

(E3)

Integrating by parts, we find that PK (k) is given by Eq. (34), where c(t) is to be substituted by

This shows that in the specified sense coefficients (D4) can be interpreted as the transmission and reflection amplitudes in the presence of the electric field.

i.e., by the unperturbed initial state, see Eq. (31).

[1] L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1945 (1964) [Sov. Phys. JETP 20, 1307 (1965)]. [2] F. H. M. Faisal, J. Phys. B 6, L89 (1973). [3] H. R. Reiss, Phys. Rev. A 22, 1786 (1980). [4] H. B. van Linden van den Heuvell and H. G. Muller, in Multiphoton Processes, edited by S. J. Smith and P. L. Knight (Cambridge University, Cambridge, England, 1988). [5] T. F. Gallagher, Phys. Rev. Lett. 61, 2304 (1988). [6] P. B. Corkum, N. H. Burnett, and F. Brunel, Phys. Rev. Lett. 62, 1259 (1989). [7] D. B. Miloˇsević, G. G. Paulus, D. Baur, and W. Becker, J. Phys. B 39, R203 (2006). [8] H. Niikura, F. Légaré, R. Hasbani, A. D. Bandrauk, M. Yu. Ivanov, D. M. Villeneuve, and P. B. Corkum, Nature (London) 417, 917 (2002). [9] H. Niikura, F. Légaré, R. Hasbani, M. Yu. Ivanov, D. M. Villeneuve, and P. B. Corkum, Nature (London) 421, 826 (2003). [10] J. Itatani, J. Levesque, D. Zeidler, H. Niikura, H. Pépin, J. C. Kieffer, P. B. Corkum, and D. M. Villeneuve, Nature (London) 432, 867 (2004). [11] T. Kanai, S. Minemoto, and H. Sakai, Nature (London) 435, 470 (2005). [12] M. Lein, Phys. Rev. Lett. 94, 053004 (2005). [13] S. Baker, J. S. Robinson, C. A. Haworth, H. Teng, R. A. Smith, C. C. Chirilˇa, M. Lein, J. W. G. Tisch, and J. P. Marangos, Science 312, 424 (2006). [14] M. Lein, J. Phys. B 40, R135 (2007). [15] T. Morishita, A.-T. Le, Z. Chen, and C. D. Lin, Phys. Rev. Lett. 100, 013903 (2008). [16] M. Okunishi, T. Morishita, G. Prümper, K. Shimada, C. D. Lin, S. Watanabe, and K. Ueda, Phys. Rev. Lett. 100, 143001 (2008). [17] D. Ray et al., Phys. Rev. Lett. 100, 143002 (2008). [18] T. Morishita, A.-T. Le, Z. Chen, and C. D. Lin, New J. Phys. 10, 025011 (2008). [19] X. X. Zhou, Z. Chen, T. Morishita, A.-T. Le, and C. D. Lin, Phys. Rev. A 77, 053410 (2008).

[20] A.-T. Le, T. Morishita, and C. D. Lin, Phys. Rev. A 78, 023814 (2008). [21] S. Minemoto, T. Umegaki, Y. Oguchi, T. Morishita, A.-T. Le, S. Watanabe, and H. Sakai, Phys. Rev. A 78, 061402(R) (2008). [22] S. Micheau, Z. Chen, T. Morishita, A.-T. Le, and C. D. Lin, J. Phys. B 42, 065402 (2009). [23] Z. Chen, A.-T. Le, T. Morishita, and C. D. Lin, J. Phys. B 42, 061001 (2009). [24] T. Morishita, M. Okunishi, K. Shimada, G. Prümper, K. Shimada, Z. Chen, S. Watanabe, K. Ueda, and C. D. Lin, J. Phys. B 42, 105205 (2009). [25] Z. Chen, A.-T. Le, T. Morishita, and C. D. Lin, Phys. Rev. A 79, 033409 (2009). [26] A.-T. Le, R. R. Lucchese, S. Tonzani, T. Morishita, and C. D. Lin, Phys. Rev. A 80, 013401 (2009). [27] T. Morishita, T. Umegaki, S. Watanabe, and C. D. Lin, J. Phys.: Conf. Ser. 194, 012011 (2009). [28] J. Xu, H.-L. Zhou, Z. Chen, and C. D. Lin, Phys. Rev. A 79, 052508 (2009). [29] A.-T. Le, R. R. Lucchese, and C. D. Lin, J. Phys. B 42, 211001 (2009). [30] Z. Chen, T. Wittmann, B. Horvath, and C. D. Lin, Phys. Rev. A 80, 061402(R) (2009). [31] M. Meckel et al., Science 320, 1478 (2008). [32] H. J. Wörner, H. Niikura, J. B. Bertrand, P. B. Corkum, and D. M. Villeneuve, Phys. Rev. Lett. 102, 103901 (2009). [33] O. Smirnova, Y. Mairesse, S. Patchkovskii, N. Dudovich, D. Villeneuve, P. Corkum, and M. Yu. Ivanov, Nature (London) 460, 972 (2009). [34] M. Busuladˇzić, A. Gazibegović-Busuladˇzić, D. B. Miloˇsević, and W. Becker, Phys. Rev. Lett. 100, 203003 (2008). [35] M. Okunishi, R. Itaya, K. Shimada, G. Prümper, K. Ueda, M. Busuladˇzić, A. Gazibegović-Busuladˇzić, D. B. Miloˇsević, and W. Becker, Phys. Rev. Lett. 103, 043001 (2009).

cK (t) = e−iE0 t ,

033415-26

(E4)



[36] S. Odˇzak and D. B. Miloˇsević, Phys. Rev. A 79, 023414 (2009). ˇ [37] A. Cerki´ c, E. Hasović, D. B. Miloˇsević, and W. Becker, Phys. Rev. A 79, 033413 (2009). [38] D. B. Miloˇsević, W. Becker, M. Okunishi, G. Prümper, K. Shimada, and K. Ueda, J. Phys. B 43, 015401 (2010). [39] M. V. Frolov, N. L. Manakov, T. S. Sarantseva, M. Yu. Emelin, M. Yu. Ryabikin, and A. F. Starace, Phys. Rev. Lett. 102, 243901 (2009). [40] M. V. Frolov, N. L. Manakov, T. S. Sarantseva, and A. F. Starace, J. Phys. B 42, 035601 (2009). [41] M. V. Frolov, N. L. Manakov, and A. F. Starace, Phys. Rev. A 79, 033406 (2009). [42] K. Sasaki, X. M. Tong, and N. Toshima, J. Phys. B 42, 165603 (2009). [43] C. Cornaggia, J. Phys. B 42, 161002 (2009). [44] Y. Guo, P. Fu, Z.-C. Yan, J. Gong, and B. Wang, Phys. Rev. A 80, 063408 (2009). [45] O. I. Tolstikhin, T. Morishita, and S. Watanabe, J. Phys.: Conf. Ser. 194, 032032 (2009). [46] O. I. Tolstikhin, Phys. Rev. A 74, 042719 (2006). [47] F. W. J. Olver, Asymptotics and Special Functions (Academic, New York, 1974). [48] M. V. Fedoryuk, Metod Perevala (Nauka, Moscow, 1977). [49] R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals (McGraw-Hill, New York, 1965). [50] M. Gavrila, J. Phys. B 35, R147 (2002). [51] P. B. Corkum, Phys. Rev. Lett. 71, 1994 (1993). [52] G. G. Paulus, W. Becker, W. Nicklich, and H. Walther, J. Phys. B 27, L703 (1994). [53] L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Nonrelativistic Theory) (Pergamon, Oxford, 1977). [54] X. Mu, J. Ruscheinski, and B. Crasemann, Phys. Rev. A 42, 2949 (1990). [55] W. Becker, A. Lohr, and M. Kleber, J. Phys. B 27, L325 (1994); 28, 1931 (1995).

[56] M. Lewenstein, K. C. Kulander, K. J. Schafer, and P. H. Bucksbaum, Phys. Rev. A 51, 1495 (1995). [57] A. Lohr, M. Kleber, R. Kopold, and W. Becker, Phys. Rev. A 55, R4003 (1997). [58] D. B. Miloˇsević and F. Ehlotzky, Phys. Rev. A 57, 5002 (1998); 58, 3124 (1998). [59] S. P. Goreslavskii and S. V. Popruzhenko, Phys. Lett. A 249, 477 (1998); Pis’ma Zh. Eksp. Teor. Fiz. 68, 858 (1998) [Sov. Phys. JETP Lett. 68, 902 (1998)]; Zh. Eksp. Teor. Fiz. 117, 895 (2000) [Sov. Phys. JETP 90, 778 (2000)]. [60] Z. Chen, T. Morishita, A.-T. Le, and C. D. Lin, Phys. Rev. A 76, 043402 (2007). [61] M. Busuladˇzić, A. Gazibegović-Busuladˇzić, and D. B. Miloˇsević, Laser Phys. 16, 289 (2006). [62] E. A. Solov’ev, Zh. Eksp. Teor. Fiz. 70, 872 (1976) [Sov. Phys. JETP 43, 453 (1976)]. [63] V. I. Savichev, Zh. Eksp. Teor. Fiz. 100, 1449 (1991) [Sov. Phys. JETP 73, 803 (1991)]. [64] O. I. Tolstikhin, Phys. Rev. A 77, 032711 (2008). [65] D. I. Bondar, W.-K. Liu, and G. L. Yudin, Phys. Rev. A 79, 065401 (2009). [66] C. Chester, B. Friedman, and F. Ursell, Proc. Cambridge Philos. Soc. 53, 599 (1957). [67] Handbook of Mathematical Functions, edited by M. Abramowitz and I. A. Stegun (Dover, New York, 1972). [68] A. J. F. Siegert, Phys. Rev. 56, 750 (1939). [69] R. G. Newton, Scattering Theory of Waves and Particles (Springer-Verlag, New York, 1982). [70] O. I. Tolstikhin, V. N. Ostrovsky, and H. Nakamura, Phys. Rev. A 58, 2077 (1998). [71] G. V. Sitnikov and O. I. Tolstikhin, Phys. Rev. A 67, 032714 (2003). [72] P. A. Batishchev and O. I. Tolstikhin, Phys. Rev. A 75, 062704 (2007). [73] S. Geltman, J. Phys. B 11, 3323 (1978).

033415-27