计算专题4 水星进动

前言

还有几天开学，想想我还是再写一个专题出来罢. 至于能不能写完还另说，但是总得开始.

这个专题倾向于对涉及较多知识的问题的研究，也算是一种对专题的更新想法吧. 同时我现在的写作更转向于写完一整篇再上传，这样可能效率会更高吧.

接下来我们开始计算. 我将考虑三个效应：经典附加力场（大行星影响）、狭义相对论、广义相对论.

Part 1经典附加力场（大行星影响）

这里最主要的就是木星的影响，我们使用的方法为木星星环等效. 如图所示，将木星的影响视作一个位置在木星轨道上的大圆环所施加的影响，大圆环半径即为木星轨道半径$R_J$，质量即为木星质量$M_J$，圆环质量分布是均匀的，水星轨道半径为$r$.

我们标记$\angle ABC=\theta$，$\angle ADC=\alpha$，考虑圆心角$d\theta$的一段弧对$D$点产生的引力场：

$$\begin{aligned} dg'=\frac{M_J}{2\pi}d\theta\cdot G\cdot\frac{1}{R_J^2+r^2-2R_Jr\cos\theta}\\\\ \end{aligned}$$

整个环产生的引力场：

$$\begin{aligned} \bold{g}'&=\hat{r}\int_0^{2\pi}dg\cdot\cos\alpha\\\\&=\hat{r}\int_0^{2\pi}\frac{GM_J}{2\pi}\frac{R_J\cos\theta-r}{(R_J^2+r^2-2R_Jr\cos\theta)^{3/2}}d\theta\\\\ \end{aligned}$$

这里涉及到$\cos\alpha$的计算，应用正弦定理和余弦定理：

$$\begin{aligned} \cos\alpha&=\sqrt{1-\sin^2\alpha}\\\\ \sin\alpha&=\sin\angle ADB=\frac{AB}{AD}\sin\theta\\\\ &=\frac{R_J\sin\theta}{\sqrt{R_J^2+r^2-2R_Jr\cos\theta}}\\\\ \Longrightarrow\quad\cos\alpha&=\sqrt{\frac{R_J^2\cos^2\theta+r^2-2R_Jr\cos\theta}{R_J^2+r^2-2R_Jr\cos\theta}}\\\\&=\frac{R_J\cos\theta-r}{\sqrt{R_J^2+r^2-2R_Jr\cos\theta}}\\\\ \end{aligned}$$

回到上面那个积分式，显然无法正常积分出结果. 所以我们考虑做多极展开，这也是处理平方反比力问题的常用手段（比如电偶极和电四极的问题中常用）. 认为$r/R_J$为一阶小量，展开到一阶项得到：

$$\begin{aligned} \bold{g}'&=\hat{r}\cdot\frac{GM_JR_J}{2\pi R_J^3}\int_0^{2\pi}\frac{\cos\theta-\frac{r}{R_J}}{(1-\frac{2r\cos\theta}{R_J}+\frac{r^2}{R^2})^{3/2}}d\theta\\\\ &\approx\hat{r}\cdot\frac{GM_JR_J}{2\pi R_J^3}\int_0^{2\pi}\cos\theta(1-\frac{r}{R_J\cos\theta}+\frac{3r\cos\theta}{R_J})d\theta\\\\ &=\hat{r}\cdot\frac{GM_J}{2\pi R_J^2}\int_0^{2\pi}(3\cos^2\theta-1)\frac{r}{R_J}d\theta+0\\\\ &=\frac{GM_Jr}{2R_J^3}\hat{r}\\\\ \end{aligned}$$

可见，星环等效产生的附加力为线性力，方向向外.

这个附加的引力场也可以通过引力场的Gauss定理来求，只需在环心处围出一个圆柱形的Gauss面即可.

接下来计算进动. 动力学方程可以写为：

$$\begin{aligned} \ddot{r}-r\dot{\theta}^2=-\frac{GM_S}{r^2}+\frac{GM_Jr}{2R_J^3}\,,\quad r^2\dot{\theta}=\frac{L}{m}\\\\ \end{aligned}$$

其中，$M_S$为太阳质量，$m$为水星质量，$L$是水星绕行的轨道角动量. 显然，方程等号右边的两个力有明显的大小差异，星环附加力场的影响远小于太阳引力场，所以使用摄动法求解. 设$r=a+\delta$，其中$\delta\ll a$，为摄动项，则

$$\begin{aligned} \ddot{a}+\ddot{\delta}&\approx-\frac{GM_S}{a^2}(1-\frac{2\delta}{a})+\frac{GM_Ja}{2R_J^3}(1+\frac{\delta}{a})+\frac{L^2}{m^2a^3}(1-\frac{3\delta}{a})\\\\ &=-\frac{GM_S}{a^2}+\frac{GM_Ja}{2R_J^3}+\frac{L^2}{m^2a^3}+(\frac{2GM_S}{a^3}+\frac{GM_J}{2R_J^3}-\frac{3L^2}{m^2a^4})\delta\\\\ \end{aligned}$$

注意到在零阶情况下，

$$\begin{aligned} \ddot{a}&=-\frac{GM_S}{a^2}+\frac{GM_Ja}{2R_J^3}+\frac{L^2}{m^2a^3}=0\\\\ \Longrightarrow&\quad\frac{L^2}{m^2a^3}=\frac{GM_S}{a^2}-\frac{GM_Ja}{2R_J^3}\\\\ \end{aligned}$$

可能产生问题的是，为什么这里会有$GM_Ja/2R_J^3$这一一阶项，毕竟这个方程是零阶的. 这是因为，$a$与$\delta$的区别来自于是否是微扰项，而非其二阶导数是否是小量，所以只要是稳定运动的成分统统应该归于$\ddot{a}$. 是为说明.

代入上述一阶方程可知，

$$\begin{aligned} \ddot{\delta}&=(\frac{2GM_S}{a^3}+\frac{GM_J}{2R_J^3}-\frac{3GM_S}{a^3}+\frac{3GM_J}{2R_J^3})\delta\\\\ &=(-1+\frac{2M_Ja^3}{M_SR_J^3})\frac{GM_S}{a^3}\delta\\\\ \end{aligned}$$

故，径向振动角频率

$$\begin{aligned} \omega_r&=\sqrt{\frac{GM_S}{a^3}(1-\frac{2M_Ja^3}{M_SR_J^3})}\approx\sqrt{\frac{GM_S}{a^3}}(1-\frac{M_Ja^3}{M_SR_J^3})\\\\ \end{aligned}$$

而注意到圆周运动的角频率为

$$\begin{aligned} \omega_0=\frac{L}{ma^2}=\sqrt{\frac{GM_S}{a^3}(1-\frac{M_Ja^3}{2M_SR_J^3})}\approx\sqrt{\frac{GM_S}{a^3}}(1-\frac{M_Ja^3}{4M_SR_J^3})\\\\ \end{aligned}$$

所以进动角频率为

$$\begin{aligned} \Omega=|\omega_r-\omega_0|=\frac{3M_Ja^3}{4M_SR_J^3}\sqrt{\frac{GM_S}{a^3}}\\\\ \end{aligned}$$

同时，轨道进动的方向与水星轨道运行的方向相同，因为$\omega_r<\omega_0$. 可以这样理解：当水星轨道运行一周后，径向简谐振动还没有达到平衡位置，还要再往前走一段距离才能到达平衡位置，这时轨道的近（远）日点就往轨道运行方向旋转了一个角度，也即进动.

周期进动角为

$$\begin{aligned} \Delta\theta&=\Omega\cdot T=\frac{3\pi M_Ja^3}{2M_SR_J^3}\\\\ \end{aligned}$$

Part 2 狭义相对论效应

在研究狭义相对论效应时，我们采用两种不同的方法：狭义相对论性Binet方程法和LRL矢量法.

方法1 狭义相对论性Binet方程

先做推导. 考虑狭义相对论性的水星能量，并做$v^2/c^2$的$O(2)$阶近似（当然也可以说是$v/c$的$O(4)$阶），得到

$$\begin{aligned} E_k&=\sqrt{p^2c^2+m^2c^4}-mc^2=mc^2\sqrt{\frac{v^2}{c^2(1-\frac{v^2}{c^2})}+1}-mc^2\\\\ &=mc^2\sqrt{\frac{1}{1-\frac{v^2}{c^2}}}-mc^2\\\\ &\approx\frac{1}{2}mv^2 +\frac{3}{8}mc^2(\frac{v^2}{c^2})^2\\\\ \end{aligned}$$

同时，在狭义相对论下不计引力的修正，所以势能不变，为$E_p=-\frac{GM_Sm}{r}$. 所以总能量守恒写为

$$\begin{aligned} E_0&=-\frac{GM_Sm}{r}+\frac{1}{2}mv^2+\frac{3}{8}mc^2(\frac{v^2}{c^2})^2\\\\ \Longrightarrow&\quad\frac{v^2}{c^2}\approx2\frac{E_0+\frac{GM_Sm}{r}} {mc^2}\\\\ \end{aligned}$$

上面已经推出$O(1)$下$v^2/c^2$的结果，这是为了在之后处理$(v^2/c^2)^2$时更为方便. 注意到角动量守恒，则

$$\begin{aligned} L&=\gamma mr^2\dot{\theta}=\text{const.}\\\\ \Longrightarrow\quad\dot{\theta}=\frac{L}{\gamma mr^2}&\,,\quad\dot{r}=\dot{\theta}\frac{dr}{d\theta}=\frac{L}{\gamma mr^2}\frac{dr}{d\theta}\\\\ \end{aligned}$$

令$u=1/r$，有

$$\begin{aligned} \dot{\theta}=\frac{Lu^2}{\gamma m}\,,\quad\dot{r}=-\frac{L}{\gamma m}\frac{du}{d\theta}\\\\ \end{aligned}$$

能量守恒式写为

$$\begin{aligned} \frac{1}{2}mv^2+\frac{3}{8}mc^2(\frac{v^2}{c^2})^2-\frac{GM_Sm}{r}&=E_0\\\\ \Longrightarrow\quad\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2)+\frac{3}{8}mc^2(2\frac{E_0+\frac{GM_Sm}{r}}{mc^2})^2&=E_0+\frac{GM_Sm}{r}\\\\ \Longrightarrow\quad\frac{L^2}{2m}[(\frac{du}{d\theta})^2+u^2](1-\frac{v^2}{c^2})=(E_0+\frac{GM_Sm}{r})&(1-\frac{3}{2}\frac{E_0+\frac{GM_Sm}{r}}{mc^2})\\\\ \Longrightarrow\quad\frac{L^2}{2m}[(\frac{du}{d\theta})^2+u^2]=(E_0+\frac{GM_Sm}{r})(1+&\frac{1}{2}\frac{E_0+\frac{GM_Sm}{r}}{mc^2})\\\\ \end{aligned}$$

值得注意的是，$mc^2$是一个相当大的量，使得二阶小量$(v^2/c^2)^2$在与它相乘时也成为了一个一阶小量，这就是之前Taylor展开后要保留这一项的原因.

接下来是对上述方程进行求导. 我们先全部化为$u$与$\theta$的方程，再进行这一步操作.

$$\begin{aligned} \frac{L^2}{2m}[(\frac{du}{d\theta})^2+u^2]&=(E_0+GM_Smu)(1+\frac{1}{2}\frac{E_0+GM_Smu}{mc^2})\\\\ =&(E_0+\frac{E_0^2}{2mc^2})+(GM_Sm+\frac{GM_SmE_0}{mc^2})u+\frac{G^2M_S^2m^2}{2mc^2}u^2\\\\ \Longrightarrow\quad\frac{L^2}{m}(\frac{d^2u}{d\theta^2}+u)&=GM_Sm+\frac{GM_SmE_0}{mc^2}+\frac{G^2M_S^2m}{c^2}u\\\\ \end{aligned}$$

与正常天体运动对应不同的部分是等号右边的线性项. 这使得原来解中$\cos\theta$这一项变为新的

$$\cos\sqrt{1-\frac{G^2M_S^2m^2}{L^2c^2}}\theta$$

所以产生了进动. 周期进动角为

$$\begin{aligned} \Delta\theta=(1-\frac{1}{\sqrt{1-\frac{G^2M_S^2m^2}{L^2c^2}}})\cdot2\pi\approx\frac{\pi G^2M_S^2m^2}{L^2c^2}\\\\ \end{aligned}$$

方法2 LRL矢量

LRL矢量指的是Laplace-Runge-Lenz矢量. 这个矢量的定义为

$$\begin{aligned} \bold{B}=\bold{p}\times\bold{L}-GM_Sm^2\hat{r}\\\\ \end{aligned}$$

它在纯平方反比力的情况下守恒，为一个指向近日点的常矢量. 因为它拥有永远指向近日点的特性，我们考虑用它来计算进动角.

先对LRL矢量进行求导，特别需要注意的是引力质量在狭义相对论下守恒，而惯性质量不守恒：

$$\begin{aligned} \frac{d\bold{B}}{dt}&=\frac{d\bold{p}}{dt}\times\bold{L}+\bold{p}\times\frac{d\bold{L}}{dt}-GM_Sm\hat{r}\frac{dm}{dt}-GM_Sm^2\frac{d\hat{r}}{dt}\\\\ &=-\frac{GM_Sm\hat{r}}{r^2}\times\bold{L}+0-GM_Sm\hat{r}\frac{dm}{dt}-GM_Sm^2\dot{\theta}\hat{\theta}\\\\ &=-GM_Sm\hat{r}\frac{dm}{dt}=-\frac{GM_Sm}{c^2}\frac{dE}{dt}\hat{r}\\\\ &=-\frac{GM_Sm}{c^2}(\bold{v}\cdot\bold{F})\hat{r}=-\frac{G^2M_S^2m^2}{r^2c^2}\dot{r}\hat{r}\\\\ \end{aligned}$$

另外一种分解方式：

$$\begin{aligned} \frac{d\bold{B}}{dt}=\bold{\omega}\times\bold{B}+\dot{B}\frac{\bold{B}}{|\bold{B}|}\\\\ \end{aligned}$$

即，$\bold{B}$的变化率来自于旋转和长度变化两个正交的成分. 在这里，我们显然需要的是旋转的部分，所以有投影

$$\begin{aligned} \omega&=\frac{G^2M_S^2m^2\dot{r}}{c^2Br^2}\sin\theta=\frac{G^2M_S^2m^2}{c^2B}\frac{du}{d\theta}\dot{\theta}\sin\theta\\\\ \end{aligned}$$

为了处理$du/d\theta$，先要知道的是，$|\bold{B}|=mvL$，这样我们就可以进行下一步操作：

$$\begin{aligned} L^2\frac{du}{d\theta}&=m^2r^4\dot{\theta}^2\cdot\frac{\dot{r}}{r^2}\cdot\frac{1}{\dot{\theta}}=m\dot{r}L\\\\ &=mv\sin\theta\cdot L=B\sin\theta\\\\ \end{aligned}$$

代入原来的式子，得到

$$\begin{aligned} \omega&=\frac{G^2M_S^2m^2}{c^2B}\cdot\frac{B}{L^2}\sin\theta\cdot\frac{L}{mr^2}\sin\theta=\frac{G^2M_S^2m}{c^2L}\frac{\sin^2\theta}{r^2}\\\\ \end{aligned}$$

周期进动角

$$\begin{aligned} \Delta\theta=\int_0^T\omega \text{d} t=\frac{G^2M_S^2m^2}{c^2L^2}\int_0^{2\pi}\sin^2\theta\text{d} t=\frac{\pi G^2M_S^2m^2}{c^2L^2}\\\\ \end{aligned}$$

与前面一种方法算出来的结果相同.

Part 3 广义相对论效应

（由于我正在学习广义相对论的路上，所以这里直接给出等效结论）

$$\begin{aligned} (\frac{\text{d} r}{\text{d}\theta})^2-\frac{m^2e^2r^4c^2}{L^2}+r^2(1+\frac{m^2r^2c^2}{L^2})(1-\frac{2GM_S}{c^2r})=0\\\\ \end{aligned}$$

令$u=1/r$，并求导：

$$\begin{aligned} (\frac{\text{d} u}{\text{d} \theta})^2+u^2&=\frac{m^2(e^2-1)c^2}{L^2}+\frac{2GM_Sm^2}{L^2}u+\frac{2GM_S}{c^2}u^3\\\\ \Longrightarrow\quad\frac{\mathrm{d}^2 u}{\text{d}\theta^2}+u&=\frac{GM_Sm^2}{L^2}+\frac{3GM_S}{c^2}u^2\\\\ \end{aligned}$$

仍然是一个附加势的问题，与Part 1类似，我们考虑使用摄动法. 零阶解：

$$\begin{aligned} u_0=\frac{GM_Sm^2}{L^2}(1+e\cos\theta)\\\\ \end{aligned}$$

代入方程求一阶解：

$$\begin{aligned} \frac{\mathrm{d}^2u_1}{\text{d}\theta^2}+u_1&=\frac{3GM_S}{c^2}u_0^2=\frac{3G^3M_S^3m^4}{L^4c^2}(1+2e\cos\theta+e^2\cos^2\theta)\\\\ \end{aligned}$$

注意到如下方程的通解：

$$\begin{aligned} \frac{\text{d}^2u}{\text{d}\theta^2}+u=1&\longrightarrow u=1+A\sin\theta+B\cos\theta\\\\ \frac{\text{d}^2u}{\text{d}\theta^2}+u=\cos\theta&\longrightarrow u=\frac{1}{2}\theta\sin\theta\\\\ \frac{\text{d}^2u}{\text{d}\theta^2}+u=\cos\mu\theta&\longrightarrow u=\frac{\cos\mu\theta}{1-\mu^2}\\\\ \end{aligned}$$

将$\cos^2\theta$用二倍角公式化为$\cos2\theta$和常数的组合就可应用上面的公式.

可以得到一阶解为

$$\begin{aligned} u_1=\frac{3G^3M_S^3m^4}{L^4c^2}(1+e\theta\sin\theta+\frac{1}{2}e^2-\frac{1}{6}e^2\cos2\theta)\\\\ \end{aligned}$$

产生进动效应的只有$\theta\sin\theta$项，在研究进动问题时可以只计入这一项与零阶项的影响. 这样，整体的解就变为

$$\begin{aligned} u(\theta)=\frac{GM_Sm^2}{L^2}[1+e(\cos\theta+\frac{3G^2M_S^2m^2}{L^2c^2}\theta\sin\theta)]\\\\ \end{aligned}$$

因为$G^2M_S^2m^2/L^2c^2$是小量，所以近似有

$$\begin{aligned} u(\theta)&\approx\frac{GM_Sm^2}{L^2}[1+e(\cos\theta\cos(\frac{3G^2M_S^2m^2}{L^2c^2}\theta)+\sin(\frac{3G^2M_S^2m^2}{L^2c^2}\theta)\sin\theta]\\\\ &=\frac{GM_Sm^2}{L^2}\{1+e\cos{[(1-\frac{3G^2M_S^2m^2}{L^2c^2})\theta]}\}\\\\ \end{aligned}$$

故周期进动角：

$$\begin{aligned} \Delta\theta=2\pi\cdot\frac{3G^2M_S^2m^2}{L^2c^2}=\frac{6\pi G^2M_S^2m^2}{L^2c^2}\\\\ \end{aligned}$$

此结果恰好为狭义相对论结果的6倍.

总结与反思

写完这个专题历时不长，真正花在完成内容上的时间应该只有不到一天. 但是时间上的跨度跨越了开学典礼与入学的各项手续，真正开启了大学生活之后，我逐渐意识到，时间确实相当紧张，我变得不太有时间做这些回顾性的工作，接下来的时间应该留给新知识的学习和能力的提升.

或许我之后不会再写计算专题，或许我假期有时间还会再写，但是不得不承认这段时间我的计算能力还是因为这个有一小点提升，这是有意义的. 在开学之后的时间，我可能会每周更新一次blog，算是对自己的一种督促，也让自己每周有一个总结和回顾，有一点小小的成就感.

先暂时到这里为止吧.