謝小濤 董元浩 王建康

【摘要】束縛量子系統(tǒng)的本征值精確可解模型稀少，如何利用數(shù)值方法得到系統(tǒng)的本征值是初等量子力學(xué)教學(xué)過程中不得不面對(duì)的問題之一。本文從線性諧振子模型出發(fā)，借助諧振子本征值態(tài)所建構(gòu)的希爾伯特空間，利用Matlab提供的庫(kù)函數(shù)求解任意束縛勢(shì)量子系統(tǒng)的本征值和本征函數(shù)。數(shù)值方法的講授將極大提高學(xué)生靈活應(yīng)用量子力學(xué)知識(shí)解決實(shí)際問題的能力。

【關(guān)鍵詞】諧振子 束縛態(tài) 薛定諤方程 Matlab

【中圖分類號(hào)】O413 【文獻(xiàn)標(biāo)識(shí)碼】A 【文章編號(hào)】2095-3089（2018）03-0048-02

引言

薛定諤方程束縛態(tài)離散譜的計(jì)算是初等量子力學(xué)所研究的基本問題之一。對(duì)于定態(tài)薛定諤方程本征值求解問題，由于可解的模型稀少，數(shù)值方法成為必不可少的解決方案。數(shù)值計(jì)算在量子力學(xué)教學(xué)方面如何取得突破是值得探討的問題之一，其能有效幫助學(xué)生形成更加清晰的物理圖像并能大幅提高學(xué)生初步掌握并解決實(shí)際問題的能力。本文不關(guān)注上述數(shù)值方法，而側(cè)重于如何基于線性諧振子模型本征波函數(shù)所建構(gòu)的正交完備基求解系統(tǒng)的本征值及其相應(yīng)的本征函數(shù)。正交完備基所組成的希爾伯特相關(guān)性質(zhì)在初等量子力學(xué)教學(xué)過程均會(huì)進(jìn)行詳細(xì)講解，但是相關(guān)概念依舊相當(dāng)抽象，學(xué)生難以掌握。本文所介紹的數(shù)值方法有助于學(xué)生理解相關(guān)概念。

Matlab是mathworks公司發(fā)布的集數(shù)值計(jì)算、矩陣計(jì)算、數(shù)據(jù)可視化以及系統(tǒng)建模仿真等功能的計(jì)算語(yǔ)言[1]。本文利用Matlab強(qiáng)大的矩陣運(yùn)算和繪圖功能實(shí)現(xiàn)定態(tài)薛定諤方程離散譜和相應(yīng)本征函數(shù)的數(shù)值求解。文中數(shù)值計(jì)算所使用的主要Matlab代碼已在正文中給出。

一、線性諧振子模型的本征值與本征函數(shù)

線性諧振子模型無論在經(jīng)典力學(xué)還是在量子力學(xué)均占有重要地位，其應(yīng)用領(lǐng)域非常廣泛。一維線性諧振子所滿足的定態(tài)薛定諤方程可寫為：，其中和E分別為普朗克常數(shù)、系統(tǒng)哈密頓量、振子的有效質(zhì)量、特征頻率和能量。該系統(tǒng)的本征值和相應(yīng)的本征波函數(shù)為[2]

（1）

和

（2）

其中n=0，1，2…，和。和分別為諧振子能級(jí)升降算符（或產(chǎn)生湮滅算符）。|n>為諧振子本征能量標(biāo)記為n的狀態(tài)（即Fock態(tài)）。需要指出的是態(tài)|n>具有正交歸一性，即（克羅內(nèi)克函數(shù)）。可證明的是坐標(biāo)表象的本征波函數(shù)可構(gòu)成完備的希爾伯特空間。上式中，為厄米多項(xiàng)式，其可顯式寫為：

（3）

其中，函數(shù)表示向下取整，Matlab標(biāo)準(zhǔn)庫(kù)函數(shù)floor可實(shí)現(xiàn)該運(yùn)算。由方程（2）和（3）可得到坐標(biāo)表象下的本征函數(shù)，其可由以下自定義函數(shù)Harmonic實(shí)現(xiàn)：

function psi_n=Harmonic（alpha，x，n）

xt=sqrt（alpha）*x；

N_n=（alpha^2/pi）^0.25/sqrt（2.0^n*factorial（n））；

h = zeros（1，n+1）；n_fact=factorial（n）；

for m=0：floor（n/2）

h（2*m+1）=n_fact*（-1.0）^m/（factorial（m）...

*factorial（n-2.0*m））*2^（n-2.0*m）；

end

Herfun=polyval（h，xt）；

psi_n=N_n*Herfun.*exp（-1.0*xt.^2/2.0）；

需要指出的是，厄米多項(xiàng)式采用方程（3）計(jì)算，而不是Matlab提供的庫(kù)函數(shù)hermiteH。該庫(kù)函數(shù)運(yùn)行十分緩慢，其原因在于該庫(kù)函數(shù)根據(jù)厄米多項(xiàng)式遞推關(guān)系生成高階形式。

二、計(jì)算方法

接下來我們介紹如何利用線性諧振子本征態(tài)函數(shù)所建構(gòu)的希爾伯特正交空間表示坐標(biāo)和動(dòng)量算符。由量子力學(xué)表象理論可知，坐標(biāo)和動(dòng)量算符具有矩陣形式，其矩陣元對(duì)應(yīng)于直接內(nèi)積運(yùn)算，可表示為和（n和m=0，1，2…）。基于類似的矩陣形式，人們可以進(jìn)一步的計(jì)算任意勢(shì)函數(shù)的矩陣元形式，其形如。利用式子，

，和，上述矩陣元的值可表示為

（4）

和

（5）

相應(yīng)的矩陣形式可寫為：

（6）

和

（7）

需要指出的是實(shí)際的數(shù)值計(jì)算過程中，展開的基矢不可能取無窮多個(gè)，必須采取截?cái)?。此類截?cái)鄬⒂绊懹?jì)算精度。如果指標(biāo)n且其最大值為N（截?cái)啵?，那么上述坐?biāo)與動(dòng)量算符將被截?cái)喑蔀榈木仃嚒Ｔ诮酉聛淼挠?jì)算過程中，我們采用自然單位（）。

哈密頓量的矩陣表示都可以由方程（6）和（7）計(jì)算得到，即將坐標(biāo)和動(dòng)量算符直接由其矩陣形式替換，如。實(shí)際數(shù)值計(jì)算過程中，哈密頓量的矩陣表示也將被截?cái)?，本征值及其本征函?shù)計(jì)算將被引入誤差，且計(jì)算結(jié)果的收斂性與截?cái)郚大小有關(guān)[3]。如能得到哈密頓量矩陣的本征值，則該本征值對(duì)應(yīng)該量子系統(tǒng)的能量本征值。如能獲得本征值En所對(duì)應(yīng)的本征矢量（T為轉(zhuǎn)置算符），本征值En所對(duì)應(yīng)的本征函數(shù)可由下述式子計(jì)算：

（8）

Maltab自帶的函數(shù)庫(kù)提供了矩陣的本征值和本征矢量的計(jì)算方案，如[Vec，Val]=eig（H）將直接返回H矩陣的本征矢量（Vec）和本征值（Val）。需要指出的是，Matlab軟件返回的本征值Val為矩陣形式，其對(duì)角元為我們所需的本征值，利用Matlab庫(kù)函數(shù)diag函數(shù)即可提取。另外可以證明的是，本征波函數(shù)也具有正交歸一性，即。

為驗(yàn)證該方法的可行性，我們可將該方法用于計(jì)算線性諧振子的本征值，其哈密頓量記為。帶入位置和動(dòng)量算符的矩陣形式，我們可以得到哈密頓量H的矩陣表示。該系統(tǒng)的能譜及其對(duì)應(yīng)的本征函數(shù)可由以下程序?qū)崿F(xiàn)（計(jì)算最低的60個(gè)本征值）[3]：

N=60；Mv=sqrt（1：N）；hbar=1.0；alpha=1.0；

x=1/alpha*（diag（Mv，-1）+diag（Mv，1））/sqrt（2）；

p=i*alpha*hbar*（diag（Mv，-1）-diag（Mv，1））/sqrt（2）；

H=p^2/2+x^2/2；

[Vec，Val]=eig（H）；

[Eigvalue，index]=sort（diag（Val））；

Eigvalue（1：4）

ind=index（2）；x=-10.0：0.1：10.0；psi=0.0*x；

for m=1：（N+1）

psi=psi+Vec（m，ind）*Harmonic（alpha，x，m-1）；

end

plot（x，psi）

當(dāng)取截?cái)嘀礜=60和完備基參量α=1時(shí)，數(shù)值計(jì)算結(jié)果表明：返回的本征值與理論計(jì)算結(jié)果相吻合；根據(jù)方程（8），對(duì)應(yīng)的本征函數(shù)可由返回的本征矢量計(jì)算，其值與理論計(jì)算結(jié)果相符。需要指出的是，原則上完備基中的參量α可以任意取值。實(shí)際的計(jì)算過程中，我們會(huì)發(fā)現(xiàn)：參量α對(duì)本征值的計(jì)算影響較小[5]，但是其對(duì)本征波函數(shù)的計(jì)算影響較大。較合理的選擇方案是：當(dāng)N足夠大時(shí)，坐標(biāo)算符最高次冪系數(shù)與動(dòng)量平方算符系數(shù)比值與α2接近。另外可使用的判據(jù)：取參量α，使得波函數(shù)滿足歸一化條件。

三、在本征值求解中的應(yīng)用

1.厄米量子系統(tǒng)本征值

首先，我們假定所考察的雙阱勢(shì)束縛的振子所對(duì)應(yīng)的哈密頓量為。為計(jì)算該系統(tǒng)的能譜及其相應(yīng)的本征函數(shù)，我們只需將之前的能譜計(jì)算程序第5行關(guān)于哈密頓量定義的命令代碼修改為：

k=5；H=p^2/2+（x^4-k*x^2）/2；

上述雙阱勢(shì)量子振子束縛態(tài)能譜與相應(yīng)本征波函數(shù)即可由前述代碼計(jì)算得到。圖1給出了最低的四個(gè)能量態(tài)的本征值和相應(yīng)的本征函數(shù)。

我們選取軟核勢(shì)振子作為第二個(gè)數(shù)值計(jì)算案例。軟核勢(shì)在原子物理中作為三維庫(kù)倫勢(shì)場(chǎng)在一維情形的等效勢(shì)使用。假定其所對(duì)應(yīng)的哈密頓量為：，其中g(shù)為軟核勢(shì)參量。為計(jì)算軟核勢(shì)作用下的振子能譜，同樣我們只需修改第2節(jié)程序中有關(guān)哈密頓定義部分的代碼，如下：

g=1；H=p^2/2-（eye（N+1）+g*x^2）^（-0.5）；

通過上述程序可得到該系統(tǒng)的基態(tài)能量為-0.6698。

2.非厄米量子系統(tǒng)本征值

近年來，非厄米系統(tǒng)性質(zhì)受到了廣泛的關(guān)注，尤其在時(shí)空反演（PT）對(duì)稱量子系統(tǒng)被發(fā)現(xiàn)以后[4]。PT對(duì)稱操作分別指：P（x，p）→（-x，-p）；T（x，p，i）→（-x，-p，-i）；在一定條件下，非厄米的PT對(duì)稱系統(tǒng)可具有實(shí)數(shù)本征值。正因?yàn)檫@一重要特性使得該類系統(tǒng)被廣泛研究。該類系統(tǒng)的本征值是否為實(shí)數(shù)是核心問題之一。諧振子模型也可以用于非厄米系統(tǒng)本征值的求解。為驗(yàn)證該算法的有效性，首先我們研究具有如下形式的PT對(duì)稱系統(tǒng)哈密頓量：。該系統(tǒng)的本征值為。求解上述系統(tǒng)的本征值，我們只需將第2節(jié)能譜計(jì)算程序第5行關(guān)于哈密頓量定義的命令代碼修改為：H=p^2/2+（x^2-i*2*x-eye（N+1））/2。 最低的20個(gè)本征值數(shù)值計(jì)算結(jié)果與解析解相對(duì)誤差小于10的-10次方。利用上述方法，我們計(jì)算了如下PT對(duì)稱量子系統(tǒng)的本征值：

（9）

其中β為實(shí)數(shù)。該系統(tǒng)的實(shí)數(shù)本征值隨參量β變化分布情形如圖2所示。值得指出的是：當(dāng)β=2時(shí)對(duì)應(yīng)經(jīng)典諧振子模型；當(dāng)β=4時(shí)相應(yīng)的哈密頓量為，該系統(tǒng)無束縛態(tài)。因而β≈4時(shí)系統(tǒng)本征值計(jì)算方法失效，其原因在于諧振子本征波函數(shù)所構(gòu)成的希爾伯特空間只能表示束縛態(tài)。此時(shí)需要利用其它算法進(jìn)行計(jì)算[5]。該系統(tǒng)的能譜可由以下程序?qū)崿F(xiàn)：

N=250；Mv=sqrt（1：N）；hbar=1.0；alpha=4.0；

x=1/alpha*（diag（Mv，-1）+diag（Mv，1））/sqrt（2）；

p=i*alpha*hbar*（diag（Mv，-1）-diag（Mv，1））/sqrt（2）；

for k=1：150

beta=1+k*0.03；H=p^2/2-（i*x）^beta/2；

[Vec，Val]=eig（H）；

[Eigvalue，index]=sort（diag（Val））；

hold on；

for j=1：N

if （abs（imag（Eigvalue（j）））<0.00001）

hold on；plot（beta，real（Eigvalue（j）），'b.'）；

end

end

end

四、結(jié)論

根據(jù)線性諧振子本征函數(shù)所建構(gòu)的希爾伯特空間矩陣表示，我們得到了量子系統(tǒng)哈密頓量的矩陣形式。采用Matlab軟件標(biāo)準(zhǔn)庫(kù)函數(shù)獲得了不同勢(shì)函數(shù)條件下的離散能譜及其相應(yīng)的本征函數(shù)。相關(guān)數(shù)值方法有利于學(xué)生理解并掌握量子力學(xué)表象理論，有助于提高教學(xué)效果。

參考文獻(xiàn)：

[1]張志涌.精通Matlab6.5版[M].北京：北京航空航天大學(xué)出版社，2003.

[2]周世勛.量子力學(xué)教程（第2版）[M].北京：高等教育出版社，2009：29-32，111-114.

[3]Korsch H J，Glück M.Computing quantum eigenvalues made easy [J].Eur.J.Phys，2002，23（4）：413-426.

[4]Bender C M，Boettcher S.Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry [J].Phys.Rev.Lett.，1998，80（24）：5243.

[5]Wessels G J C.A numerical and analytical investigation into non-Hermitian Hamiltonians [D]，South Africa：University of Stellenbosch，2008.

通訊作者：

謝小濤（1980—），男，湖北京山人，陜西師范大學(xué)物理學(xué)與信息技術(shù)學(xué)院教授，主要從事量子物理與量子光學(xué)的教學(xué)與研究工作。