哈希极限与天体的组成物质有关麼?或者通俗点说与星球的“硬度”有关么?

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

解出几种可解坐标情況：

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

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

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

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

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

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

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

同样地解析球平面坐标：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在这些方程中，带不同上标与下标的EA，CK，L和N都是常系数

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

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

首先来考察其某些性质：

1在这一微分方程组中，有不同n和m的未知函数UnmVnm，WnmPnm和Ψ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在均匀不可压情形为常数，所以整个微分方程是欧拉型的这些微分方程的解可表为α^α，其中α为特征值。

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

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

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

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

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

方程(23)可表达为：

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

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

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

3．把解得的内部扰动引力势的表达式代入由步骤1所得的微分方程的右边。这样可以利用边界条件(20)和(22)求得uv，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内点(xy，z)处取体元dv由式(28)可知，质元dm=ρ'dv在主星引潮力场中的势能:

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

将上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)=0W″(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)都可以近似使用