張斌 趙健 趙增秀

1)(西北核技術(shù)研究所,西安 710024)

2)(國防科技大學(xué)物理系,長沙 410073)

1 引 言

電子關(guān)聯(lián)在物理和化學(xué)的很多多體相互作用過程中發(fā)揮了重要作用[1].當(dāng)原子分子處于強(qiáng)激光場中時,電子關(guān)聯(lián)效應(yīng)最典型的一個例子是:即使強(qiáng)場不足以使得體系以順序電離的方式電離兩個電子的情況下,仍然有發(fā)生雙電離的可能.這被稱為非順序雙電離,而電子動態(tài)關(guān)聯(lián)作用被認(rèn)為是該過程的物理根源[2].當(dāng)一個電離的電子在激光場作用下與離子發(fā)生再碰撞,它有可能將另外一個電子通過“knock out”的方式電離,或者將其激發(fā)到高里德伯態(tài),然后再被強(qiáng)場電離[3,4].這兩個過程都依賴于電子關(guān)聯(lián)效應(yīng)[5].電子關(guān)聯(lián)效應(yīng)會導(dǎo)致雙電子動量譜中出現(xiàn)“V”形的手指狀結(jié)構(gòu)[6].研究關(guān)聯(lián)效應(yīng)對于非微擾強(qiáng)場物理過程的影響,需要非微擾的理論處理方法[7?9].由于計算量的問題,直接數(shù)值求解含時薛定諤方程（TDSE）的方法只能用于不超過兩個電子的多體系統(tǒng)[10,11].由于具有額外的自由度,分子在強(qiáng)場中的響應(yīng)要比原子更加復(fù)雜,這給理論上的處理帶來很大挑戰(zhàn).基于此,最簡單的分子體系——雙原子分子成為了研究分子體系基本物理過程的理想對象.對于雙原子分子,其額外的自由度為分子的核間距以及分子軸與激光電場偏振方向的夾角,其中后者可以通過脈沖激光來實現(xiàn)取向控制[12].然而,即使對于相對簡單的雙原子分子,如果分子中包含多個電子,那么理論上的處理仍然是非常困難的.需要發(fā)展能夠研究多體關(guān)聯(lián)效應(yīng)的理論方法.

作為含時Hartree-Fock(TDHF)理論[13]的推廣,人們發(fā)展了多組態(tài)的TDHF理論(MCTDHF).MCTDHF方法將多電子波函數(shù)用多個Slater行列式展開,這樣就能夠?qū)φ鎸嵉亩嚯娮硬ê瘮?shù)提供足夠好的近似.隨著組態(tài)數(shù)的增加,結(jié)果將不斷接近精確結(jié)果.因此,MCTDHF能夠比較系統(tǒng)地處理多電子交換關(guān)聯(lián)效應(yīng),而不需要像含時密度泛函方法那樣采用經(jīng)驗形式的交換關(guān)聯(lián)勢[14].不同于含時組態(tài)相互作用方法(TDCI),MCTDHF中軌道波函數(shù)也是隨時間動態(tài)演化的.TDCI方法的收斂一般需要非常多的組態(tài)數(shù)[15],而MCTDHF則對于組態(tài)數(shù)量的要求大大減小.由于具有上述優(yōu)點,MCTDHF方法在強(qiáng)場領(lǐng)域也得到越來越多的關(guān)注.Zanghellini等[16]首先將MCTDHF方程引入強(qiáng)場領(lǐng)域.隨后,人們利用一維的MCTDHF理論研究了強(qiáng)激光場作用下電子關(guān)聯(lián)效應(yīng)的影響[17?19]. 直到2011年,Hochstuhl和Bonitz[20]以及Haxton等[21]才發(fā)展出了三維的MCTDHF方法,并用于討論簡單原子分子體系(He和H2)在短遠(yuǎn)紫外(XUV)激光脈沖下的電離行為.在最新的工作中,三維的MCTDHF方法已經(jīng)被用于O2的自電離態(tài)研究[22].MCTDHF包含一組耦合的組態(tài)系數(shù)方程和軌道方程,其計算量遠(yuǎn)大于TDHF和TDDFT.由于軌道方程是非線性的,因此大都采用了較為一般的數(shù)值方法,如Runge-Kutta(RK)方法,其穩(wěn)定性較差,因此還有待于發(fā)展具有更高計算效率和穩(wěn)定性的數(shù)值方法.

為了研究電子關(guān)聯(lián)對于強(qiáng)場物理過程的影響,針對雙中心分子體系,本文發(fā)展了三維的多組態(tài)含時Hartree-Fock方法,能夠處理多電子雙原子分子在強(qiáng)激光場作用下的電子動力學(xué)行為.該方法采用橢球坐標(biāo)系以及有限元(FE)-離散變量基函數(shù)(DVR)方法[23],為在雙中心體系中應(yīng)用MCTDHF方法奠定了基礎(chǔ).進(jìn)一步,利用該方法研究了強(qiáng)XUV脈沖作用下H2分子電離過程中的電子關(guān)聯(lián)效應(yīng).

2 理論方法

2.1 MCTDHF理論

在波恩-奧本海默近似下,描述分子中電子動力學(xué)的含時薛定諤方程為[使用原子單位(a.u.)]

其中ψ(r,t)為含時多電子波函數(shù).多電子哈密頓量由單電子哈密頓量和雙電子相互作用部分組成:

其中i和j為電子的指標(biāo).在偶極近似下,用F(t)表示激光電場,采用長度規(guī)范形式下的激光與電子相互作用勢[24],則單電子哈密頓量表示為

上式右邊第一項為電子動能算符,第二項為電子與原子核的庫侖勢(其中Zk和Rk分別為第k個原子核的電荷和空間坐標(biāo),ri為電子的空間坐標(biāo)).(2)式中的雙電子算符為

在MCTDHF理論中,多電子波函數(shù)ψ(r,t)為多個Slater行列式的線性組合[16].以Nα(Nβ)表示自旋為α(β)的電子數(shù),M表示空間軌道的數(shù)量,|?p(t)?表示空間軌道.本文采用自旋限制的理論形式,即每個軌道的自旋部分和空間部分是獨立的,且?. 在二次量子化表象下,Ψ(r,t)表示為

其中Cn(t)為組態(tài)n的展開系數(shù).對于一個給定的組態(tài)空間,(5)式中的求和遍歷所有可能的個行列式.不失一般性,可假設(shè)Nα>Nβ,則自旋投影為Sz=(Nα?Nβ)/2.

在二次量子化表象下,(2)式中的可由單激發(fā)算子和雙激發(fā)算子表示為[20]

其中,分別為

這里σ和τ為自旋指標(biāo),pqrs為空間軌道指標(biāo),和分別為產(chǎn)生和湮沒算子.(6)式中,hpq和分別為[(3)式]和[(4)式]在空間軌道{?p}下的矩陣元,

通過拉格朗日變分原理,導(dǎo)出的電子運動方程可以寫為兩個耦合的方程組[20],其中描述組態(tài)系數(shù)演化的方程為

描述空間軌道演化的方程為

其中投影算子為單電子密度矩陣

dpqrs為雙電子密度矩陣,

在坐標(biāo)表象下,平均場算符為

(12)式在單組態(tài)情況下退化為TDHF形式.(11)和(12)式是相互耦合的,必須同時進(jìn)行求解.求解方法在后文進(jìn)行敘述.

2.2 數(shù)值方法

2.2.1 橢球坐標(biāo)系

如果雙原子分子的分子軸沿z軸,原子核坐標(biāo)分別在?a和a(核間距為R=2a),那么直角坐標(biāo)(x,y,z)與橢球坐標(biāo)(ξ,η,φ)的對應(yīng)關(guān)系為[25]:

其中16ξ6∞,?16η61,06φ62π.在橢球坐標(biāo)系下,電子與原子核(電荷分別為Z1和Z2)之間的庫侖勢為,

可以看出,庫侖勢的奇異性在體積元dV=a3(ξ2?η2)dξdηdφ的積分下被消除.這表明在橢球坐標(biāo)系中庫侖勢能夠被精確處理,對于雙中心體系具有天然的優(yōu)勢.

2.2.2 FE-DVR

空間軌道采用FE與DVR相結(jié)合的方法進(jìn)行離散.?(ξ,η,φ;t)在FE-DVR/DVR基函數(shù)下展開為

其中,m為磁量子數(shù),nm為展開式中最大的磁量子數(shù).對于ξ坐標(biāo),我們用了FE-DVR基函數(shù).以ξmax表示ξ坐標(biāo)的截斷值(ξmax要足夠大), 首先將區(qū)間[71,ξmax]分成ne個子區(qū)間ξ1=1<ξ2<ξ3<···<ξne?1<ξne=ξmax. 每一個區(qū)間l(定義在區(qū)間[ξl,ξl+1])通過廣義的高斯積分方法進(jìn)行離散[26].用nξ表示所有區(qū)間上基函數(shù){fi(ξ)}的總數(shù)量.另一方面,對于坐標(biāo)η(從?1到1),利用一組DVR基函數(shù)即可,采用nη點的高斯-勒讓德方法進(jìn)行離散[27].關(guān)于FE-DVR/DVR基函數(shù)更具體的介紹見文獻(xiàn)[23].

f和g都滿足正交歸一條件,即?fi|fi′?=δii′,?gj|gj′?=δjj′.三維基函數(shù)總個數(shù)為nb=nξnη(2nm+1). 任何局域勢在FE-DVR/DVR基下都是對角化的,可以用一個矢量表示.對于一般的基函數(shù)(如高斯或Slater函數(shù)[1]),雙電子積分的計算量正比于,這是MCTDHF計算中最耗時的部分.在本文的FE-DVR方法中,計算量降低到量級,為開展 MCTDHF 計算奠定了基礎(chǔ).

2.2.3 時間演化方法

下面討論時間傳播的方法.對于(11)式中的組態(tài)系數(shù)方程,我們采用了Short Iterative Lanczos(SIL)方法[28].一般地,(12)式中的軌道方程是非線性的,需要借助更一般的微分方程求解方法.本文采用了兩種方法.第一種是RK方法.這是求解微分方程最一般的方法.為了減少傳播誤差從而提高數(shù)值穩(wěn)定性,本文采用了對于時間步長具有8階精度的RK方法[27].第二種是Bulirsch-Stoer(BS)外推方法,其基本思想是:將一個大的時間間隔逐漸細(xì)分為間隔更小的子序列,每個區(qū)間上的積分采用修正中點法,最終結(jié)果是對逐步細(xì)分的子序列的積分結(jié)果進(jìn)行外推,其中外推采用有理函數(shù)或多項式[27].本文對這兩種方法都采用了自適應(yīng)精度的處理方式:即根據(jù)設(shè)定的誤差上限動態(tài)調(diào)整時間傳播步長.在兩種特殊情況下,用自適應(yīng)步長的RK方法會有更好的效果:一是方程右邊的函數(shù)值是通過查表或插值來定義的;二是微分方程的積分區(qū)間中有奇異點.除了上面兩個缺點外,BS方法是目前以最小的計算量達(dá)到高精度求解微分方程的一個最好的方法[27].兩種方法的結(jié)果可以互相比較,以驗證傳播方法的正確性.

2.3 物理量的計算

2.3.1 分子總能量

分子總能量E(t)計算為

其中右邊第一項表示電子態(tài)的能量,第二項表示原子核之間的庫侖排斥能.將(6)式代入(19)式中右邊第一項,利用密度矩陣Dpq和dpqrs的定義得到

2.3.2 單電離概率

為了計算單電離概率,首先定義單電子密度.在坐標(biāo)表象下,單電子密度為[20]

多電子波函數(shù)ψ(r,t)的模平方在全空間的積分為

(22)式中最后一個等式利用了軌道之間的正交性.如果軌道?p(t)是歸一化的,那么積分可以簡化為∫D(r,t)dr=Nα+Nβ=Ne,即等于體系中電子總數(shù),這與D(r,t)為單電子密度的定義是一致的.在數(shù)值計算中,由于使用了吸收邊界,軌道?p(t)不再是歸一化的.單電離概率的計算方法如下:首先取定一個截斷半徑r0,大于此半徑的電子密度被認(rèn)為是電離的,那么單電離概率可以定義為

可以看出,這里參數(shù)r0的選取決定了計算結(jié)果的準(zhǔn)確性.一方面r0要足夠大,以便將束縛態(tài)波包排除在r>r0區(qū)域外;另一方面,r0又不能太大,以便在計算的時間終點處電離態(tài)波包遠(yuǎn)離原子核運動到r>r0區(qū)域內(nèi).關(guān)于r0的選取見后文.

2.3.3 雙電離概率

與單電離概率類似,為了計算雙電離概率,首先定義雙電子密度[20],

雙電子密度在全空間的積分值為

同樣利用截斷半徑r0,雙電離概率計算為

注意在(23)和(26)式中,當(dāng)p?=q時??p(t)|?q(t)?r6r0可以不為零. 這與軌道之間保持正交性(在全空間積分)是不矛盾的.

2.3.4 自然軌道

MCTDHF理論中單電子軌道?p(t)沒有明顯的物理意義,分析中常常引入所謂的自然軌道.自然軌道的概念最早由 Per-Olov L?wdin在1955年引入,指的是對應(yīng)于一個多電子波函數(shù)的惟一的一組具有最大占據(jù)數(shù)的單電子波函數(shù)[29].首先對角化矩陣Dpq,得到本征值和相應(yīng)本征態(tài)為{pi},{χi}.那么,第i個自然軌道在坐標(biāo)表象下的波函數(shù)為

其中χij為χi的第i個分量.相應(yīng)的本征值pi即為占據(jù)數(shù),其取值范圍為0

3 結(jié)果與討論

3.1 計算參數(shù)

本文中的計算參數(shù)如下:H2原子核間距固定為R=1.4 a.u.(a=0.7 a.u.)(1 a.u.≈0.529 ?);坐標(biāo)ξ的范圍從1到一個截止值ξmax,本文設(shè)ξmax=rmax/a,其中rmax為球坐標(biāo)系下模擬盒子半徑的最大值;對于ξ和η,分別用了90個FEDVR基函數(shù)和18個DVR基函數(shù);量子數(shù)m的取值范圍為?6到6.激光電場脈沖具有Sine-Squared包絡(luò)的形式,

這里F0為峰值電場值,ω為載波頻率,T為脈沖長度.電場采用線偏振的形式,其偏振方向與分子軸的夾角記為θ.以rb表示開始采用吸收邊界的半徑值,在rb和rmax之間采用了吸收函數(shù)來防止波函數(shù)在邊界處導(dǎo)致的非物理反射.吸收函數(shù)采用如下形式:

盒子的半徑rmax必須要足夠大,使得物理系統(tǒng)不會受到吸收邊界的影響.隨著電離過程的發(fā)生,波函數(shù)也不在保持歸一.

3.2 關(guān)聯(lián)效應(yīng)對于分子基態(tài)總能量的影響

圖1給出虛時間演化時H2分子總能量隨時間的變化.組態(tài)空間分別為1σ(M=1),2σ(M=2),4σ(M=4) 和4σ1π(M=6). 這里的M表示空間軌道的個數(shù).1σ對應(yīng)于單組態(tài)情況,即Hartree-Fock(HF)方法.圖中兩條水平的線分別為HF方法(在相同格點參數(shù)下由DMTDHF程序計算[23])和精確結(jié)果(由雙電子TDSE計算[11]).可以看到,1σ的結(jié)果與DMTDHF的結(jié)果一致,表明本文方法的虛時演化在單組態(tài)情況是正確的.隨著組態(tài)數(shù)的增多,分子總能量逐漸向精確結(jié)果靠近,但在4σ1p時還沒完全收斂.隨著組態(tài)數(shù)的增多,達(dá)到收斂所需的演化時間不斷增大:在1σ時計算在t=20 a.u.(1 a.u.≈27.21 eV)左右收斂,而4σ1π時計算在t>100 a.u.才收斂.這表明組態(tài)數(shù)的增加,不但演化單位時間的計算量迅速增大,多體體系達(dá)到穩(wěn)定自洽狀態(tài)的復(fù)雜度也越大.

圖1 虛時間演化時H2分子總能量隨著演化時間的變化.組態(tài)空間分別為1σ,2σ,4σ,4σ1π;水平虛線表示HF 方法[23]和精確的雙電子TDSE計算[11]的能量值Fig.1.The total molecular energy of H2versus the imaginary propagation time,for different configuration spaces:1σ,2σ,4σ,4σ1π.The dashed horizontal lines represent the ground state molecular energies calculated by the Hartree-Fock method[23]and twoelectron TDSE method[11].

表1 由MCTDHF方法計算得到的H2分子基態(tài)的總能量隨著組態(tài)數(shù)的收斂情況.精確結(jié)果取自文獻(xiàn)[11]Table 1.Convergence of the ground state molecular energy of H2versus the configuration space from MCTDHF calculations.The exact result is taken from Ref.[11].

表1給出了H2分子基態(tài)總能量隨著組態(tài)數(shù)的收斂情況,并與準(zhǔn)確計算結(jié)果進(jìn)行了對比.為了定量衡量包含電子關(guān)聯(lián)的程度,對于某一組態(tài)空間X,定義關(guān)聯(lián)度Pc(X)如下,

其中E(1σ),E(X),Eexact分別為HF,MCTDHF(組態(tài)空間為X)和精確方法計算的能量值.關(guān)聯(lián)能的數(shù)值為Eexact?E(1σ)≈?0.0408 a.u.(?1.1 eV),占總能量的比重為3.47%.對于1σ組態(tài)空間(即HF),完全沒有關(guān)聯(lián)(Pc=0%).對于2σ組態(tài)空間,僅僅比1σ多了一個空間軌道,Pc增加到45.37%.當(dāng)組態(tài)空間為4σ1π時,已經(jīng)包含了90.48%的關(guān)聯(lián)能.隨著組態(tài)數(shù)的繼續(xù)增大,能量值收斂到精確結(jié)果的速度也變得越來越慢.當(dāng)組態(tài)空間為8σ1π時,其包含的組態(tài)數(shù)為144,遠(yuǎn)大于4σ1π的36個,然而其Pc值僅僅比后者高4.14%.要完全收斂到精確的能量值,還需要更多的組態(tài),特別是更高m量子數(shù)的態(tài)[如d(m=2)].

這里以6σ1π組態(tài)空間為例討論自然軌道.對于這個組態(tài)空間,表2中給出了8個自然軌道的占據(jù)數(shù),且本文計算的結(jié)果與文獻(xiàn)[21]的結(jié)果進(jìn)行了對比.需要注意的是,文獻(xiàn)的結(jié)果對于第6和第7個占據(jù)數(shù)結(jié)果在數(shù)值上有錯誤.這可以通過如下方式檢驗:所有占據(jù)數(shù)之和應(yīng)該等于電子的總數(shù)量(對于H2為2).本文計算的結(jié)果滿足上面這個條件.排除上述文獻(xiàn)結(jié)果數(shù)值上的錯誤,本文的結(jié)果與文獻(xiàn)結(jié)果符合得很好.從軌道占據(jù)數(shù)的大小可以看出,第1個自然軌道起主要作用,它實際上對應(yīng)于Hartree-Fock方法中的1σ軌道.其他幾個自然軌道的占據(jù)數(shù)都比較小,求和后的占據(jù)數(shù)占電子總數(shù)的百分比為1.8%.

表2 組態(tài)空間為6σ1π時H2分子基態(tài)的自然軌道占據(jù)數(shù),與文獻(xiàn)[21]的結(jié)果進(jìn)行了對比,結(jié)果由占據(jù)數(shù)從大到小的順序排列Table 2.Natural occupation numbers of H2molecule,when the configuration space is 6σ1π.The results are presented in a descending order.The reference values from Ref.[21]are also given for comparison.

3.3 不同傳播方法的對比以及電離概率的計算

本節(jié)討論電離概率的計算方法,所用的激光脈沖的長度為5個周期.在之前的工作中,我們基于TDHF方法,開發(fā)了處理雙原子電子動力學(xué)的DMTDHF程序[23].選取組態(tài)為1σ,將本文MCTDHF方法的結(jié)果與DMTDHF程序的結(jié)果進(jìn)行對比,可以測試在激光脈沖作用下實時演化時方法的正確性.圖2(a)給出了H2分子總能量隨著演化時間的變化.求解MCTDHF的軌道方程用了RK和BS兩種方法.這兩種方法給出了幾乎完全相同的結(jié)果.MCTDHF的結(jié)果與DMTDHF的結(jié)果也非常相近.DMTDHF中采用的是SIL傳播方法.結(jié)果的對比表明,程序的兩種傳播方法是正確的,能夠與DMTDHF取得一致.由于這里采用了自適應(yīng)精度的處理方法,圖2(b)給出了RK和BS中的動態(tài)時間步長.DMTDHF計算中采用的時間步長固定在0.01 a.u.(1 a.u.≈0.024 fs),因此不在圖中給出.RK方法中時間步長的變化較大,尤其是時間大于10 a.u.以后出現(xiàn)了許多振蕩.相比而言,BS方法中的時間步長則比較穩(wěn)定.這樣看來,BS方法應(yīng)該具有更好的穩(wěn)定性.RK方法相比于BS方法的優(yōu)越性在于計算量較小,在速度上一般比后者快1倍以上.實際使用中,出于計算量的考慮,一般首先使用RK方法.BS方法可以作為RK方法收斂性的驗證.

圖2 平行取向單組態(tài)(1σ)時不同傳播方法的對比:(a)H2分子總能量;(b)動態(tài)時間步長;其中TDHF為DMTDHF程序[23]計算的結(jié)果,時間步長固定為0.01 a.u.激光強(qiáng)度為1×1014W/cm2,頻率為2.0 a.u.Fig.2. Comparison of different propagation methods for parallel alignment and single configuration(1σ)of H2:(a)Total molecular energy;(b)the dynamic time steps.The TDHF result is calculated with the DMTDHF code[23],with a fixed time step of 0.01 a.u.The intensity and frequency of the laser pulses are 1×1014W/cm2and 2.0 a.u.,respectively.

對于電離概率的計算,由于(23)和(26)式依賴于半徑r0的選取,因此我們首先研究電離概率隨著r0的變化情況.圖3給出了當(dāng)計算盒子大小分別為rmax=40 a.u.和rmax=60 a.u.時的單電離概率和雙電離概率隨著r0的變化.可以看出,計算的結(jié)果對于r0依賴性很大.對于rmax=40 a.u.,電離概率在r0在30 a.u.和40 a.u.之間保持不變,這是因為這段空間存在著吸收邊界.同理對于rmax=60 a.u.,結(jié)果在50 a.u.和60 a.u.之間也是不變的.當(dāng)r0<5 a.u.時,電離概率計算結(jié)果隨著r0的增大而迅速減小.此時計算中在r>r0的區(qū)域內(nèi)包含了束縛態(tài)的成分,因此電離概率的計算是不準(zhǔn)確的.當(dāng)r0>20 a.u.時,結(jié)果也是隨著r0的增大而迅速減小的.這是因為電離的電子波包已經(jīng)遠(yuǎn)離原子核到達(dá)r0>20 a.u.的區(qū)域,過大的r0使得一部分電子波包被排除在計算之外.在r0=10 a.u.附近,電離概率計算結(jié)果變化比較平緩,可以比較好地定義電離概率.rmax=40 a.u.和rmax=60 a.u.的結(jié)果在r0=10 a.u.處幾乎重合,表明計算結(jié)果對于盒子大小已經(jīng)收斂.在下面的計算中,盒子大小都取為rmax=40 a.u..半徑r0的選取則根據(jù)上述討論的辦法,動態(tài)選擇電離曲線中間平臺區(qū)域的r0值.由于平臺區(qū)域變化比較緩慢,因此采用上述方法的電離概率結(jié)果是可靠的.同時,由于不需要等待所有電離波包都到達(dá)吸收邊界,因此大大縮短了模擬計算需要的收斂時間,從而減小了計算量.這對于利用MCTDHF計算電離概率提供了較大的便利.

圖3 由MCTDHF方法計算的平行取向單組態(tài)(1σ)時電離概率隨著截斷半徑r0的變化以及數(shù)值計算盒子大小的收斂性.實線和虛線分別為rmax=40 a.u.和rmax=60 a.u.時的電離概率;豎線表示r0=10 a.u.的位置;激光強(qiáng)度為1×1013W/cm2,頻率為1.0 a.u.Fig.3.Single and double ionization probabilities versus the cutof fradius r0from MCTDHF method,for parallel alignment and single configuration(1σ)of H2.The intensity and frequency of the XUV pulses are 1×1014W/cm2and 2.0 a.u.,respectively.The solid(dashed)lines are the results when the size of the simulation box rmax=40 a.u.(60 a.u.).The vertical line is placed at r0=10 a.u.

3.4 電子關(guān)聯(lián)效應(yīng)對于H2分子單電離以及雙電離的不同影響

當(dāng)激光頻率固定在1.0 a.u.(1 a.u.≈4.134×1016Hz)時,圖4給出了氫分子單電離和雙電離概率隨著激光強(qiáng)度的變化.激光強(qiáng)度范圍從5×1011W/cm2到2×1014W/cm2.激光脈沖的長度為5個周期.為了研究電子關(guān)聯(lián)在電離過程中的作用,同時給出了TDHF(1σ)和MCTDHF(4σ1π)的結(jié)果并進(jìn)行了比較.隨著激光強(qiáng)度的增大,單電離和雙電離概率都是單調(diào)遞增的.這里激光脈沖的中心頻率對應(yīng)的光子能量大于H2分子的電離能,但是小于H+2的電離能[30].因此單電離主要對應(yīng)于單光子電離過程,而雙電離則主要是雙光子電離過程的貢獻(xiàn).通過擬合可以發(fā)現(xiàn),單電離和雙電離分別近似符合一階和二階微擾理論,單電離和雙電離概率隨著光強(qiáng)的增大呈單調(diào)遞增,電離概率曲線的斜率分別近似為1和2.現(xiàn)在來看單組態(tài)和多組態(tài)計算結(jié)果的對比情況.對于單電離而言,4σ1π與1σ的結(jié)果符合得比較好.這表明電子關(guān)聯(lián)效應(yīng)對于單電離概率結(jié)果影響較小.因此,TDHF方法可以用于處理單電離過程,這驗證了在之前的工作中TDHF方法能夠比較好地重復(fù)出CO分子依賴于定向角度的單電離實驗結(jié)果[31].對于雙電離,多組態(tài)的結(jié)果與單組態(tài)有了較大的差別,在電離概率幅值上1σ為4σ1π的3倍左右.這表明電子關(guān)聯(lián)效應(yīng)在雙電離過程中發(fā)揮了重要的作用.

圖4 平行取向時H2分子單電離(SI)和雙電離(DI)概率隨著激光強(qiáng)度的變化.4σ1π的結(jié)果與1σ的結(jié)果進(jìn)行了對比,激光頻率為1.0 a.u.Fig.4. Single and double ionization probabilities of parallel aligned H2molecules versus the intensity of the XUV field.The results of 4σ1π (dashed lines)are compared with those of 1σ(solid lines).The frequency of XUV field is 1.0 a.u.

進(jìn)一步地,將激光強(qiáng)度固定為1×1013W/cm2,圖5給出了H2分子的單電離和雙電離概率隨著XUV脈沖頻率的變化.脈沖頻率范圍為0.6 a.u.到2.0 a.u..在這種情況下,隨著激光頻率的增大,單電離概率和雙電離概率都是單調(diào)減小的.這是由于電離截面隨XUV光子能量的增大而減小導(dǎo)致的.其中,多組態(tài)4σ1π與1σ的結(jié)果對比與圖4類似:單電離的結(jié)果基本重合,雙電離的結(jié)果則有比較大的差別.這進(jìn)一步驗證了電子關(guān)聯(lián)在單電離和雙電離過程中的不同作用.在圖4和圖5中,由TDHF方法計算得到的雙電離結(jié)果都與MCTDHF的結(jié)果有偏差.這是由于TDHF的單組態(tài)近似導(dǎo)致的.對于雙電子體系,在TDHF近似下雙電子密度dHF(r,r′)與單電子密度D(r)的關(guān)系可以表示為

將其代入雙電離概率的計算公式(26),得到

這表明在TDHF方法中,雙電離概率是直接依賴于單電離概率的.事實上,如果在圖4和圖5中由單電離結(jié)果根據(jù)(32)式畫出雙電離概率,結(jié)果與計算得到的雙電離概率完全一致.因此TDHF方法不能正確地給出雙電離過程的物理圖像.在MCTDHF方法中,隨著組態(tài)空間的增大以及組態(tài)數(shù)的增多,雙電子密度與單電子密度的相關(guān)度逐漸降低,因此雙電離結(jié)果也變得越來越偏離(32)式.本文計算的結(jié)果也支持上述結(jié)論.在圖4和圖5中,多組態(tài)的雙電離概率小于單組態(tài)的結(jié)果.在單組態(tài)方法中,雙電子被當(dāng)作平均場中的等價電子處理,即任意時刻都認(rèn)為這兩個電子等價.然而實際情況下,在單電離發(fā)生后,這兩個電子是不等價的(第二個電子束縛地更緊).這導(dǎo)致了單組態(tài)方法對于雙電離概率的高估.

圖5 平行取向時H2分子單電離(SI)和雙電離(DI)概率隨著激光頻率的變化.4σ1π的結(jié)果與1σ的結(jié)果進(jìn)行了對比,激光強(qiáng)度為1×1013W/cm2Fig.5. Single and double ionization probabilities of parallel aligned H2molecules versus the frequency of the XUV field.The results of 4σ1π (dashed lines)are compared with those of 1σ(solid lines).The intensity of XUV field is 1×1013W/cm2.

3.5 H2分子電離概率隨著取向角度的變化

尚未見利用三維的MCTDHF方法計算分子隨取向角度變化的電離概率的研究.圖6給出了H2分子單電離和雙電離概率隨著取向角度的變化.激光強(qiáng)度和頻率分別為1×1013W/cm2和1.0 a.u..根據(jù)之前的TDHF計算結(jié)果,對于800 nm的激光場,H2分子的單電離概率隨著取向角度從0?到90?是單調(diào)遞減的[23].在這里情況則不同,對于XUV脈沖,單電離概率和雙電離概率都是隨著取向角度單調(diào)遞增的.垂直取向(90?)時的電離概率大于平行取向(0?)時的電離概率,其比值p(90?)/p(0?)對于單電離和雙電離分別為1.5和2.6左右.不同激光場條件下的電離機(jī)理是不同的,這將導(dǎo)致電離概率在角度依賴上的定性差別[32].在800 nm激光場中,電離可能主要是非微擾條件下的隧穿電離過程,而在本文XUV激光場下,電離為微擾條件下的單光子和多光子電離過程.比如對于單光子電離,由于初態(tài)是σ軌道,單光子電離時在0?時末態(tài)只能是σ軌道,而在90?時末態(tài)為π軌道。二者哪一個更強(qiáng),依賴于頻率和對應(yīng)的矩陣元,在本文的激光條件下90?更強(qiáng)。

圖6 H2分子(a)單電離(SI)和(b)雙電離(DI)概率隨著取向角度的變化.4σ1π的結(jié)果與1σ的結(jié)果進(jìn)行了對比,激光強(qiáng)度為1×1013W/cm2,頻率為1.0 a.u.Fig.6.Single(a)and double(b)ionization probabilities of H2molecules versus the alignment angle.The results of 4σ1π (dashed lines)are compared with those of 1σ (solid lines).The intensity and frequency of XUV field are 1×1013W/cm2and 1.0 a.u.,respectively.

對比4σ1π和1σ的結(jié)果,主要差異為雙電離概率的幅值相差了3倍,而電離概率隨取向角度的依賴性相差不大.這可能是因為我們研究的H2分子體系比較簡單,僅有兩個電子.對于CO分子,內(nèi)殼層電子的動態(tài)極化極化會導(dǎo)致電離概率結(jié)果定性上的差別[31].對于具有更多電子的分子體系,多電子關(guān)聯(lián)也許會對電離概率的角度依賴性有更大的影響.

4 結(jié) 論

發(fā)展了三維的處理強(qiáng)激光場中的多電子雙原子分子電子動力學(xué)的多組態(tài)含時Hartree-Fock方法.該方法采用能夠精確處理雙中心庫侖勢的橢球坐標(biāo)系,以及能夠大大減小雙電子積分計算量的有限元-離散變量基函數(shù)方法.以H2為例討論了基態(tài)能量隨著組態(tài)數(shù)的收斂性.討論了自然占據(jù)軌道并與文獻(xiàn)結(jié)果進(jìn)行了對比.計算了在強(qiáng)XUV脈沖下H2分子的單電離和雙電離結(jié)果,討論了電子關(guān)聯(lián)效應(yīng)的作用.結(jié)果表明,電子關(guān)聯(lián)在單電離過程中影響不大,但是對于雙電離過程有明顯的影響.本文僅僅計算了總電離概率,而微分電離概率包含了更多的空間信息[6],電子關(guān)聯(lián)效應(yīng)的影響應(yīng)該更加明顯.對于電子關(guān)聯(lián)效應(yīng)的影響,還有待于更加深入的研究.當(dāng)電子數(shù)比較多時,MCTDHF方法中組態(tài)數(shù)呈指數(shù)形式增長,這成為限制MCTDHF方法應(yīng)用的主要障礙.后續(xù)有待于發(fā)展從完整的組態(tài)空間中選出和保留對于物理過程比較重要的組態(tài)的方法[15],以減小計算負(fù)擔(dān).

感謝大連化學(xué)物理研究所的李文亮博士關(guān)于MCTDHF方法的有益討論.

[1]Xu G X,Li L M,Wang D M,Chen M B 2009Quantum Chemistry:Fundamental Principle and Ab-initio Calculation Method(2nd Ed.)(Vol.1)(Beijing:Science Press)pp869–874(in Chinese)[徐光憲,黎樂民,王德民,陳敏伯2009量子化學(xué)——基本原理和從頭計算法(下)(第二版)(北京:科學(xué)出版社)第869—874頁]

[2]Becker W,Liu X J,Ho P J,Eberly J H 2012Rev.Mod.Phys.84 1011

[3]Zhao L,Zhang Q,Dong J W,Lü H,Xu H F 2016Acta Phys.Sin.65 223201(in Chinese)[趙磊,張琦,董敬偉,呂航,徐海峰2016物理學(xué)報65 223201]

[4]Jin F C,Wang B B 2016Acta Phys.Sin.65 224205(in Chinese)[金發(fā)成,王兵兵 2016物理學(xué)報 65 224205]

[5]Hu S X 2013Phys.Rev.Lett.111 123003

[6]Ye D F,Liu X,Liu J 2008Phys.Rev.Lett.101 233003

[7]Brabec T,Krausz F 2000Rev.Mod.Phys.72 545

[8]Posthumus J H 2004Rep.Progr.Phys.67 623

[9]Xiao X R,Wang M X,Li M,Geng J W,Liu Y Q,Peng L Y 2016Acta Phys.Sin.65 220203(in Chinese)[肖相如,王慕雪,黎敏,耿基偉,劉運全,彭良友 2016物理學(xué)報65 220203]

[10]Zhang Z,Peng L Y,Xu M H,Starace A F,Morishita T,Gong Q H 2011Phys.Rev.A84 043409

[11]Guan X,Bartschat K,Schneider B I 2011Phys.Rev.A83 043403

[12]Stapelfeldt H 2003Rev.Mod.Phys.75 543

[13]Kulander K 1987Phys.Rev.A36 2726

[14]Marques M A L,Ullrich C A,Nogueira F,Rubio A,Burke K,Gross E K U 2006Time-Dependent Density Functional Theory(Heidelberg,Berlin:Springer)pp1–13

[15]Hochstuhl D,Bonitz M 2012Phys.Rev.A86 053424

[16]Zanghellini J,Kitzler M,Fabian C,Brabec T,Scrinzi A 2003Laser Phys.13 1064

[17]Kitzler M,Zanghellini J,Jungreuthmayer C,Smits M,Scrinzi,Brabec T 2004Phys.Rev.A70 041401

[18]Caillat J,Zanghellini J,Kitzler M,Koch O,Kreuzer W,Scrinzi A 2005Phys.Rev.A71 012712

[19]Kato T,Kono H 2004Chem.Phys.Lett.392 533

[20]Hochstuhl D,Bonitz M 2011J.Chem.Phys.134 084106

[21]Haxton D J,Lawler K V,McCurdy C W 2011Phys.Rev.A83 063416

[22]Liao C T,Li X,Haxton D J,Rescigno T N,Lucchese R R,McCurdy C W,Sandhu A 2017Phys.Rev.A95 043427

[23]Zhang B,Yuan J M,Zhao Z X 2015Comput.Phys.Commun.194 84

[24]Lamb W E,Schlicher R R,Scully M O 1987Phys.Rev.A36 2763

[25]Abramowitz M,Stegun I A 1972Handbook of Mathmatical Functions(Washington:Dover)p752

[26]Rescigno T,McCurdy C 2011Phys.Rev.A62 032706

[27]Press W H,Teukolsky S A,Vetterling W T,Flannery B P 2007Numerical Recipes(3rd Ed.)(New York:Cambridge University Press)pp179–193,899–928

[28]Park T J,Light J C 1986J.Chem.Phys.85 5870

[29]Weinhold F,Landis C R 2001Chem.Educ.Res.Pract.2 91

[30]Computational Chemistry Comparison and Benchmark DataBase,National Institute of Standards and Technology http://cccbdb.nist.gov/expdata.asp[2017-12-20]

[31]Zhang B,Yuan J M,Zhao Z X 2013Phys.Rev.Lett.111 163001

[32]Zhang B,Yuan J M,Zhao Z X 2012Phys.Rev.A85 033421