next up previous contents
Next: 8 常微分方程式の数値解法:境界値問題 Up: 7 常微分方程式の数値解法:初期値問題 Previous: 1 Euler 法とRunge-Kutta 法

Subsections


2 ポリトロープガス球とLane-Emden方程式

常微分方程式の初期値問題の例として、ポリトロープガス球のLane-Emden方程式を解 いてみよう。

1 恒星の内部構造の基本

恒星の内部構造の釣り合いを考えよう。

2 Lane-Emden 方程式

状態方程式

\begin{displaymath}
P = P (\rho, T, X_i)
\end{displaymath} (53)

において、
\begin{displaymath}
P = K \rho^{1+{{1}\over{N}}}
\end{displaymath} (54)

として解く事を考える。これを指数 $N$ のポリトロープと呼ぶ。 すると、上の基礎方程式で、力学平衡を表す 静水圧平衡の式と連続の式の二本の式がエネルギーの式と分離する。 即ち静水圧平衡の式より
\begin{displaymath}
M_r =-{{r^2}\over{G}}{{1}\over{\rho}}{{dP}\over{dr}}
\end{displaymath} (55)

を得、これを連続の式に代入すれば
\begin{displaymath}
{{1}\over{r^2}}{{d}\over{dr}}
\left({{r^2}\over{\rho}}{{dP}\over{dr}}\right)
= -4\pi G \rho
\end{displaymath} (56)

を得る。 ここで
\begin{displaymath}
\rho=\rho_{\rm c} \theta^N
\end{displaymath} (57)

($\rho_{\rm c}$ は 星の中心密度) として、先の状態方程式を仮定すると
\begin{displaymath}
\left[
{{(N+1)K}\over{4\pi G}}\rho_{\rm c}^{{{1}\over{N}}-1}...
...\over{dr}}
\left(
r^2{{d\theta}\over{dr}}
\right)
= -\theta^N
\end{displaymath} (58)

を得る。
\begin{displaymath}
\alpha \equiv \left[
{{(N+1)K}\over{4\pi G}} \rho_{\rm c}^{{{1}\over{N}}-1}
\right]^{1/2}
\end{displaymath} (59)

$\alpha$ を定義し、
\begin{displaymath}
r \equiv \alpha \xi
\end{displaymath} (60)

で 無次元変数 $\xi$ を導入すると
\begin{displaymath}
{{1}\over{\xi^2}}{{d}\over{d\xi}}\left(\xi^2{{d\theta}\over{d\xi}}\right)
=
-\theta^N
\end{displaymath} (61)

を得る。これを 指数 $N$ の Lane-Emden equation と呼ぶ。

星の中心 $\xi=0$

\begin{displaymath}
\theta = 1
\end{displaymath} (62)


\begin{displaymath}
{{d\theta}\over{d\xi}}=0
\end{displaymath} (63)

である。また $\theta=0$ になる点が星の表面に相当する。

3 Lane-Emden方程式の解の中心近傍での振舞

Lane-Emden 方程式 (7.35) は、${\xi}=0$が、確定特異点になっている。 この様な場合、特異点の近傍での解のLaurent級数展開を考える必要がある。 いま、$\xi=0$のまわりの形式解を
\begin{displaymath}
\theta(\xi) = 1 + a_2 \xi^2 + a_4 \xi^4 + \cdots
\end{displaymath} (64)

の形に仮定して、(7.35)式の両辺に代入すれば、左辺は、
\begin{displaymath}
{\rm LHS} = 6 a_2 + 20 a_4 \xi^2 + \cdots
\end{displaymath} (65)

一方、右辺は、
\begin{displaymath}
{\rm RHS} = - ( 1 + N a_2 \xi^2 + \cdots )
\end{displaymath} (66)

となる。 両辺の${\xi^0}$$\xi^2$の係数を比較して、 $a_2 = -1/6$$a_4 = N/120$ を得る。よって、中心近傍での解は、
\begin{displaymath}
\theta (\xi) \simeq 1 - {{1}\over{6}} \xi^2 + {{N}\over{120}} \xi^4 + \cdots
\end{displaymath} (67)

と表される。

4 積分の方法

例えば、$y_1 = \theta$ $y_2 = \xi^2 d\theta/d\xi$ とおくと、 (7.35)式は次の二元連立一階微分方程式となる:
\begin{displaymath}
{{dy_1}\over{d\xi}} = \xi^{-2} y_2 ,
\end{displaymath} (68)


\begin{displaymath}
{{dy_2}\over{d\xi}} = -\xi^{2} y_1^N .
\end{displaymath} (69)

この方程式を、$\xi=0$ から外側に向かって、Runge-Kutta法で積分していけば良い。 中心は、明らかに確定特異点だから、初期値は、$\xi=0$ ではなく、中心からわずかに ずれた$\xi = h$で設定する。そこでの$y_i$の値は、上で求めた展開から得られる
\begin{displaymath}
y_1 \simeq 1 - {{1}\over{6}}\xi^2 + {{N}\over{120}} \xi^4 + \cdots ,
\end{displaymath} (70)


\begin{displaymath}
y_2 \simeq -{{1}\over{3}}\xi^3 + {{N}\over{30}} \xi^5 + \cdots
\end{displaymath} (71)

を使って求める。 積分を進めて、$y_1=0$になる点が表面である。その$\xi$$\xi_1$とすると、 星の半径は
$\displaystyle R$ $\textstyle =$ $\displaystyle \alpha \xi_1$  
  $\textstyle =$ $\displaystyle \left[ {{(N+1)K}\over{4\pi G}}\right]^{1/2}
\rho_{\rm c}^{(1-N)/2N} \xi_1$ (72)

となる。 中心から$\xi$より内側に含まれる質量は
$\displaystyle M(\xi)$ $\textstyle =$ $\displaystyle \int_0^{\alpha \xi} 4\pi \rho r^2 dr$  
  $\textstyle =$ $\displaystyle 4\pi \alpha^3 \rho_{\rm c}\int_0^{\xi} \xi^2 \theta^N d\xi ,$ (73)

或いは、(7.35)式を使って書き直すと、
$\displaystyle M(\xi)$ $\textstyle =$ $\displaystyle -4\pi\alpha^3 \rho_{\rm c}
\int_0^{\xi} {{d}\over{d\xi}}
\left(\xi^2 {{d\theta}\over{d\xi}}\right)
d\xi$  
  $\textstyle =$ $\displaystyle -4\pi\alpha^3 \rho_{\rm c} \xi^2 {{d\theta}\over{d\xi}}$ (74)

で与えられる。従って、星の全質量は
\begin{displaymath}
M = -4\pi\left[{{(N+1)K}\over{4\pi G}}\right]^{3/2}
\rho_{\...
...3-N)/2N}
\left(\xi^2{{d\theta}\over{d\xi}}\right)_{\xi=\xi_1}
\end{displaymath} (75)

となる。

5 問題

  1. 指数 $N=1.5$$N=3$、及び $N=4$ の場合について、Lane-Emden equation を Runge-Kutta 法を用いて数値的に解け。
  2. 太陽が $N=3$ のポリトロープで近似できるとして、その内部構造を求めよ。 但し太陽半径 $R_\odot = 6.96 \times 10^5$ km、 太陽質量 $M_\odot = 1.99 \times 10^{30}$ kg とし、内部は完全気体の状態方程式
    \begin{displaymath}
P_{\rm gas}= N k T
\end{displaymath} (76)

    に従うものとする。 ここで$N$ は単位体積当たりの全粒子数、$k$ はボルツマン定数。 今、水素、ヘリウム、それ以外の元素が質量比 $(X,Y,X_i)$ である 完全電離理想気体を考える。水素原子の質量を $m_{\rm H}$ で表わすとすると、 水素については核の数密度は $\rho X/m_{\rm H}$、電子の数密度は同じく $\rho X/m_{\rm H}$ 、ヘリウムについては核の数密度は $\rho Y/(4m_{\rm H})$、 電子の数密度は $\rho Y/(2m_{\rm H})$、それら以外の元素については、それぞれの 元素の atomic weight を $A_i$、atomic number を $Z_i$ とすると、 核の数密度は $\rho X_{i}/(A_{i}m_{\rm H})$、 電子の数密度は $\rho X_{i} Z_i/(A_{i}m_{\rm H})$ だから、 全粒子数 $N$
    \begin{displaymath}
N = {{\rho}\over{m_{\rm H}}}
\left( 2X + {{3}\over{4}}Y + \sum_i {{Z_i + 1}\over{A_i}}X_i \right)
\end{displaymath} (77)

    となる。 重元素については $(Z_i + 1)/A_i \sim 1/2$ だから、 水素とヘリウム以外の元素の質量比を$Z$ $(Z=1-X-Y=\sum_i X_i)$ と表わす事にすれば、
    \begin{displaymath}
N= {{\rho}\over{\mu m_{\rm H}}}
\end{displaymath} (78)

    であり、ガス圧は
    $\displaystyle P_{\rm gas}$ $\textstyle =$ $\displaystyle {{\rho}\over{\mu m_{\rm H}}} k T$ (79)
      $\textstyle =$ $\displaystyle {{\cal R}\over{\mu}}\rho T$ (80)

    ここに$\mu$
    \begin{displaymath}
\mu^{-1} = 2X +{{3}\over{4}}Y + {{1}\over{2}}Z
\end{displaymath} (81)

    で、平均分子量と呼ぶ。また ${\cal R}\equiv k/m_{\rm H} = 8.314510\times 10^3$ J$\cdot$kg$^{-1}\cdot$K$^{-1}$。 ( $k = 1.380658 \times 10^{-23}$ J$\cdot$K$^{-1}$ $m_{\rm H} = 1.672623\times 10^{-27}$ kg)。 太陽の水素、ヘリウム、金属の質量比はそれぞれ $X=0.75$$Y=0.23$$Z=0.02$ とする。
  3. 白色矮星の内部では、電子は縮退し、その縮退圧が星を支えるようになる。 電子が完全に縮退した場合、 分布関数は
    \begin{displaymath}
n_{\rm e}(p) dp = \left\{
\begin{array}{ll}
{{2}\over{h^3}}4...
...for $p < p_0$} \\
0 & \mbox{for $p > p_0$}
\end{array}\right.
\end{displaymath} (82)

    全電子数$n_{\rm e}$
    \begin{displaymath}
n_{\rm e} = \int_0^{p_0} n_{\rm e}(p) dp =
{{8\pi}\over{3h^3}}p_0^3
\end{displaymath} (83)

    だから、与えられた電子数に対しては $p_0$
    \begin{displaymath}
p_0=\left( {{3\pi^3}\over{8\pi}}n_{\rm e}\right)^{1/3}
\end{displaymath} (84)

    $p_0 \ll m_{\rm e}c^2 = 0.51 {\rm MeV}$ であれば、 $v = p/m_{\rm e}$として、電子の縮退による圧力は
    $\displaystyle P$ $\textstyle =$ $\displaystyle {{1}\over{3}}\int_0^{p_0} p v_p n(p) dp$  
      $\textstyle =$ $\displaystyle {{1}\over{3}}\int_0^{p_0} {{8\pi}\over{m_{\rm e} h^3}}p^4 dp$  
      $\textstyle =$ $\displaystyle {{8\pi}\over{15 m_{\rm e} h^3}}p_0^5$  
      $\textstyle =$ $\displaystyle {{h^2}\over{20 m_{\rm e}}}\left({{3}\over{\pi}}\right)^{2/3}
n_{\rm e}^{5/3}$ (85)
      $\textstyle =$ $\displaystyle {{h^2}\over{20 m_{\rm e}}}\left({{3}\over{\pi}}\right)^{2/3}
N_0^{5/3} \left({{\rho}\over{\mu_{\rm e}}}\right)^{5/3}$  
      $\textstyle =$ $\displaystyle 1.004 \times 10^{7}
\left({{\rho}\over{\mu_{\rm e}}}\right)^{5/3}
\quad {\rm Pa}$ (86)

    となり、$N=1.5$のポリトロープに相当する。 ここで、 $\mu_{\rm e}^{-1} = (X + {{1}\over{2}} Y + {{1}\over{2}} Z) =
{{1}\over{2}}(1+X)$。 白色矮星の状態方程式を(7.60)式とした時、その半径と質量はどれぐらいにな るか。中心密度は $\rho_{\rm c} = 10^9 {\rm kg}\cdot{\rm m}^{-3}$を目安とせよ。 また、半径の質量依存性はどうなるか。

    相対論的縮退の場合は、

    $\displaystyle P$ $\textstyle =$ $\displaystyle {{hc}\over{8}}\left({{3}\over{\pi}}\right)^{1/3}
N_0^{4/3} \left({{\rho}\over{\mu_{\rm e}}}\right)^{4/3}$  
      $\textstyle =$ $\displaystyle 1.243 \times 10^{10}
\left({{\rho}\over{\mu_{\rm e}}}\right)^{4/3}
\quad {\rm Pa}$ (87)

    となる。 状態方程式が(7.60)式で表される程の高密度の白色矮星では、 その半径と質量はどれぐらいにな るか。同じく中心密度は $\rho_{\rm c} = 10^9 {\rm kg}\cdot{\rm m}^{-3}$ を目安とせよ。 質量が中心密度によらない事を導け。 この事の意味する所を考察せよ。

6 参考文献

Hansen, C. J. and Kawaler, S. D. 1994, Stellar Interiors, (Springer-Verlag, New York).


next up previous contents
Next: 8 常微分方程式の数値解法:境界値問題 Up: 7 常微分方程式の数値解法:初期値問題 Previous: 1 Euler 法とRunge-Kutta 法
Jun Makino
平成15年4月17日