第 5 章　粒子積分と幾何イベント

5.1 同時刻の位置と速度を使う Boris 更新

入力状態は $(\bm{x}^n,\bm{v}^n)$ であり、速度だけ半時刻にずらして保持する状態ではない。まず予測中点

\begin{equation*} \bm{x}_{\mathrm{mid}}=\bm{x}^n+\frac{\delta t}{2}\bm{v}^n, \qquad \bm{E}_{\mathrm{mid}}=\bm{E}(\bm{x}_{\mathrm{mid}})\tag{5.1} \end{equation*}

で電場を一回評価する。周期軸では、場の評価点だけを基本セルへ写す。

速度更新は電場半 step、磁場回転、電場半 step である。

\begin{align*} \bm{v}^-&=\bm{v}^n+\frac qm\bm{E}_{\mathrm{mid}}\frac{\delta t}{2},\tag{5.2}\\ \bm{t}&=\frac qm\bm{B}_0\frac{\delta t}{2},\qquad \bm{s}=\frac{2\bm{t}}{1+\bm{t}\cdot\bm{t}},\tag{5.3}\\ \bm{v}'&=\bm{v}^-+\bm{v}^-\times\bm{t},\tag{5.4}\\ \bm{v}^+&=\bm{v}^-+\bm{v}'\times\bm{s},\tag{5.5}\\ \bm{v}^{n+1}&=\bm{v}^++\frac qm\bm{E}_{\mathrm{mid}}\frac{\delta t}{2},\tag{5.6}\\ \bm{x}^{n+1}&=\bm{x}^n+\frac{\delta t}{2} (\bm{v}^n+\bm{v}^{n+1}).\tag{5.7} \end{align*}

磁場回転は Cayley 型の直交変換となり、電場が零なら速度の大きさを保つ。ただし数値的な旋回角は $2\arctan(\omega_c\delta t/2)$ であり、速さを保存しても大きな step の旋回位相まで正確とは限らない。

標準の Boris 法は位相空間体積を保存する[13]が、本更新は電場の標本点と位置の更新が異なるので、その性質をそのまま引き継ぐとは限らない。滑らかな静電場では位置・速度とも 2 次精度である。ただし衝突、境界作用、バッチ間の場の更新を含む全体の精度は、別に確認する必要がある。

$\delta t$ は、一 step の旋回角 $|q|\lVert\bm{B}_0\rVert\delta t/m$ 、電場が空間的に変化する長さを粒子が通過する時間、三角形や箱の面の近くを横切る時間のいずれも解像する大きさに選ぶ。三角形との交差は一 step の chord で判定するので、一 step 中に軌道が強く曲がる場合や細かい形状を通過する場合は、 $\delta t$ を小さくして衝突位置、命中要素、吸収と escape の統計を収束させる。 $\delta t$ を半分にするときは $N_{\max}$ を 2 倍にし、一粒子を追跡できる時間 $N_{\max}\delta t$ を揃えて比較する。

5.2 三角形との最初の交点

候補軌道の chord を

\begin{equation*} \bm{x}(\lambda)=\bm{x}^n+\lambda(\bm{x}^{n+1}-\bm{x}^n), \qquad 0\le\lambda\le1\tag{5.8} \end{equation*}

とする。三角形の頂点 $\bm{a},\bm{b},\bm{c}$ に対して

\begin{equation*} \bm{a}+u(\bm{b}-\bm{a})+v(\bm{c}-\bm{a})=\bm{x}(\lambda)\tag{5.9} \end{equation*}

を Möller–Trumbore 法で解き[9]、 $u\ge0$ 、 $v\ge0$ 、 $u+v\le1$ 、 $0\le\lambda\le1$ を満たす交点の最小 $\lambda$ を採用する。法線の符号による back-face 除去は行わず、両側から衝突を検出する。

候補の絞り込みには三角形の AABB を用いる。小規模な mesh では全要素を線形に調べ、大きい mesh では一様格子と三次元 DDA[1]によって chord が通る cell を順にたどり、その cell に登録された三角形だけを精密に判定する。

5.3 計算箱の作用と event の発生順

mesh との交点と箱の面の通過の $\lambda$ を比較し、先に起こるものだけを処理する。丸め誤差の範囲で同時なら mesh への吸収を優先する。

箱の面の作用	生存粒子への処理
開放	escape、または外部の電位障壁で帰還か escape かを判定
鏡面反射	法線速度を反転し、残り時間を再積分
再配置反射	反射に加え、面内の位置を一様に選び直す
周期	反対面へ移し、速度を保って残り時間を再積分

再配置反射は、外部を移動して帰還する粒子の面内位置を縮約する選択肢である。一つの面での反射では面内の 2 軸の位置を選び直し、辺や角で複数の面に同時に達したときは、達した面に含まれない軸だけを選び直す。光電子の上部境界で水平位置を randomize する Zimmerman らの設定に着想を持つ[22]。外部軌道の飛行時間や自己無撞着シースを再現するものではない。

event での速度は chord 方向を使い、その速さを離散仕事

\begin{equation*} |\bm{v}_{\mathrm{event}}|^2=|\bm{v}^n|^2+ 2\frac qm\bm{E}_{\mathrm{mid}}\cdot (\bm{x}_{\mathrm{event}}-\bm{x}^n)\tag{5.10} \end{equation*}

に整合させる。chord の $\lambda$ を残り時間の割合にも使うため、強い加速・旋回時の真の交差時刻を解いているわけではない。電場が零なら、磁場の強さによらず event の前後で速さを保つ。右辺が正にならない退化した event では、粒子を推定で分類せずに計算を停止する。外部シースを接続した上端面の通過も同じ chord の交差で判定する。図5.1に、場の評価点と、候補 chord 上で event を選ぶ手順を示す。

同時刻 Boris 更新と最初の幾何 event — **図 5.1**　同時刻 Boris 更新と幾何 event の模式図。左では予測中点で電場を標本化し、更新後の速度から候補位置を作る。右の橙色破線は候補 chord のうち最初の吸収より後の部分であり、粒子をそこまで進めない。真の曲線軌道と交差時刻の近似は $\delta t$ の収束で確認する。

5.4 周期画像と未解決軌道

周期場の画像の層数と、衝突で調べる画像の範囲は独立に決める。衝突では、chord の AABB と基本セルの mesh の AABB が重なる画像を列挙し、命中した画像を基本セルの要素へ対応づけて電荷を堆積する。場の評価点を基本セルへ写すだけでは、セル境界の手前で周期画像の表面へ衝突する event を検出できない。

追跡上限に達した粒子は未解決として集計し、吸収や escape と区別する。未解決の割合が大きい計算は、吸収電流が収束した結果とみなさない。