中性子星の最大質量へのハイペロンの影響
平成 28 年 2 月
福井大学 工学部 物理工学科(学籍番号 12380253)
佐々木 光明目次
第 1 章 はじめに

2

第 2 章 中性子星の基礎
2.1 中性子星の天文学的な観測事実 . . . . . . . . . . . . . . . . . . . . . .
2.2 中性子星の構造 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
3

第 3 章 ハイペロンについて

5

第 4 章 相対性理論
4.1 特殊相対性理論の基礎となる仮説 . . . . . . . . . . . . . . . . . . . . .
4.2 Lorentz 変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Lorentz 収縮 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 一般相対性理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6
6
7
7
8

第 5 章 核物質の状態方程式

10

第 6 章 TOV 方程式
6.1 TOV 方程式の導出 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
11

第7章
7.1
7.2
7.3

15
15
19
20

結果および考察
中心密度と半径の関係 . . . . . . . . . . . . . . . . . . . . . . . . . . .
質量とエネルギー密度の関係 . . . . . . . . . . . . . . . . . . . . . . .
質量と半径の関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 8 章 まとめ

21

謝辞

22

参考文献

23

付 録 A MeV および fm の単位で表した物理定数

24

付 録 B 本研究で使用した 4 つの状態方程式の詳細

27付 録 C 自作 Python スクリプト
C.0.1 config.py . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.0.2 eos.py . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.0.3 tov.py . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28
28
30
32

1

第 1 章 はじめに中性子星とは超新星爆発の後にできる超高密度の天体である。超新星爆発と聞くと星の最後に起こる爆発であるとか、この爆発によってブラックホールができることがあるということは知っている人も多いかもしれない。中性子星とはその超新星爆発の後に残され、主成分が中性子であることからそう呼ばれる天体である。これは言い換
えれば一つの巨大な原子核ということができる。
この中性子星については 2010 年にこれまで観測された中でもっとも質量の大きい中性子星が発見された。これにより従来の核物質の状態方程式のモデルを見直すことになったため、原子核物理、宇宙物理の両方で今注目を浴びている天体と言える。そこ
で個人的に興味のあった宇宙物理学と配属先の研究室の専門分野である原子核物理の両方の分野にまたがった研究対象である中性子星を研究テーマとすることにする。
先述したとおり、2010 年にこれまででもっとも質量の大きい中性子星が観測されたことで、核物質の状態方程式の見直しやハイペロンの影響について議論がなされていることから本研究ではその点に注目する。ここでハイペロンとは S クォークを含むバ
リオンのことで、密度が高くなると弱い相互作用により核子がハイペロンに変わると予想されている。本研究では最大質量を TOV 方程式、状態方程式の 2 つによって求める。そしてそれぞれの状態方程式が中性子星の半径、最大質量にどう影響を及ぼすか、
またハイペロンを考慮した場合としない場合での違いについて研究を行う。

2

第 2 章 中性子星の基礎
2.1

中性子星の天文学的な観測事実

以下での中性子星の天文学的な発見の経緯について主に文献 [1] の内容を整理して記す。
光のビームを回転させ、我々から見ると点滅しているように見える星をパルサーと呼ぶ。このパルサーは今では可視光線でも点滅していることが知られている。しかし最初は可視光線ではなく電波で点滅を繰り返す星として 1967 年にケンブリッジ大学のアント二―・ヒューイッシュ(電波源のシンチレーションの研究者)、そしてその指導を受けた院生のジョスリン・ベル-バーネル (Bell-Burnell.B.) によって発見された。この最初に発見されたパルサーは 1.337 秒周期で電波を放射している。
この中性子星の存在が予言されたのはパルサー発見の 30 年以上前である。まず、1932
年にチャドウィックが原子核の反応で中性子を発見する。そして中性子が発見されるとその 2 年後の 1934 年にバーデ (Baade.W.) とツビッキ (Zwicky.F.) は中性子ばかりでできた高密度の星、「中性子星」が超新星爆発の後にできるだろうとその存在を論文のなかで予言する。1939 年にはオッペンハイマー (Oppenheimer.J) とヴォルコフ (Volkoﬀ.G.M.)
は中性子星があるとしてその質量や半径を計算、さらに中性子星の磁場や自転速度を推定し中性子星の出すエネルギーの計算をした。
このようなパルサーは宇宙全体にあると考えられる。宇宙の中では星の形成と超新星爆発が繰り返され、たくさんの中性子星が銀河の中で誕生し、たまっていく。アン
ドロメダ銀河にはまだ見つかっていないが大マゼラン銀河にはパルサーが見つかっている。私たちの銀河のなかのパルサーの分布についてみてみると、地球のまわりにたくさんあって、遠くに行くほどまばらになっている。これは遠くに行くほどパルサー
からの電波が弱くなるので見つかりにくくなるからである。

2.2

中性子星の構造

中性子星は半径が約 10[km]、質量が 1〜2[M⊙ ] ある。(M⊙ :太陽質量=1.99×1030 [kg])
恒星が超新星爆発をすることによって星の外部は吹き飛び超新星残骸となり、残った内部が圧縮されて高密度天体、すなわち中性子星になるので半径が小さく、質量が大きくなる。密度は中心で 1015 [g/cm3 ] と言われており、原子核の密度 2.8×1014 [g/cm3 ]を超えている。これに対して地球の密度は 5[g/cm3 ] しかない。
中性子星になる過程で半径が小さくなるが、もとの恒星が持っていた角運動量や磁束は保存されるので中性子星は 1 秒に 1000 回転するものや、表面の磁場が 1012 ガウス
あるものもある。これに対して地磁気は日本では平均して約 0.3 ガウスである。中性子3星はこれまでに約 2000 個ほど観測されており、その中でもっとも質量の大きい中性子星の質量は 2010 年に発見されたものでその質量は 1.97(±0.04)[M⊙ ][2] である。この中性子星が発見されるまでは約 1.44[M⊙ ] が最大質量で、核物質の状態方程式はこの値を参考に考えられていたたためこの中性子星の発見で多くの状態方程式のモデルが変更を余儀なくされることとなった。

4

第 3 章 ハイペロンについて
この章では本論文のタイトルにもあるハイペロンについて説明する。
物質を構成する最小単位である素粒子は大きく 2 つに分類されており、1 つはボソン(Boson) と呼ばれ力を媒介する粒子で、光子などを含むゲージ粒子 (gauge boson) である。もうひとつは物質の構成要素でありフェルミオン (Fermion) と呼ばれる。このう
ちフェルミオンはさらにクォーク（quark）とレプトン (lepton) に大きく分けられる。
レプトンは強い相互作用をしない素粒子で電子、ミュー粒子、各種ニュートリノ等からなり、クォークは強い相互作用をする素粒子である。クォークは 1〜3 世代によってさらに細かく分かれており、第 1 世代の u (up) と d (down)、第 2 世代の s (strange),c (charm)、第 3 世代 の b(bottom), t(top) に分類される。これら 3 個のクォークから構成される複合粒子をバリオン (baryon) と呼ぶ。そしてこのバリオンのうち第 2 世代クォークである s (strange) を 1 個以上含むバリオンのことをハイペロン (hyperons) と呼ぶ。
ここで Q を電荷（素電荷を単位として数える）、S をストレンジネス (s クォークは１単位のストレンジネスを持つ）とすると核子 (中性子、陽子) は次のような Q と S を持つ。
n (中性子)S=0 Q=0p (陽子)
S=0 Q=1
一方ハイペロンは含まれる S と Q の違いによって次の 6 種類に分けられる。
Λ（ラムダ粒子） S=-1
Q=0
Σ（シグマ粒子） S=-1
Q=-1,0,1
Ξ（グザイ粒子） S=-2
Q=-1,0これら 6 種類のハイペロンと陽子、中性子をあわせてバリオン 8 重項と呼ぶ。すべて
スピンは 1/2 である。後述するハイペロンを含む状態方程式とはこのバリオン 8 重項を考慮した状態方程式である。

5

第 4 章 相対性理論
この章では、相対性理論について説明する。中性子星は非常に重力が強いためその重力平衡をあらわす TOV 方程式は一般相対論を考慮した方程式になっている。本研究の中心課題は TOV 方程式の解法およびその解の性質を調べることなので、一般相対性
理論の理解は TOV 方程式の意味を把握するための、導出のキーポイントを知るのに必要な程度に留めておくことにした。
本章の執筆に際しては主として参考文献 [3] を参考にした。

4.1

特殊相対性理論の基礎となる仮説

相対性原理とは、ガリレオがうち立てた「物が動いているか止まっているかは、決めることが出来ない」と言う理論である。マイケルソンとモーリーによる有名な実験があり、その実験では地球自身の運動により光の媒質 (エーテル) との間に発生する相対速度を、地上のある方向と、それに垂直な方向への光の伝搬速度の違いを光の干渉を用いて測定することで決定するための装置を考案した。その結果、光の速度は光源
や観測者の運動とは無関係であることが確立した。
アインシュタインはこの問題の解決のために 2 つの仮説を設定し特殊相対性理論を導入した。
(a) すべての慣性座標系は同等である (相対性の原理)(b) 光の速度は光源や観測者の運動とは無関係に定まる (光速度不変の原理)
特殊相対論では１次元の時間と３次元の空間とを合わせ全部で４つの次元をもつ
Minkowski 時空を考え、すべての物理現象を時空における図形として認識し時間と空間とを一体として扱う座標変換、すなわち Lorentz 変換を導入する。

6

Lorentz 変換

4.2

Lorentz 変換とは行列



√ 1
v2
 1−v c2 √− c
L =  1− v22




−v
√ c
2
1− v2

0

√ 1

0

c

c

0
0

2

1− v2
c

0
0

1
0


0
  coshθ −sinhθ 0 00 
−sinhθ coshθ 0 0



=
0
1 0
  0
0
0
0
0 1
1

(4.1)

θ = tanh−1 vc
で与えられる下記のような時空の座標の一次変換

 
 
ct
ct′
x
 x′    ′ = L 
y 
y
′
z
z

(4.2)

である。(4.1) 式を書き換えて c → ∞、 vc → 0 の極限を考えると

  
t′
1
x′  −v   ′ = y   0
z′
0

0
1
0
0

0
0
1
0

 
0
t


0 x 
0  y 
1
z

(4.3)

となる。これを書き下すと t′ = t、x′ = x − vt、y ′ = y 、z ′ = z であってこの変換はNewton の古典力学で用いられる Galilei 変換と一致する。すなわち Lorentz 変換 (あるいは特殊相対性理論) が必要なのは vc が有限と考えられる場合であって vc を 0 と考えてよい場合には特殊相対性理論は不要となるのである。
座標系の Lorentz 変換で日常では考えられないようなことがいろいろと起きる。以下にその例をひとつだけ記す。

4.2.1

Lorentz 収縮

ある座標系における物体 (棒) の長さは、その両端の間の距離であり、その座標系に静止している物差しを用いて測られるべきものである。それでは棒が運動しているときには、棒の長さはどのように測ればよいか考えてみよう。
運動している棒の長さを静止した物差しで測るための具体的な手順は次のようになる。

7

『ある時刻における棒の両端の位置を同時に測り (位置に印を付け)、その 2 点間の距
離を物しで測る。』
もちろん、棒が静止していれば、この測り方は前のものと同じであることがわかるであろう。運動している物体の長さの測り方をこのようにすると、物体の長さはその速度によって異なって見えることになる。そこで座標系 S の上に長さ l の棒 AB がとまっているとする。棒が x 軸に平行であれば両端の点 A、B は t をパラメータとして (ct, x, y, z)、(ct, x + l, y, z) で与えられる 2 本の世界線を描く。世界線とは時空における質点運動の軌跡である。これを x 方向に運動する座標系 S ′ から見れば (4.1) 式を使うことにより

ct − vc x,
√v2
1 − c2x − vt−√,v2
1 − c2ct − vc x − vlc,
 √v2
1 − c2



y,

z


x + l − vt−√,
v2
1 − c2y,zとなるが、これらの世界線は新しいパラメータ t′ を使って
(
√

ct′ ,

x

1−

v2

− vt′ ,

c2

√

(

ct′ ,

(x + l)

v2

1−

c2

(4.4)

y,

(4.5)

)
(4.6)

z
)

− vt′ ,

y,

z

(4.7)

と表すこともできる。そこで S ′ では A、B 間の距離が

√

′

l =l

v2
1− 2
c

(4.8)

であり、棒は固有の長さ l より短くなっている。これを Lorentz 収縮という。棒が運動の方向に垂直に置かれている場合には長さは変わらない。

4.3

一般相対性理論

自然の世界は空間の 3 つの次元と時間の 1 つの次元を合わせた 4 次元の時空を作っている。特殊相対論については実験的な裏づけが十分にあるので局所的な時空は Minkowskiの時空として扱えるが、局所的な構造からさらに進んでいく場合には時空を広く一般の Riemann 空間 (多様体) として扱うことが必要である。
一般の時空には Minkowski 時空の慣性座標系に対応するような座標系は存在しない。
そのかわりすべての座標系が同等の立場をもっており、物理法則はいかなる座標変換を行ってもその表現の形が保たれるように共変性をもったテンソル方程式として書かれねばならない。(テンソルとは複数の成分を持ち、空間の座標変換に対していくつか

8

のベクトルの成分の積に対応した変換をうけるものをいう) これを一般相対性の原理という。
一般相対性理論について説明する。時空が Riemann 空間で表されるとして、そのRiemann 空間の距離の定義に用いられるテンソルである計量テンソル gik は重力場のポアソン方程式を拡張したアインシュタインの重力方程式

18πG
Rik − Rgik + Λ gik = − 4 Tik
2
c

(4.9)

で与えられる。R はスカラー曲率、Rik はリッチテンソル、R = g ik Rik 、Λは宇宙定数と呼ばれ、[1/(長さ)2 ] の次元を持ち、Λ−1/2 が宇宙の大きさに関係すると考えられる。
Λ gik の項 (宇宙項) はきわめて小さいがこれをゼロとしないのはアインシュタイン方程式の解の自由度を広くしておくためである。Tik はエネルギー・運動量テンソル (流体、弾性体あるいは電磁場などの場に対して時空の各点に対して定義される) である。弱い
重力の近似では、ニュートンの重力ポテンシャルを ϕ とすれば g00 ≃ 1 − 2 cϕ2 となる。
gik は時空の構造とともに重力をも記述している。質点の世界線はリーマン空間の測地線の式で与えられる。この空間の各点に慣性系が設定でき、その系での特殊相対性理論的な法則を一般座標系で書き換えることによって一般相対性理論的な電磁気学、その他の物理法則が与えられる。一般相対性理論は重力による時間の進み方の変化、太陽重力による光線の屈折など弱い重力による微小な効果として検証されている。

9

第 5 章 核物質の状態方程式
核物質の状態方程式の計算方法には大きく分けて 2 つあり、微視的理論に基づいたものと現象論に基づいたものがある。 微視的理論では核子間や、核子・ハイペロン
間の衝突実験から決めた相互作用を使って量子力学的に計算を行う。
現象論的理論では相対論的平均場模型と呼ばれるモデルの相互作用パラメータを有限核などのデータを再現するように決定し、そのモデルに基づいて計算する。
本研究では CompOSE[4] というサイトから 4 種類の状態方程式のデータを取得した。
そのうち APR、BHF-BBB2(以下 BHF と略す) という状態方程式は微視的理論に基づいて計算を行ったものでどちらもハイペロンを考慮していない状態方程式である。
残り 2 つの GM1-Y4, GM1-Y6(以下 Y4、Y6 と略す) という状態方程式は現象論に基づいて計算を行ったものでどちらもハイペロンを考慮している状態方程式である。
一般には、ハイペロンを考慮すると状態方程式が軟らかくなるため中性子星の最大質量が増えていかない。そこで本研究ではこのハイペロンを考慮した場合としなかった場合でどのような差が出るのか検証するというのが 1 つの大きな目的である。
付録 B に 4 つの状態方程式の詳細をまとめた。

10

第 6 章 TOV 方程式
この章では TOV 方程式の導出、およびその数値解析の手順について説明する。本章の執筆に際しては主として参考文献 [5]、[6] 参考にした。

6.1

TOV 方程式の導出

中性子が相対論的縮退 (フェルミ準位にある中性子の運動エネルギーが静止質量と同程度にまで増加すること) に近くなると、星を構成するガスの内部エネルギーは星の静止質量に比べ無視できなくなる。ビリアル定理により星の重力エネルギーの大きさは内部エネルギーと同程度であるからこれも静止質量に比べ無視できない。かくして中性子の重力平衡は一般相対論的に取り扱うことが必要になる。以下では球対称な星の重力平衡を考える。球対称な系の線素の表式は gik を計量テンソルとすると

ds2 = gik dxi dxk = eν c2 dt2 − eλ dr2 − r2 (dθ2 + sin2 θdϕ2 )(6.1)のように書ける。ここで，座標として物質に固定した固有座標をとれば，ν と λ は r だ
けの関数 ν(r) および λ(r) となる。この式から得られる曲率テンソルを Einstein の重力場の方程式（ここで Rik はリッチテンソル、R はスカラー曲率、G は万有引力定数、Tikはエネルギー・運動量テンソルをあらわしている。）

18πG
Rik − Rgik = − 4 Tik
2
c

(6.2)

に代入すると以下の 3 式が得られる。ε はエネルギー密度をあらわす。
′
8πG
1
1
−λ λ
ε
=
e
(
− 2) + 2
4
c
r
r
r

(6.3)

′
1
1
8πG
−λ ν
P
=
e
(
+ 2) − 2
4
c
r
r
r

(6.4)

8πG
1 −λ ′′ ν ′2 ν ′ − λ′ ν ′ λ′P=e (ν +
+
−
)
c4
2
2
r
2

(6.5)

これらの式によって星の構造も重力場も同時に定まるが、問題は一般相対論であるか
ら全空間にわたって計量テンソルを解く必要がある。ここで半径 r 内の質量を Mr とする。
∫
4π r 2
Mr = 2
εr dr
(6.6)
c 0

11

この式の両辺を r で微分すると

dMr
4π= 2 εr2
dr
c

(6.7)

となる。
また、(6.3) 式は積分でき、r=0 で λ が 0 となる式は次式となる。

e−λ = 1 −

2GMr
c2 r

(6.8)

(6.4) 式の両辺を r で積分し、(6.3)〜(6.5) および (6.7) を使って書き直せば以下の式を得る。
dP
G[M r + 4πr3 P/c2 ](ε + P )/c2=−drr[r − 2GM r/c2 ](6.9)

この球対称な星の重力平衡を記述する式は TOV 方程式とよばれる．球対称な場合の Einstein 方程式の取扱いは R.C.Tolman(1934) によってくわしく分析され,J.R.Oppenheimerと G.M.Volkoﬀ(1939) によってはじめて自由中性子からなる星に適用されたので 3 人の頭文字を取ってこう呼ばれる。
ここで文字を整理すると、半径は r、圧力は P 、エネルギー密度は ε、半径 r 内の質
量は Mr 、質量密度は ρ である。ここでバリオン 1 個あたりのエネルギーを εb 、バリオ
ン個数密度を nb とすると

ε = εb nb(6.10)

また内部エネルギーを U 、温度を T 、体積を V とすると

dU = T ds − P dV

(6.11)

より T = 0 のとき

P =−

dU
dV

(6.12)

ここでバリオンの個数を Nb とすると

Nb
V
U
εb =
Nb
と表せるからこれらを用いて (6.12) 式を変形すると
nb =

P =−

(6.13)(6.14)dUd(εb Nb )=−dVd( Nn b )
b

=−

dεb
dεb dnb=−1
dnb d( n1b )d( nb )

dεb
=−
dnb
12

d( n1b )/ dn = n dndε
b

2
b

b
b

(6.15)

となり、圧力 P は εb と nb で表すことができる。そして εb は nb で決まるので圧力 Pは nb で決まることがわかる。
(6.9) 式の変数は r, 未知関数は P 、Mr 、ε であるが、Mr は r の関数であり、(6.7) より ε が定まれば求まる。ε と P はともに nb で決まる。そこで (6.9) 式の左辺を変形しdP
dP dnb=dr
dnb dr(6.16)dnb
G[M r + 4πr3 P/c2 ](ε + P )/c2=−dP
drr[r − 2GM r/c2 ]dn(6.17)とすると

b（6.10) 式と (6.15) 式を代入すると
dεbdεbG[M r + 4πr3 (n2b dn)/c2 ](εb nb + n2b dn
)/c2
dnb
b
b
=−
dP
dr
r[r − 2GM r/c2 ]dn(6.18)bを得る。
まとめると、解くべき連立常微分方程式は (6.7) 式、(6.18) 式である。そして (6.7) 式でdεb
dP
用いる ε、(6.18) 式で用いられる ε、 dn、 dn
は nb の関数として状態方程式より与えらbbれる。

図 6.1: TOV 方程式に現れる量の依存関係
図 6.1 は TOV 方程式を解く際に現れる r の関数とその一階微分の依存関係についてまとめたものである。青色の矢印は状態方程式、緑色の矢印は TOV 方程式、赤色の矢
印は合成関数の微分法、茶色の矢印は (6.7) 式の関係を表す。
ある r の値における nb と Mr の値が与えられたとき、最終的にはその r の値にお
b
r
ける dn
と dM
の値が計算でき、それを使って r + ∆r における nb と Mr がそれぞれ
dr
dr

13

r
b
nb + dn∆r と Mr + dM
∆r として数値的に求まるのである。数値計算を行う際には付
drdr
録 C に記載したプログラムを実行し、計算方法として Runge-Kutta 法を採用した。

14

第 7 章 結果および考察
7.1

中心密度と半径の関係

2max massnb [1/fm3]

1.5

1

0.5

0
0

2

4

6

8

10

12

14

R[km]
図 7.1: APR の状態方程式を使用した場合の中心のバリオン個数密度と半径の関係
図 7.1 は APR という状態方程式での中性子星の中心のバリオン個数密度と半径の関係を表した図である。r=0 での nb の値を与えると、TOV 方程式により密度や圧力が動径座標 r の関数として定まる。圧力がゼロになる r がその解の表す中性子星の表面である。その r より内側の質量の総和がその中性子星の質量である。このように中心での nb の関数として中性子星の半径と質量が定まるのである。
赤い太線で描かれているのが、この状態方程式で計算した際に最大質量を示したバ
リオン個数密度と半径の関係である。このグラフから最大質量時以上のバリオン個数密度を持つときの解 (nb =1.5,2.0) のときは半径が最大質量時よりも小さくなっている
ことがわかる。
中心での密度が上がっているにも関わらず、質量も半径も小さくなっているというこ

15

とは中性子星の内側だけで重力によって潰れていくことを示している。このことから最大質量のとき以上のバリオン個数密度を持っている解は重力平衡的に不安定であり、このまま重力崩壊が進んでいくと中性子星はブラックホールになってしまうと考えられる。

2
max massnb [1/fm3]

1.5

1

0.5

0
0

2

4

6

8

10

12

14

R[km]
図 7.2: BHF の状態方程式を使用した場合の中心のバリオン個数密度と半径の関係
図 7.2 に BHF、図 7.3 に Y4、図 7.4 に Y6 で計算した場合の同様の図を表示する。図7.1 で見られた解の特徴は BHF、Y4、Y6 の状態方程式で計算した際にも見ることができた。このことから最大質量を超えると重力平衡的に不安定になるという解の特徴はTOV 方程式の解に共通する一般的な特徴であると考えられる。

16

2
max massnb [1/fm3]

1.5

1

0.5

0
0

2

4

6

8

10

12

14

R[km]
図 7.3: Y4 の状態方程式を使用した場合の中心のバリオン個数密度と半径の関係
図 7.3 ではほかのグラフには見られない半径 6.5km あたりでの線の屈曲が見られた。
本研究では状態方程式の計算方法にまで研究対象を広げていないが、この屈曲は状態方程式にハイペロンを考慮したことで得られたものであると考えられる。半径が 6.5km
程度より大きいときは傾きが小さく、半径が 6.5km より大きいときは傾きが大きい。これは半径が 6.5km あたりより内側でハイペロンが現れ始めるためだと考えられる。ハイペロンを考慮するとエネルギーが上がらなくなるため、ハイペロンが現れなくなる
半径以上では状態方程式が硬くなるようにパラメータを設定し、質量が増えるようにモデルのパラメータを調節してあるものと考えられる。

17

2
max massnb [1/fm3]

1.5

1

0.5

0
0

2

4

6

8

10

12

14

R[km]
図 7.4: Y6 の状態方程式を使用した場合の中心のバリオン個数密度と半径の関係
図 7.4 は図 7.3 と同様にハイペロンを考慮したにも関わらず Y4 で見られた線の屈曲があまり見られなかった。後述するが Y6 の状態方程式では 2M⊙ を超え、Y4 の状態方程式では超えなかった。このことから Y6 の状態方程式は 2M⊙ を超えるようにパラメーターを設定した、つまりハイペロンの影響があまり出ないように設定された状態方
程式であるため図 7.3 のような線の屈曲が見られなかったのではないかと考えられる。

18

7.2

質量とエネルギー密度の関係

2.5

M [solar mass unit]

2

1.5

1

APR (no hyp)
BHF-BBB2 (no hyp)
GM1-Y4 (with hyp)
GM1-Y6 (with hyp)max observed mass

0.5

0
14

14.5

15

15.5

16

3

log(epsilon)[g/cm ]図 7.5: 質量とエネルギー密度の関係
図 7.5 は 4 つの状態方程式での質量と中心のエネルギー密度の関係を表したものである。Ｙ 6 のグラフを見るとハイペロンを考慮した状態方程式であるにも関わらず最大
質量が 2M⊙ を超えた。ハイペロンを考慮した状態方程式は軟らかくなる、つまりエネ
ルギーが上がらなくなるので最大質量が低くなると予想していた。だが実際にはパラメータを適切に決定した状態方程式を使用することで、観測された最大質量を十分超えることが分かった。
またこの図からはハイペロンを考慮していない状態方程式である APR、BHF とハイペロンを考慮した状態方程式である Y4、Y6 では最大質量に達するときのエネルギー密度に違いが見られた。ハイペロンを考慮すると最大質量に達するときのエネルギー
密度が低くなるという傾向が見られた。

19

7.3

質量と半径の関係

2.5

M [solar mass unit]

2

1.5

1

APR (no hyp)
BHF-BBB2 (no hyp)
GM1-Y4 (with hyp)
GM1-Y6 (with hyp)max observed mass

0.5

0
4

6

8

10

12

14

16

R[km]
図 7.6: 質量と半径の関係
図 7.6 は 4 つの状態方程式での質量と半径の関係を表したものである。
この図からもハイペロンを考慮していない状態方程式である APR、BHF とハイペロンを考慮した状態方程式である Y4、Y6 では最大質量に達するときの半径に違いが見られた。ハイペロンを考慮すると最大質量に達するときの半径が大きくなるという傾向が見られた。
この図において例えば APR の状態方程式のグラフを見ると半径 10.5km のとき質量の異なる点が 3 点ほど取れることになる。これは TOV 方程式を解く際には、r=0 でのnb の値を与えることで中心での nb の関数として中性子星の半径と質量が定まるので、同じ半径でありながら質量も中心のバリオン個数密度の値も違う解が複数得られることもあるためである。この 3 点はどれも APR という状態方程式で計算した中性子星が再現しうる質量と半径の関係であるが、最大質量を超えた解、APR の状態方程式においては半径 10km 以下の解に関しては重力平衡的に不安定な解となっている。

20

第 8 章 まとめ
本研究では TOV 方程式を状態方程式と組み合わせて解く方法について詳しく考察し、さらに解の性質を調べることを行った。文献や参考書には TOV 方程式と状態方程式を実際に組み合わせる手順が書かれておらず、結果のみ記載しているものが多かった
ので本研究ではその過程に着手した。CompOSE というサイトから状態方程式のデータを取得し、付録 C に記載しているプログラムを組むことで数値計算した。
その結果から分かったことはハイペロンを考慮した状態方程式でもパラメータを適切に決定することで実際に観測された最大質量を再現することができるということ、またハイペロンを考慮した状態方程式と、していないもので最大質量に達するときの中
心のエネルギー密度と半径にそれぞれ違いが見られるということがわかった。さらに最大質量を超えた解は重力的に不安定であるという性質も確認することができた。
今後の課題としては今回は状態方程式を 4 つしか使用しなかったので、検討する状態方程式を増やしていくということが挙げられる。

21

謝辞
本論文を作成するにあたり田嶋直樹先生には終始変わらぬ親切丁寧なご指導賜りましたことに心より感謝し、御礼申し上げます。貴重な御時間を割いてまで根気強く丁寧に御指導頂いたこと、本論文作成に関しても多大な御協力を頂きました。さらには研究室では伊藤研人先輩、杉浦友章先輩にもご指導、ご意見を頂きましたことを厚く御礼申し上げます。同じゼミの草間究君にも御礼申し上げます。誰よりも遅くまで学生室に残り頑張る姿に刺激をもらい、私も最後まで頑張ることができました。また卒業論文発表会において貴重なご意見を頂いた物理工学科の先生方にも心より御礼申し
上げ、謝辞の言葉とさせて頂きます。
平成 28 年 2 月
佐々木光明

22

参考文献
[1] 柴田晋平 「宇宙の灯台 パルサー」 恒星社厚生閣 (2006).
[2] P.B. Demorest et al., Nature, 467, 1081 (2010).
[3] 平川浩正「相対論」共立出版株式会社 (1971)[4] http://compose.obspm.fr/[5] 林忠四郎他「岩波講座 現代物理学の基礎 [第 2 版]11 巻 宇宙物理学」岩波書店 (1978)[6] 高原文郎「新版 宇宙物理学」朝倉書店 (2015)

23

付 録A

MeV および fm の単位で表した物理定数SI 単位系（国際単位系）でのいくつかの物理定数の値を下記に記す。
h̄ = 1.054571596(82) × 10−34 [J s](A.1)

c = 2.99792458 × 108 [m s−1 ](A.2)

e = 1.602176462(63) × 10−19 [C](A.3)µ0 = 4π × 10−7 [N A−2 ](A.4)

ϵ0 = 1/µ0 c2 [F m−1 ](A.5)

mn c2 = 939.565379[MeV](A.6)

エネルギー単位の eV(電子ボルト) の定義は e[C] の電荷が 1[V] の電位差から得る運動エネルギーであるから

1[eV] = 1.602 × 10−19 [C] × 1[V](A.7)

が成り立つ。[C][V]=[J] なので

1[eV] = 1.602 × 10−19 [J]

(A.8)

1[eV]
= 1[無次元数]
1.602 × 10−19 [J](A.9)

が成立し、従ってとなる。
これらの値を使って、MeV および fm を単位とする場合の h̄c の値、1g の物質の静止
エネルギーの値、万有引力定数 G の値を以下で求める。計算過程は有効数字４桁で示すが、最終的な結果の値は数字を丸めることなく計算して得たものである。

24

(1)h̄c の値

h̄ = 1.054 × 10−34 [J s] ×

1[eV]
1.602 × 10−19 [J]= 6.580 × 10−16 [eV s]= 6.580 × 10−22 [MeV s]h̄ = 6.58211889(26) × 10−22 [MeV s]∴(A.10)

h̄c = 6.58 × 10−22 [MeV s] × 3.00 × 108 [m s−1 ]= 19.74 × 10−14 [MeV m]= 197.4[MeV fm]∴

h̄c = 197.3269601[MeV fm](A.11)(2)1g の物質の静止エネルギーの値
1g の物質の静止エネルギーを eV 単位で表すと
1g × c2 = 10−3 [kg] × (2.998 × 108 )2 [m2 s−2 ]= (2.998 × 108 )2 × 1013 [kg m2 s−2 ]= (2.998 × 108 )2 × 1013 [J] ×1[eV]
1.602 × 10−19 [J]= 5.610 × 1032 [eV]∴

1g × c2 = 5.6095892 × 1032 [eV](A.12)(3) 万有引力定数 G を fm、MeV、c で表す
G = 6.674 × 10−11 [m3 kg−1 s−2 ]= 6.674 × 10−11 [J m kg−2 ](A.9) 式より
1 [J MeV−1 ] =11.602 × 10−13

G = 6.674 × 10−11 [J m kg−2 ]6.674=
× 102 [MeV m kg−2 ]
1.602

25

(A.13)

ここで kg を c、MeV で表す

1[kg] = 1[kg c2 c−2 MeV MeV−1 ]=1[kg c2 ][MeV c−2 ]1[MeV]
2(3.0 × 108 )2 ( ms2 kg)=[Mev c−2 ]1[MeV]
J= (3.0 × 108 )2 [][MeV c−2 ]MeV
9.0 × 1029=
[MeV c−2 ]1.602
この値を (A.13) 式に代入し

6.674
9.0 × 1029
× 102 )([MeV c−2 ])2 [MeV m]
1.602
1.602
6.674 × (1.602)2=
× 10−56 [MeV−1 c4 m]
1.602 × (9.0)2
6.674 × 1.602
× 10−41 [MeV−1 c4 fm]=(9.0)2G=(= 1.324 × 10−42 [fm MeV−1 c4 ]∴G = 1.32374 × 10−42 [fm MeV−1 c4 ]26

(A.14)

付 録B

本研究で使用した 4 つの状態
方程式の詳細

温度はゼロ度、成分比は β 平衡での値に固定したものを使用している。クラスト (星殻) の非一様な構造は他の結果を使用している。

1. APRA18 相互作用等を使った微視的な変分法の結果考慮した粒子： 中性子、陽子、電子、ミュー粒子
A. Akmal, V.R. Pandharipande and D.G. Ravenhall, Phys. Rev. C, 58, 1804 (1998).

2. BHF-BBB2
核子間相互作用として Paris ポテンシャル等を使った微視的な Bueckner HartreeFock 計算結果考慮した粒子： 中性子、陽子、電子、ミュー粒子
M. Baldo, I. Bombaci and G. F. Burgio, Astron. Astrophys. 328, 274 (1997).
3,4. GM1-Y4, GM1-Y6
相対論的平均場模型という現象論的模型による計算結果
考慮した粒子： 中性子、陽子、電子、ハイペロン
・N.K. Glendenning and S.A. Moszkowski, Phys. Rev. Lett. 67, 2414 (1991)
・ M. Oertel, C. Providencia, F. Gulminelli, A. Raduta, arxiv.1412.4545. (2014)(2 solar mass の中性子星の発見後、改訂した)

27

付 録C
C.0.1

自作 Python スクリプト

config.py

このプログラムでは後のプログラムで使用するために記号や変数を定義している。
#!/cygdrive/c/Python34/python.exe#!/usr/bin/env python3"""
usage:
from config import G,M_solar_kg,neutron_mass,c,MeV_to_solar_mass,fm_to_km"""#from math import pi
c_SI= 2.99792458e+8# 光の速さ [m/s] (SI 単位系)
e_SI
= 1.6021766208e-19 # 電気素量 [C] (SI 単位系)
G_SI
= 6.67384e-11# 重力定数 [m^3/kg*s^2] (SI 単位系)
Ms_SI
= 1.9884e+30# 太陽質量 [kg]m_to_km= 1.0e-3# =(1 m)/(1 km), x[m]=x*m_to_km[km]fm_to_km = 1.0e-18# =(1 fm)/(1 km), x[fm]=x*fm_to_km[km]fm_to_cm = 1.0e-13# =(1 fm)/(1 cm), x[fm]=x*fm_to_cm[cm]MeV_to_eV = 1.0e+6# =(1 MeV)/(1 eV), x[MeV]=x*MeV_to_eV[eV]gram_SI= 1.0e-3# 1 グラム [kg](SI 単位系)

mn_MeV=939.565379# (自由中性子の質量)*(光速度)^2 (MeV)

G=G_SI/c_SI**2*Ms_SI*m_to_km# (重力定数)/(光速度)^2 [km/solar_mass]

MeV_to_solar_mass=MeV_to_eV*e_SI/(Ms_SI*c_SI**2) # =(1 MeV)/(1 solar_mass)*c^2
MeV_to_gram = MeV_to_eV*e_SI/(gram_SI*c_SI**2) # =(1 MeV)/(1 g)*c^2
ed_from_MeV_fm_to_g_cm = MeV_to_gram/fm_to_cm**3# ={(1 MeV)/(1 fm)^3}/{(1 g)*c^2/(1 cm)^3}

if __name__ == ’__main__’:print(’{:.10e} =c : speed of light [m/s]’.format(c_SI))
print(’{:.10e} =e : elementary charge [C]’.format(e_SI))
print(’{:.7e} =G : Gravitational constant [m^3/kg*s^2]’.format(G_SI))
print(’{:8e} = mass of the Sun [kg]’.format(Ms_SI))
print(’{:.9e} =(mass of a free neutron)*c^2 (MeV)’.format(mn_MeV))

28

print(’{:.6e} =G : (Gravitational constant)/c^2 [km/solar_mass]’.format(G))
print(’{:.10e} =1[MeV]/1[solar_mass]*c^2’.format(MeV_to_solar_mass))
print(’{:.10e} =1[MeV/fm^3]/1[g/cm^3]*c^2’.format(ed_from_MeV_fm_to_g_cm))

29

C.0.2

eos.pyこのプログラムでは状態方程式から導かれる関数を定義している#!/cygdrive/c/Python34/python.exe#!/usr/bin/env python3#!/cygdrive/c/Python34/python.exe
import numpy as np
import scipy as sp
from scipy.interpolate import InterpolatedUnivariateSpline#config より mn_MeV を取り込む
from config import mn_MeV # 中性子質量 m_nc^{2}=939 [MeV]
def eos_functions(eos_table_file_name):with open(eos_table_file_name,’r’) as fin:nbs=[0.0] # n_B : バリオン個数密度 [1/fm^3]
prs=[0.0] # 圧力 [MeV/fm^3]
eds=[0.0] # エネルギー密度 [MeV/fm^3]
for i,line in enumerate(fin):#print(i,line)
t=line.split()
if len(t) != 5:
print(’skip this data line:’,t)continue
nbs.append(float(t[1]))prs.append(float(t[3]))eds.append(float(t[4]))print(’# eos_functions: EOS data have been read from {:s}’.format(eos_table_file_name))f_pr = InterpolatedUnivariateSpline(nbs,prs)"""
圧力 (MeV/fm^3) はバリオン個数密度 [1/fm^3] の関数"""
f_eps_b_tilde = InterpolatedUnivariateSpline(nbs,eds)"""
tilde{epsilon_b} = epsilon_b/(mn_MeV *c^2)-1
epsilon_b = intenal energy per baryon [MeV]"""
def f_ed(nb):’’’’’’’’
中心のエネルギー密度はバリオン個数密度の関数
internal energy density[MeV/fm^3] as a function ofnb(baryon number density [1/fm^3])"""
energy_per_baryon=(f_eps_b_tilde(nb)+1)*mn_MeV
return energy_per_baryon*nb30

f_pr_der = f_pr.derivative()"""
dp/dn_{b}はバリオン個数密度の関数
d (pressure)/d baryon_number_density [MeV]
as a fucntion of baryon_number_density [1/fm^3]"""
return f_ed,f_pr,f_pr_der
if __name__ == ’__main__’:f_ed,f_pr,f_pr_der=eos_functions(’eos.table’)with open(’eos.g’,’w’) as fout:
for i in range(0,1501):
nb=i*0.0001# baryon number density
pr=float(f_pr(nb)) # pressure
ed=float(f_ed(nb)) # energy density
pr_der=float(f_pr_der(nb)) # d(pressure)/d(n_B)
pr2=f_pr_der.integral(0.0,nb)#print(’{:24.16e} {:24.16e} {:24.16e}’.format(nb,pr,ed),file=fout)
print(’{:16.7e} {:16.7e} {:16.7e} {:16.7e} {:16.7e}’.format(nb,pr,ed,pr_der,pr2),file=fout)"""with open(’eos.g’,’w’) as fout:
for nb,pr,ed in zip(nbs,prs,eds):
print(’{:24.16e} {:24.16e} {:24.16e}’.format(nb,pr,ed),file=fout)"""
print(’# eos.g has been created.’)31

C.0.3

tov.pyこのプログラムは上記 2 つのプログラムで定義した値をインポートし、実際に数値計算するためのものである。一番下の行で与えるバリオン個数密度の個数と範囲を変更できる。
#!/cygdrive/c/Python34/python.exe#!/usr/bin/env python3
from config import G, MeV_to_solar_mass, fm_to_km, ed_from_MeV_fm_to_g_cm#config というファイルからインポートする
from eos import eos_functions#eos というファイルから eos_functions をインポートする
import numpy as np
import scipy as sp
from scipy.integrate import ode
from math import pi
def create_f():fct1=4*pi*MeV_to_solar_mass/fm_to_km**3
def f(r,nb_M):"""
inputs:variable : radius r [km]list of the values of two functions at r:
nb_M[0] : n_b (baryon number density) [1/fm^3]nb_M[1] : M (mass of the part of the neutron star within radius < r)[Solar mass]
outputs: list of the values of two functions at r:
returned_object[0] : d n_b/d r [1/fm^3 km]
returned_object[1] : d M /d r [Solar mass/km]
working variables:
nb : baryon number density [1/fm^3] = y[0]M : mass of neutron star within radius < r [Solar mass]p : pressure [MeV/fm^3], = eos.f_pr(nb) (from EOS)
ed : energy density [MeV/fm^3], = eos.f_ed(nb) (from EOS)
dp_dnb : dp/d nb : [MeV], = eos.f_pr_der(nb) (from EOS)
dp_dr : dp/dr=-G(M+4*pi*r**3*f_pr)(f_ed+f_pr)/(r*(r-2G*M)) (TOV eq.)
[MeV/fm^3 km]"""
nb=nb_M[0]M=nb_M[1]ed=f_ed(nb)p=f_pr(nb)dp_dnb=f_pr_der(nb)if r>0.0:
dp_dr=-G*(M+fct1*r**3*p)*(ed+p)/(r*(r-2*G*M))else:
dp_dr=0.0
dM_dr=fct1*r**2*ed

32

dnb_dr=dp_dr/dp_dnb
return [dnb_dr,dM_dr]
return f"""
def check_f():
for M in [0.1,0.2,0.5,1.0,2.0]:
for nb in [0.0, 0.03,0.1,0.3,1.0]:
for r in [0.0, 1.0, 2.0, 5.0, 10.0]:a,b=f(r,[nb,M])
print(r,nb,M,a,b,’ :A’)#check_f()#exit(1)"""eos_data_file_name=’eos.APR.table’#読み込む状態方程式のファイルをここで指定する
f_ed, f_pr, f_pr_der = eos_functions(eos_data_file_name)f=create_f()solver=ode(f)solver.set_integrator(# ’vode’,method=’Adams’,with_jacobian=False,’dopri5’,# Runge-Kutta 4-5th orderatol=1.0e-10, # tolerance(error) = rtol*|y|+atol
rtol=1.0e-10)
def solution(nb_center):def print_each_r(r, nb, M):
p=float(f_pr(nb))
ed_CGS=float(f_ed(nb))*ed_from_MeV_fm_to_g_cm
print(’{:f} {:f} {:e} {:f} {:e} :P’.format(r, M, ed_CGS, nb, p))
rs.append(float(r))
nbs.append(nb)Ms.append(M)ps.append(p)eds_CGS.append(ed_CGS)
return p
print(’# n_b (r=0) = {:f} ;’.format(nb_center),’r[km], M[M_sun], ed[g/cm^3], nb[1/fm^3], p[MeV/fm^3]’)
r_ini,r_fin,dr=0.0, 100.00, 0.05
# unit=[km]nb_M_ini=[nb_center, 0.0] # n_b[1/fm^3], M[Solar mass]
solver.set_initial_value(nb_M_ini,r_ini)r, (nb, M) = r_ini,nb_M_ini
p=float(f_pr(nb))
ed_CGS=float(f_ed(nb))*ed_from_MeV_fm_to_g_cm

33

rs,nbs,Ms,ps,eds_CGS=[],[],[],[],[]
p=print_each_r(r,nb,M)
while solver.successful and solver.t < r_fin:delta_r=dr
if r<0.05:
delta_r=min(dr,0.001)
if nb<0.001:
delta_r=min(dr,0.0001)
elif nb<0.003:
delta_r=min(dr,0.0002)
elif nb<0.01:
delta_r=min(dr,0.0005)
elif nb<0.03:
delta_r=min(dr,0.001)
elif nb<0.1:
delta_r=min(dr,0.002)
elif nb<0.3:
delta_r=min(dr,0.005)
solver.integrate(solver.t+delta_r)r=solver.tnb,M=solver.yp=print_each_r(r,nb,M)
if p < 0.0:
breakwith open(’tov_APR_nb-r.dat’,’a’) as fout:#出力するファイル名をここで変更する
i=-2
print(’{:f} {:f} {:f} {:e} {:e} {:e} {:e} {:s} :S’.format(nbs[0],Ms[i],rs[i],eds_CGS[0],ps[0],nbs[i],ps[i],eos_data_file_name),file=fout)
print(’# done’)print()
for i in range(1,2401,1):
nb=i*0.001
solution(nb)’’’’’’’’’’’’’’
この値を変えることで計算する量を変える。この場合は 0.001〜2.400 までの 2400 点計算する
solution の () 内を任意の数字にすることでも計算でき、その場合は指定した数字のときバリオン個数密度のみ計算する’’’’’’

34