学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

13-NA25

天体活動現象の輻射磁気流体シミュレーション
松元亮治（千葉大学大学院理学研究科）
輻射と磁気流体の相互作用をシミュレートする輻射磁気流体コードを実装し，様々な天体現象に適用している。輻射輸送方程式の光線方向に関する 1 次モーメント式を解く
M1 法に基づく相対論的輻射磁気流体コードを用いてブラックホールへの質量供給率が球対称降着の場合の上限光度（エディントン光度）に対応する臨界降着率よりも遥かに高い場合のシミュレーションを実施し，細く絞られたジェットが噴出すること，輻射が円盤の回転軸方向に強くビーム化され，回転軸方向から観測した場合の光度がエディン
トン光度を超えることを示した。また，輻射磁気流体方程式の陰的解法の高速化を試み，
差動回転円盤の局所計算に用いられるシアリング境界条件を改良することによって15％程度高速化することができた。しかしながら，エネルギー輸送の時間スケールが力学的時間スケールに比べて極端に短くない場合の超並列計算では陰的解法よりも Super
Time Stepping(STS)を用いた陽的解法が適していることがわかった。一般相対論的輻射
磁気流体(GRRMHD)コード，輻射輸送方程式の直接解法コードの開発も進展した。

を，ＦνとＥνの勾配を関係づけた表式を用いて近1. 研究の目的と意義
本研究の目的は，輻射と磁気流体の相互作用を似的に解く Flux Limited Diffusion (FLD) 法をシミュレートする輻射磁気流体コードを超並列計用いたシミュレーションが実施されてきた。しか算機に実装して宇宙現象の大規模シミュレーショ

しながら，この手法では光学的に厚い領域の背後ンに適用することである。このため，物質によるにできる影を扱うことができないという問題があ輻射の吸収と散乱を考慮した輻射輸送方程式

る。そこで我々は昨年度までの共同研究において
輻射輸送方程式の 0 次モーメント式に加えて 1 次

1 ∂


+ n ⋅ ∇  Iν (t , r , n) = −σ ν Iν (t , r , n)
 c ∂tのモーメント式+ σ ν Sν (t , r , n) + σ ν ,s ∫ g (n, n' ) Iν (t , r , n' )dn'∂F ν
+ c 2 ∇ ⋅ Pν = −c (σ ν + σ ν ,s ) Fν∂tを磁気流体方程式と連立させて系の時間発展を調を解いて輻射流束の時間発展を同時に計算する 1

べる。ここで，Iν(t,r,n)は単位時間に単位面積を

次モーメント法（M1 法）のコードを開発してき横切って単位ベクトル n の方向に伝わる振動数νた。ここでＰνは次式で定義される輻射ストレステの輻射のエネルギー（輻射強度），σνは吸収断面

ンソルである。

1Pν (t , r ) = ∫ nn Iν (t , r , n)dnc

積，σν,s は散乱断面積，Sνは輻射の源泉関数（Source function）である。
輻射輸送方程式は空間 3 次元に振動数と光線方

ただし，このままでは方程式が閉じないため，以向を加えた 6 次元空間における輸送方程式である。 下の関係（M1 クロージャー）が成立すると仮定す
この方程式を光線方向について積分し，次式で定る。

義される輻射エネルギー密度Ｅνと輻射流束ベク
トルＦν

1Iν (t , r , n)dnc∫Fν (t , r ) = ∫ n Iν (t , r , n)dn

Eν (t , r ) =  1 − χ  3χ − 1 Pν = 
I+
nn  Eν
2
 2

2
3+ 4 fFχ=
f = ν2
5+2 4−3fcEν

本共同研究では，１）1 次モーメント法に基づを用いて記述される 0 次モーメント式

く輻射磁気流体コードの計算精度，数値的安定性∂Eν+ ∇ ⋅ Fν = σ ν (4π Sν − cEν )∂t

を高めて天体現象のシミュレーションに適用する1

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

こと，２）輻射輸送方程式を陰的に解く際にあらわれる大規模な連立方程式のロバストで並列効率の高い解法を実装すること，３）吸収・散乱過程を記述する輻射輸送方程式を光線方向について積分せずに角度依存性を残したまま直接解く解法を実装すること，を目指している。
輻射磁気流体コードを適用する具体的な対象としては，ブラックホールや中性子星などの天体の図２：ブラックホール降着円盤の熱平衡曲線。横軸は表面密度，縦軸は降着率。
．
ME はエディントン降着率。

まわりに形成される降着円盤，惑星形成の母体となる原始星円盤，太陽・恒星大気，天体からの輻射に照射される星間ガス等がある。特に，ブラッを示した（Kawashima et al. 2013）。しかしながクホール降着流については近年，球対称降着流のらこの計算では増光までの時間を決める角運動量

上限光度であるエディントン光度（輻射圧=重力と輸送率αをパラメータとして与えていたため，いなる光度）程度，あるいはそれ以上の光度で輝く

つ再増光するかを予測することができていなかっ高光度天体が次々に発見されて注目されている。

た。図３の計算例

例えば，2011 年 3 月に発見された Swift J1644

ではα=0.1 とした+57 という天体では発見時の光度は 1048erg/s をが，α=0.01 の場合

超え，電波放射も観測されたことから相対論的ジには再増光までの時間は約 500 日にェットが噴出したらしい。この現象は遠方の銀河中心に存在する 106なる。現象論的パ太陽質量程度のブラ

ラメータを導入すックホールに恒星がることなく角運動

接近して潮汐破壊さ

量輸送過程を組みれ，エディントン降込むことができる

着率を超える降着率

輻射磁気流体シミでブラックホールにュレーションを実星の残骸が落下し続けている事象と考えられる（図１）。

図１：共同研究メンバーの川島朋尚による潮汐破壊シミュレーション
例。濃淡は密度分布。

施することが急務となっている。

図３：エディントン降着率の 100 倍の質量供給がある
場合の輻射流体計算結果。
上 段 は光 度， 下段 は降 着率。亜臨界降着状態(ii)に遷
移後，再増光する。

原始惑星系円盤のように輻射の時間スケールがSwift J1644+57 からは増光後 1 年を経ても 106

力学的タイムスケールに比べて短い系では輻射輸

太陽質量のブラックホールのエディントン光度に送方程式を陽的に解く場合の時間刻みが小さくな対応する光度の X 線放射が観測されていたがりすぎるため，陰的な解法等を適用する必要があ2012 年 8 月にＸ線光度が 1 桁以上減少し，その後る。収束が速く，並列効率の高いソルバーを実装

も暗い状態が続いている。これはブラックホール

することが課題として残されている。

への降着率がエディントン降着率以上の超臨界降

各光線方向の輻射輸送方程式を直接解く解法に着状態（図２の(i)）から，亜臨界降着状態(図２のついて，昨年度までの共同研究では散乱が無視で(ii))に遷移したためと解釈できる。

きる場合について短い特性曲線法に基づくコード

我々は，FLD 法に基づく 2 次元軸対称輻射流体を実装して太陽大気のシミュレーションに適用しコードを用いて，この状態遷移をシミュレートし，てきた。散乱を考慮した輻射輸送方程式の直接解

Swift J1644+57 が再増光する可能性があること

法コードを開発することが課題になっている。
2

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

2. 当拠点公募型共同研究として実施した意義(3) 当公募型共同研究ならではという事項など(1) 共同研究を実施した拠点名および役割分担

従来の輻射磁気流体シミュレーションは流束制

共同研究実施大学：東京大学

限拡散近似（FLD 近似）を用いて実施されること共同研究者と役割分担は以下の通りである。

が多かった。この近似を一歩進め，輻射強度の光

千葉大学・大学院理学研究科

線方向に関する 1 次モーメントである輻射流束の松元亮治（研究代表者）：全体統括、

時間発展を計算する M1 法を適用している点に本

朝比奈雄太：コード実装、ジェットへの適用

共同研究の特色がある。
M1 法に基づくコードは，

工藤祐己：熱伝導モジュールの並列性能向上

輻射流体・輻射磁気流体シミュレーション分野で千葉大学・先進科学センター

世界をリードする複数のグループでほぼ同時に開花輪知幸：M1 コード実装、原始星円盤に適用

発が進められ，天体現象に適用した結果が公表さ海洋研究開発機構・地球内部ダイナミクス領域れはじめたところである。非常に競争が激しい分

廣瀬重信（副代表者）：陰的差分を用いた輻射野であるが，本共同研究を通して M1 法に基づく

磁気流体コードの並列性能向上

相対論的な輻射磁気流体コードが世界に先駆けて

東京大学・大学院理学系研究科

開発され，超並列計算機に実装されてプロダクシ

横山央明：シミュレーションエンジン実装

ョンランが始まった点に本研究最大の意義がある。

堀田英之：ダイナモ計算への適用

輻射のタイムスケールが力学的タイムスケール

飯島陽久：輻射輸送方程式解法の実装に比べて短い場合に必要になる陰的差分は天体シ

東京大学・情報基盤センター

ミュレーションの様々な分野で超並列化のネック

中島研吾，片桐孝洋（副代表者）：になっている部分である。関連する問題としては

陰的差分を用いた輻射磁気流体コードの連立

自己重力系におけるポアソンソルバー，磁場方向一次方程式解法の開発、並列性能向上に依存する非等方熱伝導等があり，問題に応じて国立天文台・理論研究部

前処理の方法，並列化方法等を最適化していく必

大須賀健：輻射輸送方程式直接解法の実装

要がある。大規模連立方程式についての豊富な知高橋博之：相対論的輻射磁気流体コード実装

識と経験があり，この分野のエキスパートである

宇宙航空研究開発機構

東京大学情報基盤センターの中島氏，片桐氏にメ

小田寛：状態遷移シミュレーションの実施

ンバーに加わっていただき，助言を得ている。

法政大学・人間環境学部

3. 研究成果の詳細と当初計画の達成状況

松本倫明：解適合格子法に基づく自己重力磁(1) 研究成果の詳細について

気流体コードの並列性能向上

3.1 M1 法に基づく特殊相対論的輻射磁気流体

九州大学・大学院理学研究院（SRRMHD）コードの適用

町田真美：降着円盤・銀河ガス円盤の大局的

我々は相対論的輻射磁気流体コードの開発を進シミュレーション実施

めてきた。このコードでは輻射輸送方程式の光線

九州産業大学・工学部

方向に関する 0 次，1 次のモーメント式

中村賢仁：熱伝導モジュールの並列性能向上
上海天文台


u j Fr j ∂ t E r + ∂ j Fr j = ρκ  4πγ B − cγ E r +


c


ij
j
2
2
 γ u E r γ ui u j Pru  u j Fr 
− ρκ s −−  γ 2 + 2 cc  c  c

川島朋尚：超臨界降着流輻射スペクトル計算(2) 共同研究分野
超大規模数値計算系応用分野
3

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

 4π u i B γ Fri u k Pri k 1iij∂ t Fr + ∂ j Pr = ρκ 
−
+
2
c
c 
c2
 c
γ F i γ 2 u i
u P ikEr − k r− ρκ s  r −
c
c
 c
+

j
jk
u i  2γ u j Fr − u j u k Pr 

c 
c2


を相対論的磁気流体方程式，マクスウェル方程式と連立させて解く。ここで，κ，κs は流体とともに動く共動座標系で定義した吸収，散乱の不透明度，γはローレンツ因子，B は黒体放射強度である。これらの方程式に加えて M1 クロージャー関係

Prij 1 − χ i j 3χ − 1 i j 
n n Er
=
δ +
2
 2


磁気散逸は Maxwell 方程式に有限の電気伝導度

図４：相対論的輻射磁気流体コードを用いた超臨界降着流のシミュレーション結果。
カラーは輻射流束，矢印は速度場を表す。

を仮定することでその効果を取り入れることがで圧よりも優勢になり，円盤は輻射圧によって支えきる。これらをまとめると流体 5 本，Maxwell 方られる。この輻射圧によってプラズマが加速され

程式 6 本，輻射 4 本の計 15 本の双曲型方程式になて軸方向に絞られた強いジェットと，より広い開り，これらを陽的に解く。本解法の詳細は口角を持ったアウトフローが形成された。ジェッTakahashi and Ohsuga (2013)に報告されている。

トとアウトフローの速度は円盤からの輻射による

このコードをブラックホールへの超臨界降着流

加速と周りの輻射場との衝突による減速の兼ね合に適用した結果を図４に示す。この計算では軸対いによって決まり，その速度は光速の 30-40%にな称性を仮定している。カラーは輻射流束，矢印はることが示された。図４に示されているように降速度場を示す。ブラックホールは 10 太陽質量を仮着流からの輻射は回転軸方向に強くビーム化され，定した。ブラックホールへのガス降着率はエディ

この方向から観測した場合の光度はエディントン

ントン降着率の 100 倍程度である。同様な計算は光度を超える。

を仮定することによって方程式は閉じる。

FLD 法 に 基 づ く 輻 射 磁気 流 体 コ ー ド を 用 い て

今回の計算ではジェットは輻射によって加速さ

Ohsuga ら（2009, 2011）によっても実施されていれ，磁場によるローレンツ力によって細く絞られるが，非相対論的な計算であるため放射抵抗が無ていることがわかった。つまり，遠心力によって視されており，また輻射によるジェットの収束機

広がろうとする力がローレンツ力と釣り合うこと

構が働かないという欠点があった。相対論的なコで細いジェット形状を維持している。理論モデル

ードを適用した本計算により，これらの効果を考によると，ジェットは円盤からの放射による放射慮することが可能になった。

抵抗によって収束すると考えられていたがこの効

超臨界降着では多量のガスが降着するために非果はマイナーであり，第０近似でジェットはロー

常に大きな重力エネルギーが解放される。また，レンツ力によって収束する，といえる。

円盤内のガス密度が高いために円盤内における放射効率が高い。このため，輻射圧がガス圧や磁気4学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

3.3 輻射輸送方程式の直接解法の実装

3.2 一般相対論的輻射磁気流体(GRRMHD)コードのこれまで用いてきた SRRMHD コードや現在開発

開発
ブラックホール近傍では一般相対論の効果が重している GRRMHD コードでは輻射場を光線方向に要になる。特にブラックホールが回転しているとついてモーメント展開して 1 次で打ち切ることでそのエネルギーによって強いアウトフローが形成

式を簡略化した。この方法は輻射輸送方程式を運されることが知られている。そこで本年度後半に動量空間で積分することによって計算コストが軽

GRRMHD コードの構築を行った。メトリックは時間減されるため，より長時間の計算が可能になると的に一定とし，以下のような Kerr-Schild 座標をいう利点がある。しかし，運動量空間を積分して用いた。

しまうために光線どうしの衝突といった非物理的な現象を引き起こしてしまうことが知られている。

ここでαは lapse 関数で重力ポテンシャルを表し，

そこで我々は光線方向への依存性は残し，振動数ι

β は shift と呼ばれるブラックホールの回転によについては積分した輻射輸送方程式を解く数値ス

る時空の引きずりを表す。解くべき式は

キームを開発している。この方法はエネルギー空間を積分しているため，輻射スペクトルは得られないが上述の光線の衝突といった問題を回避できで，それぞれ磁気流体及び輻射のエネルギー運動

る。開発中のコードによる計算例を図６に示す。

量保存，質量保存，誘導方程式である。輻射場はこれまで同様モーメント量で記述し，M-1 クロージャーを仮定すると輻射場のエネルギー運動量テ
ンソルが以下のように与えられる。

ここで

は輻射圧，

は輻射フラックスの強度の指標となる輻射４元速度，

はメトリックの反変成分である。図５にテスト計算結果を示す。

図６：左の境界面から散乱体（破線）に向けて光を照射した場合の計算結果。
（上）輻射輸送方
程式を直接解いた結果。
（下）M1 クロージャー
による結果。直接解法では光は散乱体の間を通り，一部は反射されて散乱体の前面で輻射強度が増える。M1 スキームでは入射光と反射光が衝突し，衝撃波の様な非物理的構造が作られる。
3.4 陰解法にあらわれるマトリクスソルバーの高速化
降着円盤の局所 3 次元計算に用いる FLD 法に基づく輻射磁気流体コードについて，本年度前半はまず，昨年度まで開発をしてきたマルチグリッド
図５：GRRMHD による超臨界降着円盤の大局的数値実験結果。左は輻射エネルギー密度，右はガス密度を示す。

ソルバー（MPI を用いて領域分割による並列化が施してある）について，ボトルネックとなってい5

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

る部分の抽出とその高速化について検討を行った。 クであるため判定が難しいが，少なくとも目立っ
その結果，まず，ボトルネックとなっている部分

た速度向上は見られなかった。

は境界条件のモジュールで使用されている

本年度後半には，上に述べた改良点を施したコ

mpi_allgather であることが判明した。

ードを使用して陰解法におけるクーラン数がより

本コードは，降着円盤の一部を，デカルト座標大きい問題（すなわち，輻射の時間スケールが，系(x,y,z)を用いて取り出すシアリングボックス

流体の時間スケールに比べてより小さい問題）を近似用のコードであり，そのため，動径方向であ解くことを試みた。しかし，この場合，マトリクる x 方向の境界条件は y 方向にシアした周期境界

スソルバーが収束せずに計算が破綻することがとなっている。このシアの程度は，毎ステップ変

度々見られた。主な原因は，線形方程式の解とし化するため，y 方向に領域分割した場合に，必要て得られる輻射エネルギーの値から（非線形方程な情報をどのプロセスと通信するかを判定するの式を通して）求めたガスエネルギーの値が負となは煩雑である。そのため，これまでは x 方向境界

ることであった。そこで，このような場合には，のゴーストセルの値を，
mpi_allgather を使って，ニュートン・ラプソンの反復中に限ってガスエネ

y 方向のすべてのプロセスで共有する方法を取っルギーに最低値を設定する，という対策を試みたてきた。この方法は冗長ではあるもの，上記の煩ところ，ほとんどの場合で，計算を安定に進める

雑さを避けることができる。しかし，前述のとおことが出来るようになった。一方，この方法でもり，診断の結果，この方法が問題となっているこ

うまくいかない場合は，時間刻みを小さくした上とが判明した。

で流体計算からやり直さざるを得ず，今後の課題

そのため，毎ステップ，シアの程度から y 方向として残った。

のどのプロセスと通信するかを特定した上で，
mpi_sendrecv によって y 方向の情報を共有するよ

3.5 非等方な拡散項を含む方程式解法の高速化うにした。これによって，ボトルネックとなって天体プラズマでは磁力線に沿った方向への熱伝

いた部分に関しては，必要最小限の情報通信で済導の効率が流体速度の 1000 倍を超える場合があ

むようになったが，実際の計算速度では最大 15%

り，このような硬い方程式系を解く効率的な数値

程度の高速化に留まった。

解法が必要となる。光学的に厚い領域の輻射輸送

このように，シアリングボックス近似を用いたは拡散方程式に従うため，熱伝導方程式の解法をシミュレーションの一つの問題は，x 方向の境界高速化することは輻射輸送方程式の解法を高速化

条件の複雑さにある。一方，バックグラウンドのするヒントにもなると期待される。

シア流は時間的に変化しないので，これを全体のまず，陰的解法を検討した。係数行列が非対称

流れ場から差し引くことで数値誤差を減らすことであるため，前処理付き Krylov 部分空間法の一種が出来るほか，x 方向の境界条件の一部を簡単化であり，安定性に優れたリスタート付の Flexible

することができる。そこで，本年度はこの方法の

GMRES 法を実装し，並列化を行った。行列ベクト実装も行った。実は，これまで x 方向のシミュレ

ル積の計算効率化のため，係数行列の保存形式とーションボックスサイズを大きくした際に x 方向

し ての境界近傍において数値的エネルギー散逸率の評

ELLpack-itpack format (ELL)，Compressed Sparse

価に数値誤差が大きく影響していた。しかし，バRow format (CSR)について FX10 上で速度の比較を

ックグラウンドシア流を差し引くことにより，こ行った結果に基づき，有限差分的なステンシルにの数値誤差の問題はほぼ完全に解決された。一方，

対して高速な DIA 形式を採用した。

DIAgonal

storageformat(DIA) ，

境界条件が簡単化されたことによって計算速度がGMRES 法の前処理として SSOR(対称 SOR)法及び向上したかどうかについては，流れがカオティッ

ILU(0)分解(LU 分解のうち係数行列の非ゼロ成分6学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月と同じ場所のみを計算したもの)を試みた結果，

る光線衝突による衝撃波の発生等の非物理的な現SSOR 前処理は緩和パラメータを適切にとることで， 象を生じない計算結果が得られることを確認した。
5 倍程度の収束率が改善した。マルチグリッド前

この課題は当初目的をほぼ達成することができた。

処理などのより高度な前処理を行った場合，さら【陰解法に現れるマトリクスソルバーの高速化】に数倍の改善が見込まれる。他方 ILU(0)前処理は原始惑星系円盤等では輻射拡散方程式とガス・

ゼロピボットが頻繁に表れ，分解がスムーズに行輻射エネルギー交換の非線形方程式を連立させた

えないことが分かった。

方程式を陰的に解く必要がある。本年度の目標は，陰的解法とは別のアプローチとして，Super

このマトリックスソルバーについて，これまでにTime Stepping (STS)法を検討した。STS 法では，実装してきたマルチグリッドソルバーの効率化，複数のステップをひとまとめにして安定性を考え高速化をはかることであった。シアリングボック

ることにより，放物型微分方程式の安定条件を緩ス近似を用いたコードでは，ボトルネックとなっ和させることができる。STS 法を非等方熱伝導にているのはソルバー自身というよりも，MPI で領適用した場合，SSOR 前処理付き GMRES 法のさらに域分割した場合の境界条件に伴う通信部分である

10 倍程度の改善が見込まれることが分かった。
STSことが判明した。そこで，この境界条件部分の高法は陽的な解法であるため，大規模な並列化に伴速化に取り組んだ。結果は，最大でも 15%の速度う性能低下の心配も少ない。以上の理由から，非向上ということで，検討の余地が残されている。

等方熱伝導の解法として Super Time Stepping 法

当初計画では，１）これまでに実装してきたマルチが有望であるという結論を得た。

グリッドソルバーの効率化・高速化，２）前処理として不完全 LU 分解， 反復法として BiCGSTAB,(2) 当初計画の達成状況について

GPBiCG 等のクリロフ部分空間法を用いた手法と【1 次モーメント法に基づくシミュレーションコの並列性能の比較，の 2 点を行う予定であった。

ードの改良と高光度降着流への適用】

このうち，１）については，前節に述べたように，

特殊相対論的輻射磁気流体(SRRMHD)コードを用

現時点で可能な対策はほとんど行っており，ほぼいたブラックホールへの超臨界降着流のシミュレ

100%の達成度と言って良い。一方，２）についてーションを実施し，輻射圧で加速された光速の数

は，１）に時間を取られたため，ほとんど進める

10%の速度を持つジェットが形成されることを示ことが出来なかった。

すことができた。ブラックホール降着流の状態遷

熱伝導方程式についてはエネルギー輸送の時間

移等，より長い時間スケールの計算に適用するたスケールが力学的時間スケールよりも短い場合のめには数値的な安定性を改善する必要がある。こ

非等方熱伝導の解法として前処理付き GMRES 法及の課題の達成度は 60%程度である。

び Super Time Stepping(STS)法を検討した。そのブラックホール近傍の計算の数値的安定性を高結果，STS 法が有望であるという結論を得て天体

め，また，ブラックホールの回転の効果等をとりシミュレーションに適用中である。

入れるため，一般相対論的輻射磁気流体(GRRMHD)
コードの開発を行った。コード開発は順調に進み，

4. 今後の展望
数値的安定性に改善の余地はあるものの，高光テスト計算結果が得られつつある。GRRMHD については当初計画以上に進展したと考えている。

度降着円盤に適用可能な国際競争力のあるコード【輻射輸送方程式の直接解法の実装】が実用段階に達したことから，これを用いて超臨振動数積分した輻射輸送方程式の直接解法を実界降着状態におけるジェット・アウトフロー形成，

装してテスト計算を実施し，散乱体での反射等の降着率の減少に伴う超臨界降着状態から亜臨界降計算において M1 スキームによる計算にあらわれ着状態への状態遷移過程等の 3 次元シミュレーシ
7

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

ョンを実施する計画である。

に変更することを試みる。ただし，シアリングボ

ジェット・アウトフローによる質量放出率等はックス近似を用いているコードに関しては，もし一般相対論的効果，特にブラックホールの回転に

前述のように，境界条件がボトルネックである状依存することが予想される。１）降着円盤から噴況が変わらないとすれば，１）境界条件を呼ぶ回出するジェットの速度は観測される相対論的ジェ数が出来るだけ少なくなるような手法が最も高速ットを説明できるか，２）ジェット／アウトフロ化につながる，２）マトリクス要素の境界条件部

ーのパワーは観測を説明できるか，また，そのエ分が毎ステップ変化していることから，マトリク

ネルギーの内訳(ガス，磁場，輻射)はどのようにス要素生成の高速化もポイントとなる，といったなっているのか，３)強いジェット／アウトフロー

ことが予想される。

が形成された場合，ブラックホールへの質量降着

問 題 に よ っ て は 陰 解 法 を 用 い ず ， STS(Super

率はどの程度になるか，特にその値に上限値はあTime Stepping)法に基づく陽解法を適用することるのか，について調べる必要がある。これらにはも試みていく予定である。

一般相対論効果が重要であり，GRRMHD コードを用いて初めて解決する問題である。より数値的に安5. 研究成果リスト

定な GRRMHD コードを完成させ，降着円盤の大局的(1) 学術論文（投稿中のものは「投稿中」と明記）

な数値実験を行いたい。

・Takahashi, H.R. and Ohsuga, K., A Numerical
Treatment of Anisotropic Radiation Fields

他方，M1 スキームに基づく輻射磁気流体コード
は，光線どうしの衝突のために非物理的な衝撃波

CoupledwithRelativistic

Resistive面が形成される等の弱点もある。そこで，輻射輸

Magnetofluids, Astrophysical Journal, 772,

送方程式を光線方向について積分せず，角度依存127 (16pages), 2013

性を残して解く直接解法コードの開発を進める。

・Kawashima, T., Ohsuga, K., Usui, R., Kawai,

今年度実装した，振動数について積分した輻射輸

N., Negoro, H., and Matsumoto, R., Recurrent

送方程式を解くコードの効率を高めて計算を実施

Outbursts and Jet Ejections Expected in Swift

し，M1 コードによる結果と比較して，M1 コードのJ1644+57;

精度，適用限界などを明らかにしたい。近い将来，

Supermassive Black Hole, Publ. Astron. Soc.

「京」計算機等を用いて輻射輸送直接解法に基づJapan 65, L8 (5pages), 2013

Limit-Cycle

Activities

ina・Kanno, Y., Harada, T., and Hanawa, T., Kinetic

く輻射磁気流体シミュレーションを実施すること

Scheme for Solving the M1 Model of Radiativeを視野に入れて研究を進める計画である。

Tarnsfer, Publ. Astron. Soc. Japan 65, 72

陰解法に現れるマトリクスソルバーの高速化に(19pages), 2013

関して， 本年度は，これまでに実装してきたマル
チグリッドソルバーの効率化・高速化を行うとと・Nomura, M., Ohsuga, K., Wada, K., Susa, H.,

もに，計算をロバストに進めるための処置・工夫and Misawa, T., Modeling Line-Driven Diskを行った。しかしながら，陰解法におけるクーラ

Wind for Broad Absorption Lines of Quasars,

ン数が大きい問題ではマトリクスソルバーが収束

Publ. Astron. Soc. Japan 65, 40 (10pages),

せずに計算が破綻する場合が見られた。これは，2013

ガウス―ザイデル法をスムーザとして用いている

・Takeuchi, S., Ohsuga, K., and Mineshige, S.,

現在のマルチグリッドソルバーの限界を示していClumpy Outflows from Supercritical Accretion

る。そこで，今後は，今年度取り組むことが出来

Flow, Publ. Astron. Soc. Japan 65, 88

なかった，マトリクスソルバー自身をまったく別(12pages), 2013
・Yang, X., Yuan, F., Ohsuga, K., and D. Bu, Twoの方法（すなわち，前述のクリロフ部分空間法）
8

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

Dimensional

Numerical

Simulations

T., Magnetohydrodynamic Simulations of theofSupercritical Accretion Flows Revisited,

Hard-to-Soft

Astrophys. J. 780, 79 (15pages), 2014

Accretion Flows, APPC12, The 12th Asia

・Turner, N.J., Benisty, M., Dullemond, C.P, andPacific

Hirose, S., Herbig Stars' Near-Infrared

Transition

PhysicsinBlack

Conference

of

Hole
AAPPS,

Makuhari, Chiba Japan, July 14-19, 2013

Excess: An Origin in the Protostellar Disk's

・Kawashima, T., Ohsuga, K., Usui, R., Kawai,

Magnetically-Supported Atmosphere,

N., Negoro, H., Matsumoto, R., Limit-Cycle

Astrophys. J. 780, 42 (9pages), 2014

Activities in a Tidal Disruption Event by aSupermassive Black Hole, APPC12, The 12th

(2) 国際会議プロシーディングス

Asia Pacific Physics Conference of AAPPS,・ Hanawa, T., Kanno, Y., and Harada, T.,

Makuhari, Chiba, Japan, July 14-19, 2013

Reconstruction Method for M1 Equations of・Takahashi, H.R., Numerical Study ofRadiative Transfer, Proceedings of a 7th

Supercritical Accretion Flow onto the Black

International

Hole Using Relativistic Radiation

Conferenceon

Numerical

Modeling of Space Plasma Flows (ASTRONUM

Magnetohydrodynamic code, APPC12, The 12th2012), Ed. Pogorelov, N.V., Audit, E., andAsia Pacific Physics Conference of AAPPS,

Zank, G.P., Astronomical Society of the

Makuhari, Chiba, Japan, July 14-19, 2013
・ Takahashi, H.R,

Pacific, p.233, 2013

Explicit-Implicit Scheme

for Relativistic Radiation Hydrodynamics,(3) 国際会議発表

ASTRONUM 2013, 8th International Conference・Matsumoto, R., Kawashima, T., Machida, M.,

on Numerical Modeling of Space Plasma Flows,Global Three-dimensional Magneto-

Biarritz, France, July 1st - July 5th 2013,hydrodynamic Simulations of Hard-to-Soft(invited)Transition and Jet Formation, Suzaku-MAXI・Hanawa, T., Radiative Hydrodynamical

2014, Expanding the Frontiers of the X-ray

Simulations of Protoplanetary Disks,

Universe, Ehime, Japan, February 19-22, 2014

ASTRONUM 2013, 8th International Conference・Takahashi, H.R., Numerical Study of Jets fromon Numerical Modeling of Space Plasma Flows,

Super Critical Accretion Disks, Suzaku-MAXI

Biarritz, France, July 1st – July 5th, 20132014, Expanding the Frontiers of the X-ray

・Matsumoto, R., X-ray Flares Triggered by the

Universe, Ehime, Japan, February 19-22, 2014

Tidal Disruption of a Gas Cloud Approaching・Hirose, S. et al., Thermal Equilibria of

the Galactic Center Supermassive Black Hole,

Dwarf-Nova Accretion Disks obtained by

12th International Workshop on theradiation MHD simulations, Kyoto

Interrelationship between Plasma

Mini-Workshop on Dwarf Novae and RelatedExperiments in the Laboratory and in Space

Systems, Kyoto University, October 29, 2013

(IPELS), July 1-5, 2013, Hakuba, Nagano,・ Ohsuga, K., Radiation-MHD Simulations ofJapan (invited)

Accretion Flows and Outflows around Black・Hirose, S. et al., MHD Turbulence in Accretion

Holes, Black Holes, Jets, and Outflows,

Disks with Ionization Transitions, Revealing

Katmandu, Nepal, Oct. 14-18, 2013 (invited)

Radiative

・Matsumoto, R., Ono, T., Oda, H., and Kawashima,

Processes

near

Black

Holes,

Princeton Center for Theoretical Science,9

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

年 9 月 10 日－12 日

May 3, 2013 （invited）
・Ohsuga, K., Global Simulations of・川島朋尚，小川崇之，松元亮治，銀河系中心 SgrSuper-Eddington Accretion Flows and Outflows,

A*に接近するガス雲 G2 と高温降着流の相互作

Revealing Radiative Processes Near Black

用：3 次元磁気流体シミュレーションによるフHoles, Princeton Center for Theoretical

レア予測，日本天文学会秋季年会，東北大学，

Science, May 1-May 3, 2013 (invited)2013 年 9 月 10 日
・高橋博之，大須賀健，超臨界降着流からのジェ(4) 国内会議発表

ット・アウトフローの形成，日本天文学会秋季・松元亮治，川島朋尚，Swift J1644+57 が示す巨

年会，東北大学，2013 年 9 月 12 日

大ブラックホール降着円盤の状態遷移，日本天・廣瀬重信，MHD 乱流のグリッドスケールでのエ文学会 2014 年春季年会，国際基督教大学，2014

ネルギー散逸率の評価について, 日本天文学会

年 3 月 22 日

秋季年会, 東北大学, 2013 年 9 月 10 日

・川島朋尚，松元亮治，ガス雲 G2 の衝突による銀河中心 Sgr A* 高温降着流内部での磁場増幅と(5) その他（特許，プレス発表，著書等）

降着率変化，日本天文学会 2014 年春季年会，国なし

際基督教大学，2014 年 3 月 21 日
・川島朋尚，小川崇之，小田寛，松元亮治，ブラックホール降着流状態遷移中の準周期的な磁気リコネクション，日本流体力学会年会 2013，東京農工大学，2013 年 9 月 13 日
・松元亮治，小川崇之，小野貴史，川島朋尚，銀河系中心ブラックホールへのガス雲落下の磁気
流体数値実験，日本流体力学会年会 2013，東京
農工大学，2013 年 9 月 13 日
・高橋博之，大須賀健，超臨界降着からのアウトフロー形成，日本流体力学会年会 2013，東京農工大学，2013 年 9 月 13 日
・廣瀬重信，Omer Blaes, Matt Coleman, Julian
Krolik, 佐野孝好，水素電離を伴う降着円盤におけるエネルギー輸送, 日本流体力学会年会,
東京農工大学, 2013 年 9 月 13 日
・松元亮治，小野貴史，小田寛，川島朋尚，ブラックホール候補天体の明るいハードステートにおける時間変動の起源，日本天文学会秋季年会，
東北大学，2013 年 9 月 10 日－12 日
・川島朋尚，大須賀健，薄井竜一，河合誠之，根來均，松元亮治，突発天体 Swift J1644+57 の再バースト予測 II: 潮汐破壊星による質量供給率の時間変化を考慮した輻射流体シミュレーシ
ョン，日本天文学会秋季年会，東北大学，201310