分子動力学

分子動力学（Molecular Dynamics, MD）は、原子・分子をニュートンの運動方程式に従って時間発展させることで物質の構造・輸送・反応を計算機上で追うシミュレーション手法である。核融合研究においては、プラズマと接する材料の損傷や水素・ヘリウムの取り込み、液体金属の界面挙動、強く結合した高密度プラズマの輸送係数など、原子スケールの物理が効く現象を調べる標準的な道具となっている。プラズマそのものを粒子として解くのではなく、固体・液体・分子といった凝縮系の側を原子分解能で扱う点が、流体モデルやプラズマ粒子シミュレーションとの基本的な違いである。

分子動力学とは・シミュレーションの組み立て方

MDの基本は単純である。系に含まれるすべての原子について、他の原子からの力を求め、運動方程式を微小時間刻みで積分して軌道を追跡する。得られた軌道の集合から、動径分布関数（原子間距離の統計）、拡散係数、熱伝導率、せん断粘性係数といった物理量を統計処理として抽出する。時間刻みは原子振動の周期より十分短くしなければならず、典型的にはフェムト秒程度である。このため1回の計算で追えるのはナノ秒から、大規模計算でもマイクロ秒程度の時間と、数万から数億原子程度の空間スケールに限られる。この「時間・空間の窓の狭さ」がMDの本質的な制約であり、プラズマ面材料の研究では、実験で問題になるフリュエンス（単位面積あたりの入射粒子数）に到達するまでの時間をシミュレーションで再現できないことが、結果の解釈を難しくする最大の要因となってきた。

シミュレーションセルの作り方を、プラズマ対向材料の典型例で見ると分かりやすい。タングステン結晶の表面を特定の結晶面で切り出し、その上から重水素原子を一定のエネルギー・入射角で打ち込む。原子が表面に結合したまま残るか、分子として脱離するか、内部に拡散するかは、すべて軌道から自動的に決まる。次の図は、重水素を高濃度に吸収させたタングステン(110)表面のシミュレーションセルを二つの異なるポテンシャルで計算したものである。

図では青い球がタングステン原子、桃色の球が重水素原子を表し、上半分の表層に重水素が濃く分布している様子が分かる。同じ系でも(a)と(b)で表面近くの重水素の配置が異なっており、原子間ポテンシャルの選択が結果を左右することを示している。

原子間ポテンシャル：古典・反応性・第一原理

MDの精度は、原子間の力をどう与えるかにかかっている。古典MDでは、電子の自由度をあらかじめ平均化し、原子位置だけの関数としてポテンシャルエネルギーを与える。金属や共有結合材料向けに、埋め込み原子法（EAM）やテルソフ（Tersoff）型ポテンシャルなどが使われ、数億原子の大規模計算が可能である。ただし、これらのポテンシャルは特定の結合状態を記述するよう調整されたものであり、結合の生成・切断を伴う化学反応を正しく扱えるとは限らない。反応性のある系のために開発されたのがレアックスFF（ReaxFF）に代表される反応性ポテンシャルで、結合次数を原子間距離の関数として連続的に変化させることで、分子の解離・再結合を古典MDの枠内で記述する。プラズマ中の活性酸素種と有機材料の反応など、化学的に複雑な系の解析に用いられている。

もう一つの極が第一原理MD（ab initio MD）である。原子間力を密度汎用理論などの電子状態計算から毎ステップ求めるため、ポテンシャルの調整に依存せず化学反応や電荷移動を記述できる。その代償として扱える原子数は数百から数千程度に限られる。液体リチウム中の重水素の挙動や、溶融塩ブランケット材料Li2BeF4中のトリチウムの拡散など、電荷移動や化学結合が本質的な問題では第一原理MDが用いられる。高温高密度の状態では、電子状態計算のコストがさらに増すため、軌道を明示的に解かない近似（軌道自由度を消去した密度汎関数）と量子MDを橋渡しする手法の開発も進んでいる。量子核効果が必要な場合には、核の量子性を経路積分で取り込む経路積分MDもある。要するに、系の大きさと必要な化学的精度のトレードオフの中で、古典・反応性・第一原理という三段階のポテンシャルを使い分けるのが現状である。

プラズマ対向材料：損傷・ヘリウム・水素とタングステン

核融合炉のダイバータや第一壁にはタングステンが最有力のプラズマ対向材料であり、MDの最大の応用分野の一つになっている。プラズマからの水素同位体やヘリウムの入射は、物理スパッタリング、化学スパッタリング、格子欠陥の生成、気泡の成長といった材料劣化を引き起こす。MDはこれらを原子レベルで追えるほか、スパッタされた粒子のエネルギー分布・角度分布・分子の振動回転状態といった情報も与える。この情報は重要である。というのも、表面から放出される水素分子の振動状態によって分子助起再結合の速度が桁で変わるように、放出粒子の性質はエッジプラズマの輸送やデタッチメントの予測を左右するからである。実際、水素入射による炭化水素の「迅速な化学スパッタリング」の解明はMDによるものであり、同様の機構がベリリウム表面でも働き、実験で観測されたBeD分子の放出を説明した。MDで得た分子の振動状態分布は、JETで測定されたBeDの振動温度と比較して検証されている。

ヘリウムについては、タングステン中での自己トラップと気泡成長が大きな主題である。ヘリウムはタングステン格子にほとんど固溶せず、空孔などに捕まるとさらにヘリウムを引き寄せるため、表面近くに気泡を形成する。大規模MD計算によれば、この挙動はヘリウムの入射フラックスに強く依存する。フラックスが高いと表面近くで自己トラップが起きやすく、表層に気泡の層ができる。シミュレーション時間をITER相当の低フラックスまで2桁広げ、マイクロ秒オーダーの計算を行うと、気泡が一斉に破裂する「バブルバースト」イベントが観測され、表面形態が不連続に変化する。破裂した気泡はヘリウムがより深部へ達する経路となり、深部の大きな気泡がさらに大きな表面構造を作る。同時に、この結果は従来のほとんどのMD計算がフラックスを実際より高く設定せざるを得ず、有限サイズ効果に悩まされてきたことを示した。すなわち、MDの結果を炉条件へ外挿する際には、シミュレーションが再現している条件そのものを慎重に見極める必要がある。

中性子照射に伴う損傷カスケードのMD計算では、電子・フォノン結合によるエネルギーの電子系への散逸が損傷の形態を変えることも示されている。金属中では入射粒子が原子に渡したエネルギーが一部電子系へ逃げるため、これを無視した古典MDは損傷を過大評価しがちであり、温度依存の電子・フォノン結合を組み込む改良が行われている。

液体金属の表面と界面：濡れ・すべり・トリチウム

液体リチウムをプラズマ対向材に使う「液体壁」構想では、流動する液体金属が熱を運び、蒸発したリチウムが蒸気シールドとして入射熱流を緩和する。この概念の評価には、液体リチウムと固体壁の界面での流れの境界条件の理解が不可欠である。流体力学では通常「no-slip（すべりなし）」、すなわち壁に接する液体が壁と同じ速度を持つと仮定するが、ナノスケールではこれは一つの極限にすぎず、純粋なすべりから多層のロッキング（液体が壁に固定される）まで連続的な境界条件が存在する。すべりの程度は、液体・固体間の相互作用の強さ、温度、表面粗さ、壁と液体の密度の整合性などに依存する。

図はクエット流（平行板間のせん断流）におけるすべり長さの定義と、(a) 多層ロッキングですべり長さが負になる場合、(b) no-slip、(c) 正のすべり長さを持つすべりの場合の速度分布を模式的に示したものである。MDではこのような流れを原子レベルで計算し、すべり長さを直接求めることができる。リチウムとタングステンのナノチャネルについてのMD計算では、温度、壁の格子方位、壁のひずみ、表面粗さがすべりに与える影響が調べられている。粗さのある壁では、次の図のように矩形の突起を設けたチャネルで流れを計算する。

図は下部壁に矩形の粗さ要素を持つタングステンナノチャネル内のリチウム流のスナップショットであり、赤が壁原子、青がリチウム原子を表す。この種の計算は、毛細管多孔体を用いた液体リチウムダイバータの設計に必要な、ナノスケールの流動・熱輸送データを与える。

液体リチウムの課題は水素同位体への高い反応性にある。プラズマ中の不純物を減らす効果がある一方で、トリチウムの保持が問題になる。第一原理MDによる研究では、高濃度の重水素を含む液体リチウム中で、濃度と温度の条件によっては岩塩構造のLiDが液体中で析出する液体–固体相転移が起きること、高濃度ではD2分子が形成されることが予測されている。これはトリチウム保持の評価に直接効く知見である。同様に、液体リチウムのタングステン表面での濡れ挙動、Fe–W界面での拡散、ブランケット増殖材Li2TiO3の多孔性と熱伝達など、液体金属・セラミックス界面の物性がMDで調べられている。

高密度・強結合プラズマのMD

MDは凝縮系だけでなく、プラズマそのものにも適用できる。対象はクーロン相互作用が強い、いわゆる強結合プラズマである。理想プラズマでは粒子は多数の弱い衝突を重ねて拡散するが、カップリング定数（クーロンエネルギーと熱エネルギーの比）が1を超えると、粒子は互いの位置を強く意識して局所構造を作り、通常の輸送理論が破綻する。ウォームデンスマター（温稠密物質）や慣性核融合の圧縮核、白色矮星の内部などがこの領域に当たる。

ここでMDは、摂動論に頼らず輸送係数を「測定」する役割を果たす。平衡MDで粒子軌道から速度自己相関や応力相関を計算し、グリーン・久保公式を通じて拡散係数、熱伝導率、せん断粘性係数を求める。この方法は、古典一成分プラズマ（OCP）、ヤワ型一成分プラズマ（遮蔽されたクーロンポテンシャルを用いるモデルで、ダストプラズマやレーザー冷却イオン系のモデルでもある）、二成分完全電離強結合プラズマなどに適用され、強結合領域での輸送係数の信頼できるデータを与えてきた。電荷状態が複数混在する系のMD計算では、ウォームデンスマター領域での自己拡散係数が評価されている。また、強結合な二イオンプラズマにおける相転移や、磁場中の衝突性プラズマ中の帯電ダスト粒子の異常輸送もMDで調べられている。

この分野特有の技術的課題は、長距離クーロン力の扱いと、量子効果の近似である。クーロン系では力の計算をどう打ち切るかが結果に影響するため、初期配置の最適化など数値的な工夫も提案されている。高温では電子の量子性が重要になるため、量子MDと、より安価な軌道フリーMDを状態方程式の計算で接続する試みも進んでいる。強結合プラズマのMDは、実験が困難な状態点での物性を与える基礎データ源として、核融合と天体物理の双方で役割を担っている。