从线性系统的视角看电磁势场的三类方程

从线性系统的角度出发，静电场与电磁波场的势函数都满足线性偏微分方程：

$\mathcal{L}\,u(\mathbf r) = s(\mathbf r),\tag{1}$

其中 $\mathcal{L}$ 是空间算子， $s(\mathbf r)$ 是源分布。
对应三种典型情形：

方程类型	数学形式	物理对应
拉普拉斯方程	$\nabla^2 u = 0$	无源静态势场
泊松方程	$\nabla^2 u = -\rho/\varepsilon_0$	有源静态势场
亥姆霍兹方程	$(\nabla^2 + k^2) u = -f(\mathbf r)$	有传播的时谐电磁波场

下面以亥姆霍兹方程为例，计算它的格林函数的形式

从线性系统出发：点源响应与卷积叠加

因为系统是线性的，可以借用信号系统的思想：

求解“点源”输入（δ 函数）的输出；
再用卷积叠加任意源。

于是定义格林函数 $G(\mathbf r, \mathbf r')$ ：

$\mathcal{L}\,G(\mathbf r, \mathbf r') = -\delta(\mathbf r - \mathbf r'),\tag{2}$

它正是系统对 δ 源的响应。

因此任意源的场是所有点源响应的叠加：

$u(\mathbf r) = \int G(\mathbf r,\mathbf r')\,s(\mathbf r')\,d^3r' + u_\text{hom}(\mathbf r),\tag{3}$

其中：

卷积积分 = 特解（由源决定）；
齐次解 $u_\text{hom}$ = 通解部分（由边界条件决定）。

这与常微分方程的「特解 + 通解」结构完全对应。

为了求得格林函数，引入算子的本征函数展开

算子 $\mathcal{L}$ 是线性且自伴的（self-adjoint operator，也称厄米算符，Hermitian），
因此可以对角化，存在一组完备的本征函数 $\{\phi_{lm}\}$ ：

$\mathcal{L}\,\phi_{lm} = \lambda_l \phi_{lm}.\tag{4}$

于是可形式上写出算子的逆：

$\mathcal{L}^{-1} = \sum_{l,m} \frac{|\phi_{lm}\rangle \langle \phi_{lm}|}{\lambda_l},\tag{5}$

从而得到格林函数的模态展开：

$G(\mathbf r,\mathbf r') = \sum_{l,m} \frac{\phi_{lm}(\mathbf r)\,\phi_{lm}^*(\mathbf r')}{\lambda_l}.\tag{6}$

所以求 $G$ 的本质，就是寻找算子的本征函数体系。

算子可对角化 → 分离变量

由于拉普拉斯算子在球坐标下可写为：

$\nabla^2 = \frac{1}{r^2}\frac{\partial}{\partial r}\!\left(r^2 \frac{\partial}{\partial r}\right) + \frac{1}{r^2}\nabla_\Omega^2,\tag{7}$

其中 $\nabla_\Omega^2$ 仅作用在角变量 $(\theta,\phi)$ 。

这意味着算子可以被分离为径向部分和角向部分，且它们对易（可同时对角化）。
因此可以假设：

$u(r,\theta,\phi) = R(r)\,Y(\theta,\phi),\tag{8}$

代入方程后得到：

$\frac{1}{R}\frac{d}{dr}\!\left(r^2\frac{dR}{dr}\right) + k^2r^2 = -\frac{1}{Y}\nabla^2_\Omega Y.\tag{9}$

两边分别依赖于独立变量 $r$ 和 $(\theta,\phi)$ ，因此必须分别等于同一个常数 $l(l+1)$ 。

角向方程 → 球谐函数（角动量本征态）

这里奇怪的是为什么两边必须等于同一个常数，而且这个常数还必须是一个关于整数 $l$ 的式子 $l(l+1)$ ?

首先，左右两边如果不同时等于一个常数，而是一个依赖于左边或右边自变量的因变量，那么就必定存在一个自变量使得左右两式不相等。或者说，正是因为左右两边相互独立，所以它们必须等于一个无法依赖左右两边的量，即一个常数量才行。
$l$ 为什么是一个整数，而且常数量刚好是 $l(l+1)$ ?

这里就涉及到对角向方程（角动量算子的本征方程）的求解了。我们一开始不知道这个常数量是多少，因此设为 $\lambda$ ，此时角方程写成：

$\left\{ \begin{aligned} &\nabla_\Omega^2 = \frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\!\left(\sin\theta\frac{\partial}{\partial\theta}\right) + \frac{1}{\sin^2\theta}\frac{\partial^2}{\partial\phi^2}.\\ &\nabla^2_\Omega Y=-\lambda Y \end{aligned} \right.\tag{10}$

同样的，这里的 $Y$ 还能写成 $Y(\theta,\phi) = \Theta(\theta)\,\Phi(\phi)$ ，这个 $Y$ 的本征方程函数解即为球谐函数，带入到上面角动量算子方程中

$\frac{\sin\theta}{\Theta}\frac{d}{d\theta}\!\left(\sin\theta \frac{d\Theta}{d\theta}\right) + \lambda\sin^2\theta = -\,\frac{1}{\Phi}\frac{d^2\Phi}{d\phi^2}.\tag{11}$

这里再来一次分离变量，左右两边都独立且等于 $M$ ，则

$\frac{d^2\Phi}{d\phi^2} = -M \Phi.\tag{12}$

这个方程的解太简单了，就是一个复指数函数

$\Phi(\phi) = A e^{im\phi} + B e^{-im\phi},\quad m=\pm \sqrt{M}\tag{13}$

所以我们这里干脆直接定义这个常数为 $m^2$ ，又因为物理场必须在绕一圈（ $\phi\to \phi+2\pi$ ）后保持单值，

$\Phi(\phi+2\pi)=\Phi(\phi) \Rightarrow e^{im2\pi}=1 \Rightarrow m \in \mathbb{Z}.$

于是 $m$ 是一个整数，称为方位量子数（azimuthal quantum number）,某些地方也叫做磁量子数。

为了保证正交归一性（ $\int_0^{2\pi}\Phi_m(\phi)\Phi^*_{m'}(\phi)d\phi=\delta_{mm'}$ ），仅保留单一方向的指数项 $e^{im\phi}$ 。

带回原角动量算子方程(11)中

$\frac{1}{\sin\theta}\frac{d}{d\theta}\!\left(\sin\theta\frac{d\Theta}{d\theta}\right) + \Big[\lambda - \frac{m^2}{\sin^2\theta}\Big]\Theta = 0.\tag{14}$

这个方程就是连带勒让德方程（Associated Legendre Equation），当 $m=0$ 时，为普通勒让德方程。

把变量换成 $x = \cos\theta$ ，

$(1-x^2)\frac{d^2\Theta}{dx^2} - 2x\frac{d\Theta}{dx} + \left[\lambda - \frac{m^2}{1-x^2}\right]\Theta = 0.\tag{15}$

这里需要用级数法（幂级数展开）求解。考虑到在端点 $x=\pm1$ （即 $\theta=0,\pi$ 处）存在 $\frac{m^2}{1-x^2}$ 发散项，普通幂级数（如 $\sum a_nx^n$ ）无法直接收敛，需要用到 Frobenius 方法（弗罗贝尼乌斯级数法）去分析局部幂级数解。

Frobenius 级数形式：

$\Theta(x) = (1-x)^{s}\sum_{n=0}^\infty a_n(1-x)^n,$

该思想是当方程系数含有 $1/(1-x)^2$ 或 $1/(1-x)$ 这样的奇点时，在代入级数展开后，最低阶项的系数通常给出一个“指标方程（indicial equation）”，用于确定 $s$ 的取值。这一步本质上是：“找出解在奇点附近的最强幂次因子，使得方程不再发散”。于是，在 $x\approx 1$ 处做局部分析，当 $x\to1$ 时，

$1-x^2 \approx 2(1-x).\tag{16}$

方程主导项近似为：

$(1-x^2)\Theta'' - \frac{m^2}{1-x^2}\Theta \approx 0.\tag{17}$

项	含因子	当 $x\to1$ 的幂次
$(1-x^2)\Theta''$	约 $(1-x^2)^{s-1}$	指数 $s-1$
$-2x\Theta'$	约 $(1-x^2)^{s-1/2}$	指数 $s-1/2$ ，弱得多
$-\frac{m^2}{1-x^2}\Theta$	约 $(1-x^2)^{s-1}$	同阶
$+\lambda\Theta$	约 $(1-x^2)^{s}$	更弱

因此只关注发散最大的（17）中的两个项，它们为0就可以了。

设 $\Theta \sim (1-x)^s$ ，则：

$\Theta' \sim s(1-x)^{s-1}, \quad \Theta'' \sim s(s-1)(1-x)^{s-2}.\tag{18}$

代入得：

$(1-x^2)\Theta'' \sim 2s(s-1)(1-x)^{s-1}, \quad \frac{m^2}{1-x^2}\Theta \sim \frac{m^2}{2}(1-x)^{s-1}.\tag{19}$

平衡两边主导幂次 → 系数相等：

$2s(s-1) = \frac{m^2}{2}.\tag{20}$

得到指标方程：

$s(s-1) = \frac{m^2}{4}\quad\to\quad s=\frac{1\pm\sqrt{1+m^2}}{2}.\tag{21}$

首先 $s<0$ 的根要舍弃，因为这会导致 $\Theta \sim (1-x)^s$ 发散，我本来就是要假设一个 Frobenius 级数形式的解靠 $(1-x)^s$ 去抵消掉 $s<0$ 的奇点，因此 $s=\frac{1+\sqrt{1+m^2}}{2}$ ， $\Theta \sim (1-x)^{\frac{1+\sqrt{1+m^2}}{2}}$

带入回（19）中， $(1-x)^{s-1}$ 的指数 $s-1=\frac{\sqrt{1+m^2}-1}{2}$ 需要大于等于零，就能保证发散项得到收敛。

这个指数有点复杂，不利于后面计算，我们仅需要抵消掉奇异项就行，因此挑选一个 $s-1=\frac{\sqrt{1+m^2}-1}{2}\ge 0$ 的简单的值就行。根据两边之和大于等于第三边 $\sqrt{1}+\sqrt{m^2}\ge\sqrt{1+m^2}\to|m|/2\ge\frac{\sqrt{1+m^2}-1}{2}$ ，由此得到了教科书上谜之假设，设

$\Theta(x) = (1-x^2)^{|m|/2}\,f(x)\quad, f(x)=\sum_{n=0}^\infty a_n(1-x)^n.\tag{22}$

教科书上 99% 不会告诉你为什么要这么假设，好像它天生就知道解的形式要长这样。但我这里讲清楚了为什么要这样假设，因为指标方程求解出来幂次需要拿去抵消分母，从而消除发散项。然后将这个形式解带入（15）得到修正后（消除奇异项）的勒让德方程

$(1-x^2)\frac{d^2f}{dx^2} - 2(1+|m|)x\frac{df}{dx} + [\lambda - |m|(|m|+1)]f = 0.\tag{23}$

这个时候就能使用幂级数展开的方法，找递推关系来求解。将（22）的级数表达式带入（23）并逐幂比较，得到对每个 $n$ 的递推：

$a_{n+2}=-\,\frac{\lambda-(n+|m|)(n+|m|+1)}{(n+2)(n+1)}\,a_{n}\;.\tag{24}$

由于这个多项式如果不在某一阶终止的话，考虑 $n\to\infty$ ， $a_{n+2}/a_n\sim1$ ， $f(x)$ 在 $x=1$ 处发散，与实际存在的有限物理场的解相违背，因此存在某阶 $N$ 使级数在此截断，必有

$λ=(N+|m|)(N+|m|+1).\tag{25}$

把 $l=N+|m|$ （即 $l\ge |m|$ ）记作角量子数，则

$\boxed{\;\lambda=l(l+1),\qquad m=-l,-l+1,\dots, l-1,l.}\tag{26}$

这里要记住，我们要求的是函数 $\Theta(x)$ ，它的函数解是连带勒让德多项式 $P_l^m(x)$ ，这个解与普通勒让德方程的解勒让德多项式 $P_l(x)$ 之间存在特殊关系（实际上就是提取了指标幂次并左移了 $m$ 阶）：

$\boxed{ P_l^m(x) = (-1)^m (1-x^2)^{m/2} \frac{d^m}{dx^m} P_l(x), \quad m\ge 0. }\tag{27}$

这里需要注意，定义是 $m\ge0$ ，要求 $m<0$ 的情况需要根据对称性和积分正交性得到

$\boxed{ P_l^{-m}(x) = (-1)^m \frac{(l-m)!}{(l+m)!}\,P_l^m(x),\quad m>0. }\tag{28}$

这里暂时不解释为什么是这样的，因为这需要复杂的数学证明，我们只需要记住角向方程被我们求解出来了，解是球谐函数，而球谐函数是连带勒让德多项式与方位角方向复指数函数的乘积（ $N_{l,m}$ 为归一化系数）

$\boxed{Y_l^m(\theta,\phi)=N_{l,m}\cdot P_l^m(\cos\theta)\cdot e^{im\phi}}\tag{29}$

径向方程 → 球贝塞尔函数（径向传播模态）

结合前面求解的角向方程，径向方程可以写为：

$\frac{1}{r^2}\frac{d}{dr}\!\left(r^2\frac{dR}{dr}\right) + \left[k^2 - \frac{l(l+1)}{r^2}\right]R = 0,\tag{30}$

即球贝塞尔方程，其解为：

$R(r) = A j_l(kr) + B n_l(kr),\tag{31}$

对应径向方向的振荡模式（内解与外解）。

这里奇怪的是球贝塞尔方程是啥？你凭什么说它是个球贝塞尔方程？内解和外解是什么？

首先介绍贝塞尔方程，它和球贝塞尔方程有什么关系。贝塞尔方程的起源是求解二维圆柱坐标下的亥姆霍兹方程，它最有名的应用是求解圆柱波导中光场的传播模式。

在平面波动或电磁场问题中，我们有：

$(\nabla^2 + k^2)\psi(r,\phi) = 0.\tag{32}$

在圆柱坐标 $(r,\phi)$ 中，拉普拉斯算子是：

$\nabla^2 = \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\right) +\frac{1}{r^2}\frac{\partial^2}{\partial \phi^2}.\tag{33}$

代入分离变量后得到：

$\frac{r}{R}\frac{d}{dr}\left(r\frac{dR}{dr}\right)+k^2 r^2=\frac{1}{\Phi}\frac{d^2\Phi}{d\phi^2}=m^2. \tag{34}$

同样的，右边是方位角的解，配合周期性和归一化要求， $\Phi=e^{im\phi},m\in\mathbb{Z}$

现在来看左边二维柱坐标下的径向方程

$r^2 R'' + r R' + (k^2 r^2 - m^2)R=0.\tag{35}$

对比最基本的贝塞尔方程的定义几乎完全一样

$x^2 y'' + x y' + (x^2 - \nu^2)y = 0.\tag{36}$

利用前面提到的幂级数法进行求解，设

$y(x) = \sum_{n=0}^{\infty} a_n x^{n+s}, \quad a_0\neq 0.\tag{37}$

带入到（36）中并展开

$x^2\sum a_n(n+s)(n+s-1)x^{n+s-2} +x\sum a_n(n+s)x^{n+s-1} +(x^2-\nu^2)\sum a_n x^{n+s}=0.\tag{38}$

把所有项合并成 $x^{n+s}$ 的系数：

$\sum a_n\big[(n+s)^2-\nu^2\big]x^{n+s} +\sum a_{n-2}x^{n+s}=0.\tag{39}$

得到递推关系

$a_n = -\frac{a_{n-2}}{(n+s)^2 - \nu^2}.\tag{40}$

取 $n=0$ 项，由于不存在 $a_{n-2}$ 的项

$(s^2-\nu^2)a_0=0.\tag{41}$

因为 $a_0\ne0$ ，得

$\boxed{s=\pm\nu.}\tag{42}$

这说明有两组独立的级数（对应 $x^{+\nu}$ 和 $x^{-\nu}$ 两类解）。

取 $s=\nu$ ，那么（40）为：

$a_n = -\frac{a_{n-2}}{(n+\nu)^2 - \nu^2} = -\frac{a_{n-2}}{n(n+2\nu)}. \tag{43}$

从这里我们看到，奇偶项分离成两支级数。我们通常取最简单的正则解（原点有限），即只保留偶数项：

$a_1=0,\quad a_0 \text{ 任意。}$

由（43）：

$a_2 = -\frac{a_0}{2(2+2\nu)},\quad a_4 = -\frac{a_2}{4(4+2\nu)} = \frac{a_0}{2\cdot4(2+\nu)(4+\nu)},\;\text{等等。}$

由此可以归纳出一般项：

$a_{2k} = \frac{(-1)^k a_0}{2^{2k} k!\,\Gamma(\nu+k+1)/\Gamma(\nu+1)} = \frac{(-1)^k a_0\Gamma(\nu+1)}{k!\,\Gamma(k+\nu+1)}\left(\frac{1}{2}\right)^{2k}. \tag{44}$

得到级数形式的解为

$y(x) = x^{\nu}\sum_{k=0}^{\infty} a_{2k} x^{2k} = a_0 x^{\nu}\sum_{k=0}^{\infty}\frac{(-1)^k\Gamma(\nu+1)}{k!\,\Gamma(k+\nu+1)}\left(\frac{x}{2}\right)^{2k}=a_0\Gamma(\nu+1)2^\nu\sum_{k=0}^{\infty}\frac{(-1)^k}{k!\,(k+\nu)!}\left(\frac{x}{2}\right)^{2k+\nu}.\tag{45}$

这里的 $a_0$ 需要根据归一化来求，我们记右侧级数为贝塞尔多项式 $J_n(x)$

$J_\nu(x)=\sum_{k=0}^{\infty}\frac{(-1)^k}{k!\,(k+\nu)!}\left(\frac{x}{2}\right)^{2k+\nu}.\tag{46}$

当 $x\to 0$ 时，主导项为最低阶无穷小，即 $k=0$ 项保留， $J_\nu(x\to0)=\frac{1}{\Gamma(\nu+1)}(\frac{x}{2})^\nu$ ，令 $J_0(0)=1$ （后面介绍为什么要让这个等于1）

$a_0=\frac{1}{2^\nu\Gamma(\nu+1)}\quad\Rightarrow\quad y(x)=J_\nu(x), \quad s=\nu\ \text{的解}$

取 $s=-\nu$ ，类似推导得到另一支解

$J_{-\nu}(x) = \sum_{k=0}^{\infty} \frac{(-1)^k}{k!\,\Gamma(k-\nu+1)} \left(\frac{x}{2}\right)^{2k-\nu}.\tag{47}$

这里发现当 $\nu=n$ 为整数时

$J_{-n}(x)=(-1)^n J_n(x),\tag{48}$

两个解不再独立，因此还需要定义另一个独立的解，也就是为什么会有第二类贝塞尔函数 $Y_\nu(x)$

上面是贝塞尔函数的级数定义，它还有一个更为本质，不同视角下的定义——积分定义

考虑亥姆霍兹方程在笛卡尔坐标系下平面波的解

$(\nabla^2 + k^2)\psi(r,\phi) = 0\qquad\Rightarrow\qquad \psi = e^{i k x} = e^{i k r \cos\phi}.\tag{49}$

这个解在圆柱坐标系下，平面波可以展开为各个角向模式的叠加：

$e^{i k r \cos\phi} = \sum_{m=-\infty}^{\infty} i^m J_m(k r) e^{i m \phi}. \tag{50}$

这个式子非常关键，它说明贝塞尔函数 $J_m(kr)$ 是“平面波在圆柱坐标下的模式函数”。

于是贝塞尔函数可以理解为：

它是角度平均的平面波相位的结果。

如何理解这句话呢？从形式上可以看出，（50）相当于对 $e^{ikr\cos\phi}$ 进行了傅里叶展开，每个等间隔角度对应了一个方位角频谱点，该频谱点的幅度为 $i^mJ_m(kr)$ ，它们在方位角频域上共同构成了平面波的空间频谱。

因此，当需要将笛卡尔坐标下亥姆霍兹方程的解转变到圆柱坐标下时，就需要将原本笛卡尔坐标空间的频谱转变为圆柱坐标下方位角维度上的频谱！

现在我们想知道这个方位角频谱中，各个频谱点的系数 $i^mJ_m(kr)$ 是怎么来的，根据 $\cos\alpha=\frac{1}{2}(e^{i\theta}+e^{-i\theta})$ ，令 $t=ie^{i\theta},\ z=kr$ ，（50）可看作

$e^{\frac{iz}{2}(e^{i\theta}+e^{-i\theta})}=e^{\frac{z}{2}(t-\frac{1}{t})}=\sum_{n=-\infty}^{\infty}a_n(z)t^n\tag{51}$

已知 $e^x=\sum_{m=0}^{\infty}x^m/m!$ ，（51）可展开为

$e^{\frac{z}{2}(t-\frac{1}{t})} =\sum_{m=0}^\infty \frac{1}{m!}\left(\frac{z}{2}(t-\frac{1}{t})\right)^m =\sum_{m=0}^\infty \frac{1}{m!}\left(\frac{z}{2}\right)^m\sum_{k=0}^m \binom{m}{k} t^{k}( -t^{-1})^{m-k} =\sum_{m=0}^\infty \frac{1}{m!}\left(\frac{z}{2}\right)^m\sum_{k=0}^m(-1)^{m-k}\binom{m}{k}t^{2k-m}.\tag{52}$

其中 $\binom{m}{k}$ 是组合数 $C^k_m$ ，有 $\binom{m}{k}=\frac{m!}{k!(m-k)!}$

令 $n=2k-m$ ，即 $k=\frac{m+n}{2}$ ，那么只有当 $m+n$ 为偶数时（否则 $k$ 不是整数），该项才存在

于是系数 $a_n(z)$ 为：

$a_n(z)=\sum_{m=0}^{\infty}\sum_{\substack{n=-m\\m+n\ \text{even}}}^{m} \frac{(-1)^{(m-n)/2}}{(\frac{m-n}{2})!(\frac{m+n}{2})!} \left(\frac{z}{2}\right)^m.\tag{53}$

修改索引 $p=(m-n)/2$ ，（53）可写为

$a_n(z)=\sum_{p=0}^\infty \frac{(-1)^p}{p!(p+n)!} \left(\frac{z}{2}\right)^{2p+n}.\tag{54}$

对比（46）

$J_\nu(x)=\sum_{k=0}^{\infty}\frac{(-1)^k}{k!\,(k+\nu)!}\left(\frac{x}{2}\right)^{2k+\nu}\Rightarrow a_n(z)=J_\nu(x)\Rightarrow e^{ikr\cos\phi}=\sum_{n=-\infty}^{\infty}J_n(kr)(ie^{i\phi})^n =\sum_{m=-\infty}^{\infty}i^mJ_m(kr)e^{im\phi}.\tag{55}$

由此可见，贝塞尔函数是平面波在圆柱坐标下傅里叶变换的空间频谱系数

到目前为止我们都是从级数方法去求解的，如果直接对平面波 $e^{ikr\cos\phi}$ 做空间傅里叶变化呢？

我们把（50）两边都乘上 $e^{-i m' \phi}$ ，并对 $\phi$ 从 $0$ 到 $2\pi$ 积分：

$\int_0^{2\pi} e^{i k r \cos\phi}\,e^{-i m'\phi}\,d\phi = \sum_m i^m J_m(k r)\int_0^{2\pi} e^{i(m-m')\phi}\,d\phi=2\pi i^{m'}J_{m'}(kr).\tag{56}$

得到

$J_{m'}(kr) = \frac{1}{2\pi i^{m'}} \int_0^{2\pi} e^{i(kr\cos\phi - m'\phi)}\,d\phi=\frac{1}{2\pi}\int_0^{2\pi} e^{i(kr\cos\phi - m'\phi-m'\pi/2)}\,d\phi.\tag{57}$

令 $\theta=\phi+\pi/2,\ m=m',\ x=kr$

$\boxed{ J_{m}(x) =\frac{1}{2\pi}\int_0^{2\pi} e^{i(x\sin\theta - m\theta)}\,d\theta. }\tag{58}$

这便是贝塞尔函数的积分定义。考虑到贝塞尔函数为实函数，且积分从 $0$ 到 $2\pi$ 的指数函数是周期的，利用欧拉公式分开（58）的实部并利用积分区域的平移和对称性

$\boxed{ J_{\nu}(x) =\frac{1}{\pi}\int_0^{\pi} \cos(\nu \tau-x\sin\tau)\,d\tau. }\tag{59}$

这里的结果就和教科书上的是一模一样的了，现在就可以来讨论为什么第二类贝塞尔函数要写成

$Y_\nu(x)=\frac{J_\nu(x)\cos(\nu\pi)-J_{-\nu}(x)}{\sin(\nu\pi)}.\tag{60}$

根据贝塞尔函数的积分定义，将 $-\nu$ 替换进（59）

$\boxed{ J_{-\nu}(x) =\frac{1}{\pi}\int_0^{\pi} \cos(\nu \tau+x\sin\tau)\,d\tau. }\tag{60}$

（59）和（60）在 $\nu$ 为整数时线性相关

$J_{-\nu}(x)=(-1)^\nu J_\nu(x)\quad(\nu\in\mathbb Z).\tag{61}$

为什么呢？明明被积函数一个减一个加，却还能在 $\nu$ 为整数时线性相关？
设

$I_+(x)=\int_0^\pi \cos\!\big(\nu\tau + x\sin\tau\big)\,d\tau,\qquad I_-(x)=\int_0^\pi \cos\!\big(\nu\tau - x\sin\tau\big)\,d\tau.\tag{62}$

对 $I_+(x)$ 做代换 $u=\pi-\tau$ （因此 $\sin\tau=\sin u$ 且 $d\tau=-du$ ）：

$I_+(x)=-\int_\pi^0 \cos\!\big(\nu(\pi-u)+x\sin u\big)\,du = \int_0^\pi \cos\!\big(\nu\pi-\nu u + x\sin u\big)\,du.\tag{63}$

用余弦差角公式

$\cos(A-B)=\cos A\cos B+\sin A\sin B\tag{64}$

取 $A=\nu\pi$ , $B=\nu u - x\sin u$ ，得

$I_+(x)=\cos(\nu\pi)\!\int_0^\pi \cos\!\big(\nu u - x\sin u\big)\,du \;+\;\sin(\nu\pi)\!\int_0^\pi \sin\!\big(\nu u - x\sin u\big)\,du.\tag{65}$

当 $\nu\in\mathbb Z$ 时， $\sin(\nu\pi)=0$ ， $\cos(\nu\pi)=(-1)^\nu$ 。于是

$J_{-\nu}(x)=\frac{1}{\pi}I_+(x)=\frac{(-1)^\nu}{\pi}\int_0^\pi \cos\!\big(\nu u - x\sin u\big)\,du = \frac{(-1)^\nu}{\pi}I_-(x)=(-1)^\nu J_\nu(x).\tag{66}$

因此，根据前面贝塞尔函数的级数定义可知，虽然我们求得 $s=\pm\nu$ 时两个解，但这两个解在 $\nu$ 为整数时线性相关，不能组合出所有可能的解。基于此，需要构造另一个在 $\nu$ 为整数时与 $J_\nu(x)$ 线性无关的解。

根据 $J_\nu(x)$ 的积分定义式，这是一个余弦函数，虽然内部又包含了一个正弦函数，但不妨碍我尝试将外部余弦函数改为正弦，用来作为它自己的另一个线性无关解，记

$Y_\nu(x)=-\frac{1}{\pi}\int_0^{\pi} \sin(\nu \tau-x\sin\tau)\,d\tau.\tag{67}$

对 $J_{-\nu}(x)$ 做代换 $u=\pi-\tau$ （因此 $\sin\tau=\sin u$ 且 $d\tau=-du$ ）：

$J_{-\nu}(x)=\cos(\nu\pi)\!\int_0^\pi \cos\!\big(\nu u - x\sin u\big)\,du \;+\;\sin(\nu\pi)\!\int_0^\pi \sin\!\big(\nu u - x\sin u\big)\,du=\cos(\nu\pi)J_\nu+\sin(\nu\pi)(-Y_\nu(x)).\tag{68}$

整理得到

$Y_\nu(x)=-\frac{1}{\pi}\int_0^{\pi} \sin(\nu \tau-x\sin\tau)\,d\tau=\frac{J_\nu(x)\cos(\nu\pi)-J_{-\nu}(x)}{\sin(\nu\pi)}.\tag{69}$

知道了第一类贝塞尔函数 $J_\nu(x)$ 和第二类贝塞尔函数 $Y_\nu(x)$ （诺伊曼函数），现在可以解释内解和外解是什么了。将它们画出来后可以看到，第一类贝塞尔函数在原点是有限的，这对应于介质内的**“共振”或“驻波”**解，常用于封闭边界问题（如波导、圆盘模态），“内解”——源点处无奇异、适合描述内部场。第二类贝塞尔函数即诺伊曼函数在原点是发散的，它与 $J_\nu(x)$ 线性独立，反映了外部散射场，“外解”——外部辐射或奇异点源。

观察到

$\boxed{ \begin{aligned} J_{\nu}(x) &=\frac{1}{\pi}\int_0^{\pi} \cos(x\sin\tau-\nu \tau)\,d\tau \\ Y_\nu(x)&=\frac{1}{\pi}\int_0^{\pi} \sin(x\sin\tau-\nu \tau)\,d\tau \end{aligned} } \tag{70}$

如果我给第二类贝塞尔函数加个 $i$ ，再与第一类加起来，就能拼成一个复指数形式，这正好是第三类贝塞尔函数，即汉克尔函数的由来

$\boxed{ \begin{aligned} H_\nu^{(1)}(x)&=J_\nu(x)+iY_\nu(x)=\frac{1}{\pi}\int_0^\pi e^{i(x\sin\tau-\nu\tau)}d\tau\\ H_\nu^{(2)}(x)&=J_\nu(x)-iY_\nu(x)=\frac{1}{\pi}\int_0^\pi e^{-i(x\sin\tau-\nu\tau)}d\tau \end{aligned} }\tag{71}$

为什么要将贝塞尔函数的两个线性无关的解拼起来，凑成第三类解？因为拼起来的汉克尔函数刚好描述了辐射解！

由于积分式不够显式，我们可以研究远场处 $H_\nu^{(1),(2)}$ 的渐近行为，用一个渐近的显式对该解的情况进行描述，这里需要用到驻相法（stationary phase method）

驻相法的核心思想在于：当 $x$ 很大时，相位变化极快，积分里的指数在不同 $\tau$ 取值处的相位都不同，它们在复平面上不断绕圆旋转。如果相位在不同地方变化很快，指数项 $e^{ix\Phi(\tau)}$ 的方向到处不同，于是这些向量“首尾相连”后相互抵消，积分结果趋向零。只有那些区域里——相位变化慢的地方，这些指数向量的方向几乎相同，才会相加叠加强烈。

“相位变化慢”意味着导数小：

$\frac{d\Phi}{d\tau} = 0.\tag{72}$

这些点称为 相位平稳点 (stationary points)。在这些点附近，相位函数 $\Phi(\tau)$ 的一阶变化为零，指数项的相位在局部几乎“驻定”，所以相干叠加。

对于 $H_\nu^{(1)}$ ，相位函数 $\Phi(\tau)=x\sin\tau-\nu\tau$ ，令 $\Phi'(\tau)=x\cos\tau-\nu=0$ ，当 $x\to\infty$ ，找到驻点 $\tau=\pi/2$ 。

在驻点 $\tau=\pi/2$ 附近对 $\Phi(\tau)$ 展开

$\Phi(\tau)\approx x(1-\frac{t^2}{2})-\nu(\frac{\pi}{2}+t)=x-\frac{\nu\pi}{2}-\nu t-\frac{xt^2}{2}.\tag{73}$

将（73）带入（71），得到近似积分

$H_\nu^{(1)}(x)\approx\frac{1}{\pi}e^{i(x-\nu\pi/2)}\int_{0}^{\pi} e^{-i(\nu t+xt^2/2)}dt\approx\frac{1}{\pi}e^{i(x-\nu\pi/2)}\int_{-\infty}^{\infty} e^{-i(\nu t+xt^2/2)}dt\tag{74}$

这里注意，积分上下限之所以能扩展到无穷，是因为 $x\to\infty$ 时，相位变化非常快，远离驻点的区域相位剧烈震荡，彼此相消整体趋近于0，因此主要贡献早已包含在 $[0,\ \pi]$ 区间内，不论是扩展到更大的范围甚至 $[-\infty,\ \infty]$ 结果都是一样的。（74）右边是一个标准的复高斯积分，考虑到 $\nu t$ 作为线性相移相对高阶项 $xt^2/2$ 可忽略，因此近似：

$\int e^{-i(\nu t+\frac{x t^2}{2})}\,dt \approx e^{-i\pi/4}\sqrt{\frac{2\pi}{x}}.\tag{75}$

最后得到

$\boxed{ \begin{aligned} H_\nu^{(1)}(x\to\infty)&\sim\sqrt{\frac{2}{\pi x}}e^{i(x-\nu\pi/2-\pi/4)}\\ H_\nu^{(2)}(x\to\infty)&\sim\sqrt{\frac{2}{\pi x}}e^{-i(x-\nu\pi/2-\pi/4)} \end{aligned} }\tag{76}$

$x=kr$ ，由此可知

函数	相位	远场行为	物理意义
$H_\nu^{(1)}$	$e^{+ix}$	出射波（outgoing）	向外辐射
$H_\nu^{(2)}$	$e^{-ix}$	入射波（incoming）	向内传播

搞明白了二维情况下的贝塞尔方程，再来看三维下的球贝塞尔方程，这个时候坐标系从圆柱坐标变换到了球坐标

维度	波动方程径向项	典型解	远场衰减	物理意义
2D	$r^2R''+rR'+(k^2r^2-\nu^2)R=0$	$H_\nu^{(1,2)}(kr)$	$1/\sqrt{r}$	圆柱波（line source）
3D	$r^2R''+2rR'+(k^2r^2-l(l+1))R=0$	$h_l^{(1,2)}(kr)$	$1/r$	球面波（point source）

其实这两个方程都可以化成同一个形式，对于二维贝塞尔方程，令 $R(r)=\frac{u(r)}{\sqrt{r}}$ ，对三维贝塞尔方程，令 $R(r)=\frac{u(r)}{r}$ 得：

$\begin{aligned} u''r^2+[k^2r^2-(\nu-\frac{1}{2})(\nu+\frac{1}{2})]u&=0,\qquad\text{二维贝塞尔方程}\\ u''r^2+[k^2r^2-l(l+1)]u&=0,\qquad\text{三维贝塞尔方程} \end{aligned}\tag{77}$

因此

$\boxed{ \text{球贝塞尔方程} \;\Longleftrightarrow\; \text{普通贝塞尔方程的阶 } \nu=l+\tfrac{1}{2}. }\tag{78}$

因此，球贝塞尔方程的解

$j_l(x) = \sqrt{\frac{\pi}{2x}}\,J_{l+\tfrac{1}{2}}(x), \qquad n_l(x) = \sqrt{\frac{\pi}{2x}}\,Y_{l+\tfrac{1}{2}}(x).\tag{79}$

这里的系数是如何取的，为什么是 $\sqrt{\frac{\pi}{2x}}$ ？

答：这个系数的目的是为了让球贝塞尔函数自然满足三维球面波 $e^{\pm ix}/x$ 的幅度归一化。

根据前面驻相法求出的汉克尔函数 $H_\nu^{(1,2)}$ 的渐近式，根据汉克尔函数与第一第二类贝塞尔函数的关系

$\left\{ \begin{aligned} J_\nu(x)&=\frac{1}{2}(H_\nu^{(1)}+H_\nu^{(2)})\\ Y_\nu(x)&=\frac{1}{2i}(H_\nu^{(1)}-H_\nu^{(2)}) \end{aligned} \right.\tag{80}$

得到 $J_\nu(x)$ 与 $Y_\nu(x)$ 的渐近式

$\boxed{ \begin{aligned} J_\nu(x\to\infty)&\sim\sqrt{\frac{2}{\pi x}}\cos(x-\frac{\nu\pi}{2}-\frac{\pi}{4})\\ Y_\nu(x\to\infty)&\sim\sqrt{\frac{2}{\pi x}}\sin(x-\frac{\nu\pi}{2}-\frac{\pi}{4}) \end{aligned} }\tag{81}$

因此，球贝塞尔 $j_l(x)$ 和球诺伊曼 $y_l(x)$ 的渐近式

$\boxed{ \begin{aligned} j_l(x\to\infty)&\sim CJ_{l+\frac{1}{2}}=C\sqrt{\frac{2}{\pi x}}\cos(x-\frac{(l+1/2)\pi}{2}-\frac{\pi}{4})=C\sqrt{\frac{2}{\pi x}}\sin(x-\frac{l\pi}{2})\\ y_l(x\to\infty)&\sim CY_{l+\frac{1}{2}}=C\sqrt{\frac{2}{\pi x}}\sin(x-\frac{(l+1/2)\pi}{2}-\frac{\pi}{4})=-C\sqrt{\frac{2}{\pi x}}\cos(x-\frac{l\pi}{2}) \end{aligned} }\tag{82}$

如果我想使球汉克尔函数的解与三维球面波的解 $e^{\pm ix}/x$ 保持一致，即幅度归一化，此时

$h_l^{(1,2)}(x\to\infty)\sim \frac{e^{\pm i(x-(l+1)\pi/2)}}{x}\tag{83}$

观察系数，利用关系式 $h_l^{(1,2)}=j_l\pm iy_l$

$C\sqrt{\frac{2}{\pi x}}=\frac{1}{x}\qquad\Rightarrow\qquad C=\sqrt{\frac{\pi}{2x}}\tag{84}$

由此得到了（79）

组合成球汉克尔函数（传播解）：

$h_l^{(1)}=j_l+i y_l\sim \frac{e^{\,i(x-\frac{l\pi}{2})}}{ix}\quad(\text{出射波}), \qquad h_l^{(2)}=j_l-i y_l\sim \frac{ie^{-i(x-\frac{l\pi}{2})}}{x}\quad(\text{入射波}).\tag{85}$

格林函数可以展开为完整的模态和角向基

最终：

$G(\mathbf r,\mathbf r') = ik \sum_{l=0}^\infty \sum_{m=-l}^l j_l(kr_<)\,h_l^{(1)}(kr_>)\, Y_l^m(\theta,\phi) Y_l^{m*}(\theta',\phi'),\tag{86}$

它清晰展示：

球谐函数控制角度分布；
球贝塞尔/汉克尔函数控制径向传播；
格林函数是所有本征模态的叠加；
$k\to0$ 时退化为静电势格林函数 $1/4\pi r$ 。
$r_< = \min(r,r')$ ， $r_> = \max(r,r')$

这里需要关注的几点：

我们明明求出了所有的本征函数，包括径向的第一二类贝塞尔函数、第一二类汉克尔函数，以及角向的球谐函数，为什么这里只用到了第一类贝塞尔函数和第一类汉克尔函数？
为什么第一类贝塞尔函数和第一类汉克尔函数的复共轭没有乘上？求解本征函数体系不是还要乘上复共轭么？
为什么有个系数 $ik$ ，它是用来做归一化的么？

第一个问题，为什么只保留了第一类贝塞尔函数和第一类汉克尔函数。

格林函数是满足：

$(\nabla^2 + k^2)G(\mathbf r,\mathbf r') = -\delta(\mathbf r-\mathbf r'),\tag{87}$

并且符合一个特定的物理边界条件的解。

在波动问题中，我们常选的是Sommerfeld 辐射条件：

$\lim_{r\to\infty} r\!\left(\frac{\partial G}{\partial r} - ikG\right)=0.\tag{88}$

这要求——远处的波必须是向外传播的球面波（outgoing wave）。就好比已知 $\bold{r}=\bold{r'}$ 处有个源，该源往外辐射球面波，确定了是往外传播的物理边界条件。（87）的负号代表了能量往外流失。因此保留第一类汉克尔函数，因为它描述的就是往外传播的球面波解。

对 $\bold{r}<\bold{r'}$ ：要求在原点有限因为 δ 源之外没有物理奇点，所以在 $r=0$ 附近解必须有限。 $j_l(kr)$ 在 $r\to0$ 时 $\sim r^l$ ，而 $y_l\sim r^{-(l+1)}$ 发散，因此舍弃 $y_l$ 。

第二个问题，为什么没有乘上复共轭。

球谐函数乘上复共轭是因为在把三维格林函数写成分离变量形式时，利用了球谐函数的完备展开：

$\delta(\Omega - \Omega') =\sum_{l=0}^\infty\sum_{m=-l}^{l} Y_l^m(\Omega)\,Y_l^{m*}(\Omega'). \tag{89}$

这样做是为了保证：

$\int d\Omega' \,Y_l^{m*}(\Omega')\,Y_{l'}^{m'}(\Omega') =\delta_{ll'}\delta_{mm'},\tag{90}$

从而能正确地把格林函数的角向部分投影出来。所以共轭号是由正交关系强制出来的，是为了在角度空间上“分辨不同模态”。

径向函数没有加复共轭是因为径向部分是实数函数（或定向复解）对于实数 $k$ ：

$j_l(kr)$ 是实函数；
$h_l^{(1)}(kr)$ 是复函数，但我们固定选出一个方向的物理传播（出射波）。

如果把 $h_l^{(1)}$ 和 $h_l^{(2)}$ 都放进去，那就会出现它们的复共轭关系：

$h_l^{(2)}(x) = [h_l^{(1)}(x)]^{*}.\tag{92}$

但在定义格林函数时，我们只要一种传播方向（出射），因此只选 $h_l^{(1)}$ ，而不取它的共轭。

第三个问题，为什么系数是 $ik$ 。

这个 $ik$ 是对 $(\nabla^2 + k^2)G(\mathbf r,\mathbf r') = -\delta(\mathbf r-\mathbf r')$ 两侧积分后归一化得到的。具体推导过程如下

三维 δ 函数在球坐标下的分解是：

$\boxed{ \delta(\mathbf r - \mathbf r') = \frac{1}{r^2}\,\delta(r - r')\,\delta(\Omega - \Omega'), } \tag{92}$

分母出现 $r^2$ 是体积元 $d^3r=r^2\sin\theta drd\theta d\phi=r^2drd\Omega$ 导致的，假设有：

$\delta(\mathbf r-\mathbf r') = A(r)\,\delta(r-r')\,\delta(\Omega-\Omega').\tag{93}$

将它代入积分定义：

$\int\,\delta(\mathbf r-\mathbf r')\,d^3r = \int\,A(r)\,\delta(r-r')\,\delta(\Omega-\Omega')\,r^2\,dr\,d\Omega=1\qquad\Rightarrow\qquad A(r)=\frac{1}{r^2}.\tag{94}$

其中角向 δ 用球谐函数的正交展开表示：

$\delta(\Omega - \Omega') = \sum_{l=0}^\infty \sum_{m=-l}^l Y_l^m(\theta,\phi)\,Y_l^{m*}(\theta',\phi'). \tag{95}$

展开格林函数：

$\boxed{G(\mathbf r,\mathbf r') = \sum_{l=0}^{\infty}\sum_{m=-l}^{l} g_l(r,r')\,Y_l^m(\theta,\phi)\,Y_l^{m*}(\theta',\phi').} \tag{96}$

这里 $g_l(r,r')$ 就是径向格林函数。

代入（87）：

$(\nabla^2 + k^2)G = \sum_{l,m} \left[ \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2 \frac{\partial g_l}{\partial r}\right) -\frac{l(l+1)}{r^2}g_l +k^2 g_l \right] Y_l^m(\theta,\phi)\,Y_l^{m*}(\theta',\phi') = -\,\frac{1}{r^2}\sum_{l,m}\delta(r - r')Y_l^m Y_l^{m*}. \tag{97}$

利用球谐函数的完备性和正交性，可以逐个匹配系数 $(l,m)$ 对应的方程，得到每个模式的径向方程：

$\boxed{ \frac{1}{r^2}\frac{d}{dr}\!\left(r^2\frac{dg_l}{dr}\right) +\Big(k^2-\frac{l(l+1)}{r^2}\Big)g_l(r,r') = -\,\frac{\delta(r - r')}{r^2}. } \tag{98}$
把（98）两边乘以 $r^{2}$ ，再在 $[r'-\epsilon,r'+\epsilon]$ 上积分，得到标准的跳跃条件：

$\boxed{\;\Big[r^{2}\frac{dg_l}{dr}\Big]_{r'+0}^{r'-0}=-1.\;} \tag{99}$

由于齐次方程：

$\frac{1}{r^2}\frac{d}{dr}(r^2 g_l') + \left(k^2-\frac{l(l+1)}{r^2}\right)g_l = 0.\tag{100}$

它的两个独立解是：

区域解形式物理意义

$r < r'$ $j_l(kr)$ 内解（在原点有限）

$r > r'$ $h_l^{(1)}(kr)$ 外解（出射波）

因此写成拼接形式：

$g_l(r,r') = C\, j_l(k r_<)\,h_l^{(1)}(k r_>). \tag{101}$

其中 $r_<$ 表示较小的半径， $r_>$ 表示较大的半径。 $C$ 是尚未确定的归一化常数。

现在分别计算 $r>r'$ 与 $r<r'$ 时的导数：

$\frac{dg_l}{dr}= \begin{cases} C\,k\,j_l'(kr)\,h_l^{(1)}(kr'), & r<r',\\[4pt] C\,k\,j_l(kr')\,h_l^{(1)'}(kr), & r>r'. \end{cases}\tag{102}$

求它们在 $r=r'$ 的差：

$\Big[\frac{dg_l}{dr}\Big]_{r'+0}^{r'-0} = Ck\,[j_l(kr')h_l^{(1)'}(kr') - h_l^{(1)}(kr')j_l'(kr')]. \tag{103}$

这一组括号里的量，正是朗斯基行列式（Wronskian）。

什么是朗斯基行列式，对于任意两个函数 $f(x)$ 、 $g(x)$ ，它们的朗斯基定义为：

$W\{f,g\}(x) = f(x)g'(x) - f'(x)g(x).\tag{104}$

它衡量的是这两个函数是否线性无关。如果 $W\neq0$ ，表示它们确实是独立解。在贝塞尔方程里，朗斯基行列式不仅不为零，而且是个常数。
现在要求 $j_l(kr)$ 和 $h_l^{(1)}(kr)$ 的朗斯基行列式，考虑球贝塞尔函数的 $x\to\infty$ 渐近式

$\boxed{ \begin{aligned} j_l(x\to\infty)&\sim \frac{\sin(x-\frac{l\pi}{2})}{x}\\ y_l(x\to\infty)&\sim -\frac{\cos(x-\frac{l\pi}{2})}{x} \end{aligned} }\tag{105}$

利用 $h_l^{(1)}=j_l+i\,y_l$ ：

$\begin{aligned} W_x\{j_l,h_l^{(1)}\} &= j_l(h_l^{(1)})'-j_l' h_l^{(1)}\\ &= j_l(j_l'+i y_l')-j_l'(j_l+i y_l)\\ &= i\,(j_l y_l'-j_l' y_l) = i\,W_x\{j_l,y_l\}. \end{aligned}\tag{106}$

针对不同的 $l$ 进行讨论：

$W_x\{j_l,h_l^{(1)}\} = i(j_ly'_l-j_l'y_l) =i\cdot \begin{cases} \begin{aligned} &\frac{\sin x}{x}(-\frac{\cos x}{x})'-(\frac{\sin x}{x})'(-\frac{\cos x}{x}),\quad l=4k,\ k\in\mathbb Z\\ &(-\frac{\cos x}{x})(-\frac{\sin x}{x})'-(-\frac{\cos x}{x})'(-\frac{\sin x}{x}),\quad l=4k+1,\ k\in\mathbb Z\\ &(-\frac{\sin x}{x})(\frac{\cos x}{x})'-(-\frac{\sin x}{x})'\frac{\cos x}{x},\quad l=4k+2,\ k\in\mathbb Z\\ &\frac{\cos x}{x}(\frac{\sin x}{x})'-(\frac{\cos x}{x})'\frac{\sin x}{x},\quad l=4k+3,\ k\in\mathbb Z \end{aligned} \end{cases} =\frac{i}{x^2}. \tag{105}$

发现朗斯基行列式的结果都为 $i/x^2$

代入跳跃条件确定常数 $C$

$r'^2\Big[\frac{dg_l}{dr}\Big]_{r'+0}^{r'-0} = r'^2\,Ck\,W_x\{j_l,h_l^{(1)}\}\Big|_{x=kr'} = r'^2\,Ck\,\frac{i}{(kr')^2} = \frac{C\,i}{k} = -1 \quad\Longrightarrow\quad \boxed{C = ik}.\tag{106}$

因此得到径向格林函数

$\boxed{ g_l(r,r') = ik\,j_l(k r_<)\,h_l^{(1)}(k r_>). } \tag{107}$

原格林函数的展开式（86）得证。

区域	解形式	物理意义
$r < r'$	$j_l(kr)$	内解（在原点有限）
$r > r'$	$h_l^{(1)}(kr)$	外解（出射波）

附录一些证明

首先证明（28）

**方程对称性：**方程只含 $m^2$ ，所以 $P_l^{-m}$ 与 $P_l^m$ 必须成比例。

$P_l^{-m} = C_{lm} P_l^m.\tag{108}$

正交归一化条件：
为了让所有 $P_l^m(x)$ 在 $[-1,1]$ 上满足正交性：

$\int_{-1}^{1} P_l^m(x) P_{l'}^{m}(x)\,dx = \frac{2(l+m)!}{(2l+1)(l-m)!}\,\delta_{ll'},\tag{109}$

我们必须选择那个 $C_{lm}$ 使得

$P_l^{-m}\,P_l^{m}\text{ 的积分值正确对称。}$

代入 $P_l^{-m}=C_{lm} P_l^m$ ，就会发现：

$|C_{lm}|^2 = \frac{(l-m)!}{(l+m)!}.\tag{110}$

考虑相位一致性（由递推式符号确定），就得到：

$C_{lm} = (-1)^m \frac{(l-m)!}{(l+m)!}.\tag{111}$

既如此，正交归一化条件（109）又是如何得到的呢？

想要证明（109），就得先找到普通勒让德多项式的归一化条件

$\int_{-1}^{1} P_l(x) P_{l'}(x)\,dx = \frac{2}{2l+1}\,\delta_{ll'}.\tag{112}$

先根据 Rodrigues 罗德里格斯公式（普通勒让德多项式的解析式）推导（112），再根据普通和连带勒让德多项式的递推关系得到（109）

$\left\{ \begin{aligned} P_l(x)&=\frac{1}{2^l l!}\,\frac{d^l}{dx^l}\big[(x^2-1)^l\big],\quad \text{Rodrigues}\\ P_l^m(x) &= (-1)^m (1-x^2)^{m/2} \frac{d^m}{dx^m} P_l(x), \quad m\ge 0. \end{aligned} \right.\tag{113}$

定义

$I_l=\int_{-1}^{1}[P_l(x)]^2\,dx = \frac{1}{2^l l!}\int_{-1}^{1} P_l(x)\,\frac{d^l}{dx^l}\big[(x^2-1)^l\big]\,dx. \tag{114}$

把上式右端对第二因子连续作 $l$ 次分部积分（每次都把导数“移”到 $P_l$ 上）。因为在端点 $x=\pm1$ 时 $(x^2-1)^k=0$ （对 $k\ge1$ 及其低阶导数），所有边界项都为 0，于是

$\begin{aligned} I_l=\frac{(-1)^l}{2^l l!}\int_{-1}^{1} (x^2-1)^l\,\frac{d^l P_l}{dx^l}\,dx&=\frac{1}{2^l l!}\int_{-1}^{1} (1-x^2)^l\,\frac{d^l P_l}{dx^l}\,dx\\ &=\frac{1}{2^l l!}\int_{-1}^{1} (1-x^2)^l\,\frac{d^l}{dx^l}(\frac{(2l)!}{2^l(l!)^2}x^l+\dots)\,dx\\ &=\frac{1}{2^l l!}\int_{-1}^{1} (1-x^2)^l\,l!\frac{(2l)!}{2^l(l!)^2}\,dx\\ &=\frac{(2l)!}{2^{2l}(l!)^2}\int_{-1}^{1} (1-x^2)^ldx \end{aligned}\tag{115}$

现在问题变成了求 $\int_{-1}^{1} (1-x^2)^ldx$ 的积分了，这就是个纯数学问题，不会求的重新去学高数。

直接用 Beta 函数（ $\Beta$ 函数是第一类欧拉积分， $\Gamma$ 函数是第二类欧拉积分）：

$\int_{-1}^{1}(1-x^2)^l dx =2\int_{0}^{1}(1-x^2)^l dx\overset{t=x^2}{=} \int_{0}^{1} t^{-1/2}(1-t)^l dt = B\!\left(\tfrac12,l+1\right) =\dfrac{\Gamma(\tfrac12)\Gamma(l+1)}{\Gamma(l+\tfrac32)}.\tag{116}$

用递推 $\Gamma(z+1)=z\,\Gamma(z)$ 与 $\Gamma(\tfrac12)=\sqrt{\pi}$ ：

$\begin{aligned} \Gamma\!\left(l+\tfrac{3}{2}\right) &=\Big(l+\tfrac12\Big)\Big(l-\tfrac12\Big)\cdots\Big(\tfrac12\Big)\,\Gamma\!\left(\tfrac12\right)\\ &=\left[\prod_{k=0}^{l}\Big(k+\tfrac12\Big)\right]\sqrt{\pi} =\left[\frac{1\cdot3\cdot5\cdots(2l+1)}{2^{\,l+1}}\right]\sqrt{\pi}\\ &=\frac{(2l+1)!!}{2^{\,l+1}}\sqrt{\pi}. \end{aligned}\tag{117}$

于是

$\int_{-1}^{1}(1-x^2)^l dx=\frac{\sqrt{\pi}\,\Gamma(l+1)}{\dfrac{(2l+1)!!}{2^{\,l+1}}\sqrt{\pi}} =\frac{2^{\,l+1}\,\Gamma(l+1)}{(2l+1)!!}.\tag{118}$

再用 $\Gamma(l+1)=l!$ 得

$\boxed{\,\int_{-1}^{1}(1-x^2)^l dx=\frac{2^{\,l+1}\,l!}{(2l+1)!!}\,}.\tag{118}$

恒等式：

$(2l+1)!!=\frac{(2l+1)!}{2^{\,l}\,l!}.\tag{119}$

代入可得

$I_l=\frac{(2l)!}{2^{2l}(l!)^2}\int_{-1}^{1} (1-x^2)^ldx=\frac{(2l)!}{2^{2l}(l!)^2}\cdot\frac{2^{\,l+1}\,l!}{\dfrac{(2l+1)!}{2^{\,l}\,l!}} =\frac{2}{2l+1}\tag{120}$

（112）得证。

重新整理一下，我们想从左证明右

$\left\{ \begin{aligned} &\int_{-1}^{1} P_l(x) P_{l'}(x)\,dx = \frac{2}{2l+1}\,\delta_{ll'}.\\ &P_l^m(x) = (-1)^m (1-x^2)^{m/2} \frac{d^m}{dx^m} P_l(x), \quad m\ge 0\\ &(1-x^2)P_l'' - 2xP_l' + l(l+1)P_l = 0. \end{aligned} \right. \quad\Rightarrow \int_{-1}^{1} P_l^m(x) P_{l'}^{m}(x)\,dx = \frac{2(l+m)!}{(2l+1)(l-m)!}\,\delta_{ll'}\tag{121}$

直接带入好了，定义

$I_{m,l}\;=\;\int_{-1}^1 (1-x^2)^m\,[P_l^{(m)}(x)]^2\,dx, \qquad m=0,1,2,\dots\tag{122}$

这里 $P_l^{(m)}$ 指的是 $P_l$ 的 $m$ 阶导数 $\frac{d^m}{dx^m} P_l(x)$

对其中一个 $P^{(m)}$ 做分部积分（以下省略 $l$ 下标， $P^{(k)}\equiv P_l^{(k)}$ ），取

$u=(1-x^2)^m P^{(m)},\quad dv=P^{(m)}\,dx.\tag{123}$

则

$\begin{aligned} I_m=uv\Big|^1_{-1}-\int^1_{-1}vdu&=\Big[(1-x^2)^m P^{(m)}\,P^{(m-1)}\Big]_{-1}^{1}-\int_{-1}^1 \frac{d}{dx}\big[(1-x^2)^m P^{(m)}\big]\;P^{(m-1)}\,dx\\ &=0-\int^1_{-1}\big[m(-2x)(1-x^2)^{m-1} P^{(m)} + (1-x^2)^m P^{(m+1)}\big]P^{(m-1)}dx. \end{aligned}\tag{124}$

原普通勒让德方程：

$(1-x^2)P'' - 2xP' + l(l+1)P = 0.\tag{125}$

对 $x$ 连续微分 $m-1$ 次（ $m\ge1$ ），利用莱布尼茨法则，得到一个标准的递推恒等式：

$(1-x^2) P^{(m+1)} - 2m x\,P^{(m)} + \big[l(l+1)-m(m-1)\big]\,P^{(m-1)}=0. \tag{126}$

（126）乘以 $(1-x^2)^{m-1}P^{(m-1)}$ 并带入（124）积分

$\begin{aligned} I_m&=-\int^1_{-1}\Big[m(-2x)P^{(m)} + 2mxP^{(m)}-\big[l(l+1)-m(m-1)\big]\,P^{(m-1)}\Big]P^{(m-1)}(1-x^2)^{m-1}dx\\ &=\big[l(l+1)-m(m-1)\big]\int^1_{-1} (1-x^2)^{m-1}\big[P^{(m-1)}\big]^2dx=(l+m)(l-m+1)I_{m-1}. \end{aligned}\tag{127}$

反复应用（127）：

$I_m(l)=\prod_{k=1}^m (l+k)(l-k+1)\;I_0(l) =\frac{(l+m)!}{(l-m)!}\; \frac{2}{2l+1}.\tag{128}$

由此得到

$\int_{-1}^1 (1-x^2)^m\,[P_l^{(m)}]^2 dx =\frac{2(l+m)!}{(2l+1)(l-m)!}.\tag{129}$

再根据 $P_l^m(x)=(-1)^m(1-x^2)^{m/2}P_l^{(m)}(x)$

（129）等价于

$\int_{-1}^{1} P_l^m(x)\,P_l^m(x)\,dx =\frac{2(l+m)!}{(2l+1)(l-m)!}.\tag{130}$

（109）得证。

再来计算（29）的归一化系数

格林函数展开中用这种“未归一化”的 $Y_l^m$ 时，角向 δ 函数的展开式不是

$\delta(\Omega-\Omega')=\sum_{l,m}Y_l^m(\Omega)Y_l^{m*}(\Omega'),\tag{131}$

而是多一个系数：

$\delta(\Omega-\Omega')=\sum_{l,m}\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!} Y_l^m(\Omega)Y_l^{m*}(\Omega').\tag{132}$

因此做正交积分后

$N_{l,m}=(-1)^m\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}\tag{133}$

系数怎么来的？

定义球谐函数 $Y_l^m$ 的归一化函数 $S_l$

$S_l(\Omega,\Omega')=\sum_{m=-l}^{l}Y_l^m(\Omega)Y_l^{m*}(\Omega') .\tag{134}$

令 $\Omega'=\Omega$ ，由归一化条件

$\int_{S^2} |Y_l^m|^2\,d\Omega=1 = \int_{0}^{2\pi}\!\!\int_{-1}^{1} N_{lm}^2\,[P_l^m(x)]^2\,dx\,d\phi = 2\pi\,N_{lm}^2\!\int_{-1}^{1}[P_l^m(x)]^2 dx.\tag{135}$

带入（130）得到

$\int_{S^2} |Y_l^m|^2\,d\Omega =2\pi \cdot N_{l,m}^2 \cdot \frac{2(l+m)!}{(2l+1)(l-m)!} =\;1.\tag{136}$

由此解得（133）