氢原子的跃迁偶极子矩阵元计算（束缚态之间）

贡献者： addis; 切糕糕

　　 本文使用原子单位制。

1. $z$ 方向电场，束缚态到束缚态

　　 表 1 给出了核电荷数 $Z=1$ 时的 $ \left\langle \psi_{n',l',m'} \middle| z \middle| \psi_{n,l,m} \right\rangle $，由于这是一个实对称矩阵，只给出矩阵的下半三角。当 $Z > 1$ 时把表中每个矩阵元除以 $Z$ 即可。这是因为 $\psi_{n,l,m}$ 与 $Z$ 成反比进行缩放（保持归一化），导致 $ \left\lvert \psi_{n,l,m} \right\rvert ^2$ 和 $z$ 的平均值也是如此。另外注意 $ \left\langle \psi_{n',l',m'} \middle| z \middle| \psi_{n,l,m} \right\rangle = \left\langle \psi_{n',l',-m'} \middle| z \middle| \psi_{n,l,-m} \right\rangle $。

　　 根据式 3 可以用 3j 系数表示角向积分，进而化简为式 20

　　 这可以用于计算类氢原子斯塔克效应以及跃迁率（未完成）等。

　　 Mathematica 代码如下（请自行修改矩阵尺寸和循环范围），HydrogenR 函数见类氢原子的束缚态。

DipoleZ[Z_, n1_, l1_, m1_, n2_, l2_, 
   m2_] := (-1)^
    m1 Sqrt[(2 l1 + 1) (2 l2 + 1)] ThreeJSymbol[{l1, 0}, {1, 0}, {l2, 
     0}] Integrate[
    HydrogenR[Z, n1, l1, r]\[Conjugate] HydrogenR[Z, n2, l2, 
      r] r^3, {r, 0, +\[Infinity]}] ThreeJSymbol[{l1, -m1}, {1, 
     0}, {l2, m2}];
d = ConstantArray[0, {10, 10}];
i1 = 0; i2 = 0;
For[n1 = 1, n1 <= 4, n1++, For[l1 = 0, l1 < n1, l1++,
  ++i1; i2 = 0;
  For[n2 = 1, n2 <= 4, n2++, For[l2 = 0, l2 < n2, l2++,
    ++i2;
    d[[i1, i2]] = DipoleZ[1, n1, l1, 0, n2, l2, 0];
  ]]
]];
Print[d // MatrixForm];

(* eq_SelRul_1 *)
Dipole[Z_, n1_, l1_, m1_, n2_, l2_, m2_] :=\
(-1)^m1 Sqrt[(2 l1 + 1) (2 l2 + 1)]\
  ThreeJSymbol[{l1, 0}, {1, 0}, {l2, 0}] Integrate[\
  HydrogenR[Z, n1, l1, r]\[Conjugate] HydrogenR[Z, n2, l2, r] r^3, {r,\
    0, +∞}] {(ThreeJSymbol[{l1, -m1}, {1, -1}, {l2, m2}] - \
     ThreeJSymbol[{l1, -m1}, {1, 1}, {l2, m2}])/Sqrt[\
   2], (ThreeJSymbol[{l1, -m1}, {1, -1}, {l2, m2}] + \
     ThreeJSymbol[{l1, -m1}, {1, 1}, {l2, m2}])/Sqrt[2], \
  ThreeJSymbol[{l1, -m1}, {1, 0}, {l2, m2}]}

　　 Matlab 代码如下（请自行修改矩阵尺寸和循环范围），hydrogen_Rnl 函数见类氢原子的束缚态。

% hydrogen <n1,l1,m1|z|n2,l2,m2>
% r_max is upper bound of r integral
function d_z = hydrogen_trans_dipole_z(Z,n1,l1,m1,n2,l2,m2,r_max)
if (n1 <= 0 || l1 < 0 || l1 >= n1 || abs(m1) > l1)
    error('illegal n1,l1,m1');
end
if (n2 <= 0 || l2 < 0 || l2 >= n2 || abs(m2) > l2)
    error('illegal n2,l2,m2');
end
if (abs(l1-l2) ~= 1 || m1 ~= m2)
    d_z = 0; return;
end
integrand = @(r) hydrogen_Rnl(Z,n1,l1,r) ...
    .* hydrogen_Rnl(Z,n2,l2,r) .* r.^3;
% radial integral
% change of variable r = tan(x) to avoid infinite upperbound
I_r = integral(integrand, 0, r_max);
if isnan(I_r) || isinf(I_r)
    error('integral failed!');
end
% angular integral
I_th = (-1)^m2 * sqrt((2*l1+1)*(2*l2+1)) ...
    * ThreeJ(l1,0,1,0,l2,0) * ThreeJ(l1,-m1,1,0,l2,m2);
d_z = I_r * I_th;
end

　　 以下函数需要 parallel toolbox 加速，如果没有可能会较慢（parfor 自动变为 for）。

% hydrogen <n1,l1,m1|z|n2,l2,m2> matrix
% n_l(i,:) = [n, l] for each basis
function [d_z, n_l] = hydrogen_trans_dipole_z_mat(Z, m, n_max, r_max)
if nargin == 0
% === params ===
Z = 1;
m = 0;
n_max = 5;
% ==============
end
Ndim = (n_max-abs(m))*(n_max-abs(m)+1)/2;
n_l = zeros(Ndim, 2);
n1l1_n2l2_k1k2 = zeros(Ndim*(Ndim-1)/2, 6);
k1 = 0; k = 0;
for n1 = 1:n_max
    for l1 = m:n1-1
        k1 = k1 + 1; k2 = 0;
        n_l(k1, :) = [n1, l1];
        for n2 = 1:n_max
            for l2 = m:n2-1
                k2 = k2 + 1;
                if (k2 >= k1)
                    continue;
                end
                k = k + 1;
                n1l1_n2l2_k1k2(k, :) = [n1, l1, n2, l2, k1, k2];
            end
        end
    end
end
if k ~= Ndim*(Ndim-1)/2 || k1 ~= Ndim
    error('error!')
end

Nlist = size(n1l1_n2l2_k1k2, 1);
d_val = zeros(1, Nlist);
parfor i = 1:Nlist
    n1 = n1l1_n2l2_k1k2(i,1); l1 = n1l1_n2l2_k1k2(i,2);
    n2 = n1l1_n2l2_k1k2(i,3); l2 = n1l1_n2l2_k1k2(i,4);
    disp(n1l1_n2l2_k1k2(i,1:4));
    d_val(i) = hydrogen_trans_dipole_z(Z,n1,l1,m,n2,l2,m,r_max);
end

d_z = zeros(Ndim, Ndim);
for i = 1:Nlist
    k1 = n1l1_n2l2_k1k2(i,5); k2 = n1l1_n2l2_k1k2(i,6);
    d_z(k1, k2) = d_val(i);
end
% copy to upper triangle
d_z = d_z + d_z.';
end

2. 好本征态列表（z 方向电场）

　　 注意第 $n$ 能级有 $n^2$ 基底。

　　 【$n=2, m=0$】

　　 【$n=3, m=0$】

　　 【$n=3, m=1$】

　　 【$n=4, m=0$】

　　 【$n=4, m=1$】

　　 【$n=4, m=2$】

　　 【$n=5, m=0$】

　　 【$n=5, m=1$】

　　 【$n=5, m=2$】

　　 【$n=5, m=3$】