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

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

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

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

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

本願発明者は、上記課題を解決するために鋭意研究を重ねた結果、数値的に解くべきプラズマの方程式を一つにできることを見出した。即ち、物理量を高周波振動成分と平均速度成分に分離し、パラメータフィッティングを併用することによって、真空の電磁界の計算と同程度の計算量（数値的に解くべき式が一つしか増えない）で、マイクロ波プラズマを考慮した数値解析が可能となる。本発明の目的は、以下の手段によって達成される。

即ち、本発明に係るマイクロ波プラズマ解析プログラム（１）は、
ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値計算するコンピュータ用プログラムであって、
コンピュータに、
数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、パラメータａ_L、ｐ_L及びＥ_L、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの
圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
前記電子密度分布から誘電率分布を求める第２の機能と、
前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第３の機能と、
前記電界分布、前記形状情報、前記メッシュ情報、前記パラメータａ_L、ｐ_L及びＥ_L、
前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、
及び
として、
の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第４の機能とを実現させ、
前記第４の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第２の機能〜前記第４の機能を順に繰り返して実現させることを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（２）は、上記のマイクロ波プラズマ解析プログラム（１）において、
前記パラメータａ_L、ｐ_L及びＥ_Lがそれぞれ、ａ_L＝０．３６７５７、ｐ_L＝２．４６４
×１０¹⁵（ｅＶ／ｓｅｃ）、Ｅ_L＝１８．５３２７（ｅＶ^aL）であることを特徴としてい
る。

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

また、本発明に係るマイクロ波プラズマ解析プログラム（４）は、
ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値計算するコンピュータ用プログラムであって、
コンピュータに、
数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
前記反応率をｋ_sとし、前記エネルギー閾値をＥ_sとし、且つ前記反応種の密度をｎ_sと
し、電子温度Ｔ_eを所定範囲で変化させて
を計算して得られた前記電子温度Ｔ_eに関するグラフに、
がフィッティングするようにパラメータａ_L、ｐ_L及びＥ_Lを決定する第２の機能と、
前記電子密度分布から誘電率分布を求める第３の機能と、
前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第４の機能と、
前記電界分布、前記形状情報、前記メッシュ情報、前記パラメータａ_L、ｐ_L及びＥ_L、
前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、
及び
として、
の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第５の機能とを実現させ、
前記第５の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第３の機能〜前記第５の機能を順に繰り返して実現させることを特徴としている。

また、本発明に係るマイクロ波プラズマ解析プログラム（５）は、上記のマイクロ波プラズマ解析プログラム（１）〜（４）の何れかにおいて、
前記反応が、ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺、ｅ＋Ｈ₂→ｅ＋Ｈ₂ ^*、ｅ＋Ｈ₂→ｅ＋２Ｈ、及びＨ₃ ⁺＋ｅ→３Ｈであり、
反応ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺の反応率が１．０×１０^-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が入力される。

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

ステップＳ５において、パラメータフィッティングを行う。即ち、ステップＳ４において入力された反応の反応率ｋ_s及びエネルギー閾値Ｅ_s、並びに反応種の密度ｎ_sを用いて
、電子の温度Ｔ_eを所定範囲で変化させて、
を計算し、得られたＴ_eに関するグラフ（図４参照）に、
がフィッティングするように、即ち、式１０のグラフを式１１で近似できるように３つのパラメータａ_L、ｐ_L、Ｅ_Lを決定する。

図４は、パラメータフィッティングを説明するためのグラフであり、「×」印は式１０で計算された値であり、曲線は、決定されたパラメータを用いた式１１の曲線である。これらのパラメータが決定されれば、式１２（電子のエネルギー輸送方程式）から、式１３を用いて電子の温度Ｔ_eは式１４によって表される。ここで｜Ｅ｜は電界強度を表す。

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

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

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

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

ステップＳ９において、電子密度分布の計算結果の収束性を判断する。具体的には、ｋ回目、ｋ−１回目の電子密度分布の計算結果をｎ_e ^k、ｎ_e ^k-1で表し、ステップＳ４で入力された基準パラメータｎ_thを用いて、
を満たせば、収束したと判断してステップＳ１１に移行し、満たしていなければステップＳ６に戻る。ここで、τの積分は放電領域での体積積分を表す。ステップＳ６に戻った場合、誘電率ε_p（ｒ）を求める際は、ステップＳ８で計算された電子密度ｎ_e（ｒ）を用いる。

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

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

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

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

以上では、ステップＳ５においてパラメータフィッティングを行う場合を説明したが、予めパラメータフィッティングを行い、パラメータａ_L、ｐ_L、Ｅ_Lを計算しておいてもよ
い。その場合には、ステップＳ４での初期値の入力において、予め求められたパラメータa_L、ｐ_L、Ｅ_Lを入力するようにすれば、ステップＳ５の処理が不要になる。特に、対象とする反応を特定種類の反応に固定してシミュレーションを繰り返す場合には、パラメータａ_L、ｐ_L、Ｅ_Lを一度だけ計算すればよいので、全体としての処理時間をより短縮するこ
とができる。

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

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

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

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

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

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

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

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

また、フィッティングパラメータには、予め計算した値である、Ｅ_L＝１８．５３２７
（ｅＶ^aL）、ｐ_L＝２．４６４×１０¹⁵（ｅＶ／ｓｅｃ）、ａ_L＝０．３６７５７を用いた。また、背景ガスの圧力は１６０Torr、温度は１０³Ｋとした。また、簡単のため、反応
種密度ｎ_sは、圧力１６０Torr、温度１０³Ｋからｎ_s＝１．５×１０²⁴とした。

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

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

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

符号の説明

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

Claims (6)

1. ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値計算するコンピュータ用プログラムであって、
   コンピュータに、
   数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、パラメータａ_L、ｐ_L及びＥ_L、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの
   圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
   前記電子密度分布から誘電率分布を求める第２の機能と、
   前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第３の機能と、
   前記電界分布、前記形状情報、前記メッシュ情報、前記パラメータａ_L、ｐ_L及びＥ_L、
   前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、
   及び
   として、
   の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第４の機能とを実現させ、
   前記第４の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第２の機能〜前記第４の機能を順に繰り返して実現させることを特徴とするマイクロ波プラズマ解析プログラム。
2. 前記パラメータａ_L、ｐ_L及びＥ_Lがそれぞれ、ａ_L＝０．３６７５７、ｐ_L＝２．４６４
   ×１０¹⁵（ｅＶ／ｓｅｃ）、Ｅ_L＝１８．５３２７（ｅＶ^aL）であることを特徴とする請
   求項１に記載のマイクロ波プラズマ解析プログラム。
3. 前記第４の機能が、前記電界分布において、電界強度が０以外の値であるメッシュに関してのみ、新たに電子密度を求め、電界強度が０のメッシュに関しては、電子密度を０とする機能であることを特徴とする請求項１又は２に記載のマイクロ波プラズマ解析プログラム。
4. ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数
   値計算するコンピュータ用プログラムであって、
   コンピュータに、
   数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
   前記反応率をｋ_sとし、前記エネルギー閾値をＥ_sとし、且つ前記反応種の密度をｎ_sと
   し、電子温度Ｔ_eを所定範囲で変化させて
   を計算して得られた前記電子温度Ｔ_eに関するグラフに、
   がフィッティングするようにパラメータａ_L、ｐ_L及びＥ_Lを決定する第２の機能と、
   前記電子密度分布から誘電率分布を求める第３の機能と、
   前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第４の機能と、
   前記電界分布、前記形状情報、前記メッシュ情報、前記パラメータａ_L、ｐ_L及びＥ_L、
   前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、
   及び
   として、
   の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第５の機能とを実現させ、
   前記第５の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第３の機能〜前記第５の機能を順に繰り返して実現させることを特徴とするマイクロ波プラズマ解析プログラム。
5. 前記反応が、ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺、ｅ＋Ｈ₂→ｅ＋Ｈ₂ ^*、ｅ＋Ｈ₂→ｅ＋２Ｈ、及びＨ₃ ⁺＋ｅ→３Ｈであり、
   反応ｅ＋Ｈ₂→２ｅ＋Ｈ₂ ⁺の反応率が１．０×１０^-14ｍ³／sec、エネルギー閾値が１５．４ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
   反応ｅ＋Ｈ₂→ｅ＋Ｈ₂ ^*の反応率が６．５×１０^-15ｍ³／sec、エネルギー閾値が１２．０ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
   反応ｅ＋Ｈ₂→ｅ＋２Ｈの反応率が１．０×１０^-14ｍ³／sec、エネルギー閾値が１０．０ｅＶ、反応種の密度が１．５×１０²⁴ｍ^-3であり、
   反応Ｈ₃ ⁺＋ｅ→３Ｈの反応率が１．０×１０^-13ｍ³／sec、エネルギー閾値が０．０ｅ
   Ｖ、反応種の密度が１．５×１０²⁴ｍ^-3であることを特徴とする請求項１〜４の何れか１項に記載のマイクロ波プラズマ解析プログラム。
6. 前記形状情報が円盤状のサセプタの形状情報を含み、
   前記サセプタの円形表面の中心からの距離をｒｍとし、前記サセプタ外部のｒ＜０．０２である領域内で、初期電子密度が２．５×１０¹⁷×[１．０＋ｃｏｓ（πｒ／０．０２
   ）]であり、前記サセプタ外部のｒ≧０．０２である領域で初期電子密度が０であること
   を特徴とする請求項５に記載のマイクロ波プラズマ解析プログラム。