数 値 的 相 対
論
柴田 大(東大総合文化)
内
容
0. イントロダクション
1. 数値的相対論において必要な要素
2. 現状:どこまで進んだのか
3. 3D計算例:連星中性子星の合体
4. 大規模・軸対称数値的相対論
5.まとめ
アニメーションは、http://esa.c.u-tokyo.ac.jp/~shibata/indexj.htmlを参照
0.イントロダクション:何故?
(A) 一般論
相対論的天体
-中性子星、ブラックホール、それらの連星、宇宙全体-
のダイナミクスを調べたいのならば、
アインシュタイン方程式を解かなくてはならない。
アインシュタイン方程式=連立非線形方程式。
残念ながら、数値的にしか解くことが出来ない。
数値的相対論
何故:Ⅱ
(B) 調べる必要があるほど重要な問題があるのか?YES
い) 重力波天文学が現実化しつつある。
TAMA,LIGO(米),GEO(独・英)が稼動中
連星の合体や星の重力崩壊により発生する
重力波の波形を理論的に予測しなくてはならない
LIGO:Hanford=感度上昇中
(B) 調べる必要があるほど重要な問題があるのか?続き
ろ) 解明すべき興味深い天体物理的問題がある。
例) ブラックホールの形成過程、
ガンマ線バーストの中心エンジン
=ブラックホール+ディスク
連星中性子星
の合体
ブラックホール
+中性子星
高速回転
重力崩壊
ブラックホール+ディスクを
作ってみせたい
(B) 調べる必要があるほど重要な問題があるのか?続き
(は) 一般相対論的現象は観測するのが難しい、珍しい
- 現実的に起こっていると予想できるが、観測できない現
象を計算機上で明らかにする。
- 現実には起こらないかもしれないが、一般相対論の特
色が暴かれうる原理的問題を明らかにする。
1.数値的相対論において必要な要素
基礎方程式
• アインシュタイン方程式
• 相対論的物質場に対する発展方程式
相対論的流体(相対論的粒子、スカラー場)
G
G  8 4 T
c

 T  0




   u  0 


アインシュタイン方程式の解き方
時間軸(任意)
3+1分解
時間的垂直
ベクトル n

時間一定の
空間面
ij
3次元計量とその微分もどき
を時間発展させる K  
ij
ij
, K ij 
i
G n  n  8 T n  n : Hamiltonian constraint 

E
 i  4e 


i
G n  k  8 T n  k : Momentum constraint 

i B  0
G  i j  8 T  i j : Evolution equation
 Ei  (  B)i  4

 ij : 3-metric, n  : timelike normal
 Bi  (  E )i
 
 
ji 


T=0で束縛方程式を解き、その後は発展方程式(双曲型)を解く。
数値的相対論に必要なコード
サイエンスを行うには必要
現実的初期条件
ゲージ条件(座標条件)
特異点を避ける、座標の歪みをただす、など
計量と物質場(流体、素過程)を発展させる
まず式の定式化、次に数値技術の開発
ブラックホール
誕生
重力波と計量から抽出
ホライズンの位置を決める
Horizon finders
ホライズン近傍で
計量の勾配がきつくなる
波動帯までグリッドが必要
巨大コンピュータが必要
ホライズン周りで特別な取り扱い
2.現状:どこまで進んだのか?
A. 現実的初期条件の与え方に関する現状
方程式
・重力場:スカラー、ベクトル楕円型方程式(束縛方程式その他)
・流体:第一積分の式(例えば、ベルヌーイの式)
A-1:相対論的回転星(大質量星の核、中性子星)
任意の回転則に対して、平衡形状を求める必要あり。
1989年小松-江里口-蜂巣以降可能。
A-2:連星中性子星
渦度ゼロの準平衡形状を用意する必要あり。
1999年瓜生、Bonazzola-Marck-Gourgoulhon以降可能。
A-3:連星ブラックホール、ブラックホール-中性子星
未だにもっともらしい初期条件は存在しない。
B.アインシュタイン発展方程式に対する定式化
Standard 3+1(ADM) : Evolve  ij , Kij 
(3metric & extrinsic curvature)
12 hyperbolic equations
(Arnold,Deser,Misner 1962; York 1979)
不安定
(発散する数値的・非物理的モードが
存在するため)
数値計算に適した形式
Modified 3+1 (Nakamura 87, Shibata-Nakamura 95, S 99):
1例
  ij   1/ 3 ij ,  =det  ij  , Fk   ij  i jk , 


1



ij
1/ 3 
K

K

,
A


K


K
ij
ij
ij
 ij



3





12 hyperbolic eqs. for  ij , Aij

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流の昔からある方法(人工粘性を付加。)
・比較的簡単にコードが作成可能。
・星の周りに大気を置く必要なし。
・強い衝撃波が発生しなければ、まあよろしい(S99,00)。
(例:連星中性子星の合体、中性子星の振動)
・強い衝撃波が発生する問題に対しては疑問。
(星の重力崩壊、超新星爆発)
(2) 最近:特性曲線に沿って解く(例;ロー法)
(Valencia(スぺ), Munchen(独)で近年研究されてきた。
Toni Font, http://www.livingreviews.org/Articles/)
重力崩壊など、衝撃波が本質的になる場合にも有効。
特殊相対論における標準的テスト
(この私の計算例では3次Roe法もどきを使用)
Riemann Shock Tube
N = 400, G = 5/3
200点だけプロット
Wall Shock
V = 0.9.
N = 400, G = 4/3
D. ゲージ条件(座標条件)
どのように時間方向を
選んでも構わない。
垂線からのずれの自由度
=シフトベクトルの自由度
自由度
(シフトベクトル)
自由度
(ラプス関数)
どのような形に時間一定の
空間面を選んでも構わない。
垂線の距離の自由度
=ラプス関数の自由度
・ Slice : Maximal slicing(ポピュラー)
-- ブラックホール形成時に特異点を避ける。
-- 座標特異点を避ける。
T
球対称計算でもチェック可能
・ Spatial gauge : Minimal distortion gauge
(Smarr&York、78)
-- 座標のゆがみやねじれを発生させない。
本質的に3D問題
Frame dragging
座標が捻じれる。
座標変換自由度(シフトベクトル)で捻じれを抑える。
赤道面上での3次元計量のある成分(xx)の
時間発展 : zero shift vector
T=0
T~P/6
P
 xx  1  2 at t ~
2
T~P/3
座標ゆがみが単調に積もる.
diverge
Minimal distortion gauge(もどき)の場合
T~0
T~P/2
T~0
Metric (xx) は単調に振動
T~3P/4
|  xx  1| 0.1 at t ~ P
E. 3次元空間に対するApparent horizon finder
面を決める= r (q,f)に対する楕円型の方程式に帰着
 2 r q ,   S  r ,q , , g  , K  
昔:軸上での境界条件が分からないので、どう解いていいか
わからない??
現在:そんなことはない。簡単に解ける(柴田-瓜生)。
(境界条件なし)
q
f f1
f f1
fずれたところの点
F. コンピュータ資源(他力本願)
・Machine : FACOM VPP5000 (NAOJ, Vector-Parallel machine)
・ 最大 = 48 PE (720GByte)
・ 実際上 = 32 PE (500GBytes)
Total mass M
・必要なメモリー:連星の合体を例として。
r
GW  ISCO
Require

2
 GM   rc 
 58  2 

c
7
GM



3/ 2
 GM 
L  GW & x  0.2  2 
 c 
 rc 
L
 290 

x
7
GM


2
3/ 2
L >> r
 rc 
& N  580 

7
GM


2
3/ 2
典型的グリッド:505*505*253 (赤道面対称)
(~120GByte、~100CPU hours、~10000ステップ)
これまで最大:633*633*317
(~250Gbyte、~200CPU hours、~10000ステップ)
N=500で十分なわけではないが、必要最小限度には到達。
M / R = 0.16 case
(R ~13km for M=1.4)
Green:
633*633*317
Z軸で観測した
連星中性子星からの
重力波の波形
Blue:
505*505*253
Dotted Curves:
293*293*147
Red : 2.5PN amp.
Time
まあまあの収束
PN計算とも一致
必要な要素:現状のまとめ
現実的初期条件
回転星、連星中性子星=OK
ゲージ条件(座標条件)
OK
計量と物質場(流体)を時間発展させる
OK
ブラックホール
誕生
だいたいOK。
重力波と計量から抽出 工夫は必要。
ホライズンの位置を決める
Horizon finders OK
ホライズン近傍で
計量の勾配がきつくなる
波動帯までグリッドが必要
巨大コンピュータが必要
ホライズン周りで特別な取り扱い
未完
未解決問題:ブラックホールの取り扱い
BH
a
特異点にヒット
しないように、
時間を止める
ラプス関数
周りは時間
が進む
勾配がきつくなる。
いずれ発散する。
(2空間面間
の距離)
R
地平線
数値計算で
取り扱えない。
解決法(Unruh)
地平線の内側を切り取ってしまえ。
切り取る
座標点が
なくなる
境界条件は? 安定な数値計算のための定式化や差分法は?
=> 初めに1BHありきならOK
(Potsdam, Cornell, Illinois)。
しかし、途中で形成する場合や2BHの場合は未解決
コードチェックについてのコメント
• 既知の解析解との比較
・特殊相対論流体:
リーマンショックチューブ問題、壁衝撃波問題。
・GR(例;S99、Fontら02)
‐ 線形重力波(Teukolsky wave)。
‐ 安定な球対称星を長時間安定に保つ。
- 球対称星の既知の固有振動(球対称、非球対称)。
- 安定な回転星を長時間安定に保つ。
- 球対称ダストの重力崩壊。
- 不安定な星の重力崩壊。
+形成するブラックホールの質量のチェック。
• 収束性のテスト:2次のコードなら、誤差も2次で収束。
(ただし、衝撃波が存在すると精度が落ちる)
振動する球対称中性子星の進化(軸対称計算)
ハミルトン束縛の破れ
N=90(33)
120(44)
180(65)
Time
EOS: P = (G1e with G  2.
シミュレーションの寿命
M  1.4Msolar)、R=14km
(N+1, N+1) grid ; N = 90, 120, 180
~ 30振動周期
NS is covered by 33, 44, 65 points.
高速剛体回転する中性子星の進化
N=240 (80)
180 (60)
120 (40)
90 (30)
Density
at r = 0
Lapse
at r = 0
~ 10 oscillations
K
Time
M = 91% of maximum mass
動径振動の固有振動数
球対称の
場合
高速剛体
回転の場合
f
Density (r=0)
チャンドラの公式(‘64)で計算した結果と一致
不安定な球対称中性子星の
ブラックホールへの重力崩壊
Mass of Apparent horizon in units of ADM mass
正しい質量に
落ち着く
N=180
120
90
Time
NS is initially covered by N/2 points
不安定な中性子星のブラックホールへの
重力崩壊(カーパラメータ a = 0.56)
正しい質量に
落ち着く
Mass of AH
N = 240 180 120.
初期条件から予想
されるカーBHの
地平線質量
時間
3.3D数値相対論:私の計算例
• 現在比較的容易に(講演者が)実行可能な対象:
連星中性子星の合体(次の大原さんの話も参照)
回転中性子星の動的安定性など
(重力崩壊、非軸対称不安定性)
• 興味はあるが手付かず:
ブラックホール-中性子星連星の合体、潮汐破壊
連星中性子星の合体 :まず動機
・合体で誕生するのは、中性子星かブラックホールか?
1.4M  1.4M  2.8M  2M ,
Maximum mass
Spherical
NS+NS=BH?
Soft EOS
~1.5M
Stiff EOS
~2.0M
Rigid Rotation +20%
~1.8M
~2.4M
Differential Rotation +>50%
>2.3M
>3.0M
・合体後、ブラックホールが誕生したならディスクは形成されるか?
?
・重力波の特徴は?特徴的周波数は?
たぶん、周波数は3kHz程度なので、普通のレーザー干渉計で検出するのは
困難。しかし、特殊なデザインの干渉計やバー型検出器のよいターゲットには?
セットアップ:現状(柴田ー瓜生、PTP107,‘02)
・状態方程式 t = 0: P = K ^G
t > 0: P = (G1e: G  2 & 2.25 Here 2.25
初期条件=準平衡状態(瓜生、谷口) (Kerr parameter)
Compactness
Total rest mass
 M / R 
M* TOTAL / M*Max (J=0)
0.12
1.16
等質量
0.14
1.35
0.16
1.50
0.12 vs 0.14
1.25
0.16 vs 0.18
1.58
M*TOTAL : Total rest-mass
M*Max
Spin parameter
m 2 / m1
J / M2
1.00
0.95
0.91
0.95
0.89
: 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
1
1
1
0.86
0.905
M/R = 0.12 equal mass case : 重い中性子星の形成
形成された重い中性子星 = 高速差動回転
Angular
velocity
剛体回転の場合の
ケプラー角速度
ブラックホール誕生時のディスク質量
=Negligible for merger of equal mass.
r = 6M.
Mass for r > 6M
~ 0%
Mass for r > 3M
~ 0.1%
Apparent
horizon
非等質量合体におけるディスクの質量
質量比およそ0.9
r = 6M.
Mass for r > 6M
~ 6%
Almost
BH
r = 3M.
Mass of r > 3M
~ 7.5%
Disk mass ~ 0.1 Solar_mass
合体後に形成される天体に関するまとめ
等質量の場合
・ Low mass cases
数秒程度は生き延びるであろう重い中性子星
・ High mass cases
ブラックホール。周りにディスクは形成されない。
10%程度質量が異なる場合
・ Low mass cases
数秒程度は生き延びるであろう重い中性子星
・ High mass cases
ブラックホール。ディスクが形成されうる。
Z軸で捕らえた重力波の波形
中性子星の
準固有振動
に伴った
重力波
crash
BH-QNM would appear
BH-QNM would appear
crash
~ 2 msec
ブラックホールが速やかに形成する場合の
予想重力波波形
未計算
合体形成物体の
準固有振動
重力波の波形のフーリエスペクトル
QPO of massive NS ~ 3kHz
M/R
Black 0.12
Blue
0.14
Green 0.16
BHのQNM
≒5-10kHz
~1kHz
重力波光度
M / R = 0.16
ブラックホール形成
M / R = 0.14
0.12 中性子星形成
重い中性子星の運命
5


dE
c
6
 GW Luminosity
 10  
dt
G
GM 2
2  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
-2
2
T
1  M   5GM 
 GW emission timescale ~
~ 10 
Sec

2 
dE / dt
 3M   Rc 
重力波放出後に角運動量を失い、重力崩壊
4.大規模・軸対称数値的相対論
何故今、軸対称に回帰?
• ポツダムのグループによって、安定に長時間シミュ
レーションするための、新しい技法が考案された
• アダプティブメッシュ法のような特殊な方法に頼らな
くても、現状のコンピュータパワーならば、1000^2
のメッシュは簡単に取ることが出来る。
• 3D計算の前に様々なテストが安価で可能。
• 重力崩壊によるブラックホールの形成、など解決さ
れていない重要な問題が多く残されている。
計算に要するメモリー
・GR code : 変数の数、約200-250
 N var  N 
Memory = 2GBytes 


 250  1000 
2
・性能の良いパソコンならN=1000も十分に可能
何が難しかったのか?
• 軸対称であれば、普通は極座標か円筒座標を使う。
 これらは座標特異点を持つ。
 No negative r, R, q
 差分を局所的に変える必要あり
 数値不安定化し易い
Z
R
さらに、正則性を保つための工夫も必要。
正則性問題
  
Axisymmetry  L  ij  0 :   





i
i
 RR  f3  f1R 2



2
2
   R ( f3  f1R ) 


3
 R  R f 2


 fi
 Rz  Rf 4

  Rf

5
 z

 zz  f 6

 
1 
e.g., 2   RR  2
R 
R
 f ( x, y, z ) : Regular func.

 should be regular

特殊な取り扱いが必要。数値不安定性が頻繁に発生
Cartoon法のレビュー
Y
不要
3 points
Solve equations only at y = 0
軸対称を仮定すれば
同等な点.
X
・ x、y、z座標を用いる = 座標特異点なし
・ 軸対称境界条件を y=+,-yにて付加
・ Total grid number = N*3*N for (x, y, z)
私の計算例:回転星の重力崩壊
• Parametric EOS (Following Mueller, Dimmelmeier, …)
P  PPolytrope  PThermal
PThermal   G Thermal  1 e Thermal
 K1  G1
PPolytrope  
G2
K

 2
e Thermal  e  e Polytrope
4
G1 ~
3
G2  2
   Nuc 

   Nuc 
G Thermal  1.5
G =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 = 0
Animation
is started here.
Qualitatively the same as Type I of Dimmelmeier et al (02).
重力波の波形
h
2
sin q
特徴的周波数=数100Hz
時間
5.まとめ
・ 空間3Dの一般相対論的シミュレーションは実行可能である。
科学的結果 - 例えば連星中性子星の合体後に誕生する
天体や重力波の波形 - を得ることができる。
・今後の問題
・コンピュータの性能に限界があるので、十分な精度の
計算は未だに難しい。(精度倍にはパワー8倍要)
=〉 AMR/FMRなどの技術開発が不可欠。
・ブラックホール形成後のシミュレーションの継続。
・ 大規模な軸対称数値的相対論が可能になった。
ターゲットは以下のような現実的問題の計算:
-星の重力崩壊によるブラックホール、中性子星の形成
-中性子星のアクリーション起源の重力崩壊
-中性子星の相転移、など
ダウンロード

PPT - 理論天文学宇宙物理学懇談会