電磁プラズマ粒子シミュレーション¶

著者:: 松本洋介（千葉大学）

宇宙空間における爆発的現象に伴う、非熱的粒子の生成メカニズムは宇宙物理学に残された最大の問題の一つであると言えます。流体近似を行った磁気流体（MHD）シミュレーションでは、このような高エネルギー粒子生成メカニズムを自己無同着に理解することができないため、無衝突プラズマの第一原理計算手法によるアプローチが代わりに大きな役割を果たします。

本章では、そのような計算手法のひとつである、 電磁プラズマ粒子シミュレーション のアルゴリズムについて説明します。

概要¶

変数の定義¶

q¶: 電荷
m¶: 静止質量
c¶: 光速
E¶: 電場
B¶: 磁場
単位系¶: pCANS ではCGS系を採用しています。

基礎方程式¶

基礎となる方程式系は、Vlasov方程式と

(1)¶ $\frac{\partial f_s({\bf x},{\bf v})}{\partial t} +{\bf v}\cdot\nabla f_s({\bf x},{\bf v})+\frac{q_s}{m_s}\left( {\bf E }+\frac{\bf v}{c}\times{\bf B} \right) \cdot\nabla_v f_s({\bf x},{\bf v}) = 0$

Maxwell方程式

(2)¶ $\frac{\partial \bf B}{\partial t} = - c \nabla \times {\bf E}$

(3)¶ $\frac{\partial \bf E}{\partial t} = + c\nabla \times {\bf B} - 4\pi \sum_s^{i,e} q_s \int\int\int {\bf v}f_s({\bf x},{\bf v}) d{\bf v}$

(4)¶ $\nabla \cdot {\bf B} = 0$

(5)¶ $\nabla \cdot {\bf E} = 4\pi \sum_s^{i,e} q_s \int\int\int f_s({\bf x},{\bf v}) d{\bf v}$

のVlasov-Maxwell方程式系です。式 (1) - (7) をカップルして解く無衝突プラズマの第一原理計算手法のうち、式 (1) の分布関数 $f({\bf x},{\bf v})$ の時間発展を直接解く手法を電磁ブラソフシミュレーション、分布関数 $f({\bf x},{\bf v})$ を粒子に代表させて解く手法を 電磁プラズマ粒子シミュレーション と呼びます。以下では、 電磁プラズマ粒子シミュレーション においてこれら方程式をどのようにして数値的に解いていくかを見ていきます。

粒子と電磁場の時間発展¶

分布関数の位相空間上での時間発展は、実空間上の移流（移動、式 (1) の第2項）と速度空間上の移流（加速、式 (1) の第3項）の二つに分けて考えることができます [1] 。粒子シミュレーションではこれを、個々の粒子の運動を以下の様にして解くことにより、表現します。

加速：

(6)¶ $\frac{{\bf u_p}^{t+\Delta t/2}-{\bf u_p}^{t-\Delta t/2}}{\Delta t} = \frac{q_p}{m_p}\left({\bf E}^t+\frac{{\bf u_p}^t}{c \gamma_p^t}\times{\bf B}^t\right)$

移動：

(7)¶ $\frac{{\bf x_p}^{t+\Delta t}-{\bf x_p}^t}{\Delta t} = \frac{{\bf u_p}^{t+\Delta t/2}}{\gamma_p^{t+\Delta t/2}}$

ここで、 $\bf x_p$ と $\bf u_p$ はそれぞれ粒子pの位置、4元速度を表し、各変数の右肩にある $t$ や、 $\Delta t$ は、時刻とシミュレーションにおける時間刻み幅を表します。また、 $\gamma_p=\sqrt{1+u_p^2/c^2}$ はローレンツ因子を表します。Maxwell方程式とは電流密度（式 (3) 右辺第2項）を通じてカップルします。電流密度は以下の様にして計算します。

(8)¶ ${\bf J}^{t+\Delta t/2} = \sum_{p=1}^{N} q_p \frac{{\bf u_p}^{t+\Delta t/2}}{\gamma_p^{t+\Delta t/2}}$

ここで、Nは対象となる粒子数を表します。上記電流密度を使うと結局、

(9)¶ $\frac{{\bf B}^{t+\Delta t}-{\bf B}^t}{\Delta t} = - c \nabla \times {\bf E}^{t+\Delta t/2}$

(10)¶ $\frac{{\bf E}^{t+\Delta t}-{\bf E}^t}{\Delta t} = + c\nabla \times {\bf B}^{t+\Delta t/2} - 4\pi{\bf J}^{t+\Delta t/2}$

が粒子シミュレーションで解くべきMaxwell方程式となります。

物理量の時間更新手順を図に示します。時間更新は、

粒子の速度を更新する（加速）
粒子の位置を更新する（移動）
半整数時刻における電流密度を、整数時刻の位置と半整数時刻の速度から計算する
電磁場を更新
1に戻る

のような手順で進めます。以下では、式 (4) と (7) を満たしながら、式 (6) - (10) をどのようにして解いていくかを具体的に見ていきます。

Footnotes

Buneman-Boris法¶

粒子の速度の時間更新は、式 (6) を元に行います。ここで ${\bf u}^{t+\Delta t/2}$ を求めるのに、 ${\bf u}^{t}$ の情報が必要なことに注意してください。従って、時間更新に工夫が必要となります。以下では、プラズマ粒子シミュレーションにおける標準的アルゴリズムである、 Buneman-Boris法 について説明します。Buneman-Boris法では、粒子の加速を、電場による加速（速度空間上における移流）と、磁場による加速（速度空間上における回転）に分けて取り扱います。

まず、以下のような変数を定義します。

(11)¶ ${\bf u_p}^{-} = {\bf u_p}^{t-\Delta t/2}+\frac{q_p}{m_p}\frac{\Delta t}{2}{\bf E}^t \\ {\bf u_p}^{+} = {\bf u_p}^{t+\Delta t/2}-\frac{q_p}{m_p}\frac{\Delta t}{2}{\bf E}^t$

この新しい変数を用いて式 (6) を書き表すと、

(12)¶ $\frac{{\bf u_p}^{+}-{\bf u_p}^{-}}{\Delta t} = \frac{q_p}{2 \gamma_p^n m_p c}\left ({\bf u_p}^{+}+{\bf u_p}^{-} \right) \times {\bf B}^t$

となります。ここで、式 (12) に $({\bf u_p}^{+}+{\bf u_p}^{-})$ との内積を取ると、右辺は0となるので、その結果 $|{\bf u_p}^{+}|=|{\bf u_p}^{-}|$ の関係が得られます。このことから、 Buneman-Boris法におけるベクトルの関係のように、 ${\bf u_p}^{-}$ から ${\bf u_p}^{+}$ を求めるには、 $\bf B$ を軸として、 ${\bf u_p}^{-}$ を角度 $\theta$ だけ回転させれば良いことがわかります。

_images/b-b.png — Buneman-Boris法におけるベクトルの関係¶

角度 $\theta$ は

(13)¶ $\tan{\frac{\theta}{2}} = \frac{|{\bf u_p}^{+}-{\bf u_p}^{-}|}{|{\bf u_p}^{+}+{\bf u_p}^{-}|} = \frac{q_pB_p^n \Delta t}{2m_p \gamma_p c} = \frac{\Omega_p \Delta t}{2}$

で与えられます。ここで、 $\Omega_p=q_pB_p^n/m_p \gamma_p c$ は粒子pのジャイロ周期を表します。これより角度 $\theta$ は、

(14)¶ $\theta = 2 \arctan{\frac{\Omega_p \Delta t}{2}} = \Omega_p \Delta t \left( 1 - \frac{(\Omega_p \Delta t)^2}{12} + ... \right)$

から求められ、時間2次精度でジャイロ運動を記述できることが分かります。

四則演算数を少なくするために、実際の計算では一度に角度 $\theta$ の回転を行わず、

(15)¶ ${\bf u_p}^{*} = {\bf u_p}^{-}+{\bf u_p}^{-} \times {\bf t} \\ {\bf u_p}^{+} = {\bf u_p}^{-}+{\bf u_p}^{*} \times {\bf s}$

の2段階に分けて解くことを行なっています。ここで ${\bf t}$ と ${\bf s}$ はそれぞれ、

(16)¶ ${\bf t} = \frac{{\bf B}^n}{|{\bf B}^n|} \tan{\frac{\theta}{2}} = \frac{q {\bf B}}{m_p \gamma_p c}\frac{\Delta t}{2}$

(17)¶ ${\bf s} = \frac{2 {\bf t}}{1+t^2}$

で与えられます（ Buneman-Boris法におけるベクトルの関係）。

以上をまとめると、粒子の加速は

電場による加速

(18)¶ ${\bf u_p}^{-} = {\bf u_p}^{t-\Delta t/2}+\frac{q_p}{m_p}{\bf E}\cdot\frac{\Delta t}{2}$

磁場による加速（回転）

(19)¶ ${\bf u_p}^{*} = {\bf u_p}^{-}+{\bf u_p}^{-} \times {\bf t} \\ {\bf u_p}^{+} = {\bf u_p}^{-}+{\bf u_p}^{*} \times {\bf s}$

電場による加速

(20)¶ ${\bf u_p}^{t+\Delta t/2} = {\bf u_p}^{+}+\frac{q_p}{m_p}{\bf E}\cdot\frac{\Delta t}{2}$

の3段階で行います。

最後に、粒子の位置を更新は、式 (7) を

(21)¶ ${\bf x_p}^{t+\Delta t} = {\bf x_p}^t+\frac{\bf u_p}{\gamma_p}^{t+\Delta t/2}{\Delta t}$

によって解くことにより行います。

Particle-in-Cell法¶

粒子の運動を解く上で（式 (18) - (20) ）、粒子の位置における電磁場の情報が必要となります。電磁プラズマ粒子シミュレーションにおいて、粒子はシミュレーション空間内の任意の位置に存在するのに対して、電磁場は、下記で示すように、空間内のグリッド点上に定義されています。また、式 (8) や下記で示す電磁場の時間発展を解く際には、逆にグリッド点での電流密度の情報が必要となります。粒子シミュレーションでは、ラグランジュ量である粒子とオイラー量である（電磁）場を、粒子を質点ではなく有限の広がりを持った 超粒子 として扱うことによって、結びつけます。

ここで1次元（ $x$ 方向）の場合を例にとって、具体的な計算を見て行きましょう。粒子の位置 $x_p$ における場 $F_p$ は、グリッド点 $i$ に定義されたの場 $F(X_i)$ の情報を用いて

(22)¶ $F_p(x_p) = \sum_i F(X_i) S(X_i-x_p) = \sum_i F_i S_i(x_p)$

より求めます。ここで、 $X_i$ はグリッド点 $i$ の位置を表します。また、 $S$ は粒子の形状関数（shape factor）と呼ばれるもので、ある粒子の位置（ $x_p$ ）からグリッド点に対して及ぼす影響度合い（重み）を表します。形状関数は、0次のnearest grid point（NGP）法、 1次の cloud in cell (CIC) 法が広く使われています。また近年では、2次のスプラインを用いた、よりノイズが少ない手法も使われるようになってきました。それぞれの手法に対応する形状関数を各手法の形状関数に示します。

NGPの形状関数は、

$\begin{equation} S_i(x_p) = \left \{ \begin{array}{cc} 1, & {\rm if}\ \left|\frac{x_p-X_i}{\Delta x}\right| \le \frac{1}{2} \\ & \\ 0, & {\rm else} \end{array} \right. \nonumber \end{equation}$

で表され、名前の通り、近接するグリッド点にのみ影響を及ぼす形状をしています。CIC法では、

$\begin{equation} S_i(x_p) = \left \{ \begin{array}{cc} 1-\frac{|x_p-X_i|}{\Delta x}, & {\rm if}\ \left|\frac{x_p-X_i}{\Delta x}\right| \le 1 \\ & \\ 0, & {\rm else} \end{array} \right. \nonumber \end{equation}$

のように線形補間を行います。2次のスプラインの方法では、

$\begin{equation} S_i(x_p) = \left \{ \begin{array}{cc} \frac{3}{4}-\left(\frac{x_p-X_i}{\Delta x}\right)^2, & {\rm if}\ |\frac{x_p-X_i}{\Delta x}| \le \frac{1}{2} \\ & \\ 0, & {\rm else} \end{array} \right. \nonumber \end{equation} \begin{equation} S_{i\pm1}(x_p) = \left \{ \begin{array}{cc} \frac{1}{2} \left(\frac{1}{2}\pm\frac{x_p-X_i}{\Delta x} \right)^2, & {\rm if}\ |\frac{x_p-X_i}{\Delta x}| \le \frac{1}{2} \\ & \\ 0, & {\rm else} \end{array} \right. \nonumber \end{equation}$

のような、なめらかな（微分が連続な）形状関数となっているため、ノイズを大きく減らすことができます。 pCANS では、1次元コードはCIC法を採用しており、多次元コードでは0−2次までの形状関数を選択できるようになっています。

グリッド点 $i$ の電流密度を求めるには、逆に粒子の位置と速度から

(23)¶ ${{\bf J}(X_i)} = \frac{1}{\Delta x}\sum_p^N q_p \frac{\bf u_p}{\gamma_p} S_i(x_p)$

のようにして求めます。以上では1次元の計算を例に取って見てきましたが、多次元計算の場合は、

(24)¶ $S_{i,j,k}({\bf x_p}) = S_i(x_p)S_j(y_p)S_k(z_p)$

の様に、1次元の形状関数を掛け合わせることで求めることができます。

以上の様にして、粒子とグリッド点に定義される場をカップルさせる手法を一般に Particle-in-Cell法 と呼びます。

電磁場の数値解法¶

計算グリッド¶

電磁場は電磁場の配置図のように、シミュレーション空間内に定義されたグリッド点に配置されます。このような入れ子（スタッガード）格子を Yee格子と呼びます。

このように電磁場を配置することにより、式 (4) を満たしつつ（後述）、空間2次精度で式 (9) , (10) の右辺を差分化して、表現することができます。例えば、式 (9) の $x$ 成分は

(25)¶ $\frac{B_{x,i+1/2,j,k}^{t+\Delta t}-B_{x,i+1/2,j,k}^t}{\Delta t} = -c \frac{E_{z,i+1/2,j+1/2,k}^{t+\Delta t/2}-E_{z,i+1/2,j-1/2,k}^{t+\Delta t/2}}{\Delta y} + c \frac{E_{y,i+1/2,j,k+1/2}^{t+\Delta t/2}-E_{y,i+1/2,j,k-1/2}^{t+\Delta t/2}}{\Delta z}$

のようにして、差分化します。

陰的解法¶

次に、マクスウェル方程式 (9) , (10) の数値解法について説明します。真空中のマクスウェル方程式の数値解法のうち、Finite-Difference Time-Domain (FDTD) 法 ( Yee, 1966 )が標準的な解法として知られていますが、 pCANS では星野が開発した陰的解法を採用しています。本手法は、

CFL条件が光速cに縛られず、安定に解ける

電場と磁場を同じ時間ステップ上に定義できる

時間・空間2次精度（FDTD法と同じ）

という特徴がある一方、

反復法の導入

といった、新たな数値アルゴリズムが必要となります。そのため、陽解法より計算コストがかかりますが、電磁プラズマ粒子シミュレーションでは Particle-in-Cell法が全計算コストの大部分を占めるため、それほど全体の計算コストを増やすことなく陰解法の安定性のメリットを享受できます [2] 。以下では本手法について説明します。

式 (9) , (10) を以下の様に書き換えます。

(26)¶ $\frac{{\bf B}^{t+\Delta t}-{\bf B}^t}{\Delta t} = - c \nabla \times {\bf E}^{t+\theta\Delta t}$

(27)¶ $\frac{{\bf E}^{t+\Delta t}-{\bf E}^t}{\Delta t} = + c \nabla \times {\bf B}^{t+\theta\Delta t} - 4\pi{\bf J}^{t+\Delta t/2}$

ここで、

(28)¶ $f^{t+\theta\Delta t}=(1-\theta) f^{t}+\theta f^{t+\Delta t}$

です。 $\theta=0.5$ の場合は式 (9) , (10) に一致し、2次精度のクランク・ニコルソン法と呼ばれます。 $\theta=1$ の場合は1次精度の後進オイラー法となります。 pCANS では $\theta=0.505$ など、0.5より少し大きな値を使うことにより、適度に後進オイラーによる減衰効果を取り入れて、短波長ノイズを落としています。式 (28) に従って式 (26) , (27) を変形すると、

(29)¶ $\delta {\bf B} = -{\theta c \Delta t} \nabla \times \delta {\bf E} - c \Delta t \nabla \times {\bf E}^{t}$

(30)¶ $\delta {\bf E} = \theta c \Delta t \nabla \times \delta {\bf B} + c \Delta t \nabla \times {\bf B}^{t} - 4 \pi \Delta t {\bf J}^{t+\Delta t/2}$

と、 $(\delta {\bf E},\delta {\bf B})=({\bf E}^{t+\Delta t}-{\bf E}^{t},{\bf B}^{t+\Delta t}-{\bf B}^{t})$ についての方程式が得られます。式 (30) に $\nabla \times$ を作用させて、 $\nabla \times \bf \delta E$ を式 (29) に代入して $\bf \delta B$ について整理すると、

(31)¶ $\left[1-\left( \theta c \Delta t \right)^2 \nabla^2 \right] {\bf \delta B} = \theta \left(c \Delta t \right)^2 \left( \nabla^2 {\bf B}^t + \frac{4\pi}{c} \nabla \times {\bf J}^{t+\Delta t/2} \right) -c \Delta t \nabla \times {\bf E}^t$

となり、 ${\bf \delta B}$ に対する方程式が得られます。ここで、右辺は既知量であることに注意してください。簡単のために、1次元計算（ $\partial/\partial y = \partial/\partial z=0$ ）を例にとってみましょう。グリッド点 $i$ における式 (31) の磁場の一成分（ $\delta B_s (s=y,z)$ , $B_x=const.$ ）に着目すると、

(32)¶ $&\delta B_{s,i} -\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 \left( \delta B_{s,i+1}-2\delta B_{s,i}+\delta B_{s,i-1} \right) \\ &= \theta \left( c \Delta t \right)^2 \left( \frac{B^t_{s,i+1}-2B^t_{s,i}+B^t_{s,i-1}}{\Delta x^2}+\frac{4\pi}{c} \nabla \times {\bf J}^{t+\Delta t/2}|_s \right) -c \Delta t \nabla \times {\bf E}^t|_s \\ &= b_{s,i}$

のように表されます。これは、

$\begin{equation} \left( \begin{array}{ccccc} . & . & . & . & . \\ -\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 & 1+2\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 & -\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 & 0 & 0 \\ 0 & -\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 & 1+2\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 & -\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 & 0 \\ 0 & 0 & -\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 & 1+2\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 & -\left( \frac{\theta c \Delta t}{\Delta x} \right)^2 \\ . & . & . & . & . \end{array} \right) \left( \begin{array}{c} . \\ \delta B_{s,i-1}\\ \delta B_{s,i} \\ \delta B_{s,i+1} \\ . \end{array} \right ) = \left( \begin{array}{c} . \\ b_{s,i-1} \\ b_{s,i} \\ b_{s,i+1} \\ . \end{array} \right) \nonumber \end{equation}$

のように整理できます。シミュレーション領域の境界条件を与えると式 (32) は閉じて、結局、連立1次方程式

(33)¶ $A{\bf x} = {\bf b}$

の形で表すことができます。式 (33) の解 ${\bf x}$ は $A$ の逆行列 $A^{-1}$ を作用させて

(34)¶ ${\bf x} = A^{-1} {\bf b}$

によって求めることができますが、通常 $A^{-1}$ を解析的に求めるには大きな数値コストが必要となるため、反復法によって解を求めます。 pCANS ではKrylov部分空間法の一種である共役勾配法（Conjugate Gradient method）と呼ばれる反復法を用いて式 (31) を解いています。特に、 $\theta c \Delta t/\Delta x<1$ の場合は、Aは対角優勢な行列になるため、反復法における収束は少ない反復回数で実現されます。 $\delta \bf B$ を求めた後、式 (30) に代入することによって、 $\delta \bf E$ を得ることができます。

Footnotes

束縛条件を満たすために¶

電磁場を数値的に解くなかで、束縛条件である式 (4), (7) は必ずしも保証されておらず、その破れは数値的破綻をもたらします。以下では、これら束縛条件を満たしながら時間発展を行う方法につい説明します。

磁場の発散なし¶

電磁場を電磁場の配置図のように配置した場合、 (9) は2次元において、以下のように差分化できます。

(35)¶ $B^{t+\Delta t}_{x,i+1/2,j} &= B^{t}_{x,i+1/2,j}-\Delta t \frac{E^{t+\Delta t/2}_{z,i+1/2,j+1/2}-E^{t+\Delta t/2}_{z,i+1/2,j-1/2}}{\Delta y} \\ B^{t+\Delta t}_{y,i,j+1/2} &= B^{t}_{y,i,j+1/2}+\Delta t \frac{E^{t+\Delta t/2}_{z,i+1/2,j+1/2}-E^{t+\Delta t/2}_{z,i-1/2,j+1/2}}{\Delta x}$

その結果、差分化した $\nabla \cdot {\bf B}=0$ 、つまり、

(36)¶ $\nabla \cdot {\bf B} = \frac{B_{x,i+1/2,j,k}-B_{x,i-1/2,j,k}}{\Delta x}+\frac{B_{y,i,j+1/2,k}-B_{x,i,j-1/2,k}}{\Delta y}$

は時間更新前後で維持されることがわかります。これは、式 (35) の左辺を式 (36) に代入すると明らかで、式 (35) の右辺の電場が相殺することにより $\nabla \cdot {\bf B}^{t+\Delta t} = \nabla \cdot {\bf B}^t$ の関係が導かれます。従って、計算の初期条件（t=0）で $\nabla \cdot {\bf B}=0$ となっていれば、 差分化した $\nabla \cdot {\bf B}$ は計算の時間発展の間も丸め誤差の範囲で維持されます。

ガウスの法則¶

式 (3) に $\nabla\cdot$ を作用させて (7) を代入すると、

(37)¶ $\frac{\partial \rho_e}{\partial t} + \nabla\cdot{\bf J} = 0$

の、電荷保存則が得られます。PIC法で電荷密度を電流密度と独立に計算すると、必ずしも式 (37) を満たしているとは限らないため、（静）電場に誤差が生じます。その誤差を少なく維持するために、式 (7) を満たすように電場を補正する手法がしばしば取られ、以下ではその補正方法を紹介します。式 (27) を解いた後の電場を ${\bf E}^{'}$ とし、ポアソン方程式

(38)¶ $\nabla^2 \phi = \nabla \cdot {\bf E}^{'}-4\pi\rho_e$

導入します。ここで、式 (38) 右辺は必ずしも $0$ でないことに注意してください。これを共役勾配法などで解くことによって $\phi$ を求め、

(39)¶ ${\bf E^{t+\Delta t}} = {\bf E}^{'} -\nabla \phi$

より電場の補正を行います。補正した電場は $\nabla \cdot {\bf E}^{t+\Delta t}=4\pi\rho_e$ であることは明らかです。

上記の電場補正方法では、PIC法により電荷密度を算出し、ポアソン方程式を解く必要が出てくるため、計算コストが増えます。これを回避するため、式 (37) を満たすような電流密度の計算方法が開発されており、この種の方法を電流密度計算における 電荷保存法 といいます。

電荷密度 $\rho_e$ を半整数グリッド点に定義すると、式 (37) は（2次元の場合）、

(40)¶ $\frac{\rho_{e,i+1/2,j+1/2}^{t+\Delta t}-\rho_{e,i+1/2,j+1/2}^t}{\Delta t} + \frac{J_{x,i+1,j+1/2}^{t+\Delta t/2}-J_{x,i,j+1/2}^{t+\Delta t/2}}{\Delta x} + \frac{J_{y,i+1/2,j+1}^{t+\Delta t/2}-J_{y,i+1/2,j}^{t+\Delta t/2}}{\Delta y} = 0$

のように差分化することができます。これは $(i+1/2,j+1/2)$ を中心としたセル内部の電荷の時間変化はセル表面から流入出する電荷フラックスによって変動すると考えられます。

まず、粒子が時刻 $t$ から $t+\Delta t$ 後に $(x_1,y_1)$ から $(x_2,y_2)$ に移動したとします。また、

$i_1 = floor\left(\frac{x_1}{\Delta x}+0.5\right) \\ j_1 = floor\left(\frac{y_1}{\Delta y}+0.5\right) \\ i_2 = floor\left(\frac{x_2}{\Delta x}+0.5\right) \\ j_2 = floor\left(\frac{y_2}{\Delta y}+0.5\right)$

とします。ここで、floor(x)関数はxを越えない最大の整数を返します。粒子が $(i_1,j_1)$ セル内に留まっている（ $i_1=i_2,\ j_1=j_2$ ）場合、粒子がグリッド点に電流として寄与する量は、粒子の始点 $(x_1,y_1)$ と終点 $(x_2,y_2)$ から以下のように計算することができます（ Eastwood, 1995 ）。

(41)¶ $J_{x,i_1,j_1-1/2} &= \frac{1}{\Delta x \Delta y} F_x S^{CIC}_{j_1-1/2}(y_h) \\ J_{x,i_1,j_1+1/2} &= \frac{1}{\Delta x \Delta y} F_x S^{CIC}_{j_1+1/2}(y_h) \\ J_{y,i_1-1/2,j_1} &= \frac{1}{\Delta x \Delta y} F_y S^{CIC}_{i_1-1/2}(x_h) \\ J_{y,i_1+1/2,j_1} &= \frac{1}{\Delta x \Delta y} F_y S^{CIC}_{i_1+1/2}(x_h)$

ここで、 $(x_h,y_h)=((x_1+x_2)/2,(y_1+y_2)/2)$ , $S^{CIC}$ はCIC法の形状関数、 ${\bf F}=(F_x,F_y)$ は

(42)¶ $F_x = q_p\frac{x_2-x_1}{\Delta t} \\ F_y = q_p\frac{y_2-y_1}{\Delta t}$

で与えられる電荷フラックスを表します。

以上では、同じセル内に留まっている場合のみを議論しましたが、一般には隣のセルに移ることも考えられます。その場合の取扱いの違いを図に示します。 Villasenor & Buneman (1992) では、始点と終点を直線で結び、図の場合、セル境界によって区切られた３つのセグメントそれぞれのセル内で上記の計算を行っています（Villasenor-Buneman法）。この場合、セル境界での交点を計算するのに分岐命令が多数必要となるため、計算コストが大きくなります。その問題を解消するために、 Umeda et al. (2003) は、図のように、途中で電荷密度が定義されるグリッド点を経由することでセル境界での交点の計算を少なくした方法を開発しています。この場合、軌道は直線とならないため、ジグザグ法と呼ばれています。一方、 Esirkepov (2001) は、電荷を（2次元の場合）q/2の2つに分けて、それぞれの経路を直線で表す手法を開発しました。これを Density Decomposition法 と言います。本手法の利点として、