基于Numerov算法的原子束縛態(tài)波數(shù)的數(shù)值求解

郭志堅， 李小明， 鄭興榮， 孫 乾

(1.隴東學(xué)院 電氣工程學(xué)院，甘肅 慶陽 745000； 2.環(huán)縣職業(yè)中等專業(yè)學(xué)校，甘肅 慶陽 745000)

隨著激光技術(shù)的飛速發(fā)展，激光與原子分子的相互作用成了熱門的研究領(lǐng)域.在強(qiáng)激光與原子分子的相互作用中，電離是最基本的過程[1-2]. 在原子強(qiáng)場電離理論的研究中，計算原子或分子的基態(tài)及激發(fā)態(tài)波函數(shù)是一項必要的工作[3-5]. 由于粒子所處的勢場通常很復(fù)雜，迄今為止能解析求解的勢場形式非常少，其相應(yīng)的波函數(shù)很難從理論上通過解薛定諤方程得到，大部分都要通過近似方法或者數(shù)值求解得到.近年來，由于計算機(jī)技術(shù)的快速發(fā)展，數(shù)值求解微分方程得到了廣泛應(yīng)用，例如，求解一階微分方程的歐拉法、求解二階微分方程的龍格庫塔法等[5-7].微觀粒子滿足的薛定諤方程是二階微分方程，本文利用Numerov算法[8-10]對原子定態(tài)薛定諤方程進(jìn)行求解，采用矩陣解法將哈密頓量簡化為三對角矩陣，然后通過矩陣對角化求解本征值及本征函數(shù)，不僅大大提升了計算效率，而且求解精度非常高. 另外，引入原子模型勢可以對不同原子的波函數(shù)進(jìn)行求解.除氫原子的波函數(shù)外，本文還求解了He、Ne、Ar原子的波函數(shù)，擴(kuò)大了其應(yīng)用范圍.

1 電子在庫侖場中的運動

在氫原子或類氫原子體系中，電子在原子核的庫侖場中運動的勢能為

(1)

這時的定態(tài)薛定諤方程為

(2)

球坐標(biāo)系下的拉普拉斯算子為

(3)

定態(tài)薛定諤方程的解為

(4)

式中:R(r)為徑向波函數(shù)，Ylm(θ,φ)為球諧函數(shù)，其中l(wèi)和m分別表示軌道量子數(shù)和磁量子數(shù).徑向波函數(shù)滿足下列方程

(5)

解得

(6)

其中:ρ=2Zr/na0,a0是玻爾最小軌道半徑，n和l表示主量子數(shù)以及軌道量子數(shù).

2 原子波函數(shù)的數(shù)值求解——Numerov算法

求解氫原子的定態(tài)波函數(shù)非常方便，但在許多情況下需要求解高激發(fā)態(tài)的波函數(shù)，用解析解處理十分不便.另外，對于多個電子體系(除類氫原子之外)，其勢能比較復(fù)雜，很難精確求解，數(shù)值解法是常用的一種方法. 徑向方程(1)可改寫為

ER(r).

(7)

u(r)=rR(r)，

(8)

則(7)式變?yōu)?/p>

(9)

我們定義一個有效勢Vl(r)，則(9)式進(jìn)一步化簡為

(10)

(10)式在形式上與一維定態(tài)薛定諤方程相似，只是其中的勢能換成了有效勢.由于(10)式為一維方程，可改寫為

(11)

根據(jù)泰勒公式可得

(12)

(13)

(12)式與(13)式相加得

(14)

對于束縛態(tài)問題，邊界條件為ψn(0)=ψn(X)=0.將最大徑向距離X分為N+1份，步長為s，即X=(N+1)s，結(jié)合(12)～(14)式得

-ψn(x+s)-ψn(x-s)+2ψn(x)+

(15)

令λn=s2[2mEn/?2]，b(x)=2+s2[2mV(x)/?2]，則(15)式改寫為

-ψn(x+s)-ψn(x-s)+b(x)ψn(x)=λnψn(x).

(16)

將ψn(s),ψn(2s),ψn(3s),…,ψn(is)簡寫為vn(1)，vn(2)，vn(3)，…,vn(i)，令b(is)=bi，得到下列方程組：

b1·vn(1)+(-1)·vn(2)+

0·vn(3)+0·vn(4)+…=λnvn(1),

(-1)·vn(1)+b2·vn(2)+

(-1)·vn(3)+0·vn(4)+…=λnvn(2),0·vn(1)+(-1)·vn(2)+

(17)

b3·vn(3)+(-1)·vn(4)+…=λnvn(3),

0·vn(1)+0·vn(2)+

(-1)·vn(3)+b4·vn(4)+…=λnvn(4),……

方程組(17)可表示為矩陣形式:

Mvn=λnvn，

(18)

其中

(19)

M為典型的三對角矩陣，通過矩陣對角化可以求出能量本征值λn及相應(yīng)的本征函數(shù)vn.至此，將解一維定態(tài)薛定諤方程問題轉(zhuǎn)化為求矩陣的特征值和特征向量問題，從而求解束縛態(tài)波函數(shù)及能級就非常方便.

3 結(jié)果與分析

在具體的數(shù)值計算過程中，為使求解方便，統(tǒng)一采用原子單位(atomic units，a.u.).對于能量，1 a.u.=27.2 eV. 為了求解更多原子的束縛態(tài)波函數(shù)，引入模型勢[5]：

(20)

除氫原子外，本文還數(shù)值求解了He、Ne、Ar的束縛態(tài)波函數(shù)，模型勢的具體參數(shù)見表1.

表1 原子的模型勢參數(shù)及電離能

4種原子的模型勢隨徑向距離改變的勢能曲線見圖1.從圖1可以看出，由于原子核與核外電子的影響， H、He、Ne、Ar的模型勢曲線有明顯差異，但隨著徑向距離增大，當(dāng)r>5 a.u.時均開始趨于0.

圖1 H、He、Ne、Ar的模型勢曲線圖

3.1 氫原子數(shù)值解與解析解

(21)

在數(shù)值求解中引入徑向波函數(shù)u(r)=rR(r)，這4個態(tài)徑向波函數(shù)的數(shù)值解和解析解見圖2.從圖2可以看出，數(shù)值解與解析解非常吻合，徑向距離在40 a.u.以內(nèi)完全吻合.這說明數(shù)值解的精度很高，采用Numerov算法是可取的.同時程序運行非?？?，說明該算法的效率很高.圖2中1s態(tài)的最大徑向距離最小，r=8 a.u.處波函數(shù)趨于0，2s、2p態(tài)中2s態(tài)范圍小于2p態(tài)，3p態(tài)徑向距離最大，大約為r=20 a.u.這是因為隨著n和l增加，電子具有更大的能量，所以能夠運動得更遠(yuǎn).

圖2 氫原子徑向波函數(shù)的解析解與數(shù)值解

電子在空間徑向出現(xiàn)的概率用徑向概率密度函數(shù)表示.電子在區(qū)間[r,r+dr]上出現(xiàn)的概率可表示為

(22)

其中：u2即徑向概率密度，u2最大值所對應(yīng)的r值為原子的最概然半徑(量子力學(xué)中最概然半徑指的是電子在空間中可能出現(xiàn)的最大球殼半徑).氫原子1s、2s、3s態(tài)及2p、3p和 4p態(tài)的徑向概率密度函數(shù)的數(shù)值計算結(jié)果見圖3，其中每條曲線的最大峰值對應(yīng)的r為原子的最概然半徑.從圖3可以看出波函數(shù)徑向概率分布的節(jié)點數(shù)滿足nr=n-l-1(不包括邊界點)，即1s、2s、3s態(tài)的節(jié)點數(shù)分別為0、1、2，而2p、3p 、4p態(tài)的節(jié)點數(shù)也分別為0、1、2.

圖3 氫原子徑向概率分布的數(shù)值解

3.2 He、Ne、Ar原子的波函數(shù)

利用模型勢給出了He原子1s態(tài)和2p態(tài)的徑向波函數(shù)(圖4).顯然，2p態(tài)比1s態(tài)的徑向距離大.

圖4 氦原子徑向波函數(shù)的數(shù)值解

圖5為Ne原子的3s態(tài)和2p態(tài)波函數(shù)，其中基態(tài)為2p態(tài)，徑向波函數(shù)節(jié)點數(shù)為0，而3s態(tài)節(jié)點數(shù)為2. 圖6為Ar原子的3p態(tài)和4s態(tài)波函數(shù)，3p態(tài)為基態(tài)，4s態(tài)為第一激發(fā)態(tài)，其中3p態(tài)節(jié)點數(shù)為1，4s態(tài)節(jié)點數(shù)為3. 根據(jù)量子理論，節(jié)點處電子出現(xiàn)的概率為0，即電子在該處不出現(xiàn).

圖5 氖原子徑向波函數(shù)的數(shù)值解

圖6 氬原子徑向波函數(shù)的數(shù)值解

將氫原子、氦原子、氖原子及氬原子的基態(tài)波函數(shù)H(1s)、He(1s)、Ne(2p)與Ar(3p)作比較，結(jié)果見圖7. 從圖7可以明顯地發(fā)現(xiàn)Ne與He的基態(tài)波函數(shù)非常相似，其趨于0的形式也保持一致；而H與Ar的基態(tài)波函數(shù)在r=1 a.u.處相交，此后其形式較為接近，在r>1 a.u.的區(qū)域，由于Ar原子基態(tài)波函數(shù)存在節(jié)點，從而有了區(qū)別，但總體而言，H與Ar的基態(tài)波函數(shù)具有較高的相似性.由表1可知, H與Ar的電離能分別為0.5 a.u.和0.579 a.u.,十分接近; 而Ne與He的電離能分別為0.793 a.u.和0.904 a.u.，也比較接近. 這說明原子的基態(tài)波函數(shù)與原子電離能的相關(guān)性很高.

圖7 H、He、Ne、Ar的基態(tài)波函數(shù)

4 結(jié)論

本文詳細(xì)介紹了氫原子束縛態(tài)徑向波函數(shù)的解析求解過程，在此基礎(chǔ)上引入Numerov算法將氫原子束縛態(tài)波函數(shù)的求解數(shù)值化，借助于FORTRAN和MATLAB軟件可以輕松求解氫原子的束縛態(tài)波函數(shù).本文將氫原子的束縛態(tài)波函數(shù)(1s、2s、2p、3p)的數(shù)值解與解析解作了對比，結(jié)果非常吻合.除此之外，我們引入了模型勢，可以對不同的原子束縛態(tài)進(jìn)行求解.本文給出了He、Ne、Ar原子的模型勢參數(shù)及電離能，數(shù)值求解了He、Ne、Ar原子的部分束縛態(tài)波函數(shù).最后，對比了H、He、Ne、Ar原子的基態(tài)波函數(shù)，發(fā)現(xiàn)H與Ar基態(tài)波函數(shù)較為相似，He與Ne基態(tài)波函數(shù)較為接近，結(jié)果表明原子基態(tài)波函數(shù)與電離能有直接關(guān)系.