洛希极限完整推导

大家好，最近很多小伙伴都看了电影流浪地球，里面涉及到了叫作洛希极限的概念。本专栏会对洛希极限进行完整推导，需要一定量的物理知识基础。另外，我在观看时发现流浪地球里的洛希极限数据好像有点问题。

洛希极限粗略计算的数据如下：

木星对地球的刚体洛希极限是在实际上是在木星内部的，刚体洛希极限大概是0.8倍木星半径(约5.6万公里)，流体洛希极限约为1.5倍木星半径(约10.5万公里)

但这和电影里的数据相差很多。电影里刚体洛希极限是890745427米，流体是1713024931米。

后来我发现原来制作组用错数据了，把太阳相对木星的洛希极限当成了木星相对地球的洛希极限弄上去了。制作组审核时应该更认真一下。不过，这也可能是为了有这样的剧情而妥协的。

洛希极限推导：

接下来我们来了解洛希极限推导，之前我在评论区简易推导了刚体的洛希极限。不过那个其实并不完整，没有考虑地球加速度导致的自转惯性离心力。在这个专栏会进行相对完整的推导，并且还将加上流体洛希极限的推导，希望大家喜欢。

卫星在不分解的情况下接近主星的极限距离取决于卫星的刚性。最极端的情况下，一个完全刚性的卫星将一直保持其形状，直到潮汐力将其分开。而相反，流动性高的卫星由于潮汐力的增加会使其变形，导致卫星被拉伸长，使得其更容易分裂。

大多数真实的卫星的洛希极限都位于这两个极端之间，一般的卫星既不是完美的刚体也不是完美的流体。

刚体洛希极限：

刚体不会考虑在到达洛希极限前潮汐力导致的变形。因此，刚体洛希极限的假设虽然不太切实际，但却可大大简化了计算。下面是不考虑卫星加速度导致的自转离心作用的。

为了确定洛希极限，需要引入一个小质量u，它位于卫星距离主星较近的一侧。

u这部分会受到两种力：卫星对u的引力和主星对u的潮汐力。卫星对u的引力可通过万有引力定律表示。

潮汐力是靠近主要卫星边缘的引力的差异导致的：

由于一般情况下：r<<R，R<d，分母可以近似调整一下，得到下面的式子。

当引力和潮汐力相互平衡时，达到了洛希极限。

因此：

由于质量可以改成密度与体积的乘积。上面的式子还可以进行调整。

消去共同因数后可以得到简易的刚体洛希极限表达式。

更精确的公式：

由于近距离卫星很可能在同步自转的近圆轨道上轨道运行，因此还需考虑旋转产生的离心力对其的影响。其表达方式如下：

将离心力加入上面潮汐力的等式中：

但由于主星的密度可能未知，我们也可以换一种表达方式。

因此，得到主星的质量并估计卫星的密度就足以计算主星相对卫星的洛希极限。

那么如果更严谨一些，考虑自转产生的动态坐标影响呢？

先大致模拟出其状态和参数。

然后表达出中心径向力的平衡。

由于刚体会维持稳定，因此得到下面的平衡等式。

将应力分量进行二阶线性展开：

将上述式表达成克罗内克函数的形式，并带入刚体的平衡等式。

作为刚体，将变形系数α=1，带入，得到：

解出几种可解坐标情况：

将y轴和z轴作为解析平面解析：

根据不同的刚体裂解程度可以分析出以下一种不同的洛希极限表达式：

将这些结果与最上面的基本结果作比较：

下图是上述几种不同结果的分裂曲线：

换个角度，我们来分析以下被环绕的天体，同样列出坐标等式：

下面是分别处于刚体和流体下两种不同环绕天体的解离程度变化：

由于β<1，我们可以将上式做近似的变换。

同样地，解析球平面坐标：

得到几种结果表达式：

消去相应系数，得到最简结果：

这样的环绕天体又产生了新的分裂曲线：

至此，刚体洛希极限基本推导完毕，已经尽量考虑到了所有能考虑的情况。计算是可以根据不同情况选用不同的公式。

流体洛希极限：

流体洛希极限考虑了卫星的变形。一个极端的情况是一颗被潮汐锁定的流体卫星绕着主星运行，任何作用在卫星上的力都会使它变形成为长椭球体。

其计算十分复杂，洛希本人提供了以下的结果：

本专栏提供的推导方法与洛希本人的不同，up主查询了洛希本人的文献，结果都是法语看不懂.....，且洛希本人使用的方法是洛希坐标法，相对比较复杂，我这里提供的办法是振动稳定性法和最小势能原理法，相对更简洁一点，但仍需要一定物理知识。

设有一质量为m，密度均匀分布为ρ，半径为a1的球体，此物体以角速度ω自转。另有一引起潮汐力的质量为M的球体，它们间的距离为R。

在自转和公转同步性假设下，有：

假设其为均匀或不均匀的不可压流体，在自转与潮汐力的作用下，它的变形后的平衡形状可表达为：

其中：

在均匀不可压情况下，可得到：

对均匀不可压介质，以其变形形态为基本态，考虑附加扰动。则线性化扰动方程为：

其中ρ0为常数，代表均匀介质密度

P为扰动压力

U、V、W分别为α、θ、φ向的扰动速度

Ψ为自引力扰动势

因为是不可压缩介质，故附有不可压缩条件：

扰动引力势满足以下方程：

对所有扰动量作展开：

精确到一阶小量，用式(4)把式(5)和(6)简化成：

把式(4)乘以sin²θ，式(10)乘以sinθ，然后把式(9)代入式(4)，(10)和(11)，即得：

对不可压缩条件(7)，在把式(9)代入后可得：

其中：

扰动引力势方程(8)在用了式(9)后为：

其中：

方程(12)，(13)，(14)，(15)和(17)是我们要解的基本方程。变量为Unm，Vnm，Wnm，Pnm，Ψnm，它们都仅为α的函数。

在这些方程中，带不同上标与下标的E，A，C，K，L和N都是常系数。

方程(12)，(13)，(14)，(15)和(17)成立的条件即要求在这些方程中带有同样项：

并且它们的系数和为0。由此可以导出Unm，Vnm，Wnm，Pnm，Ψnm，所满足的方程。显然这是一个具有无穷个未知量的无穷维线性微分方程。

首先来考察其某些性质：

1，在这一微分方程组中，有不同n和m的未知函数Unm，Vnm，Wnm，Pnm和Ψnm是互相耦合的。且容易看出，具有m为偶数的未知量相互耦合，具有m为奇数的未知量也相互耦合。但具有m为奇数的未知量是独立于m为偶数的未知量的。以后仅讨论m为偶数的情形，对m为奇数情形是不难同样讨论的。

2．科里奥利力的出现使方程(5)和(6)不再恒等。Wnm≠0的引入解决了这一问题。并且，如果变形的平衡状态（即基本态）是不对称的话，容易看出，在周向两方向传播的波也是耦合的。

特别在m=0情形，扰动量必须取如下形式：

这样在方程(12)-(18)中，当n=0时，为应用表达式(19)，应该把m和mwπ0分别用0和wπ0来替代。

这是不同于旋转星的振动稳定性情形。但相似于既考虑自转又考虑潮汐的振动稳定性情形。

3．因为f02和f2在均匀不可压情形为常数，所以整个微分方程是欧拉型的，这些微分方程的解可表为α^α，其中α为特征值。

边界条件：

由上面可以知道，应有的边界条件为：

在中心处，a=0

1.压力P约有限，

2.扰动引力势的梯度有限；

在表面处a=a1，压力为0，即Pnm=0

由于附加扰动，原来变形后的平衡态的表面产生了附加变形。因此必须考虑表面物质的位移。这一受扰动的物质迁移可以作为面质量看待．而面质量产生的扰动引力势梯度应有关系：

其中n为平衡态表面的单位外法线值，v为径向位移。从径向速度表达式，有：

从方程(23)可知，我们既应求得结构内部的扰动引力势，也应求得结构外部的扰动引力势。显然外部扰动引力势与内部扰动引力势在表面处要相等，即：

方程(23)可表达为：

结果：

从上面的讨论可知，最后归结为解在一定边界条件下的无穷个未知函数的相互耦合的无穷个微分方程。为实际计算，可以限于讨论到n=2的情形，这时共有十九个未知量：

1．从四个方程(12)-(15)可得以u，v，w和p为未知函数的常微分方程组。其非线性项系以毋来表示的。

2．分别利用球心处有限与无穷远处趋于零的条件，由方程(17)求得内、外扰动引力势的表达式。并利用条件(25)，把内、外扰动引力势的某些待定常数先确定

3．把解得的内部扰动引力势的表达式代入由步骤1所得的微分方程的右边。这样可以利用边界条件(20)和(22)求得u，v，w和p的表达式

4．把所求得的u，Ψ+和Ψ-代入边界条件(23)，就可得到决定四个待定常数的齐次线性方程组。为了得到不全为0的解，其充要条件是其系数行列式为0，这样就得到决定特征频率σ的方程，这一特征频率可以是实数，也可以是复数，主要决定于g02和g2的取值。当频率为实数时，这一原先的平衡状态为稳定的；而当频率为复数时，主要决定于其虚部。当虚部由正变负时，原先的平衡状态也就变成不稳定的了。所以频率的虚部为0的状态为临界状态。洛希极限就由这临界状态求得。

像已经知道的那样，在上面求解常微分方程时都需先解特征方程。而恰恰这些特征方程的系数包含有洛希极限项（即g02和g2）。因此要用直接方法来解这特征方程是非常困难的。

我们提出一间接方法：首先，任选M/(πR³ρ)的值。这样就可向上面所描述的过程进行求解，因为M/(πR³ρ)的任意性，一般来说特征频率σ为复数。从而改变M/(πR³ρ)的值，就可以得到不同的特征频率。这样的过程一直持续到特征频率的虚部接近于0。在这一间接方法中，重要的是初值的选取。整个过程可由计算机来实现。

最终，我们在m<<M的情形下，用本文方法得到的洛希极限为M/(πR³ρ)=0.0875，偏心率e=0.875。虽然不同的计算方法存在着一定差异，但我们可以在方程(12)-(17)中取更多的项来改进其精度。更重要的是，这一方法可以推广到非均匀介质的情形。将此数据推算便可得到洛希的结果。

那么如果考虑物体成为洛希椭球体，并进行不均匀的变换结果又会如何呢？

流体更复杂的情况：

主星的引潮力关于二星连线具有旋转对称性，故它被拉伸为长旋转椭球。这种一般发生在小天体上，如坠入木星的苏梅克-列维9号彗星。接下来，我会使用最小势能原理法，利用椭球模型来分析这种相对理想的情况。

这里先通过引力场与静电场的类比来引入一个概念: 一个质量为 m 的物体的引力自能。静电场的基本规律是库仑定律:所带电荷量分别为 q1、q2的二个点电荷之间的作用力。引力场的基本规律是牛顿 引力定律: 质量分别为 m1、m2的二个质点之间的作用力。由于二者遵守的基本规律相似，故有很多相同概念，如势能、势、场强等。对比表达式可知，若采用字符替换规则:，则可以得到引力场情况下的对应的公式而不必重新推导。一个带电体的静电自能定义为把它的各部分电荷从无限远分散的状态聚集成该带电体时外界抵抗电荷之间静电力所做的功，用 We 表示。

其中 ρe(r)是电荷密度。由于同性电荷之间总是斥力作用，故 We＞0。同理，一个物体的引力自能定义为把它的各部分质元从无限远分散的状态聚集成该物体时外界抵抗质元之间引力所做的功，用Wg表示。采用上段符号替换规则，则有：

其中ρ(r)是密度。由于质元之间总是引力作用，故Wg＜0。根据均匀带电的长旋转椭球的静电自能表达式按上边符号替换规则，可得密度均匀的长旋转椭球状物体的引力自能:

设伴星原来半径为R，质量为m，密度为ρ'且不可压缩，它经过主星附近时变成一个半长轴为a、半短轴为b的长旋转椭球。由于体积保持不变，根据式(30)，可得伴星的引力自能为：

以Gm²/R为能量单位，做Wg随e变化的曲线。结果表明，e越大，自能越大，当e=0即球形时自能最低。故由最小势能原理可知，当伴星是孤立物体时，它应该呈现球形。

密度均匀的长旋转椭球在引潮力场中的势能如下图，设主星质量为m1，密度为ρ，伴星质心O与主星质心O1距离为r。以O为原点建立坐标系，其中z轴沿O1O的方向。在此坐标系中，伴星椭球所占的空间区域V表示为：

在V内点(x，y，z)处取体元dv，由式(28)可知，质元dm=ρ'dv在主星引潮力场中的势能:

对上式积分，可得伴星在引潮力场中的势能:

做坐标变换:x=γbsinθcosφ，y=γbsinθsinφ，z=γacosθ，体积元变换dv=ab²γ²sinθdγdθdφ，则式(32)表示的区域V变换为:0≤γ≤1，0≤θ≤π，0≤φ≤2π。考虑到质量与体积密度的关系m=Vρ'，则式(34)中关于x²的积分等于:

同理：

将上3式代入式(34)得：

从上式容易看出，e越大，Wt越低，这与Wg的变化趋势相反，说明引潮力有把球形天体拉长的趋势。

平衡条件根据式(31)和式(35)，伴星系统总能量:

由上式可知，当m1、m、R、r确定后，系统总能量是e的函数:W=W(e)。根据最小势能原理，系统的稳定平衡条件是:系统在某偏心率e0处能量有极小值。只有这样，系统在受到微扰偏离平衡时，微扰消失后才会自动恢复到平衡状态。例如，以Gm²/R为能量单位，分别取λ=0.06、0.07、0.08，做W随e变化的曲线，三条曲线自上而下对应的λ是递增的。显然点A是极小值点，该点对应的e值为平衡态偏心率e0。而点B是极大值点，系统在此点很不稳定。根据数学分析知识，系统的稳定平衡，即能量取极小值的条件为：

以Gm²/R为能量单位，对式(36)进行无量纲化，这时能量用小写w表示:

其中，λ=(m1/m)(R/r)³。利用计算机作求导运算:

我们先通过上边二式，研究当主星质量m1、伴星质量m、伴星半径R确定后，伴星从远处逐渐靠近主星，即r逐渐减小的过程中，它的平衡形状偏心率的变化问题。不妨设m1等于地球质量，m等于月球质量，R等于月球半径，则λ=81.3(R/r)³。计算结果见下表。表中第1栏列出不同的r/R值，第2栏为对应的λ值的100倍，第3栏是与上边r和λ对应的、通过式(39)和方程W'(e)=0求出的数值解e0，第4栏是与上述r和λ、e0对应的、通过式(40)确定的、关于W″(e0)是否大于0的判断。数据表明，对于给定的所有距离比r/R或λ，W″(e0)均大于0，均存在平衡形状。下图绘制了表中给出的其中4个不同的平衡形状，曲线旁边的数字0.8、0.6、0.4、0.2表示对应的椭球偏心率e0。由图表可以看出，随着伴星与主星距离r的减小，即λ的增大，伴星的形状偏心率越来越大，它被拉得越来越长。不难想象，当r减小到某临界值rc时，λ增大到某临界值λc，平衡条件刚好不满足，系统即将失去平衡发生破裂。rc就是过客天体的洛希极限。设此时的形状偏心率为ec，则令式(37)第二式右边由“＞0”变为“=0”，可得临界条件：

将式(39)、式(40)代入方程组(41)，用计算机对未知量λc、ec进行数值求解，得：

根据上边的初步分析知，λ＜λc时存在平衡状态。

下面通过一阶导数W'(e)随e变化的曲线，分析系统在λ=λc附近的平衡问题。分别取λ＜λc、λ=λc、λ＞λc时，在同一图中作3条W'(e)－e曲线，重点考察每条曲线与横轴的交点，也就是考察函数W(e)取极值必要条件W'(e)=0对应的点是不是对应平衡状态。

首先可以看到，3条曲线都在原点与横轴相交，在此点对应的e=0，即对应球形。显然，从物理上判断，系统稳定在球形是不可能的，因为引潮力会把流体天体拉长;从数学上讲，图上表明3条曲线在该点的斜率皆小于0，即W″(e0)＜0，函数W(e)在此点有极大值，故无论λ取何值，原点都不是平衡状态点。

除原点之外，λ＜λc对应曲线与横轴有二个交点A、B，由图可见，曲线在点A的斜率大于0，即W″(e0)＞0，函数W(e)在此点有极小值，故点A是平衡状态。在点B的斜率小于0，即W″(e0)＜0，函数W(e)在此点有极大值，故点B不是平衡状态。所以有结论:λ＜λc时存在平衡状态;除原点之外，λ=λc对应曲线与横轴仅有一个交点———切点C;除原点之外，λ＞λc对应曲线与横轴无交点，故此时不满足平衡条件。所以λ=λc确实是系统有无平衡状态的分界点。

至于λ=λc时是否有平衡状态，由于对应曲线在点C与横轴相切，即W'(ec)=0，W″(ec)=0，需要看W(e)在点C的三阶导数是否等于0和四阶导数是否大于0。经数值计算，发现三阶导数不等于0，故此点不存在极值，更谈不上有极小值，所以此时不存在平衡状态。综上所述，λ≥λc时系统是不稳定的。由此不等式和如下关系:λ=(m1/m)(R/r)³、m1=ρ4πR³/3、m=ρ'4πR³/3、式(42)，可得洛希极限：

如果引入扁率的变量c/R，重新解方程(41)，那么ec会稍微变小一点，并得到下面的相对更加精确的结果：

至此，流体洛希极限推导完毕，计算时公式(27)、(43)、(44)都可以近似使用。