羅昕睿，余 聰

(中山大學 物理與天文學院，廣東 珠海 519082)

多方球經常被作為模型來描述氣狀恒星的內部結構. 對于靜止的多方球，利用流體靜力學平衡方程以及多方關系可以得到萊恩-埃姆登方程，通過求解該方程可以近似描述恒星的內部結構. 但在更一般的情況下，恒星處于自轉中，因此旋轉對恒星多方球模型的外邊界與內部結構的影響尤為重要.而當考慮旋轉后，就需要對萊恩-埃姆登方程進行修改，修改后的方程從一維的常微分方程變成二維的偏微分方程.當自轉速度較慢時，此時的自轉多方球可以看成是對靜止多方球的一個偏移，所以可以利用微擾展開的方法進行求解，進而對自轉較慢的多方球的外邊界及內部結構有一個較為精準的認識.

1 自轉情形下對萊登-埃姆登方程的修改

在多方球模型中，假設壓強和密度滿足下列關系:

(1)

其中，p為內部壓強，ρ為內部密度，K為常數，n為多方指數.這個關系給出了氣狀恒星內部的壓強是如何隨著內部密度的改變而改變的，并由此結合一些基本動力學方程與熱力學方程就可以計算出恒星的內部密度分布、內部質量分布、內部溫度分布、熵輪廓等基本性質，進而使我們對恒星內部結構有一個直觀的認識.盡管這個關系是假設的，但利用該關系計算出來的結果與一些恒星的標準模型作比較卻非常接近，因此這個假設是非常合理的[8].下面我們主要利用流體靜力學平衡方程，分別討論靜止情形下的多方球滿足的方程與考慮自轉后的多方球滿足的方程.

對于靜止多方球內部達到流體靜力學平衡后，會滿足下列方程：

?p=ρ?V

(2)

?2V=-4πGρ

(3)

其中，式(3)是引力勢能V滿足的泊松方程.對徑向長度r和密度ρ進行無量綱處理，并代入上述式(2)、(3)，得到

(4)

考慮自轉后，假設自轉方向是沿著直角坐標下的z軸方向，則多方球的流體靜力學平衡方程在球坐標系(r,?,φ)中的表達式為

(5)

(6)

μ=cos ?

(7)

其中，ω為多方球的自轉速度.而式(3)也可以寫為

(8)

將式(5)、式(6)、式(8)結合在一起，可得

(9)

(10)

(11)

其中，ξ為無量綱化后的徑向長度，ρcΘn為無量綱化后的密度，v為無量綱化后的自轉速度的平方. 式(10)為自轉情形下對萊恩-埃姆登方程的修改[1-3].

2 求解自轉較慢情形下的方程

假設自轉較慢情形下方程的解是對靜止多方球的解的一個微擾[1,3]，即

Θ(ξ,μ)=θ(ξ)+vΨ(ξ,μ)+…

(12)

其中θ(ξ)為靜止情形下萊恩-埃姆登方程的解，故不依賴角度.若僅考慮v的一階項，將式(12)代入式(10)，并結合式(4)可得

(13)

利用分離變量的思路[1,3]，Ψ可以展開為

(14)

其中，ψ0、ψj都是僅依賴ξ的函數.Pj(μ)是僅依賴μ的勒讓德函數，其對應的勒讓德方程如下

(15)

將展開式(14)代入式(13)，并利用式(15)進行化簡，可得ψ0、ψj滿足的方程如下

(16)

j=1,2,3,4,…

(17)

下面需要通過計算勢能，進而確定解式(14)的每一項前面的系數Aj.

將式(8)(即球坐標下的泊松方程)寫成無量綱化的形式，并代入假設解，得

(18)

仍然利用分離變量的思路，假設上述方程的解具有如下形式：

(19)

其中，U(ξ)是靜止多方球的勢能，V0、Vj是僅依賴ξ的函數，Pj(μ)是僅依賴μ的勒讓德函數.將式(19)代入式(18)，通過比較Pj(μ)前的系數，可得

(20)

(21)

(22)

利用式(20)可得

U=Rθ+cons

(23)

利用式(16)和式(21)可得

(24)

利用式(17)和式(22)可得

(25)

進而得到式(22)的一個特解：

Vj=RAjψj+cons

(26)

但式(22)的通解還需加上下面常微分方程的解:

(27)

而式(27)的解為RBjξj，Bj為常數，故式(22)的通解為

Vj=R(Ajψj+Bjξj)+cons

(28)

綜合式(23)、式(24)、式(28)，可得最終解為

(29)

此時，式(29)將勢能V與無量綱化后表征內部壓強的Θ聯(lián)系在了一起，但仍存在一個待定系數Bj，而式(5)正好是壓強與勢能滿足的關系，可以借此確定Bj.

(30)

將式(29)代入式(30)，通過比較Pj(μ)的系數，可得

Bj=0,j≠2,

(31)

進而可得

(32)

需要注意的是，式(32)表征的是多方球內部某處的勢能，多方球外部的勢能應由下列展開式所表征，即

(33)

而多方球內外的勢能大小及勢能一階導數在多方球外邊界處應保持相等，否則勢能會出現(xiàn)不連續(xù)性.由于這里考慮的多方球自轉速度較慢，多方球外邊界可近似取為靜止多方球的外邊界，設為ξ1，則由V(ξ1)=Vext(ξ1)，V′(ξ1)=V′ext(ξ1)，可得

Aj=Cj=0,j≠2

而j=2，為

(34)

(35)

解得

(36)

進而得到最終解為

Θ=θ+

(37)

由式(16)、式(17)知ψ0(ξ)和ψ2(ξ)滿足下列的常微分方程：

(38)

(39)

式(4)、式(38)、式(39)均為常微分方程，利用四階顯式龍格-庫塔方法[6,7]可數值求解上述常微分方程.

3 求解自轉多方球的外邊界

(40)

(41)

多方球內部壓強越往外越小，故θ(ξ)是一個單調遞減的函數，所以-θ′(ξ1)=|θ′(ξ1)|，進而得到自轉情形下的多方球外邊界為

(42)

n=1,v=0.018

n=1.5,v=0.018

n=2,v=0.018圖1 自轉情形下的多方球外邊界(n ≤ 2)

n=3,v=0.001 8

n=3.5,v=0.001 8

n=4,v=0.000 18圖2 自轉情形下的多方球外邊界(3≤n≤4)

圖3 自轉情形下多方球扁率隨無量綱化自轉速度的平方的變化(n≤2)

圖4 自轉情形下多方球扁率隨無量綱化自轉速度的平方的變化(3≤n≤4)

由圖可見，n越大，自轉情形下多方球扁率受自轉速度的影響越敏感，形變越明顯.可以猜測，這主要是因為，n越大，其靜止多方球模型的密度分布越為集中，如圖5所示，導致靜止多方球的外邊界附近壓強較小，進而更容易受到由自轉帶來的離心力的影響，造成形變.

圖5 不同n的靜止多方球(無量綱)密度分布

表1 不同n的靜止多方球的及

4 自轉情形下的多方球內部密度分布

以n=2,v=0.018為例，自轉情形下的多方球沿不同角度的內部密度分布如圖6、圖7所示.

圖6 自轉多方球沿不同角度的內部密度分布(無量綱)

圖7 對圖6進行局部放大

當然，不同的n對應的密度分布差異值所取的徑向長度也會有所不同，仍然選取對應兩極方向和赤道方向的密度差異最為明顯的無量綱徑向長度ξ即可.這里僅代表n=2.

圖8 密度分布差異值與(無量綱化)自轉速度的平方的關系(n=2)

可以看到，隨著自轉速度的增加，密度分布差異值也在增加，這是符合預測的.因為自轉速度越大，多方球的赤道方向上的物質受到更強的離心力，進而往外流動，并且兩極處的物質也更容易“甩”到赤道附近，使得赤道方向上的密度分布越來越松散，而兩極方向上的密度分布越來越集中.

5 求解自轉多方球的質量半徑關系

最后，我們計算自轉多方球的質量半徑關系，利用式(43)

M=2π?ρr2drdμ

(43)

將無量綱化半徑ξ及密度Θ代入，并保留v的一階項，得到從中心到任意徑向長度ξ′的質量公式：

(44)

其中

A=θn+nθn-1vψ0(ξ)-

(45)

(46)

(47)

(48)

(49)

因此，質量半徑關系為

(50)

對于靜止多方球，即v=0時，其質量半徑關系為

M0=-Gξ′2θ′(ξ′)

(51)

利用式(50)可以改寫式(49)為

(52)

若將二者畫出來，如圖9—12所示.

圖9 多方球(無量綱化)質量半徑關系(n=1,v=0.018)

圖10 多方球(無量綱化)質量半徑關系(n=2,v=0.018)

圖11 多方球(無量綱化)質量半徑關系(n=3,v=0.001 8)

圖12 多方球(無量綱化)質量半徑關系(n=4,v=0.000 18)

可以看到，在相同中心密度的情況下，自轉多方球的質量要大于靜止多方球的質量，并且由圖可以看到，半徑越大，自轉多方球與靜止多方球的質量相差越大.這是可以被理解的，在靠近核心處，離心力的作用還較小，可以被忽略，故主要的流體靜力學平衡還是引力與壓強達到平衡，所以靜止多方球與自轉多方球的質量較為接近；而越往外，離心力的作用越來越大，當大到一定程度后，離心力的作用不可被忽略，此時向內的引力需要和向外的壓強與離心力達到平衡，故需要的更大質量來提供足夠的引力，所以越往外，自轉多方球的質量與靜止多方球的質量相差越大.

6 討論求解方法的局限性

本文利用微擾法求解自轉速度較慢的多方球時較為準確，但當自轉速度較快時，微擾法將不再適用.若仍采用微擾法進行求解，求解出來的壓強為0的等高線會沿赤道裂開，也就代表了微擾法求解的多方球外邊界也沿赤道裂開，如圖13所示.本文認為當出現(xiàn)這個狀況時，微擾法將不再適用于求解自轉多方球，應該使用其它辦法進行求解.

圖13 利用微擾法求解的自轉較快多方球外邊界(n=2)