库仑势能的球谐展开

贡献者： addis

本文处于草稿阶段。

预备知识　球坐标系中的拉普拉斯方程，库仑定律

　　库仑势能的球谐展开为

\begin{equation} \frac{1}{ \left\lvert \boldsymbol{\mathbf{r}} - \boldsymbol{\mathbf{r}} ' \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}}} ') Y_{l,m}( \hat{\boldsymbol{\mathbf{r}}} )~. \end{equation}

根据式 24 ，令 $\alpha$ 为 $ \hat{\boldsymbol{\mathbf{r}}} , \hat{\boldsymbol{\mathbf{r}}} '$ 的夹角，上式也可以记为

\begin{equation} \frac{1}{ \left\lvert \boldsymbol{\mathbf{r}} - \boldsymbol{\mathbf{r}} ' \right\rvert } = \sum_{l=0}^{\infty} \frac{r_<^l}{r_>^{l+1}} P_l(\cos\alpha) = \sum_{l=0}^{\infty} \sqrt{\frac{4\pi}{2l+1}} \frac{r_<^l}{r_>^{l+1}} Y_{l,0}(\alpha, 0)~. \end{equation}

图 1：左：$1/ \left\lvert \boldsymbol{\mathbf{r}} - \hat{\boldsymbol{\mathbf{z}}} \right\rvert $ 使用式 2 展开，右：左图延 $z$ 轴的截线 $1/ \left\lvert z - 1 \right\rvert $。注意左图中 $r = 1$ 的圆周上有肉眼可见的不收敛，增大 $l$ 的个数即可。

1. 证明

　　和 “平面波的球谐展开” 中的证明类似，我们既可以直接积分得到径向函数，根据式 2 还可以假设 $ \boldsymbol{\mathbf{r}} '$ 在 $z$ 轴。也可以在球坐标中解偏微分方程。

　　除 $ \boldsymbol{\mathbf{r}} '$ 点外，$f( \boldsymbol{\mathbf{r}} ) = \frac{1}{ \left\lvert \boldsymbol{\mathbf{r}} - \boldsymbol{\mathbf{r}} ' \right\rvert }$ 符合拉普拉斯方程。球坐标系中拉普拉斯方程的通解为（式 15 ）

\begin{equation} f( \boldsymbol{\mathbf{r}} ) = f(r, \theta, \phi) = \sum_{l = 0}^\infty \sum_{m = -l}^l \left(C_{l,m} r^l + \frac{C'_{l,m}}{r^{l+1}} \right) P_l^m(\cos\theta) \mathrm{e} ^{ \mathrm{i} m\phi}~, \end{equation}

把空间划分为 $r < r'$ 和 $r > r'$ 两个区域。$r = 0$ 时函数需要有定义，所以 $r < r'$ 区间 $C'_{l,m} = 0$。另外要求 $r \to \infty$ 时 $f( \boldsymbol{\mathbf{r}} ) \to 0$，所以当 $r > r'$，$C_{l,m} = 0$。

　　接下来对比系数即可求出 $C_{l,m}$ 和 $C_{l,m}'$。

未完成：具体计算

2. 附：画图代码

　　图 1 的 Matlab 代码如下，其中 surfSph 函数见 “Matlab 球坐标中的分布图”。SphHarm 函数见 “球谐函数数值计算”。

代码 1：Ylm_Coul_exp.m

% expansion of 1/|r-r'| into Y_lm
% r' in z axis
lmax = 20;
r0 = 1;
Nr = 300; Nth = 300;
r = linspace(1e-5, 2.5, Nr);
mark = r <= r0;
th = linspace(0, 2*pi, Nth); ph = zeros(size(th));
data = zeros(Nr, Nth);
for l = 0:lmax
    % eq_PChYlm_1
    data(mark, :) = data(mark, :) + ...
        sqrt(4*pi/(2*l+1)) * (r(mark).').^l ./ r0^(l+1) ...
        .* SphHarm(l, 0, th, ph);
    data(~mark, :) = data(~mark, :) + ...
        sqrt(4*pi/(2*l+1)) * r0^l ./ (r(~mark).').^(l+1) ...
        .* SphHarm(l, 0, th, ph);
end

[R, Th] = ndgrid(r, th); Ph = zeros(size(Th));
figure;
subplot(1, 2, 1);
surfSph(R, Th, Ph, data);
caxis([0, 5]);
title(['l <= ' num2str(l)]);
colorbar; colormap jet;

subplot(1, 2, 2);
plot(r, data(:,1));
hold on; plot(r, 1./abs(r-r0), '.');
legend('球谐展开', '1/|z-z''|');
axis([0, r(end), 0, 5]);
xlabel z;
title(['l <= ' num2str(l)]);