1 BPM-Q 算法的介绍

在研究 BPM 算法时，研究者关注的是电磁场的 TE 模式和 TM 模式。我们令电磁场的传输方向是 z 方向，并且使用缓变包络近似。那么对应的 TE 模式和 TM 模式方程可以写成[4]

其中，fy ( x , y , z ) = A ( x , y , z ) exp ( - iβz)，函数 f = E，H。参数 n 是传播过程中的折射率，而 ne 则是有效折射率，k 则是波数。我们可以发现这样沿着 z 方向传播的 TE 和 TM 模式满足的方程，有对于 z 的一阶偏微商和对于 x 的二阶偏微商。这样的方程形似量子力学中描述微观粒子运动的薛定谔方程。那么接下来我们就从薛定谔方程出发，来描述如何使用 BPM-Q 算法的。

对于一维微观运动的粒子，放置在势场 U ( x , t ) 中后，满足的运动方程为

其中，波函数 ψ ( x , t ) 用以描述微观粒子的概率幅。上述的方程是与时间位置有关的偏微分方程。为了描述问题的方便，在文章中的讨论里，我们选择自然单位制，也即

。我们将适用于电磁场传输方程求解的 BPM 算法进行稍微的变形，得到适用于求解薛定谔方程的 BPM-Q 算法。图 1 为 BPM-Q 算法的流程图，用于说明如何利用波传播方法求解量子力学的演化问题。

如图 1 所示，第一步，我们输入量子初态并确定求解窗口及步长、确定势垒大小；第二步，我们求解量子初态在 t 方向的演化并对薛定谔方程进行差分处理；第三步，我们添加边界条件并从 t = 0 开始进入循环，写出对应时刻的三对角矩阵，再利用追赶法求解对应时刻的波函数，时间步长增加一步，进入判断条件，如果时间小于给定的定值则再次进入循环；如果时间大于给定的定值则跳出循环，画出波函数随时间演化图形即结果。

2 BPM-Q 算法在量子力学中的应用

在上一节中，我们介绍了 BPM-Q 算法的技术方案。接下来，我们将以量子力学中出现的几个具体的例子，来验证我们方法的可行性。

2.1 一维无限深势阱求解方案

首先，我们研究的是量子力学书中出现的最简单的例子——一维无限深势阱问题。一维无限深势阱是量子力学中的典型问题，势阱中的束缚态波函数可以解析地给出，这是我们所熟知的。当势阱的形式较为复杂时，解析求解将变得困难，这时适宜采用数值方法。在本文中我们利用 BPM-Q 算法对一维无限深势阱中运动粒子的波函数进行了数值计算，并与严格解析解相比较，表明了 BPM-Q 算法的可靠性和正确性。

在该问题中，势场在位置空间的分布满足如下

如图 2 所示，可以看到，xL 粒子运动的禁区，粒子在该区域出现的概率为零；在 0≤x≤L 的区域，粒子波函数呈现驻波形式，这是由边界条件决定的，粒子在 x~x+Δx 区域出现的概率为 |ψ(x)|2Δx。

当势场分布如式(4)所示时，我们是可以将薛定谔方程进行分离变量法解析计算出波函数的。我们令 ( x , t ) = ψ ( x ) f ( t )， 我们可以直接得到波函数和能量的解析表达式[5]

下面我们使用 BPM-Q 算法来数值计算该问题，并和解析结果进行对比。在数值计算中，我们选择的参数如下，n 是正值的整数，

是约化普朗克常数，m 是电子质量。无限深势阱中定态薛定谔方程为

我们把宽度为 L 的区域均分成相距为 Δx 的 N 个小区间，这样就可以把薛定谔方程写成差分方程形式

在讨论中，我们使用的是自然单位制，其中

和 m 都为 1，那么该表达式变成

因此，上式可表示成矩阵形式：Hψn = Eψn。例如，N=5 的二阶导数矩阵可以写成 3×3 矩阵

H 矩阵大小为 ( N - 2 ) × ( N - 2 ) 而不是 N×N，因为函数的二阶导数无法在终点进行计算，即 n=1 和 n=N。

2.2 一维无限深势阱数值模拟

我们利用 BPM-Q 算法对一维无限深势阱中运动粒子的薛定谔方程进行数值求解。如前所述，粒子不可能出现在 xL 区域，即在禁区内有 ψn(x)=0，因此我们只需要在 0≤x≤L 范围内进行求解即可。通过 Matlab 编程，我们利用上述方法数值计算给出的前五个本征函数对应的本征值如表 1 所示，求得的前三个本征函数如图 3(c)所示。

由图 3 可知，计算结果和我们预期的一致，一维无限深势阱中运动粒子的波函数是束缚定态波函数，是符合式(7)的驻波解。随着量子数 n 的增加，波函数的节点数增加，对应的能量本征值增大，呈现出离散的能谱。从表 1 可以看到，数值方法给出的本征值与解析式基本相符。我们可以看出，利用 BPM-Q 算法得到数值结果和解析结果基本一致。从而验证了我们发展的 BPM-Q 算法的可行性。

3 一维方势垒穿透问题

其次，我们讨论一维方势垒的穿透问题。一维方势垒的穿透问题具有很重要的意义，在该问题的求解中，满足薛定谔方程的微观粒子在运动中能够展示出不同于经典物理的现象。即能量低于势垒高度的粒子，也可以穿过势垒,我们称为隧道效应（Tunnel Effect）。它是由于粒子的波动性引起的，我们假设微观粒子以能量 E 从势垒的左边入射，如图 4 所示，方势垒的形式为

考虑了波函数的连续有界等边界条件后，我们可以得到波包透射系数的解析表达式为[5,6]

其中，T0 为透射系数，E0 为粒子能量。在求解该问题时，我们讨论如何利用 BPM-Q 算法求解该问题。

3.1 一维方势垒穿透求解方案

一维方势垒中粒子的波函数所满足的定态薛定谔方程为

我们将方势垒均分为 N 段，每一段的宽度均为 Δx，将薛定谔方程写为差分形式，那么表达式变为

在讨论中，我们使用的是自然单位制，其中

和 m 都为1，那么表达式变成

考虑边界条件，我们可以将薛定谔方程写为矩阵形式

其为 ( N + 1 ) × ( N + 1 ) 阶的矩阵方程，其中

当 2≤n≤N 时，

使用追赶法对该矩阵方程进行求解得到波函数 ψ。透射系数为

，为势垒右侧出现波包的概率[6]。

3.2 一维方势垒穿透数值模拟

假设一维势垒的势能函数为

我们分别使用 BPM-Q 算法和采用严格解析式求解得到此方形势垒的透射系数，结果如图 5 所示，基本一致，验证了 BPM-Q 算法的正确性。

4 含时薛定谔方程(势场与时间无关)

4.1 含时薛定谔方程求解方案

在这里我们将 BPM-Q 算法应用到含时薛定谔方程，在我们的讨论中，我们使用的是自然单位制，其中

和 m 都为 1，那么该表达式变成

接下来我们分别在 x 方向和 t 方向对上式进行离散化，使等号左边包含所有 ψ 的 l+1 项，等号右边包含所 ψ 的 l 项，其差分格式依次可表示为

化简得到

式中

则原方程可化为 Aψ = d，式中

由于矩阵 A 是三对角矩阵，因此

从上式可以解出 ψ，得到了 ψ 的所有值，也就得到了微观粒子的分布随时间的变化情况。

4.2 含时薛定谔方程数值模拟

我们假设微观粒子以能量 E 从势垒的左边入射，方势垒的形式为

接下来我们假设，当 t = 0 时，提供初值 ψ ( x , 0 ) = sin ( 0.36x )，x∈( 0 , 30 )。进行数值模拟后的结果如图 6 所示。

图 6(a)中，我们假设全反射的边界条件，取空间步长 dx = 0.04，时间步长 dt = 0.04，势垒 U = 35，势垒放在 X 坐标的 150~170 的位置，入射光波从 X 轴的 1~30 位置入射，当时间计算到 t = 100 时，我们得到 BPM-Q 算法的不同时刻的概率密度函数随着空间节点的变化关系如图 6(a)所示，其中纵轴表示空间节点，可以看出，当 t = 1 时，图中表示的是初始试探概率密度函数随着时间节点的变化关系，此时微观粒子还没有与方势垒作用;随着时间的推移，在 t = 25 左右，微观粒子与方势垒开始作用；t = 30 时，其作用效果比较明显；在 t = 40 时,可以清楚地看出，在方势垒的下面有一定的波形，说明微观粒子有一定的概率可以穿透方势垒，这正是量子隧道效应的具体体现。同时从量子力学观点来看，由于粒子的波动性,问题与波透过一层介质相似，有一部分波穿过,一部分波被反射回去。故随着粒子与方势垒的作用，图中上面区域概率密度函数的最大峰值减小，波形密集度增加，势垒区和下面区域有一部分波。按照波函数的统计诠释,粒子有一部分概率穿过势垒，有一定的概率被反射回去。

我们改变势垒的大小，当其他条件不变，势垒 U = 50 时，我们，得到 BPM-Q 算法的不同时刻的概率密度函数随着空间节点的变化关系如图 6(b)所示，在 t = 25 时，微观粒子与方势垒开始作用；t = 30 时，其作用效果比较明显；在 t = 40 时，可以清楚地看出，在方势垒的下面波形非常微弱，说明高势垒可以阻挡微观粒子穿透方势垒，从而绝大部分的波被高势垒反射回去。按照波函数的统计诠释,粒子以小概率穿过高势垒，大概率被反射回去。

我们再改变势垒的大小，当其他条件不变，势垒 U = 0 时，我们得到 BPM-Q 算法的不同时刻的概率密度函数随着空间节点的变化关系如图 6(c)所示，可以清楚的看出，在没有方势垒的情况下，粒子按照正常的轨迹在空间中行走。

通过具体计算实例，进一步证实了 BPM-Q 算法的正确性与合理性，为其在量子计算中的应用提供了理论基础。

5 含时薛定谔方程(势场与时间有关)

5.1 含时薛定谔方程求解方案

上一节讨论了当势场不变时，粒子穿过势垒时态的演变情况，而在更一般的情况下，势场也随着时间发生变化，BPM-Q 算法仍可以进行有效的计算。在我们的讨论中，我们使用的是自然单位制，其中

和 m 都为 1，那么微观粒子波函数满足方程

5.2 含时薛定谔方程数值模拟

我们假设微观粒子以能量 E 从势垒的左边入射，方势垒的形式为

我们假设初始波包的形态为：ψ ( x , 0 ) = sin ( 0.36x )，x ∈( 0 , 30 )；进行数值模拟后的结果如图 7 所示。

图 7(a)中，我们假设全反射的边界条件，时间节点 t < 40，所选区域的势垒高度为零，t > 40，势垒高度为 35，当时间计算到 t = 100 时，我们得到 BPM-Q 算法的不同时刻的概率密度函数随着时间节点的变化关系如图 7(a) 所示，其中纵轴表示空间节点，可以看出，当时间 t < 40 时，概率密度函数的变化情况与势垒为零时基本一致，微观粒子在区域内传播与反射不受干扰；随着时间的推移，在时间 t = 60 左右，在边界处反射回来的波包碰到了，会受到方势垒的影响；t = 70 时，其作用效果比较明显。此时可以清楚地看出，在 x ≥ 170 的位置有一定的波形，而在 x ≤ 150 的位置处也出现了一定的波形，表明微观粒子有一定概率透射势垒，这也是量子隧道效应的体现。

当我们继续扩大势垒高度至 U = 60 时，概率密度函数的变化情况如图 7(b) 所示。t < 40 时势垒高度为零，波包运动不受影响。t > 40 时，势垒高度 U = 60，在波包与势垒进行作用后，在 x ≤ 150 的位置处几乎没有出现波形，表明波包穿透势垒的概率变得微乎其微。

6 结语

在本论文中，我们将适用于数值求解电磁场 TE 模式和 TM 模式的 BPM 算法拓展到量子力学的研究中，并发展出适用于描述薛定谔方程的 BPM-Q 算法。我们讨论了三种不同的势场下，微观粒子的运动行为。对于不含时的一维无限深势阱和一维方势垒中，我们使用 BPM-Q 算法得到的数值结果和教材上的理论计算表达式进行对比，我们发现能够吻合的非常好。对于含时的一维势垒问题，我们也使用 BPM-Q 算法给出了波包在该势场下的运动规律。BPM-Q 算法的优点包括形式简单、便于处理非均匀结构、和具有天然的并行特点等。BPM-Q 算法可以在量子力学领域中起到重要作用。BPM-Q 算法对今后进一步开展相关的创新工作有着重要的借鉴和启迪意义。

参考文献

[1]张三慧. 大学物理学[M]. 3版. 北京: 清华大学出版社, 2018.

[2]赵凯华, 罗蔚茵. 新概念物理教程 量子物理[M]. 2版. 北京: 高等教育出版社, 2015.

[3]胡海云, 缪劲松, 冯艳全, 等. 大学物理. 第四卷, 近代物理[M]. 北京: 高等教育出版社, 2017.

[4]熊群. 有限差分光束传输法的研究[D]. 成都: 电子科技大学, 2008.

[5]曾谨言. 量子力学卷 1[M]. 5 版. 北京: 科学出版社, 2013.

[6]李明. 有限差分方法求任意形状势垒的透射系数[J]. 科技信息, 2013(21): 161.

LI M. Finite difference method for calculating transmission coefficient of arbitrary shape barrier[J]. Science and Technology Information, 2013(21): 161.

基金项目：国家自然科学基金（11974046）；北京理工大学教改项目（经典光学波传播方法展示量子力学现象）。

作者简介：陈天，男，北京理工大学副教授，主要从事物理教学科研工作，研究方向为拓扑量子计算，[email protected]。

引文格式: 张灿, 周卫康, 陈天. 波传播算法在量子力学中的应用及其数值模拟[J]. 物理与工程, 2022, 32(5): 150-156.

Cite this article: ZHANG C, ZHOU W K, CHEN T. Application of Beam Propagation Method in quantum mechanics and its numerical simulation[J]. Physics and Engineering, 2022, 32(5): 150-156. (in Chinese)

END返回搜狐，查看更多