扁率对行星自由核章动与低阶本征模耦合影响的初步研究∗

张冕黄乘利

(1中国科学院上海天文台上海200030)

(2中国科学院行星科学重点实验室上海200030)

(3中国科学院大学北京100049)

1 引言

行星不是严格意义上的刚体,它们的自转会使其形变,大多数行星的形状可以近似看作2阶椭球.1735年左右,巴黎科学院派遣了两支测量队,一支去了秘鲁,另一支去了拉普兰,试图用实验的方法测定地球椭球的扁率.1743年,克莱劳发表了他的关于地球扁率是深度函数的著名微分公式,这一公式基于流体静平衡的理论,所以至今还是基本成立的[1].除了克莱劳和拉道的经典方法外,Lichtenstein完善了Liapounoff的积分方程法,Wavre提出了一种平衡曲面的方法,Macke发展了一套现代的位能守恒方法[1].此外,Chandrasekhar[2]用惯性矩的方法详细讨论和研究了一般星体的平衡形状.

形变后的行星的平衡形状通常可以看作1阶椭球体,这会引起行星的某些物理性质发生变化.这些变化引起的现象则可以反过来检视行星的扁率.本文通过自由核章动的变化和低阶本征模的耦合来讨论这个问题.自由核章动(Free Core Nutation,FCN)是由于地核与地幔的旋转轴不一致引起的一个近周日的旋转模,它连接了天文观测与地球内部.核幔边界处的扁率对此模周期有很大的影响,根据Neuberg等[3]和黄乘利[4]的计算,发现只需把核幔边界处的扁率增加5%,就可以使原本理论计算得到的460 sd的周期,变成与观测一致的430 sd.本文则通过整体扁率的放缩,研究自由核章动周期的变化.

土星在太阳系内有着最大的扁率,它的高速旋转和大扁率都会使低阶本征模产生巨大的耦合.Vorontsov等[5]最早研究了土星在自转和扁率影响下的自由振荡,他们采用了扰动法,即把自转和扁率看作小扰动,对原始的无自转和无扁率下的结果进行修正.但是土星自转很快,其低阶的自由振荡周期与自转周期相比,已经不适合看做小扰动,且土星的大扁率使低阶的本征模强烈耦合,在此情况下使用扰动法是值得商榷的.本文采用直接计算法对这个问题进行了初步的讨论.

2 扁率对自由核章动的影响

Hough[6]于1895年利用由液态核和固体壳组成的两层均质地球模型预言了一个近周日的摆动,此即自由核章动.既然固体核是均质的,那么它在核幔边界处的动力学扁率ec与其形状扁率是一致的.利用这个模型可以得到一个退行的近周日摆动,它存在一个频率:

其中,Ω、A和Ac分别是自转角速度、整体的赤道惯性矩和内核的赤道惯性矩,ν则是自由核章动角频率.对于初步地球参考模型(preliminary reference Earth model,PREM)[7],(1)式会得到350 sd的自由核章动周期.但是Hough预测的自由核章动在当时并不被认可,直到1926年,Jeffreys[8]利用地震学资料证实了液体核的存在.Hough[6]采用的是角动量守恒法来计算自由核章动周期,这是一种经典而古老的方法,对于复杂的地球模型,现代常用的是由Smith[9−10]提出并完善的线形动量法.它把复杂动力学方程转换成常微分方程组,把需要计算的变量进行一定精度的截断,并把这些变量从扁球体计算域转换到等效的球体计算域.这些变量在地心处根据正则边界条件构造出初值,利用常微分方程组从地心积分到地表,之后在自由表面检验是否满足边界条件来判断是否为正确解.

由(1)式可以看出:自由核章动与核幔边界处的扁率有着重要的联系,Neuberg等[3]和黄乘利[4]的计算发现:只需把核幔边界处的扁率增加5%,就可以使原本用Smith方法计算自由核章动得到的460 sd,变成与观测一致的430 sd.需要强调的是:不能通过人为修改核幔边界处的扁率令理论结算结果与观测结果一致.与以往研究不同,本文通过整体放缩扁率剖面来计算扁率与自由核章动的关系.

2.1 基本方程和边界条件

在弹性各向同性的流体静平衡的固体圈层中,小位移引起的周期振荡的动力学方程可以写成[9,11−12]

其中,λ和µ为Lame参数,为单位并矢.在等熵的无粘流体圈层中,小振动的动力学方程可以写为[14]:

其中,ρ1和p1分别定义为:

其中,α为压缩波波速.液体圈层中应力张量可以写成:

其中,G为引力常数.在地心处要求各物理量必须正则.在自由表面处,需满足如下条件:

2.2 地球模型和扁率

本文采用的是初步地球参考模型[7],把地球分为12层,去掉了对自由核章动不敏感的海洋层.物理参数表示采用的是12层模型,而未知变量表示采用了3层模型:固体内核、液体外核和固体的地幔(加地壳).之后对PREM模型进行扁率改正,从Wavre的积分微分方程出发[1]:

等密度面(等势面)r=r0[1+ϵ2P2(θ)],其中,ϵ2为系数,r0为等效球面的半径,P2为2阶勒让德函数,θ为极角.解出方程(14)中的ϵ2,在1阶近似下可以得到扁率. 方程(14)可以使用Galerkin谱方法解得,使用方程(14)的优点是:不需要像Clairaut方程一样给出边界条件,只需知道行星的平均等效半径和密度剖面.这是因为Clairaut方程是方程(14)进行2次求导后的结果,丢失了信息.本文计算得到的扁率剖面与黄乘利[4]的结果是一致的,例如在核幔边界处ϵ2值都为0.00254656.

2.3 变量的表示

球坐标系中球谐矢量变量可以表示为[12]:

2.4 计算方法

本文采用了Galerkin方法来计算自由核章动,而非Smith[10]的传统方法.Seyed-Mahmoud等[15]曾用Galerkin方法解算过自转可压缩自引力流体的振荡,得到了不错的结果.本文把Galerkin方法应用到3层模型,来计算自由核章动.令未知量写成幂级数的组合:

其中,ck是系数.把试探函数χj乘到待解方程的两边,积分可以得到

其中,R是边界.上式是对球体的积分,如果求解域为扁球体,则可写成:

同时球域上的物理参数X(R)修正为扁球体上的. 加入边界条件后,解算方程(23)就能得到本征模的频率.

2.5 扁率对自由核章动的影响

具体做法为:在得到扁率剖面f后,本文不在核幔边界处修改扁率,因为这样会使物理量变化太过突兀,而是采用了全局乘以一个系数进行研究,即:

c与自由核章动周期的计算结果列于表1,可以看出扁率对地球的自由核章动影响很大.这与由Hough-Poincare模型得出的结论是一致的:Hough-Poincare模型下的自由核章动周期T可以由(1)式推导出

可以看出自由核章动周期与核幔边界处扁率的倒数成正相关.

表1 自由核章动周期与扁率放缩系数的关系Table 1 The relation of FCN period and the coeffecient c

3 大扁率对土星低阶本征模的影响

Vorontsov等人最早使用现代方法研究太阳系内类木行星的自由振荡[16].土星在太阳系内有着最大的扁率,它的高速自转和大扁率都会使低阶本征模产生巨大的耦合.Vorontsov等[5]最早研究了土星在自转和具有大扁率下的自由振荡.他们采用了扰动法,即把自转和大扁率看作小扰动,对原始的无自转和无扁率下的模型使用传统的Alterman等[17]的方法进行计算,再于结果上进行小扰动的修正.但是土星自转是高速的,其低阶的自由振荡周期与自转周期相比,已经不适合看做小扰动,且土星的大扁率使低阶的本征模强烈耦合,在此情况下使用扰动法是值得商榷的.在Vorontsov和Zharkov之后,无人讨论高速自转和大扁率下的土星低阶自由振荡.本文采用直接计算法对这个问题进行了初步的讨论.

3.1 土星模型与控制方程

本文采用了多方模型来构建土星模型,它是一种压强和密度的关系[18]:

其中,P是压强,ρ是密度,K是一个常数,常数l是多方指数,表征星体的形态:中子星可以用l在0.5到1之间的多方模型描述;恒星核、褐矮星、气态巨行星甚至类地行星都可以用l=1.5的多方模型描述;主序恒星可以用l=3的多方模型描述[18].在行星科学研究中,气态巨行星也可以用l=1的多方模型[19]描述.在具体解算有自引力、球对称和多方流体的星体内部的状态时,人们常采用一个无量纲的泊松方程:

其中,ξ是无量纲的半径,定义为r=αξ,α为常数;γ是与密度相关的量:ρ(ξ)=ρcγl(ξ),ρc是球心处的密度.(27)式的初始条件为γ(0)=1和γ′(0)=0(γ′是γ的导数). 本文取多方指数l=1.5,平均半径Rp=58242 km,总质量Mp=5.68×1026kg,自转周期Ω0=10.55 h,计算得到土星模型,其半径和密度的关系如图1[20−21].

图1 多方指数为1.5的土星密度剖面Fig.1 The Saturn’s density profile of the polytropic model with n=1.5

这样得到的是球对称的土星模型,对于2阶扁球体的土星模型,则利用(14)式进行2阶扁球体的改正.经过计算可以得到扁率剖面,在土星表面,本文多方模型的扁率为0.0918,这与实际的土星的扁率是很接近的.对于土星和木星低阶的自由振荡的计算,目前的主流做法还是采用弹性力学的手段,例如Vorontsov等[5]、Marley[22]、Gudkova等[23]和Fuller等[24]的计算方法,因为土星和木星不是完美的气态行星,它们有固体核和液体圈层.作为初步研究,本文采用了1层土星模型,采用了2.1节中等熵的无粘流体圈层中小振动的动力学方程:

方程(6)中的压缩波波速可以写为:

其中,β表示与完全绝热分层情况下密度梯度的分离程度,它正比于局部Brunt-Vaisala频率的平方:

本文采用了β=0的绝热模型,以此解出α2.

3.2 计算结果与讨论

本文计算了无自转球体、自转球体和自转椭球体下的低阶本征模,结果见表2.本文先从球体非自转的模型开始,得到了的基频为145.2 min,相较于土星的自转速度10.55 h已经不能看做小量,扰动法不是很适合于这种情况.同时得到的基频周期是117.7 min,的基频周期是103.2 min.

在有自转的扁球体模型中,扁率会使得截断进一步增长.对于模,如果只截断到,会得到基频周期为108.5 min;如果截断到,会得到基频周期为107.9 min;如果截断到,会得到基频周期为112.3 min;如果截断到,会得到基频周期为112.8 min.对于模,如果只截断到,会得到基频周期为109.3 min; 如果截断到,会得到基频周期为114.9 min; 如果截断到, 会得到基频周期为111.3 min; 如果截断到, 会得到基频周期为113.0 min.可见在大扁率下,截断对低阶本征模的影响是巨大的.

此外,Goldreich等[25]认为类木行星的自由振荡本征模的能量来自湍流对流,不同本征模的能量有两种比例关系[26].第1种是等分法,即每个模的能量都是一样的,这是强耦合的情况.第2种是每个模的能量正比于f−13/2(f是模的频率),这是弱耦合的情况.如果ft＞1(t是对流流失的特征时间),则为弱耦合,反之为强耦合.现有的结果似乎暗示土星是弱耦合[22].本文也对这种能量比例对低阶本征模的影响做了初步的计算,假设在有自转的扁球体模型中, 对于模取截断形式, 并写成:

其中,c为系数, 用来表征的振幅比. 需要说明的是: 这只是一个不严格的计算尝试.不同模的能量比可以看做振幅比的平方.c=1时,得到基频的周期为111 min;c=0.1时,可以得到基频的周期为118 min;c=100时,可以得到基频的周期为108 min.而对于有自转的球体,无论取怎样的截断,计算的结果都是同样的125.7 min,可见:在大扁率的耦合作用下,低阶本征模计算时,需要考虑不同模的能量比例,而现在的研究都没有包含这个约束条件.

表2 土星的低阶本征模Table 2The lower degree normal modes of Saturn

对于土星这样自转十分迅速、扁率很大的行星,本文采用1阶椭球近似,求得的内部形状和密度结构与真实条件会有不小的差距.对于星体的平衡形状,著名的钱德拉塞卡用惯性矩的方法做出了很多定性的讨论[2].但若要实际解算扁率剖面来进行本征模的计算,还是需要采用Clairaut(1阶)、Darwin(2阶)、Dennis(3阶甚至更高)的方法,这些方法都是令椭球内部点的重力位中某些系数为零求得的.除了Darwin的经典计算2阶扁率的方法外,Lichtenstein完善了Liapounoff的积分方程法,Macke发展了一套现代位能守恒方法,这些方法都可以计算高阶的形状[1].对于高阶形状,我们接下来拟采用很少人使用的Wavre的平衡曲面法进行研究,这种方法有很多优点:(1)类似于广义相对论中互相影响的物质分布与时空几何结构,Wavre理论可以把物理属性结构与几何分层结构完全分开.(2)Wavre理论不仅仅可以使用球谐函数,还可以使用椭圆函数,甚至放弃这些统一的特殊函数和全局坐标系,使用微分几何中的活动标架.这意味着可以抛弃传统的球谐函数(三角函数)的表示形式,例如r=R(1+c1sin2θ+c2sin4θ),而是直接用分层处的平均曲率J和分层厚度N来表示.这样做的好处是:2阶以上扁率用球谐函数表示的形状,不是流体静平衡的(Hamy-Pizzetti定理),而Wavre法配合活动标架法计算的曲面是严格平衡的.(3)Wavre理论还可以计算非流体静平衡状态下的形状,例如土星的较差自转使土星的形状发生改变,这种情况其他方法难以施展.对于高阶形状对周期等的影响，将会在接下来的工作中继续.

4 总结与展望

扁率对类地行星的自由核章动有着巨大的影响,人们可以使用观测得到的自由核章动周期来研究类地行星内部的结构.大扁率使得高速自转的类木行星的低阶本征模强烈耦合,其在球体模型与扁球体模型下计算得到的低阶本征模的频率有着很大的差距.其中不同本征模的能量分配比例会使结果有着明显的差异,但当前这方面的研究很少,本文尝试计算了不同模能量比例引起的作用.计算类木行星的低阶本征模频率,可以研究其内部结构,尤其是有丰富光环结构的土星.Fuller等[24]已经尝试用其C环中波纹变化[27]来研究土星的低阶本征模,以此来研究其内部的结构.但是Fuller等[24]采用的是一个简单的球体模型,这与实际的大扁率扁球体是有一定差距的.本研究希望将来能用扁球体模型来研究土星的C环内波纹变化与低阶本征模的关系,从而更深入地了解土星的内部结构.

[1]Moritz H.Figure of the Earth:Geodesy and the Earth’s Interior.Karlsruhe:Wichmann,1990

[2]Chandrasekhar S.Ellipsoidal Figures of Equilibrium.New Haven:Yale University Press,1969

[3]Neuberg J,Hinderer J,Zrn W.On the Complex Eigenfrequency of the“Nearly Diurnal Free Wobble”and Its Geophysical Interpretation//McCarthy D D,Carter W E.Variations in Earth Rotation.Washington:American Geophysical Union,1990:11-16

[4]黄乘利.包含海洋和大气的非刚体地球章动研究.上海:中国科学院上海天文台,1999

[5]Vorontsov S,Zharkov V.AZh,1981,58:1101

[6]Hough S S.RSPTA,1895,186:469

[7]Dziewonski A M,Anderson D L.PEPI,1981,25:297

[8]Jeffreys H.GeoJI,1926,1:371

[9]Smith M L.GeoJI,1974,37:491

[10]Smith M L.GeoJI,1977,50:103

[11]Wu X P,Wahr J M.GeoJI,1997,128:18

[12]Dahlen F A,Tromp J.Theoretical Global Seismology.Princeton:Princeton University Press,1998

[13]Rogister Y,Rochester M G.GeoJI,2004,159:874

[14]Graham G A C.Continuum Mechanics and Its Applications.Boca Raton:CRC Press,1989:797-823

[15]Seyed-Mahmoud B,Rochester M.PEPI,2006,156:143

[16]Vorontsov S V,Zharkov V N,Lubimov V M.Icar,1976,27:109

[17]Alterman Z,Jarosch H,Pekeris C L.RSPTA,1959,252:80

[18]Horedt G P.Polytropes:Applications in Astrophysics and Related Fields.Dordrecht:Kluwer Academic Publisher,2004:569-570

[19]de Pater I,Lissauer J J.Planetary Sciences.Cambridge:Cambridge University Press,2015:274-276

[20]张冕,黄乘利.天文学报,2017,58:12

[21]Zhang M,Huang C L.ChA&A,2018,42:129

[22]Marley M S.Icar,1991,94:420

[23]Gudkova T V,Zharkov V N.AstL,2003,29:674

[24]Fuller J,Lai D,Storch N I.Icar,2014,231:34

[25]Goldreich P,Kumar P.ApJ,1988,326:462

[26]Goldreich P,Kumar P.ApJ,1990,363:694

[27]Hedman M M,Nicholson P D.AJ,2013,146:12