コラム - 磁気ミラー効果のシミュレーション

カスプ型の非一様磁場中でプラズマが閉じ込められるプロセスは, 磁気ミラー領域における荷電粒子の反射や損失により説明される事が多い. 磁気ミラーの効果は, スパッタ装置などで利用されるマグネトロンプラズマの基本原理であり, 他には核融合炉や粒子加速器へも応用されている^{[参考1,2,3]}. また, 地磁気圏のオーロラ研究においても重要な効果の1つとされる^[参考4].

しかし, 従来のプロセスプラズマのシミュレーションモデルでは, 磁場だけでなく, 電場や衝突によっても速度分布が乱されるため, 磁気ミラー効果のみを可視化して精度良く評価することが困難であった. 本記事では, 磁気ミラー効果の評価を目的とした計算モデルを構築し, そのシミュレーション結果とあわせて詳細を報告する.

記載内容 (頁内リンク)

1. はじめに
2. シミュレーション設定
3. 関連用語と変数記号
4. 出力変数と評価項目
5. シミュレーション結果 - 電子1つの解析
6. シミュレーション結果 - 複数電子の解析
7. まとめ
8. 付録. 理論式の導出
9. 補注
10. 参考資料

PDF版（※動画は再生されません）

非一様磁場中の電子の運動を, 粒子法の一種である Particle-in-Cell 法（PIC法）でシミュレーションした. ツールには弊社で開発されたソフトウェアである Particle-PLUS^[参考5] を用いた. 空間は軸対称を想定してモデル化し, 2次元の計算メッシュを用いた. 磁場配置は, Particle-PLUS で使用されている基本例題（マグネトロンプラズマのチュートリアル）と同じものである. その例題のシミュレーション結果を下図に示す.

今回は, 図で示されているプラズマの部位に複数の電子を配置して, その後の運動の様子を観察した. 電子群の初速は温度1eVのマクスウェル分布^(補注a) で与えている. なお, 電場（境界条件から与えられる電場, および空間中の荷電粒子が生成する電場）や粒子間衝突による散乱は, 今回は評価の邪魔となる. それらの影響を無くすために, 空間全域の電位と衝突断面積はゼロに固定している.

ミラー効果は磁気モーメントの断熱不変性から説明される. 関連用語および変数の定義と, それらが持つ意味を本節にまとめておく. 各種理論式の導出過程は 付録 に記載する.

3.1. 変数記号

本記事内では荷電粒子の変数について次表の記号を用いる.

記号	意味
$\vec{B}$	荷電粒子が受ける磁場（局所磁場）(補注b)
$m$	質量
$\vec{v}$	速度
$\vec{v}_{/\!/}$	$\vec{v}$ のうち $\vec{B}$ に平行な成分（下図参照）
$\vec{v}_{\perp}$	$\vec{v}$ のうち $\vec{B}$ に垂直な成分（下図参照）
$\theta$	$\vec{B}$ に対する $\vec{v}$ のピッチ角（3.3節参照）
$K$	運動エネルギー
$\mu$	磁気モーメント（3.2節参照）
$R$	ミラー比（3.3節参照）

また, 本記事では, ある磁束線に沿って運動する荷電粒子の位置と, 受ける磁場が最も弱くなる位置, および最も強くなる位置, の3種類を区別して記述する場合がある. それらの位置を表す添字には ${\rm c}$ , $0$ , $1$ を使用する. 例えば, 荷電粒子が受ける磁場を表記する場合, その瞬時値を $B_{\rm c}$ , 最小値を $B_0$ , 最大値を $B_1$ と表す（付録図参照）.

3.2. 磁気モーメント

磁場中でローレンツ力を受けて旋回運動^(補注c) している荷電粒子の 磁気モーメント（Magnetic Moment）$\mu$ は次のように定義される. \begin{align} \mu := \frac{mv_\perp^2}{2B} \end{align} これは一般的な周期系における 作用変数（Action Variable）に相当し, 磁場 $B$ が旋回に対して緩やかに変化する場合は一定不変に保たれる^(補注d,e). このような量を 断熱不変量（Adiabatic Invariance）という^{[参考1,6,7]}. この不変性に着目すると, $B$ と $v_\perp$ の間には, 下図に示す関係が成り立つ.

3.3. 磁気ミラー

上図で示したような, 磁束線が漏斗状に収束している領域では $B$ が強くなる. 荷電粒子がその領域に近づくと, $\mu$ を一定に保つために $v_{\perp}$ は増大するが, 全運動エネルギー $K := \frac{1}{2}mv^2$ が一定であるため, 必然的に $v_{/\!/}$ は減少する. そして, $B$ が十分大きい場合は $v_{/\!/} = 0$ となり, その位置で旋回中心の進行方向が反転する. この現象を ミラー効果（Mirror Effect）という.

プラズマの磁気閉じ込めでは, 2つの磁気ミラーで挟んだ磁場配位がよく利用される（下図参照）. この場合, 荷電粒子の感じる磁場は中央で最も弱くなり, ミラー領域で最も強くなる. それらの磁場強度の比, \begin{align} R := \frac{B_1}{B_0} \quad (> 1) \end{align} を ミラー比（Mirror Ratio）という. これはミラー効果の強さを特徴付ける量であり, 粒子は理想的には反射を繰り返して $B/B_0 \le R$ の弱磁場領域に閉じ込められる（付録2参照）. 十分大きなミラー比が確保されない場合, 粒子は反射されずに逃げてしまう.

ここで, $\vec{v}_{/\!/}$ と $\vec{v}$ の成す角（ピッチ角）$\theta$ を, 次式のように導入する^(補注f). \begin{align} \tan \theta := \frac{v_{\perp}}{v_{/\!/}} \quad (0 \le \theta \le 90°) \end{align} このとき, 反射して閉じ込められる臨界条件の粒子に対しては, $\mu$ と $K$ の不変性から次の関係が導かれる（付録1参照）. \begin{align} R = 1/\sin ^2 \theta _0 \end{align} これは重要な結果であり, 荷電粒子の閉じ込めに必要なミラー比の大きさは, 粒子の質量や電荷に依らず, 初期位置（磁場最小点）でのピッチ角のみで決定されることを意味している.

3.4. ロスコーン

$\theta$ の小さい粒子は $v_{/\!/} = 0$ となる事ができず, ミラー領域を通過して逃げ出す（付録3参照）. ミラー領域から粒子が失われる現象を ミラー損失（Mirror Loss）という. この範囲は速度空間内では下図のように円錐状となっており ロスコーン（Loss Cone; 損失円錐）と呼ばれる. このときの荷電粒子の速度分布関数は $\theta/\theta _0 \lt 1$ の円錐部分が消失した異方性分布になり, 結果としてマクスウェル-ボルツマン分布からずれる^[参考6].

今回は扱わないが, もし衝突などによりピッチ角が変化してロスコーン内へ散乱されると, その粒子は系から失われる. 一般に, 電子はイオンに比べて高い衝突周波数を持つため, より容易に失われやすい^[参考1].

今回のシミュレーションでは, 特定の空間セル内に配置された電子の中から無作為に抽出（サンプリング）した粒子を対象として, 次の物理量の時間推移を出力した.

• 磁気モーメント $\mu$
• 運動エネルギー $K$
• 電子が受ける磁場 $\vec{B}$
• 磁場に平行な方向の速度成分 $\vec{v}_{/\!/}$
• 磁場に垂直な方向の速度成分 $\vec{v}_{\perp}$
• ピッチ角 $\theta$

ポスト処理では gnuplot^[参考8] を用いて各種グラフを描画し, 主に次の事項に着目して可視化と評価を行った.

• $\mu$ と $K$ の保存性（一定に保たれているか）
• ミラー効果による電子の閉じ込め性能（位置空間プロットによる軌跡の確認）
• ミラー反転位置における 式(6) の成立可否
• ロスコーン領域の形成（速度空間プロットによる評価）

位置空間における電子1つの軌跡を下の動画に示す. また, そのときの速度成分 $(v_{/\!/}, v_{\perp})$ およびピッチ角 $\theta$ の時間推移をその下の図に示す. 電子が $B$ の強い領域に近づくにつれて $\theta$ が大きくなり, $\theta = 90°$（$v_{/\!/} = 0\,$）となる位置で旋回中心の進行方向が反転する様子が観察された. この電子はミラー効果によって安定的に磁場内に閉じ込められている, と言える.

次に運動エネルギー $K$ と磁気モーメント $\mu$ の時間推移を下図に示す. いずれの物理量も長時間平均では一定に保たれており, 良い保存量である事が確認された.

なお, $\mu$ に関しては 1.35×10^－9 秒の短い周期で振動する成分が観測された. この振動周期は, その位置の磁場強度 $B$ から計算される電子の旋回周期に一致する. この周期よりも短い時間スケールでは, 電子の受ける $B$ が（空間的な不均一性により）急激に変化するために, それに合わせて $\mu$ も局所的に変動するものと理解できる. 断熱不変性が成立するための条件に関しては 3.2節 で触れた通りである. また, 同様の振動成分は約1060step（3.2×10^－8 秒）の周期でも観測されており, これは旋回中心の磁場方向に沿った往復振動に該当する^[参考7]. これらの振動成分は以降の議論や解析には影響を与えないため, 本記事ではこれ以上深く立ち入らない.

電子の初期位置が異なると, 受ける $\vec{B}$ が変わる. また, 旋回中心が位置する磁束線も変わる可能性があり, そうなると $R$ の値も一意に定まらくなる. これはロスコーンの正確な分析を困難にする.

そこで, 本記事では解析アプローチを2つに分けて進める. まず6.1節では, 様々な位置に配置された電子を対象として, 粒子運動論から得られる 式(6) の広範な再現性を検証する. 続く6.2節では, 磁場環境が同等とみなせるほど初期位置が近い電子のみを対象として, 詳細なロスコーン分析を試みる.

6.1. 初期位置が大きく異なる複数電子

ミラー効果で捕捉された電子サンプル10個について, 定義に従って計算されるミラー比 $R:=B_1/B_0$ と, 磁場最小位置でのピッチ角から計算される $R_0:=1/ \sin ^2 \theta _0$ の値を比較した. その結果, 両者は誤差0.5%以内の極めて高い精度で一致した. 下図の表はその計算結果の抜粋である. これらの定量的な比較により, シミュレーション上で $R=R_0$ , すなわち 式(6) の関係式が高い精度で成立していることが確認された.

6.2. 初期位置が比較的近い複数電子

本節では, 新たな電子サンプル199個を使用した. そのうち無作為に抽出した20個について, 位置空間および速度空間における軌跡を下の動画に示す. 橙色の軌跡は磁場内に捕捉された電子（8個）を表しており, 青色の軌跡は損失した電子（12個）を表している. 本シミュレーションでは Y = 1mm の位置に壁面が存在しており, そこへ到達した電子は空間から消滅する. これは損失に計上している. 動画では, 初期ピッチ角の小さい粒子が空間内で $\theta = 90°$ ($v_{/\!/} = 0\,$) の反射条件に達することができず, 壁面に衝突して消滅していく様子が観察された.

次に, 各電子のピッチ角 $\theta$ が最小となる（磁場強度が最小となる）瞬間を下図の速度空間上にプロットした. 薄い青色の点は左側のミラーで損失する電子（20個）, 濃い青色の点は右側のミラーで損失する電子（43個）, 橙色の点は最後まで捕捉される電子（136個）を表している. 本シミュレーションでは, 左右で空間の磁場強度が異なる非対称な磁場配位（左側の方が強い磁場. 2章参照）を使用しているため, ロスコーン境界は左右で個別に存在すると考えられる. 磁束線上の空間磁場強度の最小値と最大値から見積もられるロスコーン角は, 左側で36.5°, 右側で47.5°であり, 図中ではこれらのロスコーン領域の境界を点線で図示している. 点線の色は損失電子の色と同じにしてある. 実際に損失した電子は, それぞれのロスコーン領域に収まっていることが確認された.

なおロスコーン角 47.5°の境界付近では, 閉じ込めと損失の判定結果に 2°～3°程度の誤差が生じている. これは, 電子のらせん運動が有限の旋回半径を持つことに起因していると考えられ, その具体的な要因として次の2つの可能性が挙げられる.

• 旋回中心軌道のばらつき：
  同一の位置から運動を開始したとしても, 旋回の初期位相によって旋回中心の軌道（旋回中心が位置する磁束線）が電子ごとに僅かに変わる可能性がある. 磁束線が異なると, $B_0$ と $B_1$ も変わり得る.
• 境界到達時の軌道のずれ：
  らせん運動の軌道と旋回中心の軌道の違いにより, 壁面への衝突直前で受ける磁場強度が異なる可能性がある.

• シミュレーションソフト Particle-PLUS を用いて, 電場と粒子間衝突の影響を排除した系を構築することで, 磁気ミラー効果そのもののダイナミクスを詳細に解析することが可能となった.
• カスプ磁場内に捕捉された電子を対象とした検証により, 磁気モーメントの断熱不変性, および粒子運動論から導かれる 式(6) の関係性が, シミュレーション上でも高い精度で成立することが確認された.
• 損失した電子の最小ピッチ角を速度空間上にプロットすることで, ロスコーン領域の明確な可視化に成功した. シミュレーションから得られた実際のロスコーン角は, 式(6) から予測される理論値と概ね一致することが示された.

以下では各種理論式の導出を示すが, ここでは簡易的な証明に留めている. より厳密な導出過程や詳細に関しては, 教科書などの専門文献を参照されたい^[参考1,9].

なお, 式中に登場する添字 ${\rm c}$ , $0$ , $1$ の物理的な意味については 3.1節 で説明されているが, 本シミュレーション結果の図を用いて改めて示すと次のようになる.

付録1. ミラー比と最小ピッチ角との関係式(6)の導出

ある磁束線に沿ってミラー効果で反射を繰り返す荷電粒子を考える. 粒子の受ける磁場が最小・最大となる位置（添字 $0$ と $1$ ）について, $\mu$ と $K$ の不変性を立式する. \begin{align*} \frac{{v_0}^2_\perp}{B_0} &= \frac{{v_1}^2_\perp}{B_1} \\ \\ v_0^2 &= v_1^2 \\ &= {v_1}^2_{\perp} + {v_1}^2_{/\!/} \\ &= {v_1}^2_{\perp} \quad (\because \; v_{1/\!/} = 0) \end{align*} これら2つの式から $v_{1\perp}$ を変数消去すると次式が得られる. $$\frac{B_1}{B_0} = \frac{v_0^2}{{v_0}^2_\perp}$$ 左辺は $R$ の 定義式(4) に等しい. 右辺は $\sin \theta = v_\perp/v$ である事を用いて $\theta _0$ で書き換えられる. すると, $$R = \frac{1}{\sin ^2 \theta _0}$$ が得られる.

付録2. ミラー磁場による閉じ込め条件(3.3節)の導出

ある磁束線上の荷電粒子の位置（添字 ${\rm c}$ ）と, 磁場が最小となる位置（添字 $0$ ）について, $\mu$ と $K$ の不変性を立式する. \begin{align*} \frac{{v_0}^2_\perp}{B_0} &= \frac{{v_{\rm c}}^2_\perp}{B_{\rm c}} \\ \\ v_0^2 &= v_{\rm c}^2 \\ &= {v_{\rm c}}^2_{\perp} + {v_{\rm c}}^2_{/\!/} \end{align*} これら2つの式から $v_{{\rm c}\perp}$ を変数消去すると次式が得られる. $${v_{\rm c}}^2_{/\!/} = v_0^2 - \frac{B_{\rm c}}{B_0} {v_0}^2_\perp$$ ここで, ${v_{\rm c}}^2_{/\!/} \ge 0$ が実数解を持つためには, \begin{align*} v_0^2 - \frac{B_{\rm c}}{B_0} {v_0}^2_\perp & \ge 0 \\ ⇔ \quad \frac{B_{\rm c}}{B_0} & \le \frac{v_0^2}{{v_0}^2_\perp} \\ ⇔ \quad \frac{B_{\rm c}}{B_0} & \le R \end{align*} が満たされなければいけない. これはすなわち, 荷電粒子が $\mu$ を一定に保ちながら運動するためには, 受けることのできる $B_{\rm c}$ には上限が存在し, $B/B_0 > R$ の強磁場領域へは侵入できない事を意味する.

付録3. ロスコーン範囲(3.4節)の導出

ある磁束線上の荷電粒子の位置（添字 ${\rm c}$ ）と, 磁場が最大となる位置（添字 $1$ ）について, $\mu$ と $K$ の不変性を立式する. \begin{align*} \frac{{v_{\rm c}}^2_\perp}{B_{\rm c}} &= \frac{{v_1}^2_\perp}{B_1} \\ \\ v_{\rm c}^2 &= {v_1}^2_\perp \end{align*} 2つの式から $v_{1\perp}$ を変数消去すると次式が得られる. \begin{align*} \frac{v_{\rm c}^2}{{v_{\rm c}}^2_\perp} &= \frac{B_1}{B_{\rm c}} \\ &= R \frac{B_0}{B_{\rm c}} \end{align*} ここで, 添字の定義（上図参照）より $B_0 \le B_{\rm c}$ なので, $$\frac{v_{\rm c}^2}{{v_{\rm c}}^2_\perp} \le R$$ が満たされなければいけない. 両辺をピッチ角でそれぞれ表し, 角度範囲が $0 \le \theta \le 90°$ である事も考慮すると, \begin{align*} \frac{1}{\sin ^2 \theta _{\rm c}} & \le \frac{1}{\sin ^2 \theta _0} \\ ⇔ \quad \frac{\sin ^2 \theta _{\rm c}}{\sin ^2 \theta _0} & \ge 1 \\ ⇔ \quad \frac{\sin \theta _{\rm c}}{\sin \theta _0} & \ge 1 \\ ⇔ \quad \frac{\theta _{\rm c}}{\theta _0} & \ge 1 \end{align*} が得られる. これが閉じ込められる条件であるため, 逆に損失する $\theta$ の範囲は, $$ \frac{\theta}{\theta _0} \lt 1 $$ である.

1. 3次元速度ベクトルが等方性のマクスウェル-ボルツマン分布^[参考10] （Maxwell-Boltzmann Distribution^{[wikipedia↗]}）に従うとしている. これらの電子サンプル群は, シミュレーション開始時にモンテカルロ法（Monte Carlo Method^{[wikipedia↗]}）で生成している.
2. マクスウェル電磁気学における $\vec{B}$ は, 正確には磁束密度（Magnetic Flux Density^{[wikipedia↗]}）と呼ばれる量である. ただし, 本記事では磁場強度（Magnetic Field Strength）$\vec{H}$ が登場しないため, 記述の簡潔さを優先して $\vec{B}$ を単に「磁場」と呼んでいる.
3. この回転運動は, サイクロトロン運動（Cyclotron motion）, ジャイロ運動（Gyromotion）, ラーモア運動（Larmor motion）等と呼ばれる. いずれの用語も「磁場で誘起される荷電粒子の円運動」を指し, 本質的な違いは無い. 歴史的・文脈的な背景として, サイクロトロンは粒子加速器に由来する用語であり, 軌跡に着目する傾向がある. ジャイロは角運動量保存則のために回転軸が歳差運動する現象, いわゆるジャイロ効果全般を指す. ラーモアは磁気モーメントの歳差運動の発見者に因む. 読者の混乱を避けるために, 本記事ではなるべく「旋回運動」という表現に統一している.
4. 外部パラメータに依存するハミルトン系において, ハミルトニアン（Hamiltonian^{[wikipedia↗]}）の時間変化が十分に緩やかであれば, その系は断熱的であると言う. なお, 歴史的経緯から熱力学でも同じ「断熱（Adiabatic）」という語が登場するが, 異なる意味を持つため誤解を生みやすい表現ではある. ハミルトン系の断熱状態は, 現代の物理学における（Quasistatic Process）の概念に近い^[参考11].
5. ここでの「Bが緩やかに変化する」とは, 旋回周期および旋回半径に対して B の時間的・空間的に変化が十分に小さい事を意味する^[参考9].
6. 本記事ではピッチ角の定義域を 0≤θ≤ 90° としているが, 文献によっては 0≤θ≤180° で定義される場合もある. もし運動方程式に基づいて粒子の進行方向（前進・後退）まで追跡・議論する際は, 後者の方が扱いやすいであろう. 一方で, 本記事は保存則ベースの議論で完結しており, またミラー効果による粒子の捕捉・損失のダイナミクスは θ=90° を境にして対称である（3.4節 で示したロスコーンは鏡面対称である）ため, 今回は 3.3節 の定義域で十分である.

1. F. F. Chen 「プラズマ物理入門」, 内田岱二郎 訳, 丸善 (1977)
2. T. Motomura and T. Tabaru, Vacuum and Surface Science, Vol.63, No.6 (2020) 283-288
3. H. Akimune, T. Yamagishi, I. Ikeda and F. Okamoto, Kakuyugo Kenkyu, Vol.42, No.4 (1979) 472-531
4. Y. Miyoshi et al., J. Plasma Fusion Res. Vol.97, No.5 (2021) 269-275
5. 株式会社ウェーブフロント 「Particle-PLUS」, https://www.wavefront.co.jp/CAE/particle-plus/^{[リンク↗]}
6. 長倉三郎, 井口洋夫, 江沢洋, 岩村秀, 佐藤文隆, 久保亮五 「岩波 理化学辞典」, 岩波書店 (1998)
7. S. Shinji, NICT Reports, Vol.67, No.1 (2021) 105-113
8. T. Williams and C. Kelley 「gnuplot」, http://www.gnuplot.info/^{[リンク↗]}
9. 電気学会, 電気学会大学講座シリーズ 「プラズマ工学」, オーム社 (1997)
10. 株式会社ウェーブフロント, コラム 「プラズマ・流体シミュレーションの基礎」, https://www.wavefront.co.jp/CAE/particle-plus/column/plasma_simulation_basic.html^{[リンク↗]}
11. 予備校のノリで学ぶ大学の数学・物理 「断熱不変量(古典力学系)」, YouTube (2019), https://www.youtube.com/watch?v=3QBghN-QwQ8^{[リンク↗]}