JP2009174028A - マイクロ波プラズマ解析プログラム - Google Patents

Abstract

【課題】マイクロ波プラズマ分布を数値解析するプログラムを提供すること。
【解決手段】計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、反応毎のパラメータ、背景ガスの圧量及び温度、入力電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付け（Ｓ１〜Ｓ４）、電子密度分布から誘電率分布を求め（Ｓ５）、形状情報、メッシュ情報、境界条件、誘電率分布、並びに、電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式から電界分布を求め（Ｓ６）、電界分布から電子温度を求め（Ｓ７）、電界分布、電子温度、形状情報、メッシュ情報、フィッティングパラメータ、反応毎のパラメータ、背景ガスの圧力及び温度、並びに、電磁波の周波数の下で、電子輸送方程式を用いた数値計算によって電子密度分布を求め（Ｓ８）、電子密度分布の変化が収束するまで処理（Ｓ５〜Ｓ９）を繰り返す。
【選択図】図３

Description

本発明は、マイクロ波プラズマを解析するプログラムに関し、特にダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値解析することができるマイクロ波プラズマ解析プログラムに関する。

マイクロ波プラズマを用いたＣＶＤ装置の設計、合成条件の割り出し、合成速度の見積もり等には、プラズマの中の電子密度、電子温度やパワー密度等の情報が必要である。しかし、実験によりマイクロ波プラズマを観測することは困難であり、そのため、数値計算に頼っているのが現状である。

真空の電磁界を計算する方法、及びそのためのコンピュータソフトウェアは公知であり、複数のコンピュータソフトウェアが市販されている。しかし、マイクロ波プラズマを考慮した電磁界の解析機能を持つ、市販のコンピュータソフトウェアは無い。

マイクロ波プラプラズマでは、高周波電磁界の下で、電子、イオン及び中性ガスが混在している。従って、変数となる物理量には、次の(i)〜(iii)がある。
(i)電界及び磁界（２つの３次元ベクトル変数）
(ii)電子の密度、速度及び温度（２つのスカラー変数と１つの３次元ベクトル変数）
(iii)イオン及び中性ガスの、密度、速度及び温度（（２つのスカラー変数と１つの３次
元ベクトル変数）×種の数）
なお、(iii)に関しては、例えば水素の放電であれば、Ｈ、Ｈ₂、Ｈ⁺、Ｈ₂ ⁺、Ｈ₂ ^*など
のそれぞれの種について変数を導入する必要がある。本明細書において、上付き添字の＊は原子、分子の励起状態を表す。

一般には、これら全ての物理量に対応する支配方程式を解くことにより、各物理量の空間分布を決定することができる。従って、イオン及び中性ガスの種として、それぞれ一つだけを考慮するとしても、それらの定常状態における分布を求めるには、合計２１（＝６＋５＋５＋５）の偏微分方程式を連立して解かなくてはならない。

このように、マイクロ波プラズマの支配方程式は膨大且つ複雑であり、それを用いた数値計算によるシミュレーションは非常に困難である。即ち、数値的に解くべき微分方程式の数が非常に多い、数値計算における解の収束性が悪い、計算に時間が掛かる（計算コストが高い）、設定すべき未知のパラメータ量が膨大であるなどの問題がある。そして、これらのために、既存のソフトウェアでは、精度の良い計算結果を得ることができなかった。

一方、下記の非特許文献１及び非特許文献２には、ＣＶＤ装置におけるマイクロ波プラズマの電磁界分布を解析する方法が開示されている。

また、本発明者による下記非特許文献３には、物理量を高周波振動成分と平均速度成分に分離し、パラメータフィッティングを併用することによって、真空の電磁界の計算と同程度の計算量で、マイクロ波プラズマを考慮した数値解析が可能となることが開示されている。
H. Yamada, A. Chayahara, Y. Mokuno, Y. Soda, Y. Horino, and N. Fujimori, Diamond & Related Materials, 14(2005) pp.1776-1779 中川貴嗣、金載浩、鳥羽孝幸、桂井誠；電気学会論文誌Ａ、123(2003) pp.481-489 H. Yamada, A. Chayahara, Y. Mokuno, Journal of Applied Physics, 101, 063302(2007)

しかし、非特許文献１では、特定の反応に特化した式を用いて解析しているので、適用範囲が限定される。また、非特許文献１では、（６）式の左辺の微分係数の評価が粗く、特性長Ｌを新たに導入して、▽・（ｎ_eｕ_e）をｎ_e｜ｕ_e｜／Ｌ（ここで、ｕ_eは電子の流
速、即ち３次元ベクトルであり、｜ｕ_e｜は流速の絶対値である。）として評価している
こと、及び、電子の速度ｕ_eに高周波振動成分（同文献中（１）式）を採用していること
から、電界強度と電子密度は比例関係に束縛され、従って、高電子密度の場合に精度良い計算結果を得ることができない問題がある。さらに、非特許文献１では、再結合反応率ｋ_Rが電子密度の絶対値を適当な大きさとなる様に調節するためのフリーパラメータとして
扱われている。そして、先に述べた理由から計算精度が低いため、非特許文献１では、ｋ_Rは実際の値（例えば、M. T. Leu, et al., Physical Review A 8, 413 (1973)参照）よ
りも１００倍程度大きい値にしなくてはならなかった。

非特許文献２では、電磁気に関するマクスウェルの方程式に加えて、プラズマの方程式として、電子及びイオンの粒子保存の式及びエネルギー保存の式の両方を用いて数値計算を行うことが必要であり、計算負荷が大きい問題がある。

非特許文献３では、数値的に解くべきプラズマの方程式、即ちマイクロ波プラズマ中の電子のエネルギー輸送方程式に対して、３つのパラメータ（ａ_L、ｐ_L及びＥ_L）を含む試
験関数

を導入し、パラメータフィッティングを行うが、この方法ではある程度の誤差が発生することになる。また、この誤差はT_eに依存して変化するため、計算領域中で誤差が一様でなく、誤差評価が容易ではない。特に試験関数が適切でなければ、計算精度が悪くなり、場合によっては解けないことも起こり得る。

また、マイクロ波プラズマの支配方程式系を解くに当たり初期の電子密度分布を仮定する必要があるが、その与え方に関する手法が開示されていない。電子密度分布は真空領域内に局所的に存在する場合が多いが、例えば数値計算の対象とする空間全体に渡って一様な電子密度を仮定した場合には、局所的に存在する解へ収束するまでの計算コストが必要となる。また、非現実的な不安定な解へ収束することも起こり得る。

本発明は、上記した問題を解決するためになされたものであり、ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を、所望の精度で数値解析することができるマイクロ波プラズマ解析プログラムを提供することを目的とする。

本発明の目的は、以下の手段によって達成される。

即ち、本発明に係るマイクロ波プラズマ解析プログラム（１）は、ダイヤモンド合成に
使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値計算するコンピュータ用プログラムであって、
コンピュータに、
数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
電子密度分布から誘電率分布を求める第２の機能と、
前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第３の機能と、
反応種をｓで表し、電荷素量をｅ、電子質量をｍ_e、電子の衝突振動数をν_e、電子温度をＴ_e、電磁界角周波数をω、電界をＥ、反応種ｓの密度をｎ_s、並びに、電離・励起・再結合等の各種反応をＸで表し、反応Ｘのエネルギー閾値と反応率とをそれぞれｋ_s ^(X)及びＥ_s ^(X)として、前記電界分布の下で、

の方程式を数値計算によって解き、電子温度を求める第４の機能と、
前記電界分布、前記形状情報、前記メッシュ情報、前記電子温度、前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、反応種ｓのモル分率を

として、

の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第５の機能とを実現させ、
前記初期の電子密度分布が、数値計算の対象とする空間の各場所での電子密度を０と仮定した電界強度分布において、真空領域に相当する部分で電界強度が最大となる点で最大値を持つ分布であり、
前記第５の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第２の機能〜前記第５の機能を順に繰り返して実現させることを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（２）は、上記のマイクロ波プ
ラズマ解析プログラム（１）において、
前記第５の機能が、前記電界分布において、電界強度が０以外の値であるメッシュに関してのみ、新たに電子密度を求め、電界強度が０のメッシュに関しては、電子密度を０とする機能であることを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（３）は、上記のマイクロ波プラズマ解析プログラム（１）又は（２）において、
前記第４の機能で使用する方程式が、

であることを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（４）は、上記のマイクロ波プラズマ解析プログラム（３）において、
前記第４の機能で使用する方程式が、

であり、
前記第４の機能で使用する数値計算が２分法であり、該２分法の初期値である電子温度を０ｅＶ以上１０ｅＶ以下とすることを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（５）は、上記のマイクロ波プラズマ解析プログラム（３）において、
前記第４の機能で使用する方程式が、

であり、
前記第４の機能で使用する数値計算が反復法であり、該反復法の初期値である電子温度を０ｅＶ以上１０ｅＶ以下とすることを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（６）は、上記のマイクロ波プラズマ解析プログラム（３）において、
前記第４の機能で使用する方程式が、

であり、
前記第４の機能で使用する数値計算が変分法であり、電子温度が０ｅＶ以上１０ｅＶ以下の範囲で解を探索することを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（７）は、上記のマイクロ波プラズマ解析プログラム（１）〜（６）の何れかにおいて、
前記反応が、ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺、ｅ＋Ｈ₂→ｅ＋Ｈ₂ ^*、ｅ＋Ｈ₂→ｅ＋２Ｈ、及びＨ₃ ⁺＋ｅ→３Ｈであり、
反応ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺の反応率が１．０×１０^-14ｍ³／sec、エネルギー閾値が１５．４ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
反応ｅ＋Ｈ₂→ｅ＋Ｈ₂ ^*の反応率が６．５×１０^-15ｍ³／sec、エネルギー閾値が１２．０ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
反応ｅ＋Ｈ₂→ｅ＋２Ｈの反応率が１．０×１０^-14ｍ³／sec、エネルギー閾値が１０．０ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
反応Ｈ₃ ⁺＋ｅ→３Ｈの反応率が１．０×１０^-13ｍ³／sec、エネルギー閾値が０．０ｅ
Ｖ、反応種の密度が１．５×１０²⁴ｍ^-3であることを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（８）は、上記のマイクロ波プラズマ解析プログラム（１）〜（７）の何れかにおいて、
前記形状情報が円盤状のサセプタの形状情報を含み、
前記サセプタの円形表面の中心からの距離をｒｍとし、前記サセプタ外部のｒ＜０．０２である領域内で、初期電子密度が２．５×１０¹⁷×[１．０＋ｃｏｓ（πｒ／０．０２
）]であり、前記サセプタ外部のｒ≧０．０２である領域で初期電子密度が０であること
を特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（９）は、
初期の電子密度分布が、数値計算の対象とする空間の各場所での電子密度を０と仮定した電界強度分布において、真空領域に相当する部分で電界強度が最大となる点で最大値を持つ分布であることを特徴としている。

本発明のマイクロ波プラズマ解析プログラムによれば、ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を、従来よりも短時間で、且つ精度良く数値解析することができる。これによって、ダイヤモンド合成に適したＣＶＤ装置を効率的に設計することが可能となる。

また、数値計算の対象とする空間全体に渡って一様に所望の精度で電子温度を求めることができ、従来のパラメータフィッティングを行う場合よりも精度良く電子温度を求めることができる。従って、計算結果の信頼性がより高くなる。

また、初期の電子密度分布の与え方が工夫されているため、計算コストを抑制しつつより安定に現実的な解を求められる。

以下、本発明に係る実施の形態を、添付した図面に基づいて説明する。

まず、本発明の前提条件及び解くべき方程式に関して説明する。本発明が解析対象とするのは、ダイヤモンド合成のためのマイクロ波プラズマ放電である。９割程度が水素によって構成される原料ガスを比較的高いガス圧力（数Torr以上）で使用する。マイクロ波プラズマ発生装置には、図１、図２に示すような断面図の構造を有するＣＶＤ装置を想定している。これらのＣＶＤ装置は、サセプタの天面中央で放電が起こることが想定されており、各部の寸法、形状に応じて、形成されるマイクロ波プラズマ分布を解析することが要求される。

図１において、マイクロ波発生部（図示せず）で発生されマイクロ波入力部１から導入されたマイクロ波は、ウェーブガイド２を通ってアンテナ７を介して第１円筒部３及び第２円筒部５で定在波を形成する。石英板４を挟んで第１円筒部と隔離された第２円筒部５内の真空領域６には、プラズマ源となるガスが、ガス供給部（図示せず）から所定の圧力で供給されており、サセプタ８の天面中央で放電が起こり、真空領域６にプラズマが発生する。そして、サセプタ８に載置された物（図示せず）の表面にダイヤモンドを積層させる。また、図２は、図１と構造は異なるが、図１と同様にサセプタ８の天面中央で放電が起こり、真空領域中にプラズマが発生する。図２において、図１と同じ機能をする構成要素（形状は異なる）は同じ符号を付している。

従来技術に関して上記したように、一般には、イオン及び中性ガスの種の数が多くなると、解くべき方程式の数が非常に多くなり、正確に解析することは困難である。しかし、通常のＣＶＤ条件、特に本発明の対象であるダイヤモンドの合成条件においては、電子密度の空間変化に対して、イオン及び中性ガスの空間分布の時間変化は比較的緩やかであり、これらの非一様性の電子の密度、流速及び温度、並びに電磁界分布へ影響は少ない。従って、イオン及び中性ガスに関して一様性を仮定すれば、イオン及び中性ガスについて解く必要が無く、電子に関する変数のみを考慮すればよい。更に、以下に示すように、電子の流速及び温度を決定する方程式も解析的に解くことができる。

まず、電子に関する物理量を高周波振動成分と平均速度成分に分離することによって、電子の運動量輸送方程式が解析的に解け、電子の流速を他の物理量で表すことができる。

電子の運動量輸送方程式は、次式１のように表せる。

式１において、ｎ_eは電子密度、ｕ_eは電子流速（３次元ベクトル）、Ｔ_eは電子温度、
Ｅは電界（３次元ベクトル）、ｔは時間、ν_eは電子の衝突振動数、ｍ_eは電子質量である。

ここで、電子流速ｕ_eを平均速度成分ｕ_e0と高周波成分ｕ_e1との和（ｕ_e＝ｕ_e0＋ｕ_e1）で表すと、式２、式３が導出できる。ここで、通常、電子温度の勾配は電子密度勾配と比較して緩やかであることから、電子温度に起因する微分係数は無視した。

式２、３において、ωは電磁界角周波数、ｅは電荷素量である。

次に、上記非特許文献１では、定常状態における電子の質量輸送方程式を、式４で表されるとしている。

式４の右辺第１項は電子と化学種（中性原子、中性分子、イオン等）との電離反応による電子の生成を表し、右辺第２項は電子と化学種との再結合反応による電子の消滅を表す。

式４は、非特許文献１で考慮されている化学反応のみに特化した場合の式であるが、一般には、生成項及び消滅項は式５、式６のように、可能な反応の総和として表される。ここで、可能な反応には、種類の異なる反応に限らず、同じ種類の反応であっても反応のエネルギー閾値（以下、単にエネルギー閾値とも記す）や反応率が異なる反応をも含む。即ち、例えば励起反応であっても、反応に関わる原子・分子種が異なる励起反応であれば、一般にはエネルギー閾値や反応率が異なるので、添え字ｓに一意の値を割り当てて別の反応として扱う。

式５、式６において、sは反応の種類を示す添字であり、ｎ_sは反応種の密度であり、

は、それぞれ電離反応の反応率、電離反応のエネルギー閾値、再結合反応の反応率、再結合反応のエネルギー閾値を表す。消滅項については、電荷準中性条件により、電子密度とイオン密度は局所的にほぼ等しいことから、式７を用いている。

ここで、

はイオン種ｓのモル分率を表す。

非特許文献１では、電離反応及び再結合反応を、ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺とこの逆反応（
２ｅ＋Ｈ₂ ⁺→ｅ＋Ｈ₂）のみに限定し、再結合反応のエネルギー閾値を０としているため
、式４の様になっている。式４の右辺を一般形である式５及び式６で置き換えると、式８を得る。

この式は高周波振動成分を含まない平均化された質量輸送の釣り合いを示しているということが重要である。

以上によって得られた式８の左辺のｕ_eを式３の静的な平均成分ｕ_e0で置き換え、以下
の式９を得る。

次に、以上で決定した方程式を用いて数値解析する方法を、図３に示したフローチャートに従って説明する。図３は、本発明の実施の形態に係るマイクロ波プラズマ解析プログラムを示すフローチャートである。

以下に説明する処理は、コンピュータを用いて微分方程式を数値的に積分する通常の処理に、多くの点で類似する処理である。また、以下において、特に断らない限り、コンピュータの演算処理手段（以下、ＣＰＵと記す）が行う処理として記載する。ＣＰＵは、メモリをワーク領域として使用して演算などの処理を行う。ＣＰＵが、処理を行う前提条件は、外部から人によって操作手段（コンピュータ用キーボード、マウスなど）が操作されて、データとして入力される。入力されたデータは、メモリの所定領域に一時記憶され、演算処理に使用される。入力されたデータや演算処理の結果は、適宜記録部（ハードディ
スクなど）に記録される。また、以下においては、図１に示した構成のＣＶＤ装置を対象として説明する。

ステップＳ１において、数値計算の対象とする空間、及びその空間内に存在する物の形状に関する情報の入力を受け付ける。具体的には、人がコンピュータの操作手段を操作して、ＣＰＵが提供するグラフィカルユーザインタフェースを介してモニタ画面上で、図１に示した形状に関する情報を入力する。人が、各構成部（第１円筒部、第２円筒部、石英板、サセプタなど）の形状、寸法を指定し、対象とする空間中に各構成部品を配置することで、計算対象となる空間を区切る境界の位置、形状、寸法に関するデータが、メモリ上に生成される。生成されたデータは適宜、記録部（ハードディスク）に保存される。

ステップＳ２において、ステップＳ１で決定した形状をメッシュに分割する情報の入力を受け付ける。例えば、人が操作手段を操作して、各構成部の各辺を等間隔に分割する数を指定することによって、対向する分割点が線分で結ばれ、各構成部がメッシュに分割される。メッシュ分割に関する情報の入力は、各構成部のメッシュサイズを指定することで行うこともできる。

ステップＳ３において、境界条件の入力を受け付ける。具体的には、金属部分（図１では、第１円筒部、第２円筒部、ウェーブガイド、アンテナ、サセプタ）を完全導体とし、その表面に沿う方向の電界成分を０とする。また、電磁波が入力される部分（マイクロ波入力部）及び放電領域を指定する。

ステップＳ４において、計算の初期条件の入力を受け付ける。各構成部の比誘電率、電気伝導率、初期の電子密度ｎ_e ⁽ⁱⁿⁱ⁾の空間分布、入力される電磁波の周波数、投入電力、電磁波モードが入力される。また、考慮する各反応ｓの反応率ｋ_s及びエネルギー閾値Ｅ_s、並びに反応種の密度ｎ_sが入力される。また、背景ガスの圧力及び温度が入力される。
ここで、背景ガスとはプラズマ源のガスを意味し、背景ガスの圧力及び温度は、ν_eを評
価する時に使用される。その評価式は、Ｐgをガス圧、Ｔgをガス温度として、ν_ｅ＝１．４４×１０¹²×Ｐg／Ｔgであり、これは文献W. Tan and T. A. Grotjohn, Diamond and Related Materials 4, 1145 (1995)等で用いられている評価式である。さらに、繰り返し
処理の収束性を判断するための基準パラメータｎ_thが入力される。

なお、初期の電子密度ｎ_e ⁽ⁱⁿⁱ⁾の空間分布に関しては、別途、全領域でｎ_e＝０と仮定
して電界強度分布を求め、真空領域に相当する部分で電界強度が最大となった点を中心に、所定の広がり及び最大値（例えば、１×１０⁺¹⁷）を持つ分布を、初期の電子密度ｎ_e ⁽ⁱⁿⁱ⁾の空間分布とする。このように初期の電子密度ｎ_e ⁽ⁱⁿⁱ⁾の空間分布を設定すれば、初
期の電子密度ｎ_e ⁽ⁱⁿⁱ⁾として、真空領域全体において一定値（例えば０）を設定する場合よりも、後述の計算を速く収束させることができる。

ここで、考慮する反応は、電離反応、励起反応、解離反応、再結合反応であるが、上記したように同じ種類の反応であってもエネルギー閾値や反応率が異なる場合の反応をも含み、添え字ｓに一意の値を割り当てて別の反応として扱う。

ステップＳ５において、電子密度ｎ_e（ｒ）から、ステップＳ３で指定した放電領域で
の誘電率ε_p（ｒ）を求める。ここで、誘電率ε_p（ｒ）は、電子密度ｎ_e（ｒ）を用いて

で求める。後述する通り、ステップＳ５〜Ｓ８は、所定の条件を満足するまで繰り返される。条件を満足しない場合には、ステップＳ６に戻って再度式１０から誘電率ε_p（ｒ）
を求めるが、この場合は後述するステップＳ８での処理により計算された電子密度ｎ_e（
ｒ）を用いる。ステップＳ９の判定を経ない、初めてステップＳ５において誘電率ε_p（
ｒ）を求める際は、ステップＳ４において指定した初期の電子密度ｎ_e ⁽ⁱⁿⁱ⁾の空間分布から誘電率ε_p（ｒ）を求める。

ステップＳ６において、ステップＳ１〜Ｓ４で指定された固定条件（形状、境界条件）及び、ステップＳ５において求められた誘電率ε_p（ｒ）の下で、マクスウェルの方程式
を用いて、電界分布Ｅ（ｒ）を数値計算する。ここでＥ（ｒ）は、３次元の位置ベクトルｒで指定される位置における電界の３次元ベクトルである。電磁気に関するマクスウェルの方程式、及び、有限要素法などを用いて所定の境界条件の下でマクスウェルの方程式を数値積分し、所定空間内の電界分布Ｅ（ｒ）を求める方法は公知であるので、ここでは説明を省略する。

ステップＳ７において、次の式１１で表される電子のエネルギー輸送方程式を解く。

具体的には、ステップＳ６で求めた電界分布Ｅ（ｒ）を用いて、式１１を２分法によって解き、電子温度Ｔeを求める。図４は、２分法によって電子温度Ｔeを求める処理を示すグラフである。２分法は公知であるので、ここでは詳細は省略し、簡単に説明する。図４において、横軸は電子温度Ｔeであり、縦軸は式１１の右辺（Ｒ(Ｔe)とする）である。水平の点線は式１１の左辺の値（Ｌとする）を表す。異なる２つのＴeの値Ｔe1、Ｔe2に対
して、それらの中点Ｔem＝（Ｔe1＋Ｔe2）／２を求め、中点Ｔemにおける式１１の左辺の値Ｒ(Ｔem)と、式１１の右辺Ｌとの大小を比較する。そして、Ｒ(Ｔem)＜Ｌであれば、Ｔemを新たにＴe1とし、Ｒ(Ｔem)＞Ｌであれば、Ｔemを新たにＴe2とし、Ｒ(Ｔem)とＬとの差が所定値ε以下になるまで繰り返す。従って、所定値εによって解の精度を指定することができる。２分法を適用する場合のＴeの初期値は任意に設定することができるが、効
率的に計算するには、Ｔeの２つの初期値を０ｅＶ以上１０ｅＶ以下の範囲とすることが
望ましい。

ステップＳ８において、ステップＳ６で求められた電界Ｅ（ｒ）（各メッシュにおける３次元ベクトル）及びステップＳ７で求めた電子温度Ｔ_eを用いて、ステップＳ１〜Ｓ４
で指定された固定条件（形状、境界条件）の下で、上記で導出した式９の微分方程式を用いて、電子密度ｎ_e（ｒ）を数値計算によって求める。

ここで、各メッシュについて数値計算を行う前に、そのメッシュの電界強度Ｅ（ｒ）が０であるか否かを判断し、Ｅ（ｒ）＝０のメッシュについては数値計算をせずにｎ_e＝０
とし、Ｅ（ｒ）≠０のメッシュについてのみ数値計算を行う。なお、有限要素法などを用いて所定の境界条件の下で式９の微分方程式を数値積分する方法は公知であるので、ここでは説明を省略する。

ステップＳ９において、電子密度分布の計算結果の収束性を判断する。具体的には、ｋ回目、ｋ−１回目の電子密度分布の計算結果をｎ_e ^k、ｎ_e ^k-1で表し、ステップＳ４で入力された基準パラメータｎ_thを用いて、

を満たせば、収束したと判断してステップＳ１０に移行し、満たしていなければステップＳ５に戻る。ここで、τの積分は放電領域での体積積分を表す。ステップＳ５に戻った場合、誘電率ε_p（ｒ）を求める際は、ステップＳ８で計算された電子密度ｎ_e（ｒ）を用いる。

ステップＳ１０において、最新の計算結果である電子密度分布ｎ_e（ｒ）を、例えばモ
ニタ画面にコンピュータグラフィック（以下、ＣＧと記す）によって描画することによって、提示する。

最後にステップＳ１１において、終了の指示があったか否かを判断し、終了の指示があるまで、ステップＳ１〜Ｓ１０の処理を繰り返す。

以上によって、ステップＳ１０で提示されたＣＧ画像を、人が目視によって評価し、ＣＶＤ装置の形状及び寸法設計の妥当性を評価することができる。

また、ステップＳ８における数値計算は、１つの方程式を用いた処理でよいので、従来よりも高速にマイクロ波プラズマの分布を計算することが可能である。

上記した一連の処理を行うコンピュータプログラムを構成する場合、図５に示したように、誘電率を計算（ステップＳ５）するモジュール、マクスウェル方程式を用いて電界を数値計算（ステップＳ６）するモジュール、電子温度を計算（ステップＳ７）するモジュール、及び式９を用いて電子密度を数値計算（ステップＳ８）するモジュールを、それぞれ別のモジュールに構成することができる。この様に構成することによって、電子密度、電子温度及び誘電率を計算するモジュールだけを新たにプログラミングすれば、それ以外の処理に関しては、公知の電磁界解析用ソフトウェアを使用することができる。

また、上記では電子のエネルギー輸送方程式を式１１としたが、次式のように、Ｔeに
関する依存性がこれと異なる原料ガス系（Ａｒ系など）も存在する。

ここで、βsは、原料ガスに応じて適宜決定されるパラメータである。βs＝０であれば、式１２は式１１と等しくなる。βs≠０である場合には、非特許文献３に開示された試験
関数を用いてパラメータフィッティングすると、Ｔeへの依存性を試験関数で十分に表す
ことができないために、誤差が大きくなる。

また、Ｔeへの依存性がさらに複雑な場合（電離・励起・再結合等の各種反応をＸで表
し、一般的には式１３で表される）には、非特許文献３に開示された試験関数を用いてパラメータフィッティングを用いる方法では、より誤差が大きくなると予想される。

しかし、これらの式についても、２分法、反復法、変分法などの方法による数値計算を行うことで、所望の精度でＴeを求めることができる。反復法、変分法は公知であるので
、ここでは説明を省略する。反復法、変分法を用いる場合にも、効率的に計算するには、２分法を用いる場合と同様に初期値を設定することが望ましい。即ち、反復法の場合には、電子温度Ｔeの初期値を０ｅＶ以上１０ｅＶ以下の範囲とし、変分法の場合には、電子
温度Ｔeが０ｅＶ以上１０ｅＶ以下の範囲で解を探索することが望ましい。

以下に、実施例を示し、本発明の特徴をより明確にする。本実施例では、図１に示すＣＶＤ装置を想定し、本発明のマイクロ波プラズマ解析プログラムによる数値計算を行った。

図６は、各部の寸法を付した断面図である。メッシュサイズ、即ち、１つの正方形メッシュの１辺の長さは、サセプタ天面中央上方の、直径４ｃｍ高さ３ｃｍの円柱領域（図６の点線で囲んだ領域）は２ｍｍ、この領域以外は１ｃｍに設定した。

境界条件に関しては、黒色の構成部（第１円筒部、第２円筒部、アンテナ、サセプタ、ウェーブガイド）は全て完全導体とし、表面に沿う電界成分を０とし、電磁波はマイクロ波入力部から導入されるように設定した。

放電領域は、石英板の下側の第２円筒部内である。マイクロ波入力部から入力される電磁波は、周波数が２．４５ＧＨｚ、投入電力が１６００Ｗ、電磁波モードが最低次のＴＥモードとした。

初期の電子密度分布は、円盤状のサセプタの中心からの距離をｒ（ｍ）とし、サセプタ
上方の放電領域において、ｒ＜０．０２である領域では、ｎ_e ⁽ⁱⁿⁱ⁾＝２．５×１０¹⁷×[
１．０＋ｃｏｓ（πｒ／０．０２）]であり、ｒ≧０．０２の領域では、ｎ_e ⁽ⁱⁿⁱ⁾＝０で
ある分布とした。

各構成部の固定の誘電率及び電気伝導率については、石英板が比誘電率＝３．８、電気伝導率＝１ｅ^-17ｍ／Ωとし、空気領域は比誘電率＝３．８、電気伝導率＝０ｍ／Ωであ
るとした。

プラズに関わるパラメータに関しては、表１の４種類の反応を想定し、電子密度分布の収束条件として、ｎ_th＝０．０１とした。

また、背景ガスの圧力は１６０Torr、温度は１０³Ｋとした。また、簡単のため、反応
種密度ｎ_sは、圧力１６０Torr、温度１０³Ｋからｎ_s＝１．５×１０²⁴とした。

以上の条件で、上記したステップＳ５〜Ｓ９の処理を実行した結果を図７に示す。図７の左半分には、図６の破線領域の電子密度を、右半分には図６の破線領域のｐ_abを示している。何れも、明度が高いほど（白いほど）値が大きい。また、図８は、プラズマ発生領域の電界強度の分布を示す図である。図７、図８には、各値が最大となる場所を黒丸で示し、その最大値を数値で記載している。

図７、図８から、本シミュレーションによって計算されたパワー密度が、実験から見積もられるパワー密度に一致していることが分かった。また、電子密度が最大となる位置（図７参照）と電界強度が最大となる位置（図８参照）とが異なっており、上記非特許文献１では再現できなかった、高密度プラズマによって電磁界が遮蔽される効果を、本実施例で再現することができた。このように、本実施例によって得られた分布の傾向は、実際の放電の様子や、他のより大規模なシミュレーションを行った結果と一致していた。

本発明が解析対象とするＣＶＤ装置のプラズマ発生部の一例を示す断面図である。 本発明が解析対象とするＣＶＤ装置のプラズマ発生部の別の例を示す断面図である。 本発明の実施の形態に係るマイクロ波プラズマ解析プログラムを示すフローチャートである。 本発明の実施の形態に係るマイクロ波プラズマ解析プログラムにおける２分法を説明するための図である。 本発明の実施の形態に係るマイクロ波プラズマ解析プログラムの構成の一例を示すブロック図である。 本発明の実施例を示すＣＶＤ装置のプラズマ発生部の構造を示す断面図である。 数値計算の結果を示す図である。 数値計算の結果を示す図である。

符号の説明

１ マイクロ波入力部
２ ウェーブガイド
３ 第1円筒部
４ 石英板
５ 第２円筒部
６ 真空領域
７ アンテナ
８ サセプタ

Claims (9)

1. ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値計算するコンピュータ用プログラムであって、
   コンピュータに、
   数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
   電子密度分布から誘電率分布を求める第２の機能と、
   前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第３の機能と、
   反応種をｓで表し、電荷素量をｅ、電子質量をｍ_e、電子の衝突振動数をν_e、電子温度をＴ_e、電磁界角周波数をω、電界をＥ、反応種ｓの密度をｎ_s、並びに、電離・励起・再結合等の各種反応をＸで表し、反応Ｘのエネルギー閾値と反応率とをそれぞれｋ_s ^(X)及びＥ_s ^(X)として、前記電界分布の下で、
   の方程式を数値計算によって解き、電子温度を求める第４の機能と、
   前記電界分布、前記形状情報、前記メッシュ情報、前記電子温度、前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、反応種ｓのモル分率を
   として、
   の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第５の機能とを実現させ、
   前記初期の電子密度分布が、数値計算の対象とする空間の各場所での電子密度を０と仮定した電界強度分布において、真空領域に相当する部分で電界強度が最大となる点で最大値を持つ分布であり、
   前記第５の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第２の機能〜前記第５の機能を順に繰り返して実現させることを特徴とするマイクロ波プラズマ解析プログラム。
2. 前記第５の機能が、前記電界分布において、電界強度が０以外の値であるメッシュに関してのみ、新たに電子密度を求め、電界強度が０のメッシュに関しては、電子密度を０とする機能であることを特徴とする請求項１に記載のマイクロ波プラズマ解析プログラム。
3. 前記第４の機能で使用する方程式が、
   であることを特徴とする請求項１又は２に記載のマイクロ波プラズマ解析プログラム。
4. 前記第４の機能で使用する方程式が、
   であり、
   前記第４の機能で使用する数値計算が２分法であり、該２分法の初期値である電子温度を０ｅＶ以上１０ｅＶ以下とすることを特徴とする請求項３に記載のマイクロ波プラズマ解析プログラム。
5. 前記第４の機能で使用する方程式が、
   であり、
   前記第４の機能で使用する数値計算が反復法であり、該反復法の初期値である電子温度を０ｅＶ以上１０ｅＶ以下とすることを特徴とする請求項３に記載のマイクロ波プラズマ解析プログラム。
6. 前記第４の機能で使用する方程式が、
   であり、
   前記第４の機能で使用する数値計算が変分法であり、電子温度が０ｅＶ以上１０ｅＶ以下の範囲で解を探索することを特徴とする請求項３に記載のマイクロ波プラズマ解析プログラム。
7. 前記反応が、ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺、ｅ＋Ｈ₂→ｅ＋Ｈ₂ ^*、ｅ＋Ｈ₂→ｅ＋２Ｈ、及びＨ₃ ⁺＋ｅ→３Ｈであり、
   反応ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺の反応率が１．０×１０^-14ｍ³／sec、エネルギー閾値が１５．４ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
   反応ｅ＋Ｈ₂→ｅ＋Ｈ₂ ^*の反応率が６．５×１０^-15ｍ³／sec、エネルギー閾値が１２．０ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
   反応ｅ＋Ｈ₂→ｅ＋２Ｈの反応率が１．０×１０^-14ｍ³／sec、エネルギー閾値が１０．０ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
   反応Ｈ₃ ⁺＋ｅ→３Ｈの反応率が１．０×１０^-13ｍ³／sec、エネルギー閾値が０．０ｅ
   Ｖ、反応種の密度が１．５×１０²⁴ｍ^-3であることを特徴とする請求項１〜６の何れか１項に記載のマイクロ波プラズマ解析プログラム。
8. 前記形状情報が円盤状のサセプタの形状情報を含み、
   前記サセプタの円形表面の中心からの距離をｒｍとし、前記サセプタ外部のｒ＜０．０２である領域内で、初期電子密度が２．５×１０¹⁷×[１．０＋ｃｏｓ（πｒ／０．０２
   ）]であり、前記サセプタ外部のｒ≧０．０２である領域で初期電子密度が０であること
   を特徴とする請求項１〜７の何れか１項に記載のマイクロ波プラズマ解析プログラム。
9. 初期の電子密度分布が、数値計算の対象とする空間の各場所での電子密度を０と仮定した電界強度分布において、真空領域に相当する部分で電界強度が最大となる点で最大値を持つ分布であることを特徴とする、マイクロ波プラズマ解析プログラム。