氢原子电离计算（一阶微扰，Matlab）

贡献者： addis

本词条处于草稿阶段。

预备知识　跃迁概率（含时微扰）

　　本文使用原子单位制。

图 1：偶极子绝对值，线：平面波近似，点：无近似（代码：plot_hydrogen_trans_dipole.m）

1. 长度规范（平面波近似）

　　归一化的平面波和归一化的氢原子基态为（用平面波近似末态库仑函数）

\begin{equation} \left\lvert \boldsymbol{\mathbf{k}} \right\rangle = (2\pi)^{-3/2} \mathrm{e} ^{ \mathrm{i} \boldsymbol{\mathbf{k}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} }~, \qquad \left\lvert 0 \right\rangle = \frac{1}{\sqrt{\pi}} \mathrm{e} ^{-r}~. \end{equation}

长度规范下的跃迁偶极子，可以在极坐标系中积分（令 $ \hat{\boldsymbol{\mathbf{k}}} $ 方向为极轴，$ \hat{\boldsymbol{\mathbf{r}}} $ 与其夹角为 $\theta$）

\begin{equation} \left\langle \boldsymbol{\mathbf{k}} \middle| \boldsymbol{\mathbf{r}} \middle| 0 \right\rangle = \frac{ \hat{\boldsymbol{\mathbf{k}}} }{(2\pi)^{3/2}\sqrt{\pi}} \int_0^{+\infty} \int_0^\pi \mathrm{e} ^{-r} \mathrm{e} ^{- \mathrm{i} k r \cos\theta} r \cos\theta \cdot 2\pi r^2 \sin\theta \,\mathrm{d}{\theta} \,\mathrm{d}{r} ~ \end{equation}

换元，令 $u = \cos\theta$，得¹

\begin{equation} \begin{aligned} \left\langle \boldsymbol{\mathbf{k}} \middle| \boldsymbol{\mathbf{r}} \middle| 0 \right\rangle &= \frac{1}{\sqrt 2 \pi} \hat{\boldsymbol{\mathbf{k}}} \int_0^{+\infty} r^3 \mathrm{e} ^{-r} \int_{-1}^1 \mathrm{e} ^{- \mathrm{i} k r u} u \,\mathrm{d}{u} \cdot \,\mathrm{d}{r} \\ &= \mathrm{i} \frac{\sqrt2 \hat{\boldsymbol{\mathbf{k}}} }{k\pi} \int_0^{+\infty} r^2 \mathrm{e} ^{-r} \left[ \cos\left(kr\right) - \frac{1}{kr} \sin\left(kr\right) \right] \,\mathrm{d}{r} \\ &= - \mathrm{i} \frac{8 \sqrt2}{\pi} \frac{ \boldsymbol{\mathbf{k}} }{(k^2+1)^3}~, \end{aligned} \end{equation}

注意这是一个纯虚数。Matlab 代码如下：

代码 1：hydrogen_trans_dipole0.m

% hydrogen transition dipole, approximate Coulomb plane wave with plane wave
% since middle and right parts are symmetric,
%     assuming \bvec k in z+ direction
% <\bvec k|\bvec r|n,0,0> = [0, 0, <\bvec k|z|n,0,0>]
% output numerical integration and analytical result (eq_HyIon2_1)
function [dipole_z, dipole_analy_z] = hydrogen_trans_dipole0(kz, ZZ, n)
% === params ===
rmax = 20; Nr = 200; Nth = 100;
% ==============
k = [0, 0, kz];
r = linspace(0, rmax, Nr); dr = rmax / (Nr-1);
th = linspace(0, pi, Nth); dth = pi / (Nth-1);
ph = 0;
[R, Th] = ndgrid(r, th);
Z = R .* cos(Th);
Psi_n = hydrogen_Psi(ZZ, n, 0, 0, r', th, ph);
% <\bvec k|\bvec r|0>
Psi_k = 1/(2*pi)^(3/2) * exp(1i*kz.*Z);
dipole_z = sum(sum(conj(Psi_k).*Z.*Psi_n .*R.^2.*sin(Th))) ...
    * dr * dth * 2*pi;
% eq_HyIon2_1
if n == 1
    dipole_analy_z = -1i*(8*sqrt(2))/pi * kz / (dot(k,k) + 1)^3;
end
end

2. 速度规范

　　注意一阶微扰理论中的初态和末态波函数都是无微扰（无外场）情况下的，与规范无关。要计算 $ \left\langle \boldsymbol{\mathbf{k}} \middle| \boldsymbol\nabla \middle| 0 \right\rangle $，先看积分

\begin{equation} \int \exp\left(- \mathrm{i} \boldsymbol{\mathbf{k}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \right) \exp\left(-r\right) \,\mathrm{d}^{3}{r} = \frac{8\pi }{(k^2 + 1)^2}~. \end{equation}

使用算符 $ \boldsymbol\nabla $ 的反厄米性得

\begin{equation} \begin{aligned} &\int \exp\left(- \mathrm{i} \boldsymbol{\mathbf{k}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \right) \boldsymbol\nabla \exp\left(-r\right) \,\mathrm{d}^{3}{r} = -\int [ \boldsymbol\nabla \exp\left( \mathrm{i} \boldsymbol{\mathbf{k}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \right) ]^* \exp\left(-r\right) \,\mathrm{d}^{3}{r} \\ &= \mathrm{i} \boldsymbol{\mathbf{k}} \int \exp\left(- \mathrm{i} \boldsymbol{\mathbf{k}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} \right) \exp\left(-r\right) \,\mathrm{d}^{3}{r} = \mathrm{i} \frac{8 \pi \boldsymbol{\mathbf{k}} }{(k^2 + 1)^2}~, \end{aligned} \end{equation}

乘以归一化系数得

\begin{equation} \left\langle \boldsymbol{\mathbf{k}} \middle| \boldsymbol\nabla \middle| 0 \right\rangle = \mathrm{i} \frac{2\sqrt{2}}{\pi}\frac{ \boldsymbol{\mathbf{k}} }{(k^2 + 1)^2}~. \end{equation}

该式代入式 11 （$q = -1$）得微分截面为

\begin{equation} \frac{\partial \sigma}{\partial \Omega} = \frac{32}{mc\omega} \frac{k( \boldsymbol{\mathbf{k}} \boldsymbol\cdot \hat{\boldsymbol{\mathbf{e}}} )^2}{(k^2 + 1)^4} = \frac{64}{mc} \frac{k( \boldsymbol{\mathbf{k}} \boldsymbol\cdot \hat{\boldsymbol{\mathbf{e}}} )^2}{(k^2 + 1)^5}~. \end{equation}

对于质子数为 $Z$ 类氢原子有

\begin{equation} \frac{\partial \sigma}{\partial \Omega} = \frac{32 Z^5}{mc\omega} \frac{k( \boldsymbol{\mathbf{k}} \boldsymbol\cdot \hat{\boldsymbol{\mathbf{e}}} )^2}{(k^2 + Z^2)^4}~. \end{equation}

3. 两种规范对比

　　如果 $ \left\lvert \boldsymbol{\mathbf{k}} \right\rangle $ 是库仑函数（能量本征态）应该有（式 2 ）

\begin{equation} \left\langle \boldsymbol{\mathbf{k}} \middle| \boldsymbol{\mathbf{r}} \middle| 0 \right\rangle = -\frac{ \left\langle \boldsymbol{\mathbf{k}} \middle| \boldsymbol\nabla \middle| 0 \right\rangle }{m\omega_{k0}}~. \end{equation}

其中 $\omega_{k0} = k^2/2 + 1/2$，但实际上式 3 和式 6 满足

\begin{equation} \left\langle \boldsymbol{\mathbf{k}} \middle| \boldsymbol{\mathbf{r}} \middle| 0 \right\rangle = -2\frac{ \left\langle \boldsymbol{\mathbf{k}} \middle| \boldsymbol\nabla \middle| 0 \right\rangle }{m\omega_{k0}}~. \end{equation}

这说明在使用平面波近似库仑函数时，长度规范的 transition amplitude 恰好是速度规范的 2 倍，截面是四倍（待求证）。

　　教材中推导微分截面一般使用速度规范，因为速度规范的结果与实验吻合更好。

4. 使用库仑平面波

长度规范

　　理论上若把上面的平面波换成库仑平面波（库仑势能中的精确散射态），那么理论上用不同的规范结果是一样的。

\begin{equation} \begin{aligned} &\quad \langle{\psi_{ \boldsymbol{\mathbf{k}} }^{(-)}}|{z}|{\psi_{n,l,m}}\rangle = \int \psi_{ \boldsymbol{\mathbf{k}} }^{(-)}( \boldsymbol{\mathbf{r}} )^* z \psi_{nlm}(r,\theta,\phi) r^2 \,\mathrm{d}{r} \,\mathrm{d}{\Omega} \\ &= \sum_{l',m'} \frac{ \mathrm{i} ^{-l'}}{k} \mathrm{e} ^{ \mathrm{i} \sigma_{l'}} Y_{l',m'}( \hat{\boldsymbol{\mathbf{k}}} ) \int \frac{1}{r} \sqrt{\frac{2}{\pi}} F_{l'} \left(\eta, kr \right) Y_{l',m'}^*( \hat{\boldsymbol{\mathbf{r}}} ) \left(r \sqrt{\frac{4\pi}{3}} Y_{1,0}( \hat{\boldsymbol{\mathbf{r}}} ) \right) R_{n,l}(r) Y_{l,m}( \boldsymbol{\mathbf{r}} ) r \,\mathrm{d}{r} \,\mathrm{d}{\Omega} \\ &= \sqrt{\frac{8}{3}}\sum_{l',m'} \frac{ \mathrm{i} ^{-l'}}{k} \mathrm{e} ^{ \mathrm{i} \sigma_{l'}} Y_{l',m'}( \hat{\boldsymbol{\mathbf{k}}} ) \int F_{l'} \left(\eta, kr \right) R_{n,l}(r) r \,\mathrm{d}{r} \int Y_{l',m'}^*( \hat{\boldsymbol{\mathbf{r}}} )Y_{1,0}( \hat{\boldsymbol{\mathbf{r}}} )Y_{l,m}( \boldsymbol{\mathbf{r}} ) \,\mathrm{d}{\Omega} \\ &= \sqrt{\frac{2}{\pi}} \sum_{l'} \sqrt{(2l'+1)(2l+1)} \frac{ \mathrm{i} ^{-l'}}{k} \mathrm{e} ^{ \mathrm{i} \sigma_{l'}} \int_0^\infty F_{l'} \left(\eta, kr \right) R_{n,l}(r) r \,\mathrm{d}{r} \times\\ & \qquad \sum_{m'} (-1)^{m'} Y_{l',m'}( \hat{\boldsymbol{\mathbf{k}}} ) \begin{pmatrix}l'& 1& l\\ 0 & 0 & 0\end{pmatrix} \begin{pmatrix}l' & 1 & l\\ -m' & 0 & m\end{pmatrix} ~. \end{aligned} \end{equation}

其中 $\eta = -Z/k$。对 $l=0$ 有

\begin{equation} \langle{\psi_{ \boldsymbol{\mathbf{k}} }^{(-)}}|{z}|{\psi_{n00}}\rangle = \frac{1}{\sqrt{2}\pi k \mathrm{i} } \mathrm{e} ^{ \mathrm{i} \sigma_1} \cos\theta \int_0^\infty F_1 \left(\eta, kr \right) rR_{n,0}(r) \,\mathrm{d}{r} ~. \end{equation}

　　笔者没有见过该积分的解析解，长度规范的 Matlab 数值积分代码如下：

代码 2：hydrogen_trans_dipole.m

function dipole_z = hydrogen_trans_dipole(kz, ZZ, n)
eta = -ZZ/kz;
f = @(r) coulombF_sym(1, eta, kz*r) .* r .* hydrogen_Rnl(ZZ, n, 0, r);
dipole_z = exp(1i*coulomb_sigma(1,eta)) / (sqrt(2)*pi*kz*1i) ...
    * integral(f, 0, Inf, 'AbsTol',1e-10) ;
end

速度规范

未完成：速度规范数值积分，验证和长度规范一样

　　画图代码：

代码 3：plot_hydrogen_trans_dipole.m

% == params ==
n = 1; ZZ = 1;
kmin = 1; kmax = 3; Nk = 50;
% ============

kz = linspace(kmin, kmax, Nk);
dipole_z0 = zeros(1, Nk);
dipole_z0_ana = dipole_z0;
dipole_z = dipole_z0;

% plane wave approx
for i = 1:Nk
    [dipole_z0(i), dipole_z0_ana(i)]
        = hydrogen_trans_dipole0(kz(i), ZZ, n);
end

% Coulomb wave
parfor i = 1:Nk
    disp(i);
    dipole_z(i) = hydrogen_trans_dipole(kz(i), ZZ, n);
end

figure; plot(kz, abs(dipole_z0)); hold on;
plot(kz, abs(dipole_z0_ana), '--');
plot(kz, abs(dipole_z), '.');
legend({'dipole0', 'dipole0-analy', 'dipole'});
grid on;
xlabel 'k [au]'; ylabel 'abs dipole';

1. ^ 最后一步可通过 Wolfram Alpha 或 Mathematica 获得。

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