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

贡献者： 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 这么费时的办法了。所以还是老老实实把平面波展开成贝塞尔函数。