諧振子模型在量子系統(tǒng)本征值計算中的應用

謝小濤 董元浩 王建康

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

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

【中圖分類號】O413 【文獻標識碼】A 【文章編號】2095-3089（2018）03-0048-02

引言

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

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

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

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

（1）

和

（2）

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

（3）

其中，函數(shù)表示向下取整，Matlab標準庫函數(shù)floor可實現(xiàn)該運算。由方程（2）和（3）可得到坐標表象下的本征函數(shù)，其可由以下自定義函數(shù)Harmonic實現(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）；

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

二、計算方法

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

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

（4）

和

（5）

相應的矩陣形式可寫為：

（6）

和

（7）

需要指出的是實際的數(shù)值計算過程中，展開的基矢不可能取無窮多個，必須采取截斷。此類截斷將影響計算精度。如果指標n且其最大值為N（截斷），那么上述坐標與動量算符將被截斷成為的矩陣。在接下來的計算過程中，我們采用自然單位（）。

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

（8）

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

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

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

三、在本征值求解中的應用

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

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

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

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

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

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

（9）

其中β為實數(shù)。該系統(tǒng)的實數(shù)本征值隨參量β變化分布情形如圖2所示。值得指出的是：當β=2時對應經(jīng)典諧振子模型；當β=4時相應的哈密頓量為，該系統(tǒng)無束縛態(tài)。因而β≈4時系統(tǒng)本征值計算方法失效，其原因在于諧振子本征波函數(shù)所構成的希爾伯特空間只能表示束縛態(tài)。此時需要利用其它算法進行計算[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

四、結論

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

參考文獻：

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

[2]周世勛.量子力學教程（第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—），男，湖北京山人，陜西師范大學物理學與信息技術學院教授，主要從事量子物理與量子光學的教學與研究工作。