前言
It seems clear that the present quantum mechanics is not in its final form. Some further changes will be needed, just about as drastic as the changes made in passing from Bohr’s orbit theory to quantum mechanics. Some day a new quantum mechanics, a relativistic one, will be discovered, in which we will not have these infinities occurring at all. It might very well be that the new quantum mechanics will have determinism in the way that Einstein wanted.
— Paul A.M. Dirac
1922 年 2 月 7 日晚上至 8 日清早，德国中部城市气温低寒，下雪，法兰克福大学里进行的施特恩－格拉赫实验(Stern—Gerlach experiment)首次成功测量到了电中性银原子束在非均匀磁场中的双重分裂现象，明确印证了微观尺度世界的原子角动量量子化特性本质！实验是由德国物理学家瓦尔特·格拉赫(Walther Gerlach，1889—1979)独自进行的，因为奥托·施特恩(Otto Stern，1888—1969)此时正短暂地在德国北部的罗斯托克大学任教(1921 年 9 月至 1922 年 12 月)，遇假期时才会回到法兰克福与格拉赫讨论数据，并进一步修正和改善实验设计。根据格拉赫当时的一位博士生回忆说，格拉赫是一位夜猫子，喜欢晚上 9 点进实验室，一直工作到第二天早晨。但采用这种工作时间的部分原因，也可能是出于当时的实验条件相对简陋，以及实验室空间狭小的限制；而且实验过程必须持续运转数个小时以上，才能在冷凝玻璃盘上积淀出足够多数量的银原子，以便可以清楚显影。在漫漫长夜里，如果仪器运转顺利，尤其是真空系统与装置的准直都维持正常平稳，格拉赫便会在一旁阅读文献、撰写论文，或是准备上课讲义，与今日的科学研究者的专注工作形式无异。
在玻尔的模型里，电子的轨道运动相当于一个闭合电路中的电流，而一个有电流通过的闭合电路的磁性作用就相当于一个磁壳，其磁矩等于$\mu=iA$，$i$是以安培为单位的电流，电子每个周期绕轨道一周，故$i=\frac{e}{\tau}$，$\tau$是周期，单位为秒。$A$是电路包围的面积，其可以进行如下计算
\[A=\int_{0}^{2\pi}\frac{1}{2}r\cdot rd\phi=\frac{1}{2}\int_{0}^{\tau}r^2\omega dt=\frac{1}{2m}\int_{0}^{\tau}mr^2\omega dt=\frac{p_{\phi}}{2m}\tau\]其中$p_{\phi}$叫做轨道角动量，结合之前的结论得
\[\mu=\frac{e}{2m}p_{\phi}\]此式表示轨道运动产生的磁矩与轨道角动量之间的数值关系。我们知道在玻尔的模型里，$p_{\phi}=\frac{n_{\phi}h}{2\pi}$，所以
\[\mu=n_{\phi}\frac{he}{4\pi m}=n_{\phi}\mu_B, \mu_B=0.92732\times10^{-23} A\cdot m^2\]$\mu_B$是轨道磁矩的最小单元，叫做玻尔磁子。为了验证这样的猜想，如上图的实验装置图，格拉赫在电炉里加热银原子至1000$^{\circ}$C((约为$10^3$K数量级)使其蒸发，通过数个狭缝后形成细束，再经过一个垂直方向上不均匀的磁场区域，最后撞到屏幕上，银原子通过的区域都被抽成了真空。如果银原子在磁场里的取向是任意的，那么屏幕上应该有连续变化的一片黑斑，而事实却是只有两条黑斑。我们对实验原理进行一些分析，从小孔射出的银原子的速率分布函数不满足麦克斯韦速率分布律($f(\nu)=4\pi(\frac{m}{2\pi kT})^{3/2}e^{-\frac{m\nu^2}{2kT}}\nu^2$)，而是满足碰壁原子的速率分布函数
\[F(\nu)=(\frac{m}{2 kT})^2e^{-\frac{m\nu^2}{2kT}}\nu^3\]利用$\frac{\mathrm{d} F(\nu)}{\mathrm{d} \nu}=0$可得最可几速率$m\nu^2=3kT$。此式表明当温度在1000$^{\circ}$C左右时，原子平均动能大约为 0.13eV，远远没有达到银原子从基态跃迁到最近的激发态的能量。根据温度热平衡下的原子按能级的玻尔兹曼分布
\[\frac{N_0}{N_1}=\rm{exp}[-\frac{1}{kT}(E_0-E_1)]\sim 10^{26}\]在基态的银原子的数目远远大于第一激发态的数目。射出后的银原子通过在极小尺度内的非均匀磁场，一个磁矩在不均匀磁场中感受到一个力为
\[f=\mu_z\frac{\partial B}{\partial z}=\mu\frac{\partial B}{\partial z}\cos{\beta}\]其中$\beta$是磁矩与z轴的夹角。原子在穿过长$L$的不均匀磁场中偏转的距离为
\[S=\frac{1}{2}at^2=\frac{1}{2}\frac{f}{m}(\frac{L}{\nu})^2=\frac{1}{2m}\frac{\partial B}{\partial z}(\frac{L}{\nu})^2\mu\cos{\beta}\]如果仅仅是$\mu$量子化($\mu=-\sqrt{l(l+1)}\mu_B$)，而$\cos{\beta}$可以是任意的话，那么与两条黑斑不能对应。由此可知，只能有两个$\beta$值，说明银原子在磁场中只能有两个取向。值得注意的是，屏幕上的两条黑斑是略有宽度的，不是很细的线条，这是由于银原子从热炉上发出时具有一个速度的分布。
随着波尔提出的旧量子论以及爱因斯坦提出的相对论，在 1916 年，索末菲考虑了相对论效应后的氢原子能量推广了玻尔的理论，提出了量子化通则以及电子椭圆轨道运动。按照索末菲理论，施特恩－格拉赫实验的银电子轨道量子数被证实为零,所 以预期原子束穿过不均匀磁场后，不管$l$取何值，由于磁场力的作用都将被分裂成($m_l=2l+1$)奇数条子原子束，修正后的波尔理论对于施特恩－格拉赫实验中间一条银原子条的消失感到十分不适，促使科学家的思想作出改变。比如，只有角动量为半整数(1/2)。
1925 年，年龄还不到 25 岁的两个荷兰莱顿大学的学生，乔治·乌伦贝克(George Eugene Uhlenbeck，1900—1988)和塞缪尔·古德斯米特(Samuel Abraham Goudsmit，1902—1978)根据一系列的实验事实(碱金属双线和反常塞曼效应)提出大胆假设：电子不是点电荷，它除了轨道角动量外，还要自旋运动，它具有固有的自旋角动量：$s=\frac{\hbar}{2}$，其在空间中任何一个方向的投影只能取两个自旋量子数$\pm\frac{1}{2}$。这样原始的，带有机械性质的运动，就如同地球自转一样，但经过计算在切向速度将大大超过光速(若电子看做带有电荷-e,半径为$10^{-14}$cm的小球)！正是这些概念上的困难，他们的假设一开始就遭到很多人的反对，包括泡利，以致他们想把已写好的文章收回，但是他们的导师却已把稿子寄出发表，并安慰说:”你们还年轻，有些荒唐，没关系”。后来的事实证明，电子自旋的概念是微观物理学中最重要的概念。
01
为了往更深的理论内容迈进，我们先回顾一下自由粒子的薛定谔方程：
\[E\psi(\mathbf{x},t)=i\hbar\frac{\partial }{\partial t}\psi(\mathbf{x},t)=H\psi(\mathbf{x},t)=\frac{\mathbf{p}}{2m}\psi(\mathbf{x},t)=-\frac{\hbar^2}{2m}\nabla^2\psi(\mathbf{x},t)\]其中$\mathbf{p}=-i\hbar\nabla$是粒子的动量算符，$E=i\frac{\partial }{\partial t}$是粒子的能量算符。这个方程有两个最基本的特征，第一，这是一个线性偏微分方程，因此满足量子力学最基本的原理，态的线性叠加原理。第二，这个方程里只出现时间的一阶导数，这使得我们可以定义一个正定的守恒概率$\int d^3\mathbf{x}| \psi(\mathbf{x},t) |^2$。薛定谔在发表这个方程时取的是非相对论形式，意味着它不满足爱因斯坦的狭义相对论。然而，考虑到相对论以后，自由粒子的能量与动量之间的关系应该是
\[E^2=\mathbf{p}^2c^2+m^2c^4\]因此，考虑到相对论，我们就应该将自由粒子哈密顿算符取成$H=\sqrt{\mathbf{p}^2c^2+m^2c^4}$。这里涉及到算符的开方，一个可行的定义是利用泰勒展开，但是这一方面会使得方程变得丑陋，另一方面，在一般性的相对论情形下，算符$\mathbf{p}^2c^2$与$m^2c^4$的大小可能相当，这样泰勒展开就没有良好的数学定义。看起来绕过这个困难的唯一办法是，放弃开方，直接使用相对论关系，也就是说，将自由粒子薛定谔方程推广成
\[-\hbar^2\frac{\partial^2 }{\partial t^2}\psi(\mathbf{x},t)=(\mathbf{p}^2c^2+m^2c^4)\psi(\mathbf{x},t)=(-\hbar^2c^2\nabla^2+m^2c^4)\psi(\mathbf{x},t)\]这样一来我们得到的就是一个关于时间的二阶微分方程——克莱因-戈尔登方程(Klein–Gordon equation)，由瑞典理论物理学家奥斯卡·克莱因(Oskar Benjamin Klein，1894—1977)和德国人沃尔特·戈尔登于二十世纪二三十年代分别独立推导得出的。移项后两边除以$c^2$和$\hbar^2$等
\[(\square+\frac{m^2c^2}{\hbar^2})\psi(\mathbf{x},t)=0\]其中，$\square=\frac{1}{c^2}\partial_t^2-\nabla^2=\partial^{\mu}\partial_{\mu}$是达朗贝尔算符(d’Alembert operator)。与薛定谔方程不同，克莱因-戈登方程在洛伦兹变换下是不变的，因此它是相对论量子力学方程的合格候选者之一。然而不幸的是，它还有诸多缺陷，这些缺陷使其正确的物理解释变得极其难以理解。该方程大体上产生了三个基本问题。首先，最明显的是，将其方程的平面波解$\psi(\mathbf{x},t)=\rm{exp}^{\frac{i}{\hbar}(\mathbf{p}\cdot\mathbf{x}-Et)}$代入，能量是$E=\pm\sqrt{\mathbf{p}^2c^2+m^2c^4}$，存在正负两个解。经典力学不存在量子跃迁特别是自发跃迁，于是总能量守恒就足以保证经典粒子运动的稳定性。而在量子理论范畴内，存在负无穷能级的负能解是难以自圆其说的。由于量子力学中有状态跃迁的概念，特别是有自发跃迁，我们无法解释为什么电子不能跳跃到负无穷能级解中，故无法确定稳定的基态，导致真空不稳定。第二个问题是克莱因-戈尔登方程的概率密度不正定，我们列出薛定谔方程的一般形式，
\[i\hbar\frac{\partial }{\partial t}\psi(\mathbf{x},t)=-\frac{\hbar^2}{2m}\nabla^2\psi(\mathbf{x},t)+V(\mathbf{x},t)\psi(\mathbf{x},t)\]其概率密度函数的定义是$\rho(\mathbf{x},t)=\psi^{\ast }(\mathbf{x},t)\psi(\mathbf{x},t)=| \psi |^2$，将上式对时间求导可以知道
\[\frac{\partial \rho}{\partial t}=\frac{\partial \psi^{\ast }}{\partial t}\psi+\psi^{\ast }\frac{\partial \psi}{\partial t}\]通过原来的式子，我们知道
\[\frac{\partial \psi}{\partial t}=-\frac{\hbar}{2mi}\nabla^2\psi+\frac{1}{i\hbar}V\psi\] \[\frac{\partial \psi^{\ast }}{\partial t}=\frac{\hbar}{2mi}\nabla^2\psi^{\ast }-\frac{1}{i\hbar}V\psi^{\ast }\]将这两个式子带入时间导函数得
\[\frac{\partial \rho}{\partial t}=-\frac{\hbar}{2mi}(\psi^{\ast }\nabla^2\psi-\psi\nabla^2\psi^{\ast })=-\frac{\hbar}{2mi}\nabla\cdot(\psi^{\ast }\nabla\psi-\psi\nabla\psi^{\ast })\]我们假设$\mathbf{J}=\frac{\hbar}{2mi}(\psi^{\ast }\nabla\psi-\psi\nabla\psi^{\ast })$，改写上式
\[\frac{\partial \rho}{\partial t}+\nabla\cdot\mathbf{J}=0\]$\mathbf{J}$的物理意义就是在单位时间内垂直流过单位面积的概率，叫做概率流密度。这个式子乃是概率守恒定律的微分形式。在克莱因-戈尔登方程里，我们也可以进行这样的变换。将原方程左边乘以$\psi^{\ast}$，减去其本身的共轭方程，我们得到
\[\psi^{\ast}(\square+\frac{m^2c^2}{\hbar^2})\psi-\psi(\square+\frac{m^2c^2}{\hbar^2})\psi^{\ast}=\psi^{\ast}\square\psi-\psi\square\psi^{\ast}=0\]已知$\square=\frac{1}{c^2}\partial^2_t-\nabla^2$，我们利用矢量分析的知识降低阶数，并让方程两边都乘以$\frac{\hbar}{2mi}$可得
\[\frac{\hbar}{2mi}\frac{1}{c^2}\frac{\partial }{\partial t}(\psi^{\ast}\frac{\partial \psi}{\partial t}-\psi\frac{\partial \psi^{\ast}}{\partial t})-\frac{\hbar}{2mi}\nabla\cdot(\psi^{\ast}\nabla\psi-\psi\nabla\psi^{\ast})=0\]整理成连续性方程(continuity equation)标准形式，
\[\frac{\partial }{\partial t}[\frac{i\hbar}{2mc^2}(\psi^{\ast}\frac{\partial \psi}{\partial t}-\psi\frac{\partial \psi^{\ast}}{\partial t})]+\nabla\cdot\frac{\hbar}{2mi}(\psi^{\ast}\nabla\psi-\psi\nabla\psi^{\ast})=0\]所以知道，其概率密度$\rho$和概率流密度$\mathbf{J}$为
\[\rho=\frac{i\hbar}{2mc^2}(\psi^{\ast }\frac{\partial \psi}{\partial t}-\psi\frac{\partial \psi^{\ast }}{\partial t})\] \[\mathbf{J}=\frac{\hbar}{2mi}(\psi^{\ast }\nabla\psi-\psi\nabla\psi^{\ast })\]我们知道，在这里，由于方程含时间的二阶导数，概率密度$\rho$表达式包含$\frac{\partial \psi}{\partial t}$和$\frac{\partial \psi^{\ast }}{\partial t}$，这和非相对论的不同。因为求解此二阶微分方程时，波函数和其导函数的初始条件可以彼此独立给定，造成概率密度$\rho$不一定是正定的，它就失去了作为概率密度解释的资格，并且一个正定的守恒概率，是整个量子力学理论的基石之一。最后一点，也是来自于原始方程里的对时间二阶导数，其初始条件$\frac{\partial \psi}{\partial t}|_0$无直接的物理意义，并且会给系统创造出一个额外的自由度。以上缺陷说明，一般而言，克莱因-戈登方程不宜作为相对论性粒子单粒子状态波函数方程。这导致 20 世纪 20 年代末和 30 年代初暂时放弃了克莱因-戈登方程。为了解决这个问题，一位新的物理学家来了。
02
为了克服克莱因-戈登方程的负概率密度问题，英国理论物理学家保罗·狄拉克(Paul Adrien Maurice Dirac，1902-1984)于 1928 年发现了另一版本满足相对论的方程——狄拉克方程。狄拉克意识到的这一条新的道路意味着，一方面相对论当然要满足，但在另一方面方程应该仅仅只含时间的一阶导数这个特征也应该保持。由于在相对论中时间和空间可以相互转换，因此这两点结合起来的逻辑推论就是，我们应该让方程对于动量算符(即对空间的导数)也是一阶的。而要保持一阶就必须将下面的开方开出来，
\[H=\sqrt{\mathbf{pc}^2+(mc^2)^2}\]这怎么可能呢？然而狄拉克却为此想出了“人类历史上最重要的开方”，只要你放弃对普通的数的坚持，转入到奇妙的对易世界里去。下面的一些论述来自B站Up主柚柚柚子235的视频。我们都知道，下列式子成立
\[\sqrt{a^2+b^2+2ab}=\pm(a+b)\]其实开方就是一个代数关系，那么对于我们也可以假设这个开方已经开出来了，但是不知道系数
\[\sqrt{a^2+b^2+c^2}=\sigma_1a+\sigma_2b+\sigma_3c\]只要它们的系数满足
\[\left\{\begin{matrix} \sigma^2_1= \sigma^2_2= \sigma^2_3= 1\\ \sigma_1\sigma_2+\sigma_2\sigma_1=0\\ \sigma_1\sigma_3+\sigma_3\sigma_1=0\\ \sigma_2\sigma_3+\sigma_3\sigma_2=0 \end{matrix}\right.\]如果$\sigma_1\sigma_2+\sigma_2\sigma_1=2\sigma_1\sigma_2=0$满足交换律，那么$\sigma_1$或$\sigma_2$为零，这四个等式不能同时被满足。问题可能就出现在这个交换律上，我们根据上式其实知道的是
\[\left\{\begin{matrix} \sigma_1\sigma_2=-\sigma_2\sigma_1\\ \sigma_1\sigma_3=-\sigma_3\sigma_1\\ \sigma_2\sigma_3=-\sigma_3\sigma_2 \end{matrix}\right.\]这个关系多么像三个单位矩阵之间的互相作用，这三个单位矩阵就是泡利矩阵和单位矩阵$\mathbf{I}$
\[\sigma_1=\begin{pmatrix} 0 & 1\\ 1 & 0\\ \end{pmatrix}, \; \sigma_2=\begin{pmatrix} 0 & -i\\ i & 0\\ \end{pmatrix}, \; \sigma_3=\begin{pmatrix} 1 & 0\\ 0 & -1\\ \end{pmatrix}, \; I=\begin{pmatrix} 1 & 0\\ 0 & 1\\ \end{pmatrix}\]那么对于
\[\sqrt{a^2+b^2+c^2+d^2}=\alpha_1a+\alpha_2b+\alpha_3c+\beta d\]同样的道理，
\[\alpha_1=\begin{pmatrix} 0 & \sigma_1\\ \sigma_1 & 0\\ \end{pmatrix}, \; \alpha_2=\begin{pmatrix} 0 & \sigma_2\\ \sigma_2 & 0\\ \end{pmatrix}, \; \alpha_3=\begin{pmatrix} 0 & \sigma_3\\ \sigma_3 & 0\\ \end{pmatrix}, \; \beta=\begin{pmatrix} I_{2\times 2} & 0\\ 0 & -I_{2\times 2}\\ \end{pmatrix}\]那么这些怎样被狄拉克应用的呢？我们重新审视自由粒子的相对论薛定谔方程
\[i\hbar\frac{\partial \psi}{\partial t}=H\psi=\sqrt{(\mathbf{p}c)^2+(mc^2)^2}\psi=[c(\alpha_1p_1+\alpha_2p_2+\alpha_3p_3)+\beta mc^2]\psi\]其中$p_i=-i\hbar\frac{\partial }{\partial x_i}, \; i=1,2,3$是动量算符(对空间的一阶导数)，将其写成向量形式
\[i\hbar\frac{\partial \psi}{\partial t}=(-i\hbar c\mathbf{\alpha}\cdot\mathbf{\nabla}+\beta mc^2)\psi\]这里的系数满足条件
\[\left\{\begin{matrix} (\alpha_i\alpha_j+\alpha_j\alpha_i)=2\delta_{ij}\\ \alpha_i\beta+\beta\alpha_i=0\\ (\alpha_i)^2=\beta^2=I \end{matrix}\right.\]如果我们标注$\gamma^0=\beta$和$\gamma^i=\beta\alpha_i$，再方程两边乘以矩阵$\beta$并除以$c$得
\[-i\hbar\beta\partial_t\psi-i\hbar\beta\alpha_i\partial_i\psi+mc\psi=(i\hbar\gamma^0\partial_0+i\hbar\gamma^i\partial_i-mc)\psi=(i\hbar\gamma^{\mu}\partial_{\mu}-mc)\psi=0\]若再利用费曼斜线标记和自然单位制$\hbar=c=1$，狄拉克方程的极简形式就是
\[(i\gamma^{\mu}\partial_{\mu}-m)\psi=(i\partial\!\!\!/-m)\psi=0\]可以验证，$\gamma^{\mu}$是四个反对易矩阵，其关系满足\(\{ \gamma^{\mu},\gamma^{\nu} \} =-2g^{\mu\nu}\)
\[g=\begin{pmatrix} -1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ \end{pmatrix}\]而对应的波函数解也是四维列向量，称作旋量(Spinor)，$\psi=(\psi_1,\psi_2,\psi_3,\psi_4)^{T}$，其共轭转置，也就是埃尔米特共轭矩阵为四维行向量$\psi^{\dagger }=(\psi^{\ast }_1,\psi^{\ast }_2,\psi^{\ast }_3,\psi^{\ast }_4)$。
让我们先检验一下狄拉克方程的概率密度是否正定。狄拉克方程左边乘以波函数的埃尔米特共轭$\psi^{\dagger }$
\[i\hbar\psi^{\dagger }\frac{\partial \psi}{\partial t}=\frac{\hbar c}{i}\psi^{\dagger }\gamma^i\partial_i\psi+mc^2\psi^{\dagger }\gamma^0\psi\]这个式子的复共轭是
\[-i\hbar\frac{\partial \psi^{\dagger }}{\partial t}\psi=-\frac{\hbar c}{i}(\partial_i\psi^{\dagger })\gamma^{i\dagger }\psi+mc^2\psi^{\dagger }\gamma^{0\dagger }\psi\]它们之间的差是
\[\frac{\partial }{\partial t}(\psi^{\dagger }\psi)=-c[(\partial_i\psi^{\dagger })\gamma^{i\dagger }\psi+\psi^{\dagger }\gamma^i\partial_i\psi]+\frac{imc^2}{\hbar}(\psi^{\dagger }\gamma^{0\dagger }\psi-\psi^{\dagger }\gamma^0\psi)\]我们根据之前的定义很容易知道，$\gamma^{i\dagger }_i=\gamma^i$，$\gamma^{0\dagger }=\gamma^0$都满足埃尔米特矩阵。所以单值的概率密度$\rho$和四维的概率流密度 $\mathbf{j}\equiv (j^0,j^1,j^2,j^3)^{T}$ 满足 $\frac{\partial }{\partial t}\rho+\nabla \mathbf{j}=0$
\[\rho\equiv \psi^{\dagger }\psi=\sum_{4}^{i=1}\psi_i^{\ast}\psi_{\mu}\] \[\mathbf{j}^i\equiv c\psi^{\dagger }\gamma^i\psi\]我们可以看到概率密度与克莱因-戈登方程不同，是正定的。
03
在用狄拉克方程求解氢原子前，我们先通过解决一些简单的场景来熟悉一下。静止粒子(particle at rest)的狄拉克方程因没有动能算符对能量的贡献，所以为
\[i\hbar\frac{\partial \psi}{\partial t}=m_ec^2\gamma^0\psi\] \[\frac{\partial }{\partial t}\begin{bmatrix} \psi_{\uparrow}^{(+)}\\ \psi_{\downarrow}^{(+)}\\ \psi_{\uparrow}^{(-)}\\ \psi_{\downarrow}^{(-)} \end{bmatrix}=-\frac{im_ec^2}{\hbar}\begin{bmatrix} 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\\ \end{bmatrix}\begin{bmatrix} \psi_{\uparrow}^{(+)}\\ \psi_{\downarrow}^{(+)}\\ \psi_{\uparrow}^{(-)}\\ \psi_{\downarrow}^{(-)} \end{bmatrix}\]这个方程因为$\gamma^0$是个$4\times 4$的方阵，所以有四个波函数解
\[\psi_{\uparrow}^{(+)}=e^{-\frac{im_ec^2}{\hbar}t}\begin{pmatrix} 1 \\ 0 \\ 0 \\ 0\end{pmatrix}, \; \psi_{\downarrow}^{(+)}=e^{-\frac{im_ec^2}{\hbar}t}\begin{pmatrix} 0 \\ 1 \\ 0 \\ 0\end{pmatrix}, \\ \psi_{\uparrow}^{(-)}=e^{\frac{im_ec^2}{\hbar}t}\begin{pmatrix} 0 \\ 0 \\ 1 \\ 0\end{pmatrix}, \; \psi_{\downarrow}^{(-)}=e^{\frac{im_ec^2}{\hbar}t}\begin{pmatrix} 0 \\ 0 \\ 0 \\ 1\end{pmatrix}\]其中$\psi_{\uparrow}^{(+)}$，$\psi_{\downarrow}^{(+)}$，$\psi_{\uparrow}^{(-)}$，$\psi_{\downarrow}^{(-)}$分别对应的是正负能量$m_ec^2$的解。回想一下，我们已经在克莱因-戈登方程的背景下遇到过负能量，狄拉克方程的构建是为了重现自由运动粒子的克莱因-戈登方程的结果，因此也适用于静止粒子。至于负能量，会在稍后的内容里进行解释。现在我们求解自由粒子的狄拉克方程，先提出解可能符合的一种平面波的拟设(Ansatz)，就像之前定态波函数一样分离时间变量，
\[\psi(\mathbf{x},t)=u(\mathbf{x})e^{\frac{i}{\hbar}(\mathbf{p}\cdot\mathbf{x}-Et)}\]我们知道，自由粒子的狄拉克方程可以写作
\[(i\hbar\gamma^{\mu}\partial_{\mu}-mc)\psi=[-i\hbar(\gamma^0\partial_0+\gamma^1\partial_1+\gamma^2\partial_2+\gamma^3\partial_3)+mcI_{4}]\psi=0\]因为$\gamma^0=\beta$和$\gamma^i=\beta\alpha_i$，把式子写成矩阵的形式
\[\left\{-i\hbar \left [ \begin{pmatrix} I_2 & 0_2 \\ 0_2 & -I_2 \\ \end{pmatrix}\frac{\partial }{\partial t}+\sum_{i=1}^{3}\begin{pmatrix} 0_2 & \sigma_i \\ -\sigma_i & 0_2 \\ \end{pmatrix}\frac{\partial }{\partial x_i}\right ]+mc\begin{pmatrix} I_2 & 0_2 \\ 0_2 & I_2 \\ \end{pmatrix}\right\}\psi=0\]根据动量算符$\mathbf{p}=-i\hbar\nabla$和泡利矩阵与动量的标积$\sigma\cdot\mathbf{p}=\sigma_1p_x+\sigma_2p_y+\sigma_3p_z$，
\[\left [-\begin{pmatrix} I_2 & 0_2 \\ 0_2 & -I_2 \\ \end{pmatrix}i\hbar\frac{\partial }{\partial t}+c\begin{pmatrix} 0_2 & \sigma\cdot\mathbf{p} \\ -\sigma\cdot\mathbf{p} & 0_2 \\ \end{pmatrix}+mc^2\begin{pmatrix} I_2 & 0_2 \\ 0_2 & I_2 \\ \end{pmatrix}\right ]\psi=0\]方程两边乘以$c$，和左乘矩阵$\gamma^0$可得
\[\left [-\begin{pmatrix} I_2 & 0_2 \\ 0_2 & I_2 \\ \end{pmatrix}i\hbar\frac{\partial }{\partial t}+c\begin{pmatrix} 0_2 & \sigma\cdot\mathbf{p} \\ \sigma\cdot\mathbf{p} & 0_2 \\ \end{pmatrix}+mc^2\begin{pmatrix} I_2 & 0_2 \\ 0_2 & -I_2 \\ \end{pmatrix}\right ]\psi=0\]再根据之前的拟设，加上$i\hbar\frac{\partial }{\partial t}e^{-i\frac{Et}{\hbar}}=Ee^{-i\frac{Et}{\hbar}}$，把上式的第一项计算后放到等式右边
\[\left [c\begin{pmatrix} 0_2 & \sigma\cdot\mathbf{p} \\ \sigma\cdot\mathbf{p} & 0_2 \\ \end{pmatrix}+\begin{pmatrix} mc^2I_2 & 0_2 \\ 0_2 & -mc^2I_2 \\ \end{pmatrix}\right ]u(\mathbf{x})=E\begin{pmatrix} I_2 & 0_2 \\ 0_2 & I_2 \\ \end{pmatrix}u(\mathbf{x})\]假设$u=(u^L,u^S)^{T}$，再简化变换后得到
\[\begin{pmatrix} mc^2-E & c\sigma\cdot\mathbf{p}\\ c\sigma\cdot\mathbf{p} & -mc^2-E\\ \end{pmatrix}\begin{pmatrix} u^L\\ u^S \end{pmatrix}=\begin{pmatrix} 0\\ 0 \end{pmatrix}\]计算这个矩阵可以得到关系
\[u^L=\left [ \frac{c\sigma\cdot\mathbf{p}}{E-mc^2} \right ]u^S, \; u^S=\left [ \frac{c\sigma\cdot\mathbf{p}}{E+mc^2} \right ]u^L\]接着
\[u^S=\left [ \frac{c\sigma\cdot\mathbf{p}}{E+mc^2} \right ]\left [ \frac{c\sigma\cdot\mathbf{p}}{E-mc^2} \right ]u^S\]所以下列关系必须成立
\[\left [ \frac{c\sigma\cdot\mathbf{p}}{E+mc^2} \right ]\left [ \frac{c\sigma\cdot\mathbf{p}}{E-mc^2} \right ]=1\]计算可以知道体系的能量为$E=\pm\sqrt{c^2\mathbf{p}^2+m^2c^4}$。我们知道
\[\sigma\cdot\mathbf{p}=\sigma_1p_x+\sigma_2p_y+\sigma_3p_z=\begin{pmatrix} p_z & p_x-ip_y \\ p_x+ip_y & -p_z \\ \end{pmatrix}\]我们如果取正能量$E=\sqrt{c^2\mathbf{p}^2+m^2c^4}$，并且电子自旋向上的方向，即
\[u^{L(+)}_{\uparrow}=\begin{pmatrix} 1\\ 0 \end{pmatrix}\]需要通过$u^L$和$u^S$的关系计算
\[u_{\uparrow}^{S(+)}=\left [ \frac{c\sigma\cdot\mathbf{p}}{E+mc^2} \right ]u^{L(+)}_{\uparrow}=\begin{pmatrix} \frac{cp_z}{E+mc^2} & \frac{c(p_x-ip_y)}{E+mc^2} \\ \frac{c(p_x+ip_y)}{E+mc^2} & \frac{-cp_z}{E+mc^2} \\ \end{pmatrix}\begin{pmatrix} 1\\ 0 \end{pmatrix}=\begin{pmatrix} \frac{cp_z}{E+mc^2}\\ \frac{c(p_x+ip_y)}{E+mc^2} \end{pmatrix}\]所以最后，该情况下对应波函数中一部分应该是
\[u^{(+)}_{\uparrow}=\mathcal{N}\begin{pmatrix} u^{L(+)}_{\uparrow}\\ u^{S(+)}_{\uparrow} \end{pmatrix}=\mathcal{N}\begin{pmatrix} 1\\ 0\\ \frac{cp_z}{E+mc^2}\\ \frac{c(p_x+ip_y)}{E+mc^2} \end{pmatrix}\]同理有正能量下，电子自旋方向向下的波函数中一部分为
\[u^{(+)}_{\downarrow}=\mathcal{N}\begin{pmatrix} u^{L(+)}_{\downarrow}\\ u^{S(+)}_{\downarrow} \end{pmatrix}=\mathcal{N}\begin{pmatrix} 0\\ 1\\ \frac{c(p_x-ip_y)}{E+mc^2}\\ \frac{-cp_z}{E+mc^2} \end{pmatrix}\]还有负能量下的电子自旋方向上和下的两个波函数中一部分为
\[u^{(-)}_{\uparrow}=\mathcal{N}\begin{pmatrix} u^{L(-)}_{\uparrow}\\ u^{S(-)}_{\uparrow} \end{pmatrix}=\mathcal{N}\begin{pmatrix} \frac{cp_z}{E-mc^2}\\ \frac{c(p_x+ip_y)}{E-mc^2}\\ 1\\ 0 \end{pmatrix}, \; u^{(-)}_{\downarrow}=\mathcal{N}\begin{pmatrix} u^{L(-)}_{\downarrow}\\ u^{S(-)}_{\downarrow} \end{pmatrix}=\mathcal{N}\begin{pmatrix} \frac{c(p_x-ip_y)}{E-mc^2}\\ \frac{-cp_z}{E-mc^2}\\ 0\\ 1 \end{pmatrix}\]归一化系数
\[\mathcal{N}=\sqrt{\frac{|E|+mc^2}{2|E|}}\]是从下列的一个例子中求得
\[u_{\uparrow}^{(+)\dagger } \cdot u^{(+)}_{\uparrow}=\mathcal{N}^2\left[ 1+\frac{c^2\mathbf{p}^2}{(E+mc^2)^2}\right]\overset{\underset{\mathrm{def}}{}}{=}1\]留下一个小作业：继续通过计算证明，对应负能量的波函数，其概率密度仍是正的，但是它的概率流密度流动的方向跟动量算符的方向相反。狄拉克方程仍然存在着负能量态，也就意味着，它的能量本征值没有下界，可以无限小。能量最低原理告诉我们，粒子总是更喜欢待在能量更低的态上。如果一个系统的能量本征值没有下界的话，那粒子就会持续不断地往更低的能态跃迁，这样的系统就不可能稳定存在。为了解释这个现象，狄拉克注意到电子是一个费米子，遵守泡利不相容原理，它的每一个能态上只能占据一个电子，狄拉克想象所有的负能态都已经被电子占据满了，即是说真空其实是一片填满负能电子的海洋，称为狄拉克海(Dirac sea)。由于负能电子都被占满了，所以真空中的电子就不再能够向负能态跃迁了，它只能待在正能态上，这样就解决了稳定性问题。
狄拉克进一步想象负能海中的某个电子被激发到了正能态。这时候真空中就产生了两种结果，一是多出了一个正能电子，二是负能海空出了一个未被占据的态，称之为空穴。负能海填满了电子当然带负电，但作为一种时空背景，我们无法感受到这个负能海的电荷。但是，如果负能海空出了一个态，那在我们的感受中，这个空穴就应该表现为正电荷！而且由于空穴是空出来的一个电子态，所以它在负能海中移动时其惯性应该和电子一样。即是说空穴也会表现出一个质量，并且和电子的质量一样。这就是说，负能海中的空穴在所有方面的表现都和电子一样，除了电荷相反，是正电荷。因此我们完全可以把负能海中的空穴当成是一种基本粒子，由于电荷相反，狄拉克称之为正电子(Positron)。所以激发负能海的结果就是在真空中产生了一对粒子，电子和正电子。反过来，假设我们把正电子和电子放在一起，则由于正电子是空穴，所以这时候正能的电子就会向下跃迁到这个空穴上，留下一个填满的负能海真空，同时辐射出光子。我们称这个过程为电子和正电子的湮灭。所以，电子和正电子可以成对产生，也可以成对湮灭！正电子除了电荷相反之外，其它表现都和电子相同，通常我们称之为电子的反粒子。正电子这种反粒子的存在是狄拉克的一个伟大预言，而它很快就在实验中观测到了！
04
终于我们来到了氢原子的狄拉克方程。
准备工作
我们先定义一个自旋为$1/2$的粒子总角动量为其轨道角动量和自旋角动量之和：
\[J=L+S=\mathbf{r}\times\mathbf{p}+\frac{\hbar}{2}\mathbf{\sigma}\]通过之前薛定谔方程的研究，我们知道轨道角动量的本征值满足$L^2=\sum_{i=1}^{3}L_i^2=l(l+1)\hbar^2, L_z=m_l\hbar$，其中$l=0,1,2,3,…;m_l=-l,-l+1,…,+l$。这里同理，继续定义自旋角动量的本征值满足$S^2=\sum_{i=1}^{3}S_i^2=s(s+1)\hbar^2=3/4\hbar^2, S_z=m_s\hbar$，其中$s=1/2;m_s=-s,-s+1,…,s=-1/2,1/2$。那么总角动量的本质值就满足$J^2=j(j+1)\hbar^2, J_z=L_z+S_z=m_j\hbar$，其中$j=l\pm s=l\pm 1/2=1/2,3/2,…;m_j=m_l+m_s=m_l\pm 1/2$。现在，我们把外电磁场下的狄拉克方程写出来
\[H\psi=\left [c\alpha\cdot(\mathbf{p}+\frac{q_e}{c}\mathbf{A})+\beta mc^2+q_eV \right]=i\hbar\frac{\partial }{\partial t}\psi=E\psi\]我们可以将其简化为电子受氢原子核中质子库伦力下的球对称中心势场问题，那么，$q_e=-e$，$\mathbf{A}=0$，$V=Ze/r$，跟上一章的静止粒子相似
\[H_{dirac}\psi=\left [c\alpha\cdot\mathbf{p}+\beta mc^2-\frac{Ze^2}{r} \right]=i\hbar\frac{\partial }{\partial t}\psi=E\psi\] \[\begin{pmatrix} mc^2-Ze^2/r-E & c\sigma\cdot\mathbf{p}\\ c\sigma\cdot\mathbf{p} & -mc^2-Ze^2/r-E\\ \end{pmatrix}\begin{pmatrix} u^L\\ u^S \end{pmatrix}=\begin{pmatrix} 0\\ 0 \end{pmatrix}\]很明显，为了求解，我们需要在球坐标系下表示$\sigma\cdot\mathbf{p}$，因为中心势场$Ze^2/r$已经是该坐标形式下的了，而且常数$mc^2$不受坐标变换影响。根据泡利矩阵的性质可知
\[\sigma\cdot\mathbf{p}=I_2\otimes(\sigma\cdot\mathbf{p})=(\frac{\sigma\cdot\mathbf{r}}{r})(\frac{\sigma\cdot\mathbf{r}}{r})(\sigma\cdot\mathbf{p})=\frac{1}{r^2}(\sigma\cdot\mathbf{r})(\sigma\cdot\mathbf{r})(\sigma\cdot\mathbf{p})\]有一个重要的矩阵变换公式 $(\sigma\cdot\mathbf{a})(\sigma\cdot\mathbf{b})=\mathbf{a}\cdot\mathbf{b}+i\sigma\cdot(\mathbf{a}\times\mathbf{b})$，所以
\[\sigma\cdot\mathbf{p}=\frac{1}{r^2}(\sigma\cdot\mathbf{r})(\sigma\cdot\mathbf{r})(\sigma\cdot\mathbf{p})=\frac{1}{r^2}(\sigma\cdot\mathbf{r})(\mathbf{r}\cdot\mathbf{p}+i\sigma\cdot(\mathbf{r}\times\mathbf{p}))=\frac{1}{r^2}(\sigma\cdot\mathbf{r})(\mathbf{r}\cdot\mathbf{p}+i(\sigma\cdot L))\]对于右边括号里第一项的算符，易知 $\mathbf{r}\cdot\mathbf{p}=\mathbf{r}\cdot(-i\hbar\nabla)=-i\hbar(\mathbf{r}\cdot\nabla)=-i\hbar r\frac{\partial }{\partial r}$只对波函数的径向部分起作用。因为$S=\hbar\sigma/2$，根据总角动量的定义知道，该算符的平方为
\[J^2=(L+S)^2=L^2+2(L\cdot S)+S^2=L^2+\hbar(\sigma\cdot L)+S^2\]于是，对于右边括号第二项的算符，其本征值为
\[\sigma\cdot L=\frac{1}{\hbar}(J^2-L^2-S^2)=\hbar\left[j(j+1)-l(l+1)-\frac{3}{4} \right]\]这里，我们假设另外一个重要的算符 $\mathcal{K}=-(\sigma\cdot L+\hbar I_2)$，其本征值可得
\[k=\mp(j+\frac{1}{2})\hbar\]对于左边括号里的算符，可以直接转换为球坐标系，并应用欧拉公式 $e^{(\pm i\phi)}=\cos{\phi}\pm i\sin{\phi}$
\[\sigma\cdot\mathbf{r}=\sigma_xx+\sigma_yy+\sigma_zz=\begin{pmatrix} z & x-iy \\ x+iy & -z\\ \end{pmatrix}=\begin{pmatrix} \cos{\theta} & \sin{\theta}e^{-i\phi} \\ \sin{\theta}e^{+i\phi} & -\cos{\theta}\\ \end{pmatrix}\]这里该算符对所有坐标系的表达式，可以理解是从笛卡尔坐标系变换到任意球坐标系的方法，而为了之后的计算简便，我们取$\theta=0$，即$z=r$来表达这两个坐标系的对应关系，于是
\[\sigma\cdot\mathbf{r}\overset{\underset{\theta=0}{}}{=}r\begin{pmatrix} 1 & 0 \\ 0 & -1\\ \end{pmatrix}\]最终形式的氢原子狄拉克方程写成
\[\begin{pmatrix} mc^2-Ze^2/r-E & \frac{c}{r}\frac{\sigma\cdot\mathbf{r}}{r}\left[ -i\hbar r\frac{\partial }{\partial r}+i(\sigma\cdot L)\right]\\ c\frac{1}{r}\frac{\sigma\cdot\mathbf{r}}{r}\left[ -i\hbar r\frac{\partial }{\partial r}+i(\sigma\cdot L)\right] & -mc^2-Ze^2/r-E\\ \end{pmatrix}\begin{pmatrix} u^L\\ u^S \end{pmatrix}=\begin{pmatrix} 0\\ 0 \end{pmatrix}\]或者是
\[\frac{c}{r}\frac{\sigma\cdot\mathbf{r}}{r}\left[ -i\hbar r\frac{\partial }{\partial r}+i(\sigma\cdot L)\right ]\begin{pmatrix} u^S\\ u^L \end{pmatrix}=\begin{pmatrix} E+Ze^2/r-mc^2 & 0\\ 0 & E+Ze^2/r+mc^2\\ \end{pmatrix}\begin{pmatrix} u^L\\ u^S \end{pmatrix}\]最后，我们还是提出一个符合径角分离的拟设解(Ansatz)，至此准备工作已经完成。
\[u =\begin{pmatrix} u^L\\ u^S \end{pmatrix} = \begin{pmatrix} \mathcal{F}_i(r)\mathcal{O}^{(+)}_{jm_j}(\theta,\phi)\\ i\mathcal{G}_i(r)\mathcal{O}^{(-)}_{jm_j}(\theta,\phi) \end{pmatrix}\]波函数的求解
在上一节的拟设解里，第二个角向波函数其实是球谐旋量(Spherical Spinors)，或者叫泡利旋量，它的表达式为
\[\mathcal{O}_{jm_j}(\theta,\phi)=\mathcal{C}Y_{lm_l}\rho_{m_s}=\mathcal{C}Y_{l(m_j-m_s)}\rho_{m_s}=\mathcal{C}_{\frac{1}{2}}Y_{l(m_j-\frac{1}{2})}\begin{pmatrix} 1 \\ 0 \end{pmatrix} +\mathcal{C}_{-\frac{1}{2}}Y_{l(m_j+\frac{1}{2})}\begin{pmatrix} 0 \\ 1 \end{pmatrix}=\begin{pmatrix} \mathcal{C}_{\frac{1}{2}}Y_{l(m_j-\frac{1}{2})}\\ \mathcal{C}_{-\frac{1}{2}}Y_{l(m_j+\frac{1}{2})} \end{pmatrix}\]其中$\mathcal{C}$是克莱布希－高登系数(Clebsch–Gordan coefficients)，用于描述两个角动量耦合时，它们本征函数的组合系数。当取 $j=l+1/2$时，即 $s=1/2$时，根据克莱布希－高登系数定义得
\[\mathcal{O}^{(+)}_{jm_j}=\begin{pmatrix} \sqrt{\frac{l+m_j+1/2}{2l+1}}Y_{l(m_j-\frac{1}{2})}\\ \sqrt{\frac{l-m_j+1/2}{2l+1}}Y_{l(m_j+\frac{1}{2})} \end{pmatrix}\]同理，当取 $j=l’-1/2$时，即 $s=-1/2$时，注意这里 $l’=l+1$当我们取相同的 $j$,因为需要保证 $l’>j$这一关系的成立。
\[\mathcal{O}^{(-)}_{jm_j}=\begin{pmatrix} -\sqrt{\frac{l'-m_j+1/2}{2l'+1}}Y_{l'(m_j-\frac{1}{2})}\\ \sqrt{\frac{l'+m_j+1/2}{2l'+1}}Y_{l'(m_j+\frac{1}{2})} \end{pmatrix}\]该本征函数满足下列关系
\[J^2\mathcal{O}^{(\pm)}_{jm_j}=j(j+1)\hbar^2\mathcal{O}^{(\pm)}_{jm_j}\] \[J_z\mathcal{O}^{(\pm)}_{jm_j}=m_j\hbar\mathcal{O}^{(\pm)}_{jm_j}\] \[L^2\mathcal{O}^{(\pm)}_{jm_j}=l(l+1)\hbar^2\mathcal{O}^{(\pm)}_{jm_j}\] \[S^2\mathcal{O}^{(\pm)}_{jm_j}=s(s+1)\hbar^2\mathcal{O}^{(\pm)}_{jm_j}\] \[\mathcal{K}\mathcal{O}^{(\pm)}_{jm_j}=k\hbar\mathcal{O}^{(\pm)}_{jm_j}=\mp (j+\frac{1}{2})\hbar\mathcal{O}^{(\pm)}_{jm_j}\]还有一点是
\[(\frac{\sigma\cdot\mathbf{r}}{r})\mathcal{O}^{(\pm)}_{jm_j}=-\mathcal{O}^{(\mp)}_{jm_j}\]现在让我们来尝试着解一下狄拉克方程
\[\frac{c}{r}\frac{\sigma\cdot\mathbf{r}}{r}\left[ -i\hbar r\frac{\partial }{\partial r}+i(\sigma\cdot L)\right ]\begin{pmatrix} i\mathcal{G}_i(r)\mathcal{O}^{(-)}_{jm_j}\\ \mathcal{F}_i(r)\mathcal{O}^{(+)}_{jm_j} \end{pmatrix}=\begin{pmatrix} E+Ze^2/r-mc^2 & 0\\ 0 & E+Ze^2/r+mc^2\\ \end{pmatrix}\begin{pmatrix} \mathcal{F}_i(r)\mathcal{O}^{(+)}_{jm_j}\\ i\mathcal{G}_i(r)\mathcal{O}^{(-)}_{jm_j} \end{pmatrix}\]使用关系式 $k\hbar\mathcal{O}{_{jm_j}^{(+)}}=-(\sigma\cdot L+\hbar)\mathcal{O}{^{(+)}_{jm_j}}, \, k\hbar\mathcal{O}{_{jm_j}^{(-)}}=(\sigma\cdot L+\hbar)\mathcal{O}{_{jm_j}^{(-)}}$，整理得
\[\frac{c}{r}\frac{\sigma\cdot\mathbf{r}}{r}\begin{pmatrix} \left[ \hbar r\frac{\partial }{\partial r}-(k-1)\hbar \right ]\mathcal{G}_i(r)\mathcal{O}^{(-)}_{jm_j}\\ \left[ -i\hbar r\frac{\partial }{\partial r}-i(k+1)\hbar \right ]\mathcal{F}_i(r)\mathcal{O}^{(+)}_{jm_j} \end{pmatrix}=\begin{pmatrix} E+Ze^2/r-mc^2 & 0\\ 0 & E+Ze^2/r+mc^2\\ \end{pmatrix}\begin{pmatrix} \mathcal{F}_i(r)\mathcal{O}^{(+)}_{jm_j}\\ i\mathcal{G}_i(r)\mathcal{O}^{(-)}_{jm_j} \end{pmatrix}\]再使用关系式 $(\frac{\sigma\cdot\mathbf{r}}{r})\mathcal{O}{_{jm_j}^{(\pm)}}=-\mathcal{O}{_{jm_j}^{(\mp)}}$，可得
\[\frac{\hbar c}{r}\begin{pmatrix} \left[ -r\frac{\partial }{\partial r}+(k-1) \right ]\mathcal{G}_i(r)\mathcal{O}^{(+)}_{jm_j}\\ \left[ i r\frac{\partial }{\partial r}+i(k+1) \right ]\mathcal{F}_i(r)\mathcal{O}^{(-)}_{jm_j} \end{pmatrix}=\begin{pmatrix} E+Ze^2/r-mc^2 & 0\\ 0 & E+Ze^2/r+mc^2\\ \end{pmatrix}\begin{pmatrix} \mathcal{F}_i(r)\mathcal{O}^{(+)}_{jm_j}\\ i\mathcal{G}_i(r)\mathcal{O}^{(-)}_{jm_j} \end{pmatrix}\]方程两边都去除函数 $\mathcal{O}{_{jm_j}^{(+)}}$ 和 $i\mathcal{O}{_{jm_j}^{(-)}}$ 可知
\[\hbar c\begin{pmatrix} -\frac{\partial \mathcal{G}_i(r)}{\partial r}+\frac{(k-1)}{r}\mathcal{G}_i(r) \\ \frac{\partial \mathcal{F}_i(r)}{\partial r}+\frac{(k+1)}{r} \mathcal{F}_i(r) \end{pmatrix}=\begin{pmatrix} \left[ E+Ze^2/r-mc^2 \right] \mathcal{F}_i(r)\\ \left[ E+Ze^2/r+mc^2 \right] \mathcal{G}_i(r) \end{pmatrix}\]我们至此得到了两个耦合的径向方程，使用假设 $F=r\mathcal{F}$ 和 $G=r\mathcal{G}$ 进一步简化这两个方程
\[\hbar c\begin{pmatrix} -\frac{1}{r}\frac{\partial G(r)}{\partial r}+\frac{G(r)}{r^2}+\frac{kG(r)}{r^2}-\frac{G(r)}{r^2} \\ \frac{1}{r}\frac{\partial F(r)}{\partial r}-\frac{F(r)}{r^2}+\frac{F(r)}{r^2}+\frac{kF(r)}{r^2} \end{pmatrix}=\begin{pmatrix} \left[ E+Ze^2/r-mc^2 \right] \frac{F(r)}{r}\\ \left[ E+Ze^2/r+mc^2 \right] \frac{G(r)}{r} \end{pmatrix}\] \[\begin{pmatrix} \frac{\partial G(r)}{\partial r}-\frac{kG(r)}{r}\\ \frac{\partial F(r)}{\partial r}+\frac{kF(r)}{r} \end{pmatrix}=\begin{pmatrix} \frac{mc^2-E-Ze^2/r}{\hbar c} F(r)\\ \frac{mc^2+E+Ze^2/r}{\hbar c} G(r) \end{pmatrix}\]我们之前提到过的精细结构常数 $\alpha=\frac{e^2}{\hbar c}=\frac{1}{137}$，再定义 $\lambda_1=\frac{mc^2+E}{\hbar c}$，$\lambda_2=\frac{mc^2-E}{\hbar c}$ 和 $\rho=\sqrt{\lambda_1\lambda_2}r$
\[\begin{pmatrix} \frac{\partial G(r)}{\partial r}-\frac{kG(r)}{r}\\ \frac{\partial F(r)}{\partial r}+\frac{kF(r)}{r} \end{pmatrix}=\begin{pmatrix} (\lambda_1-\frac{Z\alpha}{r}) F(r)\\ (\lambda_2+\frac{Z\alpha}{r}) G(r) \end{pmatrix}\]整理可以得到
\[\begin{pmatrix} (\frac{\partial }{\partial \rho}-\frac{k}{\rho})G(r)-(\sqrt{\frac{\lambda_2}{\lambda_1}}-\frac{Z\alpha}{\rho}) F(r)\\ (\frac{\partial }{\partial \rho}+\frac{k}{\rho})F(r)-(\sqrt{\frac{\lambda_1}{\lambda_2}}+\frac{Z\alpha}{\rho}) G(r) \end{pmatrix}=0\]这里，根据玻尔对应原理，即解的普遍形式应该包含非相对论的情况，比如，当 $\rho\to\infty$ 时，两个方程的渐近形式分别为 $\frac{\partial G}{\partial \rho}-\sqrt{\frac{\lambda_2}{\lambda_1}}F=0$ 和 $\frac{\partial F}{\partial \rho}-\sqrt{\frac{\lambda_2}{\lambda_1}}G=0$。于是，令 $F=Ae^{-\rho}$ 和 $G=Be^{-\rho}$ 代入得 $A=-\sqrt{\frac{\lambda_2}{\lambda_1}}B$，这两个解就渐近地成正比的关系，可采用带 $e^{-\rho}$ 乘子的广义幂级数试探解，所以我们假设
\[F=e^{-\rho}\rho^{s}\sum_{m=0}^{\infty}a_m\rho^m=e^{-\rho}\sum_{m=0}^{\infty}a_m\rho^{s+m}\] \[G=e^{-\rho}\rho^{s}\sum_{m=0}^{\infty}b_m\rho^m=e^{-\rho}\sum_{m=0}^{\infty}b_m\rho^{s+m}\]假设里的指数项会在 $\rho$ 很大时让整个解都趋近于零。我们需要去验证当 $\rho$ 趋近于零时，如果选择合适的 $a_0$，$b_0$ 和 $s$ 的值，会有合适的解。现在我们将上面的假设带入原来的式子
\[\begin{aligned} \sum_{m=0}^{\infty}\left(-a_m \rho^{s+m}+a_m(s+m) \rho^{s+m-1}-a_m k \rho^{s+m-1}-b_m \sqrt{\frac{\lambda_2}{\lambda_1}} \rho^{s+m}+b_m Z \alpha \rho^{s+m-1}\right) & =0 \\ \sum_{m=0}^{\infty}\left(-b_m \rho^{s+m}+b_m(s+m) \rho^{s+m-1}+b_m k \rho^{s+m-1}-a_m \sqrt{\frac{\lambda_1}{\lambda_2}} \rho^{s+m}-a_m Z \alpha \rho^{s+m-1} \right) & =0 \end{aligned}\]令上列广义幂级数展开式中 $s+m$ 次幂项的系数和为零，得到系数的递推关系
\[\begin{aligned} \left(-a_m+a_{m+1}(s+m+1)-a_{m+1} k-b_m \sqrt{\frac{\lambda_2}{\lambda_1}}+b_{m+1} Z \alpha\right) & =0 \\ \left(-b_m+b_{m+1}(s+m+1)+b_{m+1} k-a_m \sqrt{\frac{\lambda_1}{\lambda_2}}-a_{m+1} Z \alpha\right) & =0 \\ -a_m+(s+m+1-k) a_{m+1}-\sqrt{\frac{\lambda_2}{\lambda_1}} b_m+Z \alpha b_{m+1} & =0 \\ -b_m+(s+m+1+k) b_{m+1}-\sqrt{\frac{\lambda_1}{\lambda_2}} a_m-Z \alpha a_{m+1} & =0 \end{aligned}\]上面的式子就是 $a_m$ 和 $b_m$ 的递归关系。因为我们的解最低次幂项是 $\rho^{s}$，需要确保没有比它更低的幂级数，所以我们看当 $m=-1$ 时的递归关系，并让 $a_m=b_m=0$，
\[\begin{array}{r} (s-k) a_0+Z \alpha b_0=0 \\ -Z \alpha a_0+(s+k) b_0=0 \\ \left(\begin{array}{cc} (s-k) & Z \alpha \\ -Z \alpha & (s+k) \end{array}\right)\left(\begin{array}{l} a_0 \\ b_0 \end{array}\right)=0 \\ s^2-k^2+Z^2 \alpha^2=0 \\ s^2=k^2-Z^2 \alpha^2 \\ s=\pm \sqrt{k^2-Z^2 \alpha^2} \end{array}\]注意到当$k$是一个非零整数时，$Z^2 \alpha^2$是一个很小的值。所以我们取正解 $s=+\sqrt{k^2-Z^2 \alpha^2}$。我们还注意到，原来假设的幂级数解不能无限长，必须在某个 $n_r$ 值处截断，这就意味着，当系数在 $n_r$ 处为零后，其后所有系数均为零。将 $a_{n_r}=b_{n_r}=0$ 代入递推关系得 $a_{n_r}=-\sqrt{\lambda_2/\lambda_1}b_{n_r}$，与我们之前的结论相符合。最后一步，我们将 $m=n_r-1$ 代入递归关系式，这里相当于我们承认 $a_{m+1}=b_{m+1}=0$
\[-a_{n_r-1}+\left(s+n_r-k\right) a_{n_r}-\sqrt{\frac{\lambda_2}{\lambda_1}} b_{n_r-1}+Z \alpha b_{n_r}=0 \\ -\sqrt{\frac{\lambda_1}{\lambda_2}} a_{n_r-1}-Z \alpha a_{n_r}-b_{n_r-1}+\left(s+n_r+k\right) b_{n_r}=0\]让第一个式子方程两边乘以 $\sqrt{\lambda_1/\lambda_2}$，然后减去第二个式子得
\[\begin{aligned} & \left(\left(s+n_r-k\right) \sqrt{\frac{\lambda_1}{\lambda_2}}+Z \alpha\right) a_{n_r}+\left(Z \alpha \sqrt{\frac{\lambda_1}{\lambda_2}}-\left(s+n_r+k\right)\right) b_{n_r}=0 \\ & \left(-\left(s+n_r-k\right) \sqrt{\frac{\lambda_1}{\lambda_2}}-Z \alpha\right) \sqrt{\frac{\lambda_2}{\lambda_1}} b_{n_r}+\left(Z \alpha \sqrt{\frac{\lambda_1}{\lambda_2}}-\left(s+n_r+k\right)\right) b_{n_r}=0 \\ & -\left(s+n_r-k\right)-Z \alpha \sqrt{\frac{\lambda_2}{\lambda_1}}+Z \alpha \sqrt{\frac{\lambda_1}{\lambda_2}}-\left(s+n_r+k\right)=0 \\ & -\left(s+n_r-k\right) \sqrt{\lambda_1 \lambda_2}-Z \alpha \lambda_2+Z \alpha \lambda_1-\left(s+n_r+k\right) \sqrt{\lambda_1 \lambda_2}=0 \\ & -2\left(s+n_r\right) \sqrt{\lambda_1 \lambda_2}+Z \alpha\left(\lambda_1-\lambda_2\right)=0 \\ & 2\left(s+n_r\right) \sqrt{\lambda_1 \lambda_2}=Z \alpha\left(\lambda_1-\lambda_2\right) \\ & 2\left(s+n_r\right) \sqrt{m^2 c^4-E^2}=2 Z \alpha E \\ & \left(s+n_r\right) \sqrt{m^2 c^4-E^2}=Z \alpha E \\ & \left(s+n_r\right)^2\left(m^2 c^4-E^2\right)=Z^2 \alpha^2 E^2 \\ & \left(s+n_r\right)^2\left(m^2 c^4\right)=\left(Z^2 \alpha^2+\left(s+n_r\right)^2\right) E^2 \\ & \frac{\left(s+n_r\right)^2}{\left(\left(s+n_r\right)^2+Z^2 \alpha^2\right)}\left(m^2 c^4\right)=E^2 \\ & E^2=\frac{m^2 c^4}{\left(1+\frac{Z^2 \alpha^2}{\left(s+n_r\right)^2}\right)} \\ & E=\frac{m c^2}{\sqrt{1+\frac{Z^2 \alpha^2}{\left(n_r+s\right)^2}}} \\ & E=\frac{m c^2}{\sqrt{1+\frac{Z^2 \alpha^2}{\left(n_r+\sqrt{\left.k^2-Z^2 \alpha^2\right)^2}\right.}}} \\ & \boxed{E=\frac{m c^2}{\sqrt{1+\frac{Z^2 \alpha^2}{\left(n_r+\sqrt{\left(j+\frac{1}{2}\right)^2-Z^2 \alpha^2}\right)^2}}}} \\ & \end{aligned}\]这里，我们定义得主量子数 $n=n_r+ | k |=n_r+j+1/2$，所以氢原子光谱的精细结构如下图
总结
我们通过结合量子力学与相对论，了解了狄拉克方程和其对氢原子光谱精细结构的计算。然而，这个理论也有它自身的缺陷，比如，还有氢原子的超精细结构，一个最明显的例子就是兰姆位移，即在狄拉克方程里 $2S_{1/2}$ 与 $2P_{1/2}$ 是简并的，只是 $k=\pm 1$的符号区别(算能级是取绝对值)，但经实验发现它们之间还是有区别的。然而，对于狄拉克方程最重要的问题是，它本身是一个“单粒子”的方程，无法描述粒子的产生与湮灭。于是，新的理论——量子场论——正在这样的追寻中悄然而生。
Reference
- 《Advanced Quantum Mechanics》 F. Schwabl, Springer-Verlag, Berlin, 3rd ed., 2005
- 《Relativistic Quantum Chemistry 》Markus Reiher Alexander Wolf, Second Edition
- Course notes of Physics 130A in University of California San Diego
- 《Modern Quantum Mechanics》 J. J. Sakurai, Cambridge University Press, 3rd ed., 2020
- 《Quantum Mechanics of One- and Two-Electron Atoms》 Hans A. Bethe , Edwin E. Salpeter, Springer US, 1977