谢小涛 董元浩 王建康

【摘要】束缚量子系统的本征值精确可解模型稀少，如何利用数值方法得到系统的本征值是初等量子力学教学过程中不得不面对的问题之一。本文从线性谐振子模型出发，借助谐振子本征值态所建构的希尔伯特空间，利用Matlab提供的库函数求解任意束缚势量子系统的本征值和本征函数。数值方法的讲授将极大提高学生灵活应用量子力学知识解决实际问题的能力。

【关键词】谐振子 束缚态 薛定谔方程 Matlab

【中图分类号】O413 【文献标识码】A 【文章编号】2095-3089（2018）03-0048-02

引言

薛定谔方程束缚态离散谱的计算是初等量子力学所研究的基本问题之一。对于定态薛定谔方程本征值求解问题，由于可解的模型稀少，数值方法成为必不可少的解决方案。数值计算在量子力学教学方面如何取得突破是值得探讨的问题之一，其能有效帮助学生形成更加清晰的物理图像并能大幅提高学生初步掌握并解决实际问题的能力。本文不关注上述数值方法，而侧重于如何基于线性谐振子模型本征波函数所建构的正交完備基求解系统的本征值及其相应的本征函数。正交完备基所组成的希尔伯特相关性质在初等量子力学教学过程均会进行详细讲解，但是相关概念依旧相当抽象，学生难以掌握。本文所介绍的数值方法有助于学生理解相关概念。

Matlab是mathworks公司发布的集数值计算、矩阵计算、数据可视化以及系统建模仿真等功能的计算语言[1]。本文利用Matlab强大的矩阵运算和绘图功能实现定态薛定谔方程离散谱和相应本征函数的数值求解。文中数值计算所使用的主要Matlab代码已在正文中给出。

一、线性谐振子模型的本征值与本征函数

线性谐振子模型无论在经典力学还是在量子力学均占有重要地位，其应用领域非常广泛。一维线性谐振子所满足的定态薛定谔方程可写为：，其中和E分别为普朗克常数、系统哈密顿量、振子的有效质量、特征频率和能量。该系统的本征值和相应的本征波函数为[2]

（1）

和

（2）

其中n=0，1，2…，和。和分别为谐振子能级升降算符（或产生湮灭算符）。|n>为谐振子本征能量标记为n的状态（即Fock态）。需要指出的是态|n>具有正交归一性，即（克罗内克函数）。可证明的是坐标表象的本征波函数可构成完备的希尔伯特空间。上式中，为厄米多项式，其可显式写为：

（3）

其中，函数表示向下取整，Matlab标准库函数floor可实现该运算。由方程（2）和（3）可得到坐标表象下的本征函数，其可由以下自定义函数Harmonic实现：

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提供的库函数hermiteH。该库函数运行十分缓慢，其原因在于该库函数根据厄米多项式递推关系生成高阶形式。

二、计算方法

接下来我们介绍如何利用线性谐振子本征态函数所建构的希尔伯特正交空间表示坐标和动量算符。由量子力学表象理论可知，坐标和动量算符具有矩阵形式，其矩阵元对应于直接内积运算，可表示为和（n和m=0，1，2…）。基于类似的矩阵形式，人们可以进一步的计算任意势函数的矩阵元形式，其形如。利用式子，

，和，上述矩阵元的值可表示为

（4）

和

（5）

相应的矩阵形式可写为：

（6）

和

（7）

需要指出的是实际的数值计算过程中，展开的基矢不可能取无穷多个，必须采取截断。此类截断将影响计算精度。如果指标n且其最大值为N（截断），那么上述坐标与动量算符将被截断成为的矩阵。在接下来的计算过程中，我们采用自然单位（）。

哈密顿量的矩阵表示都可以由方程（6）和（7）计算得到，即将坐标和动量算符直接由其矩阵形式替换，如。实际数值计算过程中，哈密顿量的矩阵表示也将被截断，本征值及其本征函数计算将被引入误差，且计算结果的收敛性与截断N大小有关[3]。如能得到哈密顿量矩阵的本征值，则该本征值对应该量子系统的能量本征值。如能获得本征值En所对应的本征矢量（T为转置算符），本征值En所对应的本征函数可由下述式子计算：

（8）

Maltab自带的函数库提供了矩阵的本征值和本征矢量的计算方案，如[Vec，Val]=eig（H）将直接返回H矩阵的本征矢量（Vec）和本征值（Val）。需要指出的是，Matlab软件返回的本征值Val为矩阵形式，其对角元为我们所需的本征值，利用Matlab库函数diag函数即可提取。另外可以证明的是，本征波函数也具有正交归一性，即。

为验证该方法的可行性，我们可将该方法用于计算线性谐振子的本征值，其哈密顿量记为。带入位置和动量算符的矩阵形式，我们可以得到哈密顿量H的矩阵表示。该系统的能谱及其对应的本征函数可由以下程序实现（计算最低的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时，数值计算结果表明：返回的本征值与理论计算结果相吻合；根据方程（8），对应的本征函数可由返回的本征矢量计算，其值与理论计算结果相符。需要指出的是，原则上完备基中的参量α可以任意取值。实际的计算过程中，我们会发现：参量α对本征值的计算影响较小[5]，但是其对本征波函数的计算影响较大。较合理的选择方案是：当N足够大时，坐标算符最高次幂系数与动量平方算符系数比值与α2接近。另外可使用的判据：取参量α，使得波函数满足归一化条件。

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

1.厄米量子系统本征值

首先，我们假定所考察的双阱势束缚的振子所对应的哈密顿量为。为计算该系统的能谱及其相应的本征函数，我们只需将之前的能谱计算程序第5行关于哈密顿量定义的命令代码修改为：

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

上述双阱势量子振子束缚态能谱与相应本征波函数即可由前述代码计算得到。图1给出了最低的四个能量态的本征值和相应的本征函数。

我们选取软核势振子作为第二个数值计算案例。软核势在原子物理中作为三维库伦势场在一维情形的等效势使用。假定其所对应的哈密顿量为：，其中g为软核势参量。为计算软核势作用下的振子能谱，同样我们只需修改第2节程序中有关哈密顿定义部分的代码，如下：

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

通过上述程序可得到该系统的基态能量为-0.6698。

2.非厄米量子系统本征值

近年来，非厄米系统性质受到了广泛的关注，尤其在时空反演（PT）对称量子系统被发现以后[4]。PT对称操作分别指：P（x，p）→（-x，-p）；T（x，p，i）→（-x，-p，-i）；在一定条件下，非厄米的PT对称系统可具有实数本征值。正因为这一重要特性使得该类系统被广泛研究。该类系统的本征值是否为实数是核心问题之一。谐振子模型也可以用于非厄米系统本征值的求解。为验证该算法的有效性，首先我们研究具有如下形式的PT对称系统哈密顿量：。该系统的本征值为。求解上述系统的本征值，我们只需将第2节能谱计算程序第5行关于哈密顿量定义的命令代码修改为：H=p^2/2+（x^2-i*2*x-eye（N+1））/2。 最低的20个本征值数值计算结果与解析解相对误差小于10的-10次方。利用上述方法，我们计算了如下PT对称量子系统的本征值：

（9）

其中β为实数。该系统的实数本征值随参量β变化分布情形如图2所示。值得指出的是：当β=2时对应经典谐振子模型；当β=4时相应的哈密顿量为，该系统无束缚态。因而β≈4时系统本征值计算方法失效，其原因在于谐振子本征波函数所构成的希尔伯特空间只能表示束缚态。此时需要利用其它算法进行计算[5]。该系统的能谱可由以下程序实现：

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

四、結论

根据线性谐振子本征函数所建构的希尔伯特空间矩阵表示，我们得到了量子系统哈密顿量的矩阵形式。采用Matlab软件标准库函数获得了不同势函数条件下的离散能谱及其相应的本征函数。相关数值方法有利于学生理解并掌握量子力学表象理论，有助于提高教学效果。

参考文献：

[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—），男，湖北京山人，陕西师范大学物理学与信息技术学院教授，主要从事量子物理与量子光学的教学与研究工作。