(57)【要約】【課題】マイクロ波プラズマ分布を数値解析するプログラムを提供すること。【解決手段】計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、フィッティングパラメータ、反応毎のパラメータ、背景ガスの圧量及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付け（Ｓ１〜Ｓ４）、電子密度分布から誘電率分布を求め（Ｓ６）、形状情報、メッシュ情報、境界条件、誘電率分布、並びに、電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いて電界分布を求め（Ｓ７）、電界分布、形状情報、メッシュ情報、フィッティングパラメータ、反応毎のパラメータ、背景ガスの圧力及び温度、並びに、電磁波の周波数の下で、一つの電子輸送方程式を用いた数値計算によって電子密度分布を求め（Ｓ８）、電子密度分布の変化が収束するまで処理（Ｓ６〜Ｓ９）を繰り返す。【選択図】図３
2008031515.tif 000002
【特許請求の範囲】
【請求項１】
ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値計算するコンピュータ用プログラムであって、
コンピュータに、
数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、パラメータａL、ｐL及びＥL、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの
圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
前記電子密度分布から誘電率分布を求める第２の機能と、
前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第３の機能と、
前記電界分布、前記形状情報、前記メッシュ情報、前記パラメータａL、ｐL及びＥL、
前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、
【数１】
2008031515.tif 000027
及び
【数２】
2008031515.tif 000028
として、
【数３】
2008031515.tif 000029
の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第４の機能とを実現させ、
前記第４の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第２の機能〜前記第４の機能を順に繰り返して実現させることを特徴とするマイクロ波プラズマ解析プログラム。
【請求項２】
前記パラメータａL、ｐL及びＥLがそれぞれ、ａL＝０．３６７５７、ｐL＝２．４６４
×１０15（ｅＶ／ｓｅｃ）、ＥL＝１８．５３２７（ｅＶaL）であることを特徴とする請
求項１に記載のマイクロ波プラズマ解析プログラム。
【請求項３】
前記第４の機能が、前記電界分布において、電界強度が０以外の値であるメッシュに関してのみ、新たに電子密度を求め、電界強度が０のメッシュに関しては、電子密度を０とする機能であることを特徴とする請求項１又は２に記載のマイクロ波プラズマ解析プログラム。
【請求項４】
ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数
値計算するコンピュータ用プログラムであって、
コンピュータに、
数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
前記反応率をｋsとし、前記エネルギー閾値をＥsとし、且つ前記反応種の密度をｎsと
し、電子温度Ｔeを所定範囲で変化させて
【数４】
2008031515.tif 000030
を計算して得られた前記電子温度Ｔeに関するグラフに、
【数５】
2008031515.tif 000031
がフィッティングするようにパラメータａL、ｐL及びＥLを決定する第２の機能と、
前記電子密度分布から誘電率分布を求める第３の機能と、
前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第４の機能と、
前記電界分布、前記形状情報、前記メッシュ情報、前記パラメータａL、ｐL及びＥL、
前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、
【数６】
2008031515.tif 000032
及び
【数７】
2008031515.tif 000033
として、
【数８】
2008031515.tif 000034
の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第５の機能とを実現させ、
前記第５の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第３の機能〜前記第５の機能を順に繰り返して実現させることを特徴とするマイクロ波プラズマ解析プログラム。
【請求項５】
前記反応が、ｅ＋Ｈ2→２ｅ＋Ｈ2+、ｅ＋Ｈ2→ｅ＋Ｈ2*、ｅ＋Ｈ2→ｅ＋２Ｈ、及びＨ3+＋ｅ→３Ｈであり、
反応ｅ＋Ｈ2→２ｅ＋Ｈ2+の反応率が１．０×１０-14ｍ3／sec、エネルギー閾値が１５．４ｅＶ、反応種の密度が１．５×１０24ｍ-3であり、
反応ｅ＋Ｈ2→ｅ＋Ｈ2*の反応率が６．５×１０-15ｍ3／sec、エネルギー閾値が１２．０ｅＶ、反応種の密度が１．５×１０24ｍ-3であり、
反応ｅ＋Ｈ2→ｅ＋２Ｈの反応率が１．０×１０-14ｍ3／sec、エネルギー閾値が１０．０ｅＶ、反応種の密度が１．５×１０24ｍ-3であり、
反応Ｈ3+＋ｅ→３Ｈの反応率が１．０×１０-13ｍ3／sec、エネルギー閾値が０．０ｅ
Ｖ、反応種の密度が１．５×１０24ｍ-3であることを特徴とする請求項１〜４の何れか１項に記載のマイクロ波プラズマ解析プログラム。
【請求項６】
前記形状情報が円盤状のサセプタの形状情報を含み、
前記サセプタの円形表面の中心からの距離をｒｍとし、前記サセプタ外部のｒ＜０．０２である領域内で、初期電子密度が２．５×１０17×[１．０＋ｃｏｓ（πｒ／０．０２
）]であり、前記サセプタ外部のｒ≧０．０２である領域で初期電子密度が０であること
を特徴とする請求項５に記載のマイクロ波プラズマ解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、マイクロ波プラズマを解析するプログラムに関し、特にダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値解析することができるマイクロ波プラズマ解析プログラムに関する。
【背景技術】
【０００２】
マイクロ波プラズマを用いたＣＶＤ装置の設計、合成条件の割り出し、合成速度の見積もり等には、プラズマの中の電子密度、電子温度やパワー密度等の情報が必要である。しかし、実験によりマイクロ波プラズマを観測することは困難であり、そのため、数値計算に頼っているのが現状である。
【０００３】
真空の電磁界を計算する方法、及びそのためのコンピュータソフトウェアは公知であり、複数のコンピュータソフトウェアが市販されている。しかし、マイクロ波プラズマを考慮した電磁界の解析機能を持つ、市販のコンピュータソフトウェアは無い。
【０００４】
マイクロ波プラプラズマでは、高周波電磁界の下で、電子、イオン及び中性ガスが混在している。従って、変数となる物理量には、次の(i)〜(iii)がある。
(i)電界及び磁界（２つの３次元ベクトル変数）
(ii)電子の密度、速度及び温度（２つのスカラー変数と１つの３次元ベクトル変数）
(iii)イオン及び中性ガスの、密度、速度及び温度（（２つのスカラー変数と１つの３次
元ベクトル変数）×種の数）
なお、(iii)に関しては、例えば水素の放電であれば、Ｈ、Ｈ2、Ｈ+、Ｈ2+、Ｈ2*など
のそれぞれの種について変数を導入する必要がある。本明細書において、上付き添字の＊は原子、分子の励起状態を表す。
【０００５】
一般には、これら全ての物理量に対応する支配方程式を解くことにより、各物理量の空間分布を決定することができる。従って、イオン及び中性ガスの種として、それぞれ一つだけを考慮するとしても、それらの定常状態における分布を求めるには、合計２１（＝６＋５＋５＋５）の偏微分方程式を連立して解かなくてはならない。
【０００６】
このように、マイクロ波プラズマの支配方程式は膨大且つ複雑であり、それを用いた数値計算によるシミュレーションは非常に困難である。即ち、数値的に解くべき微分方程式の数が非常に多い、数値計算における解の収束性が悪い、計算に時間が掛かる（計算コストが高い）、設定すべき未知のパラメータ量が膨大であるなどの問題がある。そして、これらのために、既存のソフトウェアでは、精度の良い計算結果を得ることができなかった。
【０００７】
一方、下記の非特許文献１及び非特許文献２には、ＣＶＤ装置におけるマイクロ波プラズマの電磁界分布を解析する方法が開示されている。
【非特許文献１】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、
前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、
【数１】
2008031515.tif 000003
及び
【数２】
2008031515.tif 000004
として、
【数３】
2008031515.tif 000005
の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第４の機能とを実現させ、
前記第４の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第２の機能〜前記第４の機能を順に繰り返して実現させることを特徴としている。
【００１３】
また、本発明に係るマイクロ波プラズマ解析プログラム（２）は、上記のマイクロ波プラズマ解析プログラム（１）において、
前記パラメータａL、ｐL及びＥLがそれぞれ、ａL＝０．３６７５７、ｐL＝２．４６４
×１０15（ｅＶ／ｓｅｃ）、ＥL＝１８．５３２７（ｅＶaL）であることを特徴としてい
る。
【００１４】
また、本発明に係るマイクロ波プラズマ解析プログラム（３）は、上記のマイクロ波プラズマ解析プログラム（１）又は（２）において、
前記第４の機能が、前記電界分布において、電界強度が０以外の値であるメッシュに関してのみ、新たに電子密度を求め、電界強度が０のメッシュに関しては、電子密度を０とする機能であることを特徴としている。
【００１５】
また、本発明に係るマイクロ波プラズマ解析プログラム（４）は、
ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を数値計算するコンピュータ用プログラムであって、
コンピュータに、
数値計算領域の形状情報、誘電率、電気伝導率、メッシュ情報、及び境界条件、反応毎の反応率、エネルギー閾値及び反応種の密度、背景ガスの圧力及び温度、入力される電磁波の周波数、投入電力及び電磁波モード、並びに、初期の電子密度分布の入力を受け付ける第1の機能と、
前記反応率をｋsとし、前記エネルギー閾値をＥsとし、且つ前記反応種の密度をｎsと
し、電子温度Ｔeを所定範囲で変化させて
【数４】
2008031515.tif 000006
を計算して得られた前記電子温度Ｔeに関するグラフに、
【数５】
2008031515.tif 000007
がフィッティングするようにパラメータａL、ｐL及びＥLを決定する第２の機能と、
前記電子密度分布から誘電率分布を求める第３の機能と、
前記形状情報、前記メッシュ情報、前記境界条件、前記誘電率分布、並びに、前記電磁波の周波数、投入電力及び電磁波モードの下で、マクスウェル方程式を用いた数値計算によって電界分布を求める第４の機能と、
前記電界分布、前記形状情報、前記メッシュ情報、前記パラメータａL、ｐL及びＥL、
前記反応率、前記エネルギー閾値、前記反応種の密度、前記背景ガスの圧力及び温度、並びに、前記電磁波の周波数の下で、
【数６】
2008031515.tif 000008
及び
【数７】
2008031515.tif 000009
として、
【数８】
2008031515.tif 000010
の方程式を用いた数値計算によって新たに電子密度分布を求めて前記電子密度分布とする第５の機能とを実現させ、
前記第５の機能によって求められた電子密度分布の変化が収束するまで、前記コンピュータに前記第３の機能〜前記第５の機能を順に繰り返して実現させることを特徴としている。
【００１６】
また、本発明に係るマイクロ波プラズマ解析プログラム（５）は、上記のマイクロ波プラズマ解析プログラム（１）〜（４）の何れかにおいて、
前記反応が、ｅ＋Ｈ2→２ｅ＋Ｈ2+、ｅ＋Ｈ2→ｅ＋Ｈ2*、ｅ＋Ｈ2→ｅ＋２Ｈ、及びＨ3+＋ｅ→３Ｈであり、
反応ｅ＋Ｈ2→２ｅ＋Ｈ2+の反応率が１．０×１０-14ｍ3／sec、エネルギー閾値が１５．４ｅＶ、反応種の密度が１．５×１０24ｍ-3であり、
反応ｅ＋Ｈ2→ｅ＋Ｈ2*の反応率が６．５×１０-15ｍ3／sec、エネルギー閾値が１２．０ｅＶ、反応種の密度が１．５×１０24ｍ-3であり、
反応ｅ＋Ｈ2→ｅ＋２Ｈの反応率が１．０×１０-14ｍ3／sec、エネルギー閾値が１０．０ｅＶ、反応種の密度が１．５×１０24ｍ-3であり、
反応Ｈ3+＋ｅ→３Ｈの反応率が１．０×１０-13ｍ3／sec、エネルギー閾値が０．０ｅ
Ｖ、反応種の密度が１．５×１０24ｍ-3であることを特徴としている。
【００１７】
また、本発明に係るマイクロ波プラズマ解析プログラム（６）は、上記のマイクロ波プラズマ解析プログラム（５）において、
前記形状情報が円盤状のサセプタの形状情報を含み、
前記サセプタの円形表面の中心からの距離をｒｍとし、前記サセプタ外部のｒ＜０．０
２である領域内で、初期電子密度が２．５×１０17×[１．０＋ｃｏｓ（πｒ／０．０２
）]であり、前記サセプタ外部のｒ≧０．０２である領域で初期電子密度が０であること
を特徴としている。
【発明の効果】
【００１８】
本発明のマイクロ波プラズマ解析プログラムによれば、ダイヤモンド合成に使用されるＣＶＤ装置内で発生するマイクロ波プラズマの分布を、従来よりも短時間で、且つ精度良く数値解析することができる。これによって、ダイヤモンド合成に適したＣＶＤ装置を効率的に設計することが可能となる。
【発明を実施するための最良の形態】
【００１９】
以下、本発明に係る実施の形態を、添付した図面に基づいて説明する。
【００２０】
まず、本発明の前提条件及び解くべき方程式に関して説明する。本発明が解析対象とするのは、ダイヤモンド合成のためのマイクロ波プラズマ放電である。９割程度が水素によって構成される原料ガスを比較的高いガス圧力（数Torr以上）で使用する。マイクロ波プラズマ発生装置には、図１、図２に示すような断面図の構造を有するＣＶＤ装置を想定している。これらのＣＶＤ装置は、サセプタの天面中央で放電が起こることが想定されており、各部の寸法、形状に応じて、形成されるマイクロ波プラズマ分布を解析することが要求される。
【００２１】
図１において、マイクロ波発生部（図示せず）で発生されマイクロ波入力部１から導入されたマイクロ波は、ウェーブガイド２を通ってアンテナ７を介して第１円筒部３及び第２円筒部５で定在波を形成する。石英板４を挟んで第１円筒部と隔離された第２円筒部５内の真空領域６には、プラズマ源となるガスが、ガス供給部（図示せず）から所定の圧力で供給されており、サセプタ８の天面中央で放電が起こり、真空領域６にプラズマが発生する。そして、サセプタ８に載置された物（図示せず）の表面にダイヤモンドを積層させる。また、図２は、図１と構造は異なるが、図１と同様にサセプタ８の天面中央で放電が起こり、真空領域中にプラズマが発生する。図２において、図１と同じ機能をする構成要素（形状は異なる）は同じ符号を付している。
【００２２】
従来技術に関して上記したように、一般には、イオン及び中性ガスの種の数が多くなると、解くべき方程式の数が非常に多くなり、正確に解析することは困難である。しかし、通常のＣＶＤ条件、特に本発明の対象であるダイヤモンドの合成条件においては、電子密度の空間変化に対して、イオン及び中性ガスの空間分布の時間変化は比較的緩やかであり、これらの非一様性の電子の密度、流速及び温度、並びに電磁界分布へ影響は少ない。従って、イオン及び中性ガスに関して一様性を仮定すれば、イオン及び中性ガスについて解く必要が無く、電子に関する変数のみを考慮すればよい。更に、以下の（１）、（２）に示すように、電子の流速及び温度を決定する方程式も解析的に解くことができる。
【００２３】
（１）電子に関する物理量を高周波振動成分と平均速度成分に分離することによって、電子の運動量輸送方程式が解析的に解け、電子の流速を他の物理量で表すことができる。
【００２４】
電子の運動量輸送方程式は、次式１のように表せる。
【００２５】
【数９】
2008031515.tif 000011
式１において、ｎeは電子密度、ｕeは電子流速（３次元ベクトル）、Ｔeは電子温度、
Ｅは電界（３次元ベクトル）、ｔは時間、νeは電子の衝突振動数、ｍeは電子質量である。
【００２６】
ここで、電子流速ｕeを平均速度成分ｕe0と高周波成分ｕe1との和（ｕe＝ｕe0＋ｕe1）で表すと、式２、式３が導出できる。ここで、通常、電子温度の勾配は電子密度勾配と比較して緩やかであることから、電子温度に起因する微分係数は無視した。
【００２７】
【数１０】
2008031515.tif 000012
式２、３において、ωは電磁界角周波数、ｅは電荷素量である。
【００２８】
（２）パラメータフィッティングを併用することによって、電子のエネルギー輸送方程式を解析的に解くことができ、電子温度を他の物理量で表すことができる。
【００２９】
上記非特許文献１では、定常状態における電子の質量輸送方程式を、式４で表されるとしている。
【００３０】
【数１１】
2008031515.tif 000013
式４の右辺第１項は電子と化学種（中性原子、中性分子、イオン等）との電離反応による電子の生成を表し、右辺第２項は電子と化学種との再結合反応による電子の消滅を表す。
【００３１】
式４は、非特許文献１で考慮されている化学反応のみに特化した場合の式であるが、一般には、生成項及び消滅項は式５、式６のように、可能な反応の総和として表される。ここで、可能な反応には、種類の異なる反応に限らず、同じ種類の反応であっても反応のエネルギー閾値（以下、単にエネルギー閾値とも記す）や反応率が異なる反応をも含む。即ち、例えば励起反応であっても、反応に関わる原子・分子種が異なる励起反応であれば、一般にはエネルギー閾値や反応率が異なるので、添え字ｓに一意の値を割り当てて別の反応として扱う。
【００３２】
【数１２】
2008031515.tif 000014
式５、式６において、sは反応の種類を示す添字であり、ｎsは反応種の密度であり、
【数１３】
2008031515.tif 000015
は、それぞれ電離反応の反応率、電離反応のエネルギー閾値、再結合反応の反応率、再結合反応のエネルギー閾値を表す。消滅項については、電荷準中性条件により、電子密度とイオン密度は局所的にほぼ等しいことから、式７を用いている。
【００３３】
【数１４】
2008031515.tif 000016
ここで、
【数１５】
2008031515.tif 000017
はイオン種ｓのモル分率を表す。
【００３４】
非特許文献１では、電離反応及び再結合反応を、ｅ＋Ｈ2→２ｅ＋Ｈ2+とこの逆反応（
２ｅ＋Ｈ2+→ｅ＋Ｈ2）のみに限定し、再結合反応のエネルギー閾値を０としているため
、式４の様になっている。式４の右辺を一般形である式５及び式６で置き換えると、式８を得る。
【００３５】
【数１６】
2008031515.tif 000018
この式は高周波振動成分を含まない平均化された質量輸送の釣り合いを示しているということが重要である。
【００３６】
以上の（１）、（２）によって得られた式８の左辺のｕeを式３の静的な平均成分ｕe0
で置き換え、以下の式９を得る。
【００３７】
【数１７】
2008031515.tif 000019
次に、以上で決定した方程式を用いて数値解析する方法を、図３に示したフローチャートに従って説明する。図３は、本発明の実施の形態に係るマイクロ波プラズマ解析プログラムを示すフローチャートである。
【００３８】
以下に説明する処理は、コンピュータを用いて微分方程式を数値的に積分する通常の処理に、多くの点で類似する処理である。また、以下において、特に断らない限り、コンピュータの演算処理手段（以下、ＣＰＵと記す）が行う処理として記載する。ＣＰＵは、メモリをワーク領域として使用して演算などの処理を行う。ＣＰＵが、処理を行う前提条件
は、外部から人によって操作手段（コンピュータ用キーボード、マウスなど）が操作されて、データとして入力される。入力されたデータは、メモリの所定領域に一時記憶され、演算処理に使用される。入力されたデータや演算処理の結果は、適宜記録部（ハードディスクなど）に記録される。また、以下においては、図１に示した構成のＣＶＤ装置を対象として説明する。
【００３９】
ステップＳ１において、数値計算の対象とする空間、及びその空間内に存在する物の形状に関する情報の入力を受け付ける。具体的には、人がコンピュータの操作手段を操作して、ＣＰＵが提供するグラフィカルユーザインタフェースを介してモニタ画面上で、図１に示した形状に関する情報を入力する。人が、各構成部（第１円筒部、第２円筒部、石英板、サセプタなど）の形状、寸法を指定し、対象とする空間中に各構成部品を配置することで、計算対象となる空間を区切る境界の位置、形状、寸法に関するデータが、メモリ上に生成される。生成されたデータは適宜、記録部（ハードディスク）に保存される。
【００４０】
ステップＳ２において、ステップＳ１で決定した形状をメッシュに分割する情報の入力を受け付ける。例えば、人が操作手段を操作して、各構成部の各辺を等間隔に分割する数を指定することによって、対向する分割点が線分で結ばれ、各構成部がメッシュに分割される。メッシュ分割に関する情報の入力は、各構成部のメッシュサイズを指定することで行うこともできる。
【００４１】
ステップＳ３において、境界条件の入力を受け付ける。具体的には、金属部分（図１では、第１円筒部、第２円筒部、ウェーブガイド、アンテナ、サセプタ）を完全導体とし、その表面に沿う方向の電界成分を０とする。また、電磁波が入力される部分（マイクロ波入力部）及び放電領域を指定する。
【００４２】
ステップＳ４において、計算の初期条件の入力を受け付ける。各構成部の比誘電率、電気伝導率、初期の電子密度ｎe(ini)の空間分布、入力される電磁波の周波数、投入電力、電磁波モードが入力される。また、考慮する各反応ｓの反応率ｋs及びエネルギー閾値Ｅs、並びに反応種の密度ｎsが入力される。また、背景ガスの圧力及び温度が入力される。
ここで、背景ガスとはプラズマ源のガスを意味し、背景ガスの圧力及び温度は、νeを評
価する時に使用される。その評価式は、Ｐgをガス圧、Ｔgをガス温度として、νｅ＝１．４４×１０12×Ｐg／Ｔgであり、これは文献W. Tan and T. A. Grotjohn, Diamond and Related Materials 4, 1145 (1995)等で用いられている評価式である。さらに、繰り返し
処理の収束性を判断するための基準パラメータｎthが入力される。
【００４３】
ここで、考慮する反応は、電離反応、励起反応、解離反応、再結合反応であるが、上記したように同じ種類の反応であってもエネルギー閾値や反応率が異なる場合の反応をも含み、添え字ｓに一意の値を割り当てて別の反応として扱う。
【００４４】
ステップＳ５において、パラメータフィッティングを行う。即ち、ステップＳ４において入力された反応の反応率ｋs及びエネルギー閾値Ｅs、並びに反応種の密度ｎsを用いて
、電子の温度Ｔeを所定範囲で変化させて、
【数１８】
2008031515.tif 000020
を計算し、得られたＴeに関するグラフ（図４参照）に、
【数１９】
2008031515.tif 000021
がフィッティングするように、即ち、式１０のグラフを式１１で近似できるように３つのパラメータａL、ｐL、ＥLを決定する。
【００４５】
図４は、パラメータフィッティングを説明するためのグラフであり、「×」印は式１０で計算された値であり、曲線は、決定されたパラメータを用いた式１１の曲線である。これらのパラメータが決定されれば、式１２（電子のエネルギー輸送方程式）から、式１３を用いて電子の温度Ｔeは式１４によって表される。ここで｜Ｅ｜は電界強度を表す。
【００４６】
【数２０】
2008031515.tif 000022
【００４７】
ステップＳ６において、電子密度ｎe（ｒ）から、ステップＳ３で指定した放電領域で
の誘電率εp（ｒ）を求める。ここで、誘電率εp（ｒ）は、電子密度ｎe（ｒ）を用いて
【数２１】
2008031515.tif 000023
で求める。後述する通り、ステップＳ６、Ｓ７、及びＳ８は、所定の条件を満足するまで繰り返される。条件を満足しない場合には、ステップＳ６に戻って再度１５式から誘電率εp（ｒ）を求めるが、この場合は後述するステップＳ８での処理により計算された電子
密度ｎe（ｒ）を用いる。ステップＳ９の判定を経ない、初めてステップＳ６において誘
電率εp（ｒ）を求める際は、ステップＳ４において指定した初期の電子密度ｎe(ini)の
空間分布から誘電率εp（ｒ）を求める。
【００４８】
ステップＳ７において、ステップＳ１〜Ｓ４で指定された固定条件（形状、境界条件）及び、ステップ６において求められた誘電率εp（ｒ）の下で、マクスウェルの方程式を
用いて、電界分布Ｅ（ｒ）を数値計算する。ここでＥ（ｒ）は、３次元の位置ベクトルｒで指定される位置における電界の３次元ベクトルである。電磁気に関するマクスウェルの方程式、及び、有限要素法などを用いて所定の境界条件の下でマクスウェルの方程式を数値積分し、所定空間内の電界分布Ｅ（ｒ）を求める方法は公知であるので、ここでは説明を省略する。
【００４９】
ステップＳ８において、ステップＳ７で求められた電界Ｅ（ｒ）（各メッシュにおける３次元ベクトル）を用いて、ステップＳ１〜Ｓ４で指定された固定条件（形状、境界条件
）の下で、上記で導出した式９の微分方程式を用いて、電子密度ｎe（ｒ）を数値計算に
よって求める。
【００５０】
【数２２】
2008031515.tif 000024
ここで、電子温度Ｔeは、ステップＳ７で求められた電界Ｅ（ｒ）を用いて、式１４に
よって計算される。また、各メッシュについて数値計算を行う前に、そのメッシュの電界強度Ｅ（ｒ）が０であるか否かを判断し、Ｅ（ｒ）＝０のメッシュについては数値計算をせずにｎe＝０とし、Ｅ（ｒ）≠０のメッシュについてのみ数値計算を行う。なお、有限
要素法などを用いて所定の境界条件の下で式９の微分方程式を数値積分する方法は公知であるので、ここでは説明を省略する。
【００５１】
ステップＳ９において、電子密度分布の計算結果の収束性を判断する。具体的には、ｋ回目、ｋ−１回目の電子密度分布の計算結果をｎek、ｎek-1で表し、ステップＳ４で入力された基準パラメータｎthを用いて、
【数２３】
2008031515.tif 000025
を満たせば、収束したと判断してステップＳ１１に移行し、満たしていなければステップＳ６に戻る。ここで、τの積分は放電領域での体積積分を表す。ステップＳ６に戻った場合、誘電率εp（ｒ）を求める際は、ステップＳ８で計算された電子密度ｎe（ｒ）を用いる。
【００５２】
ステップＳ１０において、最新の計算結果である電子密度分布ｎe（ｒ）を、例えばモ
ニタ画面にコンピュータグラフィック（以下、ＣＧと記す）によって描画することによって、提示する。
【００５３】
最後にステップＳ１１において、終了の指示があったか否かを判断し、終了の指示があるまで、ステップＳ１〜Ｓ１０の処理を繰り返す。
【００５４】
以上によって、ステップＳ１０で提示されたＣＧ画像を、人が目視によって評価し、ＣＶＤ装置の形状及び寸法設計の妥当性を評価することができる。
【００５５】
また、ステップＳ８における数値計算は、１つの方程式を用いた処理でよいので、従来よりも高速にマイクロ波プラズマの分布を計算することが可能である。
【００５６】
以上では、ステップＳ５においてパラメータフィッティングを行う場合を説明したが、予めパラメータフィッティングを行い、パラメータａL、ｐL、ＥLを計算しておいてもよ
い。その場合には、ステップＳ４での初期値の入力において、予め求められたパラメータaL、ｐL、ＥLを入力するようにすれば、ステップＳ５の処理が不要になる。特に、対象とする反応を特定種類の反応に固定してシミュレーションを繰り返す場合には、パラメータａL、ｐL、ＥLを一度だけ計算すればよいので、全体としての処理時間をより短縮するこ
とができる。
【００５７】
上記した一連の処理を行うコンピュータプログラムを構成する場合、図５に示したよう
に、誘電率を計算（ステップＳ６）するモジュール、マクスウェル方程式を用いて電界を数値計算（ステップＳ７）するモジュール、及び式９を用いて電子密度を数値計算（ステップＳ８）するモジュールを、それぞれ別のモジュールに構成することができる。この様に構成することによって、電子密度及び誘電率を計算するモジュールだけを新たにプログラミングすれば、それ以外の処理に関しては、公知の電磁界解析用ソフトウェアを使用することができる。
【実施例】
【００５８】
以下に、実施例を示し、本発明の特徴をより明確にする。本実施例では、図１に示すＣＶＤ装置を想定し、本発明のマイクロ波プラズマ解析プログラムによる数値計算を行った。
【００５９】
図６は、各部の寸法を付した断面図である。メッシュサイズ、即ち、１つの正方形メッシュの１辺の長さは、サセプタ天面中央上方の、直径４ｃｍ高さ３ｃｍの円柱領域（図６の点線で囲んだ領域）は２ｍｍ、この領域以外は１ｃｍに設定した。
【００６０】
境界条件に関しては、黒色の構成部（第１円筒部、第２円筒部、アンテナ、サセプタ、ウェーブガイド）は全て完全導体とし、表面に沿う電界成分を０とし、電磁波はマイクロ波入力部から導入されるように設定した。
【００６１】
放電領域は、石英板の下側の第２円筒部内である。マイクロ波入力部から入力される電磁波は、周波数が２．４５ＧＨｚ、投入電力が１６００Ｗ、電磁波モードが最低次のＴＥモードとした。
【００６２】
初期の電子密度分布は、円盤状のサセプタの中心からの距離をｒ（ｍ）とし、サセプタ上方の放電領域において、ｒ＜０．０２である領域では、ｎe(ini)＝２．５×１０17×[
１．０＋ｃｏｓ（πｒ／０．０２）]であり、ｒ≧０．０２の領域では、ｎe(ini)＝０で
ある分布とした。
【００６３】
各構成部の固定の誘電率及び電気伝導率については、石英板が比誘電率＝３．８、電気伝導率＝１ｅ-17ｍ／Ωとし、空気領域は比誘電率＝３．８、電気伝導率＝０ｍ／Ωであ
るとした。
【００６４】
プラズに関わるパラメータに関しては、表１の４種類の反応を想定し、電子密度分布の収束条件として、ｎth＝０．０１とした。
【００６５】
【表１】
2008031515.tif 000026
また、フィッティングパラメータには、予め計算した値である、ＥL＝１８．５３２７
（ｅＶaL）、ｐL＝２．４６４×１０15（ｅＶ／ｓｅｃ）、ａL＝０．３６７５７を用いた。また、背景ガスの圧力は１６０Torr、温度は１０3Ｋとした。また、簡単のため、反応
種密度ｎsは、圧力１６０Torr、温度１０3Ｋからｎs＝１．５×１０24とした。
【００６６】
以上の条件で、上記したステップＳ６〜Ｓ９の処理を実行した結果を図７に示す。図７の左半分には、図６の破線領域の電子密度を、右半分には図６の破線領域のｐabを示している。何れも、明度が高いほど（白いほど）値が大きい。また、図８は、プラズマ発生領域の電界強度の分布を示す図である。図７、図８には、各値が最大となる場所を黒丸で示し、その最大値を数値で記載している。
【００６７】
図７、図８から、本シミュレーションによって計算されたパワー密度が、実験から見積もられるパワー密度に一致していることが分かった。また、電子密度が最大となる位置（図７参照）と電界強度が最大となる位置（図８参照）とが異なっており、上記非特許文献１では再現できなかった、高密度プラズマによって電磁界が遮蔽される効果を、本実施例で再現することができた。このように、本実施例によって得られた分布の傾向は、実際の放電の様子や、他のより大規模なシミュレーションを行った結果と一致していた。
【図面の簡単な説明】
【００６８】
【図１】本発明が解析対象とするＣＶＤ装置のプラズマ発生部の一例を示す断面図である。
【図２】本発明が解析対象とするＣＶＤ装置のプラズマ発生部の別の例を示す断面図である。
【図３】本発明の実施の形態に係るマイクロ波プラズマ解析プログラムを示すフローチャートである。
【図４】本発明の実施の形態に係るマイクロ波プラズマ解析プログラムにおけるパラメータフィッティングを説明するための図である。
【図５】本発明の実施の形態に係るマイクロ波プラズマ解析プログラムの構成の一例を示すブロック図である。
【図６】本発明の実施例を示すＣＶＤ装置のプラズマ発生部の構造を示す断面図である。
【図７】数値計算の結果を示す図である。
【図８】数値計算の結果を示す図である。
【符号の説明】
【００６９】
１    マイクロ波入力部
２    ウェーブガイド
３    第1円筒部
４    石英板
５    第２円筒部
６    真空領域
７    アンテナ
８    サセプタ
【図１】
2008031515.tif 000035
【図２】
2008031515.tif 000036
【図３】
2008031515.tif 000037
【図４】
2008031515.tif 000038
【図５】
2008031515.tif 000039
【図６】
2008031515.tif 000040
【図７】
2008031515.tif 000041
【図８】
2008031515.tif 000042