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

贡献者： addis

1. 算符拆分

预备知识 1　氢原子的含时薛定谔方程（球坐标），算符的指数函数、波函数传播子

　　本文使用原子单位制。在实际的程序中，我们可以把演化子 $ \exp\left(- \mathrm{i} H \Delta t\right) $ 拆成 3 项。虽然这么做会引入一定的误差（$H_0$ 和 $ \boldsymbol{\mathbf{\mathcal E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} $ 不对易），但是却大大提高了效率

\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} H_1 \Delta t\right) \exp\left(- \mathrm{i} H_0 \frac{\Delta t}{2}\right) + \mathcal{O}\left(\Delta t^3 \right) ~. \end{equation}

　　例如对于线偏振光（式 7 ），在每个时间步长 $\Delta t$ 中，我们可以把波函数先根据方程

\begin{equation} H_0 \psi_{l} = \mathrm{i} \frac{\partial \psi_{l}}{\partial t} ~. \end{equation}

演化 $\Delta t/2$，其中 $t_{mid}$ 取这段时间的中点。再对每个 $r$ 根据方程

\begin{equation} \sum_{l'} \mathcal E(t_{mid})rF_{l, l'} \psi_{l'} = \mathrm{i} \frac{\partial \psi_{l}}{\partial t} ~, \end{equation}

演化 $\Delta t$。最后再根据式 2 演化 $\Delta t/2$。具体演化算法有多种，将在下面介绍。

　　至于相邻两步之间产生的 $ \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) $，取决于所使用的算法这么做以后是否会引入额外误差（例如 Crank-Nicolson 算法就不宜这么做）。

2. 网格和演化算法

预备知识 2　Crank-Nicolson 算法（一维）

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

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

　　演化可以并使用 Crank-Nicolson 算法Crank-Nicolson 算法（一维）演化。但是 Lanczos 算法效率更高，而且可以实时判断误差改变步长。

　　拆分后的每个算符（矩阵）演化的算法可以一样或不一样。

3. 电场演化的直接计算

预备知识 3　平面波的球谐展开

　　事实上，注意到 $ \exp\left( \mathrm{i} q \boldsymbol{\mathbf{\mathcal E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \,\mathrm{d}{t} \right) $ 不过是一个普通的平面波函数而不是微分算符，所以我们只需要把它和波函数相乘：$ \exp\left( \mathrm{i} q \boldsymbol{\mathbf{\mathcal E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \,\mathrm{d}{t} \right) \Psi$。为了使相乘后的函数仍然具有式 4 的形式，可以先根据式 1 对其进行分波展开（$ \exp\left( \mathrm{i} q \boldsymbol{\mathbf{\mathcal E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \,\mathrm{d}{t} \right) = \sum_{l'',m''}R_{l'',m''}Y_{l'',m''}$）。所以演化后的每个分波的径向波函数就是

\begin{equation} \begin{aligned} \psi_{l,m} &= r \left\langle Y_{l,m} \middle| \exp\left( \mathrm{i} q \boldsymbol{\mathbf{\mathcal E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \,\mathrm{d}{t} \right) \Psi \right\rangle = \left\langle Y_{l,m} \middle| \sum_{l'',m''} R_{l'',m''} Y_{l'',m''} \sum_{l',m'} \psi_{l',m'} Y_{l',m'} \right\rangle \\ &= \sum_{l',m'} \left[\sum_{l'',m''}R_{l'',m''} \left\langle Y_{l,m} \middle| Y_{l'',m''} \middle| Y_{l',m'} \right\rangle \right] \psi_{l',m'}~. \end{aligned} \end{equation}

现在和上文一样，可以用式 19 把三个球谐函数的积分变为两个 3j 符号的乘积，再根据选择定则，排除两个求和中等于零的项。另外如果电场只沿 $z$ 方向，那么式中所有 $m$ 为零。

　　式 4 相当于对每个 $r$ 处不同分波的波函数进行一个矩阵乘法，但每个 $r$ 处的矩阵是不同的。

　　一种看似可能的近似方法是把 $ \exp\left( \mathrm{i} q \boldsymbol{\mathbf{\mathcal E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \,\mathrm{d}{t} \right) $ 展开为前两三项。但这样一来其实和直接求 $ \boldsymbol{\mathbf{F}} $ 矩阵进而 $ \boldsymbol{\mathbf{F}} ^2$ 没有什么区别了。事实上保留两三项并不稳定，不然也不需要用 Crank-Nicolson 或者 Lanczos 这么费时的办法了。所以还是老老实实把平面波展开成贝塞尔函数。

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