氦原子数值解 TDSE 笔记

             

预备知识 广义球谐函数

1. 波函数展开

   使用角向基底为两个电子总轨道角动量和 $z$ 分量的本征态,即广义球谐函数 $\mathcal{Y}_{l_1,l_2}^{L,M}( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2)$.

   总波函数所在的空间可以看做角向空间和径向空间的张量积空间.现在有了角向空间的基底,总波函数就可以在该基底上展开

\begin{equation} \left\lvert \Psi \right\rangle = \sum_\lambda \left\lvert R_\lambda \right\rangle \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}

   由张量积空间中的结论,我们可以求哈密顿算符 $H$ 关于角向基底的 “矩阵元” $H_{\lambda, \lambda'}$,每个矩阵元是径向空间中的一个算符.列出薛定谔方程的 “矩阵形式”,就得到了一组 couple 的径向波函数的薛定谔方程

\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)$ 是单个电子的哈密顿算符,对应对角矩阵
\begin{equation} H_i = K_{ri} + \frac{L_i^2}{2m r_i^2} - \frac{2}{r_i} \qquad (i = 1,2) \end{equation}
其中第二项是对角矩阵是因为 $\mathcal Y$ 基底是 $L_i^2$ 的本征函数.$V_{12}$ 是两电子之间的库仑势能,$V_{Fi}$ 是电场对每个电子的作用势能.

2. 波函数的对称性

   我们一般假设电子自旋态恒为 singlet(反对称),所以总波函数满足粒子交换对称 $\Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \Psi( \boldsymbol{\mathbf{r}} _2, \boldsymbol{\mathbf{r}} _1)$.由式 5 可以导出分波的对称性

\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 系数都是实数,这必定是一个实数矩阵.选择定则如下

   爱华的论文中将式 11 用 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{\phi_i}|{ \left[-\frac{1}{2m} \frac{\mathrm{d}^{2}}{\mathrm{d}{r}^{2}} + \frac{l(l+1)}{2mr} \right] }|{\phi_j}\rangle \\ &= -\frac{1}{2m} \langle{\phi_i}|{ \frac{\mathrm{d}^{2}}{\mathrm{d}{r}^{2}} }|{\phi_j}\rangle + \delta_{i,j}\frac{l(l+1)}{2mr_i^2} \end{aligned} \end{equation}
其中 $\phi_i$ 是一维 FEDVR 格点在 $r_i$ 处的正交归一基底,$r_{max}$ 是计算 $T_l$ 时最后一个 FE 的边界.虽然式 15 中含有 $r_{max}$,但结果却与 $r_{max}$ 无关,原因是 $T_l$ 也和 $r_{max}$ 相关,两项中的 $r_{max}$ 抵消了.

4. 电场作用项

   对每个电子,

\begin{equation} V_F = \boldsymbol{\mathbf{E}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} = E_x x + E_y y + 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} \begin{aligned} & \left\langle \mathcal{Y}_{l_1 l_2}^{LM} \middle| z_1 \middle| \mathcal{Y}_{l_1' l_2'}^{L'M'} \right\rangle = 2\sqrt{\frac{\pi}{3}} r \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'} r \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 \iint 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)} r \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)} r \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_1 \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}\\ &\times \begin{bmatrix}1 & l'_1 & l_1\\ 0 & 0 & 0\end{bmatrix} \begin{bmatrix}0 & l'_2 & l_2\\ 0 & 0 & 0\end{bmatrix} \begin{bmatrix}1 & L' & L\\ 0 & 0 & 0\end{bmatrix} \left\{\begin{matrix}1 & l'_1 & l_1\\ 0 & l'_2 & l_2\\ 1 & L' & L\end{matrix}\right\} \end{aligned} \end{equation}
\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 = r\sqrt{9(2l'_1+1)(2l'_2+1)(2L'+1)/(4\pi)^2}\\ &\times \begin{bmatrix}0 & l'_1 & l_1\\ 0 & 0 & 0\end{bmatrix} \begin{bmatrix}1 & l'_2 & l_2\\ 0 & 0 & 0\end{bmatrix} \begin{bmatrix}1 & L' & L\\ 0 & 0 & 0\end{bmatrix} \left\{\begin{matrix}0 & l'_1 & l_1\\ 1 & l'_2 & l_2\\ 1 & L' & L\end{matrix}\right\} \end{aligned} \end{equation}

5. 非线性偏振电场

   对非线性偏振电场,式 17 的前两项不为零,但只和 $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}

6. 束缚态

   能量较低的束缚态可以轻易用虚时间演化得到,即把 $ \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 ).

7. 波函数演化

   演化(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$ 对易,我们可以将它们独立演化

\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 防止反射.)

8. 分波收敛分析

   使用 $l_{1,max}, l_{2,max}, L_{max}$ 来生成分波列表有时候会有较大的浪费,会出现许多概率极小的分波.最好的办法是把各个分波的概率画成图,从图中可以轻易看出分波是否收敛以及是否有浪费.

图
图 1:$l_{1,max}, l_{2,max}, L_{max} = 5$,可以看出 $L$ 较小时需要更大 $l_1, l_2$,而 $L$ 较大时浪费严重.理想情况是每个 $L$ 中最小概率的分波都差不多.

9. 时间步长误差分析以及变步长演化

   首先要注意的是一个电磁场周期至少要有若干个步长,可以在程序中检查.

   对于任何演化算法,一个最简单的方法就是每隔一段时间,使用 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)$

   对于 Lanzos,有标准的误差公式,也可以和 $ \left\lVert \Delta \Psi' - \Delta \Psi \right\rVert / \left\lVert \Delta \Psi \right\rVert $ 对比看看哪个更准.

10. 分波误差分析(原创)

   我们截取有限个分波时,在分波演化(式 29 )的过程中,会引入新的误差.截取以后,$V_F, V_{12}$ 仍然是厄米矩阵,所以他们的传播子同样是幺正的,所以不能用波函数的模长来估计误差.所以可以用另一种方法来估计误差.误差可以每隔几步估算一次,令 $ \boldsymbol{\mathbf{V}} '_F, \boldsymbol{\mathbf{V}} '_{12}$ 为两倍基底数量的矩阵,那么可以计算 $ \boldsymbol{\mathbf{(}} \mathrm{i} V'_F\Delta t) \boldsymbol{\mathbf{v}} $ 比基底多出来的几个分波,对每个分波的绝对值进行求和.

图
图 2:分波误差分析

11. 束缚态

   理论上可以数值解总哈密顿矩阵的本征值,但实际上这个矩阵大到几乎不可能解出(事实上我们从来不会把总哈密顿矩阵直接储存,而是储存 split operator).我们使用虚时间演化(imaginary time propagation)

   虚时演化的初始态我用的是两个 He+ 的基态相乘(实数波函数),和爱华略有不同,不过这影响不大.由于我们的 split operator 都是实数矩阵,使用虚时间后其指数函数也是实数矩阵,所以虚时间演化得到的基态波函数也一定会是纯实数. 我们只需用 $L = 0, M = 0$ 的所有 partial waves 即可.

12. 双电子近似

   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$ 越小屏蔽势能的半径就越小,越接近于氦原子,电离势能越大.

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

广告位

投放详情

         

© 小时科技 保留一切权利