数値 的 相対論

柴田 大(東大総合文化)内容

0. イントロダクション
1. 数値的相対論において必要な要素
2. 現状：どこまで進んだのか
3. ３Ｄ計算例：連星中性子星の合体
4. 大規模・軸対称数値的相対論
５．まとめ
アニメーションは、http://esa.c.u-tokyo.ac.jp/~shibata/indexj.htmlを参照

０．イントロダクション：何故？
（Ａ） 一般論
相対論的天体−中性子星、ブラックホール、それらの連星、宇宙全体−のダイナミクスを調べたいのならば、アインシュタイン方程式を解かなくてはならない。

アインシュタイン方程式＝連立非線形方程式。

残念ながら、数値的にしか解くことが出来ない。

数値的相対論何故：Ⅱ（Ｂ） 調べる必要があるほど重要な問題があるのか？ＹＥＳい） 重力波天文学が現実化しつつある。
ＴＡＭＡ，ＬＩＧＯ（米），ＧＥＯ（独・英）が稼動中
連星の合体や星の重力崩壊により発生する重力波の波形を理論的に予測しなくてはならない

ＬＩＧＯ：Ｈａｎｆｏｒｄ＝感度上昇中

（Ｂ） 調べる必要があるほど重要な問題があるのか？続きろ） 解明すべき興味深い天体物理的問題がある。
例） ブラックホールの形成過程、ガンマ線バーストの中心エンジン＝ブラックホール＋ディスク連星中性子星の合体

ブラックホール＋中性子星

高速回転
重力崩壊

ブラックホール＋ディスクを作ってみせたい（Ｂ） 調べる必要があるほど重要な問題があるのか？続き（は） 一般相対論的現象は観測するのが難しい、珍しい- 現実的に起こっていると予想できるが、観測できない現象を計算機上で明らかにする。
- 現実には起こらないかもしれないが、一般相対論の特色が暴かれうる原理的問題を明らかにする。

１．数値的相対論において必要な要素
基礎方程式• アインシュタイン方程式• 相対論的物質場に対する発展方程式
相対論的流体（相対論的粒子、スカラー場）

G
Gµν = 8π 4 Tµνcµ∇ µTν = 0µ∇ µ ρ u = 0 ()アインシュタイン方程式の解き方
時間軸（任意）３＋１分解
時間的垂直ベクトル n３次元計量とその微分もどきを時間発展させる K ∝ γ&ijij(γ , K )

時間一定の空間面ijijiGµν n µ nν = 8π Tµν n µ nν : Hamiltonian constraint E∇
 i = 4πρ e iGµν n γ k = 8π Tµν n γ k : Momentum constraint ∇i B = 0Gµν γ iµγ νj = 8π Tµν γ iµγ νj : Evolution equation
 E& i = (∇ × B )i − 4π ji &µγ ij : 3-metric, n : timelike normal
 Bi = −(∇ × E )iµ νµ ν

Ｔ＝０で束縛方程式を解き、その後は発展方程式（双曲型）を解く。

数値的相対論に必要なコード
サイエンスを行うには必要

現実的初期条件

ゲージ条件（座標条件）
特異点を避ける、座標の歪みをただす、など
計量と物質場（流体、素過程）を発展させる
まず式の定式化、次に数値技術の開発
ブラックホール
誕生

重力波と計量から抽出
ホライズンの位置を決める
Horizon finders
ホライズン近傍で計量の勾配がきつくなる

波動帯までグリッドが必要
巨大コンピュータが必要
ホライズン周りで特別な取り扱い

２．現状：どこまで進んだのか？

Ａ. 現実的初期条件の与え方に関する現状
方程式
・重力場：スカラー、ベクトル楕円型方程式（束縛方程式その他）
・流体：第一積分の式（例えば、ベルヌーイの式）
Ａ-１：相対論的回転星（大質量星の核、中性子星）
任意の回転則に対して、平衡形状を求める必要あり。
１９８９年小松−江里口−蜂巣以降可能。
Ａ-２：連星中性子星
渦度ゼロの準平衡形状を用意する必要あり。
１９９９年瓜生、Bonazzola-Marck-Gourgoulhon以降可能。
Ａ-３：連星ブラックホール、ブラックホール-中性子星
未だにもっともらしい初期条件は存在しない。

Ｂ．アインシュタイン発展方程式に対する定式化Standard 3+1(ADM) : Evolve (γ ij , K ij )(3metric & extrinsic curvature)
12 hyperbolic equations（Ａｒｎｏｌｄ，Ｄｅｓｅｒ，Ｍｉｓｎｅｒ 1962； Ｙｏｒｋ 1979）

不安定(発散する数値的・非物理的モードが存在するため)

数値計算に適した形式
Modified 3+1 (Nakamura 87, Shibata-Nakamura 95, S 99):
 γ%ij = γ −1/ 3γ ij , γ =det (γ ij ) , Fk = δ ij ∂ iγ% jk , 1
ij−1/ 3 %
 K ij − γ ij K  K = K ijγ , Aij = γ312 hyperbolic eqs. for γ% , A%１例(ijij)+5 auxiliary nonhyperbolic eqs. for (γ , K , Fi )+ appropriate use of constraints equations and det (γ%ij ) =1
for rewriting the evolution equations
安定
その他にも、Kidder – Scheel - Teukolsky, Scheel – Lindblom，Alucubierreらなどの重要な仕事あり。世界的に活発な研究領域である。
真貝君のポスター参照。

C. 相対論的流体(1) Wilson流の昔からある方法(人工粘性を付加。)
・比較的簡単にコードが作成可能。
・星の周りに大気を置く必要なし。
・強い衝撃波が発生しなければ、まあよろしい（Ｓ９９，００）。
（例：連星中性子星の合体、中性子星の振動）
・強い衝撃波が発生する問題に対しては疑問。
（星の重力崩壊、超新星爆発）(2) 最近:特性曲線に沿って解く（例；ロー法）(Valencia（スぺ）, Munchen（独）で近年研究されてきた。
Toni Font, http://www.livingreviews.org/Articles/)
重力崩壊など、衝撃波が本質的になる場合にも有効。

特殊相対論における標準的テスト(この私の計算例では３次Ｒｏｅ法もどきを使用)
Riemann Shock Tube

N = 400， Γ = 5/3
２００点だけプロット

Wall ShockV = 0.9.
N = 400， Γ = 4/3

D. ゲージ条件（座標条件）
どのように時間方向を選んでも構わない。

垂線からのずれの自由度=シフトベクトルの自由度

自由度（シフトベクトル）

自由度（ラプス関数）

どのような形に時間一定の空間面を選んでも構わない。

垂線の距離の自由度＝ラプス関数の自由度

・ Slice : Maximal slicing（ポピュラー）-- ブラックホール形成時に特異点を避ける。
-- 座標特異点を避ける。
Ｔ

球対称計算でもチェック可能・ Spatial gauge : Minimal distortion gauge(Smarr＆York、７８)-- 座標のゆがみやねじれを発生させない。
本質的に３Ｄ問題

βkVk

Frame dragging

Vk

座標が捻じれる。

座標変換自由度（シフトベクトル）で捻じれを抑える。

t=0.0t=4.8

赤道面上での３次元計量のある成分(γxx)の時間発展 : zero shift vector
T=0T~P/6

t=8.7

t=12.9P
γ xx − 1 ≈ 2 at t ~2
T~P/3
座標ゆがみが単調に積もる.

diverge

Minimal distortion gauge（もどき）の場合
t=4.7

t=12.5

T~0

T~P/2

T~0

Metric (γxx) は単調に振動
T~3P/4t=19.4| γ xx − 1|< 0.1 at t ~ P
t=26.9

E. ３次元空間に対するApparent horizon finder
面を決める＝ r (θ,φ)に対する楕円型の方程式に帰着∆ 2 r (θ , ϕ ) = S  r ,θ , ϕ , g µν , K µν 
昔：軸上での境界条件が分からないので、どう解いていいかわからない？？
現在：そんなことはない。簡単に解ける（柴田-瓜生）。
（境界条件なし）
θφ= φ1φ= φ1+πφ=πずれたところの点F. コンピュータ資源（他力本願）
・Machine : FACOM VPP5000 (NAOJ, Vector-Parallel machine)
・ 最大 = 48 PE （720GByte）
・ 実際上 = 32 PE (500GBytes)
Total mass M
・必要なメモリー:連星の合体を例として。
r GM   rc λGW ≤ λISCO ≈ 58  2  c   7GM 2Require⇒

3/ 2 GM L ≥ λGW & ∆x ≤ 0.2  2  c  rc L≥ 290 ∆x7GM

23/ 2L >> r rc & N ≥ 580 
7GM

23/ 2

典型的グリッド：５０５＊５０５＊２５３ （赤道面対称）(〜120GByte、〜100CPU hours、〜10000ステップ）
これまで最大：６３３＊６３３＊３１７（〜250Gbyte、〜200CPU hours、〜10000ステップ）

N=500で十分なわけではないが、必要最小限度には到達。

M / R = 0.16 case(R ~13km for M=1.4)
Green：
６３３＊６３３＊３１７

Ｚ軸で観測した連星中性子星からの重力波の波形

Blue：５０５＊５０５＊２５３
Dotted Curves：
２９３＊２９３＊１４７
Red ： 2.5PN amp.

Time

まあまあの収束ＰＮ計算とも一致

必要な要素：現状のまとめ
現実的初期条件

回転星、連星中性子星＝ＯＫ
ゲージ条件（座標条件）
ＯＫ

計量と物質場（流体）を時間発展させる
ブラックホール誕生

ＯＫ
だいたいＯＫ。
重力波と計量から抽出 工夫は必要。

ホライズンの位置を決める
Horizon finders ＯＫ
ホライズン近傍で計量の勾配がきつくなる波動帯までグリッドが必要
巨大コンピュータが必要

ホライズン周りで特別な取り扱い
未完未解決問題：ブラックホールの取り扱い
ＢＨ

α

特異点にヒットしないように、時間を止める

ラプス関数

周りは時間が進む勾配がきつくなる。
いずれ発散する。

（２空間面間の距離）

Ｒ
地平線

数値計算で取り扱えない。

解決法（Ｕｎｒｕｈ）
地平線の内側を切り取ってしまえ。
切り取る座標点がなくなる境界条件は？ 安定な数値計算のための定式化や差分法は？
=> 初めに１BHありきならOK（Potsdam, Cornell, Ｉllinois）。
しかし、途中で形成する場合や２BHの場合は未解決

コードチェックについてのコメント• 既知の解析解との比較
・特殊相対論流体：
リーマンショックチューブ問題、壁衝撃波問題。
・ＧＲ（例；Ｓ９９、Fontら０２）‐ 線形重力波（Teukolsky wave）。
‐ 安定な球対称星を長時間安定に保つ。
- 球対称星の既知の固有振動（球対称、非球対称）。
- 安定な回転星を長時間安定に保つ。
- 球対称ダストの重力崩壊。
- 不安定な星の重力崩壊。
＋形成するブラックホールの質量のチェック。
• 収束性のテスト：２次のコードなら、誤差も２次で収束。
（ただし、衝撃波が存在すると精度が落ちる）

振動する球対称中性子星の進化(軸対称計算）
ハミルトン束縛の破れ
N=90(33)

120(44)180(65)Time
EOS: P = (Γ−1)ρε with Γ = 2.
シミュレーションの寿命
Μ = 1.4Μ(solar)、R=14km(N+1, N+1) grid ; N = 90, 120, 180~ 30振動周期
NS is covered by 33, 44, 65 points.

高速剛体回転する中性子星の進化

Density
at r = 0

N=240 (80)180 (60)120 (40)90 (30)Lapse
at r = 0~ 10 oscillations

KTime

M = 91% of maximum mass

動径振動の固有振動数
球対称の場合高速剛体回転の場合

fDensity (r=0)
チャンドラの公式( ６４)で計算した結果と一致

不安定な球対称中性子星のブラックホールへの重力崩壊
Mass of Apparent horizon in units of ADM mass
正しい質量に落ち着く
N=180120

90

Time
NS is initially covered by N/2 points

不安定な中性子星のブラックホールへの重力崩壊（カーパラメータ a = 0.56）
正しい質量に落ち着く

Mass of AH
N = 240 180 120.

初期条件から予想されるカーＢＨの地平線質量時間

３．３D数値相対論：私の計算例• 現在比較的容易に(講演者が)実行可能な対象：
連星中性子星の合体（次の大原さんの話も参照）
回転中性子星の動的安定性など（重力崩壊、非軸対称不安定性）• 興味はあるが手付かず：
ブラックホール−中性子星連星の合体、潮汐破壊

連星中性子星の合体 ：まず動機・合体で誕生するのは、中性子星かブラックホールか？
1.4M + 1.4M = 2.8M > 2M ,NS+NS=BH?

Maximum massSoft EOS

Stiff EOS

Spherical

~1.5M

~2.0M

~1.8M

~2.4M

>2.3M

>3.0M

Rigid Rotation+20%

Differential Rotation +>50%

・合体後、ブラックホールが誕生したならディスクは形成されるか？

？
・重力波の特徴は？特徴的周波数は？
たぶん、周波数は３ｋＨｚ程度なので、普通のレーザー干渉計で検出するのは困難。しかし、特殊なデザインの干渉計やバー型検出器のよいターゲットには？

セットアップ：現状（柴田ー瓜生、ＰＴＰ１０７， ０２）
・状態方程式 t = 0： P = K ρ^Γt > 0： P = （Γ−1)ρε： Γ = 2 & 2.25 (Here 2.25)
初期条件＝準平衡状態(瓜生、谷口) (Kerr parameter)
Compactness

Total rest mass

( M / R )∞M* TOTAL / M*Max (J=0)

Spin parameterm 2 / m1

J / M2

1.16

1.00

1

1.35

0.95

10.16

1.50

0.91

10.12 vs 0.14

1.25

0.95

0.86

0.16 vs 0.18

1.58

0.89

0.905

0.12
0.14

等質量

M*TOTAL : Total rest-mass
M*Max: Maximum rest-mass of spherical star in isolation

Note (M/R) = 0.14 & 0.16 mean
R = 15km & 13km if M = 1.4 Solar mass

M/R = 0.12 equal mass case : 重い中性子星の形成

形成された重い中性子星 ＝ 高速差動回転

Angularvelocity

剛体回転の場合のケプラー角速度

ブラックホール誕生時のディスク質量＝Negligible for merger of equal mass.

r = 6M.
Mass for r > 6M
~ 0%
Mass for r > 3M
~ 0.1%
Apparenthorizon非等質量合体におけるディスクの質量
質量比およそ0.9

r = 6M.
Mass for r > 6M
~ 6%

AlmostBH

r = 3M.
Mass of r > 3M~ 7.5%

Disk mass ~ 0.1 Solar_mass

合体後に形成される天体に関するまとめ
等質量の場合・ Low mass cases

数秒程度は生き延びるであろう重い中性子星・ High mass cases
ブラックホール。周りにディスクは形成されない。

１０％程度質量が異なる場合
・ Low mass cases
数秒程度は生き延びるであろう重い中性子星
・ High mass cases
ブラックホール。ディスクが形成されうる。

Z軸で捕らえた重力波の波形
中性子星の準固有振動に伴った重力波

crash

BH-QNM would appear

BH-QNM would appear
crash~ 2 msecブラックホールが速やかに形成する場合の予想重力波波形

未計算

合体形成物体の準固有振動

重力波の波形のフーリエスペクトルQPO of massive NS ~ 3kHz
Ｍ／Ｒ
Black ０．１２
Blue０．１４
Green ０．１６

ＢＨのＱＮＭ≒５−１０ｋＨｚ~1kHz

重力波光度M / R = 0.16

ブラックホール形成

M / R = 0.140.12

中性子星形成重い中性子星の運命
5dE−6  c • GW Luminosity≈ 10  dt
GGM 22  5GM • Potential energy W ~= 0.1Mc 2 2R Rc  5GM • Kinetic energy T ~ (0.1- 0.2)W ~ 10 Mc 2 Rc
M : Mass, R : Typical radius of massive NS-22

T−1  M   5GM ~ 10 Sec
⇒ GW emission timescale ~2 dE / dt
 3M   Rc 
重力波放出後に角運動量を失い、重力崩壊

４．大規模・軸対称数値的相対論
何故今、軸対称に回帰?
• ポツダムのグループによって、安定に長時間シミュレーションするための、新しい技法が考案された• アダプティブメッシュ法のような特殊な方法に頼らなくても、現状のコンピュータパワーならば、１０００＾２のメッシュは簡単に取ることが出来る。
• ３Ｄ計算の前に様々なテストが安価で可能。
• 重力崩壊によるブラックホールの形成、など解決されていない重要な問題が多く残されている。

計算に要するメモリー
・GR code ： 変数の数、約２００−２５０

 N var  N Memory = 2GBytes  250  1000 2・性能の良いパソコンならＮ＝１０００も十分に可能何が難しかったのか？
・ 軸対称であれば、普通は極座標か円筒座標を使う。
⇒ これらは座標特異点を持つ。
⇒ No negative r, R, θ
⇒ 差分を局所的に変える必要あり⇒ 数値不安定化し易い
ZR

さらに、正則性を保つための工夫も必要。

正則性問題 ∂ Axisymmetry ⇒ Lξ γ ij = 0 : ξ = ∂ϕ
ii

γ RR = f3 + f1R 222
γ ϕϕ = R ( f3 − f1R ) 3
γ Rϕ = R f 2
⇒ fi = f ( x, y, z ) : Regular func.
γ Rz = Rf 4γ = Rf
5 ϕz
γ zz = f 6

γ ϕϕ 
1 e.g., 2  γ RR − 2  should be regular
R R 
特殊な取り扱いが必要。数値不安定性が頻繁に発生

Cartoon法のレビュー
Y
不要

3 points

Solve equations only at y = 0
軸対称を仮定すれば同等な点.

X・ ｘ、ｙ、ｚ座標を用いる ＝ 座標特異点なし
・ 軸対称境界条件を y=+,-∆yにて付加
・ Total grid number = N＊3＊N for (x, y, z)私の計算例：回転星の重力崩壊• Parametric EOS (Following Mueller, Dimmelmeier, …)

P = PPolytrope + PThermal
PThermal = ( Γ Thermal − 1) ρε Thermal K1 ρ Γ1PPolytrope = Γ2 K 2 ρε Thermal = ε − ε Polytrope4
Γ1 ~3

Γ2 ≥ 2ρ ≤ ρ Nuc ρ ≥ ρ Nuc Γ Thermal = 1.5

Γ =4/3 & ρ ~ 1.e10 g/cc の回転平衡形状を初期条件

Collapse from a rigidly rotating
initial condition with central density ~1e10 g/cc
Density
at r = 0

At t = 0,T/W = 9.e-3
ρ (r=0) = 1.e10
M = 1.49 Solar
J/M^2 = 1.14

Lapse
at r = 0Animationis started here.

Qualitatively the same as Type I of Dimmelmeier et al (02).

重力波の波形h+2sin θ

特徴的周波数＝数１００Ｈｚ

時間

５．まとめ・ 空間3Dの一般相対論的シミュレーションは実行可能である。
科学的結果 − 例えば連星中性子星の合体後に誕生する天体や重力波の波形 − を得ることができる。
・今後の問題・コンピュータの性能に限界があるので、十分な精度の計算は未だに難しい。（精度倍にはパワー８倍要）＝〉 ＡＭＲ／ＦＭＲなどの技術開発が不可欠。
・ブラックホール形成後のシミュレーションの継続。
・ 大規模な軸対称数値的相対論が可能になった。
ターゲットは以下のような現実的問題の計算：−星の重力崩壊によるブラックホール、中性子星の形成−中性子星のアクリーション起源の重力崩壊−中性子星の相転移、など