羅昕睿，余 聰

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

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

1 自轉(zhuǎn)情形下對萊登-埃姆登方程的修改

在多方球模型中，假設(shè)壓強和密度滿足下列關(guān)系:

(1)

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

對于靜止多方球內(nèi)部達到流體靜力學(xué)平衡后，會滿足下列方程：

?p=ρ?V

(2)

?2V=-4πGρ

(3)

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

(4)

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

(5)

(6)

μ=cos ?

(7)

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

(8)

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

(9)

(10)

(11)

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

2 求解自轉(zhuǎn)較慢情形下的方程

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

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

(12)

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

(13)

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

(14)

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

(15)

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

(16)

j=1,2,3,4,…

(17)

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

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

(18)

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

(19)

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

(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為常數(shù)，故式(22)的通解為

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

(28)

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

(29)

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

(30)

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

Bj=0,j≠2,

(31)

進而可得

(32)

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

(33)

而多方球內(nèi)外的勢能大小及勢能一階導(dǎo)數(shù)在多方球外邊界處應(yīng)保持相等，否則勢能會出現(xiàn)不連續(xù)性.由于這里考慮的多方球自轉(zhuǎn)速度較慢，多方球外邊界可近似取為靜止多方球的外邊界，設(shè)為ξ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]可數(shù)值求解上述常微分方程.

3 求解自轉(zhuǎn)多方球的外邊界

(40)

(41)

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

(42)

n=1,v=0.018

n=1.5,v=0.018

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

n=3,v=0.001 8

n=3.5,v=0.001 8

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

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

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

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

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

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

4 自轉(zhuǎn)情形下的多方球內(nèi)部密度分布

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

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

圖7 對圖6進行局部放大

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

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

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

5 求解自轉(zhuǎn)多方球的質(zhì)量半徑關(guān)系

最后，我們計算自轉(zhuǎn)多方球的質(zhì)量半徑關(guān)系，利用式(43)

M=2π?ρr2drdμ

(43)

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

(44)

其中

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

(45)

(46)

(47)

(48)

(49)

因此，質(zhì)量半徑關(guān)系為

(50)

對于靜止多方球，即v=0時，其質(zhì)量半徑關(guān)系為

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

(51)

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

(52)

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

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

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

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

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

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

6 討論求解方法的局限性

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

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