課題 4 報告

牧野淳一郎
東京工業大学理工学研究科
理学研究流動機構(4/1 から「地球生命研究所」)
HPCI 戦略プログラム分野 5 全体シンポジウム（2013/3/5)

話の順番• 課題 4 の話 (あんまり去年と変わらないところが多い)• その辺に関係するサイエンス・計算機スキーム等の話

最終成果• ダークマターハローシミュレーション、銀河形成シミュレーション、及び星形成、惑星形成シミュレーションについて、「京」及び近い将
来のスーパーコンピューターの性能を有効に生かすことができるアル
ゴリズム、シミュレーションコードを開発すること• 開発したシミュレーションコードを使った大規模シミュレーションで、– ダークマターハローの最小サイズハローからの形成過程– 銀河形成過程、特に、銀河の多様性の起源を明らかにすること• 星形成、惑星形成過程についても研究を進め、特に系外惑星の多様性の起源の理解を進めること

最終目標 — 「数値目標」(1)• ダークマターハローシミュレーションでは、理論ピーク性能の 25-50の性能を「京」全体に近いスケールの並列計算で実現する• ダークマターハローシミュレーションでは、大空間スケール (100Mpc程度)、小スケール (100-1000pc) の両方について、少なくとも 81923 、可能であれば 163843 粒子程度の分解能でのシミュレーションを複数実行する。
• 銀河形成については、理論ピークの 10-20% の性能を「京」全体の1/10 程度の規模の利用で実現する• 銀河形成については、流体粒子で 108 以上、質量分解能 (SPH 粒子質量) で 100-1000 太陽質量程度の計算を複数行ない、銀河の多様性の起源を解明する。また、可能であれば流体粒子で 109 以上の粒子数で、銀河群、銀河団スケールのシミュレーションを行う。

最終目標 — 「数値目標」(2)• 星形成・惑星形成シミュレーションについては、流体部分については銀河形成シミュレーションコード、また微惑星からの集積過程についてはダークマターシミュレーションコードの成果を利用したソフトウェア開発を進める。

年次計画• H23: ダークマター、銀河形成コードについて、「京」向けの最適化• H24 まずダークマターシミュレーションについて大規模計算を行う。
銀河形成についてはコードの改良、小規模なテスト計算を進める。星・惑星形成についてはコード開発を進める。
• H25 以降 ダークマターシミュレーションを継続する。銀河形成については大規模計算を進める。H26 前後から星・惑星形成についてもプロダクトランに入る。
課題 4 でやりたいこと• ダークマターハローから惑星形成までの (主に) 重力による構造形成の(輻射・磁場が重要でないというわけではない)• 空間・時間的にアダプティブな大規模並列コードでのシミュレーションによる理解

年次計画の達成状況とか
今年度:• まずダークマターシミュレーションについて大規模計算を行う:計算終
わってないが始めてはいる• 銀河形成についてはコードの改良、小規模なテスト計算を進める。:
輻射入りコードの「京」むけ最適化はだいぶ進んだ• 星・惑星形成についてはコード開発を進める。:研究員採用したくらい、、、

実施体制
組織

• 東工大 牧野、井田 — 斎藤、馬場• 筑波大学 梅村 — 石山、岡本、長谷川• 千葉大学 松元• 国立天文台 富阪、小久保、大須賀 — 高橋
テーマと人ダークマターハロー: 石山、牧野、梅村
銀河形成: 斎藤、馬場、岡本、長谷川、牧野、梅村星・惑星形成: 富阪、小久保、井田
ブラックホール成長: 高橋、大須賀
今年度末辺り異動: 岡本、斎藤
来年度新任: 台坂淳子 (惑星形成)

やっているテーマ• ダークマターハロー: 石山君がするので省略• 銀河形成: 今日はこっち関係の話を少し• 星・惑星形成: 今日は省略• ブラックホール成長: 高橋君がするので省略
作業部会資料

銀河形成の計算の話

銀河形成の計算の話をしようと思ってたんですがちょっと気分変えて怪しい話を。

銀河力学と気候変動• 導入: スベンスマーク仮説と銀河-太陽相互作用• 銀河のスパイラル構造• 太陽の銀河内運動を遡る• まとめ

導入: スベンスマーク仮説と銀河- 太陽相互作用• スベンスマーク仮説• 銀河渦状肢と太陽スベンスマーク仮説
Svensmark 2007 から
基本的には、銀河宇宙線の地球までふってくる量が増えると雲が増えて寒くなる、という話
宇宙線が増えるメカニズム• 地球磁場の変化• 太陽風の変化• 銀河宇宙線自体の変化: 近傍での超新星爆発の増加とか

長周期の気候変動• 1.4 億年くらいの周期の気候変動がある(牧野はよく知らないので詳しくは知っている人に聞いて下さい)• 長周期の起源: 地球の内部や軌道運動ではなさそう (???)• 銀河の渦状肢を通過すると宇宙線は増えるのでは？

銀河系と太陽• 渦状肢は定常密度波で、太陽の位置の円運動とは違う角速度で動いている• なので、太陽がほぼ周期的に渦状肢を横切る• 渦状肢のところでは、星間ガスが圧縮されて活発な星形成が起きている• 宇宙線が多くなっていて寒冷化

いくつもの疑問• そもそも星形成が高いくらいで本当に寒冷化なんかするのか？(今日はこの話はしません。すみません)• 渦状肢って本当に一定のパターンのがあるのか？
• 太陽の運動はどんなふうなのか？

銀河の渦巻構造の最近の理解• 国立天文台ではスベンスマーク仮説をなんとかしようと思ってたわけでは全くなく、純粋に銀河円盤の力学、渦巻構造の起源の研究をしていた (馬場、斎藤、和田、牧野辺り)• 観測でも、VLBI 観測で銀河円盤の中の星の運動が少しみえてきた30 年くらい前から教科書にのってる理論— 定常密度波理論• 左の図のように、「何故か」星は楕円軌道をそろって回っていて、その軸が半径によってずれると「仮定」する• そう仮定すると時間たっても消えない渦巻構造はできる• 「何故か」は不明なんだけど、、、• 他の説も色々あるけどこれが正しいという決定的なものはなかった• 数学的には色々議論できるのは密度波理論

最近の観測• VLBI (超長基線干渉計) で遠くの星の距離と動く速度を三角測量で決める• 結構円運動からはずれて動いている• 定常密度波理論はこういうことはな
いはず
はて？

最近の (まあちょっと前の、、、) シミュレーション(馬場他 2009)
animation123)
空間構造と、円運動からのずれ

左:星の分布右: 冷たい星間ガス

観測と比べると結構それっぽいものができるようになってきた。

「スベンスマーク仮説」は？
• 渦状肢は定常ではないし、ケプラー速度と違う「パターン速度」があるわけではない• なので、「太陽と渦状肢の周期的遭遇」はない• 本当のところはどうか、我々の銀河系に近い (「太陽」から観測すると大域構造が非常に近い) シミュレーションモデル
で太陽に近い速度の星の運動を過去にさかのぼってみた。

結果:銀河系と太陽

青:寒冷期
寒冷期には銀河中心に近い？これはフリーパラメータなしで位相まで一致。

太陽周りの環境変化• 上: 紫外線での明るさ (星形成率を表す)• 中: 星間ガスの密度• 下: 超新星発生率• これらが高いと寒冷期になる？

まとめ• スベンスマーク仮説: 銀河渦状肢との周期的遭遇で気候変動• 現代的な銀河円盤シミュレーションではこういうことはおこらない• が、太陽のエピサイクル運動による銀河中心からの距離変化が周期的環境変化を起こしている• この周期は気候変動の周期と実際に一致していて、関係している可能性はある。
• バーとの相互作用も考えると 6-10 億年スケールでの変動もありえる• 「銀河古気候学」みたい感じのことができるかも。

SPH のもっと怪しい話• 去年のシンポジウムでは斎藤君が DISPH の話をした• これはこれで色々上手くいくが問題も残る• さらに拡張した定式化を考えた• まだあんまりテストはしてないSaitoh and Makino 2012
従来の SPH1
新しい SPH1
従来の SPH2
新しい SPH2
接触不連続を扱える SPH の定式化 (Saitoh and Makino, まだ論文通ってない、、、)

SPH と接触不連続、KH 不安定
Agertz et al (MN 2007, 380, 963)
• SPH と Grid コードで、「Brob test」の答が全然違う• もっと簡単な Kelvin-Helmholtz 不安定 (計算は 3 次元)
でも全然違う• SPH だめじゃんどれくらい違うか (1)• 周りよりつめたい (温度1/10、、密度 10 倍) ガスの球を超音速で動かす• 上から 3 個は Grid
• 下の 2 つは Gasoline(下は 10M 粒子)• SPH では 境界での不安
定が起きないで、冷たい流体が固まりのまま。

どれくらい違うか (2)SPH では KH 不安定が起きない。

どれくらい違うか (3)2 流体の境界面で妙な隙間ができる。このため力が働かない？

何故そんなことが起きるか？
基本的には、 SPH の運動方程式自体の問題、ということになっているらしい。
SPH の基本式おさらい
密度の推定ρ(x) =∑

mj W (x − xj ),(1)f (x′)W (x − x′)dx′.

(2)jある物理量 f の推定⟨f ⟩(x) =∫SPH の基本おさらいのつづき (1)f の微分: ⟨∇f ⟩ = ∇⟨f ⟩ で、以下の恒等式
1=

∑mjj

1
ρ(x)

W (x − xj ).

(3)を使って、さらにもうちょっと近似して⟨∇f ⟩(x) ∼∑jmjf (xj )ρ(xj )∇W (x − xj ).

(4)

SPH の基本おさらいのつづき (2)運動方程式は − ρ1 ∇P を計算する。この時に恒等式

1ρ

∇P =PPρ2∇ρ + ∇ 2 .
ρ

(5)を使って対称化すると
v̇i = −になる。

∑jm
j

Pi
ρ2i



+Pj ∂
ρ2j

∂xiW (xi − xj ),(6)密度不連続面での振る舞い
通常の SPH では、ここまでの変形の 2 箇所で ρ の微分可能性を仮定している。以下の 2 つの「恒等式」である
1∑1 = mjW (x − xj ).
(7)
ρ(x)
jP
1P
∇P = 2 ∇ρ + ∇ 2 .
(8)ρρρ
SPH でカーネル推定した密度は滑らか。このため、• 接触不連続の低密度側では、密度を過大推定、高密度側では過小推定する• その結果、圧力、その空間微分もデタラメになる。結果として粒子が再配置される

対策「根本的」な理由:
ρ は滑らかだけど u (内部エネルギー) はジャンプがあるまま。
この観点では、u を滑らかにすればよい。色々提案あり。

• u にもカーネル推定した量を使う• u を拡散させる (人工熱伝導)• 質量密度でない密度 (数密度とか) を使う
それぞれ、それなりにうまくいくケースもある。

新しい提案 — 思想
圧力が本来変わってないのに、密度が不連続なだけでおかしなことが起こるのは何故か？
物理量 (とその微分) の推定式に密度を使うから:⟨f ⟩(x) =∑mj f (xj )W (x − xj ).

(9)ρ(xj )
ここでやっていることは、本質的には体積要素 dx を mj /ρ(xj )で置き換えているだけ。
j

粒子の占める体積の推定さえできれば別に何を使ってもいいはず (但し、色々論文をみてもこれ以外の方法をやっているものはみたことがない)新しい提案 — 原理
質量密度の代わりに何を使うか？
気体 (理想気体) は状態方程式 P V = nRT で規定される。
ここにはそもそも質量密度とかない。右辺は本質的には熱エネ
ルギー (に比熱による係数かけたもの)内部エネルギー密度 (結局圧力と同じ) を使えばどうか？
粒子当りの内部エネルギーは今の SPH でも時間発展させる量なので、単位体積当りの内部エネルギー、つまり圧力の空間分布は質量密度を使わなくても計算できる。
接触不連続では圧力は (もちろん) 連続なので、変なことは起きないかもしれない。
これはこれでいいけど、非理想気体は？

非理想気体への一般化の思想
状態方程式が理想気体と違う時にどうすればいいか？

• 基本的な思想は接触不連続で自然に連続なものを連続にして、それを体積推定に使う。
• そうすると、自然な熱力学的量はそもそも圧力 p であって内部エネルギー密度 u じゃない
• よく考えると、別に p でなくても、例えば p だけの単調な関数であればいい。例えば、y = pα みたいなの。
内部エネルギーの代わりに、「圧力の適当な関数に粒子当りの体積を掛けた量」yV を基本的変数にしてみる。理想気体の場合は比例係数が変わるだけで同じ式になるように気をつける。

圧力をさらに一般化する意図• 圧力は変化幅が大きい。これは strong shock で問題になる。つまり、低圧力側の量が計算にきいてこないので、色々おかしなことになる。
• 圧力に依存するが変化幅の小さい量にすればこの問題は回避できそうな気がする• 表面でも問題が少ないはず。
• 本当かね?まあとりあえず式を導出してみる非理想気体への一般化 (1)y は、 p だけに依存する単調で微分可能な関数とする。
粒子当りの yV 項を Zi と書く。この量の空間密度を
y=∑j

Zj W (x − xj ).

(10)で定義する。そうすると、他の物理量のカーネル推定は⟨f ⟩(x) =∑Zj f (xj )jy(xj )⟨∇f ⟩(x) =∑Zj f (xj )jy(xj )W (x − xj ),

(11)

空間微分は

∇W (x − xj ).

(12)一般化 (2)エネルギー方程式が∑ pi Zi Zj
U̇i =(v i − v j )∇W (xi − xj ).
yi yj
j

(13)

運動方程式

miv̇i = −∑ZiZjjyi yj(pi + pj )∇W (xi − xj ).

(14)

Z の方程式Żi = Zi(αγi − 1)

∑Zj

jyj

(v i − v j )∇W (xi − xj ).

(15)

一応これで方程式系は閉じるが、これだと U と Z がコンシステントにならない。
エネルギー方程式も積分して、内部エネルギーの値がそっちになるように Z を補正することを考える

Z の補正と y の具体的な形
Z, y を求めるのに以下の反復をする。
Zi,new = f (yi,old, Ui)∑yi,new =Zj,new W (x − xj ).

(16)
(17)

j

最初の y は Z を時間積分して求めたものを使う。
そもそも y に何を使うか？

• なんとなく自然な形は y = pα で、 α は正で 1 より小さい (1/5 とか 1/10) に取ること？
• 意味: 変化幅を小さくする• 実験必要、、、

と思ってたわけですが、• ここで極限を考える (とにかく極限でどうなってるか理解しろ!というのが先人の教え)• 色々ありえるが、 α = 0 の極限を。
• そうすると y = 1• つまり、∑j

Zj W (x − xj ) = 1.

(18)

この式はなに？
∑Zj W (x − xj ) = 1.

(19)

j

• Zi は粒子があるところでは恒等的に “1” になる量を与える。
• ということは、「体積要素そのもの」ということ？体積と思うなら Viと書くことにしよう。というわけで、∑Vj W (xi − xj ) = 1

j• そもそも Vi は求まるか？安定か？ → OK(今日は省略)• これで SPH の式はでてくるか？ → OK
• それはなんかいいことがあるか？ → 多分ある

(20)

これで SPH の式はでてくるか？
特に問題なくそれらしい式がでてくる。
エネルギー方程式は多分U̇i =∑jpiViVj (v i − v j )∇W (xi − xj ).

(21)

運動方程式はmiv̇i = −∑jViVj (pi + pj )∇W (xi − xj ).

(22)というようになる。V の時間発展も Z と同様。

V̇i = −Vi∑jVj (v i − v j )∇W (xi − xj ).

(23)

SPH のところのまとめ• 従来の SPH の問題点である、接触不連続での異常な振舞いは、密度の代わりに圧力を基本的変数にすることで解消できる• 但し、これは、強い衝撃波や、圧力 0 の表面で問題。また、非理想気体の扱いが不明• 非理想気体で、圧力の任意の関数を基本的変数に取る形に拡張した• さらに、「圧力の任意の関数」ではなくて「定数」を使う SPH を導出した。
• 書いてみると「もっとも自然」な形。なぜ今まで知られてなかったのか？
• これからテスト。

まとめ• ダークマターハロー計算: コード開発は大変素晴らしい性能がでるところまで進んだ。 2012 年ゴードンベル賞• 実際の計算はまあその的• (実際に使っている人は結構疲弊している感じが、、、)• 他のサブテーマの実装・チューニング状況は色々• 惑星形成のコード実装を進めたい