ダークマターとは
夜空に煌めく無数の星々や，我々が住む地球の親 星である太陽は銀河系（天の川銀河）のメンバであ る．銀河系は半径数万光年程度のディスク状に分布 した数千億個の星々からなると考えられている．銀 河が多数群集まっているのが銀河団であり，それら はさらに超銀河団という集団構造を形成している． 超銀河団同士もフィラメントやシート状に分布する 銀河でつながっていて，宇宙の大規模構造を形成す る．このように宇宙は階層的な構造をしている． こういった我々の目で観測できる天体は，原子や 分子のような「バリオン」と呼ばれる物質で構成さ れている．最新の宇宙物理学は，バリオン物質は宇 宙の構成要素全体の4～5% 程度の割合でしか存在 しないことを明らかにした．残り25% 程度は「ダ ークマター」と呼ばれる重力のみ作用する物質，残 りは宇宙の加速膨張にかかわる「ダークエネルギー」 である．宇宙誕生からおよそ140億年間にわたって， どのようにダークマターが重力的な構造形成の主要 な役割を果たしてきたかを俯瞰してみよう． 宇宙誕生直後は，ダークマターは空間全体にほと んど一様に存在していた．しかしいたるところに，密 度が少しだけ高いところ，低いところ（密度揺らぎ） が存在した．ある程度高密度な領域では，重力によ って宇宙膨張に打ち勝つことができ収縮していく．こ うしてできる安定なダークマター天体が，ダークマタ ーハローである．そしてハローの中心の深い重力ポテ ンシャルの中にバリオンが集まる．バリオンが高密度 に凝縮することで原始星ができ，ハローが合体を繰 り返すことで，現在観測される銀河のような多種多 様な天体を形成していったと考えられている☆ 1．図 -1 にシミュレーションによって再現された，宇宙のダー クマター大規模構造が形成し進化する様子を示した． ダークマターの存在は，さまざまな観測で間接的 に確かめられている．間接的というのは，星や銀河 のように光を介して直接観測できないが，“見える 物質”による重力だけでは説明できない観測結果が 数多く存在するという意味である．1つの例として， 回転曲線問題というものがある．銀河の星々は銀河 中心に対しておおよそ回転運動をしている．回転速 度はその星の位置より内側にある物質の総質量によ っていて，物質が多いほど速度も大きい．ところが， 観測された星の回転速度は，星の総質量から推定さ れるものより大きい．これは見えない物質，つまり ダークマターが存在することを強く示唆している． 逆に回転速度からダークマターの量を見積もること ができ，銀河系は半径数十万光年にわたって広がっ た，1兆太陽質量程度の，楕円状のハローが作る巨 大な重力場の中心に存在していると考えられている． 理論的にも，バリオンの密度ゆらぎの大きさだけ では，現在の宇宙の構造を作るには不十分であるこ とが分かっている．よってダークマターは天体形成 に必要不可欠であり，逆にダークマターの分布がバ リオンの集まり方，つまり銀河や銀河団などの分布 を決定する．その分布には宇宙初期の密度揺らぎの 性質についてや，宇宙年齢，各物質成分の割合など の情報が刻み込まれている1）． それでは宇宙の構造形成，進化に重要な役割を担 ってきたダークマターとは一体何だろうか ? 実は ☆ 1 観測される銀河は，バリオン質量の数倍から数十倍にもおよぶダー クマターハローの中心に存在していると考えられている．石山智明
筑波大学計算科学研究センター「京」 の威力で宇宙の正体に迫る
─ダークマターの超大規模シミュレーション─
8
ダークマターの素粒子としての正体はよく分かって おらず，宇宙物理学的にも，素粒子物理学的にも重 要な未解決問題の1つである．有力な候補として「超 対称性粒子」である「ニュートラリーノ」が提案さ れている．ニュートラリーノは，自己対消滅をしガ ンマ線を放出するという特徴がある．このガンマ線 を検出することで，間接的にダークマターを検出し ようとする試みが世界中でなされているが，決定的 な結果は得られていない．
スーパーコンピュータ ・ シミュレ
ーション
宇宙の構造形成過程は，重力の非線形性が本質的 に重要な役割を果たしているため，そのダイナミク スを研究するには数値シミュレーションが非常に有 用である．宇宙空間の物理現象は，地上で実験して 検証することが不可能であり，また質量，空間，時 間スケールが非常に広大である．したがって，スー パーコンピュータによる大規模シミュレーションが 必要不可欠である．時にシミュレーションは「理論 の望遠鏡」とも呼ばれ，宇宙の研究において非常に 重要な位置を占めている． 特にダークマターは直接観測できないので，ダー クマター構造形成の研究は，シミュレーションによ って大きく進められてきた．手法としては，重力多 体シミュレーション，いわゆる N 体シミュレーシ ョンがよく用いられる．これは対象を N 個の粒子 としてモデル化し，粒子間にはたらく相互作用重力 を計算し，運動方程式を解いて，粒子の時間発展を 追うシミュレーション手法である．ダークマターの ほかに，星団や銀河，ブラックホールなど適用範囲 は非常に幅広い．ダークマターの場合，個々の粒子 は位相空間上のダークマターを離散化したものをあ らわす．シミュレーション粒子がダークマター素粒 子に相当するわけではない． 宇宙初期の密度揺らぎを持ったダークマターの分 布からスタートし，重力相互作用による時間発展 を追うものを，特に宇宙論的シミュレーションと 呼ぶ．無限遠方からの重力を取り入れるために周 期境界条件が用いられる．初期密度揺らぎは宇宙 背景放射の観測などから非常によく制限されてい る．対象とする天体によって，扱う領域はさまざま である．銀河の空間分布を調べたい場合は，観測 できる宇宙の全領域にせまる領域を扱う．銀河の 微細構造を調べたり，我々の銀河の形成過程を研 究するときは1つのハローだけを扱うことが多い． シミュレーションに用いる粒子数はきわめて重要 である．どれだけ広い宇宙をシミュレーションでき るかが決まり，銀河分布の精密さに影響する．また どれだけ微細なダークマター構造が分解できるかが 決まり，それは銀河系自身のモデルや，ダークマタ ー対消滅シグナルの予測の正確さに直結する． 図 -1 シミュレーションによって再現された，ダークマター構造の進化の様子．明るさはダークマターの空間密度を表し，明るいところ は密度が高い．宇宙初期（一番左のパネル）に存在した微小な密度揺らぎが，時間経過とともに（左から右）重力により発達していく．● 一番右のパネルは宇宙誕生から約 140 億年後の現在の宇宙の姿を表している．図中右下の数値は空間スケールを表し，1Mpc はおよそ 300 万光年である．いたるところに存在する楕円状の高密度領域それぞれが個々のダークマターハローである．重力はすべての粒子ペア間にはたらくが，粒子数 Nに対して，ペア数は N2に比例する．そのまま 計算すると，粒子数を10倍にしたときに計算量は 100倍となる．これでは粒子数を増やすのは非常に 困難である．そこでツリー法2）のような近似法が よく使われる．基本的なアイディアは，近傍の粒子 からの力は直接計算するが，十分遠方の粒子群から の力はまとめて多重極展開からの力として計算する というものである．ツリー法では計算量を O (N log N) に減らすことができ，粒子数を大きく増やすこ とが可能となる． 宇宙論的シミュレーションではツリー法をい わ ゆ る Particle-Mesh（PM） 法 と 組 み 合 わ せ た TreePM 法がよく用いられる（図 -2）．この方法では， ツリー法で計算するのは，ある粒子から見てカット オフ半径内☆ 2にある粒子からの力のみである．カ ットオフ半径より外側は， PM 力からの寄与のみと なる．カットオフの導入により，通信量を大きく減 らすことができる．PM 法では，まず粒子の位置と 質量から一様メッシュ状の密度場を計算し，FFT とポアソン方程式から一様メッシュ状のポテンシャ ルを計算する．そして差分と補間から，各粒子位置 ☆ 2 許容できる誤差と計算量の観点から，典型的には PM で用いるメッ シュ幅の 3 倍程度に設定する． 宇宙論的シミュレーションでは，前述のように大 規模化が必要不可欠であるため，TreePM 法を用 いた大規模並列コードの開発が世界のいくつかの グループで進められてきた．しかし，これまで約 1,000並列度程度までに対応した並列化が報告され ているくらいであった．我々のグループでも C++ による OpenMP+MPI 並列コード，"GreeM" を開 発してきた．「京」のような数千並列以上の大規模 環境での代表的なボトルネックには， （1）構造の時間発展によるロードインバランス （2） PM 法における，メッシュ構造の大規模全対全 通信☆ 3などがあるが，我々は新しい並列化手法を 開発し，8万ノードにもおよぶ「京」のフルシステ ム上で効率良く動作するアプリケーションの開発 に成功した．開発にあたっては，2014年2月現在， 理化学研究所計算科学研究機構粒子系シミュレータ 研究チームの似鳥啓吾研究員，同チームリーダーの 牧野淳一郎教授の協力を得た．
●
●
構造の発達によるロードインバランス
並列化は，シミュレーション空間を並列数で分割 し，各ノードに割り当てることで行う．宇宙論的シ ミュレーションでは図 -1から分かるように，時間 発展とともにいたるところに高密度なダークマター ハローが発達する．よって領域分割が空間一様で時 間不変である場合，構造が発達するとロードバラン スが非常に悪化する．多くのコードでは，morton ordering や peano
hirbert curve のような空間充填曲線に沿って粒子 を並べ，並列数で均等に分割して領域分割を行って いる．この方法ではロードバランスを均等にしやす い反面，領域形状が複雑化し，ステップ間での粒子 交換や，ツリー法による力の計算に必要な，他のノ ードに存在するツリー構造の通信粒度が細かくなる 傾向にある． ☆ 3 通信にかかわる全ノードから全ノードへのデータの相互交換． Tree 図 -2 TreePM 法の概念図．粒子分布を 8●分木で表現する（左パ ネル）．ツリー法で計算するのは，ある粒子から見てカットオフ 半径（緑色の円）内にある粒子からの力のみである．それより外 側は PM●のみで計算される．ツリー法では，十分に離れた粒子群 からの力はまとめて多重極展開として計算される（緑丸）．右下 のグラフは粒子間距離の関数としての，各力の大きさである．緑 線がツリー法で，青線が PM で計算される力である．これらの合 計がニュートン重力（破線）である．
そこで我々のコードでは，再帰的多段分割法に基 づいた3次元領域分割を採用した．3軸それぞれの 分割数を決め，まず x 軸方向に計算時間が均等に なるように分割する．そして分割された空間ごとに 独立に，y 軸方向，z 軸方向に分割していく．結果 領域形状は直方体となり，通信粒度が小さくなるこ とを防ぐことができる．領域分割の例を図 -3●に示 す．構造が発達する個所で細分化されている様子が 分かる．またこの領域分割は「京」のようなトーラ スネットワークにマップしやすく，通信最適化を促 進しやすいというメリットがある． 通信量を抑えるために，領域分割はサンプル粒子 に対して行っている．各プロセスは適当なサンプリ ングレートで粒子をサンプルし，ルートノードがそ れらを集める．そして領域再分割を行い，結果を全 体に送りかえす．このとき，各ノードは1つ前のス テップの計算時間に比例するように，サンプリング レートを決める．前のステップで少し計算負荷が大 きかったノードは，サンプルされる粒子数も少し大 きくなる．領域再設定は，各ノードでのサンプル粒 子の数が均等になるように行われる．したがってサ ンプル粒子数が大きかったノードは，次のステップ で計算負荷が減るように領域再設定が行われる． これらの方法により，各ステップで動的に，ロー ドバランスを調整しながら領域分割を行っている． さらにサンプリングによる浮動を最小限にするため に，過去数ステップの領域形状を平均化している． こうして，最小限の通信で8万ノードでも計算負荷 をほぼ均等に保つことに成功した．
●
●
FFT 計算のための，メッシュ構造の大規模
全対全通信
PM で力を計算するとき FFT を実行するが，我々 のコードでは1次元スラブ並列にのみ対応した FFTW ライブラリ☆ 4，または3次元並列に対応し た富士通製の FFT ライブラリを使用している．い ずれの場合でも，FFT 用の分割されたメッシュデ ータの構造は，ノード間で均一でなければならない． ところが領域形状は直方体であり，図 -3でも見ら れるように，3軸不等の度合いがかなり大きくなる こともある．したがって各ノードの担当領域を覆う ローカルなメッシュは，ノード間でサイズや形状が 大きく異なるが，それを均一なスラブ分割，また は3次元分割になるようにデータ構造を変換しなけ ればならない（FFT で計算されたメッシュ上のポ テンシャルに対しても同じ変換が必要）．計算全体 の中で FFT 自体の計算量はたいしたことはないが， このデータ構造を変換する通信は大規模で，非常に 不均一な全対全通信であるため，数万並列規模だと ボトルネックになり得る． 我々は全体をいくつかのグループに分け，まずグ ループ内での全対全通信を行ってからグループ間 で通信するというように，2段階に分けて通信を行 うことによって性能を大きく改善することができ た．基本的に FFT を実行するノード数は，全体の ノード数より小さい（特に1次元スラブ並列にのみ 対応した FFTW を用いる場合）．まず，各グルー プは FFT を実行するノード全体を内包できるサイ ズ，形状にし，FFT 用のメッシュ全体をグループ 内で分散して持つようにする．そして，グループの メンバノードが持つ部分メッシュを，グループ内全 ☆ 4 http://www.fftw.org/ 図 -3 ●領域分割の例．2 次元で 8 × 8 分割している．右上のパネ ルは中心部を拡大したものであり，高密度なハロー部分で細かく 分割されていることが分かる．を持っているものの，中身はグループが覆う領域か らの寄与のみなので不完全である．次にグループ間 通信を行い，各メッシュ上の値を足し合わせること で，全体の完全な FFT 用メッシュを生成する．こ うすることで，この通信部分を数倍高速化すること に成功した3）． これらの並列化手法を実装することで，「京」の フルシステム上でもほぼ完璧なロードバランスを実 現した． さらに粒子間の相互作用重力を高速に計算するこ とが重要である．富士通コンパイラ提供の組込み関 数を用いて手動 SIMD 化を行い，さらに力を受け る粒子とおよぼす粒子両方のループで，手動ループ アンロールを行った．最適な展開数を見つけるため に，展開数を指定すれば自動でコードが生成される ようにした．我々の実装では，1相互作用の計算に は51演算必要であり，FMA 部と非 FMA 部の割 合が1:1（17 命令ずつ）であるので，理論ピーク性 能は CPU ピーク性能の75% となる．この部分だ けの単純なベンチマークテストでは，72% 程度の 実行効率（理論ピークの96%）を達成した（この 部分は似鳥氏の実装による）． 図 -4に GreeM のストロングスケーリングの測定 結果を示した．粒子数1 兆のシミュレーションでは， 「京」のフルシステムまでほぼ理想的にスケールし ていることが分かる．最終的な実効性能は，2兆ダ ークマター粒子のシミュレーションを，「京」のほ ぼ全システムを用いて，5.67ペタフロップスの実効 性能（実行効率55%）で達成した．この成果が認 められ，「ハイ・パフォーマンス・コンピューティ ングに関する国際会議 SC 12」（2012年11月，米 国ソルトレークシティー開催）において，ゴードン・ ベル賞を単独受賞した3）． ゴードン・ベル賞ファイナリストには，ピーク性 能が20ペタフロップス（「京」の約2倍）の「セコ イア」（米国ローレンス・リバモア国立研究所）を 使い，同様のダークマターシミュレーションで14 ペタフロップス（実行効率69%）を達成した米国 のグループがあった4）．ところが我々のグループの アプリケーションが実際の計算速度で上回り，ピー ク性能が約半分の「京」を用いたにもかかわらず， 1粒子あたり2.4倍の速さでシミュレーションをす ることが可能であった（米国と同じ計算機を用い た場合は5倍近く速い）．こういった点が評価され， ゴードン・ベル賞の受賞につながったと考えられる． 我々のアプリケーションが圧倒的に速いのはいく つか理由があるが，通信部分のアルゴリズムを工夫 し並列性能を上げたのに加え，優れた重力計算アル ゴリズムを用いているのが効いている．米国のグル ープのアプリケーションでもツリー法を用いている． 我々のものは，ある粒子が受ける力はあわせて数千 個程度の質点（粒子と重心）からの力として計算す れば，必要な計算精度としては十分であった．とこ ろが，米国のグループは，我々より5～6倍多い数 を必要としていた．つまり我々のアプリケーション は，遠方の粒子群のまとめ方が優れており，そもそ も必要な演算量が1/6～1/5で済んでいたのであ る．一方，力の計算の部分は浮動小数点演算の密度 が高いため，実効効率は計算速度とは逆に米国のグ ループの方が高くなっている．だがこれは無駄な計 算をしているとも言える．単に大規模なスーパーコ ンピュータ上で実行効率を高めるのではなく，より 優れたアルゴリズムを実装し，実アプリケーション 1 10 10 2 10 3 10 4 10 5 2,048 3 4,096 3 8,192 3 10,240 3 ノード数 1ス テ ッ プの計算時間 （ 秒 図 -4 「京」上での GreeM のストロングスケーリング．横軸が使 用したノード数（= 並列数．8 倍が CPU コア数），縦軸が 1● ステ ップにかかった計算時間を表す．破線は理想的なスケーリングが 達成された場合の性能である．カラーの違いは評価に用いた粒子 数の違いである（2,0483●〜 10,2403粒子）．
としての性能が優れていたことが評価され，我々の 劣勢を覆した勝利につながったと考えられる．
「京」 による世界最大規模のシミュ
レーション
「京」コンピュータの非常に強力な演算能力を最 大限活かすことのできるシミュレーションコードは 完成した．そして我々は2種類の世界最大規模の ダークマターシミュレーションを達成した（図 -5）． 左パネルは現在の宇宙の大規模構造を5,500億粒子 を用いてシミュレーションしたものである．右パネ ルは宇宙初期に形成するダークマター構造を690億 粒子を用いてシミュレーションしたものである．い ずれも，それぞれのスケールのシミュレーションで は世界最大粒子数である． 両者は似ているが，いくつか大きな違いがあるこ とが分かる．宇宙初期に形成するダークマター構造 は，現在に比べ滑らかな構造を示していることが分 かる．フィラメント構造の中や，ハロー周囲に存在 する小さいハローの数が少ないことも分かる． これらの違いは一体どこからきたのであろうか ? それは対象とするスケールの違いである．ハローは 密度揺らぎが成長してできるが，ダークマター素粒 子が固有に持っている自由運動によって，あるスケ ール以下の密度揺らぎは減衰し，そのスケール以下 のハローは形成しない．ダークマター素粒子として， 質量が100GeV 程度のニュートラリーノを考えた 場合，対応するスケールは10-6太陽質量程度である． 左パネルのシミュレーションはそれよりはるかに大 きい銀河や大規模構造の形成を対象としている一方， 右パネルのシミュレーションはこのスケールを対象 としているので，できるハローの数が少なく，構造 が滑らかになる．なお両スケールは質量で20桁以 上の差があるので，両方を同時にシミュレーション することは，エクサスケールをもってしても不可能 である． こういった大規模シミュレーションから何が分か るのだろうか ? まず，左パネルのような宇宙の大 規模構造形成シミュレーションからは，銀河や銀河 中心の超巨大ブラックホールの空間分布を予測する ことができる．超巨大ブラックホールは銀河と共進 化していると考えられているが，その形成，進化過 程は謎に包まれている．シミュレーションから得ら 図 -5 「京」による世界最大規模のダークマターシミュレーション．左パネルは，現在の宇宙の大規模構造を 5,500 億粒子 を用いてシミュレーションしたものであり，およそ 50 億光年にわたる空間の銀河分布を再現している．白枠内は拡大図で あり，右下のパネルは銀河団サイズ（約 1015太陽質量）のハローである．右パネルは，宇宙初期に形成するダークマター構 造を 690 億粒子を用いてシミュレーションしたものである．宇宙誕生からおよそ 1 億年後の宇宙であり，右下のパネルは約 10-3●太陽質量のハローである．ルの性質を理論的に予言することができる．「京」 以前のシミュレーションでは，分解能と体積の両方 が足りず，銀河の階層的形成を追いつつ，最新の観 測データと比較できるような広大な領域にかけて， 銀河とブラックホールの分布両方を予測することが できなかった．今回「京」の威力によりそれがはじ めて可能となり，共進化過程の解明への手がかりな どが得られるようになった． 右パネルのシミュレーションからは，銀河系が存 在するハローの微細構造に関する手がかりが得られ る．ここ十数年の研究により，ハローは中心が高密 度であり，さらにその中には図 -5でも見られるよ うな無数のより小さいハローが存在することが分か ってきた．また最近の研究は，小さいハローほど内 部質量の中心集中度が大きく，大きいものと比べて 多量に存在することを示唆している．したがってこ の高密度な最小のハローが，地球近傍に存在し得る かどうかが，ニュートラリーノの探査とその詳細な 性質を解明するために重要である．なぜなら，ダー クマター対消滅起源のガンマ線検出のために，地球 からできるだけ近くて，ダークマターが集中して存 在すると考えられる領域に観測衛星を向けたいから である．だが「京」以前は最小のハローのみを正し くシミュレーションできただけであり，最小のハロ ーがより大きいハローの中でどう分布するかが分か らなかった．今回のシミュレーションでは最小ハロ ーが合体してできるハローを追うことができ，銀河 サイズまで成長したとき，そのなかでの最小ハロー の分布をある程度推定することが可能となった．
エクサ時代の構造形成シミュレーシ
ョン
きたるエクサスケールの時代では，我々は「京」 の100倍の演算能力を手に入れることができるだ ろう．粒子数が100倍になれば，さらに色々面白 いサイエンスを展開できるようになる．たとえば， とはいえ，宇宙の一部分である．宇宙のサイズにせ まる領域をシミュレーションできるようになること で，我々の宇宙の構造形成モデルをより精密に決め られるようになるだろう．また銀河系のハローのよ り微細な構造を明らかにすることができ，ダークマ ター探査のための，より具体的なモデルを提案する ことができるようになるだろう． またダークマターシミュレーションの規模を上げ るだけではなく，バリオンを同時にシミュレーショ ンすることも重要である．今後は，今回開発したア プリケーションを，圧縮性流体のための代表的な計 算法である SPH 法5）などのより複雑な物理を扱う アルゴリズムと組み合わせ，星形成，銀河形成など， 宇宙の構造形成過程に関する科学的成果の創出に向 けていきたいと考えている． 参考文献 1） 宇宙論 II : 宇宙の進化，シリーズ現代の天文学，日本評論社 （2007）．2） Barnes, J. and Hut, P. : A Hierarchical O (N log N) force
-calculation Algorithm, Nature, Vol.324, pp.446-449 (Dec.
1986).
3） Ishiyama, T., Nitadori, K. and Makino, J. : 4.45 pflops
Astrophysical N-body Simulation on K Computer : The
Gravitational Trillion-body Problem. In Proceedings of the
International Conference on High Performance Computing, Networking, Storageand Analysis, SC'12, pp.5:1-5:10, Los
Alamitos, CA, USA, IEEE Computer Society Press (2012).
4） Habib, S., Morozov, V., Finkel, H., Pope, A., Heitmann, K., Kumaran, K., Peterka, T., Insley, J., Daniel, D., Fasel, P.,
Frontiere, N. and Luki´c, Z. : The Universe at Extreme Scale
: Multi-petaflop Sky Simulation on the bg/q. In Proceedings
of the International Conference on High Performance Computing, Networking, Storage and Analysis, SC'12,
pp.4:1-4:11, Los Alamitos, CA, USA, IEEE Computer
Society Press (2012). 5）シミュレーション天文学，シリーズ現代の天文学，日本評論 社 （2007）． （2014 年 3 月 1 日受付） 石山智明 [email protected] 1982 年生．2005 年東京大学教養学部広域科学科広域システム科学系 卒．2007 年同大学院修士課程修了．2010 年同大学院博士課程修了．学 術博士．日本学術振興会特別研究員を経て，2011 年より筑波大学計算 科学研究センター研究員．宇宙の構造形成，HPC の研究に従事．2012 年 ACM Gordon Bell 賞受賞．天文学会．
