利用重歸一化Numerov 方法研究超冷雙原子碰撞

白凈 謝廷

1)(大連職業(yè)技術(shù)學(xué)院電氣電子工程學(xué)院,大連 116035)

2)(中國科學(xué)院大連化學(xué)物理研究所,分子反應(yīng)動力學(xué)國家重點實驗室,大連 116023)

采用重歸一化Numerov 算法求解關(guān)于超低溫雙原子碰撞問題的非含時薛定諤方程組.以39K-133Cs 碰撞為例,研究了超低溫下雙原子Feshbach 共振的性質(zhì).結(jié)果表明,重歸一化Numerov 算法可以很精確地描述超冷條件下碰撞過程.與改進(jìn)的logarithmic derivative 算法相比,在同等參數(shù)條件下,重歸一化Numerov 方法在計算效率上雖然有一定劣勢,但在大格點步長參數(shù)范圍內(nèi)有著更好的穩(wěn)定性.提出重歸一化Numerov 和logarithmic derivative 算法相結(jié)合的傳播方法,在保證結(jié)果精度的同時大大減少了計算時間.此項算法也可以應(yīng)用于求解任意溫度下的兩體碰撞耦合薛定諤方程組.

1 引言

在過去幾十年里超冷分子的制備和應(yīng)用經(jīng)歷了長足發(fā)展.在研究超低溫原子和分子氣體相互作用的實驗中,Feshbach 共振扮演著一個主要角色[1].由于在超低溫下碰撞粒子對的動能極低,這使得我們可以通過磁場調(diào)節(jié)碰撞粒子不同通道之間的能量差.當(dāng)粒子對在入射通道的散射態(tài)和閉通道中束縛態(tài)能量發(fā)生簡并時,即為Feshbach 共振發(fā)生,表征原子間相互作用的散射長度會發(fā)散至無窮大[2].同時,入射通道的散射波函數(shù)在短程處增強,碰撞粒子對會迅速轉(zhuǎn)換為弱束縛分子態(tài)[3].由于不同通道間能量差和粒子的性質(zhì)有關(guān),Feshbach 共振只發(fā)生在一些特定強度的磁場下.許多實驗小組利用磁場來研究和控制超冷氣體的動力學(xué)和結(jié)構(gòu)性質(zhì)[4].

當(dāng)前有幾個研究小組發(fā)展了不同的計算方法來預(yù)測磁場中Feshbach 共振的性質(zhì),主要分為以下3 種:1)量子虧損理論(quantum defect theory).通過把相互作用勢的散射長度作為輸入?yún)?shù)并忽略精確的短程相互作用勢,Gao[5]發(fā)展了解析的量子虧損模型,使其可求解超冷碰撞原子對的散射性質(zhì)及磁誘導(dǎo)Feshbach 共振.2)漸近束縛態(tài)模型(asymptotic bound state model).Tiecke 等[6]利用單三重態(tài)的弱束縛態(tài)的能量本征值作為輸入?yún)?shù),發(fā)展了漸近束縛態(tài)模型,使其可以預(yù)測磁誘導(dǎo)Feshbach共振的位置及寬度.3)密耦方法(close-coupling method).通過數(shù)值求解完整哈密頓的薛定諤方程組,可以獲得精確的散射信息[7].這是一種沒有使用任何近似的方法,因此也常被用來解釋實驗觀測結(jié)果.

對于精確求解密耦方程,人們發(fā)展了幾種不同的數(shù)值方法,如對數(shù)導(dǎo)數(shù)方法(logarithmic deriviative,LOGD)[8]、R 矩陣方法[9]、龍格-庫塔法[10]、Numerov方法[11]等.其中,LOGD 被認(rèn)為是最穩(wěn)定也很高效的方法,常用來作為數(shù)值求解耦合薛定諤方程組的標(biāo)準(zhǔn).Numerov 方法可以給出散射態(tài)和束縛態(tài)的具體波函數(shù),多用于單通道問題.由于超低溫下磁場耦合了不同的束縛態(tài)通道,Numerov 方法傳播的波函數(shù)會發(fā)散至無窮大,并不適合用來研究共振過程[12],為此人們發(fā)展了更加穩(wěn)定的重歸一化Numerov方法 (renormalized Numerov,RN)[13].Colavecchia 等[14]利用改進(jìn)的RN 方法研究了溫度大于1 K 時的三原子分子碰撞反應(yīng)散射,發(fā)現(xiàn)了RN 方法在研究碰撞解離時比LOGD 方法具有更好的收斂性.Karman 等[15]推導(dǎo)了可以使用變步長的RN 方法,比較了不同算法的精度,研究了溫度為1 K 時的分子-分子散射碰撞,并給出了不同通道的散射波函數(shù).Blandon 等[16]將映射格點方法和RN 方法相結(jié)合,改進(jìn)了計算效率.然而關(guān)于RN 方法在超低溫下雙原子散射碰撞應(yīng)用并未被詳細(xì)研究.在這項工作中,使用RN 方法來研究超低溫雙原子碰撞過程.

對于堿金屬雙原子分子體系,處于絕對基態(tài)的超冷KCs 分子具有較大的永久電偶極矩,在低溫下有很強的長程偶極相互作用,是許多實驗組的研究熱點[17,18].本文以39K-133Cs 兩體碰撞為例,研究了超低溫下的Feshbach 共振性質(zhì),并與LOGD 方法作比較.結(jié)果表明RN 方法可以很精確地給出散射信息,并且在大格點步長情況下有著比LOGD更好的穩(wěn)定性.在此基礎(chǔ)上,提出RN 和LOGD 相結(jié)合的傳播方法,在保證精度的同時可以大大減少計算時間.此項工作為研究散射碰撞提供了額外一種可靠的數(shù)值方法.

2 計算方法

2.1 外磁場中堿金屬雙原子碰撞Schr?dinger方程

外磁場中超低溫雙原子碰撞的總哈密頓可以表示為

式中,μ是雙原子體系的折合質(zhì)量,R是原子核間距,l是雙原子體系的軌道角動量量子數(shù),θ和φ分別表示空間極化角和方位角.(1)式右邊第一項表示相對運動動能,第二項是離心勢能.(1)式右邊第三項是原子間的相互作用勢,可表示為

式中,VS(R)表示碰撞原子對在總自旋態(tài)S的絕熱相互作用,MS是S沿著磁場方向的投影量子數(shù),這里假設(shè)磁場方向和空間量子z軸平行.

式中,μ0表示玻爾磁子,B是磁場強度,μα(β)表示原子α(β)的核磁矩,Iα(β)和分別表示原子α(β)的核自旋及其在磁場方向的投影量子數(shù).

式中,γα(β)表示原子α(β)的超精細(xì)相互作用常數(shù),Sα(β)是原子α(β)的電子自旋量子數(shù).

定義完全非耦合基矢

式中,ψ(R)是體系的總波函數(shù)矩陣形式,V(R)是勢能相互作用矩陣,E是碰撞粒子對的總能量矩陣.求解此方程組,可以獲得散射碰撞信息等參量.

2.2 多通道重歸一化Numerov 方法

Numerov 方法是一種可以用于求解一維薛定諤方程的有效方法.在數(shù)值計算中,通常把原子核間距R離散為等間距的格點,然后求解不同格點位置下的波函數(shù).對于單通道情況,其計算核心是三點迭代關(guān)系式[13]

這里ψn表示在第n個格點處的波函數(shù);Tn=—(μh2/6)(Vn—E),其中h為格點之間的間隔,Vn表示勢能在第n個格點的值,E是入射通道的碰撞能.

通過對(6)式作兩次變換,可得RN 傳播算法.具體如下:在第一步變換中定義Fn=(1—Tn)ψn,并代入(6)式中,可得

式中Un=(2+10Tn)/(1—Tn).第二步變換是將三點迭代關(guān)系替換為兩點迭代.首先定義比率Qn=Fn+1/Fn,將其代入(7)式中,則得兩點迭代關(guān)系式

進(jìn)一步地,可得波函數(shù)在第n個格點處的導(dǎo)數(shù)

式中An=1/2(1—2Tn)/(1—Tn).由于Qn定義的是前后相鄰格點的波函數(shù)比例,因此被稱作重歸一化方法.注意上面給出的單通道LOGD 參數(shù)是在RN方法中的迭代關(guān)系式,不同于文中所說的LOGD 方法.對于超低溫雙原子散射碰撞,通常入射通道和其他多個通道耦合在一起,(10)式應(yīng)記為矩陣形式.為了完整性,寫出多通道LOGD 參數(shù)的矩陣形式

式中,I為單位方陣,An=1/2(I— 2Tn)(I—Tn)—1.在數(shù)值計算中,一般設(shè)置初始格點Qn=1=1020×I,將(10)式由短程向長程傳播,可以得到y(tǒng)矩陣在第n個格點的值.當(dāng)原子核間距Rn足夠大時,碰撞原子對的波函數(shù)可表示為[19]

式中,K是和Rn無關(guān)的反應(yīng)矩陣,J和N都是對角矩陣.對于開通道情況,J和N的矩陣元可表示為Riccati-Bessel 函數(shù)[20]

對于閉通道情況,其矩陣元可以表示為修正的第一和第二類球Bessel 函數(shù)[20]

式中kj和lj分別表示第j個通道的波矢和分波量子數(shù).對(12)式求一階導(dǎo)再除以(12)式,可得K矩陣的表達(dá)式

從而得到S矩陣[21]

式中i 表示虛數(shù)單位.

在獲得S矩陣后,就可以得到散射長度等信息.對于s 波通道來說,其散射長度a可用下式來計算[22]

式中,δ(k)表示第j個通道的相移,Sjj(k)是S矩陣中第j個通道的對角矩陣元.

3 結(jié)果與討論

本部分研究了39K-133Cs 超低溫兩體碰撞.由于碰撞能極低,高階分波對散射的貢獻(xiàn)很小,可以只考慮最低的入射s 波通道.雖然偶極-偶極相互作用可以耦合s 和d 波,d 波可能誘導(dǎo)新s 波共振[23],但是本文關(guān)注的是數(shù)值方法的準(zhǔn)確性和計算效率,因此,在計算中偶極相互作用被忽略.

采用Ferber 等[24]在2013 年由光譜數(shù)據(jù)擬合的基單重態(tài)和三重態(tài)的勢能面數(shù)據(jù).圖1 展示了39K-133Cs 雙原子體系的基態(tài)勢能相互作用.在這個體系中能量最低的電子態(tài)解離極限是39K(4s)+133Cs(6s),相關(guān)的電子態(tài)為基單重態(tài)X1Σ+和a3Σ+,對應(yīng)的散射長度分別是—18.37a0和82.24a0,a0是玻爾半徑.在求解耦合微分方程組時選定短程起點為R=3a0,長程終點為R=2000a0,每個格點間的間隔固定為0.002a0,碰撞能為1 nK.在這樣的低能區(qū)域,s 波散射性質(zhì)和碰撞能無關(guān).

圖1 39K-133Cs 體系基單重態(tài)和三重態(tài)相互作用勢Fig.1.Ground electronic potential curves of singlet and triplet states in 39K-133Cs.

由于39K 和133Cs 原子最外層各有1 個電子,所以它們的電子自旋量子數(shù)都是1/2,核自旋量子數(shù)分別為3/2 和5/2,在磁場中能量最低的超精細(xì)態(tài)通道為|fK=1,mfK=1〉+|fCs=3,mfCs=3〉,這里f和mf分別表示單個原子的總自旋量子數(shù)及其在磁場方向的投影.在磁場中,碰撞粒子對的總自旋投影量子數(shù)Mf=mfK+mfCs守恒.當(dāng)入射通道為s 波|fK=1,mfK=1〉+|fCs=3,mfCs=3〉時,它是唯一開通道,并和另外7 個Mf=4 的閉通道通過塞曼超精細(xì)相互作用耦合在一起.圖2 為入射通道為|fK=1,mfK=1〉+|fCs=3,mfCs=3〉時的s 波散射長度,可以看到散射長度在某些磁場位置展示了發(fā)散的特性,這表明發(fā)生了共振.因為改變磁場強度同時也改變了開通道和閉通道間的能量差,從而使得閉通道的束縛態(tài)可以和開通道散射態(tài)發(fā)生簡并.圖中展示的5 個理論共振位置也被其他研究小組預(yù)測到[22],而實驗上人們只發(fā)現(xiàn)了位于341.9 G (1 G=10—4T)和421.6 G 處的共振信號[25],另外3 個共振并沒有被觀測到.需要指出的是,正無窮大的散射長度表明原子間相互作用是強排斥的,反之是強吸引的.在實驗上需要制備某些種類的超冷分子時,人們總是調(diào)節(jié)磁場由強吸引區(qū)域到強排斥區(qū)域,這樣才會得到弱束縛的超冷分子,然后通過受激拉曼絕熱通道技術(shù)轉(zhuǎn)移至絕對基態(tài)就可以獲得絕對基態(tài)的超冷分子[1].在單個共振附近,散射長度a可以表示為[26]

圖2 在入射通道計算的磁場0—1000 G 范圍內(nèi)的s 波散射長度,紅線標(biāo)記錄5 個共振位置Fig.2.Calculated s-wave scattering length in the incoming channel with field range 0 —1000 G.The five resonant positions are labeled with red lines.

式中,abg是背景散射長度;ΔB是共振寬度,定義為散射長度為0 和無窮大時磁場強度之差;B0是共振位置,即散射長度為無窮大的磁場值.(18)式中比較重要的參數(shù)是共振寬度,它是衡量共振強度的參量.由于實驗裝置中磁場強度具有一定的不確定度,因此寬共振在實驗上容易被發(fā)現(xiàn),也更有應(yīng)用價值[27].根據(jù)(18)式,對計算的散射長度在所有共振附近進(jìn)行了擬合,并把參數(shù)列于表1 中.作為對比,利用Manolopoulos[28]發(fā)展的改進(jìn)的定步長LOGD 方法得到的結(jié)果也被列入表中.

從表1 可以看到,由兩種方法得到的共振參數(shù)幾乎一致,這表明了RN 方法的計算精度很高,可以和LOGD 方法相比擬.此外,也可以看到大部分s 波共振普遍偏窄(＜1 G),這是由于39K-133Cs體系的單三重態(tài)的散射長度差異較小(—18.37a0和82.24a0),不同通道之間的耦合較弱導(dǎo)致的.

表1 兩種計算方法得到的共振位置、共振寬度和背景散射長度參量Table 1. Calculated resonant positions,widths and background scattering lengths in the incoming channelby both methods.

表1 兩種計算方法得到的共振位置、共振寬度和背景散射長度參量Table 1. Calculated resonant positions,widths and background scattering lengths in the incoming channelby both methods.

接下來比較了兩種方法的計算效率.圖3(a)給出了利用兩種方法計算散射長度的程序運行時間隨格點數(shù)的變化.由于單一磁場強度下兩種方法計算散射長度的時間均很短,存在較大誤差,為了比較的準(zhǔn)確性,這里列出了計算100 個磁場強度下散射長度的總時間.計算中所有的參數(shù)都選取一致,并且都使用了固定步長方法.格點數(shù)變化范圍為2×105到2×106,對應(yīng)的格點步長范圍為10—2a0至10—3a0.可以看到LOGD 方法在計算效率上有一定優(yōu)勢,然而在少格點數(shù)情況下差異縮小,這主要是由于在RN 方法中(8)式中的Un是精確計算得來的,而改進(jìn)的LOGD 方法使用了參考勢的解析解代替精確解[28],在對格點循環(huán)過程中少了一步矩陣相乘運算.如果采用和LOGD方法的近似處理方法,原則上這兩者的計算時間應(yīng)該很接近,這也是當(dāng)前RN 方法可以繼續(xù)優(yōu)化的一方面.圖3(b)給出了磁場為1 G 時隨總格點數(shù)變化的散射長度.可以看到LOGD 方法在少格點數(shù)情況下計算誤差增大,而RN 方法在少格點數(shù)(大格點步長)情況下表現(xiàn)的更穩(wěn)定.這個現(xiàn)象在任意磁場強度下都可以被觀測到,表明使用RN 方法時可以采用較少的格點數(shù)來減少計算時間.

圖3 (a)采用兩種方法計算100 個磁場強度下散射長度的計算時間隨格點數(shù)的變化;(b)在磁場為1 G 時散射長度隨格點數(shù)的變化Fig.3.(a)Computational time of scattering length at 100 magnetic fields as a function of number of sectors by using two methods;(b)the scattering lengths as a function of number of sectors at 1 G.

第2 節(jié)介紹了RN 的計算基礎(chǔ)是泰勒展開,高階展開項被忽略(h＞6),也被稱為“解跟隨”方法,適合用于勢能變化比較劇烈的短程區(qū)域.而改進(jìn)的LOGD 方法使用了近似的參考勢,是一種“勢能跟隨”方法,適合勢能變化比較平緩的長程區(qū)域.根據(jù)上面的計算結(jié)果提出一種新傳播方法:在雙原子核間距短程處使用RN 方法,在長程處使用變步長LOGD 方法.這種方法可以利用兩者各自的優(yōu)點,從而最大程度地減少計算時間.表2 列出使用上述新傳播方法得到的5 個共振位置總誤差.這里設(shè)定短程和長程接合點為R=20a0,即R＜20a0范圍使用RN 算法,R＞20a0范圍使用LOGD 算法,注意更大的接合點會降低計算效率,而更小的接合點會引起較大的計算誤差.可以看到當(dāng)RN 方法在短程傳播步長小于 0.008a0時,5 個共振位置的總誤差小于1 G,仍可以得到較為精確的共振位置.和圖3(a)相比計算總步數(shù)大大減少,計算效率很高.預(yù)期此傳播方法應(yīng)用在復(fù)雜體系可以減少計算時間,擴(kuò)大研究體系范圍.

表2 采用短程RN 和長程變步長LOGD 相結(jié)合的傳播方法得到的共振位置總誤差,這里兩種方法的接合點在R=20a0.第一列ΔR 表示在RN 方法中使用的固定格點步長,第二列NRN 表示RN 方法在短程傳播的總步數(shù),第三列NLOGD 表示在變步長LOGD 方法在長程傳播的總步數(shù),第四列表示使用兩者結(jié)合方法產(chǎn)生的所有共振位置誤差的絕對值總和Table 2. Total errors of resonant positions by using the method combining RN method in short range with LOGD with variable step size in long range,where the connected point is located at 20a0.The first column,ΔR,represents the fixed step size in RN method.The second and third columns,NRN and NLOGD,denote the steps propagated in RN and LOGD methods,respectively.The last column is the sum of errors for all resonant positions.

表2 采用短程RN 和長程變步長LOGD 相結(jié)合的傳播方法得到的共振位置總誤差,這里兩種方法的接合點在R=20a0.第一列ΔR 表示在RN 方法中使用的固定格點步長,第二列NRN 表示RN 方法在短程傳播的總步數(shù),第三列NLOGD 表示在變步長LOGD 方法在長程傳播的總步數(shù),第四列表示使用兩者結(jié)合方法產(chǎn)生的所有共振位置誤差的絕對值總和Table 2. Total errors of resonant positions by using the method combining RN method in short range with LOGD with variable step size in long range,where the connected point is located at 20a0.The first column,ΔR,represents the fixed step size in RN method.The second and third columns,NRN and NLOGD,denote the steps propagated in RN and LOGD methods,respectively.The last column is the sum of errors for all resonant positions.

4 結(jié)論

本文利用RN 和LOGD 數(shù)值方法計算了K-Cs體系超低溫下的散射性質(zhì).結(jié)果表明RN 方法可以用于精確求解超低溫下的兩體散射過程,并且在少格點數(shù)情況下的穩(wěn)定性允許人們使用更大的格點步長,從而減少計算時間.另外,提出了短程固定步長RN 和長程變步長LOGD 相結(jié)合的傳播方法以提高計算效率.此方法對于求解任意溫度下的兩體碰撞問題同樣有效.