Subsections

• 1 恒星内部での輻射輸達
• 2 対流によるエネルギー輸送
  • 1 対流に対する安定性
  • 2 対流によるエネルギーフラックス
  • 3 対流温度勾配
  
  
• 3 おおよその値
• 4 HR図と主系列
• 5 輻射圧と臨界光度
• 6 質量光度関係
• 7 現実的な恒星内部構造の解法の方針
• 8 ステップ１：中心からの積分
• 9 問題：中心近傍での解の振舞
• 10 ステップ２：表面からの積分
• 11 問題：表面近傍の温度構造
• 12 ステップ３：中心からの積分と表面からの積分のフィッティング

3 現実的な恒星内部構造モデル

境界値問題の例として、 恒星内部を記述する以下の方程式を現実的な仮定の下に解く事を考えよう。 基本方程式は、$M_r$を独立変数に選べば、

1. 静水圧平衡
   

   $${{dP}\over{dM_r}} = - {{GM_r}\over{4\pi r^4}}$$ (15)

   
   

2. 連続の式
   

   $${{dr}\over{dM_r}} = {{1}\over{4\pi r^2\rho}}$$ (16)

   
   

3. エネルギー保存式
   

   $${{dL_r}\over{dM_r}} = \varepsilon$$ (17)

   
   

4. エネルギー輸達
   

   $${{dT}\over{dM_r}} = -{{GM_r T}\over{4\pi r^4 P}} \nabla$$ (18)

   
   
   但しここに
   

   $$\nabla \equiv \min \left[ {{3}\over{16\pi acG}}{{\kappa L_r P}\over{M_r T^4}}, \nabla_{\rm ad} \right]$$ (19)

   
   

で与えられる。これに、補助的な関係式

1. 状態方程式
   

   $$P = P (\rho, T, X, Y)$$ (20)

   
   

2. 吸収係数
   

   $$\kappa = \kappa (\rho, T, X, Y)$$ (21)

   
   

3. 核反応生成率
   

   $$\varepsilon = \varepsilon(\rho, T, X, Y)$$ (22)

   
   

を考慮する事によって、 (8.15) - (8.18)式は、 $(P, r, L_r, T)$を変数とする四元連立常微分方程式となる。 これは 境界値問題で、 境界条件は 星の中心 $M_r =0$ で

$$r=0$$ (23)

$$L_r=0$$ (24)

星の表面 $Mr=M$ で

$$P=0$$ (25)

$$T=0$$ (26)

である。

1 恒星内部での輻射輸達

恒星内部で両面の温度が $T$ と $T+dT$ の層を考える。 この層に加わるネットの輻射圧を $dP_{\rm rad}$ とする。 輻射平衡を考えると、 この輻射圧に伴う運動量は物質に吸収される。 今単位質量単位断面積 当たりの吸収係数を $\kappa$ とすると、 単位面積単位時間当たり $F_{\rm rad}$ の輻射エネルギーが入射して $dr$ 進む間に、 $F_{\rm rad}\kappa\rho dr/c$ の運動量が吸収される。 従って、

$$-dP_{\rm rad} = F_{\rm rad} \kappa\rho dr/c$$ (27)

或いは

$$P_{\rm rad} \equiv {{a}\over{3}}T^4$$ (28)

だから

$$F_{\rm rad} = -{{ac}\over{3\kappa\rho}}{{dT^4}\over{dr}}$$ (29)

を得る。

$$L_r \equiv 4\pi r^2 F_{\rm rad}$$ (30)

として

$${{dT}\over{dr}} = {{-3}\over{4ac}} {{\kappa\rho}\over{T^3}} {{L_r}\over{4\pi r^2}}$$ (31)

を得る。

2 対流によるエネルギー輸送

1 対流に対する安定性

静水圧平衡にある層で、bubble を 断熱的に $dr$ だけ平衡位置からずらせたとする。 最初の平衡位置における bubble の圧力と密度をそれぞれ $P_0$、$\rho_0$ とする。これらは周囲の圧力、密度に等しい。 $dr$ ずれた周囲の圧力、密度は $P_0 + (dP/dr)dr$、 $\rho_0 + (d\rho/dr)dr$ である。 ずれた位置での bubble の圧力 $P^*$ はその周囲の圧力に等しい。一方、密度 $\rho^*$ は 断熱変化により決まる。即ち

$$\rho^* = \rho_0 (P^*/P_0)^{(1/\gamma)}$$ (32)

bubble が元の平衡位置に戻るかどうかは $\rho^*$ と周囲の密度の大小による。 戻るのが対流安定、どんどんずれていくのが対流不安定である。 従って、線形の範囲で、対流不安定の条件は

$${{d\ln\rho}\over{dr}} - {{1}\over{\gamma}}{{d\ln P}\over{dr}} > 0$$ (33)

2 対流によるエネルギーフラックス

対流不安定の時、bubble の温度超過分は

$$dT = \left( 1- {{1}\over{\gamma}}\right) {{T}\over{P}}{{dP}\over{dr}} dr - {{dT}\over{dr}} dr$$ (34)

$$\equiv \Delta \nabla T dr$$ (35)

単位質量当たりの熱エネルギー超過分はこれに $c_p \rho$ を掛けて与えられ、 単位面積単位時間当たりのエネルギーフラックスは更に速度を掛けたものである。

$$F_{\rm conv} = \Delta \nabla T dr c_p \rho v$$ (36)

同様に密度超過分は

$$d\rho = {{\rho}\over{T}}\Delta \nabla T dr$$ (37)

で与えられるが、移動によってなされる平均仕事量は これに重力加速度と移動距離を掛けた上で平均したものである。これが運動エネルギー の元であるから、

$${{1}\over{2}}\rho v^2 = {{\rho}\over{T}}\Delta\nabla T dr {{G M_r}\over{r^2}}{{1}\over{2}} dr$$ (38)

が成り立つ。 $dr$ として、対流で混ざる代表的長さ $l$ を採ると、エネルギーフラックスは

$$F_{\rm conv} = c_p \rho \left({{GM_r}\over{T r^2}}\right)^{1/2} (\Delta \nabla T)^{3/2} {{l^2}\over{4}}$$ (39)

で与えられる。

3 対流温度勾配

対流核の場合、温度超過は温度勾配に比べて無視出来るほど小さい。そこで 温度勾配は、

$${{dT}\over{dr}} = \left( 1-{{1}\over{\gamma}} \right) {{T}\over{P}}{{dP}\over{dr}}$$ (40)

で与えられる。

3 おおよその値

微分方程式を解く前に、解のおおよその値に目星をつけよう。 基本方程式において、左辺の微分、例えば $dP/dM_r$ を単なる割り算 $P/M_r$ に置き換える。そのうえで $M_r \sim M/2$、 $r \sim R/2$、 $\rho \sim \rho_{\rm c}$、 $T \sim T_{\rm c}$、 $P \sim P_{\rm c}$ として恒星の物理諸量の大雑把な値を求めよ。

1. 静水圧平衡
   

   $${{P_{\rm c}}\over{M}} \sim - {{GM/2}\over{4\pi (R/2)^4}}$$ (41)

   
   
   これより
   

   $$P_{\rm c} \sim {{2}\over{\pi}} {{GM^2}\over{R^4}}$$ (42)

   
   

2. 連続の式
   

   $${{R}\over{M}} \sim {{1}\over{4\pi (R/2)^2 (\rho_{\rm c}/2)}}$$ (43)

   
   
   これより
   

   $$\rho_{\rm c} \sim {{2}\over{\pi}}{{M}\over{R^3}}$$ (44)

   
   

3. 状態方程式

   ガス圧が主である場合を想定すると
   

   $$P_{\rm c} \sim {{\cal R}\over{\mu}}\rho_{\rm c}T_{\rm c}$$ (45)

   
   
   これより
   

   $$T_{\rm c} \sim {{1}\over{{\cal R}/\mu}} {{P_{\rm c}}\over{\rho_{\rm c}}}$$

   $$\sim {{1}\over{{\cal R}/\mu}} {{GM}\over{R}}$$ (46)

   
   
   ただし $G = 6.6726\times 10^{-11} {\rm N}\cdot{\rm m}^2\cdot{\rm kg}^{-2}$。
4. エネルギー輸達

   エネルギーは輻射によって運ばれていると想定すると
   

   $${{T_{\rm c}}\over{M}} \sim - {{3}\over{64\pi^2 ac}} {{\bar{\kappa}L/2}\over{(T_{\rm c}/2)^3 (R/2)^4}}$$ (47)

   
   
   これより
   

   $$L \sim {{\pi^2ac}\over{3\bar{\kappa}}} {{T_{\rm c}^4 R^4}\over{M}}$$

   $$\sim {{\pi^2}\over{3\bar{\kappa}}}{{acG^4}\over{({\cal R}/\mu)^4}} M^3$$ (48)

   
   

今、太陽を想定して $M = M_{\odot} \sim 2\times 10^{30} {\rm kg}$、 $R = R_{\odot} \sim 7\times 10^{8} {\rm m}$、 ${\cal R}/\mu \sim 1.7\times 10^8$ を代入すると

$$T_{\rm c} \sim 10^7 \quad {\rm K}$$ (49)

$$P_{\rm c} \sim 10^{15} \quad {\rm Pa}$$ (50)

$$\rho_{\rm c} \sim 4\times 10^4 \quad {\rm kg}\cdot{\rm m}^{-3}$$ (51)

を得る。光度については、

$$\bar{\kappa} \sim \kappa_0 Z (1+X) {{\rho_{\rm c}/2}\over{(... ...2)^{3.5}}} \sim 8\times 10^3 \quad {\rm m}^2\cdot{\rm kg}^{-1}$$ (52)

だから

$$L \sim 10^{24} \quad {\rm J}\cdot{\rm s}{-1}$$

$$= 10^{24} \quad {\rm W}$$ (53)

を得る。

4 HR図と主系列

$$\bar{\kappa} L \sim {{{\pi}^2}\over{6}}{{acG^4}\over{({\cal R}/\mu)^4}} M^3$$ (54)

$${{\bar{\kappa}L}\over{\bar{\kappa}L_{\odot}}} \sim \left( {{M}\over{M_{\odot}}} \right)^3$$ (55)

$$\sigma T_{\rm eff}^4 = {{L}\over{4\pi R^2}}$$ (56)

$$\left( {{T_{\rm eff}}\over{T_{{\rm eff},\odot}}} \right)^4 =... ...right)^{-2} \sim \left( {{L}\over{L_{\rm\odot}}} \right)^{1/3}$$ (57)

5 輻射圧と臨界光度

輻射圧 $P_{\rm rad}$ は

$$P_{\rm rad} = {{a}\over{3}} T^4$$ (58)

で与えられる。恒星の内部ではガス圧の他にこの輻射圧が無視できない。 温度の依存性がガス圧に比べて高いので、低密度高温状態ではむしろ輻射圧が主た る圧力になる。 輻射圧の変化は 輻射輸達の方程式から

$${{dP_{\rm rad}}\over{dM_r}} = -{{1}\over{16\pi^2 c}}{{\kappa L_r}\over{r^4}}$$ (59)

これから

$$L_r = - {{16\pi^2cr^4}\over{\kappa}}{{dP_{\rm rad}}\over{dM_r}}$$

$$= - {{16\pi^2cr^4}\over{\kappa}}(1-\beta){{dP}\over{dM_r}}$$

$$= {{4\pi cGM_r}\over{\kappa}}(1-\beta)$$ (60)

これにより、輻射光度には或る上限値がある事が判る。即ち、

$$L \leq L_{\rm crit} \equiv {{4\pi cGM_r}\over{\kappa}}$$ (61)

等号は$\beta=0$で成立するが、 $\beta=0$の場合は輻射圧だけで支える極限の場合である。

6 質量光度関係

星の光度は質量の何乗ぐらいで変わっているだろうか？ 星の質量と光度の関係を 導いてみよう。 臨界光度に近い高温大質量星の場合、 $1-\beta \simeq 1$ ゆえ

$$L \sim \left( {{L_{\rm crit}}\over{M}} \right ) M$$ (62)

となり、 光度は星の質量の1乗に比例する。

一方ガス圧が主になる中質量星、低質量星では、前に導いたように、

$$L \sim {{2\pi^2}\over{3}}{{acG^4}\over{({\cal R}/\mu)^4}}{{M^3}\over{\kappa}}$$ (63)

となり、吸収係数に依存する。中質量星の場合には吸収係数としては 電子散乱が主たる寄与をしており、$\kappa$ は 温度によらない。 従って、光度は星の質量の3乗に比例する。 一方、低質量星では 吸収係数へは Kramers の項の寄与が大きい。 そこで

$$\kappa \sim \rho T^{-3.5}$$

$$\sim {{M}\over{R^3}} \left( {{M}\over{R}} \right)^{-3.5}$$

$$\sim M^{-2.5} R^{0.5}$$ (64)

よって

$$L \propto R^{-0.5}M^{5.5}$$ (65)

となり、 光度は星の質量の5.5乗に比例する。

7 現実的な恒星内部構造の解法の方針

恒星内部を記述する ４本の 微分方程式を、中心及び表面でのそれぞれ２本ずつ計４本の境界条件とともに、数値 的に解く。 中心と表面からそれぞれ数値積分をはじめ、ある適当な所で両者が フィットするように する。 星の中心から外に向かって 積分していくとしたら二つの物理量 $(P_{\rm c}, T_{\rm c})$ がパラメーターとして残る。 星の表面から内に向かって 積分していくとしたら二つの物理量 $R, L$ がパラメーターとして残る。 解がフィットするようなパラメーター $(P_{\rm c}, T_{\rm c},R, L)$ を捜す事になる。

星としては 2 - 15 $M_{\odot}$ の主系列星を考えることにしよう。 そこで 状態方程式としては理想気体と輻射を考えて

$$P = P_{\rm gas} + P_{\rm rad}$$

$$= {{\cal R}\over{\mu}}\rho T + {{a}\over{3}} T^4$$ (66)

ここに$\mu$ は

$$\mu^{-1} = 2X +{{3}\over{4}}Y + {{1}\over{2}}Z$$ (67)

で、平均分子量。また ${\cal R}\equiv k/m_{\rm H}$ で、 $a$ は

$$a= 7.565\times 10^{-16} \quad {\rm J}\cdot{\rm m}^{-3}\cdot{\rm K}^{-4}$$ (68)

で radiation constant。 吸収係数としては Kramers + 電子散乱 を考え

$$\kappa = \kappa_0 Z (1+X)\rho T^{-3.5} + 0.019(1+X)$$ (69)

$$\kappa_0 = 4.34\times 10^{24} \quad {\rm m}^2\cdot{\rm kg}^{-1}$$ (70)

とする。 核反応生成率は CNO-cycle を考え

$$\varepsilon_{{\rm CNO}} = \varepsilon_0 \rho X X_{{\rm CNO}} T^{16}$$ (71)

$$\varepsilon_0 = 8 \times 10^{-118} \quad {\rm W}\cdot{\rm kg}^{-1}$$ (72)

とする。

8 ステップ１：中心からの積分

$M=3M_{\odot}$、$M=5M_{\odot}$、$M=7M_{\odot}$、$M=10M_{\odot}$ の星の Zero-Age Main-Sequence 時のモデルを計算してみよう。 星が Zero-Age Main-Sequence に達した時は、それに先立つ林-フェーズにおける 星全体に及ぶ対流の結果、化学組成は一様になっていると考えられる。 そこで化学組成一様な星のモデルを計算する。 $X=0.75$, $Z=0.02$ としよう。 最初のステップとして中心からの積分を行う。 二つの物理量 $(P_{\rm c}, T_{\rm c})$ がパラメーターとして残る。 積分は Runge-Kutta 法で行う。 中心近傍での物理量は

$$M_r \simeq {{4\pi}\over{3}}\rho_{\rm c} r^3$$ (73)

$$P \simeq P_{\rm c} - {{2\pi}\over{3}} G \rho_{\rm c}^2 r^2$$ (74)

$$L_r \simeq {{4\pi}\over{3}} \rho_{\rm c} r^3 \varepsilon_{\rm c}$$ (75)

$$T \simeq \left\{ \begin{array}{ll} T_{\rm c} - {{\kappa_{\rm... ...}\rho_{\rm c}^2 r^2 & \mbox{if convective} \end{array}\right.$$ (76)

パラメーター $(P_{\rm c}, T_{\rm c})$ の現実的な値としては、 前回計算を行ったポリトロープ球から得られる値を用いると良いだろう。 （後で行う表面からの積分の解と fitting point で解が フィットするようなパラメーターを捜さなければならないので、パラメーターを任意 に変えられるようなプログラムを書く必要がある 。）

9 問題：中心近傍での解の振舞

中心近傍での解の振舞を議論せよ。

10 ステップ２：表面からの積分

次のステップとして表面からの積分を行う。 二つの物理量 $(L, R)$ がパラメーターとして残る。独立変数として $M_r$ を用いるのが良いだろう。 積分は 前回同様 Runge-Kutta 法で行う。

今考えているような $M=3M_{\odot}$、$M=5M_{\odot}$、$M=7M_{\odot}$、$M=10M_{\odot}$ の星では、エ ネルギーは中心近傍では対流で運ばれ、外側では輻射で運ばれる。 そこで radiative envelope を考える。 すると

$${{dP}\over{dT}} = {{16 \pi ac GM_r}\over{3L_r}}{{T^3}\over{\kappa}}$$ (77)

だが、表面近傍では $M_r \simeq M$ $L_r = L$ なので、上の式は

$${{dP}\over{dT}} \simeq {{16 \pi ac GM}\over{3L}}{{T^3}\over{\kappa}}$$ (78)

となる。ここで opacity を

$$\kappa = \tilde{\kappa_0} P^{\alpha} T^{\gamma}$$ (79)

とすると、上の式は更に

$$P^{\alpha} dP = {{A}\over{\tilde{\kappa_0}}} T^{3-\gamma} dT$$ (80)

となるから、 この式を積分して

$$P^{1+\alpha} - P_{{\rm e}}^{1+\alpha} = {{A}\over{\tilde{\ka... ...lpha}\over{4-\gamma}} (T^{4-\gamma} - T_{{\rm e}}^{4-\gamma})$$ (81)

を得る。但し、ここで

$$A \equiv {{16 \pi ac GM}\over{3L}}$$ (82)

であり、$P$、$T$ の表面での値をそれぞれ $P_{{\rm e}}$、$T_{{\rm e}}$ とした。 もし $1+\alpha > 0$ 且つ $4-\gamma > 0$ であれば、上記 $P$-$T$ 関係は 表面での値 $P_{{\rm e}}$、$T_{{\rm e}}$ に殆ど依存しなくなる。 実際 Kramers' opacity では $\alpha = 1$、 $\gamma = -4.5$ でこの条件を満たす。 そこで、ここでは、星の表面 ($M_r = M$)では $P = 0$、 $T = 0$ とする。 Kramers' opacity を採用すると、

$$P \sim \sqrt{ {{1}\over{4.25}}{{A}\over{\tilde{\kappa_0}}} } T^{4.25}$$ (83)

となる。そこで表面領域では 解くべき方程式は

$${{dP}\over{dM_r}} = - {{GM_r}\over{4\pi r^4}}$$ (84)

$${{dr}\over{dM_r}} = {{{\cal R} T}\over{4\pi r^2\mu P}}$$ (85)

$${{dL_r}\over{dM_r}} = 0$$ (86)

$${{dT}\over{dM_r}} = -{3{\kappa_0 \mu}\over{64\pi^2 ac {\cal R}}} Z(1+X) L T^{-7.5} P R^{-4}$$ (87)

となる。 これらを $M_r = M$ から内側に積分していく。 初期値として、$M_r = M$ で $P = T = 0$。そこでの $L$ 及び $r (\equiv R)$ が 二つの自由パラメーターになる。

パラメーター $(L, R)$ の現実的な値としては、 ポリトロープ球から得られる値を用いると良いだろう。 （いずれにせよ、前回行った中心からの積分の解と fitting point で解が フィットするようなパラメーターを捜さなければならないので、パラメーターを任意 に変えられるようなプログラムを書く必要がある 。）

11 問題：表面近傍の温度構造

表面近傍では温度は

$$T \simeq {{1}\over{4.25}} {{\mu}\over{{\cal R}}} GM \left({{1}\over{r}} - {{1}\over{R}}\right)$$ (88)

で良く近似できる事を導け。

12 ステップ３：中心からの積分と表面からの積分のフィッティング

中心と表面からそれぞれ数値積分をはじめ、ある適当な所で両者が フィットするようにする。 中心からの積分は 二つの物理量 $(P_{\rm c}, T_{\rm c})$ をパラメーターとして行った。 この時、 ある適当なフィッティングポイント $($ 例えば $M_r=0.7M )$ での ４変数の値 $($ これらを $T_{\rm i}$，$P_{\rm i}$，$r_{\rm i}$，$L_{\rm i}$ と記そう。$)$ は $(P_{\rm c}, T_{\rm c})$ の関数になっている。

$$T_{\rm i} = T_{\rm i} (P_{\rm c}, T_{\rm c})$$ (89)

$$P_{\rm i} = P_{\rm i} (P_{\rm c}, T_{\rm c})$$ (90)

$$r_{\rm i} = r_{\rm i} (P_{\rm c}, T_{\rm c})$$ (91)

$$L_{\rm i} = L_{\rm i} (P_{\rm c}, T_{\rm c})$$ (92)

同様にして、 ある適当なフィッティングポイント $($ 例えば $M_r=0.7M )$ での 表面から積分してきた ４変数 $($ これらを $T_{\rm e}$，$P_{\rm e}$，$r_{\rm e}$，$L_{\rm e}$ と記そう。$)$ の値は $(L, R)$ の関数になっている。

$$T_{\rm e} = T_{\rm e} (L, R)$$ (93)

$$P_{\rm e} = P_{\rm e} (L, R)$$ (94)

$$r_{\rm e} = r_{\rm e} (L, R)$$ (95)

$$L_{\rm e} = L_{\rm e} (L, R)$$ (96)

両者が一致する必要が条件は

$$T_{\rm i} = T_{\rm e}$$ (97)

$$P_{\rm i} = P_{\rm e}$$ (98)

$$r_{\rm i} = r_{\rm e}$$ (99)

$$L_{\rm i} = L_{\rm e}$$ (100)

だが、任意に選んだ $(P_{\rm c}, T_{\rm c})$, $(L, R)$ に対してはこの条件は満たされない。この条件を満たす $(P_{\rm c}, T_{\rm c})$, $(L, R)$を $(\bar{P_{\rm c}}, \bar{T_{\rm c}})$, $(\bar{L},\bar{R})$ とし、その時の フィティングポイントでの $(T, P, r, L)$ の値を $(\bar{T}, \bar{P}, \bar{r}, \bar{L})$ とし、 前に用いた $(P_{\rm c}, T_{\rm c})$, $(L, R)$ との差を$\delta$ で表わすことにする。 今、 前に用いた $(P_{\rm c}, T_{\rm c})$, $(L, R)$ がかなり $(\bar{P_{\rm c}}, \bar{T_{\rm c}})$, $(\bar{L},\bar{R})$ に近いとすると、 フィティングの条件は、$\delta$ の一次の範囲では

$$T_{\rm i}(P_{\rm c}, T_{\rm c}) + {{\partial T_{\rm i}}\ove... ...}} \delta L + {{\partial T_{\rm e}}\over{\partial R}} \delta R$$ (101)

$$P_{\rm i}(P_{\rm c}, T_{\rm c}) + {{\partial P_{\rm i}}\ove... ...}} \delta L + {{\partial P_{\rm e}}\over{\partial R}} \delta R$$ (102)

$$r_{\rm i}(P_{\rm c}, T_{\rm c}) + {{\partial r_{\rm i}}\ove... ...}} \delta L + {{\partial r_{\rm e}}\over{\partial R}} \delta R$$ (103)

$$L_{\rm i}(P_{\rm c}, T_{\rm c}) + {{\partial L_{\rm i}}\ove... ...}} \delta L + {{\partial L_{\rm e}}\over{\partial R}} \delta R$$ (104)

となる。 ここで微分係数 ${{\partial T_{\rm i}}\over{\partial P_{\rm c}}}$ などは、フィティングポイント での $T_{\rm i}$の値のパラメーター$P_{\rm c}$に対する依存性で、これらは 数値的に調べる事が出来る。従って、上の４本の式は、 $(\delta P_{\rm c}, \delta T_{\rm c}, \delta L, \delta R)$ を４つの未知変数とする代数方程式と見做す事が出来る 。

これを解いて最初に用いたパラメーター $(P_{\rm c}, T_{\rm c}, L, R)$ に対する補正量 $(\delta P_{\rm c}, \delta T_{\rm c}, \delta L, \delta R)$ を求める。このプロセスを繰り返し、収束するまで行えば フィティングする解が求められた事になる。