目次✓第１回目 高エネルギー天体物理の基礎「星の進化、超新星爆発の標準理論」✓第２回目 シミュレーション研究最前線「爆発天体現象のエンジンは解明されたか？」✓第３回目 ニュートリノ輻射流体数値計算法
ガンマ線バーストのエンジン✓第４回目 マルチメッセンジャー天文学に向けて（重力波・ニュートリノ・多波長電磁波観測）今日の前半部のお題§３-2 ニュートリノ輻射＋流体計算法
流体計算（数値流体解法）
ニュートリノ輻射輸送計算§ 3-3 ガンマ線バーストの中心駆動源

§３-１
超新星の動的（ダイナミカル）な進化を追うための数値輻射流体力学

内容• 流体力学基礎方程式の導出（done）• ランキン・ユゴニオ(Rankine-Hugoniot)
関係• Riemann 問題• ニュートリノ輻射輸送流体力学基礎方程式のまとめ

未知数

１本

３１
運動量の式

３１

エネルギーの式

１

１ 状態方程式

１１Polytropic
ポリトロープ状態方程式

１１gi=-∇iΦ

１ ポアソン方程式
方程式の数と未知数の数が同じで方程式系が閉じている。

「方程式の数と未知数の数が同じで方程式系が閉じている」⇒初期条件を与えれば、その後の時間発展が（原理的には）分かるはず。

☆ただ方程式が非線形で、複雑。
解析解はごく限られた問題でしか得られない。

したがって、流体の運動を追うためには、数値計算が不可欠。
☆更に偏微分方程式を数値的に解かなくてはならない。
いかに精度良く、数値的に流体の時間発展を追うか？
数値流体力学（Computational Fluid Dynamics:CFD)

内容• 流体力学基礎方程式の導出• Riemann問題に向けて☆Rankine Hugoniot関係導出☆具体例：Strong shock３次元シミュレーション例

密度、速度、圧力

ある瞬間、各点各点で密度、圧力、速度が不連続に与えられたとき、その後の各物理量の時間発展を数値的に追うこと、これがCFDのココロ。

もっと単純化すると、二つの状態の異なるガスをいれた容器のふたを取ったとき、その後のガスのダイナミクスが分かれば良い。

？

このような問題設定をRiemann問題と呼ぶ。

Riemann問題を解くための準備（１）
一次元保存系の流体方程式

コンパクトに以下のように書く。

Riemann問題を解くための準備（１）に対して、を定義すると、と書ける。
任意の関数Φにたいして以下の関係が成り立つ。

部分積分して以下の関係が得られる。（無限遠でGはゼロ）領域Ωを二つの領域に分けて考える。

部分積分する。

グリーンの法則より、面積分に変換。

不連続面の速度を用いて不連続面の法線方向のベクトルはと書けるを思い出すと、nと内積を取って（衝撃波静止系）では以下の関係が成り立つ、

注：

衝撃波
Rankine-Hugoniot関係と呼ぶ。

３番目のエネルギーの式から移項して、１番目の質量保存の式からMを使って、２番目に代入すると１番目の質量保存の式から、比体積を定義してこの曲線を衝撃波断熱曲線、
Rankine-Hugoniot断熱曲線と呼ぶ。
内部エネルギーは状態方程式を用いて、（ここでは圧力、密度依存とする）二番目の流束保存に代入すると、Hugoniot曲線は二つのパラメーターで書ける。

Hugoniot 曲線から読み取れる重要な性質（１）初期の状態が分かっているとき、衝撃波背面の物理量はただ一つのパラメータにより指定される。
例えば、ｐ２が与えられると、τ2は断熱曲線より与えられ、そうすると、Mがわかり、そうすると、u1,u2も分かる。
（２）とを結ぶ直線（Rayleigh直線）の傾きが、流束Mを与える。

曲線曲線☆Question!
初期の状態から、Hugoniot
曲線に載っている点だったら、どんな点にでも状態が遷移することが可能か？ NO!

エントロピーは以下のように書ける。

衝撃波通過後はエントロピーが増加すべし。

１（衝撃波前）

２（衝撃波後方）

衝撃波
さらにだから、より、M2＞０より、τ1>τ2,密度は衝撃波背面で上がるべき。
だったので、であることが分かる。

内容• 流体力学基礎方程式の導出• Riemann問題に向けて☆Rankine Hugoniot関係の導出☆具体例：Strong shockStrong shock(強い衝撃波）の場合で静止（V1 = 0）しているガスに、速度Dで衝撃波が飛び込んできたとする。この場合衝撃波背面の物理量を求めたい。

今、?

V1 = 0

衝撃波静止系に移る

u2 = D – v2に注意！

u1 = Du1 = D?

u2 = D – v2に注意！

RH関係の３つの式は以下と同値。
音速
マッハ数＞１ 超音速（supersonic)＜１ 亜音速(subsonic)

今、強い衝撃波の場合を考える。

☆この式の導出
宇宙流体力学（坂下志郎、池内了）など参照＞＞１

V2は実験室系、u2は衝撃波静止系の量であることに注意！！！

衝撃波背面の量が、衝撃波前面の量と衝撃波の速度で表せた。
（強い衝撃波問題が解けた）Riemann問題

Riemann問題？

未知数勘定CD
左右の密度（２）、速度、圧力、RSの速度
CDの速度、FSの速度(7未知数）⇔RH
Reverse shock(RS)
7本（内訳
RS：3
FS:3
ＣＤ：1)

接触不連続面（圧力、速度連続）

圧力
速度密度(CD)

Forwardshock(FS)

初期条件によって解の構造が異なる：Shock tube テスト

密度不連続

速度不連続

圧力不連続

エネルギー不連続もっとも一般的な場合
青線：解析解

密度連続

速度不連続

圧力連続

エネルギー連続

密度連続

速度連続

圧力不連続

エネルギー不連続Adaptive mesh refinement 法(適合格子法）の一例（Kifonidius et al. 2003 A&A)

差分の式にする。

差分の式にする。

これらを求めるために、リーマン問題を解く！
精密・面倒

Godnov法（厳密解）Roe 法（振幅の滑らかな変化を無視）
HLLD・C法（fast/slowを区別しない,alfven)
HLL法（中間状態を一つだけ）
簡単・粘性大§3-2
超新星シミュレーションにおけるニュートリノ輸送計算法

多次元シミュレーションでニュートリノ加熱駆動爆発をおこすのに必要な要素✓ニュートリノ輻射パート :Neutrino heating depends onneutrino luminosities, spectra,
and angular distributions.

heatingHeating-dominatedCooling-dominatedcooling

PNS

Stalled shock~200km※ちなみに

SN環境では、
ニュートリノはmassless
E = cpf (t,r,, , E, p ,  p )E R (t,r,, , E)   d p d p f
E R (t,r,,  )   dE d p d p f✓流体パート :“MGMA”(6 dimensional problem)“MG”(Multi energy-Group：エネルギー群）) or IDS(isotropic diffusion source approximation)“Gray (no energy-dependence)”（詳細については第３回目で）自転, 磁場を正確に扱うには3D計算が不可欠!
✓ニュートリノ分布関数

ごく中心部を除いてCourtesy:Multidimensional neutrino transport coupled to 3D MHD hydro(hopefully with general relativity) is required. 熱分布(Fermi-Dirac分布)
住吉光介さんではない。
⇒
少なくとも、Energy依存性を残した
輻射輸送計算が不可欠(multi-energy group輸送）（エネルギー多群輸送）解くべきはBoltzmann equation

左辺：ニュートリノ数の変動
右辺：衝突項(Collisional term)反応によるニュートリノ数の変動✓超新星コードにマイクロ物理を入れるロードマップ

Easy(1)Simple deleptonization：「Ye処方」（Ｙｅ（ρ）で電子捕獲を扱う）(2)Neutrino leakage scheme：「ニュートリノ漏れ出し法」（ニュートリノ冷却のみ近似的に扱う）(3)Simple deleptonization + Light-bulb treatment：「ライトバルブ近似」（ニュートリノ加熱を手で入れる）(4)Single energy flux-limited diffusion:「グレイ輸送」(5) Multi-energy flux-limited diffusion:「マルチグループ輸送+FLD」
Isotropic Diffusion Source Approximation（ＩＤＳＡ）(5’) Ray-by-ray Boltzmann transport :「マルチグループ輸送+α」(6) Multi-angle Boltzmann transport：「マルチエネルギー・マルチアングル輸送」
Hard

解くべきはBoltzmann equation✓超新星コードにマイクロ物理を入れるロードマップ
6 months
Only if your hydro is robust

Easy
Simple deleptonization

ニュートリノクーリングのみ入る

Neutrino leakage scheme

Simple deleptonization + Light-bulb treatment
Single energy flux-limited diffusion
Multi-energy flux-limited diffusion,
Isotropic Diffusion Source Approximation

ニュートリノ加熱が「手」で入る。

✓ニュートリノ加熱・冷却がパラメータなしに入る✓ ただ陰的解法が不可欠（implicit scheme）

Ray-by-ray Boltzmann transport

Multi-angle Boltzmann transport
Hard> 1 year (or a life work!)✓Ｙｅ 処方
球対称ボルツマン輸送計算から得られたYe(ρ）のfitting formula を使う。
(Pros) Very easy to implement.
(often employed in GR simulations:Dimmelmeier et al. (2007), Ott et al. (2008))(Cons) バウンス前しか使えない。
親星が違う時の妥当性。
球対称のシステムのみ有効。
✓Neutrino leakage scheme (NLS) (van Riper 1981, Bludman et al. 1981,

Rosswog & Liebendoerfer 2003, Kotake et al. 2003)

Iron CoreNeutrino Sphere・ Outside the neutrino sphere, neutrinos are assumed to vanish
instantaneously.
assuming
・ Inside the neutrino sphere, the beta equilibrium is assumed.

✓Ｙｅ 処方
球対称ボルツマン輸送計算から得られたYe(ρ）のfitting formula を使う。
(Pros) Very easy to implement.
(often employed in GR simulations:Dimmelmeier et al. (2007), Ott et al. (2008))(Cons) バウンス前しか使えない。
親星が違う時の妥当性。
球対称のシステムのみ有効。
✓Neutrino leakage scheme (NLS) (van Riper 1981, Bludman et al. 1981,

Rosswog & Liebendoerfer 2003, Kotake et al. 2003)

Iron Core
Neutrino Sphere・ Outside the neutrino sphere, neutrinos are assumed to vanish
instantaneously.
assuming

Cons: The neutrino heating cannot be treated.
・ Inside theimportantneutrino sphere,
the beta equilibrium is assumed.
Pros : The NLS can reproduce
featuresobtained in a Boltzmann simulation, like evolution ofneutrino energy, luminosity, etc, to some extent.
The scheme is valid only before the neutrino heating becomes
important (~50 ms after bounce).

Roadmap to implement the supernova microphysics to your code!
6 months
Only if your hydro is robust

Easy
Simple deleptonization

The neutrino coolingcan be taken into account.

Neutrino leakage scheme
Matthias parametrization + Light-bulb treatment
Single energy flux-limited diffusion
Multi-energy flux-limited diffusion,
Isotropic Diffusion Source Approximation
Ray-by-ray Boltzmann transport

The neutrino heating,but only parametrically.

✓The neutrino heating/coolingwithout any parameters.
✓ Implicit schemes are needed
to treat accurately
the matter-neutrino coupling
as well as to solve theBoltzmann equation !

Multi-angle Boltzmann transport
Hard> 1 year (or a life work!)✓Matthias’ recipe + light-bulb method(Scheck et al. 04, Ohnishi et al. 07,
Murphy & Burrows 08, Nordhaus et al. 10)

ライト・バルブ法;

M

原子中性子星から照らされるニュートリノ光度を手で与える。
(type I simulationでよくつかわれる）StandingshockLν
PNS
１０ｋｍ
Neutrino temperature
given by hand.

(Pros): Qualitative effects of neutrino heatingon convection or hydrodynamic instability
can be studied.

(Cons): The neutrino heating is completely
an input parameter !

Roadmap to implement the supernova microphysics to your code!
6 months
Only if your hydro is robust

Easy
Simple deleptonization

The neutrino coolingcan be taken into account.

Neutrino leakage scheme
Matthias parametrization + Light-bulb treatment
Single energy flux-limited diffusion
Multi-energy flux-limited diffusion,
Isotropic Diffusion Source Approximation
Ray-by-ray Boltzmann transport

The neutrino heating,but only parametrically.

✓The neutrino heating/coolingwithout any parameters.
✓ Implicit schemes are needed
to treat accurately
the matter-neutrino coupling
as well as to solve theBoltzmann equation!

Multi-angle Boltzmann transport
Hard> 1 year (or a life work!)

Why implicit scheme（陰解法） is needed ?
absorption

j: emissivity
λ： absorptivityemission

Note that

If solved explicitly,,,

The implicit treatment is stable !

ボルツマン方程式の右辺：「衝突項」の計算法

Emission and absorption<O(100) MeV＝
等エネルギー散乱 (M_{nuc}~1GeV>>Eν）

非弾性散乱 inelastic

ペア反応

(Bruenn 1985,ApJS)
nnu_eemission

j: emissivity
λ： absorptivity
The matrix element:

See e.g., Bjorken & Donellep

(Bruenn 1985,ApJS)
nnu_eemission

j: emissivity
λ： absorptivity
The matrix element:

eSee e.g., Bjorken & Donell✓ボルツマン方程式の右辺は頑張って計算できる。
✓ボルツマン方程式を解くためには、もうひとつのキーポイントが・・

pボルツマン方程式の階層構造
Intensityの角度モーメント量✓ゼロ次 ⇒ radiation energy(scalar)
外力が加わっていない場合✓１次⇒ radiation flux（輻射フラックス）(vector)✓2 次⇒ radiation pressure (輻射圧）

Specific intensityを定義する。

(tensor)ボルツマン方程式の角度モーメント量輻射が単位立体角、単位エネルギー、ある方向、dA、dtを通過する確率になるどこまで行っても閉じない。

✓ Closure することが必要(1) Flux limited diffusion法（流速制限法）(2) M1 closure(3) Variable Eddington factor法

ボルツマン方程式の解法(1/3)

Bruenn et al. 1985, Livne et al. 2005,
Kotake et al. 2006

✓流速制限法（Flux-Limited-Diffusion)
輻射は等方的と近似、等方からのずれの一次を入れる

Fickの法則から以下を仮定
Bowers&Wilson(1982)

Flux limiter

Obergaulinger and Janka (2010)

ボルツマン方程式の解法(2/3)✓M1 closure (1次のmomentの式まで解く）
階層性を閉じるためを仮定する。
Eddington tensorを以下の形に仮定する。

Freestreaming
limit

自由伝
搬

Semitransparent

遷移領域

Diffusionlimit

拡散極限（Audit et al. 2002)

Eddington factor
Flux factor☆Diffusion limitを取ると、fν →0pν →1/3, Pν(輻射圧）→ 1/3 Eν☆streaming limitを取ると、fν →1pν →1, Pν(輻射圧）→ Eν

FLD、M1共に二つの極限状態を内挿したフォーマリズムボルツマン方程式の解法(3/3)

VE; Rampp & Janka (1998)
Sn: Livne et al. (04), Hubney & Burrows(06)✓Variable Eddington factor method✓Rayにそって実際に輸送方程式を解き、Ｉを決める(ray-method)✓Intensityが求められたら、諸量が分かる✓SN法 (Spectral ordinates method)

Momentを取らず、直接ボルツマン方程式を角度方向まで差分を取って解く。
（formalismはstraightforwadだが、とにかく演算量が大変）ボルツマン方程式の解法(3/3)

ApJ (2009)

VE; Rampp & Janka (1998)
Ray-trace 計算の一例
Sn: Livne et al. (04), Hubney & Burrows(06)✓Variable Eddington factor method✓Rayにそって実際に輸送方程式を解き、Ｉを決める(ray-method)✓Intensityが求められたら、諸量が分かる✓SN法 (Spectral ordinates method)

Momentを取らず、直接ボルツマン方程式を角度方向まで差分を取って解く。
（formalismはstraightforwadだが、とにかく演算量が大変）ニュートリノ輸送計算法の長所と短所
Transparent
Diffusive regime Semi-transparent
regimeregime

Boltzmannsolver(VEF,SN)

Insufficientangularresolution

Truncation errors
in flux2成分近似ボルツマン formalism (Basel-Tokyo)
Flux-limited
Diffusion,
M1 closure
Ray-tracing
methodflux-factor ismodel-dependent

Suffer fromshort mean freepath

flux-factor
is unknown

Limited byreaction

The ideal argorthm combines the three orange fields !

２成分近似ボルツマン輸送法
The neutrino distribution function “f”
を２成分に分解する

Liebendoerfer et al. (09), ApJ

流体素片:trapped part – opaque region:streaming part – transparent region
Boltzmann equation is then ,The two components are evolved
separately ,- is the (diffusion) term,
which converts trapped into streaming part and vice verse.
-is determined between the free-streaming and the β-equilibrium limit via akind of flux limiter.

まとめ：超新星シミュレーション

imMultiEnergy

3D

TYPEII (Full simulations)

2D

1D

Trsan

D
rtop

nsois
en

BoltzmannBoltzmann
MGFLD

Ott+08

果てしなき旅は続く・・

Marek+09

Burrows+07

2CB
Suwa et al. +10
Liebendoerfer et al. 09

GrayOhnishi+07

Iwakami+08,09

TYPE I:

Nordhaus

Adiabatic

Blondin+03

1D

2D

Blondin+07

3D
Spacial Dimensions

吉報：3D,GRみんな
爆発には良さそう?

§３-3
超新星とガンマ線バースト（以降の話の筋）✓ガンマ線バースト（Gamma-Ｒay Burst):
ガンマ線がバースト的に観測される天体現象✓過去４0年に渡って起源が不明✓超新星との相関が報告される（10年前）
普通の超新星じゃない極超新星✓極超新星のエンジンの理解⇔恒星進化論の統一的解明

観測のまとめ>ミリsec

最も明るい天体～1051erg/sec

光度

赤方偏移

GRB
～1000 events/yr
等方に到来、宇宙論的起源
平均エネ～200 keV, 非熱的
継続時間10-3s～103s :

残光
Ｘ線
光学
電波時間

フォトンカウント数

ＧＲＢの時間変動

Fishman&Meegan(95)

時間LongShort

GRBの継続時間

Long burst、Short burstの二極分布それぞれタイムプロファイルが全然違う。ソースは何？

色々あるのですが現在の理解
Long
Short

GRBの継続時間

中性子星合体uncertain

超新星爆発

certain謎の天体現象GRBと始めて相関が観測された天体：超新星 だった！
Galama et al (1998) Nature

イベントレート（per galaxy)

RIbc～ 10-2 yr-1
RGRB～0.5(250) Gpc-3yr-1
SN1998bwが観測されたMatsubayashi et al.
06

超新星とGRBが無相関とすると、10年で同期する確率：
GRB980425 観測のerror circle

Rcomb～ 10-７

一方、数年の観測で相関が実際観測された。

GRBと超新星は優位な相関がある！

SN1998bw ： Hypernova(極超新星）
SN 1998bw

SN 1987A=大問題：1051 ergの爆発すら難しいのにどうやって十倍に？
E ~ 30×1051ergs

E ~ 1×1051ergs

GRB（Long)と相関がある超新星は一般にenergetic

Nomoto et al. (2002)

ふつうの超新星と極超新星は何が違うのか？
ちがうブランチに移る物理は何なのか？

「ブラックホール・回転・磁場」が鍵e -e +
Meszaros&Rees’９２νν

アニメーション

エッセンス：狭い領域に莫大なエネルギー注入✓ニュートリノ✓磁場

コラプサー・シナリオ「Collapsar」：failed supernova

ニュートリノ駆動ガンマ線バースト
Goodman 1987,Asano & Fukuyama 2000 (ApJ)θ✓加熱率∝(1 – cos θ)
円盤の対称性を考えても軸状が一番温められる✓ ニュートリノ加熱率

ニュートリノ加熱率
反応率

断面積：

Energetics 的には十分
ただし降着円盤の温度、半径に非常に敏感
運動学：
分布関数：

⇒ブラックホール周りの時空での構造を決める
一般相対論的流体計算⇒ニュートリノ加熱率
マルチアングル（フルの）ニュートリノ輻射輸送

最近の結果✓フルの一般相対論的計算＋ニュートリノ漏れ出し法

Sekiguchi & Shibata(2011)

重力崩壊⇒ブラックホール形成⇒降着円盤の成長

1054 erg!
明るさは十分

落下時間

落下時間/加熱時間✓特殊相対論計算＋Ray-tracing ニュートリノ加熱

Harikae et al. (2010)

磁場駆動ガンマ線バースト
エネルギー源：ブラックホールの回転エネルギー
Blandford-Znajek process (MNRAS,179,433,1977)
Kerr BHから引き抜ける回転エネルギー

無次元Kerr parameter
Maximally rotating Kerr BH（Schwarzshild半径を光速で回転する場合）

BZ processの継続時間(Thorne, Membrane パラダイムより）

回転しているBHの時空：Kerr 時空と ＢＺプロセス
H.K. Lee et al, astroph/9906213

Boyer-Lindquist 座標

Ergosphere:

Kerr parameter
K= 0の条件からでは、non-static observer.
Frame-dragging とよぶBHの角運動量Jと質量Mの比

Event horizon:となる条件から

PhD thesis by A. Mueller

Blandford-Znajek process のOrder評価
・Kerr BHから引き抜ける回転エネ

・Kerr BH の表面積
from・Kerr BH のentropy(面積定理）
無次元Kerr parameter・Kerr BH のirreducible mass

Maximally rotating Kerr BH（Schwarzshild半径を光速で回転する場合）

・BZで引き抜ける回転エネはのefficiencyがかかる。

・Kerr BH の全質量
BZ processの継続時間

Recent simulations of BZ processes

Nagataki (2010,2011)✓ 2D GRMHD (fixed metric)
a=0

a=0.5

Faster is Better ！

a=0.9

a=0.95

まとめニュートリノ重力崩壊

磁場コアバウンス
爆発

中性子星形成恒星進化論の最終段階を統一的に理解すること
まとめ数値シミュレーションで連続的に再現することが不可欠。
今後、高エネルギー数値天文学のグランドチャレンジ
ニュートリノ

磁場

時間発展

ブラックホール形成のダイナミクス
重力崩壊

コアバウンス

極超新星のメカニズム
γ線バーストの形成
中性子星形成爆発
磁場

ニュートリノ加熱

フォールバック

降着円盤、ブラックホール形成

極超新星→γ線バースト

恒星進化論の最終段階を統一的に理解することまとめ
数値シミュレーションで連続的に再現することが不可欠。
今後、高エネルギー数値天文学のグランドチャレンジ
ニュートリノ

磁場

時間発展

ブラックホール形成のダイナミクス
重力崩壊

コアバウンス

極超新星のメカニズム
γ線バーストの形成
中性子星形成爆発
磁場✓「ニュートリノ輸送」+「一般相対論」シミュレーションが必要
ニュートリノ加熱✓ 現在、それぞれのフェイズに特化した計算が進んでいる✓ 今後よりコンシステントな計算へ移行しつつある✓ フォールバック
少なくともExa-scaleの計算資源が必要！
降着円盤、ブラックホール形成 極超新星→γ線バースト