贡献者： addis
预备知识 广义球谐函数
，Lanczos 算法
，
氢原子球坐标薛定谔方程数值解
本文使用原子单位制。
1. 波函数展开
使用角向基底为两个电子总轨道角动量和 $z$ 分量的本征态，即广义球谐函数 $\mathcal{Y}_{l_1,l_2}^{L,M}( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2)$。
总波函数所在的空间可以看做角向空间和径向空间的张量积空间。现在有了角向空间的正交归一基底（内积由式 12 定义为 $\int \dots \,\mathrm{d}{\Omega_1} \,\mathrm{d}{\Omega_2} $），总波函数就可以在该基底上展开（$ \left\lvert R_\lambda \right\rangle $ 并不是基底也非正交归一，只相当于展开系数）
\begin{equation}
\left\lvert \Psi \right\rangle = \sum_\lambda \left\lvert R_\lambda \right\rangle \otimes \left\lvert \mathcal{Y_\lambda} \right\rangle ~.
\end{equation}
其中 $\lambda$ 是将所有的 $(l_1,l_2,L,M)$ 组合排序后的序号。$ \left\lvert R_\lambda \right\rangle $ 是二维径向波函数
\begin{equation}
\left\lvert R_\lambda \right\rangle = \frac{1}{r_1 r_2} \psi_\lambda(r_1, r_2)~.
\end{equation}
即总波函数为
\begin{equation}
\Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \sum_{L, M, l_1, l_2} \frac{1}{r_1 r_2} \psi_{l_1, l_2}^{L, M}(r_1, r_2)\mathcal{Y}_{l_1, l_2}^{L, M}( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2)~.
\end{equation}
由于总波函数的内积定义为 $\int \dots r_1^2 \,\mathrm{d}{r_1} \,\mathrm{d}{\Omega_1} r_2^2 \,\mathrm{d}{r_2} \,\mathrm{d}{\Omega_2} $，所以径向空间的内积定义为 $\int \dots r_1^2 \,\mathrm{d}{r_1} r_2^2 \,\mathrm{d}{r_2} $。
由张量积空间中的结论，我们可以求哈密顿算符 $H$ 关于角向基底的 “矩阵元” $H_{\lambda, \lambda'}$，每个矩阵元是径向空间中的一个算符。
\begin{equation}
H_{\lambda, \lambda'} = \left\langle \mathcal{Y}_{\lambda} \middle| H \middle| \mathcal{Y}_{\lambda'} \right\rangle ~.
\end{equation}
列出薛定谔方程的 “矩阵形式”，就得到了一组耦合的二维径向波函数的薛定谔方程
\begin{equation}
\sum_{\lambda'} H_{\lambda, \lambda'} \left\lvert R_{\lambda'} \right\rangle = \mathrm{i} \frac{\partial}{\partial{t}} \left\lvert R_{\lambda} \right\rangle ~,
\end{equation}
总哈密顿可以表示为
\begin{equation}
H = H_1 + H_2 + V_{12} + V_{F1} + V_{F2}~.
\end{equation}
其中 $H_i \ \ (i = 1, 2)$ 是单个电子的哈密顿算符（
式 11 ），$V_{12}$ 是两电子之间的库仑势能，$V_{Fi}$ 是电场对每个电子的作用势能。
$H_i$ 不耦合分波，对应对角矩阵
\begin{equation}
\begin{aligned}
\left\langle \mathcal{Y}_{\lambda} \middle| H_i \middle| \mathcal{Y}_{\lambda'} \right\rangle &= \delta_{\lambda,\lambda'} \left(K_{ri} + \frac{L_i^2}{2 r_i^2} - \frac{2}{r_i} \right) \\
&= \delta_{\lambda,\lambda'} \left(-\frac{1}{2m} \frac{\partial^{2}}{\partial{r_i}^{2}} + \frac{l_i^2}{2m r_i^2} - \frac{2}{r_i} \right) \qquad (i = 1,2)~,
\end{aligned}
\end{equation}
其中第二项 $L_i$ 变为 $l_i$ 是因为 $ \left\lvert \mathcal Y_\lambda \right\rangle $ 基底是 $L_i^2$ 的本征函数。剩下的哈密顿项 $V_{12} + V_{F1} + V_{F2}$ 都会耦合分波。
2. 波函数的对称性
我们一般假设电子自旋态恒为 singlet（反对称），所以总波函数满足粒子交换对称 $\Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \Psi( \boldsymbol{\mathbf{r}} _2, \boldsymbol{\mathbf{r}} _1)$。由式 7 可以导出分波的对称性
\begin{equation}
\psi_{l_1, l_2}^{L, M}(r_1, r_2) = (-1)^{l_1 + l_2 - L} \psi_{l_2, l_1}^{L, M}(r_2, r_1)~.
\end{equation}
理论上用该式可以在程序中减少一半重复的分波，并且分波的 coupling 演化也可以得到精简，但由于编程麻烦暂时没有实现。
另外可以证明（其实是数值计算发现的）$l_1 + l_2 - L$ 为奇数的分波（奇分波）和 $l_1 + l_2 - L$ 为偶数的分波（偶分波）在线偏振电场中是独立演化的，也就是 $V_{12}$ 和 $z_1, z_2$ 不会把奇偶性不同的分波耦合起来。所以如果氦原子从基态（只包含偶分波：$l_1 = l_2, L = 0$）开始在线偏振电场中演化，那么奇数分波的项将始终为零，可以在程序中去除。
3. 电子相互作用项
\begin{equation}
\begin{aligned}
V_{12} = \frac{1}{ \left\lvert \boldsymbol{\mathbf{r}} _2 - \boldsymbol{\mathbf{r}} _1 \right\rvert } &= 4\pi \sum_{l=0}^{\infty} \frac{1}{2l+1} \frac{r_<^l}{r_>^{l+1}} \sum_{m = -l}^l Y_{l,m}^*( \hat{\boldsymbol{\mathbf{r}}} _1) Y_{l,m}( \hat{\boldsymbol{\mathbf{r}}} _2)\\
&= 4\pi \sum_{l=0}^{\infty} \frac{(-1)^l}{\sqrt{2l + 1}} \frac{r_<^l}{r_>^{l+1}} \mathcal{Y}_{l,l}^{0,0} ( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2)~.
\end{aligned} \end{equation}
这就相当于将张量积空间中的线性算符分解成径向和角向算符的张量积之和。第二步使用了
式 9 ，即
\begin{equation}
\mathcal{Y}_{l,l}^{0,0} ( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2) = \frac{(-1)^l}{\sqrt{2l+1}} \sum_{m = -l}^l Y_{lm}^*( \hat{\boldsymbol{\mathbf{r}}} _1) Y_{lm}( \hat{\boldsymbol{\mathbf{r}}} _2)~.
\end{equation}
要计算矩阵元 $ \left\langle \mathcal Y_\lambda \middle| V_{12} \middle| \mathcal Y_{\lambda'} \right\rangle $，就要对 6 个球谐函数的乘积的线性组合做两次角向积分。使用式 12 可以将上式表示为六个 CG 系数相乘的线性组合（二重求和）。
\begin{equation}
\begin{aligned}
&\quad \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| Y_{lm}^*( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{lm}( \hat{\boldsymbol{\mathbf{r}}} _2) \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle \\
&= \sum_{m_1',m_2'}\sum_{m_1,m_2} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ m_1' & m_2' & M'\end{bmatrix} \times\\
&\qquad \iint Y_{l_1m_1}^*( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l_2m_2}^*( \hat{\boldsymbol{\mathbf{r}}} _2)Y_{l m}^*( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l m}( \hat{\boldsymbol{\mathbf{r}}} _2)Y_{l_1' m_1'}( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l_2' m_2'}( \hat{\boldsymbol{\mathbf{r}}} _2) \,\mathrm{d}{\Omega_1} \,\mathrm{d}{\Omega_2} \\
&=\sum_{m_1',m_2'}\sum_{m_1,m_2} (-1)^{m_1 + m_2 + m} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ m_1' & m_2' & M'\end{bmatrix} \times\\
&\qquad \int Y_{l_1, -m_1}( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l,-m}( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l_1' m_1'}( \hat{\boldsymbol{\mathbf{r}}} _1) \,\mathrm{d}{\Omega_1} \int Y_{l_2, -m_2}( \hat{\boldsymbol{\mathbf{r}}} _2)Y_{l m}( \hat{\boldsymbol{\mathbf{r}}} _2)Y_{l_2' m_2'}( \hat{\boldsymbol{\mathbf{r}}} _2) \,\mathrm{d}{\Omega_2} \\
&=\frac{2l+1}{4\pi}\sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)} \sum_{m_1',m_2'}\sum_{m_1,m_2} (-1)^{m_1 + m_2 + m} \times\\
&\qquad \begin{pmatrix}l_1 & l & l_1'\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l_1 & l & l_1'\\ -m_1 & -m & m_1'\end{pmatrix} \begin{pmatrix}l_2 & l & l_2'\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l_2 & l & l_2'\\ -m_2 & m & m_2'\end{pmatrix} \times\\
&\qquad \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ m_1' & m_2' & M'\end{bmatrix} ~.
\end{aligned} \end{equation}
现在试图对 $m$ 求和并使用选择定则（3j 符号第二行之和为 0），得（注意 $m$ 的范围）该求和中最多只有一项不为 0，所以
\begin{equation}
\begin{aligned}
& \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| \mathcal{Y}_{l,l}^{0,0} \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle
=\\
& \delta_{M,M'} (-1)^l \frac{\sqrt{2l+1}}{4\pi} \sqrt{(2l_1+1)(2l_1'+1)(2l_2+1)(2l_2'+1)}\times\\
&\sum_{m_1',m_2'}\sum_{m_1,m_2} (-1)^{m_1' + m_2} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ m_1' & m_2' & M\end{bmatrix} \begin{pmatrix}l_1 & l & l_1'\\ 0 & 0 & 0\end{pmatrix} \times\\
& \begin{pmatrix}l_1 & l & l_1'\\ -m_1 & m_1-m_1' & m_1'\end{pmatrix} \begin{pmatrix}l_2 & l & l_2'\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l_2 & l & l_2'\\ -m_2 & m_2-m_2' & m_2'\end{pmatrix} ~.
\end{aligned} \end{equation}
其中对 $m_1, m_2$ 的求和要求
\begin{equation}
m_1 + m_2 = m_1' + m_2' = M~, \qquad
\left\lvert m_1'-m_1 \right\rvert \leqslant l~, \qquad
\left\lvert m_2'-m_2 \right\rvert \leqslant l~.
\end{equation}
由于 CG 系数都是实数，这必定是一个实数矩阵。选择定则如下
- 根据总角动量守恒（$1/r_{12}$ 算符与 $L, L_z$ 算符对易），$L \ne L'$ 或 $M \ne M'$ 时，矩阵元为零。
- 根据 parity CG coefficients，当 $l_1' + l + l_1$ 或 $l_2' + l + l_2$ 为奇数时，上式等于 0。
- 根据 $\Pi \mathcal{Y} = (-1)^{l_1 + l_2} \mathcal{Y}$（式 5 ），考虑到 $1/r_{12}$ 具有偶宇称，所以 $l_1 + l_2$ 与 $l_1' + l_2'$ 的奇偶必须不同。
爱华的论文中将式 12 用 9j 符号表示，但仅限于 $M \ne M'$。然而经过数值验证，我发现应该可以将爱华的公式拓展为
\begin{equation}
\begin{aligned}
\left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| \mathcal{Y}_{ll}^{00} \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle
= &\frac{2l+1}{4\pi} \delta_{L,L'}\delta_{M,M'} \sqrt{(2l_1'+1)(2l_2'+1)(2L+1)}\times\\
& \begin{bmatrix}l & l_1' & l_1\\ 0 & 0 & 0\end{bmatrix}
\begin{bmatrix}l & l_2' & l_2\\ 0 & 0 & 0\end{bmatrix}
\left\{\begin{matrix}l & l_1' & l_1\\ l & l_2' & l_2\\ 0 & L & L\end{matrix}\right\} ~.
\end{aligned} \end{equation}
注意当 $M = M'$ 时，结果与 $M$ 无关。
泊松积分
按照 FEDVR 基底的思路，我们会自然地想把 $V_{12}$ 算符展开后的 $ \left\langle r_i r_j \middle| r_<^l/r_>^l \middle| r_{i'} r_{j'} \right\rangle $ 用高斯 Lobatto 积分来计算，然而这么做并不精确，因为 $r_<^l/r_>^l$ 用有限阶多项式展开会有误差。使用泊松积分法可以减小计算误差，且仍然保持对角矩阵的形式，即
\begin{equation}
\langle{r_i r_j}|{\frac{r_<^l}{r_>^{l+1}}}|{r_{i'} r_{j'}}\rangle \propto \delta_{i,i'} \delta_{j,j'}~.
\end{equation}
根据我的猜测（有空验证一下！），泊松积分积分的结果应该和直接用数值积分的结果是一样的。
\begin{equation}
\langle{r_i r_j}|{\frac{r_<^l}{r_>^{l+1}}}|{r_{i'} r_{j'}}\rangle = \delta_{i,i'} \delta_{j,j'} \left[\frac{2l+1}{r_i r_j \sqrt{\omega_i \omega_j}} (2 \boldsymbol{\mathbf{T}} _l)^{-1}_{i,j} + \frac{r_i^l r_j^l}{r_{max}^{2l+1}} \right] ~.
\end{equation}
推导参考 McCurdy, Solving the three-body Coulomb breakup problem using exterior complex scaling 式 68。其中 $ \boldsymbol{\mathbf{T}} _l$ 是单电子动能算符的矩阵
\begin{equation}
\begin{aligned}
(T_l)_{ij} &= \langle{r_i}|{ \left[-\frac{1}{2m} \frac{\mathrm{d}^{2}}{\mathrm{d}{r}^{2}} + \frac{l(l+1)}{2mr} \right] }|{r_j}\rangle \\
&= -\frac{1}{2m} \langle{r_i}|{ \frac{\mathrm{d}^{2}}{\mathrm{d}{r}^{2}} }|{r_j}\rangle + \delta_{i,j}\frac{l(l+1)}{2mr_i^2}~.
\end{aligned}
\end{equation}
其中 $ \left\lvert r_i \right\rangle $ 是一维 FEDVR 格点在 $r_i$ 处的
正交归一基底，$r_{max}$ 是计算 $T_l$ 时最后一个 FE 的边界。虽然
式 16 中含有 $r_{max}$，但结果却与 $r_{max}$ 无关，原因是 $T_l$ 也和 $r_{max}$ 相关，两项中的 $r_{max}$ 抵消了。
4. 电场作用项（长度规范）
长度规范中，对每个电子，
\begin{equation}
V_F = \boldsymbol{\mathbf{\mathcal E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} = {\mathcal E}_x x + {\mathcal E}_y y + {\mathcal E}_z z~
\end{equation}
与氢原子的情况相同，有
\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 =2\sqrt{\frac{\pi}{3}} rY_{1,0}~.
\end{equation}
对线偏振电场我们只需要最后一项，以第一个电子为例，令
\begin{equation}
F^{(z1)} = \cos\theta_1~.
\end{equation}
\begin{equation} \begin{aligned}
& \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| \cos\theta_1 \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle
= 2\sqrt{\frac{\pi}{3}} \sum_{m_1',m_2'}\sum_{m_1,m_2} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2 & L'\\ m_1' & m_2' & M'\end{bmatrix} \times\\
&\qquad \iint Y_{l_1m_1}^*( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l_2m_2}^*( \hat{\boldsymbol{\mathbf{r}}} _2)Y_{10}( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l_1' m_1'}( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l_2' m_2'}( \hat{\boldsymbol{\mathbf{r}}} _2) \,\mathrm{d}{\Omega_1} \,\mathrm{d}{\Omega_2} \\
&= 2 \sqrt{\frac{\pi}{3}} \delta_{l_2, l_2'} \sum_{m_1,m_2, m_1'} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ m_1' & m_2 & M'\end{bmatrix} \times\\
&\qquad \int Y_{l_1m_1}^*( \hat{\boldsymbol{\mathbf{r}}} _1) Y_{10}( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l_1' m_1'}( \hat{\boldsymbol{\mathbf{r}}} _1) \,\mathrm{d}{\Omega_1} \\
&= \delta_{l_2, l_2'} \sqrt{(2l_1 + 1)(2l_1'+1)} \sum_{m_1,m_2, m_1'} (-1)^{m_1} \times\\
&\qquad \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2 & L'\\ m_1' & m_2 & M'\end{bmatrix} \begin{pmatrix}l_1 & 1 & l_1'\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l_1 & 1 & l_1'\\ -m_1 & 0 & m_1'\end{pmatrix} \\
&= \delta_{l_2, l_2'} \delta_{M, M'} \sqrt{(2l_1 + 1)(2l_1'+1)} \sum_{m_1,m_2} (-1)^{m_1} \times\\
&\qquad \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2 & L'\\ m_1 & m_2 & M\end{bmatrix} \begin{pmatrix}l_1 & 1 & l_1'\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l_1 & 1 & l_1'\\ -m_1 & 0 & m_1\end{pmatrix} ~.
\end{aligned} \end{equation}
最后一步使用了 3j 符号中 $m_1 + m_2 + m_3 = 0$，以及 CG 系数中 $m_1 + m_2 = M$ 的性质。两个 delta 函数验证了对第一个电子作用的 $z$ 方向电场不会影响第二个电子的角动量以及总角动量的 $z$ 分量。推导中还陆续得出了 $m_1' = m_1, m_2' = m_2$ 的结论，即该电场同样不会影响每个电子的角动量 $z$ 分量（所以对 $m_1', m_2'$ 的求和消去了）。最后对 $m_1,m_2$ 的求和中，要求 $m_1+m_2=M$，且 $ \left\lvert m_1 \right\rvert < \min \left\{l_1, l'_1 \right\} $。
同理，对第二个电子有
\begin{equation} \begin{aligned}
& \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| z_2 \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle
= \delta_{l_1, l_1'} \delta_{M, M'} \sqrt{(2l_2 + 1)(2l_2'+1)} r \sum_{m_1,m_2} (-1)^{m_2} \times\\
&\qquad \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1 & l_2' & L'\\ m_1 & m_2 & M\end{bmatrix} \begin{pmatrix}l_2 & 1 & l_2'\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l_2 & 1 & l_2'\\ -m_2 & 0 & m_2\end{pmatrix} ~.
\end{aligned} \end{equation}
在爱华的毕业论文中，$M = M' = 0$ 时用 9j 符号表示为
\begin{equation}
\begin{aligned} & \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| z_i \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle =
r\sqrt{9(2l'_1+1)(2l'_2+1)(2L'+1)/(4\pi)^2}\\
& \quad\times \begin{bmatrix}\delta_{i,1} & l'_1 & l_1\\ 0 & 0 & 0\end{bmatrix}
\begin{bmatrix}\delta_{i,2} & l'_2 & l_2\\ 0 & 0 & 0\end{bmatrix}
\begin{bmatrix}1 & L' & L\\ 0 & 0 & 0\end{bmatrix}
\left\{\begin{matrix} \delta_{i,1} & l'_1 & l_1\\ \delta_{i,2} & l'_2 & l_2\\ 1 & L' & L\end{matrix} \right\} \quad (i=1,2)~,
\end{aligned}
\end{equation}
：$L' - L = \pm 1$，$M' = M$，$l_1+l_2-L$ 奇偶性不变。这就是为什么从基态（$l_1=l_2$，$L=0$）的电离计算可以排除所有 $l_1+l_2-L$ 为奇数的分波。
但由于单光子只有 $l=1$，所以 $l_1+l_2-L$ 奇偶性不变的条件还可以进一步变为 $l_1,l_2$ 其中一个改变 $\pm 1$。
5. 电场作用项（速度规范）
特别注意：在速度规范下即使只考虑从基态的单电离，也需要好几个 $L$，因为电矢势不为零时，波函数比起长度规范叠乘了一个平面波，而这个平面波需要更多分波才能展开。
类比氢原子（式 25 ），每两个分波之间的耦合从一个数变为一个算符
\begin{equation}
F^{(v1)}_{\lambda,\lambda'} = F^{(z1)}_{\lambda,\lambda'} \frac{\partial}{\partial{r_1}} + \frac{F^{(vz1)}_{\lambda,\lambda'}}{r_1}~.
\end{equation}
其中 $F^{(z1)}_{\lambda,\lambda'}$ 在
式 20 中定义。第二项
\begin{equation}
F^{(vz1)} = -\cos\theta_1 -\sin\theta_1 \frac{\partial}{\partial{\theta_1}} ~.
\end{equation}
为了方便暂且假设所有 $M=0$，有
\begin{equation} \begin{aligned}F^{(vz1)}_{\lambda,\lambda'} &= \left\langle \mathcal{Y}_{l_1 l_2}^{L} \middle| F^{(vz1)} \middle| \mathcal{Y}_{l_1' l_2'}^{L'} \right\rangle \\
&= \delta_{l_2, l_2'} \sum_{m_1,m_2, m_1'} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & 0\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ m_1' & m_2 & 0\end{bmatrix} \times\\
& \qquad\int Y_{l_1m_1}^*( \hat{\boldsymbol{\mathbf{r}}} _1) F^{(vz1)}( \hat{\boldsymbol{\mathbf{r}}} _1)Y_{l_1' m_1'}( \hat{\boldsymbol{\mathbf{r}}} _1) \,\mathrm{d}{\Omega_1} \\
&= \delta_{l_2, l_2'} \sum_{m_1} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & -m_1 & 0\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ m_1 & -m_1 & 0\end{bmatrix} \times\\
& \qquad (\delta_{l_1',l_1+1} l'_1\mathcal C_{l_1,m_1} - \delta_{l_1,l'_1+1} l_1\mathcal C_{l'_1,m_1})~.
\end{aligned} \end{equation}
积分 $\mathcal C_{l,m}$ 可以参考
式 18 。同理，
\begin{equation} \begin{aligned}F^{(vz2)}_{\lambda,\lambda'} &= \left\langle \mathcal{Y}_{l_1 l_2}^{L} \middle| F^{(vz2)} \middle| \mathcal{Y}_{l_1' l_2'}^{L'} \right\rangle \\
&= \delta_{l_1, l_1'} \sum_{m_1,m_2, m_2'} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & 0\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ m_1 & m'_2 & 0\end{bmatrix} \times\\
& \quad\int Y_{l_2m_2}^*( \hat{\boldsymbol{\mathbf{r}}} _2) F^{(vz2)}( \hat{\boldsymbol{\mathbf{r}}} _2)Y_{l_2' m_2'}( \hat{\boldsymbol{\mathbf{r}}} _2) \,\mathrm{d}{\Omega_1} \\
&= \delta_{l_1, l_1'} \sum_{m_2} \begin{bmatrix}l_1 & l_2 & L\\ -m_2 & m_2 & 0\end{bmatrix} \begin{bmatrix}l_1' & l_2' & L'\\ -m_2 & m_2 & 0\end{bmatrix} \times\\
& \qquad (\delta_{l_2',l_2+1} l_2'\mathcal C_{l_2,m_2} - \delta_{l_2,l'_2+1} l_2 \mathcal C_{l'_2,m_2})~.
\end{aligned} \end{equation}
选择定则： 经过数值验证，速度规范的选择定则同样可以保证 $l_1+l_2-L$ 奇偶性不变。
对称性： 经过数值验证，$F_{(l_1,l_2,L),(l_1',l_2',L')}^{(vz1)} = F_{(l_2,l_1,L),(l_2',l_1',L')}^{(vz2)}$。
6. 非线性偏振电场
对非线性偏振电场，式 18 的前两项不为零，但只和 $z$ 分量的结果大同小异（只有最后一个 3j 系数变为两项，再除以 $\sqrt{2}$ 以归一化）
\begin{equation} \begin{aligned}
& \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| x_1 \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle
= \frac{\delta_{l_2, l_2'} \delta_{M, M'}}{\sqrt2} \sqrt{(2l_1 + 1)(2l_1'+1)} r \times\\
&\qquad \sum_{m_1,m_2} (-1)^{m_1} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2 & L'\\ m_1 & m_2 & M\end{bmatrix} \times\\
&\qquad \begin{pmatrix}l_1 & 1 & l_1'\\ 0 & 0 & 0\end{pmatrix} \left[ \begin{pmatrix}l_1 & 1 & l_1'\\ -m_1 & -1 & m_1\end{pmatrix} - \begin{pmatrix}l_1 & 1 & l_1'\\ -m_1 & 1 & m_1\end{pmatrix} \right] ~.
\end{aligned} \end{equation}
\begin{equation} \begin{aligned}
& \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| y_1 \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle
= \frac{ \mathrm{i} \delta_{l_2, l_2'} \delta_{M, M'}}{\sqrt2} \sqrt{(2l_1 + 1)(2l_1'+1)} r \times \\
& \qquad \sum_{m_1,m_2} (-1)^{m_1} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1' & l_2 & L'\\ m_1 & m_2 & M\end{bmatrix} \times\\
&\qquad \begin{pmatrix}l_1 & 1 & l_1'\\ 0 & 0 & 0\end{pmatrix} \left[ \begin{pmatrix}l_1 & 1 & l_1'\\ -m_1 & -1 & m_1\end{pmatrix} + \begin{pmatrix}l_1 & 1 & l_1'\\ -m_1 & 1 & m_1\end{pmatrix} \right] ~.
\end{aligned} \end{equation}
\begin{equation} \begin{aligned}
& \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| x_2 \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle
= \frac{\delta_{l_1, l_1'} \delta_{M, M'}}{\sqrt2} \sqrt{(2l_2 + 1)(2l_2'+1)} r\\
&\qquad \sum_{m_1,m_2} (-1)^{m_2} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1 & l_2' & L'\\ m_1 & m_2 & M\end{bmatrix} \times\\
&\qquad \begin{pmatrix}l_2 & 1 & l_2'\\ 0 & 0 & 0\end{pmatrix} \left[ \begin{pmatrix}l_2 & 1 & l_2'\\ -m_2 & -1 & m_2\end{pmatrix} - \begin{pmatrix}l_2 & 1 & l_2'\\ -m_2 & 0 & m_1\end{pmatrix} \right] ~.
\end{aligned} \end{equation}
\begin{equation} \begin{aligned}
& \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| y_2 \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle
= \frac{ \mathrm{i} \delta_{l_1, l_1'} \delta_{M, M'}}{\sqrt2} \sqrt{(2l_2 + 1)(2l_2'+1)} r\\
&\qquad \sum_{m_1,m_2} (-1)^{m_2} \begin{bmatrix}l_1 & l_2 & L\\ m_1 & m_2 & M\end{bmatrix} \begin{bmatrix}l_1 & l_2' & L'\\ m_1 & m_2 & M\end{bmatrix} \times\\
&\qquad \begin{pmatrix}l_2 & 1 & l_2'\\ 0 & 0 & 0\end{pmatrix} \left[ \begin{pmatrix}l_2 & 1 & l_2'\\ -m_2 & -1 & m_2\end{pmatrix} + \begin{pmatrix}l_2 & 1 & l_2'\\ -m_2 & 0 & m_1\end{pmatrix} \right] ~.
\end{aligned} \end{equation}
7. 束缚态
能量较低的束缚态可以轻易用虚时间演化得到，即把 $ \exp\left(- \mathrm{i} H \Delta t\right) $ 替换为 $ \exp\left(- H \Delta t\right) $，这相当于把 $\Delta t$ 变为虚数 $- \mathrm{i} \Delta t$，因此得名。氦原子基态只包含 $L = 0$ 的分波，这要求 $l_1 = l_2 = 0, 1, 2, \dots$（参考图 2 ）。
注意氦原子基态只需要由 $L=0$ 的分波展开。另外要特别注意，千万不要对束缚态的分波进行截断，这样会产生类似冲击波的噪音，其强度可以媲美双电离。
图 1：束缚态分波截断噪音（氦原子单电离 110eV，速度规范整体 Lanczos 算法）
图 2：束缚态分波截断噪音 2（缩小时间步长，增加 Krylov 阶数）
图 3：完美波形（Lanczos 算法 11eV 单电离。Krylov 阶数 15，时间步长 0.025，Locatto 节点数 6，等间距 FE 1 au，分波配置 3,3,3）
8. Crank-Nicolson 波函数演化
对薛定谔方程式 5 做 Crank-Nicolson 演化（propagation）同样使用 split operator
\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_{int}\Delta t\right) \exp\left(- \mathrm{i} H_0\frac{\Delta t}{2}\right) + \mathcal{O}\left(\Delta t^3 \right) ~.
\end{equation}
只是这里的 $H_{int} = V_{12} + V_F$ 是所有 $ \left\langle \mathcal{Y} \middle| \ \middle| \mathcal{Y} \right\rangle $ 作用后为非对角（算符）矩阵的项，即电子的相互作用和电场对每个电子的作用。
$H_0 = H_1 + H_2$ 由于 $H_1$ 和 $H_2$ 对易，我们可以将它们独立演化而不产生误差（式 3 ）：
\begin{equation}
\exp\left(- \mathrm{i} H_0 \frac{\Delta t}{2}\right) = \exp\left(- \mathrm{i} H_1 \frac{\Delta t}{2}\right) \exp\left(- \mathrm{i} H_2 \frac{\Delta t}{2}\right) ~.
\end{equation}
也就是对每个 partial wave 的二维网格，用类似氢原子的方法演化每一列再演化每一行。
根据以上的推导可以看出 $H_{int}$ 在展开成径向算符与角向算符的张量积的线性组合后，径向算符都是关于 $r$ 的函数，这对 FEDVR 基底来说都是对角矩阵（即使使用泊松积分）。所以在演化 $ \exp\left(- \mathrm{i} H_{int}\Delta t\right) $ 时，每个二维 FEDVR 格点都可以独立进行，给大量并行提供了可能性。总结起来，程序中的三维波函数矩阵的每个方向的演化都可以每列独立进行。
在演化 $ \exp\left(- \mathrm{i} H_{int}\Delta t\right) $ 的时候，爱华的做法是进一步 split 成三个 operators
\begin{equation}
\exp\left(- \mathrm{i} H_{int}\Delta t\right) = \exp\left(- \mathrm{i} V_F\frac{\Delta t}{2}\right) \exp\left(- \mathrm{i} V_{12} \Delta t\right) \exp\left(- \mathrm{i} V_F\frac{\Delta t}{2}\right) ~.
\end{equation}
这是因为 $V_F = V_{F1} + V_{F2}$ 是含时的，和电场分量 $E_z$ 成正比。
爱华进一步分割的原因是，他的代码中预先储存了 $V_F$ 和 $V_{12}$ 对角化后的本征值和本征矢矩阵，如果直接对 $H_{int}$ 对角化，那么就需要每一步都重新对角化一次。预先储存所需的内存和分波数的平方成正比，然而当我们需要很多分波的时候，内存就远远不够了。为了解决内存问题，我径向和角向都使用 Crank-Nicolson 演化，而且每个时间步长对每个格点重新生成 $V_F + V_{12}$ 矩阵并一次演化不需要分成三次。对于波函数绝对值低于某个阈值的格点，我们甚至可以忽略不演化。
另外一点，如果只关心 single ionization，可以用一个长方形网格（见 Ossiander 2017 Nat. Phys.）。但要注意在短边加入 absorber 防止反射。）
9. Lanczos 波函数演化
Lanczos 算法是 Pazourek 组用的。注意他们并没有使用算符拆分。一口气把整个哈密顿作用在波函数上。它们使用等间距 FE，Krylov 基底的数量在 15 左右。通过调整时间步长来固定误差估计。
10. 分波收敛分析
使用 $l_{1,max}, l_{2,max}, L_{max}$ 来生成分波列表有时候会有较大的浪费，会出现许多概率极小的分波。最好的办法是把各个分波的概率画成图，从图中可以轻易看出分波是否收敛以及是否有浪费。
图 4：$l_{1,max}, l_{2,max}, L_{max} = 5$，可以看出 $L$ 较小时需要更大 $l_1, l_2$，而 $L$ 较大时浪费严重。理想情况是每个 $L$ 中最小概率的分波都差不多。
11. 时间步长误差分析以及变步长演化
首先要注意的是一个电磁场周期至少要有若干个步长，可以在程序中检查。
对于任何演化算法，一个最简单的方法就是每隔一段时间，使用 dt/N（$N = 2$ 或 $4$）为步长，演化 $N$ 次然后和 dt 一次的结果做对比，计算 $ \left\lVert \Delta \Psi' - \Delta \Psi \right\rVert / \left\lVert \Delta \Psi \right\rVert $ 用于衡量相对误差。其中 $\Delta\Psi = \Psi(t+dt) - \Psi(t)$。然而这对 Crank-Nicolson 似乎不准。
对于 Lanzos，有标准的误差公式，也可以和 $ \left\lVert \Delta \Psi' - \Delta \Psi \right\rVert / \left\lVert \Delta \Psi \right\rVert $ 对比看看哪个更准。
12. 分波误差分析（原创）
我们截取有限个分波时，在分波演化（式 34 ）的过程中，会引入新的误差。截取以后，$V_F, V_{12}$ 仍然是厄米矩阵，所以他们的传播子同样是幺正的，所以不能用波函数的模长来估计误差。所以可以用另一种方法来估计误差。误差可以每隔几步估算一次，令 $ \boldsymbol{\mathbf{V}} '_F, \boldsymbol{\mathbf{V}} '_{12}$ 为两倍基底数量的矩阵，那么可以计算 $ \boldsymbol{\mathbf{(}} \mathrm{i} V'_F\Delta t) \boldsymbol{\mathbf{v}} $ 比基底多出来的几个分波，对每个分波的绝对值进行求和。
图 5：分波误差分析
13. 束缚态
理论上可以数值解总哈密顿矩阵的本征值，但实际上这个矩阵大到几乎不可能解出（事实上我们从来不会把总哈密顿矩阵直接储存，而是储存 split operator）。我们使用虚时间演化（imaginary time propagation）。
虚时演化的初始态我用的是两个 He+ 的基态相乘（实数波函数），和爱华略有不同，不过这影响不大。由于我们的 split operator 都是实数矩阵，使用虚时间后其指数函数也是实数矩阵，所以虚时间演化得到的基态波函数也一定会是纯实数。
我们只需用 $L = 0, M = 0$ 的所有 partial waves 即可。
14. 双电子近似
Argon 的电子数为 $N = 18$，基态电离势能为 $I_p = 15.6 \,\mathrm{eV} $（$0.5733 \,\mathrm{au} $），使用等效势能
\begin{equation}
V_{eff} = \frac{-2 - (N-2) \mathrm{e} ^{-r/a}}{r} + V_0 \mathrm{e} ^{-r/b}~,
\end{equation}
满足 $V_{eff}(r\to\infty) \to -2/r$，当 $V_{eff}(r\to 0) \to -N/r$。参数 $a$ 越小屏蔽势能的半径就越小，越接近于氦原子，电离势能越大。
15. Streaking 的分波数
Pazourek 组同时实现了长度和速度规范。[Pazo12] 的 TDSE solver 主要参考是 [Feis08]，这篇讨论的是 sequential 双电离，cross section 使用了最多 $L_{max},l_{1,max},l_{2,max}$ 最大 $(4,9,9)$，且长度规范和速度规范一致。
