用Python入门量子力学

求解薛定谔方程

无限深方势阱是量子力学中非常经典的一个模型，其特点是，在势阱内势能为0，而势阱外，势能为无穷大。对应其薛定谔方程及边界条件为

$$\begin{array}{rl}& \left[-\frac{{\hslash }^2}{2m}\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}+V(x)\right]\psi (x)=E\psi (x)\\ & V(x)=\{\begin{array}{rl}0,& x\in [0,L]\\ \infty ,& \mathrm{others}\end{array}\end{array}$$

其中$\psi (x)$为波函数，$E$为能量本征值，若只考虑$V(x)=0$时的情况，则波函数的通解为

$$\psi \left(x\right)=C_1{e}^{-\frac{\sqrt{2}\sqrt{m}x\sqrt{-E}}{\hslash }}+C_2{e}^{\frac{\sqrt{2}\sqrt{m}x\sqrt{-E}}{\hslash }}$$

代码如下

from sympy import symbols, Function, Eq, diff, dsolve
from sympy import print_latexx = symbols('x')  # 空间坐标
L = symbols('L', positive=True)  # 势阱宽度
E = symbols('E')  # 能量
m = symbols('m', positive=True)  # 质量
hbar = symbols('hbar', positive=True)  # 普朗克常数
psi = Function('psi')  # 波函数
V = symbols('V')V = 0  # 无限深势阱内势能为 0
s = Eq(-hbar**2 / (2 * m) * diff(psi(x), x, 2) + V * psi(x), E * psi(x))res = dsolve(s, psi(x))
print_latex(res)

下面考虑边界条件，令$\psi (0)=0,\psi (L)=0$，则薛定谔方程的波函数和本征能量分别为

$${\psi }_n(x)=C_1\sin \left(\frac{\pi nx}{L}\right),E_n=\frac{{\pi }^2{\hslash }^2{n}^2}{2L^{2}m}$$

由于波函数${\psi }_n(x)$受到本征值的影响，sympy在计算时，会默认选择一个最简单的数值，若不事先设置$E_n$，则最终得到的$\psi (x)$将为0，故而其代码为

from sympy import pin = symbols('n', integer=True, positive=True)
hbar = symbols('hbar', positive=True)
E = n**2 * pi**2 * hbar**2 / (2 * m * L**2)schrodinger_eq = Eq(-hbar**2 / (2 * m) * psi(x).diff(x, 2), E * psi(x))res = dsolve(schrodinger_eq, psi(x), ics={psi(0): 0, psi(L): 0})
print_latex(res)

谐振子模型

一维谐振子用经典的眼光来看就是可以沿着一个方向振动的弹簧，其薛定谔方程为

$$-\frac{{\hslash }^2}{2m}\frac{{\mathrm{d}}^2\psi (x)}{\mathrm{d}{x}^2}+\frac{1}{2}m{\omega }^2{x}^2\psi (x)=E\psi (x)$$

其哈密顿算符可表示为

$$\hat{H}=-\frac{{\hslash }^2}{2m}\frac{{\text{d}}^2}{\text{d}{x}^2}+\frac{1}{2}m{\omega }^2{x}^2$$

其本征态的能级可表示为

$$E_n=(n+\frac{1}{2})\hslash \omega$$

当$n=0$时，仍有$\frac{1}{2}\hslash \omega$，体现了量子力学与经典力学的差别。

由于谐振子模型的通解中存在埃尔米特多项式，故而无法简单地通过dsolve推导出来，但sympy中提供了【qho_1d】模块，封装了谐振子的一些特性

from sympy import print_latex
from sympy.abc import m, x, omega
from sympy.physics.qho_1d import E_n, psi_n
print_latex(E_n(x, omega))
print_latex(psi_n(0, x, m, omega))

其中，【E_n】函数即为一维谐振子能量。【psi_n】则返回波函数${\psi }_n$，其输入四个参数分别是$n$，坐标$x$，质量$m$以及角频率$\omega$，示例代码的输出结果为

$$\frac{\sqrt[4]{m\omega }{e}^{-\frac{m\omega x^2}{2\hslash }}}{\sqrt[4]{\hslash }\sqrt[4]{\pi }}$$

sympy中海提供了三维谐振子模块【sho】，其中【E_nl】函数可返回谐振子能量，【R_nl】函数可返回谐振子波函数，这两个函数的前两个参数均为$n,l$表示主量子数和轨道量子数，示例如下

from sympy.physics.sho import R_nl, E_nl
from sympy.abc import r, nu, l
print_latex(E_nl(0, 0, 1))
print_latex(R_nl(0, 0, 1, r))

氢原子波函数

作为量子力学第一个成功的应用，氢原子的数学本质，是薛定谔方程在中心立场中的解，薛定谔方程和中心力场的势能$V(r)$表示如下

$$-\frac{{\hslash }^2}{2\mu }{\nabla }^2\psi +V(r)\psi =E\psi ,V(r)=-\frac{e^2}{r\pi {ϵ}_{0}r}$$

使用球坐标，将Laplace算子展开

− ℏ 2 2 μ r 2 { ∂ ∂ r ( r 2 ∂ ∂ r ) + 1 sin ⁡ 2 θ [ sin ⁡ θ ∂ ∂ θ ( sin ⁡ θ ∂ ∂ θ ) + ∂ 2 ∂ ϕ 2 ] } ψ − e 2 r π ϵ 0 r = E ψ -\frac{\hbar^2}{2\mu r^2}\{\frac{\partial}{\partial r}(r^2\frac{\partial}{\partial r})+\frac{1}{\sin^2\theta}[\sin\theta\frac{\partial}{\partial\theta}(\sin\theta\frac{\partial}{\partial\theta})+\frac{\partial^2}{\partial\phi^2}]\}\psi-\frac{e^2}{r\pi\epsilon_0r}=E\psi −2μr2ℏ2​{∂r∂​(r2∂r∂​)+sin2θ1​[sinθ∂θ∂​(sinθ∂θ∂​)+∂ϕ2∂2​]}ψ−rπϵ0​re2​=Eψ

其解为径向函数和球谐函数的乘积

$$\psi (r,\theta ,\varphi )=R_{nl}(r)Y_{lm}(\theta ,\varphi )$$

$Y_{lm}$是球谐函数，其定义为

$$Y_n^m(\theta ,\varphi )=\sqrt{\frac{2n+1}{4\pi }\frac{(n-m)!}{(n+m)!}}{e}^{im\theta }{P}_n^m(\cos \varphi )$$

其中$P_n^m$是连带勒让德多项式，定义为

$$\begin{gathered} P_n^m=(-1{)}^m(1-x^2{)}^{m/2}\frac{{\text{d}}^m}{\text{d}{x}^m}{P}_n(x) \\ {P}_n(x)=\sum _{k=0}^{\infty }\frac{(-n{)}_k(n+1{)}_k}{(k!{)}^2}(\frac{1+x}{2}{)}^k \\ {n}_k=\{\begin{array}{rlr}& 1& k=0\\ & n(n+1)\cdots (n+k-1)& k>1\end{array} \end{gathered}$$

$R_{nl}$可以继续展开为

$$R_{nl}=\sqrt{(\frac{2}{na}{)}^3\frac{(n-l-1)!}{2n[(n+l)!{]}^3}}{e}^{-\frac{r}{na}}(\frac{2r}{na}{)}^l{L}_{n-l-1}^{2l+1}(\frac{2r}{na})$$

在上面两个公式中，主要出现了两个特殊函数，在scipy中均有实现。其中$Y_{lm}$是球谐函数，对应【sph_harm】函数；$L_{n-l-1}^{2l+1}$是拉盖尔多项式，对应【assoc_laguerre】函数。

刨除这两个令人费解的符号之后，波函数表达式中就只剩下一些数值变量了，其中$a$是玻尔半径0.53$\text{A˚}$，在写代码是带入0.53，表示整个模型以埃米为单位。

$n,l,m$则为不同本征函数的本征值，一般在物理学中被称为量子数

• $n$ 主量子数，为正整数，表示电子能级或其所处壳层。
• $l$ 角量子数，取值范围是$0\sim n-1$的整数，决定电子云的形状
• $m$ 磁量子数，取值范围是$-l\sim l$，代表电子在特定亚壳层中的具体轨道。

【hydrogen】模块提供了与氢原子相关的函数

• E_nl(n) 返回Hatree能量
• E_nl_dirac(n, l) 相对论修正后的Hatree能量。
• R_nl(n, l, r) 氢原子的径向波函数
• Psi_nlm(n, l, m, r, phi, theta) 氢原子波函数$|n,l,m⟩$

这三个函数都有一个默认参数Z=1，表示原子个数；r, phi, theta为空间坐标。

E_nl_dirac另有两个变量，spin_up为自旋，默认为True,表示$⟨\uparrow |$，False即为$|\downarrow ⟩$。c为原子单位制下的光速，默认为精细结构常数的倒数$137.036$。原子物理中常用Hatree作为能量单位，$1E_h=4.3597\times 10^{-18}J$，约等于基态氢原子的电子势能。

通过下面的代码，可以打印出不同能级下电子的波函数

from sympy.physics.hydrogen import Psi_nlmr, phi, theta = symbols("r phi theta", positive=True)
Z=Symbol("Z", positive=True, integer=True, nonzero=True)print("|$nlm$|$\|\psi_{nlm}\\big>$|")
print("|:-|:-|")
for n in range(6):for l in range(n):for m in range(-l,l+1):md = latex(Psi_nlm(n,l,m,r,phi,theta,Z))print(f"|$\\psi_{{{n},{l},{m}}}$|${md}$", )

部分常见波函数的电子云如图所示。