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