氢原子球坐标薛定谔方程数值解

             

预备知识 薛定谔方程,类氢原子的定态波函数,Wigner 3j 符号

   本文使用原子单位制.我们希望求解氢原子在电场中的薛定谔方程

\begin{equation} -\frac{1}{2m} \boldsymbol{\nabla}^2 \Psi + V( \boldsymbol{\mathbf{r}} , t)\Psi = \mathrm{i} \frac{\partial}{\partial{t}} \Psi \end{equation}
虽然最直观的方法是使用直角坐标,但计算效率很低.实际中一般使用球坐标系,用球谐函数展开波函数.如果哈密顿算符是轴对称的(例如线偏振电场),那么只需要 $m = 0$ 的球谐函数,即勒让德多项式.
\begin{equation} \Psi( \boldsymbol{\mathbf{r}} , t) = \frac{1}{r}\sum_{l,m} \psi_{l,m}(r) Y_{l,m}( \boldsymbol{\mathbf{r}} ) \end{equation}
其中 $\psi_{l,m}(r)$ 是约化径向波函数(scaled radial wave function),式中每一项叫做一个分波(partial wave)

   在薛定谔方程中,一种简单的做法是使用以下势能表示库仑势能以及电磁场对电子的作用

\begin{equation} V( \boldsymbol{\mathbf{r}} , t) = -\frac{Z}{r} - q \boldsymbol{\mathbf{E}} (t) \boldsymbol\cdot \boldsymbol{\mathbf{r}} \end{equation}
在原子单位中电子电荷 $q = -1$,该势能使用了长度规范,此外我们还可以使用速度规范速度规范,也是等效的.在长度规范中,我们近似认为原子的长度远小于电磁波的波长,所以电磁波不随位置变化.

1. 线偏振光

   若我们取电场极化方向为 $ \hat{\boldsymbol{\mathbf{z}}} $,则角动量 $L_z$ 是一个守恒量.假设初始波函数关于 $ \hat{\boldsymbol{\mathbf{z}}} $ 轴对称,那么在波函数的整个演化(propagation)过程中,我们只需要 $m=0$ 的球谐函数展开波函数,即

\begin{equation} \Psi( \boldsymbol{\mathbf{r}} , t) = \frac{1}{r}\sum_{l'} \psi_{l'}(r) Y_{l', 0}( \hat{\boldsymbol{\mathbf{r}}} ) \end{equation}
另外薛定谔方程中 $V( \boldsymbol{\mathbf{r}} , t) = \boldsymbol{\mathbf{E}} (t) \boldsymbol\cdot \boldsymbol{\mathbf{r}} = E(t) z$,进而可以用球谐函数表示(式 20
\begin{equation} V( \boldsymbol{\mathbf{r}} , t) = \sqrt{\frac{4\pi}{3}} E(t) r \cdot Y_{1,0}( \hat{\boldsymbol{\mathbf{r}}} ) \end{equation}
以上两式代入薛定谔方程,再把每一项与 $Y_{l,0}( \hat{\boldsymbol{\mathbf{r}}} )$ 做内积(放入积分 $\int Y_{1,0}^*( \hat{\boldsymbol{\mathbf{r}}} ) \,\mathrm{d}{\Omega} $ 中)可得一系列耦合的方程
\begin{equation} H_0 \psi_{l}(r) + \sum_{l' = 0}^{\infty} E(t)rF_{l, l'} \psi_{l'}(r) = \mathrm{i} \frac{\partial \psi_{l}}{\partial t} \quad (l=0,1,\dots) \end{equation}
其中无场哈密顿算符为
\begin{equation} H_0 = -\frac{1}{2m} \frac{\mathrm{d}^{2}}{\mathrm{d}{r}^{2}} -\frac{Z}{r} + \frac{l(l+1)}{2mr^2} \end{equation}
矩阵 $ \boldsymbol{\mathbf{F}} $ 就是跃迁偶极子矩阵的角向积分(式 17
\begin{equation} \begin{aligned} F_{l,l'} &= \sqrt{\frac{4\pi}{3}} \left\langle Y_{l,0} \middle| Y_{1,0} \middle| Y_{l',0} \right\rangle \\ &= \sqrt{(2l+1)(2l'+1)} \begin{pmatrix}l & 1 & l'\\ 0 & 0 & 0\end{pmatrix} ^2\\ &= \sqrt{\frac{2l+1}{2l'+1}} \begin{bmatrix}l & 1 & l'\\ 0 & 0 & 0\end{bmatrix} ^2 \end{aligned} \end{equation}
可见,当没有外场的时候每一个项(即每一个分波)都可以独立演化,而电场将不同的分波耦合起来.根据氢原子的选择定则 $\Delta l = \pm 1$(式 6 ),矩阵 $ \boldsymbol{\mathbf{F}} $ 中除了两条副对角线上的元都为零.另外由式 2 易得 $ \boldsymbol{\mathbf{F}} $ 是一个对称矩阵.

2. 任意偏振光

   将所有 $(l,m)$ 按某种顺序排列,例如

\begin{equation} (0, 0),\ (1,-1),\ (1,0),\ (1,1),\ (2,-2), \dots \end{equation}
并将他们编号为 $\lambda = 1,2, \dots$,那么可以把式 2 记为
\begin{equation} \Psi( \boldsymbol{\mathbf{r}} , t) = \frac{1}{r}\sum_\lambda \psi_\lambda(r) Y_\lambda( \boldsymbol{\mathbf{r}} ) \end{equation}
势能函数为 $V( \boldsymbol{\mathbf{r}} , t) = \boldsymbol{\mathbf{E}} (t) \boldsymbol\cdot \boldsymbol{\mathbf{r}} = E_x x + E_y y + E_z z$.其中 $x,y,z$ 可以用球谐函数表示为(式 20
\begin{equation} x = \sqrt{\frac{2\pi}{3}} r (Y_{1,-1} - Y_{1,1}) \qquad y = \mathrm{i} \sqrt{\frac{2\pi}{3}} r (Y_{1,-1}+Y_{1,1}) \qquad z = \sqrt{\frac{4\pi}{3}} rY_{1,0} \end{equation}

   式 6 的耦合方程拓展为

\begin{equation} \begin{aligned} H_0 \psi_{\lambda}(r) + \sum_{\lambda'} r \left[E_x(t) F_{\lambda, \lambda'}^{x} + E_y(t) F_{\lambda, \lambda'}^{y} + E_z(t) F_{\lambda, \lambda'}^{z} \right] \psi_{\lambda'}(r) \\ = \mathrm{i} \frac{\partial \psi_{\lambda}}{\partial t} \quad (\lambda=0,1,\dots) \end{aligned} \end{equation}
三个耦合矩阵分别为
\begin{equation} \begin{aligned} F_{\lambda,\lambda'}^{x} = \sqrt{\frac{2\pi}{3}} E_x \left( \left\langle Y_{l,m} \middle| Y_{1,-1} \middle| Y_{l',m'} \right\rangle - \left\langle Y_{l,m} \middle| Y_{1,1} \middle| Y_{l',m'} \right\rangle \right) \end{aligned} \end{equation}
\begin{equation} \begin{aligned} F_{\lambda,\lambda'}^{y} = \mathrm{i} \sqrt{\frac{2\pi}{3}} E_y \left( \left\langle Y_{l,m} \middle| Y_{1,-1} \middle| Y_{l',m'} \right\rangle + \left\langle Y_{l,m} \middle| Y_{1,1} \middle| Y_{l',m'} \right\rangle \right) \end{aligned} \end{equation}
\begin{equation} \begin{aligned} F_{\lambda,\lambda'}^{z} = \sqrt{\frac{4\pi}{3}} E_z \left\langle Y_{l,m} \middle| Y_{1,0} \middle| Y_{l',m'} \right\rangle \end{aligned} \end{equation}
其中(式 17
\begin{equation} \left\langle Y_{l,m} \middle| Y_{1,m_1} \middle| Y_{l',m'} \right\rangle = (-1)^m\sqrt{\frac{3(2l+1)(2l'+1)}{4\pi}} \begin{pmatrix}l & 1 & l'\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l & 1 & l'\\ -m & m_1 & m'\end{pmatrix} \end{equation}

3. 任意势能

   如果需要加上任意势能项 $V'( \boldsymbol{\mathbf{r}} , t)$,首先需要进行分解

\begin{equation} V'( \boldsymbol{\mathbf{r}} , t) = \sum_{l,m} V'_{l,m}(r, t) Y_{l,m}( \hat{\boldsymbol{\mathbf{r}}} ) \end{equation}
那么 coupling 矩阵元为(式 17
\begin{equation} \left\langle Y_{l'',m''} \middle| V'( \boldsymbol{\mathbf{r}} , t) \middle| Y_{l',m'} \right\rangle = \sum_{l,m} V'_{l,m}(r, t) \left\langle Y_{l'',m''} \middle| Y_{l,m} \middle| Y_{l',m'} \right\rangle \end{equation}
\begin{equation} \left\langle Y_{l'',m''} \middle| Y_{l,m} \middle| Y_{l',m'} \right\rangle = (-1)^{m''}\sqrt{\frac{(2l''+1)(2l+1)(2l'+1)}{4\pi}} \begin{pmatrix}l''& l& l'\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l'' & l & l'\\ -m'' & m & m'\end{pmatrix} \end{equation}
注意矩阵元取决于 $r$.在程序中,可以把 $ \left\langle Y_{l'',m''} \middle| Y_{l,m} \middle| Y_{l',m'} \right\rangle $ 表示为三维数组 YYY(λ'', λ', λ),然后在每个 $r$ 格点对 λ 加权求和.

   对 $m = 0$ 的对称情况,$ \left\langle Y_{l'',m''} \middle| Y_{l,m} \middle| Y_{l',m'} \right\rangle $ 在 $l = 0$ 时是一个对角矩阵,$l = 1$ 时只有两个 1-副对角线不为零,$l = 2$ 时只有对角线和两个 2-副对角线不为零,以此类推.左上角的三角形也会等于零,见图 2

   在算符拆分时,若把 $V'$ 矩阵对角线上的元合并到 $H_0$ 中很可能会减小误差.

4. 算符拆分

   在实际的程序中,我们可以把演化子 $ \exp\left(- \mathrm{i} H \Delta t\right) $ 拆成 3 项.虽然这么做会引入一定的误差,但是却大大提高了效率

\begin{equation} \exp\left(- \mathrm{i} H \Delta t\right) = \exp\left(- \mathrm{i} H_0 \frac{\Delta t}{2}\right) \exp\left(- \mathrm{i} F \Delta t\right) \exp\left(- \mathrm{i} H_0 \frac{\Delta t}{2}\right) + \mathcal{O}\left(\Delta t^3 \right) \end{equation}
也就是说,在每个时间步长 $\Delta t$ 中,我们可以把波函数先根据
\begin{equation} H_0 \psi_{l}(r) = \mathrm{i} \frac{\partial \psi_{l}}{\partial t} \end{equation}
演化 $\Delta t/2$,其中 $t_{mid}$ 取这段时间的中点.再对每个 $ \boldsymbol{\mathbf{r}} $ 根据
\begin{equation} \sum_{l'} E(t_{mid})rF_{l, l'} \psi_{l'}(r) = \mathrm{i} \frac{\partial \psi_{l}}{\partial t} \end{equation}
演化 $\Delta t$.最后再根据式 21 演化 $\Delta t/2$.

   理论上也可以把演化子划分成两项 $ \exp\left(- \mathrm{i} F \Delta t\right) \exp\left(- \mathrm{i} H_0 {\Delta t}\right) $ 但据说式 20 更精确.至于相邻两步之间产生的 $ \exp\left(- \mathrm{i} H_0 {\Delta t}/{2}\right) \exp\left(- \mathrm{i} H_0 {\Delta t}/{2}\right) $ 是否可以合并为 $ \exp\left(- \mathrm{i} H_0 {\Delta t}\right) $,取决于所使用的算法这么做以后是否会引入更大误差.

5. 网格和演化算法

预备知识 Crank-Nicolson 算法(一维)

  

未完成:以下内容应该放在一维薛定谔方程里面讲解

   可以使用二维数组储存波函数,每一列(或行)是一个分波的 $\psi_l(r)$.径向网格可以使用等间距网格,但 FEDVR 网格效率要更高.

   演化可以并使用 Crank-Nicolson 算法 Crank-Nicolson 算法(一维)演化.但是 Lanczos 算法效率更高.

致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择会员制,大量广告,内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 10 元,我们一个星期内就能脱离亏损, 并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

广告位

投放详情

         

© 小时科技 保留一切权利