現在、世界中では8m級の望遠鏡が稼動し、様様な銀河探査計画が進行または計画されている。今後人類は大量の観測データを得ることになるわけだが、近傍宇宙の観測量ならば自明な解釈を与えらるとができた観測量でも、深宇宙では常に選択効果の影が付きまとい、その解釈が恣意的なものになってしまう恐れがある。特に今まで銀河の形成と進化を定性的に理解する上で大いに役に立ってきた単純なモデルを用いる場合、その適用限界には注意を払う必要がある。では従来より客観的な解釈を与えるためにはどうしたらよいのか。観測データの定量的解釈を行っていくには今後はできる限り仮定を省いたモデル、つまり、大規模数値シュミレーションがもっとも強力な道具となることは間違いないであろう。

　銀河の観測量にはたくさんあるのだが、なかでも我々が興味を持っているのは銀河の形態がどのようにしてできたのかという問題である。この問題は銀河天文学が幕を開けるとともに常に問われ続けられてきた問題であるにもかかわらず、未だにその明確な描像は得られていない。これは銀河形態が、銀河の環境、内部力学的構造、星形成史など、銀河より小さいスケールから大きいスケールの現象が複雑に絡み合って決定されるからである。この問題を銀河の分布から銀河の内部構造まで再現するシミュレーションを行うことによって解くために、私は質量分解能、重力分解能の双方に於いて高いダイナミックレンジを有する数値計算コードを開発することにした。

　まず、銀河形態が判別できる計算をする為に必要な質量及び重力ダイナミックレンジを、シミュレーションで作られた銀河ハローの典型的な力学的緩和時間と銀河ハローの典型的な年齢を比較することにより見積もった。ここで、質量ダイナミックレンジというのは系全体の質量に対する粒子の質量の比を指し、粒子の質量が全て等しい場合は粒子数となる。一方、重力ダイナミックレンジとは系全体のスケールに対しての重力の計算がNewton則に合う最も小さいスケールの比を表す。上記の見積もりから質量ダイナミックレンジは109、重力ダイナミックレンジは105を実現する数値計算コードが必要であることが分かった。

　そこで、まずParticle-Mesh(PM)法を用いた計算コードを用いることを考慮したが、PM法には大きな欠点がある。それは、重力ダイナミックレンジが小さいことである。PM法を用いた典型的なシミュレーションでは使用する格子数と粒子数を同程度にする。この場合、粒子数が109なら、重力ダイナミックレンジは〜103となる。この値は要求されている値より2桁も小さい。PM法の重力ダイナミックレンジを上げる方法としてParticle-Particle Particle-Mesh(P3M)法が今日広く使われている。P3M法では、ある粒子にかかる重力を近傍粒子からの力と遠方粒子からの力とに分離して、前者を直接総和法で、後者をPM法で分けて解いている。しかし、ダークハローが形成され系全体の密度コントラストがあがってくると、近傍粒子間の直接総和計算部が近傍粒子数の自乗に比例して急速に増大するため、大きな密度コントラストを持つ系の計算を実行することが困難になる。よって、CDMモデルのように小さい質量スケールほど初期揺らぎが大きい場合、同じ粒子数でも計算領域が小さいものほど計算時間を要する。

　このような問題をすべて克服しているのがAdaptive mesh refinement(AMR)を用いたPM法である。PM法の重力分解能、即ち、γ-2則に合う最も小さいスケール、はPM法で用いられる格子間隔によって決まっている。PM法ではこの格子は等間隔の立方格子である。しかし、AMRを用いれば、より高い重力分解能が要求される領域に対し局所的に細かい格子を配置することによって、高い重力ダイナミックレンジを得ることができる。更に、各粒子の軌道計算の時間発展の間隔を各粒子が属する格子間隔に比例する形で変化させることによって、格子を空間的に分割するだけでなく時間方向にも分割することができる。私はこのAMRを用いたPMコードを開発するとともに、種々のテストを繰り返すことによって、このコードが高い重力ダイナミックレンジをもつことを示した。

　一方、N体コードではダークマター成分の時間発展を追うことができるのだが、直接的に観測されるのはガス成分とそれから生じる星の成分である。その為、数値シミュレーションで直接的に観測データと比較するためには、N体コードと数値流体コードを組み合わせる必要がある。このAMRを用いたN体コードには自然にEuler的な数値流体コードを取り込むことができる。Euler的な数値流体法はPM法と同様に格子間隔で分解能が決まってしまうのだが、AMRを使うことによって局所的に分解能をあげることができる。また、N体成分と流体成分の重力計算は同時に解くことが可能であるAMRを用いたEuler的コードは天文以外の分野も含め様々な分野で開発が行われているが、AMRを宇宙論的な問題に用いた数値シミュレーションコードは未だ少ない。我々は2次のGodunov法に基づいたEuler的な数値流体コードにAMRを施したコードを開発した。

　流体計算コードの完成に先立ち、現在はN体コードの計算結果の解析も進めている。まず、シミュレーションの結果から銀河団や銀河のダークハローを抽出するハロー探索のコードを作成した。このコードではFriends-of-friends法を用いている。また、ハロー探索を各時刻でのN体シミュレーションの出力結果に用い、ある時刻とその前の時刻でのハローの同定を行うと、現在の、即ちz=0でのダークハローがどのような合体を経て今日に至ったかを示すマージングツリーと呼ばれる樹状図を作ることができる。このマージングツリーを用いると銀河形成の準解析的モデルを通して、本来質量分布しか表していないN体計算の計算結果に、実際の銀河の各種パラメータを空間分布とともに付与することができるようになる。これは即ち、数値的な銀河カタログを作成できるということである。このような数値銀河カタログは今後の深宇宙探査、広域探査、赤方偏移探査といった観測データを解釈するうえで基本的なデータになると考えられる。

　2001年1月より国立天文台ではベクトル並列型計算機VPP5000を中心とした新スーパーコンピューターシステムが稼動し始めた。この計算機の性能を十分に発揮するためには上で述べたAMRコードをベクトル化しなければならない。AMRコードはループ内でアクセスされるデータの種類から次の三つに分類される。即ち、粒子ループ、格子ループ、粒子−格子ループである。このうち、粒子ループと格子ループのベクトル化は容易であるが、粒子−格子ループのベクトル化には工夫を要する。まず、同じ格子に入っている粒子は連結リストでつながっているが、粒子−格子ループにおいてこれを一つ一つ手繰ってしまうと(深さ優先走査)ベクトル化できないのだが、これを幅優先走査に変換することによってベクトル化することができた。また、粒子の質量を格子に付与する個所は従来ベクトル化することができないと考えられていたが、上記の走査順序変換に加えループ分割を行うとベクトル化することができることを示した。

　更にスーパーコンピュータの性能を引き出すにはベクトル化に加え並列化を行う必要がある。その為には各階層毎に各プロセッサでの負荷が均等になるような格子分配を実装した。分割をする際はMorton順序と呼ばれる方法で格子を並べ替え、整列された格子を各プロセッサに配分する。こうするとあるプロセッサの担当する格子が比較的集まるように分配することができる。これらに加え、階層間通信、同一階層間通信、粒子−格子間通信を実装し、世界初のAMRを用いたN体コードの分散並列コードは完成した。

　我々はこのコードを用いて銀河団の形成シミュレーションを行った。その結果、粒子の分布は他のグループの計算結果とおよそ一致するのだが、銀河団ハローの中心付近での力学構造が顕著に異なる結果となった。我々のコードは他のグループと比べ時間分解能が1,2桁高くなっている。そこで、我々は改めて時間分解能の低い計算を行ったが、こちらは他のグループの結果と一致した。このことは、宇宙論的シミュレーションにおいて、ダークハローの力学構造を精密に計算するためには従来よりも短い時間刻みで計算しなければならないことを示唆している。

　我々はAMRを用いた宇宙論的シミュレーションコードの開発を行った。このコードは期待通りの性能を出しており、今後はこのコードを用いて数値銀河カタログを作成し、実際の観測データと比較をしながら、銀河生誕の秘密を解き明かしていきたいと考えている。

　本論文は、宇宙論的な構造形成の数値シミュレーションのための新しい高速な計算アルゴリズムを提案し、さらに国立天文台の並列ベクトル型計算機VPP-5000上に実装し、計算精度・速度を従来の計算コードと比較したものである。宇宙論的な初期ゆらぎからの重力不安定による構造形成の研究、特に銀河や銀河団の分布や構造の観測との比較のためには構造形成の非線形段階を追跡できる数値シミュレーションが有効なツールになる。CDM宇宙論の元では重力場は無衝突物質とみなせるダークマターに支配されているので、基本的な構造形成を追うためにはダークマターのみによるいわゆるN体シミュレーションを行えばよい。しかし、銀河、銀河団の観測と比較するには準解析的なモデルと組み合わせてダークハローを銀河と解釈しなおすか、あるいはバリオンの進化を流体力学を解いて計算し、さらに輻射輸送や星形成過程をモデル化することで銀河の観測可能な構造、光度、スペクトル等の情報をシミュレーションから直接に引き出す必要がある。本論文で提案される計算法はこのうち後者の方向を目指すものである。

　従来、このような宇宙論的N体、あるいはN体+流体のシミュレーションには基本的には2種の方法が使われてきた。一つは重力ポテンシャルは規則格子上で解くものである。この場合にはソルバーとしてFFTが利用でき、極めて高速な計算が行える。しかし、空間分解能が格子サイズに制限されるため、宇宙膨張から切り離された銀河等構造の発展を追跡するのは困難であるという大きな欠点がある。格子サイズ以下は粒子間重力を直接計算することで分解能を上げる方法もあるが、この場合には空間構造が発達し高密度領域ができると計算時間が急激に増大する。もうひとつはツリー法である。この場合には、分割を適応的に行うことで空間構造が発達した場合でも大きな計算時間の増大はないが、そのかわりそもそも一様に近い場合でも計算量が多い。

　本論文で提案される方法は、ツリー法のような再帰的な空間分割を行いながらもポアソン方程式は格子上で解くことで高速な計算コードを実現したものである。

　本論文は4章からなっている。第1章では、FFT、ツリー法などの従来の計算方法が概観され、本研究のアプローチが述べられる。

　第2章では計算スキームの詳細が述べられる。適応格子を使う計算法自体は本論文が初めてというわけではないが、従来の方法のほとんどは基本格子の中の高密度領域と同定された領域を覆う直方体領域を再分割するものであった。本論文では、基本的には各メッシュ毎に他とは独立に分割する(あるいは分割をやめて上のレベルに戻す)かどうかを決定する。分割の基準は単に1メッシュ内の粒子数が設定した限界値を超えるかどうかである。このために、特別な手続きを必要としないで自動的に空間構造の変化に対応して動的に格子構造が変化する。ポアソン方程式はFFTではなくマルチグリッドを使って解く。メッシュサイズが変わる境界では、外側のメッシュで値を内挿することで小さいメッシュのために必要な境界条件を与える。さらに、メッシュサイズ毎に違う時間刻みを採用する階層的時間刻みについても新しい方法が考案されている。

　第3章では、ベクトル化および並列化の方法と、テスト計算による速度・精度評価についてまとめられている。本論文では、再帰的な複雑なデータ構造を採用しているにもかかわらず、計算量の多い各部で巧妙なアルゴリズムの書き換えを行うことで極めて高いベクトル化効率が実現されている。また、並列化については分割レベル毎に別の空間分割を行うことで粒子系のシミュレーションで並列効率を制限する主な要因となることが多いロードバランスの不均衡を回避できている。このような新しいアプローチをとった結果、実現された計算速度は非常に優れたものである。VPP5000単一プロセッサでの速度が35万粒子/秒であり、文献にあるなかでこれを上回るのは3000プロセッサを超えるASCI Redを使ったツリー法の計算のみである。計算機の速度をLinpackベンチマークでの速度を使って規格化した比較では、ツリー法および他の適応格子の実装に比べて10倍ないし100倍高速という結果が得られており、これは非常に高く評価できるものである。並列化効率についても利用できたCPU台数の範囲で満足できる結果がえられている。

　最後に、サンタバーバラクラスター比較プロジェクトにのっとった計算をすることで他のコードの計算結果との比較を行っている。時間刻み一定の場合についてはよい一致が得られており、コードの正当性は確認されたと判断できる。

　第4章は結論であり、主に今後の展望がまとめられている。

　以上要約するに、本論文は宇宙論的N体計算のための適応格子を使った新しい計算アルゴリズムをベクトル並列計算機上に実装し、従来の計算コードに比べて画期的な高速化を実現したものである。今後この計算コードを使うことで従来は不可能であった大粒子数のシミュレーションが可能になり、大規模構造の形成から銀河形成までを同時に解くことで観測との精密な比較が可能になると期待できる。本論文の一部は吉井譲および森正夫との共同研究であるが、論文提出者が主体となってアルゴリズム開発及び実装・検証を行ったものであり、論文提出者の寄与が十分であると判断できる。したがって、委員会は全員一致で本論文提出者に博士(理学)の学位を授与できると認める。