计算专题4 水星进动

本文最后更新于 2024年8月16日 晚上

计算专题4 水星进动

前言

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

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

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

Part 1经典附加力场(大行星影响)

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

我们标记ABC=θ\angle ABC=\thetaADC=α\angle ADC=\alpha,考虑圆心角dθd\theta的一段弧对DD点产生的引力场:

dg=MJ2πdθG1RJ2+r22RJrcosθ\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}

整个环产生的引力场:

g=r^02πdgcosα=r^02πGMJ2πRJcosθr(RJ2+r22RJrcosθ)3/2dθ\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α\cos\alpha的计算,应用正弦定理和余弦定理:

cosα=1sin2αsinα=sinADB=ABADsinθ=RJsinθRJ2+r22RJrcosθcosα=RJ2cos2θ+r22RJrcosθRJ2+r22RJrcosθ=RJcosθrRJ2+r22RJrcosθ\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/RJr/R_J为一阶小量,展开到一阶项得到:

g=r^GMJRJ2πRJ302πcosθrRJ(12rcosθRJ+r2R2)3/2dθr^GMJRJ2πRJ302πcosθ(1rRJcosθ+3rcosθRJ)dθ=r^GMJ2πRJ202π(3cos2θ1)rRJdθ+0=GMJr2RJ3r^\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面即可.

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

r¨rθ˙2=GMSr2+GMJr2RJ3,r2θ˙=Lm\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}

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

a¨+δ¨GMSa2(12δa)+GMJa2RJ3(1+δa)+L2m2a3(13δa)=GMSa2+GMJa2RJ3+L2m2a3+(2GMSa3+GMJ2RJ33L2m2a4)δ\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}

注意到在零阶情况下,

a¨=GMSa2+GMJa2RJ3+L2m2a3=0L2m2a3=GMSa2GMJa2RJ3\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}

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

代入上述一阶方程可知,

δ¨=(2GMSa3+GMJ2RJ33GMSa3+3GMJ2RJ3)δ=(1+2MJa3MSRJ3)GMSa3δ\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}

故,径向振动角频率

ωr=GMSa3(12MJa3MSRJ3)GMSa3(1MJa3MSRJ3)\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}

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

ω0=Lma2=GMSa3(1MJa32MSRJ3)GMSa3(1MJa34MSRJ3)\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}

所以进动角频率为

Ω=ωrω0=3MJa34MSRJ3GMSa3\begin{aligned} \Omega=|\omega_r-\omega_0|=\frac{3M_Ja^3}{4M_SR_J^3}\sqrt{\frac{GM_S}{a^3}}\\\\ \end{aligned}

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

周期进动角为

Δθ=ΩT=3πMJa32MSRJ3\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方程

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

Ek=p2c2+m2c4mc2=mc2v2c2(1v2c2)+1mc2=mc211v2c2mc212mv2+38mc2(v2c2)2\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}

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

E0=GMSmr+12mv2+38mc2(v2c2)2v2c22E0+GMSmrmc2\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)O(1)v2/c2v^2/c^2的结果,这是为了在之后处理(v2/c2)2(v^2/c^2)^2时更为方便. 注意到角动量守恒,则

L=γmr2θ˙=const.θ˙=Lγmr2,r˙=θ˙drdθ=Lγmr2drdθ\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/ru=1/r,有

θ˙=Lu2γm,r˙=Lγmdudθ\begin{aligned} \dot{\theta}=\frac{Lu^2}{\gamma m}\,,\quad\dot{r}=-\frac{L}{\gamma m}\frac{du}{d\theta}\\\\ \end{aligned}

能量守恒式写为

12mv2+38mc2(v2c2)2GMSmr=E012m(r˙2+r2θ˙2)+38mc2(2E0+GMSmrmc2)2=E0+GMSmrL22m[(dudθ)2+u2](1v2c2)=(E0+GMSmr)(132E0+GMSmrmc2)L22m[(dudθ)2+u2]=(E0+GMSmr)(1+12E0+GMSmrmc2)\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}

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

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

L22m[(dudθ)2+u2]=(E0+GMSmu)(1+12E0+GMSmumc2)=(E0+E022mc2)+(GMSm+GMSmE0mc2)u+G2MS2m22mc2u2L2m(d2udθ2+u)=GMSm+GMSmE0mc2+G2MS2mc2u\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θ\cos\theta这一项变为新的

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

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

Δθ=(111G2MS2m2L2c2)2ππG2MS2m2L2c2\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矢量. 这个矢量的定义为

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

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

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

dBdt=dpdt×L+p×dLdtGMSmr^dmdtGMSm2dr^dt=GMSmr^r2×L+0GMSmr^dmdtGMSm2θ˙θ^=GMSmr^dmdt=GMSmc2dEdtr^=GMSmc2(vF)r^=G2MS2m2r2c2r˙r^\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}

另外一种分解方式:

dBdt=ω×B+B˙BB\begin{aligned} \frac{d\bold{B}}{dt}=\bold{\omega}\times\bold{B}+\dot{B}\frac{\bold{B}}{|\bold{B}|}\\\\ \end{aligned}

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

ω=G2MS2m2r˙c2Br2sinθ=G2MS2m2c2Bdudθθ˙sinθ\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θdu/d\theta,先要知道的是,B=mvL|\bold{B}|=mvL,这样我们就可以进行下一步操作:

L2dudθ=m2r4θ˙2r˙r21θ˙=mr˙L=mvsinθL=Bsinθ\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}

代入原来的式子,得到

ω=G2MS2m2c2BBL2sinθLmr2sinθ=G2MS2mc2Lsin2θr2\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}

周期进动角

Δθ=0Tωdt=G2MS2m2c2L202πsin2θdt=πG2MS2m2c2L2\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 广义相对论效应

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

(drdθ)2m2e2r4c2L2+r2(1+m2r2c2L2)(12GMSc2r)=0\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/ru=1/r,并求导:

(dudθ)2+u2=m2(e21)c2L2+2GMSm2L2u+2GMSc2u3d2udθ2+u=GMSm2L2+3GMSc2u2\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类似,我们考虑使用摄动法. 零阶解:

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

代入方程求一阶解:

d2u1dθ2+u1=3GMSc2u02=3G3MS3m4L4c2(1+2ecosθ+e2cos2θ)\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}

注意到如下方程的通解:

d2udθ2+u=1u=1+Asinθ+Bcosθd2udθ2+u=cosθu=12θsinθd2udθ2+u=cosμθu=cosμθ1μ2\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}

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

可以得到一阶解为

u1=3G3MS3m4L4c2(1+eθsinθ+12e216e2cos2θ)\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}

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

u(θ)=GMSm2L2[1+e(cosθ+3G2MS2m2L2c2θsinθ)]\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}

因为G2MS2m2/L2c2G^2M_S^2m^2/L^2c^2是小量,所以近似有

u(θ)GMSm2L2[1+e(cosθcos(3G2MS2m2L2c2θ)+sin(3G2MS2m2L2c2θ)sinθ]=GMSm2L2{1+ecos[(13G2MS2m2L2c2)θ]}\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}

故周期进动角:

Δθ=2π3G2MS2m2L2c2=6πG2MS2m2L2c2\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,算是对自己的一种督促,也让自己每周有一个总结和回顾,有一点小小的成就感.

先暂时到这里为止吧.


计算专题4 水星进动
https://physnya.top/2024/08/16/计算专题4 水星进动/
作者
菲兹克斯喵
发布于
2024年8月16日
更新于
2024年8月16日
许可协议