Weak-field asymptotic theory of tunneling ionization including the first ...

PHYSICAL REVIEW A 91, 063410 (2015)

Weak-field asymptotic theory of tunneling ionization including the first-order correction terms: Application to molecules Vinh H. Trinh,1 Vinh N. T. Pham,1 Oleg I. Tolstikhin,2 and Toru Morishita1 1

Department of Engineering Science, The University of Electro-Communications, 1-5-1 Chofu-ga-oka, Chofu-shi, Tokyo 182-8585, Japan 2 Moscow Institute of Physics and Technology, Dolgoprudny 141700, Russia (Received 23 March 2015; published 15 June 2015) The weak-field asymptotic theory (WFAT) of tunneling ionization including the first-order correction terms is validated for molecules by comparison with accurate calculations of molecular Siegert states in a static electric field. Both fundamental observables related to tunneling ionization, namely, the ionization rate and transverse momentum distribution of the ionized electrons, are considered. This complements our previous study of atoms [V. H. Trinh et al., Phys. Rev. A 87, 043426 (2013)]. Similarly to the atomic case, the first-order terms essentially improve the agreement between the WFAT and accurate results in a wide interval of fields up to the onset of over-the-barrier ionization. This establishes the WFAT including the first-order correction terms as an appealing alternative to accurate calculations in the tunneling regime. In addition to demonstrating the quantitative performance of the WFAT, the theory is shown to be helpful for understanding the field and orientation dependencies of the observables. In particular, we show that the first-order terms account for a deviation of the shape of the orientation dependence of the ionization rate of a molecule from that of the ionizing orbital as the field grows, which has important implications for strong-field molecular imaging techniques. This prediction of the WFAT is confirmed by comparison with time-dependent calculations. DOI: 10.1103/PhysRevA.91.063410

PACS number(s): 32.60.+i, 33.80.Rv, 42.50.Hz

I. INTRODUCTION

In the presence of an external static electric field, bound states of atoms and molecules turn into Siegert states (SSs) represented by the eigensolutions of the stationary Schrödinger equation satisfying outgoing-wave boundary conditions. The real and imaginary parts of the complex SS eigenvalue define the Stark-shifted energy and ionization rate of the state, respectively. The SS eigenfunction defines the transverse (with respect to the direction of the field) momentum distribution (TMD) of the ionized electrons. These are the main properties of atomic and molecular systems characterizing their interaction with a static electric field. The current interest to the theoretical methods of evaluating these properties is motivated by their applications in strong-field physics [1]. For example, the SS eigenvalue and TMD amplitude are required for calculating photoelectron momentum distributions produced by intense low-frequency laser pulses within the adiabatic theory [2]. The knowledge of the ionization rate is a prerequisite for the analysis of experimental photoelectron and harmonic spectra generated by such pulses [3–10]. The TMD determines the lateral width of photoelectron momentum distributions in circularly polarized laser fields [11], etc. We have developed a powerful method to calculate oneelectron SSs in the combination of an atomic [12] or molecular [13,14] potential and an electric field. This method is capable of treating both the tunneling and the over-the-barrier ionization regimes. For neutral systems in the ground state, a transition between the two regimes occurs at a critical field Fc ∼ 0.1 a.u., which corresponds to intensities I ∼ 3.5 × 1014 W/cm2 . At higher intensities, accurate but rather laborious calculations of the type reported in Refs. [12–14] are indispensable for obtaining theoretical predictions applicable at the quantitative level. Meanwhile, at lower intensities, an approximate analytical solution of the SS eigenvalue problem is possible. In the tunneling regime, F < Fc , the energy of a 1050-2947/2015/91(6)/063410(18)

SS has power-series dependence on F and can be obtained using perturbation theory [15]. The characteristics related to tunneling ionization are exponentially small in F and can be found by asymptotic methods in the form of asymptotic expansions in F . The leading-order term in the expansion of the ionization rate from spherically symmetric atomic potentials was obtained in Refs. [15–18]. The weak-field asymptotic theory (WFAT) introduced in Ref. [19] generalizes these early results to molecular potentials without any symmetry. The WFAT yields asymptotic expansions not only for the ionization rate, but also for the TMD amplitude and other properties of SSs needed for applications. Recently, this theory was generalized to many-electron atoms and molecules treated in the frozen-nuclei approximation [20]. The effects of nuclear motion also can be accounted for within the WFAT [21,22]. The WFAT is much simpler in implementation than exact calculations, but its accuracy is limited. For example, the wellknown leading-order asymptotic formula for the ionization rate of hydrogen in the ground state [15] overestimates the rate at F = 0.05 a.u. (which is well below Fc ≈ 0.125 a.u. for this system) by 68%. This is too large an error for a quantitative theory. The result can be improved by incorporating higherorder terms in the asymptotic expansion. The first-order correction for the ionization rate of an arbitrary state of hydrogen was derived in Ref. [23]. The asymptotic formula including this term underestimates the rate for the ground state at the same field by 7%; thus, the error is reduced by a factor of 10 and becomes tolerable for most of applications. It is very tedious to derive the second-order, let alone higher-order, terms analytically. At the same time, the example of hydrogen, for which such terms were found numerically [24], shows that their inclusion considerably improves the result only near the upper boundary Fc of the tunneling regime; at weaker fields, the first-order correction accounts for the major part of the difference between the leading-order approximation and exact results. This suggests that truncating the asymptotic

063410-1

©2015 American Physical Society

TRINH, PHAM, TOLSTIKHIN, AND MORISHITA


expansions of the WFAT beyond the first-order terms is a sensible compromise between accuracy and computational labor. Within the one-electron WFAT [19], the first-order correction terms for SSs in arbitrary potentials were derived in Ref. [25]. The result of Ref. [23] for the ionization rate of hydrogen emerges from this general theory as a particular case. A necessary step in establishing the WFAT including the first-order terms as an appealing alternative to accurate numerical calculations in the tunneling regime is to compare the results of the two approaches for a number of atomic and molecular systems. In Ref. [25], the theory was validated by calculations for several spherically symmetric potentials describing noble-gas atoms in the single-active-electron approximation. It was shown that the first-order terms essentially improve the agreement between the WFAT and accurate results for the ionization rate and TMD over a wide interval of fields up to the onset of over-the-barrier ionization. In this paper we present similar calculations for one-electron molecular ions H2 + and HeH2+ treated in the frozen-nuclei approximation. In addition to validating the WFAT for molecules, we also discuss some qualitative effects accounted for by the first-order terms. In particular, the molecular case differs from the atomic case considered in Ref. [25] by an additional dependence of the observables on the orientation of the molecule with respect to the field. Thus, while the tunneling ionization rate of an atom in a given field is a number, and including the first-order terms in the asymptotic expansion just makes this number closer to its exact value, for a molecule the rate depends on its orientation and the first-order terms can modify the shape of this dependence as the field grows. The possibility of such a modification of the shape of the orientation dependence of molecular ionization rates has important implications for strong-field molecular imaging techniques. We show that the WFAT enables one to quantify the effect and analyze the underlying physical mechanisms. The paper is organized as follows. In Sec. II we summarize basic equations of the theory of SSs in an electric field [12,19]. In Sec. III we summarize formulas needed to implement the WFAT [19,25]. Some details of the present numerical procedure are discussed in the Appendix. In Sec. IV we present and discuss the results of our calculations. The WFAT results are compared with accurate calculations of SSs in a static field, but we also give examples of comparison with time-dependent calculations of tunneling ionization in low-frequency laser fields. Section V concludes the paper. II. SIEGERT STATES IN AN ELECTRIC FIELD

Within the single-active-electron and frozen-nuclei approximations adopted in this work, the interaction of a molecule with an external static uniform electric field F = F ez , F 0, is described by the stationary Schrödinger equation for the active electron (atomic units are used throughout) 1 − 2 + V (r) + F z − E ψ(r) = 0. (1) The potential V (r) representing the interaction with the molecular ion implicitly depends on its shape determined by the internuclear configuration and orientation with respect to

the field. The only assumption regarding the potential is Dn Z − 2 + O(r −3 ), (2) r r where Z is the total charge of the ion, D is its dipole moment, and n = r/r. We solve Eq. (1) in parabolic coordinates defined by [15] V (r)|r→∞ = −

ξ = r + z, 0 ≤ ξ < ∞, η = r − z, 0 ≤ η < ∞, y ϕ = arctan , 0 ≤ ϕ < 2π. x The solutions are sought in the form fν (η) ν (ξ,ϕ), ψ(r) = η−1/2

(3a) (3b) (3c)

(4)

ν

where eimϕ ν (ξ,ϕ) = φν (ξ ) √ 2π

(5)

and φν (ξ ) are the solutions of the eigenvalue problem d d m2 Eξ Fξ2 ξ − +Z+ − − βν φν (ξ ) = 0, (6a) dξ dξ 4ξ 2 4 φν (ξ )|ξ →0 ∝ ξ |m|/2 , φν (ξ )|ξ →∞ = 0, ∞ φnξ m (ξ )φn ξ m (ξ ) dξ = δnξ n ξ .

(6b) (6c)

0

The parabolic channel functions (5) are identified by the multiindex ν = (nξ ,m), nξ = 0,1, . . . , m = 0, ± 1, . . . ,

(7)

where m is the azimuthal quantum number and nξ enumerates the different solutions to Eqs. (6) for a given m. In the numerical solution of Eq. (1) [12–14] it is more efficient to use an adiabatic (dependent on η) basis in the expansion (4), but in constructing the asymptotic solution for F → 0 [19,25] it is more convenient to use the diabatic (independent of η) basis (5). The functions (5) are orthonormal with respect to the inner product ∞ 2π ei(m −m)ϕ dξ dϕ = δνν , (8)

ν | ν ≡ φν (ξ )φν (ξ ) 2π 0 0 where ν = (n ξ ,m ). Substituting Eq. (4) into Eq. (1), one can obtain a set of coupled equations defining the unknown coefficient functions fν (η). The SSs are given by the solutions to these equations satisfying the regularity boundary conditions at η = 0 and the outgoing-wave boundary conditions at η → ∞, 1/2 3/2 iEη1/2 21/2 fν iF η . (9) + fν (η)|η→∞ = exp (F η)1/4 3 F 1/2 The solutions exist only for a discrete set of generally complex values of E, so this is an eigenvalue problem. The real-energy eigensolutions of Eq. (1) for F = 0 representing bound states of the unperturbed system can be constructed in the same way, but with zero asymptotic boundary condition fν (η)|η→∞ = 0 instead of Eq. (9). We assume that the eigenfunctions in this case are real, which is always possible to achieve for bound

063410-2

WEAK-FIELD ASYMPTOTIC THEORY OF TUNNELING . . .


states. For F > 0, the bound states turn into the SSs. The SS eigenvalues presented in the form i E=E− (10) 2 define the energy E and ionization rate of the state. The SS eigenfunctions are normalized by 1 ∞ ∞ 2π 2 2 ψ (r)(ξ + η)dξ dηdϕ = 1, ψ (r)dr = 4 0 0 0 (11) where the integral should be regularized by deforming the integration path in η from the real semiaxis into a contour in the complex η plane [12]. The normalization condition (11) defines the coefficients fν in Eq. (9) up to a common sign. These coefficients, together with the parabolic channel functions (5), define the TMD of the ionized electrons [19] 2

2 k⊥ 8π 2 ,ϕk , P (k⊥ ) = fν ν (12) F ν F where k⊥ = (kx ,ky ) = (k⊥ cos ϕk ,k⊥ sin ϕk ) is the transverse momentum. The energy eigenvalue E and the asymptotic coefficients fν are the main characteristics of a SS related to observables. An efficient numerical procedure to calculate these quantities in a wide interval of fields was developed in Ref. [12] for axially symmetric potentials and in Refs. [13,14] for arbitrary potentials without any symmetry.

combinations of the two degenerate states [15]. We choose them to be real. In our geometry, one of them is even (+) and the other is odd (−) with respect to the reflection y → −y, 1 + (r ) = √ [ψn|M| (r ) + ψn−|M| (r )] ψn|M| 2 cos |M|ϕ = fn|M| (r ,θ ) √ , π 1 − (r ) = √ [ψn|M| (r ) − ψn−|M| (r )] ψn|M| i 2 sin |M|ϕ = fn|M| (r ,θ ) √ . π

(14a)

(14b)

In the calculations below, we consider only states with M = 0 (σ states) and |M| = 1 (π states). The σ states are even and do not have nodes in ϕ . The even and odd π states are denoted by π + and π − , respectively; they have nodal planes which for β = 0 coincide with the yz and xz planes, respectively. The spatial structure of the σ and π ± states is illustrated in Sec. IV. Let E0 and ψ0 (r ) denote the energy and wave function of the unperturbed initial bound state in the molecular frame. By applying the standard perturbation theory [15], in the presence of a nonzero field for this state we obtain E = E0 − μz F − 12 αzz F 2 + O(F 3 ), ψ(r ) = ψ0 (r ) + ψ1 (r )F + O(F 2 ),

(15a) (15b)

where

III. WEAK-FIELD APPROXIMATIONS

μz = μz cos β,

In the weak-field limit, F → 0, the energy E of a given SS can be obtained from perturbation theory [15], while the ionization rate , coefficients fν , and TMD P (k⊥ ) can be found using the WFAT [19,25]. In this section we summarize formulas needed to implement these approximations.

(16a)

αzz = αx sin2 β + αz cos2 β,

ψ1 (r ) = −χx (r ) sin β + χz (r ) cos β,

(16b) (16c)

and A. Perturbation theory

In the calculations below we consider diatomic molecules. ˆ = (x ,y ,z ) ≡ (x ,x ,x ) denote Let r = (x,y,z) and r = Rr 1 2 3 the Cartesian coordinates of the active electron in the laboratory and molecular frames, respectively, where Rˆ is an Euler rotation [26] from the laboratory to the molecular frame. By our convention, the z axis is directed along the electric field, the internuclear axis z lies in the xz plane of the laboratory frame, and the y and y axes coincide. Then the different orientations of the molecule with respect to the field are described by a single angle β, 0 ≤ β ≤ π , defining the rotation Rˆ from z to z about the y = y axis. For diatomic molecules the potential V (r) is axially symmetric about the internuclear axis z . The unperturbed (F = 0) solutions to Eq. (1) can be characterized by the projection M of the electronic angular momentum onto this axis. Such solutions as functions of spherical coordinates in the molecular frame, r = (r ,θ ,ϕ ), have the form

eiMϕ ψnM (r ) = fn|M| (r ,θ ) √ . 2π

(13)

The states with M = 0 are doubly degenerate since their energy En|M| does not depend on the sign of M. The correct zeroth-order wave functions are given by certain linear

μxi = − ψ0 |xi |ψ0 , αxi = −2 χxi (r ) =

| ψnM |x |ψ0 |2 i , E − E 0 n|M| nM=0

ψnM |x |ψ0 i ψnM (r ). E − E 0 n|M| nM=0

(17a) (17b)

(17c)

In Eqs. (16a) and (16b), we have taken into account that the components of the permanent dipole moment vector perpendicular to the internuclear axis vanish, μx = μy = 0, and the static dipole polarizability tensor in the molecular frame is diagonal, αxi xj = αxi δij , where αxi is the polarizability in the direction of the xi axis. The summations in Eqs. (17b) and (17c) run over the complete spectrum of the unperturbed Hamiltonian excluding the initial state. More specific forms of these formulas for σ and π ± states and some details of their implementation in the present calculations are given in the Appendix. Note that while for implementing the leading-order WFAT [19] one needs only the unperturbed energy E0 and wave function ψ0 (r ) of the initial state and its dipole moment μxi , to implement the first-order correction terms [25] one additionally needs the polarizabilities αxi and the first-order

063410-3



distortion (16c) of the wave function. Thus, all terms retained in expansions (15) are needed for the following.

and (0)

gν = η1/2−βν

The WFAT for the one-electron problem was developed in Ref. [19]. In this work, the leading-order terms in the asymptotic expansions of fν and for F → 0 were obtained. The first-order correction terms in the expansions were derived in Ref. [25]. Here we summarize formulas used in the present calculations. The parabolic channel eigenvalues and eigenfunctions defined by Eqs. (5)–(8) can be found using perturbation theory, (18a)

(1) 2 ν (ξ,ϕ) = (0) ν (ξ,ϕ) + ν (ξ,ϕ)F + O(F ). (18b)

The leading-order terms in these expansions are given by [19]

|m| + 1 (0) , (19a) βν = Z − nξ + 2 nξ ! (0) 1/2 |m|/2 −ξ/2 L(|m|) (ξ ), (19b) φν (ξ ) = (ξ ) e (nξ + |m|)! nξ where =

2|E0 |,

(20)

and L(α) n (x) are the generalized Laguerre polynomials [27]. The function φν(0) (ξ ) should be substituted into Eq. (5) for φν (ξ ) to obtain (0) ν (ξ,ϕ). The explicit expressions for the first-order terms in Eqs. (18) are given in Ref. [25]. The asymptotic expansion for the coefficient fν in Eq. (9) has the form F 1 (0) fν = fν 1 + Aν F ln 2 2 4 1 + (Bν + 2i Im aν − iπ Aν )F + O(F 2 ) , (21) 2 where the leading-order term is given by [19]

(0) 1/2 gν 4 2 βν / iπ iπβν(0) 3 (0) + − μz − , fν = 1/2 exp 2 F 4 3F

− (C2 η2 + C1 η + Cl ln η + C0 ),

(25a)

2βν(1) 3β (0)2 γm − 2μz βν(0) − ν5 , (23a) − 3 μ2 μz 4μz βν(0) 10 + 18γm + 3γm2 − Bν = −αzz − z + 2 + 3 24 3 (0) (0)2 (9 − 6γm )βν (49 + 2γm )βν − − 4 4 8 5 (0)3 (0)4 3β β + ν 6 − ν 7 + 2 Re aν . (23b) 2 8

Aν = −

Here we use the notation (24)

(25b)

where Rˆ is the Euler rotation defined above. The coefficients C2 , . . . ,C0 appearing in Eq. (25b) are given in Ref. [25]. The components of the dipole moment μz and polarizability αzz in Eq. (23b) refer to the laboratory frame and hence depend on the orientation of the molecule; they are related to the components in the molecular frame given by Eqs. (17a) and (17b) by the rotation Rˆ between the frames. The techniques to extract the asymptotic coefficients gν [Eq. (25a)] from an unperturbed molecular orbital were developed for linear [28,29] and nonlinear [30] molecules. To obtain the coefficients aν [Eq. (25b)], one needs to know the first-order distortion (16c) of the initial orbital. The techniques to calculate aν for spherically symmetric atomic potentials were developed in Ref. [25]; in the present work these techniques are generalized to molecular potentials. Within the WFAT, the total ionization rate defined by Eq. (10) is given by

=

ν + O( 2 ), ν = |fν |2 , (26) ν

where ν is the partial rate for ionization into parabolic channel ν. Substituting here Eq. (21), we obtain the asymptotic expansion F 2 2 2

ν = |Gν | Wν (F ) 1+Aν F ln 2 + Bν F +O(F ln F ) , 4 (27) where Gν and Wν (F ) are the structure and field factors given by Gν = e−μz gν ,

(28)

4 2 2Z/−2nξ −|m|−1 2 3 . (29) exp − Wν (F ) = 2 F 3F

and the coefficients of the first-order terms are [25]

1 − m2 − 2Dz 4

ˆ (0) ν ψ0 (Rr)η→∞ ,

aν = gν−1 η eη/2 (1)

ˆ + (0) ˆ × ν ψ0 (Rr) ν ψ1 (Rr) η→∞

(22)

γm =

e

1/2−βν(0) /

B. Weak-field asymptotic theory

βν = βν(0) + βν(1) F + O(F 2 ),

/ η/2

In principle, expansion (21) allows one to retain in Eq. (27) also terms of orders O(F 2 ln2 F ) and O(F 2 ln F ). However, we share a view that logarithm is a “constant” rather than a function (meaning that it is a too slowly varying function), and hence these terms can be retained only simultaneously with the term O(F 2 ), which is not available in Eq. (21). Note that the dependence of ν on the orientation of the molecule is contained in the structure factor Gν and correction coefficients Aν and Bν , while the field factor Wν (F ) does not depend on the orientation. Substituting Eq. (27) into Eq. (26), in the leading-order approximation one should retain only the contribution from the dominant ionization channel corresponding to the lowest power of F in Eq. (29). The same applies to the asymptotic expansion of the TMD, which can be

063410-4



obtained by substituting Eqs. (18b) and (21) into Eq. (12). We give here only the final formulas needed for the following. In the leading-order approximation, the total ionization rate and TMD for any state are given by [19]

as = (2 − δm0 )|G0m |2 W0m (F )

(30)

and 4π (0) 2 φ (s) Pas (k⊥ ) = (2 − δm0 ) as F 0m 2 cos (mϕk ), even states, × sin2 (mϕk ), odd states,

(31)

where m is the azimuthal quantum number of the dominant ionization channel and s=

2 k⊥ . F

(32)

Here we have taken into account that ν does not depend on the sign of m, and hence there are two identical terms in Eq. (26) if |m| > 0. Equation (31) for even states with m = 0 coincides with the well-known result obtained in Refs. [31],

k 2 4π Pas (k⊥ ) = as exp − ⊥ . (33) F F In the first-order approximation, when the correction terms in Eq. (21) are retained, one should additionally include into the expansions for and P (k⊥ ) the leading-order terms from the next-to-the-dominant channels, since they have the same power of F ; for more details on this point, see Ref. [25]. The formulas become more involved and we give them only for the states considered below. For σ states at any orientation angle β and π + states at β > 0, the dominant channel is ν = (0,0). The ionization rate and TMD in this case are given by

F F

= W00 (F ) G200 1 + A00 F ln 2 + B00 F + 2 |G01 |2 , 4 2 and P (k⊥ ) =

(34)

4π F 1 W00 (F )e−s G200 1 + A00 F ln 2 + B00 + 3 [6 + 4μz − 4(μz + 1)s − s 2 ] F F 4 4 √ √ F F 3/2 . − {(1 − s)G00 G10 − 2s[Re(G01 eiϕk )]2 + 2sG00 Re(G02 e2iϕk )} 2 − sG00 (2 Im a00 − π A00 ) Re(G01 eiϕk ) 2 (35)

−

+

For π states at any β and π states at β = 0, the dominant channels are ν = (0, ± 1). In this case we have

F F ,

= W01 (F ) 2 |G01 |2 1 + A01 F ln 2 + B01 F + 2 |G02 |2 + G210 4 4 2 and P (k⊥ ) =

(36)

4π F 1 W01 (F )e−s 4s[Re(G01 eiϕk )]2 (1 + A01 F ln 2 + B01 + 3 [2(4μz + 9) − 2(2μz + 3)s − s 2 ] F F 4 4 √ √ F + {(1 − s)G10 − 4s Re(G01 eiϕk ){24 2 Im(G01 eiϕk ) Im(a01 ) + 3 2(2 − s) Re(G11 eiϕk ) + 6s Re(G03 e3iϕk )} 12 2 3/2 √ √ √ 2iϕk 2 F iϕk 2iϕk F . (37) + 2s Re(G02 e )} + 2π A01 s Re(G01 e ){(1 − s)G10 + 2s Re(G02 e )} 4 2 2

Here we have used that the coefficients Gnξ 0 are real, because the unperturbed bound-state wave function is assumed to be real. Equations (34)–(37) coincide with Eqs. (30) and (31) if one drops the first-order terms. The asymptotic expansions given above hold for F < Fc , where Fc is a boundary between the tunneling and overthe-barrier ionization regimes. This field can be estimated as [19]

Fc ≈

4

. 8|2Z − (|m| + 1)|

(38)

The condition F Fc guarantees that the first-order terms in the expansions are much smaller than the leading-order terms.

IV. RESULTS AND DISCUSSION

In this section, we illustrate the role of the first-order correction terms in applications of the WFAT to molecular systems. We consider one-electron diatomic molecular ions described by the soft-core potential V (r) = −

Z1 (r − r1 )2 +

−

Z2 (r − r2 )2 +

,

(39)

where Zi and ri = (zi sin β,0,zi cos β), i = 1,2, are the charges and positions of the nuclei in the laboratory frame, and is the softening parameter. The softening of the Coulomb singularities of the potential is not essential for the WFAT, but is dictated by limitations of the current implementation of our accurate computational method [13,14]. This potential satisfies Eq. (2) with Z = Z1 + Z2 and D = Z1 r1 + Z2 r2 . We

063410-5



A. Validation of the WFAT by comparison with accurate calculations of Siegert states in H2 +

For H2 + , the potential (39) does not change under a reflection z → −z in the molecular frame. Because of this symmetry, all the observables are invariant under the transformation β → π − β, so it is sufficient to consider orientations in the interval 0◦ β 90◦ . For all states of this nonpolar molecule μz = 0. Hence, Aν does not depend on β, while Bν depends on β only through αzz and aν ; see Eqs. (23). We consider the same three lowest states 1sσ and 2pπ ± as in Refs. [13,14]. A comparison of the exact results for the energy E of these states with the predictions of perturbation theory [Eq. (15a)] was reported in Ref. [13], so we do not discuss the behavior of E here. 1. The ground-state H2 + (1sσ )

The field-free energy of the ground 1sσ state in the present soft-core model E0 = −0.962 366 is slightly higher than the corresponding energy −1.102 634 for the pure Coulomb potential with = 0. The polarizabilities calculated using Eqs. (A2a) and (A2b) are αx = 2.8775 and αz = 5.9095, respectively [32]. At all orientations β, the dominant ionization channel is ν = (0,0). The correction coefficients (23) for this channel were calculated using Eqs. (A2c) and (A2d). We have A00 = −0.715; the dependence of B00 on β is shown in the top

1sσ, ν = (0, 0)

-7.5 -8.0 -8.5

Bν (a.u.)

consider in detail several states of a homonuclear system H2 + , but also discuss briefly a heteronuclear system HeH2+ with a permanent dipole moment (Sec. IV B 2). All the calculations for H2 + were performed with Z1 = Z2 = 1, z1 = 1, and z2 = −1; hence, Z = 2 and D = 0, which corresponds to the equilibrium internuclear distance R = 2 for the ground 1sσ state, and = 0.09. With these parameters, Eq. (39) coincides with the potential used in Refs. [13,14]. The results for HeH2+ were obtained with Z1 = 1, Z2 = 2, z1 = 3.112, and z2 = −0.778, which corresponds to the equilibrium internuclear distance R = 3.89 for the lowest bonding 2pσ state, and = 0.1. In this case Z = 3 and D = 1.556ez . The same potential was used in Ref. [19]. In Sec. IV A, we present calculations for three states of different symmetry of H2 + . The WFAT results for the asymptotic coefficients fν , ionization rate , and TMD P (k⊥ ) as functions of the field F and orientation angle β obtained by implementing equations of Sec. III are compared with accurate results obtained by solving Eq. (1) using the method and program developed in Refs. [12–14]. The goal of these calculations is to validate the WFAT and demonstrate the effect of the first-order terms on the different observables. In Sec. IV B, we discuss two situations where the first-order terms not only improve quantitative agreement with the accurate results, but are essential for their qualitative understanding. Finally, in Sec. IV C, we give examples of the comparison of the predictions of the WFAT with the results of time-dependent calculations available in the literature. In the following, for brevity, the leading-order WFAT results and those including the first-order correction terms are denoted by WFAT(0) and WFAT(1), respectively, and the results of accurate numerical calculations are referred to as “exact.”

-36.65 -36.70

2pπ − , ν = (0, 1)

-36.75 -40 -60 -80 -100

2pπ + , ν = (0, 0) 0

30

β (deg)

60

90

FIG. 1. Correction coefficients Bν [Eq. (23b)] for the dominant ionization channel in each of the three states of H2 + considered in Sec. IV A as functions of the orientation angle β.

panel of Fig. 1. The boundary of the over-the-barrier ionization regime estimated from Eq. (38) is Fc = 0.18. We first consider the asymptotic coefficient f00 for the dominant ionization channel defined by Eq. (9). Figure 2 illustrates the behavior of this coefficient as a function of the field F at three representative orientations β. The exact results (solid lines) are compared with the WFAT(1) results (dashed lines) obtained from Eq. (21). In order to eliminate a rapid variation of f00 by many orders of magnitude in the interval of F considered, which is caused by the last term in the exponent in Eq. (22), and thus to facilitate the comparison, we plot the (0) . The upper (black lines) and lower (red lines) ratio f00 /f00 halves of each panel in Fig. 2 show the real and imaginary parts of this ratio, respectively. At all orientations, the ratio obviously approaches 1 as F → 0. For the WFAT(1) results this follows immediately from Eq. (21). The exact results stop at some nonzero F , where /E ∼ 10−10 . The present calculations based on double-precision arithmetics cannot be continued to smaller F , where and f00 become too small, because of the round-off errors. However, by extrapolating the exact results to F = 0 one can expect that the ratio approaches unity, which confirms the prediction of the WFAT. Moreover, the difference between the exact and WFAT(1) results is seen to decrease ∝ F 2 as F → 0, as it should be, according to Eq. (21). Note that this holds for both the real and the imaginary parts of the ratio, and hence Eq. (21) describes correctly not only the absolute value of f00 , but also its phase. Because of the symmetry of the 1sσ state, it depends on β, for which of the higher channels fν = 0. For example, at β = 0◦ , the SS eigenfunction is axially symmetric about the z axis, and therefore all fν with m = 0 vanish. At β = 90◦ , the eigenfunction is even with respect to the reflection x → −x,

063410-6



β = 0◦

β = 45◦

β = 90◦

1.0

f00 /f00

(0)

0.8 0.6

exact WFAT(1)

0.4

exact WFAT(1)

0.2 0.1 0.0 0

0.05

0.10

0.15

0.20 0

0.05

0.10

0.15

0.20 0

0.05

0.10

0.15

0.20

F (a.u.) FIG. 2. (Color online) Ratios of the asymptotic coefficient f00 in Eq. (9) for the dominant ionization channel in H2 + (1sσ ) taken at three representative orientation angles β to its leading-order WFAT value f00(0) [Eq. (22)] as functions of the field. Solid lines, exact results; dashed lines, WFAT(1), the WFAT results including the first-order correction terms [Eq. (21)]. The upper (black lines) and lower (red lines) halves of each panel show the real and imaginary parts of the ratio, respectively.

and hence fν have nonzero values only for channels with even m. In the interval of F considered, the coefficient fν for the next-to-the-dominant channel is more than four orders of magnitude smaller than f00 . We next consider the ionization rate. In the top row of Fig. 3, the exact results for as a function of F at the same three orientations as in Fig. 2 are compared with the WFAT(0) and WFAT(1) results obtained from Eq. (30) with m = 0 and Eq. (34), respectively. The WFAT(0) consistently overestimates the rate, while the WFAT(1) is much closer to the exact results up to a field ∼Fc , where the right-hand side of Eq. (34) becomes negative. To eliminate the variation of the rate by many orders of magnitude and thus facilitate comparison of the different results on a linear scale, in the

Γ (a.u.)

β = 0◦ 100 10-2 10-4 10-6 10-8 10-10

bottom row of the figure we show the ratio / as . Because of the round-off errors mentioned above, the exact calculations cannot be continued to smaller F . However, even the available results unambiguously confirm Eq. (34). The correction terms in this equation, not included in Eq. (30), essentially improve the agreement with the exact results throughout the tunneling regime F < Fc , where the WFAT is expected to apply. Note that at all β the ratio approaches unity from below as F → 0. This is explained by the fact that the term with B00 dominates over the other correction terms in Eq. (34) (the term with A00 becomes dominant only at very small F ), and this term is negative; see the top panel of Fig. 1. The dependence of the ionization rate on the orientation angle β for several representative values of F is shown in

β = 45◦

β = 90◦

exact W FAT(0) W FAT(1)

1.0

Γ/Γas

0.8 0.6 0.4 0.2 0.0 0

0.1

0.2

0

0.1

0.2

0

0.1

0.2

F (a.u.) FIG. 3. (Color online) (Top panels) Ionization rates of H2 + (1sσ ) at three representative orientation angles β as functions of the field. (Bottom panels) Ratios of the rates to their leading-order WFAT values as . Solid (black) lines, exact results; dotted (blue) lines, WFAT(0), the leading-order WFAT results [Eq. (30) with m = 0]; dashed (red) lines, WFAT(1), the WFAT results including the first-order correction terms [Eq. (34)]. 063410-7


Z

Field F = 0.05

4.5

6x10-9

4.0

4x10-9

3.5

2x10-9

3.0

y F = 0.09

2.5 F = 0.11

x

x

2.0 1.5 1.0

0 10 20 30 40 50 60 70 80 90

Orientation angle β (deg)

FIG. 4. (Color online) (Left) Illustration of the unperturbed wave function of H2 + (1sσ ). The molecular frame (x ,y ,z ) is obtained from the laboratory frame (x,y,z) by a rotation through the angle β about the y = y axis. The dominant ionization channel for this state is ν = (0,0), independently of β. (Right) Ionization rates divided by the field factor W00 (F ) [Eq. (29)] as functions of the orientation angle β for three representative values of the field. Solid (black) lines, exact results; dotted (blue) line, WFAT(0), the leading-order WFAT results [Eq. (30) with m = 0] (in this approximation the ratio is |G00 |2 and does not depend on F ); dashed (red) lines, WFAT(1), the WFAT results including the first-order correction terms [Eq. (34)].

P (k⊥ , ϕk = 0) (a.u.)

exact WFAT(0) WFAT(1)

β

Γ/W00 (F)

Z


β = 0◦


F = 0.05

0 4 3 2 1 0 4x10-9

F = 0.15

β = 45◦

2x10-9 F = 0.05

0 3 2 1 0

F = 0.15

3x10-9

β = 90◦

2x10-9 1x10-9

F = 0.05

0 2 1

Fig. 4. To bring the results for the different fields to the same scale, we plot the rate divided by the field factor (29) for the dominant ionization channel. As follows from Eq. (30), within the WFAT(0) the ratio is given by the structure factor G00 squared. For the present nonpolar molecule μz = 0 and G00 coincides with the coefficient g00 in the asymptotic tail of the unperturbed bound-state wave function defined by Eq. (25a). Thus, the WFAT(0) results for the ratio present a property of the ionizing orbital and do not depend on F . The exact results should approach the WFAT(0) curve as F → 0. This is indeed the case. The WFAT(1) results do depend on F and are in much better agreement with the exact results. The WFAT(0) correctly reproduces the shape of the dependence of on β, and hence this dependence is determined by the asymptotic coefficient g00 . This results from the fact that the correction terms in Eq. (34) for the present case only weakly depend on β. The WFAT(1) incorporating these terms essentially improves the magnitude of at all orientations. For all values of F considered, the rate maximizes at β = 0◦ and minimizes at β = 90◦ ; that is, the electron is more easily ionized in parallel than in perpendicular orientation. This conclusion agrees with the results of the time-dependent calculations of the ionization yield from H2 + (1sσ ) in intense laser fields [33–36]. We now turn to the discussion of the TMD. We consider P (k⊥ ) defined by Eq. (12) as a function of polar coordinates (k⊥ ,ϕk ) in the plane of the transverse momentum k⊥ . For the 1sσ state, the TMD is exactly isotropic (that is, does not depend on ϕk ) at β = 0◦ . According to the WFAT, it should remain almost isotropic at β = 0◦ , because the dominant ionization channel has m = 0. Figure 5 shows the cuts of the TMD along the ray ϕk = 0 for several representative values of β and F . The WFAT(0) and WFAT(1) results are obtained from Eq. (31) for even states with m = 0 and Eq. (35), respectively. As follows

F = 0.15

0 0

0.2

0.4

0.6

0.8

k⊥ (a.u.) FIG. 5. (Color online) Cuts of the TMDs [Eq. (12)] for H2 + (1sσ ) at several representative orientation angles β and fields F along the ray ϕk = 0. Solid (black) lines, exact results; dotted (blue) lines, WFAT(0), the leading-order WFAT results [Eq. (31) with m = 0]; dashed (red) lines, WFAT(1), the WFAT results including the firstorder correction terms [Eq. (35)].

from Eq. (33), in the weak-field limit the TMD as a function of k⊥ should have a Gaussian shape with the width ∝ F 1/2 . This is qualitatively confirmed by the exact results; in particular, the width of the TMD obviously grows with F . However, the magnitudes of the WFAT(0) and exact results are different, and this difference grows with F . The correction terms included in Eq. (35) make the agreement with the exact results much better. A more subtle feature is the appearance of an anisotropy of the TMD (that is, its dependence on ϕk ) at β = 0◦ , which is not accounted for by the leading-order approximation (33). To illustrate this feature, we present in Fig. 6 the cuts of the TMDs shown in Fig. 5 (except for the ones at β = 0◦ , which are isotropic) along the circle k⊥ = 0.1. To emphasize the anisotropic part, we subtract from P (k⊥ ,ϕk ) its average over ϕk defined by 2π 1 ¯ P (k⊥ ) = P (k⊥ ,ϕk )dϕk . (40) 2π 0 Figure 6 compares the exact and WFAT(1) results for the difference. At β = 45◦ , the dominant correction to Eq. (33)

063410-8


P (k⊥ = 0.1, ϕk ) − P¯ (k⊥ = 0.1) (a.u.)

1x10-11

exact W FAT(1)

-1x10-11 0.02

corrections in Eq. (35) depend on ϕk as cos 2ϕk , which is again in accordance with the exact results. Thus, Eq. (35) correctly describes the shape of the anisotropic part of the TMD and relates it to the symmetry of the state. It also reproduces its magnitude; as expected, the agreement with the exact results becomes worse as F grows.

β = 45◦

F = 0.05

0


F = 0.15

2. Excited states H2 + (2 pπ ± )

0 -0.02 5x10-12

β = 90◦

F = 0.05

0 -5x10-12 0.005

F = 0.15

0 -0.005

-π

-π/2

0

ϕk (rad)

π/2

π

FIG. 6. (Color online) Cuts of the TMDs shown in Fig. 5 along the circle k⊥ = 0.1. To emphasize the anisotropic part of the TMDs, their average P¯ (k⊥ ) over ϕk [Eq. (40)] is subtracted. Solid (black) lines, exact results; dashed (red) lines, WFAT(1), the WFAT results including the first-order correction terms [Eq. (35)]. For this state, the exact TMD at β = 0◦ is isotropic for all fields. Within the WFAT(0), the same holds at any β.

Γ (a.u.)

at the value of k⊥ considered comes from the last term in Eq. (35). This term is ∝ cos ϕk , which explains the shape of the dependence of P (k⊥ ,ϕk ) on ϕk in this case. However, the coefficient G01 in this term vanishes at β = 90◦ . All the other

β = 0◦

100 10-2 10-4 10-6 10-8 10-10 10-12

Here we consider two excited states 2pπ ± with |M| = 1. They are degenerate in the absence of the field. Their field-free energy in the present soft-core model is E0 = −0.418 947 [37], while in the pure Coulomb potential with = 0 it is −0.428 772. It is convenient to begin the discussion with the odd 2pπ − state, because in this case the same ionization channel ν = (0,1) remains dominant at all orientation angles β. The polarizabilities in this state calculated using Eqs. (A4a) and (A4b) are αx = 19.2303 and αz = 23.4056, respectively [32]. The correction coefficients (23) were calculated using Eqs. (A4c) and (A4d). We have A01 = −1.580; the dependence of B01 on β is shown in the middle panel of Fig. 1. The critical field estimated from Eq. (38) is Fc = 0.040. Figure 7 shows the ionization rate of H2 + (2pπ − ) as a function of F for the same orientations as in Fig. 3. The WFAT(0) and WFAT(1) results are obtained from Eq. (30) with m = 1 and Eq. (36), respectively. The results look quite similar to those in Fig. 3. At all orientations, the WFAT(0) overestimates the rate while the WFAT(1) is in much better agreement with the exact results up to a field ∼Fc , where the right-hand side of Eq. (36) becomes negative. The bottom row of the figure presents the ratio / as on a linear scale. It is clearly seen that the first-order correction terms in Eq. (36) essentially improve the WFAT results. Figure 8 presents the dependence of the ionization rate on β for several values of F . In contrast to the 1sσ state, which is axially symmetric about the internuclear z axis (see

β = 45◦

β = 90◦


1.0

Γ/Γas

0.8 0.6 0.4 0.2 0.0 0

0.01

0.02

0.03 0.04 0

0.01 0.02

0.03

0.04 0

0.01 0.02

0.03

0.04

F (a.u.) FIG. 7. (Color online) Same as in Fig. 3, but for the 2pπ − state of H2 + . WFAT(0), the leading-order WFAT results [Eq. (30) with m = 1]; WFAT(1), the WFAT results including the first-order correction terms [Eq. (36)]. 063410-9



z 0.6

10x10-8

0.5

5x10-8

z

Field F = 0.015

y

H


0.4

F = 0.020

0.3

F = 0.030

0.2

Γ/W 01 (F )

β

x 0.1

x

0 10 20 30 40 50 60 70 80 90

Orientation angle β (deg)

FIG. 8. (Color online) Same as in Fig. 4, but for the 2pπ − state of H2 + . The dominant ionization channel for this state is ν = (0,1), independently of β. WFAT(0), the leading-order WFAT results [Eq. (30) with m = 1] (in this approximation the ratio is 2|G01 |2 and does not depend on F ); WFAT(1), the WFAT results including the first-order correction terms [Eq. (36)].

the left panel in Fig. 4), the 2pπ − state has a node in the plane xz including in the direction of the field. This results in a suppression of tunneling ionization from this state; within the WFAT, this suppression is accounted for by an additional power of F in the field factor (29) for the dominant ionization channel. The orientation dependence of the ratio /W01 (F ) is flatter than that for the 1sσ state (see the right panel in Fig. 4). This is explained by the fact that the unperturbed 2pπ − wave function is almost axially symmetric about the y axis, and therefore all the coefficients in Eq. (36) only weakly depend on β. The WFAT(1) results more rapidly depart from the exact results as F grows, but still give more accurate predictions for the ionization rate than the WFAT(0) results. For the 2pπ − state, P (k⊥ ,ϕk ) turns to zero at ϕk = 0 and ±π . Figure 9 shows the cuts of the TMD along the ray ϕk = π/2 for several values of β and F . The WFAT(0) and WFAT(1) results are obtained from Eq. (31) for odd states with m = 1 and Eq. (37), respectively. Although the shape of the cuts as functions of k⊥ is different from that for the 1sσ state, the situation in Fig. 9 is similar to that in Fig. 5. At all orientations and fields considered, the WFAT(0) correctly describes the shape but overestimates the magnitude of the TMD, while the WFAT(1) is in much better agreement with the exact results over a wide interval of F . The dependence of P (k⊥ ,ϕk ) on ϕk is illustrated in Fig. 10. Here we present the cuts of the TMDs shown in Fig. 9 along the circle k⊥ = 0.2. As seen from Eq. (31), for the present state the TMD is anisotropic even in the leading-order approximation of the WFAT, so there is no need to subtract its average over ϕk as we did above for the 1sσ state. In all the cases, the first-order corrections in Eq. (37) essentially improve the agreement between the WFAT and exact results. We next consider the even 2pπ + state. At β = 0◦ for any F , this state coincides with the 2pπ − state rotated about the z axis by π/2. In this case the dominant ionization channel is ν = (0,1). The rotation does not affect the observables, apart from a shift by π/2 of the angular argument ϕk of the TMD P (k⊥ ,ϕk ). Thus, the results for the 2pπ + state at β = 0◦ coincide with the

P (k⊥ , ϕk = π/2) (a.u.)

H


β = 0◦ F = 0.015

0 15x10-4 10x10-4 5x10-4

F = 0.022

0 10x10-8

β = 45◦

5x10-8 F = 0.015

0 10x10-4 5x10-4

F = 0.022

0 8x10-8 6x10-8 4x10-8 2x10-8

β = 90◦ F = 0.015

0 10x10-4 5x10-4

F = 0.022

0 0

0.1

0.2

0.3

0.4

0.5

k⊥ (a.u.) FIG. 9. (Color online) Same as in Fig. 5, but for the 2pπ − state of H2 + and the cuts are made along the ray ϕk = π/2. WFAT(0), the leading-order WFAT results [Eq. (31) with m = 1]; WFAT(1), the WFAT results including the first-order correction terms [Eq. (37)].

corresponding results for the 2pπ − state, and it is sufficient to consider nonparallel orientations with β > 0◦ . At any nonzero β, the dominant ionization channel is ν = (0,0) (the change of the dominant channel is discussed in Sec. IV B 1). From Eqs. (A3a) and (A3b) we obtain the polarizabilities αx = 93.6687 and αz = 23.4056 [32]. Using Eqs. (A3c) and (A3d) we find A00 = −10.127; the dependence of B00 on β is shown in the bottom panel of Fig. 1. The critical field estimated from Eq. (38) is Fc = 0.028. Figure 11 compares the exact and WFAT results for the ionization rate of H2 + (2pπ + ) as a function of F . The WFAT(1) curve for the ratio / as at β = 45◦ is seen to approach 1 from above as F → 0. This is explained by the fact that the term with Aν in Eq. (34) becomes dominant as F → 0, and this term is positive, because the coefficient A00 is negative. In fact, the same applies to all three states of H2 + considered here at all β, but this behavior is not seen in the bottom rows of Figs. 3 and 7 because the coefficients Aν for the dominant channel in those states have much smaller values than in the present case. Again, the WFAT(1) is always closer to the exact results than the WFAT(0), but the convergence at weak fields for the present state is somewhat slower than in the previous cases. The dependence of the ionization rate on β for several values of F is shown in Fig. 12. There are three main features to be noticed from the figure. First, the exact rates become

063410-10

WEAK-FIELD ASYMPTOTIC THEORY OF TUNNELING . . . F = 0.015

z

β = 0◦


β Field

2x10-8 0 10x10-4 5x10-4


z

0.15

F = 0.022

y

H

F = 0.015

H

P (k⊥ = 0.2, ϕk ) (a.u.)

◦

β = 45

F = 0.025

4x10-8

x

2x10-8 0 10x10-4

4x10-8

0.00 0 10 20 30 40 50 60 70 80 90

Orientation angle β (deg) F = 0.022

FIG. 12. (Color online) Same as in Figs. 4 and 8, but for the 2pπ + state of H2 + . The dominant ionization channel for this state is ν = (0,1) at β = 0◦ and (0,0) at β = 0◦ . WFAT(0), the leading-order WFAT results [Eq. (30) with m = 0] (in this approximation the ratio is |G00 |2 and does not depend on F ); WFAT(1), the WFAT results including the first-order correction terms [Eq. (34)].

5x10-4 0 6x10-8

0.05

x

F = 0.015

0.10

F = 0.020

0 6x10-8

0.20

Γ/ W00 (F )

6x10-8 4x10-8


β = 90◦

F = 0.015

2x10-8 0 10x10-4

F = 0.022

5x10-4 0

-π

-π/2

π/2

0

ϕk (rad)

π

FIG. 10. (Color online) Cuts of the TMDs shown in Fig. 9 along the circle k⊥ = 0.2. Solid (black) lines, exact results; dotted (blue) lines, WFAT(0), the leading-order WFAT results, [Eq. (31) with m = 1]; dashed (red) lines, WFAT(1), the WFAT results including the first-order correction terms [Eq. (37)].

Γ (a.u.)

β = 45◦ 100 10-2

β = 90◦


10-4 10-6 10-8

1.0

Γ/Γas

0.8 0.6 0.4 0.2 0.0 0

0.01

0.02

0.03

0.04 0

0.01

0.02

0.03

0.04

F (a.u.) FIG. 11. (Color online) Same as in Figs. 3 and 7, but for the 2pπ + state of H2 + . WFAT(0), the leading-order WFAT results [Eq. (30) with m = 0]; WFAT(1), the WFAT results including the first-order correction terms [Eq. (34)]. At β = 0◦ , the results for this state coincide with those shown in the left column of Fig. 7.

very small (but remain nonzero) at β = 0◦ . Within the WFAT, this is explained by the fact that, on the one hand, the partial ionization rate for the dominant channel (0,0) vanishes at β = 0◦ , because the structure factor G00 vanishes. On the other hand, the partial rate for the next-to-the-dominant channel (0,1), that has a nonzero value at β = 0◦ , is suppressed by an additional power of F in the field factor (29). The WFAT(1) results from Eq. (34) account for the contributions from both these channels and restore nonzero values of the rates. Second, the exact rates have a shallow minimum at β = 90◦ and a maximum at some intermediate β whose position shifts towards 90◦ as F decreases. Such a maximum in the ionization yield was also found in time-dependent calculations for the present system in an intense laser field [33]. The WFAT(0) result for the ratio /W00 (F ) is |G00 |2 and does not depend on F . This curve has a maximum at β = 90◦ . The WFAT(1) results reproduce the position of the maximum of the exact rates and converge to the WFAT(0) curve as F → 0. The third feature is a rapid variation of the shape of the orientation dependence of the rate with F , much more rapid than in the previous cases. This is explained by a large value of the difference between the polarizabilities αx and αz in the present case, resulting in a strong dependence of B00 on β; see the bottom panel in Fig. 1. We finally consider the TMD for the 2pπ + state at sufficiently large β, where channel (0,0) is dominant. In this case, the TMD is almost isotropic, as for the 1sσ state. Figure 13 presents the cuts of the TMD along the ray ϕk = 0. To make the anisotropic part visible, we again subtract from P (k⊥ ,ϕk ) its average P¯ (k⊥ ) defined by Eq. (40). The cuts of the difference along the circle k⊥ = 0.2 are shown in Fig. 14. These results look similar to those in Figs. 5 and 6. To summarize this section, the results presented confirm that the WFAT quantitatively describes all the observables at sufficiently weak fields. The first-order correction terms are shown to essentially improve the agreement with the exact results and make the WFAT applicable up to the onset of over-the-barrier ionization.

063410-11

TRINH, PHAM, TOLSTIKHIN, AND MORISHITA 6x10-6

B. Case studies

β = 45◦


4x10-6


2x10-6

Having demonstrated the role of the first-order terms in the quantitative performance of the WFAT, here we give examples of situations where these terms are essential for understanding the results at a qualitative level.

F = 0.015

P (k⊥ , ϕk = 0) (a.u.)

0 0.06

1. A change of the dominant ionization channel in H2 + (2 pπ + ) at small β

0.04 0.02

F = 0.022

0 10x10-6

β = 90◦

5x10-6

F = 0.015 0 0.10 0.05

F = 0.022

0 0

0.1

0.2

0.3

0.4

k⊥ (a.u.) FIG. 13. (Color online) Same as in Figs. 5 and 9, but for the 2pπ + state of H2 + . The cuts are made along the ray ϕk = 0. WFAT(0), the leading-order WFAT results [Eq. (31) with m = 0]; WFAT(1), the WFAT results including the first-order correction terms [Eq. (35)].

exact W FAT(1)

5x10-8

A peculiarity of the 2pπ + state in comparison with the other two states discussed in Sec. IV A is that a change of the dominant ionization channel occurs at small β. This results from the fact that the structure factor G00 turns to zero at β = 0◦ , because of the symmetry of the state. As a consequence, for a given F , at sufficiently large β channel (0,0) is dominant, but at smaller β channel (0,1) becomes dominant. The change of the dominant channel occurs at a critical angle βc (F ) = O(F 1/2 ) [13,28], where the contributions from the two channels become comparable. The effect of the interplay of the contributions from these two channels on the shape of the TMD at β ≈ βc (F ) was discussed in Ref. [14]. Here we illustrate the effect on the behavior of the ionization rate. The left panel in Fig. 15 presents the leading-order contributions to the square brackets in Eq. (34) from channels (0,0) and (0,1) as functions of β. The former term does not depend on F ; the latter one is shown for a sufficiently weak field F = 0.02 < Fc . The two terms become equal at β ≈ 6◦ , which can be taken as the value of βc (F ) for F = 0.02. The right panel shows the same ratio / as as in the bottom row of Fig. 11, but at smaller values of β. One can see that in this case the ratio as a function of F behaves very differently. At a given β, it attains a maximum and then approaches 1 from above as F decreases. The WFAT(1) curve for β = 6◦ passes through the value of 2 at F ≈ 0.02, which agrees with the point where the

-5x10-8 2x10-3

β = 45◦

F = 0.015

4.0 0.030 0.025

0 -2x10-3

|G00 |2 F |G01 |2 2κ 2

F = 0.022


0.020

5x10-9

3.0

β = 5◦ β=6

◦

1.0

β = 90◦

F = 0.02

0.005

0.5

1x10-4

0.000 0

0 -1x10-4

-π

5

10

β (deg) F = 0.022

-π/2

0

π/2

1.5

β = 20◦

0.010 F = 0.015

2.5 2.0

0.015

0 -5x10-9

3.5

Γ/Γas

P (k⊥ = 0.2, ϕk ) − P¯ (k⊥ = 0.2) (a.u.)

0

π

ϕk (rad) FIG. 14. (Color online) Cuts of the TMDs shown in Fig. 13 along the circle k⊥ = 0.2. To emphasize the anisotropic part of the TMDs, their average P¯ (k⊥ ) over ϕk [Eq. (40)] is subtracted. Solid (black) lines, exact results; dashed (red) lines, WFAT(1), the WFAT results including the first-order correction terms [Eq. (35)].

15

20 0

0.01

0.02

0.03

0.0 0.04

F (a.u.)

FIG. 15. (Color online) (Left) Leading-order contributions to the square brackets in Eq. (34) for the ionization rate of H2 + (2pπ + ) from channels ν = (0,0) (solid black line) and ν = (0,1) for F = 0.02 (dashed green line) as functions of β. (Right) Ratio of the rate to its leading-order WFAT value as (as in the bottom row of Fig. 11, but at smaller β). Solid (black) lines, exact results; dotted (blue) line, WFAT(0), the leading-order WFAT results [Eq. (30) with m = 0]; dashed (red) lines, WFAT(1), the WFAT results including the firstorder correction terms [Eq. (34)].

063410-12



β = 0◦ E (a.u.)

-0.9 -1.0 exact PT

-1.1 -1.2 100 10-2 10-4 10-6 10-8 10-10 10-12 10-14

Γ/Γas

2. An avoided crossing in HeH2+ (2 pσ )

It is important to demonstrate the performance of the WFAT also for a polar molecule. To this end, we consider the lowest bonding state 2pσ of HeH2+ . Time-dependent calculations for this system in intense laser fields were reported in Ref. [38] and references therein. The field-free energy in the present soft-core model is E0 = −0.939 749, while for = 0 it is −1.045 349. The electronic dipole moment in the molecular frame is μz = −2.202 and the polarizabilities are αx = 5.6887 and αz = 12.4554. The dominant ionization channel at all β is (0,0) and the critical field is Fc = 0.095. We restrict our treatment of this molecule to considering the energy eigenvalue (10) at the parallel (β = 0◦ ) and antiparallel (180◦ ) orientations. Figure 16 compares the exact results with weak-field approximations. In the top row, the energy is compared with the predictions of perturbation theory [Eq. (15a)]. At β = 0◦ , the results are in close agreement with each other up to F ≈ 0.07, where the exact energy rapidly changes its slope and goes down at stronger fields. This is explained by the fact that the unperturbed state for F = 0 is localized on proton in the upper potential well near r = z1 ez (see Fig. 5 in Ref. [19]). As F grows, the energy of this state goes up, while the energy of the lower of the n = 2 parabolic states localized on α particle in the lower potential well near r = z2 ez goes down. The energies of the two states pass through an avoided crossing at F ≈ 0.07, which is seen in the top left panel in Fig. 16. This does not happen if the state is initially localized in the lower well, as is the case at β = 180◦ , and perturbation theory in this case works well over the whole interval of F considered. In the middle row, the exact results for the ionization rate are compared with the predictions of the WFAT [Eq. (30) with m = 0 and Eq. (34)]. To facilitate the comparison on a linear scale, the ratio / as is shown in the bottom row of the figure. The avoided crossing at F ≈ 0.07 is seen to manifest itself also in the field dependence of the ionization rate at β = 0◦ . At this orientation A00 = −6.9473. The term with this coefficient in Eq. (30) is the dominant correction at small F , which explains why the WFAT(1) curve in the bottom left panel attains a maximum and approaches 1 from above as F → 0. On the other hand, for β = 180◦ we have A00 = 0.8925. The positive value of this coefficient explains why the WFAT(1) curve in the bottom right panel approaches 1 from below as F decreases. The exact results confirm this difference in the behavior of the ratio / as at the two orientations. The difference was noticed already in

β = 180◦

-0.8

Γ (a.u.)

two curves in the left panel cross. This means that at this β the leading-order contribution from channel (0,1) in Eq. (34) dominates over the first-order corrections to the contribution from channel (0,0). As β grows, the relative role of the last term in Eq. (34) decreases, and we return to the situation where the term with B00 dominates and the ratio approaches 1 from below, as in Fig. 11. The exact results in Figs. 11 and 15 confirm this prediction of the WFAT. We emphasize that such an agreement between the exact and WFAT results would not be possible to achieve if one retained only the leading-order terms for each of the two channels in Eq. (34), omitting the first-order correction terms containing the coefficients A00 and B00 .


1.2 1.0 0.8 0.6 0.4 0.2 0.0 0

0.02 0.04 0.06 0.08 0.1 0

0.02 0.04 0.06 0.08 0.1

F (a.u.) FIG. 16. (Color online) Energy (top row), ionization rate (middle row), and the ratio / as (bottom row) for the 2pσ state of HeH2+ as functions of field at two orientations β. Solid (black) lines, exact results; dashed-dotted (green) lines, second-order perturbation theory results [Eq. (15a)]; dotted (blue) lines, WFAT(0), the leading-order WFAT results [Eq. (30) with m = 0]; dashed (red) lines, WFAT(1), the WFAT results including the first-order correction terms [Eq. (34)].

Ref. [19], but it could not be explained there without having the first-order correction terms in Eq. (34). It should be noted that the 2pσ state of HeH2+ is unstable against a nonadiabatic electronic transition to the 1sσ state and subsequent dissociation with the lifetime of several nanoseconds [39]. The field ionization rate has the same value at F ≈ 0.045 and 0.056 at β = 0◦ and 180◦ , respectively, which corresponds to typical amplitudes of intense laser fields of current interest. This means that theoretical simulations of the interaction of HeH2+ (2pσ ) with such fields must include the nonadiabatic electronic and nuclear dynamics into account. C. Comparison with time-dependent calculations

The orientation dependence of ionization yields from molecules in intense low-frequency laser fields is in the focus of many of the current studies in strong-field physics. It is generally believed that the ionization yield reflects the symmetry and spatial structure of the unperturbed ionizing orbital, which is based on the molecular Ammosov-Delone-Krainov (MO-ADK) formula [40] and justifies various molecular imaging techniques. Indeed, the ionization yields are often directly compared with the shape of molecular orbitals in both experimental [41,42] and theoretical [33,43] papers. However, many of the recent studies [44–48] revealed deviations from such a naive picture, and the understanding of the underling physics is challenging among the strong-field physics

063410-13


2pπ u+

2pπ g+ 2 × 10

-2

2 1 × 10

-2

0 -2 -4 -2

0

x (a.u.)

2

3 2 ×10-3

0

-3 -6

0

4

4 ×10

-3

6

z (a.u.)

z (a.u.)

4

-4


-6

-3

0

3

x (a.u.)

0

6

WFAT(0) WFAT(1) F F = 0.0056 0

β (deg)

F = 0.0015

60

120

F = 0.0224

180

F = 0.0061

FIG. 17. (Color online) (Top row) Electronic density |ψ0 (x,y = 0,z)|2 in the unperturbed 2pπu+ (Fc = 0.028) and 2pπg+ (Fc = 0.0076) states of H2 + . The other panels show polar plots of the WFAT(0) and WFAT(1) ionization rates for these states as functions of β at F = 0.2Fc (middle row) and F = 0.8Fc (bottom row). The rates in each plot are normalized to have the same maximum value.

community [49]. In this section, we address this issue in the example of H2 + . The interaction of this system with intense laser pulses was studied extensively by solving the time-dependent Schrödinger equation (TDSE) [33–35]. In particular, large discrepancies between the TDSE results and the MO-ADK formula were found for the 2pπu+ and 2pπg+ states [50]. The first (ungerade) of these states coincides with the 2pπ + state considered in Secs. IV A 2 and IV B 1; for brevity, we have omitted the subscript u in the above discussions. The field-free energy of the second (gerade) state in the present soft-core model is E0 = −0.224 302. The polarizabilities in this state are αx = −1147.9976 and αz = 438.1831, and the critical field is Fc = 0.0076. The dominant ionization channel for both states is ν = (0,0). According to the WFAT, the ionization rate in the limit F → 0 is given by the leading-order asymptotic formula (30). For the present system and states we have m = 0 and μz = 0. Therefore, the orientation dependence of the ionization rate is determined by the square of the coefficient g00 in the asymptotic tail of the unperturbed ionizing orbital ψ0 (r) defined by Eq. (25a). In other words, the orientation dependence of the ionization rate indeed provides an image of the electronic density in the asymptotic region. This is illustrated for the two states under consideration in the top and middle rows in Fig. 17. The top row shows the electronic density. The middle row shows polar plots of the WFAT ionization rates at a sufficiently weak field F = 0.2Fc . To focus on the shape of the dependence on β, the WFAT(0) and WFAT(1) results in each plot are normalized to have the same maximum values.

Only the interval 0◦ β 180◦ is discussed below, since the plots are symmetric with respect to the vertical z axis. It is seen that the first-order terms play only a minor role at this field and the shapes of the orientation dependence of the rates resemble that of the orbitals. In particular, for the 2pπu+ state, the rate has a nodal line at parallel geometries (β = 0◦ and 180◦ ), reflecting the nodal line of the orbital along the z axis, and a one-lobe shape with the maximum at β = 90◦ , as does the electronic density. Similarly, for the 2pπg+ state, the rate has two nodal lines reflecting the nodes of the orbital along the z axis and in the xy plane and a two-lobe shape with maxima at β = 45◦ and 135◦ reflecting the maxima of the electronic density. However, for stronger fields the situation becomes different. The bottom row in Fig. 17 shows the WFAT rates for F = 0.8Fc . The WFAT(0) results remain unchanged in this normalized representation, but the shapes of the WFAT(1) results considerably deviate from that of the orbitals. In the following, we focus on the two main features. For the 2pπu+ state, a local minimum of the rate is formed at β = 90◦ ; we have already discussed this feature in Sec. IV A 2; see Fig. 12. For the 2pπg+ state, the maxima of the lobes are shifted closer to β = 90◦ . We emphasize that these modifications of the shape of the orientation dependence of the ionization rate do not reflect any features of the shape of the ionizing orbital. They represent strong-field effects accounted for by the first-order correction terms in Eq. (34), and hence their explanation should be sought in the physical origin of these terms. Before we turn to such an analysis, let us compare the predictions of the WFAT with the TDSE results from Ref. [33]. The comparison for the two states is presented in the bottom panels in Figs. 18 and 19. The TDSE results show the ionization yields calculated for some particular laser pulse parameters (see captions to the figures). The WFAT results show the ionization rates calculated for a field equal to the peak value of the laser field. As is usually the case in comparisons with experimental results, here we are interested only in the shape of the dependence on the orientation angle β, so all the results are normalized to have the same unity value at their maxima. For both states, the WFAT(1) predictions are seen to be in very good agreement with the TDSE results. The two main features mentioned above, namely, the appearance of the local minimum at 90◦ and the maximum at 55◦ for the 2pπu+ state, as well as the shift of the maximum to 55◦ for the 2pπg+ state seen in the TDSE results, are well reproduced by the WFAT(1). The remaining difference can be attributed to the nonadiabatic effects and the difference between the present soft-core model (39) and the purely Coulomb potential with = 0 used in Ref. [33]. The WFAT(0), on the other hand, does not reproduce these features correctly. We note that the MO-ADK model produces results similar to the WFAT(0) and also fails to reproduce these features [50]. To explain the TDSE results in the framework of the WFAT(1), we analyze the different contributions to Eq. (34) in more detail. For the present nonpolar system A00 does not depend on β. Moreover, for the states under consideration the term with G01 is negligible except at very small β. Therefore, the modification of the shape of the orientation dependence of the rate as F grows discussed above results mainly from the term with B00 . This coefficient depends on β only through the

063410-14

Ionization yield (arb. units) WFAT coefficients (a.u.)



30 0

−καzz 2a00 B00

−30 −60 −90 1.0 0.8 0.6

TDSE W FAT(0) W FAT(1)

0.4 0.2 0.0 0

10

20

30

40

50

60

70

80

90

β (deg)

Ionizationyield (arb.units) WFAT coefficients (a.u.)

FIG. 18. (Color online) The orientation dependence of WFAT coefficients and ionization yield for H2 + (2pπu+ ). (Top) The terms −αzz and 2a00 on the right-hand side of Eq. (23b) and the total coefficient B00 . (Bottom) The ionization yield obtained in timedependent calculations [33] for the laser intensity I = 1013 W/cm2 and frequency ω = 0.043 (solid line with circles) and the WFAT ionization rates for the corresponding peak field F = 0.0169 (broken lines). The results in the bottom panel are normalized to have the same unity value at the maximum.

800

−καzz a00 B00

600 400 200 0 −200 −400 1.0 0.8 0.6 0.4

TDSE W FAT(0) W FAT(1)

0.2 0.0 0

10

20

30

40

50

60

70

80

90

β (deg) FIG. 19. (Color online) Same as in Fig. 18, but for H2 + (2pπg+ ). The time-dependent results [33] are for the laser intensity I = 1012 W/cm2 and frequency ω = 0.012, and the WFAT results are for the corresponding peak field F = 0.005 33.

terms −αzz and 2a00 in Eq. (23b) (a00 is real in the present case). The first of these terms accounts for the second-order Stark shift of the energy of the ionizing state, and the second one represents the first-order distortion of its wave function. The orientation dependence of these terms as well as that of the total coefficient B00 is shown in the top panels in Figs. 18 and 19. For both states, the main contribution to B00 comes from the Stark shift term −αzz . This term demonstrates strong dependence on β because of a large difference between the polarizabilities αx and αz ; see Eq. (16b). The distortion term 2a00 gives a smaller but also non-negligible contribution. Only by taking both these terms into account can the agreement between the WFAT(1) and TDSE results seen in the figures be achieved. We mention that the importance of the second-order Stark shift in strong-field ionization of molecules was recently recognized [51,52]. To account for this effect, it was suggested to simply incorporate the shift into the ionization potential within the MO-ADK model [53] or into the strong-field approximation [54]. However, as is clear from the above discussion, such an approach accounts neither for the distortion of the ionizing orbital, and hence it does not apply to situations where the term with a00 in Eq. (23b) plays a more prominent role, nor for the other first-order terms in Eq. (34). We emphasize that to correctly describe the deviation of the shape of the orientation dependence of the ionization rate in strong fields from that of the ionizing orbital, the consistent theory WFAT(1) including all the first-order terms in the asymptotic expansion (34) is required.

V. CONCLUSIONS

In this work, to validate the WFAT including the firstorder correction terms [19,25] for molecules, its results are compared with the results of accurate calculations of molecular SSs in an electric field [12–14]. Three states of different symmetry of a nonpolar molecule H2 + and one state of a polar molecule HeH2+ are considered. In all the cases, the first-order terms are shown to essentially improve the agreement between the WFAT and accurate results over a wide interval of fields up to the onset of over-the-barrier ionization. This not only demonstrates the good quantitative performance of the WFAT, but also confirms the high accuracy of our numerical method [12–14]. These results, together with similar calculations for atoms reported in Ref. [25], establish the WFAT including the first-order correction terms as an appealing alternative to laborious accurate calculations in the tunneling regime F < Fc . In addition to the quantitative improvement of the WFAT results, the first-order terms are shown to be helpful for a qualitative analysis of the field and orientation dependencies of the ionization rate and TMD, as well as the dependence of the TMD on the direction of the transverse momentum. These dependencies are determined by the interplay of the different terms in Eqs. (34)–(37). In particular, in the example of the 2pπ + state of H2 + at small β it is demonstrated that a rapid variation of the shape of the field and orientation dependence of the ionization rate can occur in situations where there is no a single dominant ionization channel. A similar behavior of the TMD in this case was discussed in Ref. [14]. In the example

063410-15



of the 2pσ state of HeH2+ it is shown that the first-order terms reflect the presence of an avoided crossing of molecular SSs localized on the different nuclei. One of the most sound implications of the present study for applications in strong-field physics concerns the orientation dependence of molecular ionization rates. We showed that the general belief that the shape of this dependence reproduces the shape of the unperturbed ionizing orbital, which lies in the foundation of various molecular imaging techniques and is supported by the leading-order approximation of the WFAT, holds only in the weak-field limit. As the field grows, the first-order correction terms can qualitatively modify the shape of the orientation dependence of the ionization rate. The physical origin of such a modification roots in the Stark shift and distortion of the ionizing orbital. This prediction of the WFAT is confirmed by comparison with time-dependent calculations of ionization yields in intense low-frequency laser fields for two states of H2 + [33]. Finally, we mention that there still remain discrepancies between theory and experiment [44–47] regarding the direction of the maximum ionization yield of CO, CO2 , and OCS in intense laser fields. The strong-field effects represented by the first-order correction terms of the WFAT should be considered together with effects of nuclear motion [21,22] and many-electron effects [20] as a possible source of such discrepancies.

we have

μz = − fn|M| r cos θ fn|M| .

For a σ state ψn0 , the other quantities needed are given by fn1 r sin θ fn0 2 αx = , (A2a) En 1 − En0 n fn0 r cos θ fn0 2 α =2 , En 0 − En0 n =n z

(A2c)

fn0 r cos θ fn0 1 fn0 (r ,θ ) √ . χz (r ) = En0 − En 0 2π n =n

(A2d)

+ For an even π state ψn1 , we find 0 fn r sin θ fn1 2 1 fn2 r sin θ fn1 2 αx = + , En 0 − En1 2 n En 2 − En1 n

(A3a) αz = 2

fn1 r cos θ fn1 2 n =n

APPENDIX: IMPLEMENTATION OF PERTURBATION THEORY FOR σ AND π ± STATES

To implement formulas of Sec. III A, we solve Eq. (1) for F = 0 by employing a single-center expansion in the molecular frame. A similar procedure was used in Ref. [13]. However, here we use the direct product of two discrete-variable representation (DVR) basis sets in r and θ constructed from the Laguerre and Gegenbauer polynomials [55], respectively. The use of the Laguerre-DVR basis in r eliminates the issue of convergence with respect to the radius of a finite spherical box used in Ref. [13]. Instead, there appears a “soft boundary” controlled by a scaling factor relating r to the argument of the basis functions. A proper choice of this factor accelerates convergence. This basis is also more suitable for extracting the asymptotic coefficients from the wave functions at η → ∞ needed to implement the WFAT. This procedure yields a complete set of the eigensolutions of Eq. (1) with energies En|M| and wave functions (13). Then it is straightforward to calculate the matrix elements in Eqs. (17) for the different symmetries of the initial state by replacing ψ0 with the corresponding function of the form (14). For all symmetries

(A2b)

fn1 r sin θ fn0 cos ϕ χx (r ) = fn1 (r ,θ ) √ , En0 − En 1 2π n

ACKNOWLEDGMENTS

This work was supported by the Grants-in-Aid for scientific research (A), (B), and (C) from the Japan Society for the Promotion of Science. O.I.T. acknowledges the support from the Russian Foundation for Basic Research (Grant No. 14-0292110) and the Ministry of Education and Science of Russia (State Assignment No. 3.679.2014/K).

(A1)

En 1 − En1

,

(A3b)

fn0 r sin θ fn1 1 χx (r ) = fn0 (r ,θ ) √ En1 − En 0 2 π n

fn2 r sin θ fn1 cos 2ϕ + fn2 (r ,θ ) √ , En1 − En 2 2 π n

fn1 r cos θ fn1 cos ϕ fn1 (r ,θ ) √ . χz (r ) = En1 − En 1 π n =n

− Similarly, for an odd π state ψn1 , we obtain 2 2 1 fn r sin θ fn1 αx = , 2 n En 2 − En1

fn1 r cos θ fn1 2 αz = 2 , En 1 − En1 n =n

(A3c)

(A3d)

(A4a)

(A4b)

fn2 r sin θ fn1 sin 2ϕ fn2 (r ,θ ) √ , χx (r ) = En1 − En 2 2 π n

(A4c)

fn1 r cos θ fn1 sin ϕ fn1 (r ,θ ) √ . χz (r ) = En1 − En 1 π n =n

(A4d)

The summations in these formulas run over all n , including the discretized continuum states with En |M| > 0.

063410-16



[1] F. Krausz and M. Ivanov, Attosecond physics, Rev. Mod. Phys. 81, 163 (2009). [2] O. I. Tolstikhin and T. Morishita, Adiabatic theory of ionization by intense laser pulses: Finite-range potentials, Phys. Rev. A 86, 043417 (2012). [3] J. Itatani, J. Levesque, D. Zeidler, H. Niikura, H. Pépin, J. C. Kieffer, P. B. Corkum, and D. M. Villeneuve, Tomographic imaging of molecular orbitals, Nature (London) 432, 867 (2004). [4] T. Morishita, A.-T. Le, Z. Chen, and C. D. Lin, Accurate retrieval of structural information from laser-induced photoelectron and high-order harmonic spectra by few-cycle laser pulses, Phys. Rev. Lett. 100, 013903 (2008). [5] S. Minemoto, T. Umegaki, Y. Oguchi, T. Morishita, A.-T. Le, S. Watanabe, and H. Sakai, Retrieving photorecombination cross sections of atoms from high-order harmonic spectra, Phys. Rev. A 78, 061402(R) (2008). [6] S. Haessler, J. Caillat, W. Boutu, C. Giovanetti-Teixeira, T. Ruchon, T. Auguste, Z. Diveki, P. Breger, A. Maquet, B. Carré, R. Ta¨ıeb, and P. Salières, Attosecond imaging of molecular electronic wavepackets, Nat. Phys. 6, 200 (2010). [7] W. Li, A. A. Jaroń-Becker, C. W. Hogle, V. Sharma, X. Zhou, A. Becker, H. C. Kapteyn, and M. M. Murnane, Visualizing electron rearrangement in space and time during the transition from a molecule to atoms, Proc. Natl. Acad. Sci. USA 107, 20219 (2010). [8] C. Wang, M. Okunishi, R. R. Lucchese, T. Morishita, O. I. Tolstikhin, L. B. Madsen, K. Shimada, D. Ding, and K. Ueda, Extraction of electron-ion differential scattering cross sections for C2 H4 by laser-induced rescattering photoelectron spectroscopy, J. Phys. B 45, 131001 (2012). [9] J. Mikosch, A. E. Boguslavskiy, I. Wilkinson, M. Spanner, S. Patchkovskii, and A. Stolow, Channel- and angle-resolved above threshold ionization in the molecular frame, Phys. Rev. Lett. 110, 023004 (2013). [10] H. Ohmura, N. Saito, and T. Morishita, Molecular tunneling ionization of the carbonyl sulfide molecule by double-frequency phase-controlled laser fields, Phys. Rev. A 89, 013405 (2014). [11] I. Petersen, J. Henkel, and M. Lein, Signatures of molecular orbital structure in lateral electron momentum distributions from strong-field ionization, Phys. Rev. Lett. 114, 103004 (2015). [12] P. A. Batishchev, O. I. Tolstikhin, and T. Morishita, Atomic Siegert states in an electric field: Transverse momentum distribution of the ionized electrons, Phys. Rev. A 82, 023416 (2010). [13] L. Hamonou, T. Morishita, and O. I. Tolstikhin, Molecular Siegert states in an electric field, Phys. Rev. A 86, 013412 (2012). [14] V. N. T. Pham, O. I. Tolstikhin, and T. Morishita, Molecular Siegert states in an electric field. II. Transverse momentum distribution of the ionized electrons, Phys. Rev. A 89, 033426 (2014). [15] L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Nonrelativistic Theory) (Pergamon Press, Oxford, 1977). [16] B. M. Smirnov and M. I. Chibisov, The breaking up of atomic particles by an electric field and by electron collisions, Zh. Eksp. Teor. Fiz. 49, 841 (1965) [,Sov. Phys. JETP 22, 585 (1966)]. [17] S. Yu. Slavyanov, Application of the method of standard comparison problems to perturbations of the Coulomb field. The discrete spectrum, Probl. Mat. Fiz. 4, 125 (1970) [English translation: Topics in Mathematical Physics (Consultants Bureau, New York, London, 1971), Vol. 4, P. 113].

[18] T. Yamabe, A. Tachibana, and H. J. Silverstone, Theory of the ionization of the hydrogen atom by an external electrostatic field, Phys. Rev. A 16, 877 (1977). [19] O. I. Tolstikhin, T. Morishita, and L. B. Madsen, Theory of tunneling ionization of molecules: Weak-field asymptotics including dipole effects, Phys. Rev. A 84, 053423 (2011). [20] O. I. Tolstikhin, L. B. Madsen, and T. Morishita, Weak-field asymptotic theory of tunneling ionization in many-electron atomic and molecular systems, Phys. Rev. A 89, 013421 (2014). [21] O. I. Tolstikhin, H. J. Wörner, and T. Morishita, Effect of nuclear motion on tunneling ionization rates of molecules, Phys. Rev. A 87, 041401(R) (2013). [22] O. I. Tolstikhin and L. B. Madsen, Retardation effects and the Born-Oppenheimer approximation: Theory of tunneling ionization of molecules revisited, Phys. Rev. Lett. 111, 153003 (2013). [23] R. J. Damburg and V. V. Kolosov, An asymptotic approach to the Stark effect for the hydrogen atom, J. Phys. B 11, 1921 (1978); ,A hydrogen atom in a uniform electric field. III, 12, 2637 (1979). [24] H. J. Silverstone, E. Harrell, and C. Grot, High-order perturbation theory of the imaginary part of the resonance eigenvalues of the Stark effect in hydrogen and of the anharmonic oscillator with negative anharmonicity, Phys. Rev. A 24, 1925 (1981). [25] V. H. Trinh, O. I. Tolstikhin, L. B. Madsen, and T. Morishita, First-order correction terms in the weak-field asymptotic theory of tunneling ionization, Phys. Rev. A 87, 043426 (2013). [26] A. R. Edmonds, Angular Momentum in Quantum Mechanics (Princeton University Press, Princeton, NJ, 1957). [27] M. Abramowitz and I. A. Stegun (eds.), Handbook of Mathematical Functions (Dover, New York, 1972). [28] L. B. Madsen, O. I. Tolstikhin, and T. Morishita, Application of the weak-field asymptotic theory to the analysis of tunneling ionization of linear molecules, Phys. Rev. A 85, 053404 (2012). [29] L. B. Madsen, F. Jensen, O. I. Tolstikhin, and T. Morishita, Structure factors for tunneling ionization rates of molecules, Phys. Rev. A 87, 013406 (2013). [30] L. B. Madsen, F. Jensen, O. I. Tolstikhin, and T. Morishita, Application of the weak-field asymptotic theory to tunneling ionization of H2 O, Phys. Rev. A 89, 033412 (2014). [31] A. M. Perelomov, V. S. Popov, and M. V. Terent’ev, Ionization of atoms in an alternating electric field, Zh. Eksp. Teor. Fiz. 50, 1393 (1966) [,Sov. Phys. JETP 23, 924 (1966)]; ,Ionization of atoms in an alternating electric field: II, Zh. Eksp. Teor. Fiz. 51, 309 (1966) [,Sov. Phys. JETP 24, 207 (1967)]. [32] The values of the polarizabilities given in Ref. [13] by mistake include the factor 1/2 from Eq. (15a) and hence are two times smaller than the present values. [33] G. Lagmago Kamta and A. D. Bandrauk, Imaging electron molecular orbitals via ionization by intense femtosecond pulses, Phys. Rev. A 74, 033415 (2006). [34] D. A. Telnov and Shih-I. Chu, Ab initio study of the orientation effects in multiphoton ionization and high-order harmonic generation from the ground and excited electronic states of H+ 2, Phys. Rev. A 76, 043412 (2007). [35] T. K. Kjeldsen, L. A. A. Nikolopoulos, and L. B. Madsen, Solving the m-mixing problem for the three-dimensional timedependent Schrödinger equation by rotations: Application to strong-field ionization of H+ 2 , Phys. Rev. A 75, 063427 (2007).

063410-17



[36] B. Zhang, J. M. Yuan, and Z. X. Zhao, Alignment-dependent ionization of H+ 2 : From multiphoton ionization to tunneling ionization, Phys. Rev. A 85, 033421 (2012). [37] The value of this energy given in Ref. [13] contains a misprint. [38] X. B. Bian and A. D. Bandrauk, Phase control of multichannel molecular high-order harmonic generation by the asymmetric diatomic molecule HeH2+ in two-color laser fields, Phys. Rev. A 83, 023414 (2011). [39] I. Ben-Itzhak, I. Gertner, O. Heber, and B. Rosner, Experimental evidence for the existence of the 2pσ bound state of HeH2+ and its decay mechanism, Phys. Rev. Lett. 71, 1347 (1993). [40] X. M. Tong, Z. X. Zhao, and C. D. Lin, Theory of molecular tunneling ionization, Phys. Rev. A 66, 033402 (2002). [41] A. S. Alnaser, S. Voss, X.-M. Tong, C. M. Maharjan, P. Ranitovic, B. Ulrich, T. Osipov, B. Shan, Z. Chang, and C. L. Cocke, Effects of molecular structure on ion disintegration patterns in ionization of O2 and N2 by short laser pulses, Phys. Rev. Lett. 93, 113003 (2004). [42] M. Meckel, D. Comtois, D. Zeidler, A. Staudte, D. Paviˇcić, H. C. Bandulet, H. Pépin, J. C. Kieffer, R. Dörner, D. M. Villeneuve, and P. B. Corkum, Laser-induced electron tunneling and diffraction, Science 320, 1478 (2008). [43] S.-F. Zhao, J. Xu, C. Jin, A.-T. Le, and C. D. Lin, Effect of orbital symmetry on the orientation dependence of strong field tunnelling ionization of nonlinear polyatomic molecules, J. Phys. B 44, 035601 (2011). [44] J. L. Hansen, L. Holmegaard, J. H. Nielsen, H. Stapelfeldt, D. Dimitrovski, and L. B. Madsen, Orientation-dependent ionization yields from strong-field ionization of fixed-in-space linear and asymmetric top molecules, J. Phys. B 45, 015101 (2012). [45] D. Paviˇcić, K. F. Lee, D. M. Rayner, P. B. Corkum, and D. M. Villeneuve, Direct measurement of the angular dependence of ionization for N2 , O2 , and CO2 in intense laser fields, Phys. Rev. Lett. 98, 243001 (2007). [46] I. Thomann, R. Lock, V. Sharma, E. Gagnon, S. T. Pratt, H. C. Kapteyn, M. M. Murnane, and W. Li, Direct measurement of the angular dependence of the single-photon ionization of aligned N2 and CO2 , J. Phys. Chem. A 112, 9382 (2008).

[47] J. Wu, L. P. H. Schmidt, M. Kunitski, M. Meckel, S. Voss, H. Sann, H. Kim, T. Jahnke, A. Czasch, and R. Dörner, Multiorbital tunneling ionization of the CO molecule, Phys. Rev. Lett. 108, 183001 (2012). [48] J. Yao, G. Li, X. Jia, X. Hao, B. Zeng, C. Jing, W. Chu, J. Ni, H. Zhang, H. Xie, C. Zhang, Z. Zhao, J. Chen, X. Liu, Y. Cheng, and Z. Xu, Alignment-dependent fluorescence emission induced by tunnel ionization of carbon dioxide from lower-lying orbitals, Phys. Rev. Lett. 111, 133001 (2013). [49] R. Murray, M. Spanner, S. Patchkovskii, and M. Y. Ivanov, Tunnel ionization of molecules and orbital imaging, Phys. Rev. Lett. 106, 173001 (2011). [50] S.-F. Zhao, C. Jin, A.-T. Le, T. F. Jiang, and C. D. Lin, Determination of structure parameters in strong-field tunneling ionization theory of molecules, Phys. Rev. A 81, 033423 (2010). [51] T. Brabec, M. Côté, P. Boulanger, and L. Ramunno, Theory of tunnel ionization in complex systems, Phys. Rev. Lett. 95, 073001 (2005). [52] D. Dimitrovski, M. Abu-samha, L. B. Madsen, F. Filsinger, G. Meijer, J. Küpper, L. Holmegaard, L. Kalhøj, J. H. Nielsen, and H. Stapelfeldt, Ionization of oriented carbonyl sulfide molecules by intense circularly polarized laser pulses, Phys. Rev. A 83, 023405 (2011). [53] H. Li, D. Ray, S. De, I. Znakovskaya, W. Cao, G. Laurent, Z. Wang, M. F. Kling, A. T. Le, and C. L. Cocke, Orientation dependence of the ionization of CO and NO in an intense femtosecond two-color laser field, Phys. Rev. A 84, 043429 (2011). [54] D. Dimitrovski, C. P. J. Martiny, and L. B. Madsen, Strong-field ionization of polar molecules: Stark-shift-corrected strong-field approximation, Phys. Rev. A 82, 053404 (2010). [55] O. I. Tolstikhin and C. Namba, CTBC—A Program to Solve the Collinear Three-Body Coulomb Problem: Bound States and Scattering Below the Three-Body Disintegration Threshold, Research Report NIFS-779, National Institute for Fusion Science, Toki, Japan, 2003, available at http://www.nifs.ac.jp/report/nifs779.html

063410-18