贡献者： addis
预备知识 类氢原子的定态波函数
，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}
在薛定谔方程中，一种简单的做法是使用以下势能表示库仑势能以及电磁场对电子的作用
\begin{equation}
V( \boldsymbol{\mathbf{r}} , t) = -\frac{Z}{r} - q \boldsymbol{\mathbf{E}} (t) \boldsymbol\cdot \boldsymbol{\mathbf{r}} ~.
\end{equation}
在原子单位中电子电荷 $q = -1$，该势能使用了长度规范
，此外我们还可以使用速度规范速度规范
，也是等效的。在长度规范中，我们近似认为原子的长度远小于电磁波的波长，所以电磁波不随位置变化。本文使用长度规范。
虽然求解方程最直观的方法是使用直角坐标，但计算效率很低。实际中一般使用球坐标系，用球谐函数展开波函数（参考球坐标系中的定态薛定谔方程）。
\begin{equation}
\Psi( \boldsymbol{\mathbf{r}} , t) = \frac{1}{r}\sum_{l,m} \psi_{l,m}(r, t) Y_{l,m}( \boldsymbol{\mathbf{r}} )~.
\end{equation}
其中 $\psi_{l,m}(r, t)$ 是
约化径向波函数（scaled radial wave function），式中每一项叫做一个
分波（partial wave）。如果哈密顿算符是关于 $z$ 轴对称的（例如线偏振电场），且初始波函数也轴对称，那么波函数将始终保持轴对称。这时只需要 $m = 0$ 的球谐函数，即勒让德多项式（见
式 3 ）。
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, t) 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$，进而可以用球谐函数表示（
式 26 ）
\begin{equation}
V( \boldsymbol{\mathbf{r}} , t) = E(t) r \cos\theta = \sqrt{\frac{4\pi}{3}} E(t) r \cdot Y_{1,0}( \hat{\boldsymbol{\mathbf{r}}} )~.
\end{equation}
以上两式代入
式 1 ，再把每一项与 $Y_{l,0}( \hat{\boldsymbol{\mathbf{r}}} )$ 做内积（放入积分 $\int Y_{1,0}^*( \hat{\boldsymbol{\mathbf{r}}} ) \,\mathrm{d}{\Omega} $ 中）可得一系列耦合的方程
\begin{equation}
H_0 \psi_{l} + E(t)r\sum_{l' = 0}^{\infty} F_{l, l'} \psi_{l'} = \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}} $ 就是跃迁偶极子矩阵的角向积分（
式 19 ）
\begin{equation}
\begin{aligned}
F_{l,l'} &= \left\langle Y_{l,0} \middle| \cos\theta \middle| Y_{l',0} \right\rangle = \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$，那么可以把
式 3 记为
\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$ 可以用球谐函数表示为（
式 26 ）
\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}
的耦合方程拓展为
\begin{equation}
\begin{aligned}
H_0 \psi_{\lambda}(r) + r \sum_{\lambda'} \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}} \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}} \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} = \left\langle Y_{l,m} \middle| \cos\theta \middle| Y_{l',m'} \right\rangle
= \sqrt{\frac{4\pi}{3}} \left\langle Y_{l,m} \middle| Y_{1,0} \middle| Y_{l',m'} \right\rangle ~.
\end{aligned}
\end{equation}
其中（
式 19 ）
\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}
这在 “氢原子的跃迁偶极子
” 也有出现。当 $m_1=0$ 时，使用（
式 20 ）令
式 15 为
\begin{equation}
\mathcal C_{l,m} = \left\langle Y_{l,m} \middle| \cos\theta \middle| Y_{l+1,m} \right\rangle = \sqrt{\frac{(l+1)^2-m^2}{(2l+1)(2l+3)}}~.
\end{equation}
\begin{equation}
\boldsymbol{\mathbf{F}} ^{z}
= \begin{pmatrix}
0 & \mathcal C_{0m} & 0 & 0 & \dots\\
\mathcal C_{0m} & 0 & \mathcal C_{1m} & 0 &\dots\\
0 & \mathcal C_{1m} & 0 & \mathcal C_{2m} &\dots\\
\vdots & 0 & \mathcal C_{2m} & 0 & \dots\end{pmatrix} ~.
\end{equation}
3. 任意含时势能
如果要给式 2 加上一个额外的势能项 $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}
那么耦合矩阵元为（
式 19 ）
\begin{equation}
F'_{\lambda,\lambda'}(r, t) = \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}
在程序中，可以把 $ \left\langle Y_{l'',m''} \middle| Y_{l,m} \middle| Y_{l',m'} \right\rangle $ 表示为三维数组
F1(λ'', λ', λ)，然后在每个 $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 ）。耦合薛定谔方程变为
\begin{equation}
H_0 \psi_{l} + E(t)r\sum_{l' = 0}^{\infty} F_{l, l'} \psi_{l'} + \sum_{l'=0}^\infty F'_{l,l'}(r, t)\psi_{l'} = \mathrm{i} \frac{\partial \psi_{l}}{\partial t} \quad (l=0,1,\dots)~.
\end{equation}
在下面介绍的算符拆分中，若把 $F'$ 矩阵对角线上的元合并到 $H_0$ 中很可能会减小误差。
4. 算符拆分
在实际的程序中，我们可以把演化子 $ \exp\left(- \mathrm{i} H \Delta t\right) $ 拆成 3 项。虽然这么做会引入一定的误差（$H_0$ 和 $ \boldsymbol{\mathbf{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}
例如对于线偏振光（式 6 ），在每个时间步长 $\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'} E(t_{mid})rF_{l, l'} \psi_{l'} = \mathrm{i} \frac{\partial \psi_{l}}{\partial t} ~,
\end{equation}
演化 $\Delta t$。最后再根据
式 24 演化 $\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 算法 就不宜这么做）。
5. 网格和演化算法
预备知识 Crank-Nicolson 算法（一维）
未完成：以下内容应该放在一维薛定谔方程里面讲解
可以使用二维数组储存波函数，每一列（或行）是一个分波的 $\psi_l(r)$。径向网格可以使用等间距网格，但 FEDVR 网格效率要更高。
演化可以并使用 Crank-Nicolson 算法 Crank-Nicolson 算法（一维）演化。但是 Lanczos 算法效率更高，而且可以实时判断误差改变步长。
拆分后的每个算符（矩阵）演化的算法可以一样或不一样。
6. 电场演化的直接计算
事实上，注意到 $ \exp\left( \mathrm{i} q \boldsymbol{\mathbf{E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \,\mathrm{d}{t} \right) $ 不过是一个普通的平面波函数而不是微分算符，所以我们只需要把它和波函数相乘：$ \exp\left( \mathrm{i} q \boldsymbol{\mathbf{E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \,\mathrm{d}{t} \right) \Psi$。为了使相乘后的函数仍然具有式 3 的形式，可以先根据式 1 对其进行分波展开（$ \exp\left( \mathrm{i} q \boldsymbol{\mathbf{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{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$ 为零。
式 26 相当于对每个 $r$ 处不同分波的波函数进行一个矩阵乘法，但每个 $r$ 处的矩阵是不同的。
一种看似可能的近似方法是把 $ \exp\left( \mathrm{i} q \boldsymbol{\mathbf{E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \,\mathrm{d}{t} \right) $ 展开为前两三项。但这样一来其实和直接求 $ \boldsymbol{\mathbf{F}} $ 矩阵进而 $ \boldsymbol{\mathbf{F}} ^2$ 没有什么区别了。事实上保留两三项并不稳定，不然也不需要用 Crank-Nicolson 或者 Lanczos 这么费时的办法了。所以还是老老实实把平面波展开成贝塞尔函数。
7. 速度规范
要使用速度规范（注意仍然是位置表象而不是动量表象），$H_0$ 算符的计算是一样的，唯一不同的是把 $ \boldsymbol{\mathbf{E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} $ 换成了 $ \boldsymbol{\mathbf{A}} \boldsymbol\cdot \boldsymbol{\mathbf{p}} $。式 14 。
先把电场限制在 $z$ 方向，所以场的作用主要就是
\begin{equation}
\frac{\partial}{\partial{z}} = \cos\theta \frac{\partial}{\partial{r}} - \frac{\sin\theta}{r} \frac{\partial}{\partial{\theta}} ~.
\end{equation}
第二项只耦合不同的分波。但第一项要更为复杂，它耦合不同分波中同一有限元中的不同基底。所以 $ \exp\left(- \mathrm{i} \lambda \frac{\partial}{\partial{z}} \right) $ 需要把 $ \frac{\partial}{\partial{z}} $ 作用在整个波函数上面，然后用 lanczos 这样的整体方法来演化。
在 FEDVR 中，$ \partial/\partial r $ 可以用矩阵 $ \boldsymbol{\mathbf{D}} $ 精确表达，$\cos\theta$ 也在上文中可以表示为分波耦合矩阵。所以 $ \partial/\partial z $ 就是把这两个矩阵相乘即可。
要注意第一项的角向并不是 $ \left\langle Y_{l,m} \middle| \cos\theta \middle| Y_{l',m} \right\rangle $，而是要同时考虑径向
\begin{equation}
\begin{aligned}
&\quad r \left\langle Y_{l,m} \middle| \cos\theta \frac{\partial}{\partial{r}} \middle| \frac{\psi_{l',m}}{r}Y_{l',m} \right\rangle \\
&= \frac{\mathrm{d}}{\mathrm{d}{r}} \psi_{l',m} \left\langle Y_{l,m} \middle| \cos\theta \middle| Y_{l',m} \right\rangle
-\frac{\psi_{l',m}}{r} \left\langle Y_{l,m} \middle| \cos\theta \middle| Y_{l',m} \right\rangle ~.
\end{aligned}
\end{equation}
另外由
式 21 得
\begin{equation}
F^{(v2z)}_{l,l'} = - \left\langle Y_{l,m} \middle| \cos\theta + \sin\theta \frac{\partial}{\partial{\theta}} \middle| Y_{l',m} \right\rangle ~.
\end{equation}
\begin{equation}
\boldsymbol{\mathbf{F}} ^{(v2)}
= \begin{pmatrix}
0 & \mathcal C_{0m} & 0 & 0 & \dots\\
-\mathcal C_{0m} & 0 & 2\mathcal C_{1m} & 0 &\dots\\
0 & -2\mathcal C_{1m} & 0 & 3\mathcal C_{2m} &\dots\\
\vdots & 0 & -3\mathcal C_{2m} & 0 & \dots\end{pmatrix} ~.
\end{equation}
\begin{equation}
\left\langle Y_{l,m} \middle| \frac{\partial}{\partial{z}} \middle| Y_{l',m} \right\rangle = F_{l,l'}^z \frac{\partial}{\partial{r}} + \frac{F_{l,l'}^{(v2z)}}{r}~.
\end{equation}
所以式 27 的第一项实际上需要拆分成两项，没有导数的那个合并到第二项中去。
图 1：[Hprop]
致读者： 小时百科一直以来坚持所有内容免费，这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此，我们请求广大读者
热心打赏 ，使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 10 元，我们一个星期内就能脱离亏损， 并保证在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款， 他们的付出帮助了 99% 的读者免费获取知识， 我们在此表示感谢。
