3. プラズマシミュレーション
田中正明
ペガサスソフトウェア株式会社

mtanaka@psinc.co.jp

1 はじめに
1.1 プラズマシミュレーションの目的
一般にシミュレーションの目的とは、何のために行うのか（Needs）、どのような利益があるのか（Beneﬁt)、またそのシミュレーションにはそしてどの程度信頼できるのか（Proof)、ということが求められる。プラズシミュレーションの場合、プラズマの応用技術として光学的な応用を含むエネルギー分野、化学的、物理的な応用を含む物質・材料分野、廃棄物・廃ガス処理などの環境分野そしてロケット・衛星の推進エンジンへの応用を含む宇宙分野など、多岐にわたっているため、その目的も多種多様である。またプラズマは主に熱プラズマ、完全電離プラズマそして非平衡低温プラズマなどに分けられる。本稿では、非平衡低温プラズマを対象としたプラズマシミュレーションに関して述べる。

1.1.1 何のために行うのか (Needs)
プラズマシミュレーションを行えば何かが判る、という目的が曖昧な理由からでは、単なる時間の浪費であり、シミュレーションは意味が無いという結論に陥ることになる。得られた測定可能な物理量を説明する上で、また現象のみが分かっているとき、その現象を説明する上で、何が支配的であるかを知ることは重要である。従って、シミュレーションにより何が得られれば、その結果をどのように今後に活かすことができるか、という問題
点と方向性に関するシナリオを持つことにより、シミュレーションの目的は明確になる。

1.1.2 どのような利益があるのか (Beneﬁts)• 実験・試作コストの低減を行う• 測定が難しい現象を解明することにより予測の手がかりとなる• 物理量を数値化・可視化することにより現状の把握と仮想実験が可能となる• 運転条件を変化させたときの実験データが得られているとき、定性的ではあっても、目的とするものの最適化が可能となる。

1.1.3 どの程度信頼できるのか (Proof)
プラズマシミュレーションには気相中、構造物表面での物理的、化学的そして光学的な実現象を正確に取り入れるモデルは確かに不可能に近い。なぜなら、気相中での複雑な反応に伴う基礎データの欠如、もしくは信頼性そして取り扱う境界条件が時々刻々変化している場合が多いからである。しかし、運転条件、装置条件に適合した、理論的な背景に基づくモデルによるシミュレーションは、実現象を模擬したものであり、得られる数値も定性的には信頼に値するものである。また、物理現象が比較的明確な対象によっては定量的に信頼できる結果が得られる。

2 プラズマシミュレーションに用いられる計算手法
大別して気相中でのプラズマ（荷電粒子）とラジカル・励起種・供給ガス（中性粒子）の挙動をシミュレーションする計算手法と基板における微細形状の時間変化をシミュレーションする計算手法に分けられる。プラズマ装置内の気相シミュレーションにおける解析領域と微細形状のシミュレーションにおける解析体領域の比は 105 以上になる。また、近年常圧下での数 mm 以下の気相シミュレーションも頻繁に行われるようになってきた。

2.1 PIC/MCC 法 (荷電粒子のための粒子モデル)
PIC/MCC 法は、Particle-in-Cell(PIC) 法と Monte Carlo Collision(MCC) 法を組み合わせて装置内のプラズマ（荷電粒子）の挙動をシミュレーションする手法である。PIC（Particle-In-Cell）法はおよそ 40 年以上前に流体方程式の解法として開発された。その後
核融合プラズマにおける電磁流体方程式解法にも使用され、およそ 20 年前から弱電離プラズマ（非平衡低温プラズマ）シミュレーションに使用され、今日に至っている。計算手法の詳細は例えば Birdsall の教科書 [1]、南部 [2] などを参照して頂きたい。計算手順の概要は次のようなものである。

1. 静電位に関するポアソン方程式を解く
ポテンシャル計算のために各セルの電荷密度を求め、電位の境界条件を考慮した上で各セルのポテンシャルと電界を計算する。

2. 荷電粒子の運動方程式を時間積分する

セル境界上の電界、 磁界を粒子位置に内挿して ∆t 後の粒子の速度と位置を計算する。

3. 壁面での境界条件
壁面での荷電粒子の消滅や散乱、スパッタリング、2 次電子放出など境界上での粒子の計算を行う。

4. 中性粒子との衝突
荷電粒子と中性原子・分子間の衝突の種類 (弾性散乱、電離、励起、解離など) を衝突断面積に従った確率で決定し、衝突後の速度を計算し、電離によるイオンと電子の生成も考慮する。
この手順を、プラズマの密度分布や電位分布が (周期) 定常状態になるまで短い時間ステップ ∆t で繰り返していく。上記の 1,2 のみを繰り返して計算するのが PIC 法であり、4 がMCC 法である。PIC/MCC 法では、例えば 108 ∼ 1010 個の荷電粒子をまとめて、1 個の超粒子として扱い、この超粒子の運動を追跡する。
メッシュ分割した解析空間内に電子とイオンの超粒子を配置し、その超粒子の分布から決まる空間電荷をグリッド上に振り分ける。そして与えられた境界条件 (電極の電圧やアースの配置など) の下で (誘電率 ε が一様でない場合を考慮した) ポアソン方程式∇ · ε(r)∇φ(r) = −ρ(r)を解いて電位 φ を計算し、さらに電位の勾配から電界を計算する。荷電粒子の運動方程式は

md2 rdr=q(E+v×B),v=dt2dtである。ここで、 m: 荷電粒子の質量; E: 電界; B: 磁界; q: 荷電粒子の電荷; r: 荷電粒子の位置座標; v: 荷電粒子の速度である。 この運動方程式を、良く知られている蛙飛び法やベルレ法などで定式化して ∆t 後の位置と速度を計算すればよい。

2.2 DSMC 法 (中性粒子のための粒子モデル)
プラズマ装置内の原料ガスの流れ場などのように圧力が低い流れ場について解析する場
合、もしくは常圧下であっても解析体系が数 µm 以下となるとき、気体を連続体であるとみなして取り扱うナビエ-ストークス方程式は適用できず、Boltzmann 方程式にまでさかのぼって考えなくてはならない。ところが、この Bolzmann 方程式は複雑な非線形積分微分方程式であり、その取扱いは非常に難しい。DSMC 法 (Direct Simulation Monte Carlomethod、 直接シミュレーションモンテカルロ法) は Boltzmann 方程式を直接解くのではなく、Boltzmann 方程式の基になっているそれぞれの粒子の衝突過程を確率的に取り扱うことによって流れ場を解析するといった方法である。この手法の詳細は、南部 [2]、Brid [3]
そして南部 [4] の教科書を参照して頂きたい。

DSMC 法では、気体分子を PIC/MCC 法と同様、超粒子として取り扱う。真空装置内といえども装置内の全ての気体分子を計算機で追跡することは不可能で、例えば 1015 ∼ 1020個の気体分子をまとめて 1 個の超粒子として扱う。数万〜数 100 万個程度の超粒子を取り扱い、各粒子の位置、速度、内部状態等をメモリー内に記録し、衝突や境界の影響によってそれらの値を更新していくという手続きを繰り返す。
衝突過程については、与えられた衝突断面積に従って確率的に引き起こし、衝突後の 2粒子の速度もその衝突物理モデルに従い確率的に決定する。具体的には、流れ場を適当な大きさのセルに分割し、 同一セル内の 2 個の粒子を物理モデルに基づいた確率則に従って衝突対として選択するというものである。衝突後の粒子対の速度は、エネルギー、運動量等が保存するように決める。希薄気体の場合、流れ場の物理量の変化の代表的な長さは平均自由行程程度であることから、セルのサイズは平均自由行程程度の大きさとする。

DSMC 法により直接得られるものは、多数の粒子の位置、速度等のサンプルデータである。このデータを加工することにより、粒子の速度分布関数、密度分布、流速分布、温度分布などを得ることができる。統計誤差はサンプル数の 1/2 乗に反比例するので、時間平均のゆらぎが小さくなり、十分な精度を持つまでサンプリングを行うことになる。
またプラズマ装置内もしくは真空装置内内では、密度差が大きく異なる原子・分子が存在し、小さな密度を持つ原子・分子の挙動を正確にシミュレーションするには、従来のDSMC 法ではすべての粒子種の重みは同一であるため、膨大なサンプル数が必要であるが、Weight algorithm(南部 [2]) を採用することにより、ほぼ同数のサンプル数で行える。

2.3 ドリフト−拡散モデル（荷電粒子のための流体モデル)
ガス圧が高くなると荷電粒子が分子と衝突する際の平均自由行程が小さくなるので、PIC/MCC 法のような粒子モデルを用いようとすれば空間のメッシュ分割幅と時間進展のδt をともに相当小さくする必要があり、計算時間が実用的な範囲を超えてしまう。そのような場合には、電子やイオンの集団をそれぞれ連続体 (流体）として取り扱うモデルを用いる。
流体モデルでは、電子とイオンの各々についてマクロな量としての数密度 nα (r, t) 、流速ベクトル v α (r, t) （r は位置座標、t は時刻）を未知変数として記述される連続の式∂nα+ ∇ · Γα = Rα∂tおよび運動方程式

Γα = ±nµα E − Dα ∇nαを電位のポアソン方程式とカップルして解くことにより、プラズマ密度やポテンシャル、および電子・イオンのフラックスを求める。ここで添え字の α は荷電粒子の種類を識別する指標で、α = e は電子を、α = ion はイオンを意味する。また必要な場合には電子の平均エネルギー ¯(r, t) についてのバランス方程式∂(ne ¯) + ∇ ·∂t55 + eΓe · E = −ne ν ¯Γe ¯ − ne De ∇¯33も同時に解く。これらの式において Γα は粒子フラックス（ = nα v α ）、E は電界、µα は移動度、Dα は拡散係数、Rα は衝突電離や電子付着による生成/消滅レート、ν は電子衝
突による電子エネルギー損失の衝突周波数、e は素電荷である。
ここに記述されている運動方程式は厳密なものではなく、ドリフト−拡散モデルといわれる近似式である。厳密には、静磁場 B も考慮すると、

 mβ∂(mnv) + ∇ · (mnvv) = −∇(nkB T ) + nq(E + v × B) + mn
νβ (v β − v)
∂tm + mββとなる。ここで、添え字のない m、n そして v は荷電粒子を示し、添え字は中性粒子を示す。
ドリフト−拡散モデルは• 応力は等方的であり、スカラー圧力 p に帰着される。また温度 T と圧力の間にはp = nkB T の関係がある。
• 運動量移行に関する緩和時間を τm 、平均速度が変化する時間スケールを τv とすると τm /τv  1 である。

• 粒子速度のランダムな変動（熱速度）の大きさを vT としたとき、v  vT である。という仮定のもとに運動方程式の時間微分項と慣性項を無視したものであるが、広く使用されている。
これらの式を解くためには、移動度、拡散係数といった輸送パラメータや、電子やイオンの生成/消滅レート、電子エネルギー損失の衝突周波数などの値を与える必要があるが、これらは各種衝突の衝突周波数 ν によって定まり、その衝突周波数は速度分布関数と衝突断面積の積の (速度空間における）積分によって計算されるから、速度分布関数の型に関して何らかの仮定が必要である。
電子の場合分布関数に対する仮定は、

1. マクスウェル分布を用いる
この場合衝突周波数 /N （N は分子の密度）は電子温度の関数となるから、プラズマシミュレーションを行う前に適当な範囲の電子温度に対して、その温度に対応するマクスウェル分布を用いて ν(Te )/N を計算で求めておき、それらのテーブルを作っておく。移動度と拡散係数は運動量移行衝突周波数から計算する。

2. 一様電界中の電子スウォームの分布関数を用いる
この場合には通常２項近似のボルツマン方程式を別途解く。衝突周波数は換算電界E/N または電子の平均エネルギー ¯ の関数となる。これも事前にある範囲の E/Nに対して ν(E/N )/N を計算で求めておき、E/N またはその換算電界における電子平均エネルギーの関数としてテーブルを作っておく。
なお 1 の場合や 2 で電子平均エネルギーの関数としてテーブルが作成されている場合には、電子エネルギーバランス式を解いてセルごとの電子エネルギー（温度）を各時間ステップごとに求める必要がある。
イオンに関しては、イオンのエネルギーは電離や電子付着にほとんど影響しないので、エネルギーバランス式を解かずに温度は一定値と仮定して（例えば原子・分子の温度と同程度）連続の式と運動方程式のみを解く。移動度および拡散係数については、換算電界の関数として実験によって得られている値を用いるのが望ましいが、それらの情報が得られない場合には、衝突断面積を一定とした剛体球モデルとマクスウェル分布を仮定した平均速度から計算した値を用いる。文献 [5] にはイオンの移動度と拡散係数に関する理論と実測データがよくまとめられている。

2.4 緩和連続モデル（荷電粒子のための流体モデル)
連続の式の粒子フラックスにドリフト速度を用い、運動量保存則と電子のエネルギー保存則にそれぞれの緩和時間を取り入れた真壁等によるモデルがある。(文献 [6] [7])2.5 ハイブリッドモデル (荷電粒子のための流体モデルと粒子モデルとの混合モデル)
ハイブリッドモデルには以下のモデルが存在する。

1. 電子のモンテカルロシミュレーションを併用する
電子の速度分布関数（実際にはエネルギー分布関数）についてのみ PIC/MCC 同様の手法を用いて計算セルごとに値を求める。ただしこの場合にはプラズマ密度の空間分布によってプラズマ中の電界が変化するので、流体モデルによるプラズマ密度および電界の計算とモンテカルロ計算を数周期 ∼ 数十周期ぶんずつ交互に行いながら計算を進めていく。計算時間が多少増えるが、分布関数を仮定する必要が無いのが利点である。得られた電子エネルギー分布関数から各セルにおける電子温度、電子衝突による反応レートが求まる。（文献 [8]）2. 電子を PIC/MCC で計算しイオンは流体モデルで計算する
マグネトロンスパッタ装置におけるマグネトロンプラズマシミュレーションはPIC/MCC 法のみで計算する方法と電子のみ PIC/MCC 法を用い、イオンの生成率を基にしたイオンに関する流体モデルによりイオン密度を求める方法がある。マグ
ネトロンプラズマに限らず、電子の非線形性が強く、流体モデルで近似できない場合に採用される。（真壁 [9])3. 電子を流体モデルで近似しイオンの挙動はモンテカルロ法を用いる
ECR プラズマシミュレーションにおいて磁力線に沿って電子の流体モデルの定常解を求め、イオンの密度はモンテカルロシミュレーションにより求める方法。（文献 [10]）2.6 ICP 装置における誘導電界
次に電磁場とのカップリングに関して述べる。ICP 装置内のプラズマシミュレーションを行う場合には、電磁場についてのマクスウェル方程式から導かれる式∇2 Eθ − iωµ0 σEθ = 0を解いて誘導電界の周方向成分 Eθ を求め、プラズマに注入されるパワーを計算する。ω
はコイル電流の角振動数、µ0 は真空の透磁率、σ はプラズマの導電率である。σ はプラズ
マ (電子）の密度に依存するので、この場合もプラズマ密度の計算を数周期行うごとにこの方程式を解いて Eθ を更新する必要が有る。(文献 [11])2.7 荷電粒子シミュレーションと中性粒子シミュレーションとのカップリ
ング
中性ガス流れとのカップリングに関して述べる。プラズマ中では原料となる分子の他に電子やイオンとの衝突によって生じる励起種や原子が生成される。また一般には装置の一部から原料ガスが供給され、排気口から流出させるためガスの流れが生じ、それによって分子、原子密度の空間分布が変化する。中性粒子種の密度変化は電離や付着に影響を与えるし、原子や励起種の密度やフラックスに関する情報が必要な場合も多い。したがってできる限りガスの流れも解くことが望ましい。ガス圧が高い場合には、ガス流れに関してもDSMC 法ではなく、ナヴィエ−ストークス方程式などの連続体モデルを用いる。その解法は多く存在するが、多数の成分が存在する場合でも比較的容易に計算できるように、全ての中性ガス種をあわせてひとつの混合流体として考える。この混合流体に対して運動方程式（ナヴィエ−ストークス方程式）とエネルギー方程式を解き、各成分に関しては対流−拡散型の輸送方程式（質量保存則）を解く事により成分ごとの密度分布を求める。その際、中性粒子どうしの反応も反応速度定数を与えることによって考慮される。なお電子・イオンに関する物理量が変化する時間スケールに比べて中性ガスに関するそれはオーダーが数桁大きいので、両者の計算は切り離しておこなう。普通は荷電粒子に対する計算を数周期ぶん行ってからガス流れの計算を数 ms 程度行うことを反復する。

3 計算例
3.1 対向ターゲットマグネトロンスパッタ装置
マグネトロンスパッタ装置は、薄膜製造のための代表的な装置の一つである。マグネト
ロンスパッタ装置内で生じている物理現象は、静磁場中における DC(RF) グロー放電、高エネルギーイオンがターゲット原子を弾き出すスパッタリング現象、スパッタされた原子の輸送、基板へスパッタ原子が堆積する過程などがある。マグネトロンスパッタ装置内のいろいろな物理量をシミュレーションにより推測する。

2 次元 x-y 直交座標系で計算する。対向ターゲットのターゲット長は 12cm、ターゲット間距離は 8cm そしてターゲット厚さは 1cm とした。その他の計算条件は ガス種はアルゴンと酸素、ガス圧は共に 0.25Pa とした、また磁石の残留磁束密度 1.0T、 ターゲットの印加電圧 -400V、 ターゲットの材質 Si 、イオンがターゲットに入射した場合の 2 次電子放出係数 0.1、その他の境界に電子やイオンが入射した場合は消滅する、などである。

図1磁束密度分布。ターゲット上での、磁束密度は最大 650G 程度である。

3.1.1 静磁場解析
静磁場解析は、A − φ 法を有限要素法で定式化して計算しする。磁束密度分布を 図 1 に示す。

3.1.2 プラズマシミュレーション
先に求めた静磁場中でのプラズマの状態を、PIC/MCC 法で解析する。電子とイオンの初期数密度を 1.0×1014 /m3 として、プラズマが定常状態になるまで計算を続ける。代表的な計算結果の図を以下に示す。以下の図はターゲット – 基板近傍の領域のみ表示してある。

図2

電子数密度分布 [/m3 ]図4O2 (+) 密度分布 [V]

図3

図5Ar(+) 密度分布 [C/m3 ]ターゲットに入射する Ar(+) フラックス分布図 6ターゲットからスパッタされるフラ ッ ク ス の 分 布 。最 大 の フ ラ ッ ク ス は2.8×1019 /m2 /s.

図7

Si 原子の基板へのフラックス分布。
最大フラックスは 4.7×1018 /m2 /s.

3.1.3 スパッタ粒子の輸送シミュレーション
テスト粒子モンテカルロ法を用いて、Ar、酸素ガス中でのスパッタ粒子の輸送シミュ
レーションを行う。スパッタ粒子はターゲットから流入して Ar、酸素 ガス中を飛行し、基板に到達する。流入するスパッタ粒子の速度分布関数はスパッタリングシミュレーションから得られるものである。
計算結果のスパッタ粒子の基板へのフラックス分布から直接的に膜厚分布が推定できる。
また、基板へ到達したスパッタ粒子の速度分布関数も得られ、膜質の評価や微細なトレンチへの埋め込み特性の評価にも結びつけることが期待できる。

3.2 ドライエッチング装置
ドライエッチング装置は、半導体製造などの微細加工の際に基板に入射するイオンの異方性と中性ラジカルの等方性を制御する。基板の種類により、使用されるガス種などの運転条件と装置の種類がさまざまである。ここでは Cl2 ガスを用いた誘導結合型装置内の計算例 (文献 [12]) を示す。荷電粒子種に関してはハイブリッドモデルを用い、中性粒子種については DSMC 法を用いた。
図 8 が軸対称モデルの計算体系である。コイルの条件は、パワーが 500W 、周波数が13.56 MHz である。ガス圧は 2Pa とした。考慮した粒子種は Cl2 、Cl、正負イオンであり、弾性衝突、電離、電子付着、イオン−イオン再結合など 12 種類の衝突過程が考慮されている。ガスの流入量は Cl2 が 50 sccm である。
図 9 はコイルによって生成される誘導電界の分布である。この振動する電界から電子が加熱され、装置内にプラズマが生成され維持される。
図 10 はプラズマに吸収されるパワー密度分布、図 11 は電子の数密度分布である。中心に (∼ 1017 ）のプラズマが生成されている。また、図 12、 図 13 はそれぞれ電子のエネルギー分布関数、電子温度分布である。
図 14 は 電位分布、図 15∼ 図 17 はそれぞれ Cl2 、Cl(+) そして Cl ラジカルの数密度分布である。

図9
図8エッチング装置の計算モデル

コイルによる誘導電界 [V/cm]

図 10 パワーの吸収分布 [W/cm3 ]

図 11 電子密度分布 [/m3

図 12 電子エネルギー分布

図 13 電子温度分布 [eV]

図 14 電位分布]

図 16Cl(+) の密度分布 [/m2 /s]

図 15Cl2 (+) の密度分布 [/m3 ]図 17Cl の密度分布 [/m3 ]4 おわりに
すでに、熱 · 流体 · 構造 · 電磁波などのシミュレーション技術は、建築物 · 自動車 · 飛行機 · 電子機器などの設計に広く使われており、実験および試作コストの大幅な低減に成功している。またこの分野では、使い勝手の良い商用ソフトウェアも数多く販売されており、シミュレーション技術の普及に寄与している。
プラズマシミュレーションに関する研究も多くの研究者によって広く行われており、プラズマ現象の理解に大いに役立ってきた。しかし、未だ装置のデザインツールとして技術者に広く受け入れられるまでには至っていないのが本当のところであろう。理由としては、計算負荷が大きい、断面積データの不足、表面で起こる化学反応がよくわからない、装置内の測定データがほとんど手に入らず計算結果の精度の評価も難しい、などといった点があげられる。少なくとも原子・分子データの整備に関しては、日本独自の組織が必要であろう。

参考文献[1] C.K. Birdsall and A.B. Langdon: P lasma P hysics via Computer Simulation(McGraw-Hill, New York, 1985).
[2] Nanbu, K. IEEE Trans. Plasma Sci. Vol. 28, 971 (2000).
[3] G.A. Bird. Molecular Gas Dynamics and the Direct Simulation of Gas Flow.
Clarendon Press · Oxford, 1994.
[4] 南部健一. 原子 · 分子モデルを用いる数値シミュレーション 第 3 章モンテカルロ法の基礎. コロナ社, 1996.

[5] E.W. McDaniel and E.A. Mason: T he M obility and Dif f usion of Ions in Gases(Wiley, New York, 1973).
[6] K.Okazaki,T.Makabe and Y.Yamaguchi Appl. Phys. Lett. Vol. 54, 1742 (1989)[7] T.Makabe N.Nakano and Y.Yamaguchi Phys. Rev. A Vol. 45, 2520 (1992)[8] T.J. Sommerer and M.J. Kushner. Numerical investigation of the kinetics andchemistry of rf glow discharge plasmas sustained in He, N2, O2, He/N2/O2,He/CF4/O2, and SiH4/NH3 using a Monte Carlo-ﬂuid hybrid model. Journalof Applied Physics, Vol. 71, pp. 1654–1673, 1991.
[9] 真壁利明. プラズマエレクトロニクス 培風館, 2000.
[10] H.M. Wu, D.B. Graves and M. Kilgores. Two-dimensional simulation of compact

ECR plasma sources. Plasma Source Sci. Technol., Vol. 6, pp.231–239,1997.
[11] P.L.G. Ventzek, M. Grapperhaus, and M.J. Kushner. Investigation of electron
source and ion ﬂux uniformity in high plasma density inductively coupled etchingtools using two-dimensional modeling. J. Vac. Sci. Technol. B 12(6), pp. 3118–3137, Nov/Dec 1994.
[12] M. Shiozawa and K. Nanbu. Coupling of plasma and ﬂow in materials processing.
Thin Solid Films Vol. 457,pp.48–54,2004.