プログラムの高速化・並列化サービスの事例

＊ここに掲載するのは、クイーンズ大学ベルファストのJonathan S. Parker博士およびnAGのEdward S Smyth によるHECToRレポート「Implementation of established algorithms to extend HELIUM, Jonathan S Parker and Edward S Smyth, April 4, 2011」を要約したものです。

概要

レーザー励起されたヘリウム原子は、物質と光の量子力学的相関と量子力学的エネルギー交換の理論と実験両方に対する理想的な系です。非古典的相関（エンタングルメントとも呼びます）は、量子力学的宇宙における最も不可思議かつ非直観的な性質です。例えばエンタングルしたキュービット（entangled qubit）は量子計算の定式化における核心であり、相関した光子は量子暗号を可能にするものです。レーザー励起されたヘリウム原子の場合、高強度レーザー場により原子から放出された電子のペア（このプロセスはNSDI（非シーケンシャル二価イオン化：non-sequential double-electron ionization）と呼ばれます）の相関に着目します。NSDIは複雑に結合した原子波動関数に深く埋め込まれた捕えがたい効果ですが、相関した量子力学系のイオン化の重要な動的性質を詳細に明らかにします[1,2,3]。

全ての多電子原子の中で最も単純な原子であるヘリウムは、その3次元空間のシュレディンガー方程式が厳密解を持つ唯一の多電子原子です。レーザー励起された水素原子は同様に数値的厳密解を持ちますが、水素原子は単一電子系のため多電子相関を表すことはできません。高強度レーザー照射により励起された1電子原子の物理はこれまでによく理解されて来ました。ヘリウム原子についてはこうした研究から極めて異なる様相を持つことが明らかになりました。筆者らが得たヘリウムのシュレディンガー方程式の高度に整合的な解は全て、予測不能で驚く様な、説明が困難な結果を明らかにしました。

HELIUM[4]は、2電子原子あるいは強線形偏光レーザー場により励起されたイオンに対する非相対論的時間依存シュレディンガー方程式である、時間依存の5次元偏微分方程式を解くコードです。HELIUMの時間伝搬関数は、通常とは異なる小さな積分打切り誤差に対する要求に合わせて設計されました。この要求は、種々ありますが、NSDIのモデル化における困難さに由来するものです。例えば典型的な積分では、全NSDI効率は10^-8から10^-12までのオーダーになる可能性があり、1に規格化された波動関数では局所的な打切り誤差が10^-12より小さくなります。この要請に合わせるのは、低次数の積分スキームでは不可能です。このためHELIUMでは任意次数のArnoldi伝搬関数[4]を用います。Arnoldi伝搬関数は、この打切り誤差に対する要請がなくとも、高次数極限において積分効率の改善を示す場合があります。

このレポートは、NSDIにおける相関ペア生成中のレーザーと原子間のエネルギー交換をモデル化して解析する能力改善に向けたソフトウェア開発を議論します。この開発には、一般化されたHELIUMと、HELIUMによって生成された波動関数の新しい解析方法が含まれます。

目標と成果

新しいコードの実装は、既存のHELIUMを一般化されたレーザー分極を用いるように拡張します。新しいコードで用いられる一般化された行列要素に対して自己テストモジュールを開発します。この設計では、線形偏光場を用いる際にHELIUMと同等に振る舞う様に、新しいコードに新規機能を実行しないようにすることが可能です。新旧両方のHELIUMをテストしたところ、ｚ方向かｘ方向のどちらかに偏光した線形偏光レーザー場に対して、同等であることが実証されました。

このプロジェクトの驚くべき（かつ喜ばしい）結果の一つが、プロセッサ$コア$数を253から8001まで変化させた場合のHECToR上の性能が線形スケール以上に向上したことです。言い換えれば、253コアでの性能に比べ8001コアでの性能が8001/253倍に高速化すると予想することは可能ですが、このプログラムは4倍高速化しています。この理由は、253コアではローカル配列がL2キャッシュに収めるには大きすぎ、メモリー性能を阻害するためです。対照的に、8001コアではローカル配列は最も高速なL1キャッシュに十分収まります。

次に、OpenMP/MPIハイブリッド並列を実装して、利用対象ノード内の全コアを利用可能にしますが、コア当たりのメモリーは減少します。同時に、先述のレーザ分極実装のためにさらに大きな角運動量基底が必要になります。

次に、HELIUMの球面座標の出力を円柱座標へ変換する新しいポスト処理コードを実装しました。これはHELIUM実行と同数のプロセッサ$コア$数で実行します。ここで、HECToRでは16,000までのコア数を用いることができます。プログラムは、１および2つのcoupled spherical harmonics、および基底状態（6つのcoupled spherical harmonics）を、既知の円柱座標系での解析表現に対して検証されました。

次に、配位空間から運動量空間への終状態波動関数の変換に対して、マルチステージ・アルゴリズムを実装しました。それぞれ異なる優位点を持つ3種の方法を開発しました。新旧4種の手法を様々なコア数で検証し結果の一致を確認しました。新しい手法は既存の手法よりも極めて高速です。

運動量空間への変換コードの精度は、クーロン関数を既存のベッセル関数基底へ組み込むことで向上しました。任意次数のテイラー展開伝搬関数に基づく、新しい数値積分クーロン関数計算プログラムを開発し、正確に動作することを検証しました。この検証には、規格直交性のチェック、漸近的に正確な式との比較、およびこれとは別に独立に開発されたコード（パブリックドメインの低次Numerov法による積分コード）との比較が含まれます。

多くのテストを実行して、終状態の分布の厳密な定量解析が重要であることが判明し、このソフトウェアを用いて、理論・実験研究者の両方にとってチャレンジングな問題に適用しました。この結果驚くべきことに、3種の大きく異なる独立な方法による結果がほぼ同一であることが判りました。既存の手法はこの問題を解くことができませんでした
HELIUMによるHECToRを利用した科学研究は重要かつユニークなものです。ここでの開発はHELIUMの用途を大きく拡張します。HELIUMは実行時間が改善された現在、HECToR上のヘビーユーザとなっています。

作業項目１：HELIUMに対して交差レーザー場における計算を可能にする。

既存のHELIUMではレーザー光は線形偏光していると仮定されており、実効的にダイナミクスの自由度を減らし、問題は5次元の時間依存偏微分方程式に帰着します。ここでは最初のレーザーに直交する第2のレーザーによる分極計算を実装しました。これは既存のHELIUMで利用されている回転対称性を用いません。こうしてプログラムは、レーザー励起されたヘリウムに対する一般的な6次元の時間依存シュレディンガー方程式へ戻されることになります。

HELIUMでは、2電子波動関数の角依存性を各電子の角運動量量子数l₁,l₂、全角運動量量子数Lおよびそのｚ成分量子数Mで識別されるcoupled spherical harmonics基底関数で扱います。

時間依存性は動径関数のみに課され、球面調和関数の直交性から問題は4つの角運動量量子数の非対角項を含むハミルトニアンによる、動径成分が結合した時間依存2次元方程式に帰着します。

HELIUMは、電気双極子近似を用いて単一線形偏光レーザー場に対するこの方程式を解きます。HELIUMは動径空間（r₁,r₂）の特定の領域に対してプロセッサを割り当てます。これは動径変数に有限要素法を用いて解くため、最近接プロセッサ間のみで通信が生じるという利点があります。各プロセッサは行列要素に含まれる角積分に必要な情報を全て持ち、その際の通信は一切生じません。

当初は、相互作用ハミルトニアンに速度ゲージの利用を検討しましたが、長さゲージを用いた実験からこの方が誤りがなく正確であることが判明しました。計算効率の改善は、交差レーザー場問題においてきわめて重要です。なぜならより高次元空間状の積分では複雑性が増加するためです。特定のレーザー強度においては計算の収束のために、部分波基底関数サイズを一桁以上のオーダーで増加させなければなりません。こうした理由から長さゲージを実装しました。

交差レーザー場においては相互作用ハミルトニアンに2方向の交差項が生じるため、問題は4つの角運動量量子数ペア全てを含む6次元シュレディンガー方程式を解くことになります。既存のHELIUMは単一線形偏光に制限されていたため、角運動量量子数Mは固定（実際には定数０とする）でした。
ここで、最初のレーザー偏光方向をｚ軸、2番目のレーザー偏光方向をｘ軸方向とします。

·HELIUM X FIELDS:X交差場計算HELIUMコードのテスト

この新しい交差場コードを以降ではHELIUM X FIELDSと呼びます。開発したコードを、シングルプロセッサのワークステーションと、HECToRの253から8001コア上で検証しました。テストプログラムではコア数に依存せずに同一の結果を生成することを確認しました。また、既存のHELIUMコードとHELIUM Z FIELDSが同一の結果を生成することを確認しました。

HELIUMおよびHELIUM X FIELDSは自己テストモジュールを搭載しています。これは実行開始毎に部分波が以下の望ましい制約を満たしているかどうかを検証します：状態遷移ルールが正しいか、ハミルトニアンのパラメータ配列が正しく初期化されているか、および行列要素が正しく計算されているか。行列要素の角運動量部分は3-jシンボルと6-jシンボルの積で計算可能です[5]。このテストルーチンは複数の手法でこの結果を比較し、プログラミングエラーや精度を検証します。

HELIUM X FIELDSの性能を、コア数の関数として秒あたりに計算可能な時間ステップを計測しました。コンパイラはpgf90（-fastsse -O3）を用いました。8001コアでの計算時間は4000コアにおける時間の半分以下になりました。つまり線形スケール以上に向上しました。

線形スケール以上の改善（スーパースケールとも呼ばれます）は、並列処理で生じる可能性がある性質（利点）です。これはメッセージ通信オーバーヘッドを無視できる場合によく見られる現象です。この理由は、コアに割り当てられる配列サイズがコア数が増えると減少して、十分に大きいコア数ではデータ配列がCPUの最速のキャッシュ内に収まるためです。逆に計算に必要なデータ配列がキャッシュに収まらない場合は、メモリーからのデータフェッチに時間を取られ、FLOPS値は大きく減少します。8001コアでは、コア上のグリッドサイズは46x46であり、64KバイトのL1キャッシュと比べて十分小さいです。253コアではグリッドサイズは262x262であり、512KバイトのL2キャッシュサイズを超えています。HELIUM X FIELDSはメモリーアクセス速度に敏感で、スーパースケール性は望ましい結果です。

253コアではステップあたりの実行時間は88.2秒です。もし線形スケール性を仮定すれば、8001コアでは88.2*$253/8001$=2.79秒になります。実際には8001コア上で0.72秒です。これほどの劇的な影響は、極めて大きな通信バンド幅と同時に、コア上のグリッドをきわめて小さなくするほどの膨大なコア数によるものであり、多く見られるケースではありません。

作業項目２：HELIUMへのOpenMP/MPIハイブリッド並列の実装

最速のHELIUMバージョンはCray T3DおよびT3E上で1990秒の実行時間を記録しており、coupled spherical harmonics基底関数で並列化されています。最新バージョンでは2D動径方向グリッド上の並列化へ変更されています。動径グリッド上の並列化の利点は、より複雑な基底関数間の通信を用いず、全ての通信が最近接ハロ交換と集団総和演算のいずれかであることです。またこれは負荷バランスも改善します。基底関数間で負荷が大きく変化する場合に比べ、動径グリッド間では負荷は一定です。動径並列HELIUMはこれまで長く利用されてきており、7万コアまで極めて良好なスケール性を示します。

この作業の目的は、現状の2D動径グリッド上のMPI並列に加えて、基底関数上の並列化を実装することです。具体的な開発動機は以下の通りです：
１：より大きなコア数でスケールすることを可能にする。
２：同数コアの利用効率向上にハイブリッド並列が可能かどうかを調査する。
３：より重要な点は、HECToRのコア当たりのメモリーがフェーズ1で3GBであったのに対し、フェーズ2bでは1.33GB、将来のフェーズ3では1GBまで減少することです。この間にコードの改善を行った結果、基底関数サイズは大幅に増え、コア当たりの必要なメモリーは増加しました。コア当たりの動径風呂多くサイズを減らすことでコア当たり必要なメモリーサイズを減らすことができますが、実際上は制約があり9x9が最少限界サイズです。角運動量基底の並列化を用いずに計算はできますが、HECToRでの実行に際して、将来メモリー不足で幾つかのコアをアイドルさせなければならなくなる、つまり計算リソース浪費のリスクがあります。

基底関数並列化に際しては、MPIのような分散メモリー並列を用いず共有メモリー並列（特にOpenMP）を採用することで、通信バッファの更なるメモリー要求を回避しコード修正量を最小化します。現状のHECToRシステムは、24コアノードを搭載しGeminiインターコネクトを持つXE6システムです。24コアノードは2つの12コアチップから成ります。これは開発者視点からは、4つの6コアダイから成り、OpenMPを用いる場合全メモリーは性能を考慮すべきNUMAで共有されていると考えるほうが好ましいです。よってNUMAによる劣化を避けるためには、OpenMPスレッドは最大6までに制限し、MPIタスク内の全スレッドが同一ダイ内に配置されるようにaprunのオプションを指定するべきです。

·コードの変更

殆どの作業は、多くの実際の実行で用いられるArnoldi伝搬関数に対して行いました。並列化の追加は、作業が容易で性能に本質的ではない全体の加速と補正計算に対して行いました。修正点は以下の通りです：
1.MPIタスク間の全ての通信（ハロ交換とグローバル総和）はマスターOpenMPスレッド内で実行します。
2.Arnoldi伝搬関数には計算全てを含む並列領域が一つあります。サブルーチン呼び出し部分については、OpenMPスレッドで共有するようにローカル配列を伝搬関数の最初に宣言し、これをサブルーチンのローカル配列に用いたSMPバージョンに置き換えました。
3.基底関数内の全ての状態上の並列化を、OpenMPによるDOループ並列指示行を用いて直截的に並列化しました。スレッドの負荷バランスは崩れますが比較的小規模スケールのため、静的なループスケジューリングを用いても大きな影響はありません。
4.クリロフサブ空間内での固有値問題の計算は、スレッド毎の結果の共有を回避するために複製されました。全波動関数のスレッド分の更新に用いるために、各コアのキャッシュにローカルに結果を格納します。

·大規模基底関数

HECToR上で計算可能な基底関数サイズは、動径ブロックサイズだけで決まるのではなく、伝搬関数とそのオプションにも依存します。6次のArnoldi伝搬関数において動径ブロック数を20にして全210MPIタスクを用いる場合、ノードあたり4MPIタスク、タスク当たり6スレッドとすれば、タスク当たり16x16動径グリッドでの最大状態数は42925になります。

·同コア数上の効率

コア数が同じ場合にハイブリッド並列バージョンの性能を調査しました。HELIUMのオプションで正確に同数のコアを用いるようにするのは困難なため、簡単にノードあたり16コアとして、ダイ当たり4MPIタスクか4OpenMPスレッドのどちらかを用いました。用いたコンパイラはPGIとCrayです。

一般にハイブリッドバージョンはMPIのみのバージョンほど効率的とは言えません。MPIタスク当たり大きな動径グリッドが割り当てられている場合、キャッシュ効率の劣化がMPI通信の減少による改善を上回っています。しかしながらタスク当たり小さなグリッドの場合は、ハイブリッドの性能がより効率的な場合があります。

ハイブリッドバージョンの強スケール性を、動径グリッド数を63として2016MPIタスクを用いて測定しました。結果としてスケーラビリティは、大きなブロックサイズと大きな基底関数の場合のほうが良好でした。

作業項目３：球面座標から円柱座標へのHELIUM出力の並列変換コードの実装

HELIUMを用いたシュレディンガー方程式の数値的な解は、時間依存の6次元偏微分方程式、つまり7次元の複素数配列を表します。生成された生のデータそのままでは、有用な情報を得難い状態です。一般に、有用な物理情報を得る際に最初に行うのは、データの次元数を減らすことです。例えば、波動関数の絶対値の2乗を領域空間で全ての空間変数に関して積分をすれば、その領域に電子を見出す確率が得られます。この量は時間に依存して変化し、電子がある領域から別の領域へ移る動きを知るのに利用されます。イオン化に際しては、電子は原子のコア領域から積分領域の外へ移動します。この量の時間的な変化率からイオン化率を計算できます。

レーザー励起された2電子原子の物理量の発展を見るための最も良い方法は、2つの動径変数r_iの絶対値を変数にした確率密度関数P(r₁,r₂,t)です。HELIUMではこれは、coupled spherical harmonics基底で表現された波動関数の4つの角度成分を積分して得られます。

時間依存確率密度は1価イオン化（Pは何れかの動径変数のみが大きい場合に増加します）と2価イオン化（Pは両方の動径変数が共に大きい場合に増加します）の区別に極めて便利なものです。しかしながら多くの実験では、レーザー偏光に平行なイオン化電子の運動量成分を測定するため、この計算は上手いやり方ではありません。
このため、ｚを変数とした確率密度P(z₁,z₂)を計算する並列コード$CYLINDRICAL$を開発しました。これは、HELIUMでは線形偏光方向をｚ方向に採っているためです。

この確率密度関数の時間発展から、2価イオン化の場合に2つの電子が同方向に出ていくのか、反対方向に出ていくのかがわかります。さらに波束の追跡も可能です。しかしながらこの計算には大きなコストが掛ります。4次元の数値的な積分は、数万コアに渡ってテラバイト以上の大きさの配列で行われます。

作業項目４：新しい運動量空間終状態コードの実装

配位空間でのHELIUMの終状態を運動量空間へ変換する新たなソフトウェア実装を行います。この時点で既存コードは、大規模データと大規模コア数において、性能問題とMPI_gatherの実行が失敗するという問題を抱えていました。

数値的な積分で配位空間の波動関数ψ$r1,r2$が生成され、運動量空間の波動関数ψ$k1,k2$へ変換は（クーロンポテンシャルが無視できる範囲で）、電子の最終エネルギー分布の計算をする実験との比較のための簡便かつ正確な方法です。この変換は2次元ベッセル変換で行います。直近に行った一連の計算（HECToR上の8000から16000コアを用いた）では、計算のサイズと精度上の問題があったのは、数値積分ではなく最終ステージのベッセル変換でした。
問題が起きたのは、数値積分が局所的な空間領域$r1,r2$で演算する際に運動量空間全域が関与し、運動量積分が局所的な運動量状態|l1,l2,L,M>での演算に全空間が関与するためです。数値積分に引き続き、波動関数は全マシン上へ再配置されます。言い換えれば、波動関数を適切に組み上げるためには、各コア上にあるエネルギー波動関数の一部分は、マシン上の別の場所にある他のコアへ転送する必要があります。これをシリアル実行すると計算速度は何桁も遅くなります。8000コアで並列実行した場合はある程度許容できる時間ですが（最大で数時間）、MPIのメモリー制限により波動関数サイズは制約を受けます。

この問題に対して幾つかの改善に成功しました。HECToRフェーズ2b上において新しいコードは4倍から20倍に高速化しました。これは新しいハードウェアへ更新されたこともありますが、MPI集団通信の改善とアルゴリズムの改良によるものです。最初の提案である問題全体の並列化は、MPI alltoall集団通信のチューニングにより達成されました。

この作業項目4は、（終状態の運動量空間への返還に用いる）ベッセル関数基底よりも（クーロン力のみのハミルトニアンの固有状態である）クーロン動径関数基底を用いて、終状態に対する新たな変換を実装する事でもあります。この二つの基底関数系はクーロン力をゼロにした極限で同一になります。

クーロン関数は、任意次数のテイラー級数伝搬関数を用いて1電子シュレディンガー方程式を数値的に解くことで生成されます。クーロン関数はベッセル関数の場合と同じくらいに素早く生成可能になりました。基底関数の生成およびその終状態への射影を含む変換全体は、典型的な中規模から大規模サイズ問題において4分以内で完了します。

ここで1000コアを用いて波動関数を再配置して組み立てる問題を考えます。まず第一のベンチマークデータのパラメータを説明します。

各コアは506個の局所状態部分波を持ち、それぞれはブロックと呼ばれる32x32の複素数配列です。8001コア上では、総量66.3GBになります。このケースでは、HELIUMが実行時に持つメモリーはこの値の18倍になりますが、変換処理の第1ステージにおいて66.3GBのみがコア間で転送されます。

問題は、8001コアの全てに渡ってこのデータを再編成し、その後、エネルギー計算のためにベッセルあるいはクーロン基底でそれを表現して変換することです。これは3段階で行われます。初期のHECToR上では、明らかにメモリー不足とMPIライブラリの未成熟を要因として、計算が失敗していました。

第1ステージは8001コア上でデータを再編成します。既存の再編成手法では、各コアはその506ブロックのうちの一つを1番コアへ送信し、このコアはこれら8001ブロックを配列へ格納してディスクに書き込みます。これは506個の部分は全てに渡って繰り返されます。
第2ステージでは、シングルコアプログラムがこれら3次元配列の各々をディスクから読み込み、そのデータを2次元部分波へ再編成してディスクへこれらを書き込みます。
第3ステージでは、部分波は配位空間から、レーザーパルスで生成された終状態のエネルギー計算においてより適した表現へ変換されます。

·既存コードの改善作業

第1ステージでは新しい3手法が実装され、高速化と実行時の問題解決がなされました。

既存の手法1は簡潔な方法をとります。これはMPI集団通信の不具合回避のために開発されましたが、非常に低速でした。MPIの基盤においては、Send/Recvペアがコア間のデータ通信に用いられます。Send/Recvペアは、8001ブロックの転送においては部分波ループ内で8001回呼ばれます。コア0はこのデータをディスクに書き込み、この処理は各部分波に対して506回繰り返されます。
手法２は手法1と同じくSend/Recvループを用いますが、8001回の呼び出しは、1000回に一度集団同期バリアを呼ぶように変更されました。このやり方は劇的に実行性能を改善しました。
手法3は、上記の処理をMPI-Gather集団通信を用いて一回の処理で実行します。
手法4（Cray社のTom Edwards氏による）は、これとは別のMPI集団通信（MPI alltoall）を用い、入出力もさらに並列化します。ここで、506個の部分波の各々は一つのコアに割り当てられるため、この手法は506個のコアで実行します。

66.3GBの終状態を8000コアで計算する場合、HECToRの前モデル（pre-フェーズ2b）においてSend/Recvを用いた手法は2000秒掛かっていました（一般的な実計算では400GBを用いるのでさらに長時間になります）。

HECToR上でのMPI-alltoallによる手法4は、最も高速ですが、コア数が部分波数に等しくなければならないため、すべてのマシンで実行可能ではありません。部分波数は通常500から3000の範囲です。MPI-Gatherによる手法3は、既存の手法で用いられるコア数で動作するという利点を持ちます。

また、全ての手法の結果は同一であることが検証されました。

·クーロン動径関数への置換の実装え

クーロン動径関数の生成に話を戻します。これはベッセル関数を、束縛および非束縛状態のクーロン動径関数へ置き換える変換の実装です。クーロン関数はクーロンポテンシャル内の1電子シュレディンガー方程式の解であり、束縛状態か非束縛状態かにより異なります。一般にクーロン関数は有限差分ハミルトニアンの固有状態ですが、精度が要求される場合はHELIUMプログラムの境界条件を用いて有限差分グリッド上で計算しなくてはなりません。この要請は非束縛状態では緩和されます。近似を良くするには非束縛状態の計算は、離散的なグリッドでなく連続系として扱われ、数値積分には標準的な手法が用いられます。一方束縛状態の計算においては高い精度が本質的です。例えば、終状態波動関数から基底状態を除き、分析の対象である非束縛状態を残すことが必要ですが、これは全体の割合にして10^-7程度です。

そこで束縛状態を有限差分ハミルトニアン（巨大スパース行列）の部分固有値分解で求めます。Arnoldi-Lanczos繰り返し法を用いて、この行列から最小の固有値を持つN個の固有ベクトルを抜き出します。通常Nは50にセットされ、1電子イオン化ヘリウムのハミルトニアンの50の最低エネルギー束縛状態が生成され、ディスクに格納されます。HELIUMが終状態波動関数ψを生成した後、この50個の基底状態成分が除去されます。

既存コードは次のステップでは、波動関数ψ$r1,r2$を、積分ψ$r1,r2$exp$-ip1r1$exp$-ip2r2$dr1dr2を施し、運動量空間表現ψ$p1,p2$へ変換します。ここで、カルテシアン座標でなく角と動径変数で積分する際には、平面波でなく動径変数のベッセル関数の積分で表現されます。

新しいコードではベッセル関数の最適な置き換えとしてクーロン関数を用います。低次のNumerov法を用いたクーロン関数生成のパブリックドメインパッケージを試しましたが、信頼性に乏しく不正確でした。この代わりに、任意次数のテイラー級数積分を開発しました。新しいコードは次数を上げれば、実行性能を落とすことなく高精度の数値積分を実行できます。Numerov法を用いたシュレディンガー方程式の積分は典型的には、10^-4の局所打切り誤差を示しました。10次のテイラー積分は10^-10以下の極めて小さな局所打切り誤差を示しました。

既知の解析的標識との比較により、クーロン関数の積分における最大誤差は10^-9オーダーであることが示されました。またこのクーロン関数は直行性と期待されたノルムを持つことが示されました。また、動径距離rが大きい場所では、漸近的にクーロン関数とベッセル関数は一致します。

新しいコードの終状態分布を定量的に解析するために、理論研究者と実験研究者にとって課題となっている問題として、ヘリウムの2光子二価イオン化における非シーケンシャルな二価イオン化断面積計算に適用しました。ここでは、He⁺の束縛状態において1電子が単独でイオン化する、シーケンシャルなイオン化は無視します。このとき残る束縛電子はそのわずかな後に相関していない自由電子ペアとなってイオン化します。ここで関心があるのは、相関を持つ2つの電子がほぼ同時にヘリウム原子から自由になり放出される、非シーケンシャルなケースです。

シーケンシャルにイオン化する電子ペア生成を最小にするために、比較的弱い強度で計算を実行しました。パルスは18周期をかけて緩やかに上昇させ、30周期間一定にさせ、その後18周期かけて再びゼロにします。その後の積分では、レーザー場がない30周期において、二価イオン化した電子ペアは残されたイオンの強いクーロン場を離れ、互いに分離します。その後の解析では、電子は相互作用しないと仮定されて、各電子は1電子運動系として扱われます。この仮定は、さらに場がない30周期に電子を運動させて、相互作用が十分消滅していることを検証されます。再度物理量を計算することで相互作用の影響を評価することが可能です。

我々の興味は、レーザー原子間相互作用により非シーケンシャルに相関した電子ペアが電離するイベントが生じる確率です。このプロセスは2光子過程で、その確率は強度の2乗に比例します。そのイオン化断面積は、ここでは摂動領域にあるため、その確率を強さの2乗で割ったものに比例します。つまり断面積は強度に独立なため有用な物理量となります。

断面積の計算では束縛状態との分離が困難な問題となります。角運動量について積分した2電子の動径方向運づ両成分についてグラフを描き、その中から非シーケンシャルな2電子励起領域を選ばなくてはなりません。

ここで開発した複数の手法は良好な定量的正確性を持ち、各々は2パーセント以下の誤差で結果が一致します。

第1の手法は、上述のように運動量空間で実行されますが、前処理ステージでは主量子数9までのHe⁺の束縛状態を削除します。この束縛状態は水素原子様のHe⁺を有限差分グリッド上で解いた固有状態で、HELIUMにより同じ有限差分グリッド上で積分を行う際に用いたものと同じ境界条件を用いて計算されます。この束縛状態はあらかじめArnoldi固有値分解で計算されて、外部ファイルに格納されます。

第2の手法は、ヘリウムの終状態を、2つの電子が相互作用しない極限におけるヘリウムのハミルトニアンの固有状態の線形結合で分解するもので、さらに満足できる結果を示します。この新しい固有状態基底（これをここではクーロン状態と呼びます）は、運動量空間の終状態を構成する基底関数と微妙に異なります。クーロン状態はクーロン力が無視できるような遠方領域では運動量空間での状態と一致します。この遠方における基底関数間の同一性を利用しますが、不幸にして、束縛状態と非束縛状態が混合する領域はクーロン力を無視できない領域で生じます。クーロン状態を用いればこの問題を解決できます。クーロン基底では束縛状態と非束縛状態を区別することが可能です。結果としてクーロン基底では、終状態から束縛状態を分離する必要がなくなります。

謝辞

このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学$CSE$サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。

文献

[1]	Corkum P, 1993, Phys Rev Lett, 71, 1994
[2]	Parker JS, Doherty BJS, Taylor KT, Schultz KD, Blaga CI and DiMauro LF, 2006, Phys Rev Lett, 96, 133001
[3]	Moore LR, Parker JS, Meharg KJ, Armstrong GSJ and Taylor KT, 2008, J. Mod. Opt. 55, 2541
[4]	Smyth ES, Parker JS and Taylor KT, 1998, Comput Phys Comm, 144, 1
[5]	Rotenburg M, Bivens R, Metropolis N and Wooten J K, 1959, The 3-j and 6-j Symbols, Crosby Lockwood and Sons, Ltd. London.
[6]	http://www.hector.ac.uk/cse/documentation/Phase2b/
[7]	Baltuška A , Udem Th, Uiberacker M, Hentschel M, Goulielmakis E, Gohle Ch, Holzwarth R, Yakovlev VS, Scrinzi A, Hänsch TW and Krausz F, 2003, Nature, 421, 611
[8]	Kienberger R , Goulielmakis E, Uiberacker M, Baltuška A, Yakovlev V, Bammer F, Scrinzi A, Westerwalbesloh Th, Kleineberg U, Heinzmann U, Drescher M and Krausz F, 2004, Nature, 427, 817
[9]	Baker S , Robinson JS, Haworth CA, Teng H, Smith RA, Chirilă CC, Lein M, Tisch JWG and Marangos JP, 2006, Science, 312, 424
[10]	Uiberacker M, Uphues Th, Schultze M, Verhoef AJ, Yakovlev V, Kling MF, Rauschenberger J, Kabachnik NM, Schräder H, Lezius M, Kompa KL, Muller HG, Vrakking MJJ, Hendel S, Kleineberg U, Heinzmann U, Drescher M and Krausz F, 2007, Nature, 446, 627
[11]	Niikura H, Légaré F, Hasbani R, Ivanov MY, Villeneuve DM and Corkum PB, 2003, Nature, 421, 826
[12]	Kanai T, Minemoto S and Sakai H, 2005, Nature, 435, 470