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

jh140043-NA24
相対論的流体コードの開発と超高エネルギーガンマ線連星系への応用
岡崎 敦男（北海学園大学工学部）
概要
我々のグループは、SPH (Smoothed Particle Hydrodynamics) コードを用いて、世界に先駆けて超高エネルギーガンマ線を出す連星系の 3 次元シミュレーションを行ってきたが、シミュレーションの信頼性と効率性にいくつかの解決すべき課題があった。平成 26 年度には、それらの課題のうち、密度が大きく変化するところで計算精度が悪くなるという問題、およびパルサーと呼ばれる高速回転する中性子星からの相対論的な粒子の流れ（パルサー風）を同じ運動量流速を持つ非相対論的な流れで近似しているために広い領域での流れの形状を再現できないという問題を解決するようにコードを書き換え、それを超高エネルギーガンマ線連星系に適用することを計画し、年度前半の研究で前者の問題についてはコードの改良が終了し、テストシミュレーションも終了した。しかしその時点で、超ガンマ線連星系のシミュレーションにおいては、パルサー周囲のガス圧に応じてパルサー風が抑止される効果を取り入れることが必須であることを見いだしたので、そのようにコード改良しシミュレーションを行うことを最優先させた。結果として、パルサー風抑止効果を取り入れたシミュレーションを世界で初めて行うことができたが、当初の計画にあったコードの相対論化は準備の段階でとどまった。

1.

共同研究に関する情報(1) 共同研究を実施した拠点名
北海道大学情報基盤センター(2) 共同研究分野

氏 名（所属）役割

岡崎 敦男（北海学園大研究の統括とコード学工学部）の改良・数値計算

大宮 学（北海道大学情コードの並列化とチ報基盤センター）ューニング

河内 明子（東海大学理

数値計算

■

超大規模数値計算系応用分野学部）□

超大規模データ処理系応用分野

早崎□

超大容量ネットワーク技術分野

Astronomy and Space□

超大規模情報システム関連研究分野

Science Institute)(3) 参加研究者の役割分担
本共同研究の研究体制は以下の通り。

公威

コードの改良と数値計算

Christopher M. P.

コードの改良と数値Russell （デラウェア大計算学、アメリカ）

1(Korea2.

研究の目的と意義

2.1.

研究の背景

1. これまで用いてきた SPH コードを、密度の大きく変化するところでも精度良く計算でき、（特殊）相対論的な流体を扱え、KH 不安定性

近年、世界各地に撮像型チェレンコフ光望遠鏡を再現できるように書き換えるとともに、が建設され、ガンマ線天文学は本格的なイメージ

2. 広範囲の MPI 並列化とチューニングを行い、ング観測時代へと入った。それに伴い、多くの天体

3. 改良した相対論的コードを用いて超高エネルが超高エネルギーガンマ線（E > 100 GeV）を放射ギーガンマ線連星系の大規模数値シミュレーしていることが明らかになってきた。特に、最近でションを行い、観測量との詳細な比較を行う

は超高エネルギーガンマ線を放射する天体が、高エネルギー粒子加速の現場として注目されている。

ことである。そのためのコードの開発・改良は膨大な作業量になるので、研究に与える影響の大きさにそれらの天体の中に 超高エネルギーガンマ線連星より優先順位を付け、平成 26 年度は、コードを密系 と呼ばれる新しいクラスの天体がある。これら

度の大きく変化するところでも精度良く計算できの天体は大質量星と超強重力天体（中性子星あるるように改良（いわゆる grad-h 補正項を導入）しいはブラックホール）からなる連星系だと考えられた後、（特殊）相対論化に特化して作業を行い、そるが、ほとんどの系で超強重力天体の正体が分かっ

れを用いて超高エネルギーガンマ線連星系の大規ておらず、それらの系では相互作用や高エネルギー模数値シミュレーションを行うことを目標とする計放射機構について未解明な状態が続いている。

画を立てた 1 。

我々は 数 年 前 か ら 、3 次 元 SPH (SmoothedParticle Hydrodynamics) コードを用いて、いくつ2.3.

かの超高エネルギーガンマ線連星系に対するシミュ研究の意義
超高エネルギーガンマ線連星系における相互レーションを行い、相互作用の詳細を明らかにして作用、およびそこからの高エネルギー放射の物理きた。これらの計算はこの分野では例の少ない 3 次を明らかにするためには、3 次元動的モデルを構築

元流体シミュレーションであり価値の高いものであし、それを観測結果と比較することが必須である。

るが、結果の信頼性と計算の効率性に関していく

現在、3 次元シミュレーション行っているのは我々つかの解決すべき課題があった。前者ついては、計の他に 1 グループだけであり、シミュレーションで算全般に関わる問題として、(1) 密度が大きく変化得られたデータを用いて高エネルギー放射を計算するところで計算精度が悪くなることと、(2) 衝撃する手法についても我々が一歩先を行っている。し波領域で物質混合が起こりにくい、という問題があかし、前述のように、我々のモデルにはいくつかのり、特殊相対論的シミュレーションに関わる問題と解決すべき課題がある。本研究と来年度以降の研して、(3) パルサー風を同じ運動量流速を持つ非相究によりそれらの課題を解決すれば、モデルと観測対論的な流れで近似しているために広い領域でコリ結果との詳細な比較が可能になる。そのようにしてオリ力の影響を強く受けてしまい、衝撃波面の形状得られた 3 次元モデルは、超高エネルギーガンマが正しい形にならない (Bosch-Ramon et al. 2012)
1という問題があった。後者については、(4) MPI 化

テストシミュレーションで効果を確認した後、超ガンマ線が小規模で並列化の効率が悪いことが、より大規模なシミュレーションを行うことを困難にしていた。

2.2.

実際の作業では、後述のように、grad-h 補正項を導入し、連星系シミュレーションにおいてより緊急の対応を必要とする改良点が明らかになったため、それへの対応とそれに基づくシミュレーションを優先して行った。その結果、研究の目的

コードの相対論化作業はどこをどのように書き換えるか方

本研究課題の最終的な目的は、針を立てるところまでしか進まなかった。

2

線連星系の研究を大きく前進させるはずである。

代表的なやり方が変分原理から基礎方程式を導く

また、本研究で得られる知識や開発する手法は十分方法で、この方法を用いると自然に h の勾配を含に一般的なものなので、他の大質量連星系やガス円む項（grad-h 補正項と呼ばれる）が基礎方程式に盤を持つ様々な天体に対しても容易に適用できる。 含まれることになる（例えば、Rosswog 2009）。
さらに、我々の開発してきた計算コードは他の研

我々のコードは grad-h 項を含んでいない旧来の究グループ［デラウェア大学（アメリカ）、NASA、 定式化によっていたので、コードの変更のための具サウサンプトン大学（イギリス）、サンパウロ大学

体的な作業手順を以下のように 2 段階に分けて行（チリ）］にも利用されている。本研究で行う計算うことを計画した。

コードの開発・改良は、それらのグループの研究に第 1 段階: 非相対論的な変分原理に基づく定式化にも役立つものである。

基づき、SPH コードを改良する。主な改良点は grad-h 項を含む形への書き換えである。書

当拠点公募型共同研究として実施した意義

3.

き換えが適切に行われたことを確認するため
本研究の基礎をなすのは、SPH 法に基づく 3 次テストシミュレーションを実施する。

元流体コードである。高精度なシミュレーション

第 2 段階: 第 1 段階で用いた基礎方程式を相対論的を行うためには 100 万個以上の粒子が必要となり、なものに変更する。SPH コードの変更が終了

そのような計算を非常に小さなタイムステップでしたら、変更が適切に行われたことを確認す

何周もの連星軌道周期にわたって実施するには非るためテストシミュレーションを実施する。

常に大きな計算機資源を必要とする。そのために、JHPCN の計算機資源を利用したく、ここに申請し平成 26 年度前半（中間報告までの期間）で、第 1段階については、コードの変更、テストシミュレーた次第である。

ションとも終了した 2 。テストシミュレーションは、• Be 星星周円盤の形成と伴星との潮汐相互作用

前年度までに得られた研究成果の概要

4.

（衝撃波なし）新規課題のため、該当せず。

• 大質量連星系における恒星風の衝突（強い衝撃波。2 つの恒星風粒子の質量に大きな違い

今年度の研究成果の詳細

5.

なし）
以下に、中間報告までに得られた成果 (§ 5.1) と• 超高エネルギーガンマ線連星におけるパルその後に得られた成果 (§ 5.2) について詳述する。

サー風、星周円盤、恒星風の衝突（強い衝撃5.1.

波。3 種類の粒子に非常に大きな質量の違いGrad-h 補正項の導入とテストシミュレーシがある）

ョンという 3 つの異なるタイプの系に対して行った。

旧 来 の SPH 法 に よ る 定 式 化 で は smoothing

以下に、テストシミュレーションの結果についてlength h（粒子の大きさに対応する量）の場所に述べる。

よる違いを無視していた（h の勾配をゼロとおいていた）ため、h が大きく変化するところでは計算が不正確になっていた（例えば、Springel & Hernquist

2002）。そのため、現代的な SPH コードは h の勾2配を考慮した定式化に基づいて書かれている。その後述するように、相対論化作業を行わなかったため第 2 段階は未実施。

3

図 1. Be 星とブラックホールの連星 AGL J2241+4454

図 2. Be 星とブラックホールの連星 AGL J2241+4454のシミュレーションにおける動径方向の物理量分布：(a)のシミュレーションにおけるブラックホールへの降着grad-h 補正項を無視した計算 vs. (b) grad-h 補正項を率の位相依存性：(a) grad-h 補正項を無視した計算 vs.

含んだ計算。左のパネルは面密度の変化（4 周期毎）、(b) grad-h 補正項を含んだ計算。青実線は降着率、青破

右パネルは 16 周期後の物理量。黒線が面密度、赤線が線は平均降着率、赤実線は星周円盤の質量。位相 0 と 1

回転速度、青線が動径速度。

は近星点5.1.1.

Be 星星周円盤の形成と伴星との潮汐相互 （青線）が示されている。2 つのシミュレーションの比較から、密度が緩やかに変化している領域では作用のシミュレーション
連 星 系 の パ ラ メ ー タ と し て は 、最 近 初 め て

両者の間には差がほとんど見られないが、密度が急発 見 さ れ た Be 星 と ブ ラック ホ ー ル の 連 星 系

激に変化する領域では差が現れることが見て取れ

AGL J2241+4454 のものを使用した。Be 星とブる（右のパネルの青線で表された動径速度分布）。

ラックホールの質量はそれぞれ太陽の 13 倍と 5.3

図 2 は 2 つのシミュレーションにおけるブラック倍で、軌道周期は 60.37 日、軌道離心率 0.1 である。 ホールへの降着率を比較している。図 2(a) は grad-grad-h 補正項を含まない場合のシミュレーショ h 補正項を含まない場合、図 2(b) は grad-h 補正項ンと含む場合のシミュレーションの比較を図 1–2を含む場合である。図から、両者の結果は全体としに示す。図 1 は動径方向の（= 方位角方向に平均て似ているが、grad-h 補正項を含む場合の方がノし、z 方向に積分した）物理量の分布を示していイズが少なくなっていることが読み取れる。これもる。図 1(a) は grad-h 項を含まない場合、図 1(b)また、grad-h 補正項を考慮することにより期待さは grad-h 項を含む場合である。左のパネルでは星れる効果である。

周円盤形成の様子が 4 周期毎の面密度分布として5.1.2.

衝突恒星風のシミュレーション表されており、右のパネルでは形成開始から 16 周大質量星どうしの連星系ではそれぞれの恒星か期後の面密度（黒線）、回転速度（赤線）、動径速度ら高速で放出された物質（恒星風）が衝突し、そこ4から強い X 線が出ている。エータ・カリーナはそgrad-h の効果も無視しているので、たまたま問題のような相互作用を示す代表的な連星系（周期 5.5が起こっていなかったようだ。この問題の解決策と年、軌道離心率 0.9）である。この系に対して、2しては、短期的には粒子数を大幅に増加すると共につの恒星風の粒子の質量を同じにして、近星点近密度の低い媒質の解像度を下げることにより、ほぼ傍の軌道位相でシミュレーションを行ってみたが、 等しい質量を持つ 3 媒質で計算することが考えら

grad-h 補正項を含まない場合と含む場合で明確なれるが、長期的には密度によらない定式化 (Saitoh

差は出なかった。

& Makino 2013) を用いてコードを書き換える方が5.1.3.

良いだろう。新たに判明した大きな課題である。

超高エネルギーガンマ線連星におけるパルサー風、星周円盤、恒星風の衝突のシミュ

5.2.

レーション

研究計画の変更：パルサー風を ON/OFF するスイッチの導入と超高エネルギーガンマ線

図 3 は 、超 高 エ ネ ル ギ ー ガ ン マ 線 連 星 PSR

連星シミュレーション

B1259-63（Be 星と中性子星の連星系。周期 3.4 年、離心率 0.87）のシミュレ−ション開始からわずか 1

5.2.1.

計画変更の背景

我々は平成 25 年度から HESS J0632+057 と呼日後のスナップショットである。grad-h 補正項を含ばれる超高エネルギーガンマ線連星系のシミュレー

むシミュレーション（下図）は、grad-h 補正項を含ションに取り組んできた。この系では、Be 星と呼まないシミュレーションに比べてパルサー風の動圧ばれる星周円盤と恒星風を持つ大質量星と正体不がかなり大きくなっており、Be 星星周円盤と恒星明の強重力天体が、軌道周期が約 315 日、軌道離心風を非常に素早く吹き払っている。grad-h 補正項

率が 0.83 という非常に細長い軌道上を運動していを含んだシミュレーションでのパルサー風の動圧はる。我々は、HESS J0632+057 では Be 星星周円盤理論的に予想される値よりもかなり大きくなっておが連星軌道面に対して傾いていて、そのような円り、シミュレーションは図 3(b) の直後にクラッシュ盤とパルサー風の相互作用が遠星点付近で強くなした。それに対して、grad-h 補正項を含まないシるのだろうと考えシミュレーションを始めた。用いミュレーションでは衝撃波の形状は解析解と似たもた粒子数は数十万個程度だが、パルサー風というのになり、計算もクラッシュすることはなかった。

grad-h 補正項を含んだシミュレーションでパル高速の粒子の流れを短いタイムステップ（連星軌道周期の約 10−7 ）で約 3 軌道周期を追いかけるためサー風の動圧が強くなりすぎたのは、密度が極端にに、1 つのシミュレーションに約 3 か月を要し、平異なる 3 種類の媒質（密度が高い方から順に、Be成 26 年度に入ってからようやくシミュレーション

星星周円盤、Be 星恒星風、パルサー風）を極端にが終了した。しかし、シミュレーションを基に高エ異ならない粒子数（極端に異ならない空間解像度）で表すために、粒子に極端に異なる質量（パルサーネルギー放射の光度曲線を計算して見たところ光度の極大が近星点直後に現れ、観測された光度曲線

風粒子の質量は、恒星風粒子の質量の 1/100、星周とまったく合わなかった。

円盤粒子の質量の 1/10000）を与えたためと考えらBe 星を持つ超高エネルギーガンマ線連星系のシれる。SPH 法では、各点での物理量は近傍粒子のミュレーションは、その形状の複雑さのために、我々寄与を足し合わせたものなので、そもそもこのよのグループ以外にはまだ誰もシミュレーション ZAうに極端に異なる質量を持つ粒子を用いると媒質できていないが、HESS J0632+057 では、高エネの境界での物理量が正しく表されない。現状のコー

ルギー放射光度が遠星点付近で極大になるというドでは、物質混合がほとんどおこらないために異不可解な振る舞いを示すため（他の超高エネルギーなる質量を持つ粒子の分布が分離しており、さらに5図 3. 超高エネルギーガンマ線連星 PSR B1259-63 におけるパルサーと Be 星の相互作用の様子を表すスナップショッ
ト：(a) grad-h 補正項を無視した計算 vs. (b) grad-h 補正項を含んだ計算。近星点は x 軸の負方向にあり、連星軌道はx-y 平面内にある。左パネル：x-y 平面上の密度分布、中央パネル：x-z 平面上の密度分布、右パネル：y-z 面上の密度分布。

ガンマ線連星系では近星点付近で極大になる）、シ

ミュレーションから求めた X 線光度を示している。

ミュレーションどころかおおまかなシナリオすら下のパネルで、赤線はパルサー風を抑止する効果出されていなかった。そのような天体のシミュレーを組み込んだシミュレーション、紫色の線はその効ションを成功させることは、この分野の研究にとっ果を含まない場合のシミュレーションである。図かても我々の研究の評価にとってもとても重要なこと

ら、パルサー風が抑止される効果を取り入れたシである。このように考え、コードの相対論化作業をミュレーションは観測の傾向をおおむね再現できて一時棚上げし、この系のシミュレーションを成功さいることが分かる。ちなみに、パルサー周囲のガスせることを最優先することにした。

圧によりパルサー風が抑止される効果を組み込ん

5.2.2.

だシミュレーションは世界で初めての試みである。

コードの改良とシミュレーション結果

その結果は、平成 27 年 5 月に開かれた国際会議で検討の結果、近星点付近ではパルサー周囲のガ発表した。現在、論文を執筆中である。

スの圧力が高く、それによりパルサー風が抑止される可能性が高いことが分かり、その条件（を抑止期

6.

今年度の進捗状況と今後の展望

6.1.

今年度の進捗状況間も含めて定量的に表したもの）をコードに組み込むことにした。若干の試行錯誤の後、コードが正しく動くようになったので、数例のシミュレーショ

本研究の計画書に記載した目的は、(1) 計算コーンを行った。結果の一例を図 4 に示す。上のパネドの特殊相対論化とテスト計算、(2) 相対論的計算

ルは観測された X 線光度曲線で、下のパネルはシコードの並列化とチューニング、(3) 相対論的コー

6

表 1. 計画の進捗状況

項目

計算コードの相対論化

当初計

新規計

画の達

画の達成度

成度

(%)

(%)

40—0——

100

0100

20100とテスト計算
コードの並列化とチューニング
パルサー風抑止効果のコード化とテスト計算
超高エネルギーガンマ線連星への応用
全体とチューニング、および相対論化したコードを用いた大規模計算も手つかずのままである。それらを考慮すると、当初の計画の達成度は表 1 のように20%程度である。
パルサー風が抑止される効果についてのコード
図 4. 超高エネルギーガンマ線連星 HESS J0632+057

改良とシミュレーションについてはどちらも計画をの X 線光度曲線。上図：観測された光度曲線 (Aliu et100%達成した（表 1）。

al. 2014, Fig. 3)。下図：シミュレーションから求めた6.2.

光度曲線。赤線はパルサー風を抑止する効果を組み込んだシミュレーション、紫色の線はその効果を含まない場

今後の展望

前述のように、作業はコードの相対論化の手前で合のシミュレーション。両方の図とも横軸は軌道位相。

止まっている。それを再開する手はあるが、§ 5.1.3

位相 0 と 1 が近星点で、位相 0.5 が遠星点。

で述べたように、パルサー風、恒星風、星周円盤のドの超高エネルギーガンマ線連星系への応用の 3 つ

密度が極端に異なるため、それらを同程度の質量である。さらに、(1) は smoothing length h の勾配の粒子で表すには、非常に大きな粒子数を必要とすを考慮した定式化に基づく改良と相対論化の 2 段る。今のままで、grad-h 補正項を組み込んだ方式階に分かれる。また、計画当初は、(1) に 50%の時を単純に相対論化するだけでは計算時間が長大に間、(2) と (3) にあわせて 50%の時間をかけることなりすぎて実用にならない。そこで、計算効率を大になるだろうと考えていた。

きく上げるために、現行のコードではごく小規模に平成 26 年度前半に (1) の前半段階が終了したが、その後、我々にとってより緊急なコード改良（パルしか行っていない MPI 並列化をコード全体にわたり行うべく、コードを見直す作業を最近始めた。

サー風抑止効果の導入。§ 5.2 参照）の必要性が生

コード全体にわたる MPI 並列化作業が進展した

じたので、そちらの課題に集中して取り組んだ結後、相対論化作業に取りかかるつもりである。ま果、相対論化作業は準備の段階でとどまることにた、現行のコードで物質混合が起こりにくい問題なった。したがって、その後に予定していた並列化は、Price (2007) が提唱した artificial conductivity

7を導入することで解決しようと思う。

J0632+057”, Variable Galactic Gamma-

現在、平成 27 年度にコードの MPI 並列化と ar-ray Sources III（2015 年 5 月, ハイデルベ

tificial conductivity 導入を行い、平成 28 年度に相ルグ）対論化作業を行うことを考えている。

7.

(4) 国内会議発表

研究成果リスト

岡崎敦男 “Change of pulsar state in the

gamma-ray binary HESS J0632+057”, 日本研究課題に関連した今年度の研究成果は以下

本天文学会（2014 年 9 月、山形）の通りである。

岡 崎 敦 男 “Dynamical Modeling of the(1) 学術論文

Be/BH Binary AGL J2241+4454”, 日本
天文学会（2015 年 3 月、大阪）

Moritani, Y., Okazaki, A. T., Carciofi,(5) その他

A. C., et al., 2015, “Probing the nature of the TeV γ-ray binary HESSなし

J0632+057 by monitoring Be disk variabil参考文献

ity”, Astrophysical Journal Letters, 804:
L32

Aliu, E., 2014, Astrophysical Journal, 780: 168(2) 国際会議プロシーディングス

Bosch-Ramon1, V., Barkov, M. V., Khangulyan, D.,
Perucho, M., 2012, Astronomy & Astrophysics,

Okazaki, A. T., 2015, “Current status of our544, 59

understanding of Be disk physics”, ASP
Conference Series, Eds.

Price, D. 2007, Publications of the AstronomicalA. Sigut & C.

Society of Australia, 24, 159

Jones, 印刷中

Rosswog, S. 2009, New Astronomy Reviews, 53, 78(3) 国際会議発表

Saitoh, T. R. & Makino, J., 2013, AstrophysicalJournal, 768: 44

Okazaki, A. T., “Dynamical modelling of theSpringel, V. & Hernquist, L. 2002, Monthly Notices ofBe/Gamma-ray sources HESS J0632+057

the Royal Astronomical Society, 333, 649

and AGL J2241+4454”, Be/X-ray Binary
Workshop 2014（2014 年 7 月, バレンシア）
Okazaki, A. T., “Current status of our understanding of Be disk physics”, Bright
Emissaries: Be Stars as Messengers as StaDisk Physics（2014 年 8 月, ロンドン（カナダ））（招待講演）

Okazaki, A. T., “Quasi-cyclic variability of circumstellar disks around Oe/Bestars”, Magnetism and Variability in OStars（2014 年 9 月, アムステルダム）
O k a z a k i , A . T ., “Flip-flopping pulsar
model for the gamma-ray binary HESS8