Accurate Structure Parameters for Tunneling Ionization Rates of Gas-Phase Linear

Song-Feng Zhao(赵松峰), Jian-Ke Li(栗建科),Guo-Li Wang(王国利),Peng-Cheng Li(李鹏程),and Xiao-Xin Zhou(周效信)

College of Physics and Electronic Engineering,Northwest Normal University,Key Laboratory of Atomic and Molecular Physics and Functional Materials of Gansu Province,Lanzhou 730070,China

1 Introduction

With the rapid advancements in field-free alignment and high-resolution recoil-ion momentum spectroscopy(RIMS)techniques,[1−2]tunneling ionization of molecules in a strong infrared laser field has attracted much attention in recent years,[3]especially for the alignment dependent ionization probabilityP(θ)of a molecule fixed in space,[4−13]whereθis the angle between the molecular axis and the polarization direction of the laser’s electric field.This is because tunneling ionization of molecules is the first step in many interesting strong- field phenomena like high-order harmonic generation(HHG),high-energy above-threshold ionization(HATI),nonsequential double ionization(NSDI)and molecular imaging using its own recolliding electrons.[14−17]

Theoretically,P(θ)has been obtained by solving the time-dependent Schrödinger equation(TDSE)based on the single-active electron approximation,[18−24]performing time-dependent Hartree–Fock(TDHF)calculations[25]or using time-dependent density functionaltheory(TDDFT).[26−29]However,these ab-initio calculations are rather time-consuming and still quite challenging even for a single laser intensity. Fortunately some simpler theoretical models are developed and desirable for interpreting experimental results,such as the molecular Ammosov–Delone–Krainov(MO-ADK),[30−36]the molecular strong- field approximation (MO-SFA),[37−38]the molecular Perelomov–Popov–Terent’ev(MO-PPT),[39−41]the weak field asymptotic theory(WFAT)[42−44]and so on.

Recall that the Ammosov–Delone–Krainov(ADK)model was initially developed to study the tunneling ionization rate of atoms using quasistatic approximation.[45]Tonget al.[30]have generalized successfully the ADK model to diatomic molecules by considering the symmetry property and the asymptotic behavior of the molecular electronic wave functions.The MO-ADK has also been extended to calculate the tunneling ionization rate of nonlinear polyatomic molecules.[32,36]In the MO-ADK model,[30,32,36]the static tunneling ionization rate of molecules is given analytically and the rate depends on structure parameters,ionization potential(IP),alignment(or orientation)angles and the instantaneous electric field of the laser.Structure parameters of molecules can be extracted numerically from wave functions in the asymptotic region. We mention that the correct asymptotic tail of molecular wave functions is required to obtain the accurate structure parameters. In Tonget al.,[30]wave functions of diatomic molecules were calculated originally with the multiple-scattering theory.[46]However,molecular wave functions are more accessible from some quantum chemistry packages such as GAUSSIAN,[47]GAMESS[48]or MOLPRO.[49]Since molecular wave functions are expanded in terms of Gaussian basis functions in these packages,the extracted structure parameters from these molecular wave functions are not accurate enough because these basis functions do not have the correct behavior in the asymptotic region.[34]In our previous works,[34−35,50−52]we proposed an efficient method to fixthe asymptotic tail of molecular wave functions by solving the time-independent Schrödinger equation with B-spline functions,where molecular potentials can be constructed numerically based on the density-functional theory(DFT).Using this method,Madsenet al.determined structure parameters of the highest occupied molecular orbital(HOMO)for CO2with wave functions calculated by propagating the TDSE in imaginary time.[53]These days,the MO-ADK model has been one of the most favorable molecular tunneling ionization models in strong field physics and attosecond science.However,structure parameters are only available for a few molecules at the present time.Therefore,it is quite essential to determine systematically these accurate structure parameters for the MO-ADK model.We emphasize that the main goal of this paper is to extract and tabulate accurate molecular structure parameters for the HOMO of 123 gas-phase linear molecules by using our previous method.[34−35,50−52]

2 Theoretical Method

The theory part is separated into four subsections.Firstly,we will review briefly the MO-ADK model.Then we present how to construct numerically one-electron potentials of linear molecules using the LBαmodel and obtain molecular wave functions with the correct asymptotic tail by solving the time-independent Schrödinger equation with B-spline functions.Finally we describe how to extract structure parameters from molecular wave functions in the asymptotic region.

2.1 MO-ADK Model

Based on the MO-ADK model,[30]the tunneling ionization rate in a static field is given analytically by

whereFandZcare the field strength and the charge of the molecular residue(Zc=1 for neutral molecules),respectively.is the ionization potential for the given valence orbital.In Eq.(1),for linear molecules,B(m′)can be expressed as

whereClmare structure parameters of the ionizing orbital andRis the Euler angles of the molecular frame with respect to the laboratory fixed frame.Dml′,m(R)is the rotation matrix and

2.2 Construction of One-Electron Potentials for Linear Molecules

The one-electron potentials are required for solving the time-independent Schrödinger equation of linear molecules.In this paper,we use the LBαmodel[50,54−55]to calculate numerically the one-electron potentials consisting of electrostatic and exchange-correlation terms.

Based on the single-center expansion,the one-electron potentials of linear molecules can be expressed as

wherevl(r)is the partial radial component of the potential andPl(cosθe)is the Legendre polynomial.θeis the angular coordinate of the active electron in the molecular frame.The partial radial potential can be rewritten as

where the first two terms is the electrostatic potential and the last term describes the exchange-correlation interaction.

The partial electron-nucleus interactionvlnuc(r)can be written as

whereisums over theNaatoms in the molecule.Without loss of generality,we assume that linear molecules are aligned along thez-axis,thenvli(r)can be expressed as

The partial Hartree potentialis given by

withr＜=min(r,r′),r＞=max(r,r′).Hereal(r′)is

whereρ(r′,θ′)is the total electron density in the molecule and

Hereiruns over all theNeelectrons in the molecule.The wave functions of each occupied molecular orbital can be obtained from GAUSSIAN.[47]

In LBαmodel,the partial exchange-correlation potential is written as

where

The LDA correlation potentialVcL,σDA(r,θe)is calculated by using the Perdew–Wang representation for the correlation functionals[56]

εc(rs,ζ)is the correlation energy

where

In Eq.(15),εc(rs,0),εc(rs,1)and−αc(rs)can be given analytically

ParametersA,α1,β1,β2,β3,β4andphave been tabulated in Table 1 of Ref.[56].rs=[3/4π(ρ↑+ρ↓)]1/3andζ=(ρ↑−ρ↓)/(ρ↑+ρ↓)are density parameter and relative spin polarization,respectively.Here,ρ↑andρ↓are total density for spin-up and spin-down electrons,respectively.This potential is spin dependent.Note that sgnσis+1 forσ=↑and−1 forσ=↓.

2.3 Calculation of Molecular Wave Functions by Solving the Time-independent Schr¨odinger Equation

Once molecular potentials are created using the method described in Subsec.2.2,wave functions with the correct asymptotic behaviour can be obtained by solving the following time-independent Schrödinger equation of linear molecules

Because of the cylindrical symmetry,the wave functioncan be written as

withξ=cosθeandχeis the angular coordinate of the active electron in the molecular frame.The wave function

ψ(r,ξ)can be expanded by B-spline functions as

Here,Bi(r)andBj(ξ)are radial and angular B-splines,respectively.We take the corresponding number of radial and angular B-splinesNr=80 andNξ=30,respectively.By substituting Eqs.(4),(19)and(20)into Eq.(18)and then projecting onto theBi′(r)(1−ξ2)|m|/2Bj′(ξ)basis,we obtain the following matrix equation

with

whereEandCare energy matrix and coefficient matrix,respectively.The eigenfunctions and eigenvalues can be obtained by diagonalizing Eq.(21).

2.4 Determination of Molecular Structure Parameters of the HOMO for Gas-Phase Linear Molecules

Based on the single-center expansion,wave functions of the ionizing orbital for a linear molecule can also be expanded as

where Ylm(θe,χe)is the spherical harmonic functions.The radial functions can be calculated by

Then the accurate structure parametersClmcan be extracted by fitting these calculated radial functions to the following form in the asymptotic region

3 Results

Taking CO as an example,we show partial wave expansions,andυl(r))of the potentialV(r,θe)(see Eqs.(4)and(5))in Fig.1.

Wave functions and the corresponding ionization potentials of the ionized orbital of linear molecules can be obtained by solving Eq.(18).In Table 1,we compare ionization potentials from the present calculations using the LBαmodel with experimental vertical ionization potentials or the calculated ionization potentials with GAUSSIAN code when the experimental ones are unavailable.In the LBαmodel,αandβare two empirical parameters.In the present calculations,we takeβ=0.01 andαis optimized to obtain the accurate ionization potential by comparing with the corresponding experimental vertical ionization potential or the calculated results with GAUSSIAN code.These optimizedαvalues and molecular internuclear distances are also listed in Table 1.

Fig.1 (a)Partial electron-nucleus potential;(b)partial electron-electron repulsion potential;(c)partial exchange-correlation potential;(d)partial total potential.We use lmax=40[see Eq.(4)]for CO,but we only show here the l=0,1,2,3,4 terms for clarity.

Table 1 Internuclear distances,ionization potentials calculated with the present LBα model and experimental vertical ionization potentials for 123 selected gas-phase linear molecules.The calculated ionization energies with GAUSSIAN code for several molecules are also given when the experimental results are unavailable and α parameter used in LBα model are also listed.

Table 1 (continued)

Table 1 (continued)

Table 1 (continued)

We followed the fitting method proposed in Ref.[38]to extract molecular structure parameters in the asymp-totic region.In Fig.2,the accurate structure parameters of the ground state ofH+2at the equilibrium distance are determined by fitting the calculated radial wave function to the analytical form in Eq.(26).The extracted structure parametersC00,C20,andC40are 4.52,0.62,and 0.03,respectively.Using this fitting method,we extracted systematically accurate structure parameters for the HOMO of 123 gas-phase linear molecules and tabulated in Table 2 for future applications.We emphasize that these structure parameters are very crucial for calculating tunnelling ionization rates of gas-phase linear molecules using the MO-ADK model.[30]We mention that these structure parameters can also be used in the MO-PPT model,which works very well in a wide range covering from the multiphoton to the tunnelling ionization regimes.In the appendix,thezcoordinates of all atoms in a polar molecule are given in the Table 3.

Fig.2 Comparison of the asymptotic behavior of the present calculated radial wave function(solid lines)to the analytical form in Eq.(26)(dashed lines)for H+2.

We should address that there are some discrepancies between the present fitted structure parameters and those tabulated in previous references.[3,30,33−35,53,58,60]The discrepancies are originated from the different qualities of asymptotic tail of molecular wavefunctions.In Refs.[3]and[60],the molecular wavefunctions are obtained from the GAUSSIAN programs.[47]The asymptotic behavior of these molecular wavefunctions is not correct because the GAUSSIAN basis functions do not have the exponential decay form in the largerr.In Ref.[33],the program X2DHF based on the Hartree–Fock method was used to determine diatomic molecular wavefunctions with the correct asymptotic behavior.However,this program can not be used for polyatomic molecules.In Ref.[30],the multiple scattering theory was used to calculate molecular wavefunctions.However,we found that structure parameters obtained from the multiple scattering theory are also not accurate enough for simulating the experimental ionization probabilities of CO2.[34]In Refs.[34]–[35],[53],and[58],molecular wavefunctions are obtained by solving the time-independent Schrödinger equation(TIDSE)where the electron correlation has been neglected in the molecular potentials(to be called “exchange-only” model in the following).In our previous work,[50]the correlation interaction was included in the modified Leeuwen–Baerends(called LBα)model proposed by Schipperet al.[54]We con firmed that the binding energies with the LBαmodel are much more accurate than those using the“exchangeonly”model.In this paper,we calculate molecular wavefunctions with the correct asymptotic tail by solving the TIDSE with the B-spline functions and molecular potentials constructed numerically with the LBαmodel.In Table 1,one can see that our binding energies agree very well with the experimental ones in deed.We thus expect our structure parameters are more accurate than those obtained from other methods mentioned above.

Table 2 Fitted Clmcoefficients for the HOMO of 123 selected gas-phase linear molecules.Coefficients from earlier references are also listed for direct comparison.

Table 2 (continued)

Table 2 (continued)

Table 2 (continued)

Another issue that should be addressed is the validity of these fitted structure parameters for linear molecules.In our previous works,[34−35,40]we introduced an efficient method to improve the asymptotic behavior of molecular wave functions by solving the TIDSE with B-spline functions and extracted the accurate structure parameters for some linear molecules.By using our fitted structure parameters,the MO-ADK results are improved significantly by comparing with the TDSE and the TDDFT simulations for H+2,H2,N2,O2,F2,and C2H2(see Fig.2 in Ref.[35]and Fig.4 in Ref.[40]).These structure parameters are further examined by comparing the MO-ADK and the MO-PPT results with the experimental data(see Fig.3 in Ref.[34],Fig.6 in Ref.[35],and Fig.5 in Ref.[40]).We found that the MO-ADK and the MO-PPT models agree well with the experimental results if our fitted structure parameters are used.We also expect that these structure parameters tabulated in this paper are checked when the more accurate simulations and the experimental data become available in the near future.

4 Conclusions

In this paper,we determined systematically accurate structure parameters used in the MO-ADK and MO-PPT models for the HOMO of 123 gas-phase linear molecules.Wave functions with the correct asymptotic tail are obtained by solving the time-independent Schrödinger equation of linear molecules using B-spline functions,where molecular potentials are created numerically with the LBαmodel.We emphasize that these accurate structure parameters are very essential for calculating strong- field tunneling ionization rates of linear molecules using the MOADK and MO-PPT models.Furthermore,It is straightforward to generalize the present method for extracting accurate structure parameters of nonlinear polyatomic molecules.

Appendix

Table 3 The z coordinates(in Angstroms)of atoms for several selected gas-phase linear molecules at equilibrium distance calculated from the GAUSSIAN packages.

Table 3 (continued)

Table 3 (continued)

[1]H.Stapelfeldt and T.Seideman,Rev.Modern Phys.75(2003)543.

[2]J.Ullrich,R.Moshammer,A.Dorn,R.Dörner,L.Ph.H.Schmidt,and H.Schmidt-Böcking,Rep.Prog.Phys.66(2003)1463.

[3]C.D.Lin,X.M.Tong,and Z.X.Zhao,J.Mod.Opt.53(2006)21.

[4]H.Akagi,T.Otobe,A.Staudte,et al.,Science 325(2009)1364.

[5]I.V.Litvinyuk,K.F.Lee,P.W.Dooley,et al.,Phys.Rev.Lett.90(2003)233003.

[6]D.Paviˇci´c,K.F.Lee,D.M.Rayner,et al.,98(2007)243001.

[7]A.Staudte,S.Patchkovskii,D.Paviˇci´c,et al.,Phys.Rev.Lett.102(2009)033004.

[8]J.Wu,M.Meckel,S.Voss,et al.,Phys.Rev.Lett.108(2012)043002.

[9]A.S.Alnaser,S.Voss,X.M.Tong,et al.,Phys.Rev.Lett.93(2004)113003.

[10]I.Thomann,R.Lock,V.Sharma,et al.,J.Phys.Chem.A 112(2008)9382.

[11]P.von den Hoff,I.Znakovskaya,S.Zherebtsov,et al.,Appl.Phys.B 98(2010)659.

[12]X.Liu,C.Wu,Z.Wu,Y.Liu,Y.Deng,and Q.Gong,Phys.Rev.A 83(2011)035403.

[13]H.Liu,S.F.Zhao,M.Li,et al.,Phys.Rev.A 88(2013)061401(R).

[14]C.D.Lin,A.T.Le,Z.Chen,T.Morishita,and R.Lucchese,J.Phys.B 43(2010)122001.

[15]M.Lein,J.Phys.B 40(2007)R135.

[16]S.Haessler,J.Caillat,and P.Sali`eres,J.Phys.B 44(2011)203001.

[17]W.Becker,X.Liu,P.J.Ho,and J.H.Eberly,Rev.Modern Phys.84(2012)1011.

[18]S.Petretti,Y.V.Vanne,A.Saenz,et al.,Phys.Rev.Lett.104(2010)223001.

[19]M.Spanner and S.Patchkovskii,Phys.Rev.A 80(2009)063411.

[20]M.Abu-samha and L.B.Madsen,Phys.Rev.A 80(2009)023401.

[21]G.Lagmago Kamta and A.D.Bandrauk,Phys.Rev.A 74(2006)033415.

[22]Y.J.Jin,X.M.Tong,and N.Toshima,Phys.Rev.A 83(2011)063409.

[23]B.Zhang,J.Yuan,and Z.X.Zhao,Phys.Rev.A 85(2012)033421.

[24]S.L.Hu,J.Chen,X.L.Hao,and W.D.Li,Phys.Rev.A 93(2016)023424.

[25]B.Zhang,J.Yuan,and Z.X.Zhao,Phys.Rev.Lett.111(2013)163001.

[26]D.A.Telnov and Shih-I.Chu,Phys.Rev.A 79(2009)041401(R).

[27]X.Chu,Phys.Rev.A 82(2010)023407.

[28]X.Chu and M.McIntyre,Phys.Rev.A 83(2011)013409.

[29]T.Otobe and K.Yabana,Phys.Rev.A 75(2007)062507.

[30]X.M.Tong,Z.X.Zhao,and C.D.Lin,Phys.Rev.A 66(2002)033402.

[31]Z.X.Zhao,X.M.Tong,and C.D.Lin,Phys.Rev.A 67(2003)043404.

[32]T.K.Kjeldsen,C.Z.Bisgaard,L.B.Madsen,and H.Stapelfeldt,Phys.Rev.A 71(2005)013418.

[33]T.K.Kjeldsen and L.B.Madsen,Phys.Rev.A 71(2005)023411.

[34]S.F.Zhao,C.Jin,A.T.Le,T.F.Jiang,and C.D.Lin,Phys.Rev.A 80(2009)051402(R).

[35]S.F.Zhao,C.Jin,A.T.Le,T.F.Jiang,and C.D.Lin,Phys.Rev.A 81(2010)033423.

[36]S.F.Zhao,J.Xu,C.Jin,A.T.Le,and C.D.Lin,J.Phys.B 44(2011)035601.

[37]J.Muth-Böhm,A.Becker,and F.H.M.Faisal,Phys.Rev.Lett.85(2000)2280.

[38]T.K.Kjeldsen and L.B.Madsen,J.Phys.B 37(2004)2033.

[39]E.P.Benis,J.F.Xia,X.M.Tong,et al.,Phys.Rev.A 70(2004)025401.

[40]S.F.Zhao,L.Liu,and X.X.Zhou,Opt.Commun.313(2014)74.

[41]S.F.Zhao,A.T.Le,C.Jin,X.Wang,and C.D.Lin,Phys.Rev.A 93(2016)023413.

[42]L.B.Madsen,O.I.Tolstikhin,and T.Morishita,Phys.Rev.A 85(2012)053404.

[43]L.B.Madsen,F.Jensen,O.I.Tolstikhin,and T.Morishita,Phys.Rev.A 87(2013)013406.

[44]L.B.Madsen,F.Jensen,O.I.Tolstikhin,and T.Morishita,Phys.Rev.A 89(2014)033412.

[45]M.V.Ammosov,N.B.Delone,and V.P.Krainov,Zh.Eksp.Theor.Fiz.91(1986)2008[Sov.Phys.JETP 64(1986)1191].

[46]D.Dill and J.L.Dehmer,J.Chem.Phys.61(1974)692.

[47]M.J.Frisch,G.W.Trucks,H.B.Schlegel,et al.,GAUSSIAN 03,Revision B.04(Gaussian Inc.Pittsburgh,PA(2003).

[48]M.W.Schmidt,K.K.Baldridge,J.A.Boatz,et al.,J.Comput.Chem.14(1993)1347.

[49]H.J.Werner and P.J.Knowles,MOLPRO,Version 2002.6,A Package of Ab Initio Programs,Birmingham,UK(2003).

[50]S.F.Zhao,C.Jin,A.T.Le,and C.D.Lin,Phys.Rev.A 82(2010)035402.

[51]J.P.Wang,S.F.Zhao,C.R.Zhang,W.Li,and X.X.Zhou,Mol.Phys.112(2014)1102.

[52]S.F.Zhao,F.Huang,G.L.Wang,and X.X.Zhou,Commun.Theor.Phys.65(2016)366.

[53]M.Abu-samha and L.B.Madsen,Phys.Rev.A 81(2010)033416.

[54]P.R.T.Schipper,O.V.Gritsenko,S.J.A.van Gisbergen,and E.J.Baerends,J.Chem.Phys.112(2000)1344.

[55]Shih-I.Chu,J.Chem.Phys.123(2005)062207.

[56]J.P.Perdew and Y.Wang,Phys.Rev.B 45(1992)13244.

[57]R.Saito,O.I.Tolstilchin,L.B.Madsen,and T.Morishita,Atomic Data and Nuclear Data Tables 103-104(2015)4.

[58]M.Abu-samha and L.B.Madsen,Phys.Rev.A 82(2010)043413.

[59]J.Heslar,D.Telnov,and Shih-I.Chu,Phys.Rev.A 83(2011)043414.

[60]A.T.Le,X.M.Tong,and C.D.Lin,J.Mod.Opt.54(2007)967.